nowadays , many x - ray imaging experiments are performed dynamically , meaning that the object s inner structure is a function of time . in material science , the sample is often exposed to various external conditions ( e.g. temperature , pressure , etc . ) during the time of a scan . in medicine
dynamic tomographic imaging is focused on capturing various temporal changes under realistic conditions [ 13 ] .
good examples include the corrosion and oxidation processes in metals , crystal growth , crack healing   or fluid flow in a fixed porous body [ 2 , 7 , 8 , 9 ] .
when an object is scanned dynamically , two main sources of motion should be taken into account .
firstly , motion related to movement or change occurring within the object . for instance , in medical x - ray imaging , a patient s cyclic breathing pattern should be considered . in this case ,
due to the repeatability , this effect can be accounted for ( some motion artifacts can be eliminated through the gating procedure ) . in material science , however , changes can happen on multiple scales ( from coarse to fine ) , severely fragmented , independent and irregular .
this makes it more difficult and in many cases impossible to predict the expected motion .
one way of reducing motion artifacts is to scan faster and/or collect less projection data .
such strategies can lead to the second source of motion , those which are related to the acquisition process itself . nowadays
, synchrotron imaging can provide probably the fastest exposure times and rotation speeds [ 1012 ] .
consequently , the expected structural changes can be significantly minimized and the obtained projection data can be considered as acquired from an essentially stationary object , although it is a crude approximation in general .
the extent of structural differences between time frames depends on how fast the object s structure changes during each scan and how rapid is the acquisition hardware .
this , however , is not always possible , e.g. the imaging of fluid flow [ 2 , 3 , 7 , 8 ] .
additional difficulties are related to the limited exposure time and the speed of sample rotation . if the angle of radial integration increased ( e.g. to reduce the number of projections collected per scan ) it leads to blurring artifacts .
those artifacts , up to a certain degree , can be eliminated with state - of - the - art reconstruction algorithms which model the directional blurring effect in the projection space . in this work
we do not account for motion related issues , assuming that the exposure time is short and the object is ( almost ) stationary during the time of scan .
we refer to such an approach as four dimensional ( 4d ) tomography ( three dimensional spatial coordinates plus time evolution ) , where 3d  snapshots 
frequently , fast dynamic imaging is associated with a limited number of projections acquired per scan and the reconstruction problem is therefore under - determined and ill - conditioned . due to poor - conditioning , direct inversion methods fail to work and lead to a strong amplification of noise and loss of resolution in the reconstructed images . iterative image reconstruction ( iir ) techniques can approximate inversion by using the error - correcting refinements resulting in improved signal - to - noise ratio ( snr ) characteristics of reconstructed images .
however , the iir for undersampled data is an ill - posed problem and it requires some form of regularity applied to the solution to ensure convergence [ 14 , 15 ] . traditionally , the spatial regularity of smoothness ( 2 norm ) or gradient sparsity ( 1 norm ) are employed .
the spatial regularization , however , is not well fitted to the dynamic nature of 4d imaging since it completely discards the redundancy of the temporal information . in this paper
, we propose a spatio - temporal regularization approach for iir which is specifically adapted to various cases of dynamic imaging where the motion pattern is unknown and/or difficult to predict .
previously , there have been various attempts to improve spatial and temporal resolution in time - lapse tomography by employing a supplementary image as a prior in the reconstruction algorithm [ 79 , 16 , 17 ] .
for instance , in order to reconstruct fluid flow propagation through a fixed porous body ( e.g. a rock sample ) , a high resolution pre - scan of the solid porous object is used to improve spatial resolution [ 79 ] .
similarly , in medical imaging , a previous scan of the same patient can be also employed [ 16 , 17 ] .
although these methods can be highly effective in exploiting supplementary data , they can not be used in many other cases when a supplementary image is not available or uninformative .
the majority of time - lapse tomographic experiments in xmt are fully dynamic and without a fixed relationship to the initial state of the object .
therefore , it is difficult to extract any supplementary ( useful ) information from the data obtained prior to the scan .
various spatio - temporal regularization techniques have been developed recently and the majority of them are based on extracting information from the adjacent time frames only [ 1820 ] .
this ( local ) approach assumes that the closest time frames have a higher probability to be structurally similar to the regularized frame ( which is true for many fast dynamic experiments ) . in our approach
, we assume that the useful information ( e.g. structurally similar features ) might exist in much more distant time frames than simply adjacent ones ( depending on the experiment ) . the core problem is how to identify and use this information in very large temporal xmt data . in order to realize this we use nonlocal regularization on weighted graphs , which can be considered as a generalized case for nonlocal regularization [ 22 , 23 ] . while spatio - temporal regularization has been done before
we demonstrate a novel methodology for application to large datasets ( more than 4    10 voxels size ) characteristic of xmt .
the novelty of our method consists in the generalized form of the spatio - temporal penalty and the special data reduction technique to significantly accelerate calculations .
a speed - up of an order of magnitude is achieved which makes our method more feasible for reconstruction of large datasets .
the modelled numerical experiment is presented to demonstrate that the proposed method is able to provide better contrast and resolution than the classical nonlocal approach while performing ten times faster . to further demonstrate the applicability of our method it is applied to reconstruction of big 4d datasets characteristic of xmt .
let xnk , where n is the total number of image elements and k is the number of all time - frames , then x:=(x1t , x2t, 
similarly we define the measured projections vector as b:=(b1t , b2t,  ,bkt)t and noise component :=(1t,2t,  ,kt)t . here
therefore the system of linear equations to solve is : 
 ( 1)[a10  00a2000  ak][x1tx2txkt]=[b1tb2tbkt]+[1t2tkt ] 
  note that the block diagonal matrix amknk is time - invariant for many 4d imaging set - ups , that is , a1  =  a2  =    =  ak or a  =  i    a1 , where i is the identity and  is the kronecker product .
this , however , is not optimal for reconstruction of dynamically collected projections and more advanced scanning strategies can be employed [ 3 , 8 , 20 ] .
sparse projection matrix a:nm describes the scanning geometry , it is parallel beam in our case , however a can potentially accommodate more complex geometries , such as fan - beam or cone - beam . the regularized unconstrained least - squares ( ls ) problem can be formulated as : 
 ( 2)x=argminx{12ax - b22+r(x ) } , 
 where r 
( 
 x ) is the proposed spatio - temporal ( st ) penalty term and  is a regularization parameter which establishes a trade - off between the data and the regularization terms .
similarly to [ 7 , 8 ] , we use a splitting approach [ 23 , 25 ] to decouple data fidelity and regularization terms into two simpler sub - problems to solve : 
 ( 3){vn+1=xn-[at(axn - b)]xn+1=argminx(r(x)+2x - x022 ) ;  x0=vn+1 . 
  the first step in the forward - backward splitting ( fbs ) algorithm ( 3 )
solves the unregularized ls problem ( a gradient descent ( gd ) iteration with a time - step parameter  ) and the second is a weighted image denoising step . due to slow convergence of the gd method
there have been successful attempts to replace it with faster in convergence methods [ 7 , 8 , 18 ] .
similarly , a conjugate gradient least squares ( cgls ) algorithm [ 14 , 15 ] is used instead of the gd method in this paper . in practice , this substitution works well and similar results can be obtained with the cgls - based fbs method , however , the mathematical proof of convergence does not hold anymore .
the main value of our method is in the original form of the st penalty r  ( 
 x ) which is outlined in the next section .
by taking a nonlocal regularization approach based on weighted graphs   we can rewrite the second minimization problem in ( 3 ) as : 
 ( 4)x=argminx(x ) ,  (x)=1pvvx(v)p+2x - x022 , 
 where v is a set of vertices v of a weighted graph g  =  ( v , e ,  ) , p  ]0 , +    [ is a smoothness degree ,  is a similarity weight function defined on edges e and   x  ( v )   is the weighted local variation of x over the graph .
the -weighted p - laplace term in ( 4 ) is convenient to operate with since various choices of p lead to well - known nonlocal filters used in image processing [ 22 , 23 ] .
for instance , p  =  2 leads to the weighted linear heat diffusion on graphs and p  =  1 to the mean curvature or total variation ( tv ) seminorm on graphs .
two vertices u and v are said to be adjacent if the edge ( u , v )    e connects them , the weight between u and v is denoted by   ( u , v ) :   ( u , v )    e ,   ( u , v )  =    ( v , u ) . also   ( u , v )  > 0 if u    v and 0 otherwise . notation u 
the vertex u    b  ( v ) , where b  ( v ) is a cubic box ( the searching space ) size of sx    sy    sk .
potentially , the box size can be extended to the 4d case ( sx    sy    sz    sk ) , however due to significant computational costs involved it is currently infeasible for large xmt data .
therefore , our regularization is currently in 3d ( x , y+time ) and it has similarities with the video denoising problem .
the major benefit of this realization is that it is relatively fast and can be programmed in a 
slice - by - slice  manner which does not require lots of computer memory for reconstruction and is easily parallelizable .
the local variation of the weighted gradient operator   x  ( v )   is defined by 
 ( 5)x(v)=uv(u , v)(x(v)-x(u))2 . 
 
the weighted p - laplace operator at a vertex v is defined as 
 ( 6)px(v)=1puv(u , v)(x(v)-x(u ) ) , 
 where 
 ( 7)(u , v)=(u , v)(x(v)p-2+x(u)p-2 ) . 
 
one can see that for the case p  =  2 :   ( u , v )  =    ( u , v ) and for p  =  1 it is : 
 ( 8)(u , v)=(u , v)(1x(v)+1x(u) ) . 
 
both cases are convex and therefore lead to the convex minimization problem ( 4 ) , which guarantees a unique solution x under a proper choice of parameter . in this paper we use p  =  1 which provides a better regularization performance than the linear model p  =  2 .
note that the nonlocal tv seminorm ( 5 ) is smoothed by a small constant to ensure differentiability : x(v)=uv(ux(v))2+2 .
the non - convex case p  <  1 is not considered in this paper .
therefore for p    1 the system ( 4 ) has a unique solution written as : 
 ( 9)2px(v)+(x - x0)=0 ,  vv , 
 which is equivalent to the following system of algebraic equations : 
 ( 10)(+uv(u , v))x(v)-uv(u , v)x(u)=x0(v ) . 
 
the linearized gauss - jacobi iterative method can be used to solve the system ( 10 ) which also can be considered as the discrete euler - lagrange equation of minimization problem ( 4 ) [ 2123 ] .
let t be an iteration step , then the fixed point ( fp ) iteration is used to solve ( 10 ) : 
 ( 11)xt+1(v)=x0(v)+uvt(u , v)xt(u)+uvt(u , v ) ,  vv . 
  in our case ( p  =  1 )
,   ( u , v ) as taken as in ( 8) . the iteration procedure ( 11 ) can be stopped using criterion 
 x  - 
 x    <   , where  is a small constant . in our experiments
we use only one iteration of ( 11 ) mainly due to regularization performance restrictions .
the regularization on graphs for time - lapse reconstruction involves a search for similar vertices in a box b  ( v ) .
similar to nonlocal video denoising , a quadratic patch f size of rx    ry is used to provide more robust estimation of similarity between vertices u and v ( calculation of euclidian distance on patches [ 22 , 23 ] ) . the computational complexity of the proposed algorithm can go up to o(n ) and therefore becomes infeasible for large xmt datasets . here
we propose a simple approach to significantly accelerate the regularization process on graphs by considering sparsity in the temporal domain .
we developed an accelerated regularization on graphs ( arg ) technique based on two observations : a ) the temporal data are likely to be sparse and redundant ( some features remain stationary or change only a little over time ) ; b ) two vertices are likely to be similar ( structurally ) if the difference in intensity between them is lower than the maximum level of noise present in the signal .
we use these two observations to modify and accelerate the classical regularization algorithm on graphs ( rg ) .
we say that the sx , y size of the searching box b  ( v ) must be larger if a feature at the vertex v is changing ( e.g. moving ) quickly over time and sx , y of b  ( v ) must be smaller if the feature is almost stationary ( e.g. a uniform region or an edge ) .
by reducing b  ( v ) for temporally non - active ( stationary ) features we drastically reduce our search space . in order to decide if the feature is stationary or not and if the vertex u is similar to v we use the observation ( b ) .
while minimizing the cost function ( 2 ) , we intentionally do not impose any positivity constraints on our solution and therefore 
 x will have negative values on every n - iteration . the maximum noise level present in our reconstructed image
