Core Imaging Library - Part I: a versatile Python framework for tomographic imaging

We present the Core Imaging Library (CIL), an open-source Python framework for tomographic imaging with particular emphasis on reconstruction of challenging datasets. Conventional filtered back-projection reconstruction tends to be insufficient for highly noisy, incomplete, non-standard or multi-channel data arising for example in dynamic, spectral and in situ tomography. CIL provides an extensive modular optimization framework for prototyping reconstruction methods including sparsity and total variation regularization, as well as tools for loading, preprocessing and visualizing tomographic data. The capabilities of CIL are demonstrated on a synchrotron example dataset and three challenging cases spanning golden-ratio neutron tomography, cone-beam X-ray laminography and positron emission tomography. This article is part of the theme issue ‘Synergistic tomographic image reconstruction: part 2’.


Introduction
It is an exciting time for computed tomography (CT): existing imaging techniques are being pushed beyond current limits on resolution, speed and dose, while new ones are being continually developed [1].Driving forces include higher-intensity X-ray sources and photon-counting detectors enabling respectively fast time-resolved and energy-resolved imaging.In situ imaging of evolving processes and unconventional sample geometries such as laterally extended samples are also areas of great interest.Similar trends are seen across other imaging areas, including transmission electron microscopy (TEM), positron emission tomography (PET), magnetic resonance imaging (MRI), and neutron imaging, as well as joint or multicontrast imaging combining several such modalities.
Critical in CT imaging is the reconstruction step where the raw measured data is computationally combined into reconstructed volume (or higher-dimensional) data sets.Existing reconstruction software such as proprietary programs on commercial scanners are often optimised for conventional, high quality data sets, relying on filtered back projection (FBP) type reconstruction methods [2].Noisy, incomplete, non-standard or multi-channel data will generally be poorly supported or not at all.
In recent years, numerous reconstruction methods for new imaging techniques have been developed.In particular, iterative reconstruction methods based on solving suitable optimisation problems, such as sparsity and total variation regularisation, have been applied with great success to improve reconstruction quality in challenging cases [3].This however is highly specialised and time-consuming work that is rarely deployed for routine use.The result is a lack of suitable reconstruction software, severely limiting the full exploitation of new imaging opportunities.
This article presents the Core Imaging Library (CIL) -a versatile open-source Python library for processing and reconstruction of challenging tomographic imaging data.CIL is developed by the Collaborative Computational Project in Tomographic Imaging (CCPi) network and is available from https: //www.ccpi.ac.uk/CIL as well as from [4], with documentation, installation instructions and numerous demos.
CIL aims to combine the best of the two worlds of tomography and optimisation software in a single easy-to-use, highly modular and configurable Python library.Particular emphasis is on enabling a variety of regularised reconstruction methods within a "plug and play" structure in which different data fidelities, regularisers, constraints and algorithms can be easily selected and combined.The intention is that users will be able to use the existing reconstruction methods provided, or prototype their own, to deal with noisy, incomplete, non-standard and multi-channel tomographic data sets for which conventional FBP type methods and proprietary software fail to produce satisfactory results.In addition to reconstruction, CIL supplies tools for loading, preprocessing, visualising and exporting data for subsequent analysis and visual exploration.
CIL easily connects with other libraries to further combine and expand capabilities; we describe CIL plugins for ASTRA [6], TIGRE [7] and the CCPi-Regularisation (CCPi-RGL) toolkit [16], as well as interoperability with the Synergistic Image Reconstruction Framework (SIRF) [17] enabling PET and MR reconstruction using CIL.
We envision that in particular two types of researchers might find CIL useful: • Applied mathematicians and computational scientists can use existing mathematical building blocks and the modular design of CIL to rapidly implement and experiment with new reconstruction algorithms and compare them against existing state-of-the-art methods.They can easily run controlled simulation studies with test phantoms and within the same framework transition into demonstrations on real CT data.• CT experimentalists will be able to load and pre-process their standard or non-standard data sets and reconstruct them using a range of different state-of-the-art reconstruction algorithms.In this way they can experiment with, and assess the efficacy of, different methods for compensating for poor data quality or handle novel imaging modalities in relation to whatever specific imaging task they are interested in.CIL includes a number of standard test images as well as demonstration data and scripts that make it easy for users of both groups to get started using CIL for tomographic imaging.These are described in the CIL documentation and we also highlight that all data and code for the experiments presented here are available as described under Data Accessibility.
This paper describes the core functionality of CIL and demonstrates its capabilities using an illustrative running example, followed by three specialised exemplar case studies.Section 2 gives an overview of CIL and describes the functionality of all the main modules.Section 3 focuses on the optimisation module used to specify and solve reconstruction problems.Section 4 presents the three exemplar cases, before a discussion and outlook are provided in Section 5. Multi-channel functionality (e.g. for dynamic and spectral CT) is presented in the part II paper [18] and a use case of CIL for PET/MR motion compensation is given in [19]; further applications of CIL in hyperspectral X-ray and neutron tomography are presented in [20] and [21].

