Core Imaging Library -- Part II: Multichannel reconstruction for dynamic and spectral tomography
Evangelos Papoutsellis, Evelina Ametova, Claire Delplancke, Gemma Fardell, Jakob S. Jørgensen, Edoardo Pasca, Martin Turner, Ryan Warr, William R. B. Lionheart, Philip J. Withers
CCore Imaging Library – Part II: Multichannel reconstruction fordynamic and spectral tomography
Evangelos Papoutsellis , Evelina Ametova , Claire Delplancke Gemma Fardell Jakob S. Jørgensen , Edoardo Pasca Martin Turner Ryan Warr William R. B. Lionheart Philip J. Withers Henry Royce Institute, Department of Materials, The University of Manchester,Oxford Road, Manchester, M13 9PL, United Kingdom Scientific Computing Department, STFC, UKRI, Rutherford Appleton Laboratory,Harwell Campus, Didcot OX11 0QX, United Kingdom Department of Mathematics, The University of Manchester, Oxford Road, Manchester,M13 9PL, United Kingdom Department of Mathematical Sciences, University of Bath, Claverton Down, Bath,BA2 7AY, United Kingdom Department of Applied Mathematics and Computer Science, Technical University ofDenmark, Richard Petersens Plads, Building 324, 2800 Kgs. Lyngby, Denmark Laboratory for Applications of Synchrotron Radiation, Karlsruhe Institute ofTechnology, Kaiserstr. 12, 76131, Karlsruhe, Germany Research IT Services, The University of Manchester, Oxford Road, Manchester M139PL, United Kingdom ∗ Corresponding author: [email protected]
Abstract
The newly developed Core Imaging Library (CIL) is a flexible plug and play library for to-mographic imaging with a specific focus on iterative reconstruction. CIL provides building blocksfor tailored regularised reconstruction algorithms and explicitly supports multichannel tomographicdata. In the first part of this two-part publication, we introduced the fundamentals of CIL. Thispaper focuses on applications of CIL for multichannel data, e.g., dynamic and spectral. We formalisedifferent optimisation problems for colour processing, dynamic and hyperspectral tomography anddemonstrate CIL’s capabilities for designing state of the art reconstruction methods through casestudies and code snapshots.
Over recent years there has been a growing interest in Dynamic Computed Tomography (DCT), andSpectral X-ray CT thanks to the technological advancements on detector speed and sensitivity and onmultichannel photon-counting detectors (PCDs), as depicted in the EPSRC Tomography roadmap, [1].In DCT the aim is to reconstruct a series of images and depict the complete spatiotemporal responseof the scanned object. These temporal variations may occur because the composition/structure evolved,e.g. corrosion, or the object is subject to external input, e.g. compression, or the object moved duringthe scanning process.In Spectral X-ray CT (SCT), using a pixelated energy sensitive detector, it is possible to collect n energy specific radiographs, where n is the number of energy channels. As a result for any voxel inthe system it is possible to reconstruct the profile of attenuation coefficient as a function of energy, orconversely to create a tomogram corresponding to each energy bin. Since each chemical element has acharacteristic attenuation profile this provides a fingerprint of the elements in each voxel. This fingerprintis especially clear for attenuation spectra that include the energies corresponding to the X-ray absorption1 a r X i v : . [ phy s i c s . m e d - ph ] F e b dges or K-edges for the elements concerned, because there is an abrupt change in attenuation on eitherside of the edge, [2]. Moreover, in pulsed neutron imaging data, sharp edges can also be imaged, [3]. Inthis case the edges, i.e., Bragg edges, correspond to abrupt increases in the transmitted spectrum, whenthe energy is below that possible for Bragg diffraction out of the beam for each diffraction peak providingunique fingerprints corresponding to different crystal structures. In general for spectral imaging, becausethe signal is allocated to a number of energy bins rather than accumulated to give a single image, theenergy resolved data acquired usually suffer from low signal-to- noise ratio, acquisition artifacts andangular undersampling making tomographic image reconstruction difficult.The scope of this paper is to present the capabilities of the Core Imaging Library (CIL) ( ), [4], of the Collaborative Computational Project in Tomographic imaging (CCPi)for multichannel tomography. It allows to reconstruct higher quality images and ensure more accuratespatiospectral K-edge identification, see for instance [5, 6], where novel reconstruction methods are in-troduced for lab-based hyperspectral CT and Neutron Tomography respectively. CIL is an open source,object-oriented library (primarily written in Python) for processing tomographic data. We can read,load, preprocess, reconstruct, visualise and analyse tomographic data from different applications, e.g.,X-ray Computed Tomography, X-ray Laminography, Neutron tomography (NT) and Positron EmissionTomography (PET). Outline of the paper : In the first section, we give a brief overview of CIL and introduce notationfor the optimisation framework, necessary for the reader to make the transition from mathematicalformulation to code. Then, we consider a simple exemplar case study involving two simple 3-channelimaging tasks, i.e., colour denoising and inpainting. For the first task, the aim is to solve the TotalVariation (TV) denoising problem, using the Fast Gradient Projection (FGP) algorithm, [7]. For theinpainting problem, we use the Total Generalised Variation (TGV) regularisation and solve it using thePrimal-Dual Hybrid Gradient algorithm (PDHG), see [8]. In the following sections, we consider two realtomography applications, namely dynamic X-ray and hyperspectral CT. For the dynamic study in section4, we illustrate how one can configure a reconstruction geometry for a tomographic setup in CIL, runstandard non-regularised iterative algorithms and explore different regularisers over the spatiotemporaldomain. In the final example, we deal with hyperspectral 4D tomographic data. We pre-process thedata to eliminate ring artifacts and use a stochastic version of the PDHG, namely the Stochastic PDHGalgorithm, [9], under TV regularisation to reconstruct the data with different coupling between spatialand spectral dimensions.
In part I, [4], we described the main building blocks of CIL namely: cil.framework , cil.optimisation , cil.processors , cil.io and cil.utilities . We illustrated the basic usage of CIL data structures, as ap-plied to a number of X-ray CT cases with different geometries, e.g., parallel, cone and laminographyand also different modalities such as NT and PET. CIL wraps a number of third-party libraries, us-ing cil.plugins , to perform various operations required for CT reconstruction. For instance, we canuse the Astra-Toolbox, [10], to perform forward and back projection steps and filtered back projection(FBP) reconstruction and can use the CCPi-Regularisation Toolkit to employ several regularisers witha CPU/GPU hardware acceleration. In addition, CIL is designed such that the data structures of theSynergistic Image Reconstruction Framework ( SIRF ), [11], from the Collaborative Computational Projectin Synergistic Reconstruction for Biomedical Imaging (CCP-SynerBi), , can beused for PET and Magnetic Resonance Imaging (MRI) reconstruction.
The cil.optimisation framework contains three structures, namely
Function , Operator and
Algorithm that formalise a generic optimisation problem for imaging applications as u ∗ = arg min u ∈ X f ( Ku ) + g ( x ) ≡ arg min x ∈ X n (cid:88) i =1 f i ( K i ( x )) + g ( x ) . (1)We let X , Y two finite dimensional vector spaces, K : X → Y a linear operator with operator norm (cid:107) K (cid:107) = max {(cid:107) Kx (cid:107) Y : (cid:107) x (cid:107) X ≤ } and proper, convex functions f : Y → R , g : X → R . Note that in R := R ∪ {∞} Y = Y × . . . Y n , n ≥ f ( y ) := f ( y , . . . , y n ) = (cid:80) ni =1 f i ( y i ) which results to the right-side formulation in (1).In the following case studies, using different definitions for the triplet ( K , f , g ), we can expressoptimisation problems for several imaging tasks. For example, in denoising, we let K be the identityoperator, in inpainting it is a mask operator that encodes missing pixel information, while it is a projectionoperator for tomography. The functions f , g allow us to define a fidelity term, that measures the distancebetween the acquired data b and the forward-projected reconstruction image as well as well as a regulariser,which enforce a certain regularity on x . If the noise follows a Gaussian distribution an appropriate choicefor the fidelity term is (cid:107) Kx − b (cid:107) . In the case of impulse noise, the L norm (cid:107) Kx − b (cid:107) leads to moreefficient restorations and for Poisson noise the Kullback-Leibler divergence (cid:82) Kx − b log Kx is the mostsuitable choice. The choice of the regulariser, e.g., Tikhonov, TV and TGV, favours minimisers of (1)with certain geometric features and is usually weighted by positive parameters to control the influencebetween data fidelity and regularisation terms.In order to find an approximate solution for minimisation problems of the (1) form, we use differentCIL Algorithm for smooth and non-smooth objective functions such as the Conjugate Gradient LeastSquares (CGLS), Simultaneous Iterative Reconstruction Technique (SIRT) and proximal type algorithms,which are extensively used in this papers, such as the FGP and the PDHG, SPDHG algorithms. In theFGP algorithm, we require that the function g has a proximal operator defined asprox τg ( u ) := arg min v (cid:107) v − u (cid:107) + τ g ( v ) , (2)which has a ”simple” closed form solution or it can be computed efficiently numerically. Also, we assumethat f is continuously differentiable and has Lipschitz continuous gradient L . On the other hand, in thePDHG algorithm, we allow functions f and g to be non-differentiable and express (1) into a saddle pointproblem , min u ∈ X max z ∈ Y (cid:104) Ku, z (cid:105) − f ∗ ( z ) + g ( u ) , (3)where f ∗ denotes the convex conjugate of f . Under this setup, PDHG can decouple the initial problem(1) into two simple problems, using as before the proximal operators of g and in addition the proximaloperator of f ∗ , prox τf ∗ ( u ) := arg min v (cid:107) v − u (cid:107) + τ f ∗ ( v ) . (4) We begin our first demonstration, with a case study in a colour imaging framework i.e., a vector-valuedimage that has just 3 channels: Red, Green and Blue. Our test data is a high resolution double rainbow image taken from a smartphone of 1194 × { ( i, j ) | ≤ i < M, ≤ j < N, M = 1194 , N = 1353 } be a rectangular discretised grid representing our image domain and define an RGB colour image u as u : Ω → R , u = ( u , u , u ) , where u k ∈ R M × N , k = 1 , , We start with one of the most fundamental and well-studied problems in image processing, that is imagedenoising. In the pioneering work, Rudin, Osher and Fatemi (ROF), [12], introduced the TV regulariserto tackle image denoising for grayscale images u : Ω → R . Given a noisy image b corrupted by additiveGaussian noise, they solve the following optimisation problem u ∗ = arg min u (cid:107) b − u (cid:107) + α TV( u ) , (5) Image can be downloaded from https://github.com/vais-ral/CIL-data (cid:96) , norm of the gradient. If the gradientis Du := ( D y u, D x u ), where D y , D x denote the finite differences along the y and x directions respectively,we write TV( u ) = (cid:107) Du (cid:107) , = M,N (cid:88) i,j (cid:0) | ( D y u, D x u ) | (cid:1) i,j = M,N (cid:88) i,j (cid:0)(cid:113) ( D y u ) + ( D x u ) (cid:1) i,j . (6)The above definition can be extended for vector-valued images, which results to a vectorial version of thetotal variation. The gradient for the RGB case is now Du = ( Du , Du , Du ), where for each k = 1 , , Du k := ( D y u k , D x u k ). The Vectorial Total Variation (VTV), see [13], is defined asVTV( u ) := (cid:107) Du (cid:107) , = M,N (cid:88) i,j (cid:0) | ( Du , Du , Du ) | (cid:1) i,j = M,N (cid:88) i,j (cid:18)(cid:113) ( | Du | ) i,j + | Du | + | Du | (cid:19) i,j = M,N (cid:88) i,j (cid:18)(cid:118)(cid:117)(cid:117)(cid:116) (cid:88) k =1 ( D y u k ) + ( D x u k ) (cid:19) i,j (7)For both grayscale and coloured images, we can setup the regularisers (6) and (7) using the TotalVariation function in CIL. The optimisation problem that we solve for the colour denoising is similar to (5) butusing the VTV regulariser, i.e., u ∗ = arg min u (cid:107) b − u (cid:107) + α VTV( u ) , (8)where b is shown in Figure (1b). One can observe that (8) is in fact the proximal operator (2) with τ = 1 . b . We solve (8), using the Fast Gradient Projection (FGP) algorithm that is contained inthe proximal method of the TotalVariation function.
Proximal TV Denoising
VTV = 0.15*TotalVariation(max_iteration=2000)proximalVTV = VTV.proximal(noisy_data, tau=1.0)
As it is observed noise is eliminated, while preserving the edges of the image. However, total variation isknown for promoting piecewise constant reconstructions leading to images with blocky structures. Thisis called the staircasing effect and becomes apparent in smooth regions, see for instance the area aroundthe rainbow in Figure (1c).
Given an image where a specific region is unknown, the task of image inpainting is to recover the missingregion from the known part of the image. A popular application of inpainting is in art restoration, wheredamaged or missing areas are repainted, i.e., filled, based on the surrounding context. We let
D ⊂
Ω bea subdomain of Ω, i.e., the inpainting domain , where no data are known and missing information shouldbe interpolated. In this example, our input image is shown in Figure (1d), where in addition to Salt &Pepper noise, missing pixels from a repeated text are incorporated. A suitable data fidelity term for thiskind of noise distribution is the L norm that acts on the domain Ω \ D .To overcome the staircasing artifacts that TV promotes, we employ a higher-order regulariser namelythe Total Generalised Variation introduced in [14]. We let α, β > α,β ( u ) = min w α (cid:107) Du − w (cid:107) , + β (cid:107)E w (cid:107) , , (9)where E denotes the symmetrised gradient operator defined as E w = ( Dw + Dw T ). The optimisationproblem above provides a way of balancing between the first and second derivative of an image u . Inparticular, one expects that in the neighbourhood of edges, the second derivative D u is relatively ”large”,hence a reasonable choice is to let w = 0 in (9) and recover the total variation regulariser. On the otherhand, D u is relatively small in smooth regions of an image and w = Du is a proper condition forthe minimisation problem (9). Under this format, edges are preserved, as in the TV regulariser, andadditionally piecewise smooth structures are promoted.4 a) Ground Truth (b) noisy_data (Gaussian Noise) (c)
VTV denoising( α = 0 . (d) noisy_data (Salt and Pepper Noisewith missing pixels) (e) TGV inpainting( α, β ) = (0 . , . (f ) Absolute difference—(e) - (a)—
Figure 1:
Colour Processing: 1st row) Total Variation denoising. 2nd row) Total Generalised Variationinpainting and denoising. Regularising parameters are optimised based on the Structural Similarity IndexMeasure.The minimisation problem under the TGV regulariser and the L norm as a data fidelity term is thefollowing: u ∗ = arg min u (cid:107)M u − b (cid:107) + TGV α,β ( u ) ⇔ ( u ∗ , w ∗ ) = arg min u,w (cid:107)M u − b (cid:107) + α (cid:107) Du − w (cid:107) , + β (cid:107)E w (cid:107) , , (10)where the M is a diagonal operator with 1 in diagonal elements corresponding to pixels in Ω \ D and 0in D . In CIL, we use the MaskOperator that accepts as an input a 2D boolean array, i.e., mask . Sincewe have a colour image, we employ the
ChannelwiseOperator to encode the effect of missing pixels tothe RGB channels. In order to solve (10), we use the PDHG algorithm, where the first step is to express(10) in the general form of (1). Let u = ( u, w ) ∈ X and define an operator K : X → Y as K = M O D −IO E ⇒ K u = K (cid:20) uw (cid:21) = M uDu − w E w = y y y = y ∈ Y , (11)where O , I denote the zero and identity operators respectively. We continue with the definition of thefunctions f and g . The function f is a separable function that contains the three terms in (10) and isdefined as f ( y ) := f ( y , y , y ) = f ( y ) + f ( y ) + f ( y ) , where ,f ( y ) := (cid:107) b − y (cid:107) , f ( y ) := α (cid:107) y (cid:107) , , f ( y ) := β (cid:107) y (cid:107) , , (12)and g ( u ) = g ( u, w ) = O ( u ) ≡ f ( K u ) + g ( u ) = f (cid:18) M uDu − w E w (cid:19) = f ( M u ) + f ( Du − w ) + f ( E w )= (cid:107)M u − b (cid:107) + α (cid:107) Du − w (cid:107) , + β (cid:107)E w (cid:107) , , BlockOperator K and it is filled row-wise. Theelements are the GradientOperator D , the IdentityOperator I , the SymmetrisedGradientOperator E , the ChannelwiseOperator M and two ZeroOperator O . The separable function in (12) can beexpressed by the BlockFunction f , whose elements are the L1Norm , and two
MixedL21Norm functions.Finally, g is the ZeroFunction function. We choose the
PDHG algorithm to solve such an optimisationproblem. Without any user input, CIL will by default use primal/dual stepsizes σ, τ with σ = 1 . τ = . σ (cid:107) K (cid:107) that satisfy στ (cid:107) K (cid:107) < verbose=2 . However, to speed up convergence it may be needed to change such default values. PDHG: TGV-L Inpainting
K11 = ChannelwiseOperator(MaskOperator(mask), 3, dimension='append')K21 = GradientOperator(ig)K32 = SymmetrisedGradientOperator(K21.range)K12 = ZeroOperator(K32.domain, ig)K22 = -IdentityOperator(K21.range)K31 = ZeroOperator(ig, K32.range)K = BlockOperator(K11, K12, K21, K22, K31, K32, shape=(3,2))f1 = L1Norm(b=noisy_data)f2 = 0.5*MixedL21Norm()f3 = 0.2*MixedL21Norm()f = BlockFunction(f1, f2, f3)g = ZeroFunction()pdhg = PDHG(f=f, g=g, operator=K, max_iteration=2000,update_objective_interval=500)pdhg.run(verbose=2)
The TGV reconstruction is presented in Figure (1e), where there are no staircasing issues and mostof the repeated text is eliminated. We observe that the inpainting process behaves quite well when thebackground is relatively smooth, e.g., sky. However, in regions with specific textures, such as trees, leavesand grass, TGV inpainting could not completely restore the missing pixels, see the absolute differencebetween the ground truth and reconstruction in Figure (1f).In the above two examples, the optimal regularising parameters are chosen as maximising the Struc-tural Similarity Index Measure, [15], since the ground truth image is available. However, in the followingstudies, we deal with real tomographic data and the ”best” reconstructions among different parametersare based on visual inspection.
Our next case study concerns a Dynamic X-ray Computed Tomography application. In DCT, given animage sequence of acquired data, e.g., sinograms, the main goal is to reconstruct the temporal changesin the distribution of the linear attenuation coefficient. The aim of this section is to demonstrate to thereader the capabilities of CIL in a spatiotemporal tomographic setup. We show how one can specify therequired geometries and the corresponding projection operator. In addition, we compare reconstructionmethods using different regularisation setups applied to spatiotemporal data. In the following we employan open-access dataset available from [16].
Description : The sample is an agarose-gel phantom perfused with a liquid contrast agent in a 50 mlFalcon test tube, [17]. The aim of this experiment is to simulate diffusion of liquids inside plant stems,namely the flow of iodine-based contrast agents used in high resolution tomographic X-ray imaging ofplants. After the agarose solidified, five cavities were made into the gel and filled with 20% sucrosesolution to guarantee the diffusion by directing osmosis to the gel body.
Acquisition : Each scan was acquired in 4.5 minutes with intermissions of approximately 15 minutesbetween consecutive measurements. In total, the acquisition process lasted about 3 hours leading to67 sinograms one for each time state. The acquired sinograms are pre-processed using Lambert-Beernegative logarithm conversion.
Dataset : Every measurement consists of 360 projections obtained from a cone beam microCT-scanner. In this case only the central slice is used, resulting in
2D fanbeam geometry . Depending onthe detector binning, we can choose between a 512x512 or 256x256 reconstruction volume. For thisexperiment, our reconstruction volume is 256x256 of 17 time frames. Additional metadata informationare provided to setup the cone beam geometry.In order to begin our iterative reconstructions, we first need to setup our geometries, i.e., the
ImageGeometry , the
AcquisitionGeometry and the
ProjectionOperator A acting over all channels.Note that our image domain is a 3D spatiotemporal volume, i.e.,Ω = { ( i, j, t ) : 0 ≤ i < M, ≤ j < N, ≤ t < T, M = N = 256 , T = 17 } . The following code defines the geometries and the projection operator:
Geometries + Operator ag = AcquisitionGeometry.create_Cone2D(source_position=[0, -distanceSourceOrigin],detector_position=[0, distanceOriginDetector]).set_panel(num_pixels=282, pixel_size=0.2).set_channels(num_channels=17).set_angles(angles, angle_unit='radian').set_labels(['channel', 'angle', 'horizontal'])ig = ag.get_ImageGeometry()A = ProjectionOperator(ig, ag, 'gpu')
In [4], we presented in detail how one can setup the CGLS algorithm and Tikhonov regularisation forsingle channel X-ray tomography. In this section, we apply the same regularisation methods but fordynamic data. We solve the following problems: u ∗ = arg min u (cid:107) Au − b (cid:107) , (Least Squares) (13) u ∗ = arg min u (cid:107) Au − b (cid:107) + α (cid:107) Lu (cid:107) , (Tikhonov) , (14)where, b is the dynamic sinogram ( data ) of 17 time-frames (channels), see Figure 2. The second termin (14) acts as a smooth regulariser, where the linear operator L , can be either an identity or a gradientoperator. In the case of L = D , we offer the user two different modes for the GradientOperator , wherefinite differences are computed only along the spatial dimensions or for both the spatial and channel di-mensions. Therefore, if
L = GradientOperator(ig,correlation) , we take into account the derivativesacross every direction in a 3D volume, i.e., Du = ( D t u, D y u, D x u ), if correlation= ' SpaceChannels ' ,whereas if correlation= ' Space ' ) , we take into account only the derivatives across the spatial dimen-sions, i.e., Du = ( D y u, D x u ). For (14), we use correlation= ' SpaceChannels ' . The code snippet inorder to setup the above problems in CIL is identical to the one presented in [4], hence it is omitted.For CGLS reconstruction that is performed channelwise, noise is still present and can be eliminatedwith intense blurring artifacts spatially using Tikhonov regularisation, see Figure 3. Additionally, withthe spatiotemporal coupling enforced in (14), noise reduction is observed along the temporal direction, asdepicted in the dynamic profiles in Figure 4. However, one needs to consider applying different regularisersalong the temporal and spatial dimensions to obtain a better quality reconstruction. regularisation We continue our reconstructions with a more advanced choice of regularisers and explain how the usercan penalise differently the spatial and temporal domains. Motivated by the section 3, we extend thetotal variation regulariser in a spatiotemporal setting and for Du = ( D t u, D y u, D x u ), we write,TV( u ) = (cid:107) Du (cid:107) , = M,N,T (cid:88) i,j,t (cid:0)(cid:113) ( D t u ) + ( D y u ) + ( D x u ) (cid:1) i,j,t . (15)7 igure 2: Sinogram data of the gel-phantom (17 time-frames in total) for 3 different time-frames.
Figure 3:
Tomographic reconstructions of the gel phantom, using the minimisation problems (13), (14),solved by the
CGLS algorithm and (16), (17), solved by the
PDHG algorithm. Only 3 time-frames from atotal of 17 are shown. Regularising parameters are α = 0 . α = 0 .
001 (Spatiotemporal TV)and ( α, β )=(4 · − , 50) (TV - L ). 8 igure 4: Temporal variations in the attenuation of the gel-phantom for specific ROIs (blue, green)shown in the top left reconstruction of Figure 3.Under this isotropic coupling between space and time, the finite differences along the directions t, y and x are penalised equally with a single regularising parameter, promoting piecewise constant structures inthe spatiotemporal volume by solving u ∗ = arg min u> (cid:107) Au − b (cid:107) + α TV( u ) , (Spatiotemporal TV) . (16)However, it is natural to assume that the temporal domain does not contain edges and a smootherregulariser would be a more appropriate choice. Therefore, we decompose the spatiotemporal volumeand at the same time, apply different regularisation terms, i.e, total variation penalty in the spatialdomain and a (cid:107) · (cid:107) norm for the temporal domain. We solve u ∗ = arg min u> (cid:107) Au − b (cid:107) + α (cid:107) ( D y u, D x u ) (cid:107) , + β (cid:107) D t u (cid:107) , (TV-L ) . (17)In order to demonstrate the flexibility of the CIL, we choose to solve the optimisation problems (16), (17)using the PDHG algorithm but with different formations of the triplet ( K , f , g ). For (16), we setup an implicit PDHG algorithm, where one of the proximal operator does not have a closed form solution butrather an inner iterative solver is used, see [18]. We let f be the squared L norm for the data fidelityterm and the operator K for the projection operator A . As in section 4, we let g be the TotalVariation function with a nonnegative constraint using lower=0.0 . (Implicit) PDHG: Spatio-Temporal TV K = Af = 0.5*L2NormSquared(b=data)g = 0.001*TotalVariation(max_iteration=100, lower=0.0,correlation='SpaceChannels')
For (17), we use an explicit
PDHG setup, where the proximal operator of f ∗ and g have closed formsolutions, see [19]. We consider the FiniteDifferenceOperator along the temporal, y and x directionsand use the (cid:107)·(cid:107) , norm for the spatial gradient Du = ( D y u, D x u ) and the (cid:107)·(cid:107) for the temporal gradient D t . Then, we define a BlockOperator K to wrap all the three gradients, as well as the projection9perator. In the same manner, we use the BlockFunction f to wrap all three terms appeared in (17).Finally, to obtain nonnegative reconstructions, we use the IndicatorBox for the function g . A comparisonbetween implicit and explicit cases of PDHG is beyond the scope of this paper and we refer the reader to[20] for a detailed review. For both implicit and explicit cases the setup of the PDHG algorithm is identicalto the one used in section 3. (Explicit) PDHG: TV - L Dy = FiniteDifferenceOperator(ig, direction='horizontal_y')Dx = FiniteDifferenceOperator(ig, direction='horizontal_x')spatialGradient = BlockOperator(Dy, Dt, shape=(2,1))Dt = FiniteDifferenceOperator(ig, direction='channel')K = BlockOperator(A, spatialGradient, Dt, shape=(3,1))f1 = 0.5*L2NormSquared(b=data)f2 = 4e-4*MixedL21Norm()f3 = 50/2*L2NormSquared()f = BlockFunction(f1, f2, f3)g = IndicatorBox(lower=0.0)
In the last two rows in Figure 3, we present the Spatiotemporal TV and TV-L reconstructionsrespectively. Compared to the Tikhonov regularisation, spatiotemporal TV can eliminate the noise whilepreserving the edges in the spatial domain. In addition, staircasing artifacts are also apparent especiallyaround the five cavities in the gel, see for instance the zoomed regions in Figure 3. However, there is nosignificant difference between these two reconstructions with respect to the temporal activity, see purpleand blue curves in Figure 4. This is because, temporal variations are relatively small compared to spatialones, hence, they are not penalised when an isotropic coupling is enabled.One the other hand, using the TV-L regulariser, we are able to reduce the staircasing around thefive cavities as well as in the center of the tube. This is also obvious from the temporal variations forthe TV-L (red curve) of the green ROI in Figure 4. Therefore, with the proposed decoupling undera spatial edge preserving prior and a smooth temporal penalty, a better quality reconstruction can beobtained which reflects the physical characteristics of this dynamic tomography experiment. For our final case study, we demonstrate how CIL can provide several tools to pre-process, reconstructand analyse hyperspectral data. In hyperspectral X-ray computed tomography, we have the ability todetect the position and energy of each transmitted X-ray photon in a single acquisition scan. Thiswill allow us to identify materials in the sample based on the element characteristic absorption edges,i.e., K-edges. Here, we use a mineralised ore sample from a gold-rich hydrothermal vein and identifyelements of gold and lead. We first pre-process the data to correct from ring artifacts and then follow asimilar regularisation strategy as in section 4 in terms of reconstruction. This time, we compare differentformulations on the total variation regulariser and apply the SPDHG algorithm to reconstruct efficientlylarge acquisition data.
Description:
The sample is extracted from a hydrothermal vein from the Leopard Mine, Silobela,Zimbabwe. It contains materials such as pyrite (FeS ), quartz (SiO ), gold (Au) and minor amounts ofgalena (PbS), chalcopyrite (CuFeS ) and bornite (Cu FeS ). Acquisition: ◦ . The total scan time was 7.5 hrs. The acquired 4Draw sinogram data has three spatial dimension and one spectral dimension, i.e., 120 projections anglesof 80 ×
400 pixels with 800 energy bins ranging from 1.82 keV to 186.07 keV.10 ataset:
All the files of this study are freely available and can be downloaded from [22]. It containsa) 4D hyperspectral (energy-resolved) X-ray CT projection data, b) flat-field data, c) energy in keV forevery energy bin and d) the cone-geometry setup.We load the initial raw data which contain all 800 energy channels recorded during the acquisitionprocess. It is possible to consider the full energy range, however for demonstration purposes we chooseto examine an energy interval of 80 channels that encompasses the K-edges of gold (Au, 80,725 keV) andlead (Pb, 88.005 keV).
In Figure 5, one can observe that the sinogram is heavily affected by vertical stripes. Typically, someelements in the flat panel detectors behave abnormally due to a wrong calibration or beam instability.These stripes result to the well-known ring artifact phenomenon. We eliminate these artifacts usinga
RingRemover processor. It is wavelet-based algorithm, [23], acting in the sinogram domain and canhandle both 2D/3D
ImageData and
AcquisitionData as well as multichannel data. In the latter case,the algorithm is applied separately for each channel. It accepts the level of decomposition L , the wavelettype and the damping factor σ . In Figure 5, we compare the raw sinogram and the sinogram afterthe RingRemover processor. In the code below, we present the
RingRemover processor applied to the cropped_sinogram of 80 channels. The cropped_sinogram is obtained using the
Slicer processor andselecting the [318 , RingRemover cropped_sinogram = Slicer(roi={'channel': (318, 398)})(raw_data)b = RingRemover(decNum=4, wname='db25', sigma=1.5)(cropped_sinogram)
Figure 5:
Slice through the (4D) sinogram for the hyperspectral ore sample dataset, corresponding to84.37 keV, for the 40th vertical slice, before and after the RingRemover processor.
For our first reconstruction method, we use the
SIRT algorithm for the raw data before and after the
RingRemover processor. The SIRT algorithm is a particular case of a weighted least-squares problemand in addition, it accepts certain convex constraints such as a nonnegative constraint. We enforce thisconstraint using lower=0.0 . The
SIRT algorithm is applied channelwise for each 3D dataset. As shownin the code below, we first need to extract the acquisition and image geometries, i.e., ig3D and ag3D respectively for the 3D dataset and define our projection operator. We also need to define an imagegeometry for the final 4D SIRT reconstruction sirt4D_reconstruction . For every channel, we run100 iterations and its solution sirt3D_reconstruction is filled to the corresponding channel i of the sirt4D_reconstruction . Note that for every channel the SIRT algorithm is initialised with the solutionof the previous sirt3D_reconstruction . We compare the two reconstructions in Figure 6.11 IRT reconstruction ag3D = b.geometry.subset(channel=0)ig3D = ag3D.get_ImageGeometry()A = ProjectionOperator(ig3D, ag3D)ig = b.geometry.get_ImageGeometry()sirt4D_reconstruction = ig.allocate()x0 = ig3D.allocate()constraint = IndicatorBox(lower=0)sirt3D_reconstruction = SIRT(max_iteration = 100)for i in range(ig.channels):sirt3D_reconstruction.set_up(initial=x0, operator=A,constraint=constraint,data=data.subset(channel=i))sirt3D_reconstruction.run(verbose=0)sirt4D_reconstruction.fill(sirt3D_reconstruction.solution, channel=i)x0.fill(sirt3D_reconstruction.solution)
Figure 6:
SIRT reconstructions in orthogonal slices for the sinograms in Figure 5.
For the dynamic CT application, in section 4, we used two different formulations for the total variationregularisation, namely spatiotemporal TV and TV-L . Here, we follow the same strategy but this timewe allow piecewise constant structures both in the spatial and spectral dimensions. This enables us topreserve edges in the spatial domain and absorption K-edges in the energy spectrum. We first considerthe total variation regulariser extended to a 4D volume, where the gradient Du = ( D e u, D z u, D y u, D x u )is coupled isotropically. We solve u ∗ = arg min u> (cid:107) Au − b (cid:107) + α (cid:88) (cid:113) ( D e u ) + ( D z u ) + ( D y u ) + ( D x u ) . (Spatiospectral TV)(18)This regulariser combines together the spatial and energy variations which are penalised by a singleregularising parameter. Similarly to (15), if for example the finite difference along the energy is negligiblecompared to the spatial gradients, the regularisation is not active for this direction. An alternative setup,12hat produces better reconstructions results, is to enforce separate regularisation with respect to theenergy and spatial gradients, i.e., D e u and ( D z u, D y u, D x u ) respectively. Then, we solve u ∗ = arg min u> (cid:107) Au − b (cid:107) + α (cid:107) ( D z u, D y u, D x ) (cid:107) , + β (cid:107) D e u (cid:107) . (3D + spectral) TV (19)Although, we can follow exactly the same non-smooth setup presented in the previous section to solvethe above problems, one has to perform forward and backward operations of the projection operator A for all the sinogram bins per iteration. These operations are computationally expensive, especially forsuch large dataset. In order to overcome this problem, we employ the SPDHG algorithm, where theabove operations are applied to a randomly selected subset of the data in every iteration. It has beenused to different clinical imaging applications, such as PET, [24] and motion estimation/correction inPET/MR, [25], and produce significant computational improvements compared to the PDHG algorithm.The setup of the SPDHG algorithm is very similar to the PDHG algorithm. The only difference is thatwe need to specify the number of subsets, as well as a list of probabilities for each subset, see Algorithm1. For our hyperspectral reconstruction, we split the 120 projections of our sinogram data g = ( g i ) ni =1 ,into n = 10 subsets of 12 angles each. We can configure different sampling patterns, namely uniform and balanced sampling. The first choice assigns an equal probability for every subset n , i.e., p i = n . In thelatter setting, we choose a uniform probability for the data fidelity term, i.e., p i = m , 1 ≤ i ≤ m = n − p n = for the regulariser. In the following, we choose the balanced sampling approach and refer thereader to [24] for a comparison between these two sampling patterns. Algorithm 1
SPDHG
Inputs: n = γ = 1 . ρ = 0 .
99, probability p i , i = 1 , . . . , n Stepsizes: S i = γ ρ (cid:107) A i (cid:107) , T = min i γ ρp i (cid:107) A i (cid:107) , i = 1 , . . . , n Initialize: x , z = z , y Update: x k , y k , z k for k ≥ do x k +1 = prox TG ( x k − T z k ) Select i ∈ [ n ] at random with probability p i y k +1 i = prox S i F ∗ i ( y ki + S i A i x k +1 ) ∆ z = A Ti ( y k +1 i − y ki ) z k +1 = z k + ∆ z z k +1 = z k +1 + ∆ zp i end for For each subset ( g i ) ni =1 , 1 ≤ i ≤ n , we extract the corresponding acquisition geometries and create alist called list_geometries . Next, using a list of ProjectionOperator A = ( A i ) ni =1 that share the same ImageGeometry , we wrap them as the
BlockOperator K . Finally, the GradientOperator is appended.
Subsetting operators
A_i = [ProjectionOperator(ig, ageom) for ageom in list_geometries]A_i.append(GradientOperator(ig)K = BlockOperator(*A_i)
Similarly, we create a list of data fidelity terms, fsubsets , using the
L2NormSquared for each of theacquisition data ( g i ) ni =1 and wrap them as the BlockFunction f . For the spatiospectral TV, we use BlockFunction to wrap the list of data fidelity terms and the standard (cid:107) · (cid:107) , norm. However, CIL doesnot provide a build-in Function for the (3D+spectral) TV regulariser. In order to implement a splittingas in (19), we use an additional
BlockFunction that contains the last two terms, namely the (cid:107) · (cid:107) normof the finite difference along the spectral direction and (cid:107) · (cid:107) , for the spatial gradients, weighted by β and α regularising parameters respectively. 13 patiospectral TV vs (3D+spectral) TV subsets = 10fsubsets = [0.5*L2NormSquared(b=g[i]) for i in range(subsets)]f = BlockFunction(*Fsubsets, alpha*MixedL21Norm())splitTV = BlockFunction(alpha*MixedL21Norm(), beta*L1Norm())f = BlockFunction(*Fsubsets, splitTV)g = IndicatorBox(lower=0.0) Finally, the code to setup the
SPDHG algorithm for both (18), (19) minimisation problems is:
Setup SPDHG prob = [1/(2*subsets)]*subsets + [1/2]spdhg = SPDHG(f=f, g=g, operator=K, max_iteration=500, prob=prob)spdhg.run()
Figure 7:
SIRT (after
RingRemover ) and Spatiospectral TV, (3D+spectral) TV reconstructions usingthe SPDHG algorithm for energies: (left) below the gold K-edge, (middle) after the gold K-edge butbefore the lead K-edge and (right) above the lead K-edge. Regularising parameters are α = 0 .
002 for theSpatiospectral TV and ( α, β )=(0.001, 0.2) for the (3D+spectral) TV.14e compare the SIRT, Spatiospectral TV and (3D+spectral) TV reconstructions for three differentenergy channels, see Figure 7. We observe that for both SIRT and Spatiospectral TV reconstructions,suffer from significant noise artifacts spatially. In fact, it is difficult to distinguish between noisy voxelsand the signal from the gold and lead in the sample. Also, same effect is observed in Figure 8 for the twospectra for ROIs centered on the gold (ROI ) and lead (ROI ) containing materials. The first spectralplot shows a step change in attenuation corresponding to the Au K-edge (80.725 keV) whilst the secondspectrum shows a Pb K-edge (88.005 keV). The spectral noise can be reduced further by choosing ahigher regularising parameter for the Spatiospectral TV, but then small features will be lost spatiallydue to loss of contrast. Therefore, using the (3D+energy) TV regulariser, we can find an optimal pair ofparameters and balance the strength between TV regularisation in space and spectral directions. In thisway, we obtain higher quality reconstructions with better contrast and less noise. Figure 8:
Energy plots for SIRT, Spatiospectral TV and (3D+spectral) TV reconstructions for selectedROIs for the gold (ROI ) and lead (ROI ) shown in Figure 7. Multichannel CT imaging opens up many new possibilities in material and life sciences. Multichannel CTis intrinsically ’photon-hungry’ because detected photons are shared between multiple energy or time bins.Therefore acquired tomographic datasets typically do not provide sufficient information for high qualityreconstruction using traditional FBP-type algorithms. The absence of efficient reconstruction routinescapable of handling noisy and/or undersampled multichannel data hamper scientific applications of thetechnique.The inverse problem framework provides methods to treat these challenging multichannel CT datathrough iterative reconstruction with suitable regularisation which efficiently exploits prior knowledgeand inter-channel correlation. CIL implements essential building blocks, which can explicitly supportmultichannel reconstruction and four dimensional datasets. Here, we have demonstrated the potentialof CIL for multichannel CT data with three representative case-studies. Starting with a simple colourdenoising and inpainting problem, we illustrated the ability to incorporate various regularisation tech-niques, such as classical TV, vectorial TV and TGV. We also outlined how a conventional formulation ofiterative reconstruction through the optimisation framework is mapped onto CIL objects. In the secondcase study, we exploited spatiotemporal TV and a combination of spatial TV with smoothing over thetemporal direction. The later regularisation allowed us to balance the noise level in spatial and temporal15irection and to enhance reconstruction quality compared to both spatiotemporal TV and other con-ventional reconstruction techniques. To highlight the flexibility of CIL we constructed both explicit andimplicit PDHG to solve the reconstruction problem. Finally, in the last case study we endeavoured toreconstruct an energy-resolved X-ray CT dataset with high energy resolution. The dataset suffered fromheavy ring artifacts. Hence we used CIL preprocessing tools to reduce ring artifacts. We followed thesame regularisation strategy as in the second case study, i.e. a combination of spatial TV with smoothingover spectral direction, but this time we used a stochastic version of PDHG algorithm to speed-up largescale CT reconstruction. Regularisation aided better identification of K-edges in the energy-resolvedX-ray CT dataset.The ability to incorporate and balance various regularisation terms in the reconstruction routine isa promising approach to treat noisy and undersampled multichannel CT data. It is widely understoodthat one of the main challenges of iterative reconstruction is the high computational cost compared totraditional FBP-type methods. To resolve this CIL wraps hardware accelerated libraries to performcostly forward- and back-projection steps and calculate proximal of regularisation term. However thereis a continuous effort to improve performance and resolve computational bottlenecks. In terms of sup-ported imaging modalities, we can currently handle any tomographic modality which can be describedby the Beer-Lambert law and provide interoperability with Synergistic Image Reconstruction Frameworkenabling positron emission tomography and magnetic resonance imaging reconstruction using CIL. Wecontinue to enrich the library of available algorithms, regularisers, pre- and post-processing tools alongwith supported imaging models and available back-ends.
Data accessibility
Python scripts to reproduce the results for all the case studies are provided at https://github.com/vais-ral/CIL-Demos/tree/papers/Papers/CIL2/ . Author Contribution
EPap carried out all the experiments, wrote the manuscript and developed the CIL software. JJ andEPas conceived the CIL software, assisted with the experiments and developed the CIL software. EAand GF assisted with the experiments and developed the CIL software. CD co-developed the SPDHGalgorithm. RW assisted with hyperspectral case study and contributed to the CIL software. MT, WL andPW helped conceptualise and roll out the CIL software. All authors critically revised the manuscript,gave final approval for publication and agree to be held accountable for the work performed therein.
Competing interests
The author(s) declare that they have no competing interests.
Funding
The work presented here was funded by the EPSRC grants “A Reconstruction Toolkit for Multichan-nel CT” (EP/P02226X/1), “CCPi: Collaborative Computational Project in Tomographic Imaging”(EP/M022498/1 and EP/T026677/1). We acknowledge the EPSRC for funding the Henry MoseleyX-ray Imaging Facility through grants (EP/F007906/1, EP/F001452/1, EP/I02249X/1, EP/M010619/1,EP/F028431/1, and EP/M022498/1) which is part of the Henry Royce Institute for Advanced Materialsfunded by EP/R00661X/1. JSJ was partially supported by The Villum Foundation (grant no. 25893).EA was partially funded by the Federal Ministry of Education and Research (BMBF) and the Baden-W¨urttemberg Ministry of Science as part of the Excellence Strategy of the German Federal and StateGovernments. WRBL acknowledges support from a Royal Society Wolfson Research Merit Award. PJWacknowledges support from the European Research Council grant No. 695638 CORREL-CT. CD ac-knowledges support from the EPSRC grant EP/S026045/1 “PET++: Improving Localisation, Diagnosisand Quantification in Clinical and Medical PET Imaging with Randomised Optimisation”.16 cknowledgements
This work made use of computational support by CoSeC, the Computational Science Centre for ResearchCommunities, through CCPi.
References [1] EPSRC X-Ray Tomography Roadmap 2018;. https://epsrc.ukri.org/files/research/epsrc-x-ray-tomography-roadmap-2018/ .[2] Kruger RA, Riederer SJ, Mistretta CA. Relative properties of tomography, K -edge imaging, andK -edge tomography. Medical Physics. 1977 May;4(3):244–249. Available from: https://doi.org/10.1118/1.594374 .[3] Santisteban JR, Edwards L, Fitzpatrick ME, Steuwer A, Withers PJ, Daymond MR, et al. Strainimaging by Bragg edge neutron transmission. Nuclear Instruments and Methods in Physics ResearchSection A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2002 Apr;481(1-3):765–768. Available from: https://doi.org/10.1016/s0168-9002(01)01256-6 .[4] Jørgensen JS, Ametova E, Burca G, Fardell G, Papoutsellis E, Pasca E, et al. Core Imaging Library– Part I: a versatile software framework for tomographic image processing. 2020;Preprint. Availablefrom: http://arxiv.org/abs/2102.04560 .[5] Warr R, Ametova E, Fardell G, Handschuh S, Jørgensen JS, Papoutsellis E, et al. Elemental mappingby low-count reconstruction for lab-based hyperspectral computed tomography. 2020;In preparation.[6] Ametova E, Burca G, Fardell G, Jørgensen JS, Papoutsellis E, Pasca E, et al. Phase discriminatingneutron tomography using advanced reconstruction methods. 2020;In preparation.[7] Beck A, Teboulle M. Fast Gradient-Based Algorithms for Constrained Total Variation Image De-noising and Deblurring Problems. IEEE Transactions on Image Processing. 2009;18(11):2419–2434.Available from: http://dx.doi.org/10.1109/TIP.2009.2028250 .[8] Chambolle A, Pock T. A First-Order Primal-Dual Algorithm for Convex Problems with Applicationsto Imaging. Journal of Mathematical Imaging and Vision. 2011;40(1):120–145. Available from: http://dx.doi.org/10.1007/s10851-010-0251-1 .[9] Chambolle A, Ehrhardt MJ, Richt´arik P, Sch¨onlieb CB. Stochastic primal-dual hybrid gradientalgorithm with arbitrary sampling and imaging applications. SIAM Journal on Optimization.2018;28(4):2783–2808. Available from: https://doi.org/10.1137/17M1134834 .[10] van Aarle W, Palenstijn WJ, Cant J, Janssens E, Bleichrodt F, Dabravolski A, et al. Fast and flexibleX-ray tomography using the ASTRA toolbox. Optics Express. 2016 Oct;24(22):25129. Availablefrom: https://doi.org/10.1364/oe.24.025129 .[11] Ovtchinnikov E, Brown R, Kolbitsch C, Pasca E, da Costa-Luis C, Gillman AG, et al. SIRF: Syner-gistic Image Reconstruction Framework. Computer Physics Communications. 2020 Apr;249:107087.Available from: https://doi.org/10.1016/j.cpc.2019.107087 .[12] Rudin L, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Phys-ica D: Nonlinear Phenomena. 1992;60:259–268. Available from: http://dx.doi.org/10.1016/0167-2789(92)90242-F .[13] Duran J, Moeller M, Sbert C, Cremers D. Collaborative Total Variation: A General Frameworkfor Vectorial TV Models. SIAM Journal on Imaging Sciences. 2016;9(1):116–151. Available from: https://doi.org/10.1137/15M102873X .[14] Bredies K, Kunisch K, Pock T. Total Generalized Variation. SIAM Journal on Imaging Sciences.2010;3(3):492–526. Available from: http://dx.doi.org/10.1137/090769521 .[15] Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: From error visibility tostructural similarity. IEEE Transactions on Image Processing. 2004;13:600–612. Available from: http://dx.doi.org/10.1109/TIP.2003.819861 .1716] Heikkil¨a T, Help H, Meaney A. Gel phantom data for dynamic X-ray tomography. Zenodo; 2020.Available from: https://zenodo.org/record/3696817 .[17] Heikkil¨a T, Help H, Meaney A. Gel phantom data for dynamic X-ray tomography; 2020. https://arxiv.org/abs/2003.02841 . Available from: http://dx.doi.org/10.5281/zenodo.3696817 .[18] Anthoine S, Aujol JF, Boursier Y, Melot C. Some proximal methods for Poisson intensity CBCTand PET. Inverse Problems & Imaging. 2012;6:565. Available from: https://doi.org/10.3934/ipi.2012.6.565 .[19] Sidky EY, Jørgensen JS, Pan X. Convex optimization problem prototyping for image reconstructionin computed tomography with the Chambolle-Pock algorithm. Phys Med Biol. 2012;10:3065–3091.Available from: https://doi.org/10.1088/0031-9155/57/10/3065 .[20] Burger M, Sawatzky A, Steidl G. First Order Algorithms in Variational Image Processing. In:Splitting Methods in Communication, Imaging, Science, and Engineering. Springer InternationalPublishing; 2016. p. 345–407. Available from: https://doi.org/10.1007/978-3-319-41589-5_10 .[21] Egan CK, S D M J, Wilson MD, Veale MC, Seller P, Beale AM, et al. 3D chemical imaging inthe laboratory by hyperspectral X-ray computed tomography. Scientific Reports. 2015;5(1):15979.Available from: https://doi.org/10.1038/srep15979 .[22] Warr R, Egan C, Papoutsellis E, Jørgensen JS, Ametova E, Cernik R, et al.. Hyperspectral X-rayCT data set of mineralised ore sample with Au and Pb deposits. Zenodo; 2020. Available from: https://doi.org/10.5281/zenodo.4157615 .[23] M¨unch B, Trtik P, Marone F, Stampanoni M. Stripe and ring artifact removal with combinedwavelet – Fourier filtering. Opt Express. 2009 May;17(10):8567–8591. Available from: https://doi.org/10.1364/OE.17.008567 .[24] Ehrhardt MJ, Markiewicz P, Sch¨onlieb CB. Faster PET reconstruction with non-smooth priorsby randomization and preconditioning. Physics in Medicine & Biology. 2019 nov;64(22):225019.Available from: https://doi.org/10.1088/1361-6560/ab3d07https://doi.org/10.1088/1361-6560/ab3d07