can be simply found as   =  |  min ( 
 x )  | .
we also say that the trustworthy intensity differences between vertices u and v are bounded above by l , where l    ( 0 , 1 ] is an empirical constant ( we use l  =  0.4 for all our experiments ) . to establish the optimal searching box size bo 
( v ) for v  ( i , j , k )  ,  k  =  1 , 2 ,  ,
k we calculate the following measure of temporal sparsity : 
 ( 12)s(v)=s(v)+1 ,  if  |xn(v(i , j , k))-xn(u(i , j , l))|ln ,  l=1,2,  ,k , 
 where x is a quadratic average ( mean ) of intensity values in 
 x in order to reduce the influence of noise in estimation . the sparsity measure s    ( 0 , k  -  1 ] ( 12 ) , where s  =  0 corresponds to the most temporally non - sparse v - region and s  =  k  -  1 to the most temporally sparse v - region . in fig .
1 ( middle ) we demonstrate the calculated measure 
 s for time frame k  =  4 ( left ) of the noisy video sequence of k  =  20 frames .
the level of temporal sparsity can be clearly seen ( dynamic features , like ship and birds , are not sparse , while stationary features are ) .
then , if we select the lower size of the box bo  ( v ) as sxl = syl and the upper size as sxu = syu , s  ( v ) can be linked to bo  ( v ) ( see fig .
1 ( right ) ) . we divide the whole interval min(sk )  -  max(sk ) into an equal number of d steps ( we use d  =  10 in our experiments ) and linearly assign values for bo  ( k ) in a way that min(sk ) corresponds to s and max(sk ) to s respectively .
the simple reasoning here is that the larger 
 s values should be related to the smaller sizes of the searching window and vice versa .