Overview of CIL
CIL is developed mainly in Python and binary distribution is currently via Anaconda.Instructions for installation and getting started are available at https://www.ccpi.ac.uk/CIL as well as at [4].The present version 21.0 consists of six modules, as shown in Fig. 1.CIL is open-source software released under the Apache 2.0 license, while individual plugins may have a different license, e.g.ccpi.plugins.astra is GPLv3.In the following subsections the key functionality of each CIL module is explained and demonstrated, apart from ccpi.optimisation which is covered in Section 3.
As a running example (Fig. 2) we employ a 3D parallel-beam X-ray CT data set from Beamline I13-2, Diamond Light Source, Harwell, UK.The sample consisted of a 0.5 mm aluminium cylinder with a piece of steel wire embedded in a small drilled hole.A droplet of salt water was placed on top, causing corrosion to form hydrogen bubbles.The data set, which was part of a fast time-lapse experiment, consists of 91 projections over 180 • , originally acquired as size 2560-by-2160 pixels, but provided in [22] downsampled to 160-by-135 pixels.

Data readers and writers
Tomographic data comes in a variety of different formats depending on the instrument manufacturer or imaging facility.CIL currently supplies a native reader for Nikon's XTek data format, Zeiss' TXRM format, the NeXus format [23] if exported by CIL, as well as TIFF stacks.Here "native" means that a CIL AcquisitionData object incl.geometry (as described in the following subsection) will be created by the CIL reader.Other data formats can be read using e.g.DXchange [24] and a CIL AcquisitionData object can be manually constructed.CIL currently provides functionality to export/write data to disk in NeXus format or as a TIFF stack.
The steel-wire dataset is included as an example in CIL.It is in NeXus format and can be loaded using NEXUSDataReader.For example data sets in CIL we provide a convenience method that saves the user from typing the path to the datafile: Load steel-wire example dataset from cil.utilities.dataexampleimport SYNCHROTRON_PARALLEL_BEAM_DATA data = SYNCHROTRON_PARALLEL_BEAM_DATA.get()

Data structures, geometry and core functionality
CIL provides two essential classes for data representation, namely AcquisitionData for tomographic data and ImageData for reconstructed (or simulated) volume data.The steel-wire dataset was read in as an AcquisitionData that we can inspect with: At present, data is stored internally as a NumPy array and may be returned using the method as_array().AcquisitionData and ImageData use string labels rather than a positional index to represent the dimensions.In the example data, 'angle', 'vertical' and 'horizontal' refer to 91 projections each with vertical size 135 and horizontal size 160.Labels enable the user to access subsets of data without knowing the details of how it is stored underneath.For example we can extract a single projection using the method get_slice with the label and display it (Fig. 2

left) as
Extract single projection and display as image show2D(data.get_slice(angle=0),cmap='inferno', origin='upper-left') where show2D is a display functiontter in cil.utilities.display.show2D displays dimension labels on plot axes as in Fig. 2; subsequent plots omit these for space reasons.
Both ImageData and AcquisitionData behave much like a NumPy array with support for: • algebraic operators + , -, etc., • relational operators > , >= , etc., • common mathematical functions like exp, log and abs, mean, and • inner product dot and Euclidean norm norm.This makes it easy to do a range of data processing tasks.For example in Fig. 2 (left) we note the projection (which is already flat-field normalised) has values around 0.7 in the background, and not 1.0 as in typical well-normalised data.This may lead to reconstruction artifacts.A quick-fix is to scale the image to have background value ca.1.0.To do that we extract a row of the data toward the top, compute its mean and use it to normalise the data: Normalise data by mean over vertical slice of data data = data / data.get_slice(vertical=20).mean()Where possible in-place operations are supported to avoid unnecessary copying of data.For example the Lambert-Beer negative logarithm conversion can be done by: In-place mathematical operations data.log(out=data)data *= -1 The first line creates a default 3D parallel-beam geometry with a rotation axis perpendicular to the beam propagation direction.The second and third lines specify the detector dimension and the angles at which projections are acquired.Numerous configuration options are available for bespoke geometries; this is illustrated in Section 4.2, see in particular Fig. 9, for an example of cone-beam laminography.Similarly, ImageGeometry holds the geometric specification of a reconstructed volume, including numbers and sizes of voxels.

Preprocessing data
In CIL a Processor is a class that takes an ImageData or AcquisitionData as input, carries out some operations on it and returns an ImageData or AcquisitionData.Example uses include common preprocessing tasks such as resizing (e.g.cropping or binning/downsampling) data, flat-field normalization and correction for centre-of-rotation offset, see Table 1 for an overview of Processors currently in CIL.
We will demonstrate centre-of-rotation correction and cropping using a Processor.Typically it is not possible to align the rotation axis perfectly with respect to the detector, and this leads to well-known centre-of-rotation reconstruction artifacts.CIL provides different techniques to estimate and compensate, the simplest being based on cross-correlation on the central slice.First the Processor instance must be created; this is an object instance which holds any parameters specified by the user; here which slice to operate on.Once created the Processor can carry out the processing task by calling it on the targeted data set.All this can be conveniently achieved in a single code line, as shown in the first line below.
Afterwards, we use a Slicer to remove some of the empty parts of the projections by cropping 20 pixel columns on each side of all projections, while also discarding the final projection which is a mirror image of the first.This produces data90.We can further produce a subsampled data set data15 by using another Slicer, keeping only every sixth projection.

Auxiliary tools
This module contains a number of useful tools: • dataexample: Example data sets and test images such as the steel-wire dataset1 .
• display: Tools for displaying data as images, including the show2D used in the previous section and other interactive displaying tools for Jupyter notebooks.• noise: Tools to simulate different kinds of noise, including Gaussian and Poisson.
• quality measures: Mathematical metrics Mean-Square-Error (MSE) and Peak-Signal-to-Noise-Ratio (PSNR) to quantify image quality against a ground-truth image.Some of these tools are demonstrated in other sections of the present paper; for the rest we refer the reader to the CIL documentation.

CIL Plugins and interoperability with SIRF
CIL allows the use of third-party software through plugins that wrap the desired functionality.At present the following three plugins are provided: • cil.plugins.ccpiregularisation This plugin wraps a number of regularisation methods from the CCPi-RGL toolkit [16] as CIL Functions.• cil.plugins.astra:This plugin provides access to CPU and GPU-accelerated forward and back projectors in ASTRA as well as the filtered back-projection (FBP) and Feldkamp-Davis-Kress (FDK) reconstruction methods for parallel and cone-beam geometries.• cil.plugins.tigre:This plugin currently provides access to GPU-accelerated cone-beam forward and back projectors and the FDK reconstruction method of the TIGRE toolbox.Furthermore, CIL is developed to be interoperable with the Synergistic Image Reconstruction Framework (SIRF) for PET and MR imaging [17].This was achieved by synchronising naming conventions and basic class concepts: • sirf : Data structures and acquisition models of SIRF can be used from CIL without a plugin, in particular with cil.optimisation one may specify and solve optimisation problems with SIRF data.An example of this using PET data is given in Section 4.3.We demonstrate here how the cil.plugins.astraplugin, or cil.plugins.tigreplugin interchangeably, can be used to produce an FBP reconstruction of the steel-wire dataset using its FBP Processor.To compute a reconstruction we must specify the geometry we want for the reconstruction volume; for convenience, a default ImageGeometry can be determined from a given AcquisitionGeometry.The FBP Processor can then be set up and in this instance we specify for it to use GPU-acceleration, and then call it on the data set to produce a reconstruction: Set up and run GPU-accelerated FBP algorithm from ASTRA plugin data15.reorder(order='astra')ag = data15.geometryig = ag.get_ImageGeometry()recon = FBP(ig, ag, device='gpu')(data15) The first line permutes the underlying data array to the specific dimension order required by cil.plugins.astra,which may differ from how data is read into CIL.Reconstructions for both the 90-and 15-projection steel-wire datasets are seen in Fig. 3, with notable streak artifacts in the subsampled case, as is typical with few projections.

Reconstruction by solving optimisation problems
FBP type reconstruction methods have very limited capability to model and address challenging data sets.For example the type and amount of noise cannot be modelled and prior knowledge such as nonnegativity or smoothness cannot be incorporated.A much more flexible class of reconstruction methods arises from expressing the reconstructed image as the solution to an optimisation problem combining data and noise models and any prior knowledge.
The CIL optimisation module makes it simple to specify a variety of optimisation problems for reconstruction and provides a range of optimisation algorithms for their solution.