finally , the following condition is used to calculate nonlocal weights in ( 11 ) : 
 ( 13)n(u , v)={exp(-f(xn(v))-f(xn(u))h2),if|xn(v)-xn(u)|ln0,otherwise . 

 
one can see that we reduce the data space substantially by first creating a spatially variant search box bo and next we remove vertices which do not fit our assumptions of similarity using the rule ( 13 ) . as we will demonstrate later
, this approach provides an acceleration by an order of magnitude without impairment of reconstruction quality .
moreover , the averaging across dissimilar features in arg is much reduced and leads to improved resolution .
by taking a nonlocal regularization approach based on weighted graphs   we can rewrite the second minimization problem in ( 3 ) as : 
 ( 4)x=argminx(x ) ,  (x)=1pvvx(v)p+2x - x022 , 
 where v is a set of vertices v of a weighted graph g  =  ( v , e ,  ) , p  ]0 , + 
  [ is a smoothness degree ,  is a similarity weight function defined on edges e and   x  ( v )   is the weighted local variation of x over the graph .
the -weighted p - laplace term in ( 4 ) is convenient to operate with since various choices of p lead to well - known nonlocal filters used in image processing [ 22 , 23 ] .
for instance , p  =  2 leads to the weighted linear heat diffusion on graphs and p  =  1 to the mean curvature or total variation ( tv ) seminorm on graphs .
two vertices u and v are said to be adjacent if the edge ( u , v )    e connects them , the weight between u and v is denoted by   ( u , v ) :   ( u , v )    e ,   ( u , v )  =    ( v , u ) . also   ( u , v )  > 0 if u    v and 0 otherwise .
the vertex u    b  ( v ) , where b  ( v ) is a cubic box ( the searching space ) size of sx    sy    sk .
potentially , the box size can be extended to the 4d case ( sx    sy    sz    sk ) , however due to significant computational costs involved it is currently infeasible for large xmt data .
therefore , our regularization is currently in 3d ( x , y+time ) and it has similarities with the video denoising problem .
the major benefit of this realization is that it is relatively fast and can be programmed in a 
slice - by - slice  manner which does not require lots of computer memory for reconstruction and is easily parallelizable .
the local variation of the weighted gradient operator   x  ( v )   is defined by 
 ( 5)x(v)=uv(u , v)(x(v)-x(u))2 . 
 
the weighted p - laplace operator at a vertex v is defined as 
 ( 6)px(v)=1puv(u , v)(x(v)-x(u ) ) , 
 where 
 ( 7)(u , v)=(u , v)(x(v)p-2+x(u)p-2 ) . 
 
one can see that for the case p  =  2 :   ( u , v )  =    ( u , v ) and for p  =  1 it is : 
 ( 8)(u , v)=(u , v)(1x(v)+1x(u) ) . 
 
both cases are convex and therefore lead to the convex minimization problem ( 4 ) , which guarantees a unique solution x under a proper choice of parameter . in this paper we use p  =  1 which provides a better regularization performance than the linear model p  =  2 .
note that the nonlocal tv seminorm ( 5 ) is smoothed by a small constant to ensure differentiability : x(v)=uv(ux(v))2+2 .
the non - convex case p  <  1 is not considered in this paper .
therefore for p    1 the system ( 4 ) has a unique solution written as : 
 ( 9)2px(v)+(x - x0)=0 ,  vv , 
 which is equivalent to the following system of algebraic equations : 
 ( 10)(+uv(u , v))x(v)-uv(u , v)x(u)=x0(v ) . 
 
the linearized gauss - jacobi iterative method can be used to solve the system ( 10 ) which also can be considered as the discrete euler - lagrange equation of minimization problem ( 4 ) [ 2123 ] .
let t be an iteration step , then the fixed point ( fp ) iteration is used to solve ( 10 ) : 
 ( 11)xt+1(v)=x0(v)+uvt(u , v)xt(u)+uvt(u , v ) ,  vv . 
  in our case ( p  =  1 ) ,   ( u , v ) as taken as in ( 8) . the iteration procedure ( 11 ) can be stopped using criterion 
 x  - 
 x    <   , where  is a small constant . in our experiments
we use only one iteration of ( 11 ) mainly due to regularization performance restrictions .
the regularization on graphs for time - lapse reconstruction involves a search for similar vertices in a box b  ( v ) .
similar to nonlocal video denoising , a quadratic patch f size of rx    ry is used to provide more robust estimation of similarity between vertices u and v ( calculation of euclidian distance on patches [ 22 , 23 ] ) .
the computational complexity of the proposed algorithm can go up to o(n ) and therefore becomes infeasible for large xmt datasets . here
we propose a simple approach to significantly accelerate the regularization process on graphs by considering sparsity in the temporal domain .
we developed an accelerated regularization on graphs ( arg ) technique based on two observations : a ) the temporal data are likely to be sparse and redundant ( some features remain stationary or change only a little over time ) ; b ) two vertices are likely to be similar ( structurally ) if the difference in intensity between them is lower than the maximum level of noise present in the signal .
we use these two observations to modify and accelerate the classical regularization algorithm on graphs ( rg ) .
we say that the sx , y size of the searching box b  ( v ) must be larger if a feature at the vertex v is changing ( e.g. moving ) quickly over time and sx , y of b  ( v ) must be smaller if the feature is almost stationary ( e.g. a uniform region or an edge ) .
by reducing b  ( v ) for temporally non - active ( stationary ) features we drastically reduce our search space . in order to decide if the feature is stationary or not and if the vertex u is similar to v we use the observation ( b ) .
while minimizing the cost function ( 2 ) , we intentionally do not impose any positivity constraints on our solution and therefore 
 x will have negative values on every n - iteration . the maximum noise level present in our reconstructed image can be simply found as   =  |  min ( 
 x )  | .