Operators
The ccpi.optimisation module is built around the generic linear inverse problem where A is a linear operator, u is the image to be determined, and b is the measured data.In CIL u and b are normally represented by ImageData and AcquisitionData respectively, and A by a LinearOperator.
The spaces that a LinearOperator maps from and to are represented in attributes domain and range; these should each hold an ImageGeometry or AcquisitionGeometry that match with that of u and b, respectively.Reconstruction methods rely on two essential methods of a LinearOperator, namely direct, which evaluates Av for a given v, and adjoint, which evaluates A * z for a given z, where A * is the adjoint operator of A. For example, in a LinearOperator representing the discretised Radon transform for tomographic imaging, direct is forward projection, i.e., computing the sinogram corresponding to a given image, while adjoint corresponds to back-projection.
Table 2 provides an overview of the Operators available in the current version of CIL.It includes imaging models such as BlurringOperator for image deblurring problems and mathematical operators such as IdentityOperator and GradientOperator to act as building blocks for specifying optimisation problems.Operators can be combined to create new Operators through addition, scalar multiplication and composition.
The bottom two row contains ProjectionOperator from both cil.plugins.astraand cil.plugins.tigre,which wraps forward and back-projectors from the ASTRA and TIGRE toolboxes respectively, and can be used interchangeably.A ProjectionOperator can be set up simply by Create ProjectionOperator from image and acquisition geometries A = ProjectionOperator(ig, ag) and from the AcquisitionGeometry provided the relevant 2D or 3D, parallel-beam or cone-beam geometry employed; in case of the steel-wire dataset, a 3D parallel-beam geometry.

Algebraic iterative reconstruction methods
One of the most basic optimisation problems for reconstruction is least-squares minimisation, where we seek to find the image u that fits the data the best, i.e., in which the norm of the residual Au − b takes on the smallest possible value; this u we denote u and take as our reconstruction.
The Conjugate Gradient Least Squares (CGLS) algorithm [25] is an algebraic iterative method that solves exactly this problem.In CIL it is available as CGLS, which is an example of an Algorithm object.The following code sets up a CGLS algorithm -inputs required are an initial image, the operator (here ProjectionOperator from cil.plugins.astra),the data and an upper limit on the number iterations to run -and runs a specified number of iterations with verbose printing: At this point the reconstruction is available as myCGLS.solutionand can be displayed or otherwise analysed.The object-oriented design of Algorithm means that iterating can be resumed from the current state, simply by another myCGLS.runcall.
As imaging operators are often ill-conditioned with respect to inversion, small errors and inconsistencies tend to magnify during the solution process, typically rendering the final least squares u useless.CGLS exhibits semi-convergence [26] meaning that in the initial iterations the solution will approach the true underlying solution, but from a certain point the noise will increasingly contaminate the solution.The number of iterations therefore has an important regularising effect and must be chosen with care.
CIL also provides the Simultaneous Iterative Reconstruction Technique (SIRT) as SIRT, which solves a particular weighted least-squares problem [27,9].As with CGLS, it exhibits semi-convergence, however tends to require more iterations.An advantage of SIRT is that it admits the specification of convex constraints, such as a box constraints (upper and lower bounds) on u; this is done using optional input arguments lower and upper: Set up and run SIRT algorithm with bounds on pixel values mySIRT = SIRT(initial=x0, operator=A, data=b, max_iteration=1000, \ lower=0.0,upper=0.09)mySIRT.run(200, verbose=1) In Fig. 4 we see that CGLS reduces streaks but blurs edges.SIRT further reduces streaks and sharpens edges to the background; this is an effect of the nonnegativity constraint.In the steel wire example data the upper bound of 0.09 is attained causing a more uniform appearance with sharper edges.

Tikhonov regularisation with BlockOperator and BlockDataContainer
Algebraic iterative methods like CGLS and SIRT enforce regularisation of the solution implicitly by terminating iterations early.A more explicit form of regularisation is to include it directly in an optimisation formulation.The archetypal such method is Tikhonov regularisation which takes the form where D is some operator, the properties of which govern the appearance of the solution.In the simplest form D can be taken as the identity operator.Another common choice is a discrete gradient implemented where the 0 corresponds to the range of D. We can use the CGLS algorithm to solve Eq. ( 4) but we need a way to express the block structure of Ã and b.This is achieved by the BlockOperator and BlockDataContainer of CIL: Set up Tikhonov regularisation for CGLS using BlockOperator and BlockDataContainer GradientOperator automatically works out from the ImageGeometry ig which dimensions are available and sets up finite differencing in all dimensions.If two or more dimensions are present, D will in fact be a BlockOperator with a finite-differencing block for each dimension.CIL supports nesting of a BlockOperator inside another, so that Tikhonov regularisation with a Gradient operator can be conveniently expressed.In Fig. 5 (left) Tikhonov regularisation with the GradientOperator is demonstrated on the steel-wire sample.Here, α governs the solution smoothness similar to how the number of iterations affects CGLS solutions, with large α values producing smooth solutions.Here α = 1 is used as a suitable trade-off between noise reduction and smoothing.
The block structure provides the machinery to experiment with different amounts or types of regularisation in individual dimensions in a Tikhonov setting.We consider the problem where we have different regularising operators D x , D y , D z in each dimension and associated regularisation parameters α x , α y , α z .We can write this as the following block least squares problem which can be solved  L 1 -norm: by CGLS: where 0 x , 0 y and 0 z represent zero vectors of appropriate size.In Fig. 5 we show results for D x , D y and D z being finite-difference operators in each direction, achieved by the FiniteDifferenceOperator.We show two choices of sets of regularisation parameters, namely α x = α y = 30, α z = 0.1 and α x = α y = 0.1, α z = 60.We see in the former case a large amount of smoothing occurs in the horizontal dimensions due to the larger α x and α y parameters, and little in the vertical dimension, so horizontal edges are preserved.In the latter case, opposite observations can be made.
Such anisotropic regularization could be useful with objects having a layered or fibrous structure or if the measurement setup provides different resolution or noise properties in different dimensions, e.g., for non-standard scan trajectories such as tomosynthesis/laminography.

Smooth convex optimisation
CIL supports the formulation and solution of more general optimisation problems.One problem class supported is unconstrained smooth convex optimisation problems, Here f is a differentiable, convex, so-called L-smooth function, that is its gradient , ∀u 1 , u 2 for some L > 0 referred to as the Lipschitz parameter.CIL represents functions by the Function class, which maps an ImageData or AcquisitionData to a real number.Differentiable functions provide the method gradient to allow first-order optimisation 2 .An overview of Function types currently in CIL is provided in Table 3.Another example using a smooth approximation of non-smooth total variation regularisation will be given in Section 4.1.

Non-smooth convex optimisation with simple proximal mapping
Many useful reconstruction methods are formulated as non-smooth optimisation problems.Of specific interest in recent years has been sparsity-exploiting regularisation such as the L 1 -norm and total variation (TV).TV-regularisation for example has been shown capable of producing high-quality images from severely undersampled data whereas FBP produces highly noisy, streaky images.A particular problem class of interest can be formulated as where f is L-smooth and g may be non-smooth.This problem can be solved by the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [28,29], which is available in CIL as FISTA.FISTA makes use of f being smooth by calling f.gradient and assumes for g that the so-called proximal mapping, for a positive parameter τ is available as g.proximal.This means that FISTA is useful when g is "proximable", i.e., where an analytical expression for the proximal mapping exists, or it can be computed efficiently numerically.A simple, but useful case, for FISTA is to enforce constraints on the solution, i.e., require u ∈ C, where C is a convex set.In this case g is set to the (convex analysis) indicator function of C, i.e., In CIL we can express such a function using a BlockOperator, as also used in the Tikhonov example, and a BlockFunction, which essentially holds a list of Function objects.
Here we demonstrate this setup by using PDHG to solve the TV-regularised least-squares problem.As shown in [33] this problem can be written in the required form as In CIL this can be written succinctly as (with a specific choice of regularisation parameter): Set up and run PDHG for TV-regularised least-squares problem alpha = 0.02 F = BlockFunction(L2NormSquared(b=b), alpha*MixedL21Norm()) K = BlockOperator(A, GradientOperator(ig)) G = ZeroFunction() myPDHG = PDHG(f=F, operator=K, g=G, max_iteration=10000) myPDHG.run(5000,verbose=2) Figure 7 shows the resulting steel-wire dataset reconstruction which appears identical to the result of FISTA on the same problem (Fig. 6), and as such validates the two algorithms against each other.
CIL Algorithms have the option to save the history of objective values so the progress and convergence can be monitored.PDHG is a primal-dual algorithm, which means that the so-called dual maximisation problem of Eq. ( 12), which is referred to as the primal problem, is solved simultaneously.In PDHG the dual objective values are also available.The primal-dual gap, which is the difference between the primal and dual objective values, is useful for monitoring convergence as it should approach zero when the iterates converge to the solution.
Figure 7 (right) compares the primal objective, dual objective and primal-dual gap history with the objective history for FISTA on the same problem.The (primal) objectives settle at roughly the same level, again confirming that the two algorithms achieve essentially the same solution.FISTA used fewer iterations, but each iteration took about 25 times as long as a PDHG iteration.The dual objective is negative until around 3000 iterations, and the primal-dual gap is seen to approach zero, thus confirming convergence.CIL makes such algorithm comparisons straightforward.It should be stressed that the particular convergence behavior observed for FISTA and PDHG depends on internal algorithm parameters such as step sizes for which default values were used here.The user may experiment with tuning these parameters to obtain faster convergence, for example for PDHG the primal and dual step sizes may be set using the inputs sigma and tau.
In addition to PDHG a stochastic variant SPDHG [34] that can sometimes accelerate reconstruction substantially by working on problem subsets is provided in CIL as SPDHG; this is demonstrated in the Part II article [18].
An overview of all the algorithms currently supplied by CIL is provided in Table 4.