we also say that the trustworthy intensity differences between vertices u and v are bounded above by l , where l    ( 0 , 1 ] is an empirical constant ( we use l  =  0.4 for all our experiments ) . to establish the optimal searching box size bo  ( v ) for v  ( i , j , k )  ,  k  =  1 , 2 ,  , k we calculate the following measure of temporal sparsity : 
 ( 12)s(v)=s(v)+1 ,  if  |xn(v(i , j , k))-xn(u(i , j , l))|ln ,  l=1,2,  ,k , 
 where x is a quadratic average ( mean ) of intensity values in 
 x in order to reduce the influence of noise in estimation . the sparsity measure s    ( 0 , k  -  1 ] ( 12 ) , where s  =  0 corresponds to the most temporally non - sparse v - region and s  =  k  -  1 to the most temporally sparse v - region . in fig .
1 ( middle ) we demonstrate the calculated measure 
 s for time frame k  =  4 ( left ) of the noisy video sequence of k  =  20 frames .
the level of temporal sparsity can be clearly seen ( dynamic features , like ship and birds , are not sparse , while stationary features are ) .
then , if we select the lower size of the box bo  ( v ) as sxl = syl and the upper size as sxu = syu , s  ( v ) can be linked to bo  ( v ) ( see fig .
1 ( right ) ) . we divide the whole interval min(sk )  -  max(sk ) into an equal number of d steps ( we use d  =  10 in our experiments ) and linearly assign values for bo  ( k ) in a way that min(sk ) corresponds to s and max(sk ) to s respectively .
the simple reasoning here is that the larger 
 s values should be related to the smaller sizes of the searching window and vice versa .