Exemplar studies using CIL
This section presents 3 illustrative examples each demonstrating different functionality of CIL.All code and data to reproduce the results are provided, see Data Accessibility.

Neutron tomography with golden-angle data
This example demonstrates how CIL can handle other imaging modalities than X-ray, a non-standard scan geometry, and easily compare reconstruction algorithms.Contrary to X-rays, neutrons interact with atomic nuclei rather than electrons that surround them, which yields a different contrast mechanism, e.g., for neutrons hydrogen is highly attenuating while lead is almost transparent.Nevertheless, neutron data can be modelled with the Radon transform and reconstructed with the same techniques as X-ray data.
A benchmarking neutron tomography dataset (Fig. 8) was acquired at the IMAT beamline [35,36] of the ISIS Neutron and Muon Source, Harwell, UK.The raw data is available at [37] and a processed subset for this paper is available from [38].The test phantom consisted of an Al cylinder of diameter 22 mm with cylindrical holes holding 1mm and 3mm rods of high-purity elemental Cu, Fe, Ni, Ti, and Zn rods.186 projections each 512-by-512 pixels in size 0.055 mm were acquired using the non-standard golden-angle mode [39] (angular steps of 1 2 ( √ 5−1)•180 • = 111.24... • ) rather than sequential small angular increments.This was to provide complete angular coverage in case of early experiment termination and to allow experimenting with reconstruction from a reduced number of projections.An energy-sensitive micro-channel plate (MCP) detector was used [40,41] providing raw data in 2332 energy bins per pixel, which were processed and summed to simulate a conventional white-beam absorption-contrast data set for the present paper.Reconstruction and analysis of a similar energy-resolved data set is given in [21].
We use TIFFStackReader to load the data, several Processor instances to preprocess it, and initially FBP to reconstruct it.We compare with TV-regularisation, Eq. ( 11), solved with MixedL21Norm and PDHG using α = 1 and 30000 iterations, and further with a smoothed variant of TV (STV) using SmoothMixedL21Norm.The latter makes the optimisation problem smooth, so it can be solved using GD, using the same α and 10000 iterations.
The sinogram for a single slice is shown in Fig. 8 along with FBP, TV and STV reconstructions and a horizontal line profile plot as marked by the red line.The FBP reconstruction recovers the main sample features, however it is contaminated by noise, ring artifacts and streak artifacts emanating from the highest-attenuating rods.The TV and STV reconstructions remove these artifacts, while preserving edges.We see that the STV approximates the non-smooth TV very well; this also serves to validate the reconstruction algorithms against one another.