finally , the following condition is used to calculate nonlocal weights in ( 11 ) : 
 ( 13)n(u , v)={exp(-f(xn(v))-f(xn(u))h2),if|xn(v)-xn(u)|ln0,otherwise . 

 
one can see that we reduce the data space substantially by first creating a spatially variant search box bo and next we remove vertices which do not fit our assumptions of similarity using the rule ( 13 ) . as we will demonstrate later
, this approach provides an acceleration by an order of magnitude without impairment of reconstruction quality .
moreover , the averaging across dissimilar features in arg is much reduced and leads to improved resolution .
in this section we present some numerical results with the method given in section 2 . in [ 7 , 8 ]
we have shown that the st regularization which uses more than just adjacent time - frames is able to outperform spatial and some spatio - temporal approaches . when more time - frames are used in regularization
therefore in this paper we focus on two issues : first is the quantitative and qualitative performance estimation of the arg method ( shorter abbreviation for cgls - arg ) in comparison to the normal rg ( cgls - rg ) method .
we provide an omp - c function which is mex - wrapped in matlab environment to perform optimized calculation for arg and rg techniques .
we note that it is a 3d version ( x , y + time ) which is well suited for parallel beam slice - by - slice reconstruction . for smaller data sizes
, the 4d implementation can be also feasible and it should provide images with even better snr and resolution . for calculations we used intel xeon cpu 2.26 ghz with 2 processors ( 4 threads each ) . to test arg and rg approaches we created a dynamically evolving phantom from high quality tomographic measurements ( the phantom is available here ) .
the high resolution data of a mouse tibia sample was acquired at the diamond - manchester branchline ( i13 ) of the diamond light source ( dls ) . in order to simulate the dynamic behaviour of the stationary sample we take 10 vertical slices of 3d reconstructed volume and shift them horizontally using an equidistant step of 15 pixels ( see fig .
three time frames ( k  =  1 , 5 , 10 ) of the mouse tibia phantom are shown in fig .
the total size of the phantom is 400    400    10 
( k ) pixels . to avoid the  inverse crime  of reconstructing on the same grid where projection data was simulated
, we used a higher resolution of the phantom on a 800    800 isotropic pixel grid to generate projections with a strip kernel .
gaussian noise was applied to the projection data ( 5% of a true signal maximum ) .
reconstructions were calculated on a smaller 400    400 isotropic pixel grid and with a linear projection model .
we used 180 and 90 projection angles in [ 0 ,  ) angular interval ( assuming parallel beam geometry ) to reconstruct the phantom . the root mean square error ( rmse )
is calculated as 
 ( 14)rmse(x , x)=e(x - x)2 , 
 where 
 x is an exact image sequence and x is a reconstructed one . before reconstructing the data we found optimal parameters for each method by optimizing rmse for each parameter while keeping other parameters fixed ( see table 1 ) . in fig .
4 we show reconstructions ( for time frame k  =  5 of the phantom in fig .
2 ( middle ) ) with various methods for two cases of limited noisy projections ( 180 and 90 projections ) . the main aim here is to show that the accelerated arg method can be competitive ( quantitatively ) with the classical rg method , while much faster computationally . to reconstruct the data we used various sizes ( sx  =  sy ) of the searching window bsx , sy , sk , while sk  =  9 remains constant ( we used nine temporal frames for every time - frame in regularizion ) . because our dynamic phantom ( see fig .
2 ) has significant vertical shifts in time , one should use larger search window sizes ( sx  =  sy ) to track the similar edge .
the rmse was calculated for all reconstructed phantoms and given in tables 2 and 3 .
one can see that the reconstructions with the cgls method are noisy ( see fig .
the regularization on graphs ( rg ) is able to reduce noise significantly and sharpen important features .
however , the small search window reconstructions ( method rg - b9 , 9 , 9 ) have a distinctive blocky appearance ( insufficient amount of temporal information collected ) . with a larger search window ( rg - b43 , 43 , 9 ) , reconstructions are much smoother and edges are well preserved .
unfortunately , the computational time is very large for rg - b43 , 43 , 9 ( see tables 2 , 3 ) and therefore the method is impractical . the computation time is drastically reduced with the proposed method ( arg - bo,9 ) while the error is the lowest .
notably , the quality of reconstruction is very similar to the original method ( rg - b43 , 43 , 9 ) , however some small and thin features are better preserved than with rg - b43 , 43 , 9 .
the effect of under and over estimation for various methods can be clearly seen in a 1d plot shown in fig .
5 ( by arrows we show the most erroneous regions in the profile recovered by methods rg - b9 , 9 , 9 and rg - b43 , 43 , 9 ) .
notably , the rg - b43 , 43 , 9 method smooths small features while the proposed method recovers them much better ( see the parabola shaped region in fig .
we conclude that the proposed method ( arg - bo ) is more feasible for larger datasets due to performance and quality - wise it is similar or better than the classical nonlocal regularization approaches .
the main purpose of developing such a spatio - temporal regularization technique ( see section section 2 ) is to apply it to time - lapse xmt data . due to high resolution detectors and rapid acquisition
( 
 x ) is very important for iterative minimization of ( 2 ) . in order to perform fast projection - backprojection operations we use gpu - accelerated open - source package : the astra toolbox .
the proposed regularizer is implemented in c , however it can be also realized on a gpu which can provide substantial acceleration .
the following time - lapse experiment was performed on the i13 beamline of dls , where multiple fast 3d scans of time - varying sample ( melting - freezing ice - cream ) were obtained .
the ice - cream microstructure consists of three phases , air cells ( gas ) , ice crystals ( solid ) and unfrozen matrix ( liquid ) ( see fig . 6 , 7 ) .
the aim of the experiment was to understand how phase content varies with temperature variation . during the thermal cycling the inner structure of a sample is changing ( crystals melt or coarsen ) .
since dynamic changes are easily traceable ( controlled sample cooling and heating procedure ) , this dataset is a perfect case study for our reconstruction technique .
additionally , since our reconstruction method searches for some structural correlations in time it significantly enhances the edges between different phases . due to poor contrast between ice - crystals and air bubbles one has to recover the edges of the unfrozen matrix as accurately as possible . in fig .
6 one can see reconstructions with three different methods : fbp ( filtered back - projection ) , cgls and the proposed arg method .
the time - lapse data size is 2k    2k    1k    6  ( k ) voxels , and therefore 6 large 3d volumes are reconstructed .
we perform 40 iterations of cgls method with one fp iteration ( 11 ) for the arg penalty .
all parameters for the arg method except  and h are selected the same as in table 1 .
the computation time of one arg iteration ( 11 ) to regularize 6 slices ( 2k    2k    6  ( k ) ) is 170 seconds ( all available time frames are used ) .
the classical rg implementation is not used here due to long computation times ( more than 1500 seconds for one iteration ) . from the reconstructions
it can be seen that the direct method ( fbp ) and cgls both deliver unquantifiable reconstructions due to high noise levels and ring artifacts present .
that the arg method is able to suppress noise significantly while also enhance edges between different phases .
to test arg and rg approaches we created a dynamically evolving phantom from high quality tomographic measurements ( the phantom is available here ) .
the high resolution data of a mouse tibia sample was acquired at the diamond - manchester branchline ( i13 ) of the diamond light source ( dls ) . in order to simulate the dynamic behaviour of the stationary sample we take 10 vertical slices of 3d reconstructed volume and shift them horizontally using an equidistant step of 15 pixels ( see fig .
three time frames ( k  =  1 , 5 , 10 ) of the mouse tibia phantom are shown in fig .
the total size of the phantom is 400    400    10 
( k ) pixels . to avoid the  inverse crime  of reconstructing on the same grid where projection data was simulated
, we used a higher resolution of the phantom on a 800    800 isotropic pixel grid to generate projections with a strip kernel .
gaussian noise was applied to the projection data ( 5% of a true signal maximum ) .
reconstructions were calculated on a smaller 400    400 isotropic pixel grid and with a linear projection model .
we used 180 and 90 projection angles in [ 0 ,  ) angular interval ( assuming parallel beam geometry ) to reconstruct the phantom . the root mean square error ( rmse )
is calculated as 
 ( 14)rmse(x , x)=e(x - x)2 , 
 where 
 x is an exact image sequence and x is a reconstructed one . before reconstructing the data we found optimal parameters for each method by optimizing rmse for each parameter while keeping other parameters fixed ( see table 1 ) . in fig .
4 we show reconstructions ( for time frame k  =  5 of the phantom in fig .
2 ( middle ) ) with various methods for two cases of limited noisy projections ( 180 and 90 projections ) . the main aim here is to show that the accelerated arg method can be competitive ( quantitatively ) with the classical rg method , while much faster computationally . to reconstruct the data we used various sizes ( sx  =  sy ) of the searching window bsx , sy , sk , while sk  =  9 remains constant ( we used nine temporal frames for every time - frame in regularizion ) . because our dynamic phantom ( see fig .
2 ) has significant vertical shifts in time , one should use larger search window sizes ( sx  =  sy ) to track the similar edge .
the rmse was calculated for all reconstructed phantoms and given in tables 2 and 3 .
one can see that the reconstructions with the cgls method are noisy ( see fig .
the regularization on graphs ( rg ) is able to reduce noise significantly and sharpen important features .
however , the small search window reconstructions ( method rg - b9 , 9 , 9 ) have a distinctive blocky appearance ( insufficient amount of temporal information collected ) . with a larger search window ( rg - b43 , 43 , 9 ) , reconstructions are much smoother and edges
unfortunately , the computational time is very large for rg - b43 , 43 , 9 ( see tables 2 , 3 ) and therefore the method is impractical .
the computation time is drastically reduced with the proposed method ( arg - bo,9 ) while the error is the lowest .
notably , the quality of reconstruction is very similar to the original method ( rg - b43 , 43 , 9 ) , however some small and thin features are better preserved than with rg - b43 , 43 , 9 . the effect of under and over estimation for various methods can be clearly seen in a 1d plot shown in fig .
5 ( by arrows we show the most erroneous regions in the profile recovered by methods rg - b9 , 9 , 9 and rg - b43 , 43 , 9 ) .
notably , the rg - b43 , 43 , 9 method smooths small features while the proposed method recovers them much better ( see the parabola shaped region in fig . 5 ) . from those numerical experiments
we conclude that the proposed method ( arg - bo ) is more feasible for larger datasets due to performance and quality - wise it is similar or better than the classical nonlocal regularization approaches .
the main purpose of developing such a spatio - temporal regularization technique ( see section section 2 ) is to apply it to time - lapse xmt data . due to high resolution detectors and rapid acquisition
( 
 x ) is very important for iterative minimization of ( 2 ) . in order to perform fast projection - backprojection operations we use gpu - accelerated open - source package : the astra toolbox .
the proposed regularizer is implemented in c , however it can be also realized on a gpu which can provide substantial acceleration .
the following time - lapse experiment was performed on the i13 beamline of dls , where multiple fast 3d scans of time - varying sample ( melting - freezing ice - cream ) were obtained .
the ice - cream microstructure consists of three phases , air cells ( gas ) , ice crystals ( solid ) and unfrozen matrix ( liquid ) ( see fig . 6 , 7 ) .
the aim of the experiment was to understand how phase content varies with temperature variation . during the thermal cycling the inner structure of
a sample is changing ( crystals melt or coarsen ) . since dynamic changes are easily traceable ( controlled sample cooling and heating procedure ) , this dataset is a perfect case study for our reconstruction technique .
additionally , since our reconstruction method searches for some structural correlations in time it significantly enhances the edges between different phases . due to poor contrast between ice - crystals and air bubbles one
has to recover the edges of the unfrozen matrix as accurately as possible . in fig .
6 one can see reconstructions with three different methods : fbp ( filtered back - projection ) , cgls and the proposed arg method .
the time - lapse data size is 2k    2k    1k    6  ( k ) voxels , and therefore 6 large 3d volumes are reconstructed .
we perform 40 iterations of cgls method with one fp iteration ( 11 ) for the arg penalty .
all parameters for the arg method except  and h are selected the same as in table 1 .
the computation time of one arg iteration ( 11 ) to regularize 6 slices ( 2k    2k    6  ( k ) ) is 170 seconds ( all available time frames are used ) .
the classical rg implementation is not used here due to long computation times ( more than 1500 seconds for one iteration ) . from the reconstructions
it can be seen that the direct method ( fbp ) and cgls both deliver unquantifiable reconstructions due to high noise levels and ring artifacts present .
that the arg method is able to suppress noise significantly while also enhance edges between different phases .
in this paper , a generalized spatio - temporal regularization technique for time - lapse tomographic reconstruction is introduced .
the proposed method has been successfully applied to modelled and real dynamic data resulting in noiseless images with enhanced structures ( edges are sharpened while uniform areas are smoothed ) .
the significant benefits of the method are the substantial decrease in reconstruction time ( becomes feasible for large xmt dynamic datasets ) and the possibility of adding more ( potentially all the available ) temporal information into the regularization process .
the future work includes generalization to 4d regularization and application to other important time - lapse experiments .
 x - ray imaging applications in medical and material sciences are frequently limited by the number of tomographic projections collected . 
 the inversion of the limited projection data is an ill - posed problem and needs regularization . 
 traditional spatial regularization is not well adapted to the dynamic nature of time - lapse tomography since it discards the redundancy of the temporal information . in this paper 
 , we propose a novel iterative reconstruction algorithm with a nonlocal regularization term to account for time - evolving datasets . 
 the aim of the proposed nonlocal penalty is to collect the maximum relevant information in the spatial and temporal domains . with the proposed sparsity seeking approach in the temporal space , the computational complexity of the classical nonlocal regularizer 
 is substantially reduced ( at least by one order of magnitude ) . 
 the presented reconstruction method can be directly applied to various big data 4d ( x , y , z+time ) tomographic experiments in many fields . 
 we apply the proposed technique to modelled data and to real dynamic x - ray microtomography ( xmt ) data of high resolution . compared to the classical spatio - temporal nonlocal regularization approach , 
 the proposed method delivers reconstructed images of improved resolution and higher contrast while remaining significantly less computationally demanding . 