Non-standard acquisition: X-ray laminography
This example demonstrates how even more general acquisition geometries can be processed using CIL, and how cil.plugins.ccpiregularisation allows CIL to use GPU-accelerated implementations of regularising functions available in the CCPi-RGL toolkit [16].Furthermore, unlike the examples up to now, we here employ the ProjectionOperator provided by the TIGRE plugin, though the ASTRA plugin could equally have been used.
Laminography is an imaging technique designed for planar samples in which the rotation axis is tilted relative to the beam direction.Conventional imaging of planar samples often leads to severe limited-angle artifacts due to lack of transmission in-plane, while laminography can provide a more uniform exposure [42].In Transmission Electron Microscopy (TEM) the same technique is known as conical tilt.
An experimental laminography setup in the so-called rotary configuration was developed [43] for Nikon micro-CT scanners in the Manchester X-ray Imaging Facility.Promising reconstructions of a planar LEGO-brick test phantom were obtained using the CGLS algorithm.Here we use CIL on the same data [44] to demonstrate how TV-regularisation and non-negativity constraints can reduce inherent laminographic reconstruction artifacts.CIL allows the specification of very flexible scan configurations.The cone-beam laminography setup of the LEGO data set provides an illustrative case for demonstrating The data consists of 2512 projections of 798-by-574 pixels sized 0.508 mm in a 360 • cone-beam geometry.We load the data with NikonDataReader and preprocess with a couple of Processor instances to prepare it for reconstruction.For reconstruction we use the GPU-accelerated cone-beam ProjectionOperator from ccpi.plugin.tigreand FISTA to solve Eq. ( 8) for the unregularised leastsquares problem (LS) and non-negativity constrained TV-regularised least-squares (TVNN).For TVNN we use FBP_TV from cil.plugins.ccpiregularisation which implements a GPU-accelerated version of g TV , which is faster than, but otherwise equivalent to, using the native CIL TotalVariation.The full 3D volume is reconstructed for LS and TVNN, and Fig. 10 shows a horizontal and vertical slice through both.
The LEGO bricks are clearly visualised in all reconstructions.The LS reconstruction has a haze in the horizontal slice (top left), which in the vertical slice (bottom left) is seen to amount to smooth directional streaks known to be inherent for laminography; in particular horizontal edges are heavily blurred.On the other hand, fine details in the horizontal plane are preserved, for example the text "LEGO" seen on several knobs to the right.TVNN (right) reduces the haze and streaks substantially with the LEGO bricks displaying a uniform gray level and the horizontal edges in the vertical slice completely well-defined.However, some fine details are lost, including the "LEGO" text, which is a commonly observed drawback of TV-regularisation.Depending on the sample and application, this may or may not be an issue, and if necessary more sophisticated regularisers such as Total Generalised Variation (TGV) could be explored (a CIL example with TGV is given in the Part II article [18]).
As shown, CIL can process very general scan configurations and allows easy experimentation with different reconstruction methods, including using third-party software through plugins.

PET reconstruction in CIL using SIRF
SIRF (Synergistic Image Reconstruction Framework) [17] is an open-source platform for joint reconstruction of PET and MRI data developed by CCP-SyneRBI (formerly CCP-PETMR).CIL and SIRF have been developed with a large degree of interoperability, in particular data structures are aligned to enable CIL algorithms to work directly on SIRF data.As an example we demonstrate here reconstruction of the NEMA IQ Phantom [45], which is a standard phantom for testing scanner and reconstruction performance.It consists of a Perspex container with inserts of different-sized spheres, some filled with liquid with higher radioactivity concentration than the background, others with "cold" water (see [45] for more details).This allows assessment of resolution and quantification.
A 60-minute PET dataset [46] of the NEMA IQ phantom was acquired on a Siemens Biograph mMR PET/MR scanner at Institute of Nuclear Medicine, UCLH, London.Due to poor data statistics in PET a Poisson noise model is normally adopted, which leads to using the Kullback-Leibler (KL) divergence as data fidelity.We compare here reconstruction using the Ordered Subset Expectation Maximisation (OSEM) method [47] available in SIRF without using CIL, and TV-regularised KL divergence minimisation using CIL's PDHG algorithm with a KullbackLeibler data fidelity (KLTV).Instead of a CIL Operator a SIRF AcquisitionModel represents the forward model, and has all necessary methods to allow its use in CIL algorithms.
Figure 11 shows horizontal slices through the 220 × 220 × 127-voxel OSEM and KLTV reconstructions and vertical profile plots along the red line.In both cases the inserts are visible, but OSEM is highly affected by noise.KLTV reduces the noise dramatically, while preserving the insert and outer phantom edges.This may be beneficial in subsequent analysis, however a more detailed comparative study should take post-filtering into account.The purpose of this example was to give proof of principle of prototyping new reconstruction methods for PET with SIRF, using the generic algorithms of CIL, without needing to implement dedicated new algorithms in SIRF.Another example with SIRF for PET/MR motion compensation employing CIL is given in [19].

Summary and outlook
We have described the CCPi Core Imaging Library (CIL), an open-source library, primarily written in Python, for processing tomographic data, with particular emphasis on enabling a variety of regularised reconstruction methods.The structure is highly modular to allow the user to easily prototype and solve new problem formulations that improve reconstructions in cases with incomplete or low-quality data.We have demonstrated the capability and flexibility of CIL across a number of test cases, including parallel-beam, cone-beam, non-standard (laminography) scan geometry, neutron tomography and PET using SIRF data structures in CIL.Further multi-channel cases including temporal/dynamic and spectral tomography are given in [18].
CIL remains under active development with continual new functionality being added, steered by ongoing and future scientific projects.Current plans include: • adding more algorithms, functions, and operators to support an even greater set of problems, for example allow convex constraints in smooth problems; • adding more pre-/postprocessing tools, for example to handle beam hardening; • adding templates with preselected functions, algorithms, etc. to simplify solving common problems such as TV regularisation; • further integrating with other third-party open-source tomography software through the plugin capability; • introducing support for nonlinear problems, such as polarimetric neutron spin tomography [48] and electron strain tomography [49]; and • developing support for multi-modality problems.CIL is developed as open-source on GitHub, and questions, feature request and bug reports submitted as issues are welcomed.Alternatively, the developer team can be reached directly at CCPI-DEVEL@ jiscmail.ac.uk.CIL is currently distributed through the Anaconda platform; in the future additional modes of distribution such as Docker images may be provided.Installation instructions, documentation and training material is available from https://www.ccpi.ac.uk/cil as well as at [4], as are GitHub repositories with source code that may be cloned/forked and built manually.In this way users may modify and contribute back to CIL.
Finally we emphasize that a multitude of optimization and regularization methods exist beyond those currently implemented in CIL and demonstrated in the present article.Recent overviews are given for example by [50,51,52,3] with new problems and methods constantly being devised.CIL offers a modular platform to easily implement and explore such methods numerically as well as apply them directly in large-scale imaging applications.the neutron data.GF carried out the laminography case study and developed the CIL software.EPap carried out the PET case study and developed the CIL software.EPas conceived of and developed the CIL software and interoperability with SIRF.KT contributed to the PET case study, interoperability with SIRF and development of the CIL software.RW assisted with case studies and contributed to the CIL software.MT, WL and PW 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.

Figure 1 :
Figure 1: Overview of CIL module structure and contents.The cil.plugins module contains wrapper code for other software and third-party libraries that need to be installed separately to be used by CIL.

Figure 7 :
Figure 7: PDHG reconstruction of 15-projection 3D steel-wire dataset.Left two: TV-regularisation with α = 0.02, reproduces same result as FISTA in Fig. 6 on same case and parameter choice, thus validating algorithms against each other.Colour range [-0.01,0.11].Right: Objective value histories (log-log) for FISTA and PDHG on TV-regularisation problem.Both algorithms reach the same (primal) objective value, FISTA taking fewer but slower iterations.The primal-dual gap for PDHG (difference between primal and dual objectives) approaches zero indicating convergence.

Figure 9 :
Figure 9: CIL AcquisitionGeometry and ImageGeometry illustrated for the laminography cone-beam setup.Configurable parameters are shown in the legend.Parallel-beam geometry and 2D versions are also available.CIL can illustrate ImageGeometry and AcquisitionGeometry instances as in this figure using show_geometry.

Figure 10 :
Figure 10: Slices through 3D reconstruction of laminography LEGO sample.Left, top/bottom: LS reconstruction using FISTA, horizontal/vertical slice at yellow line.Right: Same using TVNN, in which laminography artifact are suppressed while edges preserved.

Table 1 :
Processors currently available in CIL.ImageGeometry and AcquisitionGeometry objects available in the attribute geometry of ImageData and AcquisitionData.AcquisitionGeometry will normally be provided as part of an AcquisitionData produced by the CIL reader.It is also possible to manually create AcquisitionGeometry and ImageGeometry from a list of geometric parameters.Had the steel-wire dataset not had geometry information included, we could have set up its geometry with the following call:

Table 2 :
Operators in CIL; and Operators from cil.plugins.astraand cil.plugins.tigre in bottom two rows.
Name Description BlockFunction Separable sum of multiple functions ConstantFunction Function taking the constant value OperatorCompositionFunction Compose function f and operator A: f (Ax) IndicatorBox Indicator function for box (lower/upper) constraints KullbackLeibler Kullback-Leibler divergence data fidelity L1Norm
i.e., where the composite function f (K•) can be written as a separable sum