MARS-MD: rejection based image domain material decomposition

This paper outlines image domain material decomposition algorithms that have been routinely used in MARS spectral CT systems. These algorithms (known collectively as MARS-MD) are based on a pragmatic heuristic for solving the under-determined problem where there are more materials than energy bins. This heuristic contains three parts: (1) splitting the problem into a number of possible sub-problems, each containing fewer materials; (2) solving each sub-problem; and (3) applying rejection criteria to eliminate all but one sub-problem's solution. An advantage of this process is that different constraints can be applied to each sub-problem if necessary. In addition, the result of this process is that solutions will be sparse in the material domain, which reduces crossover of signal between material images. Two algorithms based on this process are presented: the Segmentation variant, which uses segmented material classes to define each sub-problem; and the Angular Rejection variant, which defines the rejection criteria using the angle between reconstructed attenuation vectors.


Introduction
Spectral CT is a variant of x-ray computed tomography (CT) which uses additional measurements of photon energy to enable identification and quantification of different materials simultaneously.Energy information can be obtained through a multitude of ways.Clinical dual-energy CT is typically performed using either dual-source [1], kVp switching [2][3][4], or dual-layered detectors [5,6].More recently there has been development of pre-clinical systems which utilise energy discriminating photon counting detectors to provide up to eight channels of energy information [7][8][9][10].The MARS (Medipix All Resolution System) is a spectral CT system based on the Medipix3RX photon counting detector.This detector has eight energy bins which include one arbitrated, four charge summed, three single pixel counters.Information about these different counting modes can be found elsewhere [11].Recent MARS systems typically only utilize the arbitration and/or charge-summed counters for spectral analysis due to their superior energy resolution, yielding 4-5 unique energy bins for each measurement.
The process of converting energy information into material information is known as either basis or material decomposition (MD).This was first demonstrated to be possible by Alzarez and Macovski through the conversion of dual-energy measurements into photo-electric and Compton cross-section components [12].Later this was extended to estimating effective density and atomic number [13,14], and material specific densities [15][16][17].There are three different ways in which MD can be performed: pre-reconstruction, post-reconstruction, and (simultaneous) joint-reconstruction.Each of these have their own advantages and disadvantages.Pre-and joint-reconstruction allows for modelling of the x-ray beam poly-chromaticity, which in turn avoids beam hardening artefacts.To do this however requires use of non-linear optimization techniques which can be computationally intensive.Post-reconstruction MD has the advantage that it can be performed quickly, however it is susceptible to any beam hardening artefacts within the reconstructed data.The impact of this can be minimized by, using synchrotron x-ray sources, or applying beam-hardening correction techniques, and in some cases of measuring within narrow energy ranges.
One issue that is common to all three approaches is the linear dependencies between the attenuation distributions for materials commonly of interest in medical imaging.It was originally assumed that only two materials without K-edges could be decomposed due to the primary dependence on the photo-electric and Compton interactions in the diagnostic imaging range [18].Bornefalk later showed that the attenuation for the set of elements Z < 20 can theoretically be decomposed into four independent components [19].Despite this, the decomposition of many materials simultaneously can still prove challenging if any of them are close in effective atomic number.This paper outlines two post-reconstruction MD algorithms that has been routinely used in MARS small animal scanners, and the heuristic they are based on to minimize the impact of this basis linear dependence issue.We also describe the image domain formulation of post-reconstruction MD problem.Use of the presented algorithms is demonstrated on three multi-material phantoms.In addition, references are given to several publications describe the use and results of MARS-MD in different pre-clinical applications.

Image domain Material Decomposition
MD in the image domain processes reconstructed energy bin volumes into material image volumes.Each voxel of the reconstructed image provides an estimate of the attenuation in that voxel for each energy bin.That voxel's attenuation for a given energy bin can be parametrized by a linear combination of basis functions Using this, the attenuation for all energy bins can be written as the matrix equation where b = µ(E) is a column vector containing the attenuation of the voxel in each energy bin; ] is a matrix whose columns containing values of the basis functions at the given energies; and b = [a 1 a 2 . . .a N ] T is the column vector containing the basis function coefficients.Here we use E to represent the vector of energy bins with a slight abuse of notation.
If each energy bin is capturing a monochromatic signal then the elements of E would be just the energies of the corresponding x-rays.When working with polychromatic energy bins the elements of E are more abstractly thought of as observable energy intervals.It is not uncommon for this to be defined using effective energies to simplify the description of each energy bin.The image domain MD problem is simply the inversion of this matrix equation (x = M −1 b ) for every voxel in the volume.
Numerous approaches have been taken for solving this problem, many of which are based on a linear least squares estimates [15,20].Variations include application of different constraints such as volume conservation [21], mass conservation [17], or non-negativity [22]; each voxel being treated either independent or dependent of its neighbours [23]; and the representation of a voxel being restricted to a subset of basis functions [15,24].

MARS-MD Heuristic Process
The MARS-MD algorithms are based on a simple heuristic which is designed to be adaptable to different situations.Both algorithms presented below are based on the same three steps: given a set of N basis functions, (1) construct a number of feasible sub-problems from the basis set, (2) calculate the solution to every sub-problem, and (3) apply some rejection criteria to eliminate all but one sub-problem solution.
When determining what set of sub-problems to use we first choose the maximum number k of basis functions that will be allowed to describe a voxel (in practice this might be the maximum number of materials per voxel).We then consider all k ≤ N combinations of basis functions and further eliminate any that are non-physical and/or undesirable.A sub-problem is then formulated to find a solution to Eqn. 2.2 for each remaining sub-set of basis functions.Since each of these are calculated independently, different constraints can be applied to different sub-problems.

MARS-MD Algorithms
There are two versions of the MARS-MD algorithm following this heuristic that have been used with MARS systems -the Segmentation variant and the Angular Rejection variant.These two algorithms have been known by different names in different places.In the MARS system automated image processing chain they are known as MARS-MD v1.0 and MARS-MD v1.1 respectively.Each of these methods only considers solutions with at most two materials per voxel.For each voxel, the MARS-MD algorithms typically decompose reconstructed effective linear attenuation onto a set of mass attenuation bases for different materials.The results of the algorithms are produced as estimated densities of the given basis materials.

Sub-problem selection:
This uses a segmentation procedure we have previously published [24] to categorize each voxel to contain either air, soft tissue, or dense material.Combinations of materials within a voxel are then restricted based on its class.Air voxels are considered to be empty and set to the zero solution.Soft tissue voxels are treated as containing a combination of lipid and water obeying volume conservation.Dense material voxels which contain either a single dense material, or a single dense material in combination with water.

Sub-problem solution:
Eqn. 2.2 is solved for each basis combination using a non-negative linear least squares estimate.
The volume constraint is applied to the lipid-water sub-problem as an additional equation included within the matrix equation.To do this the lipid and water bases used in the matrix equation need to be expressed as linear attenuation.

Solution rejection:
The linear least squares error ||b − Mx * || 2 is calculated for each sub-problem solution x * .The solution with the smallest error is selected as the final solution for the respective voxel.

Sub-problem selection:
Each voxel is treated as possibly containing either a combination of lipid and water, a single dense material, or a single dense material in combination with water.These are the same basis combinations as used for the Segmentation variant, but without the additional layer of exclusion from the segmentation.

Sub-problem solution:
Eqn. 2.2 is solved for each basis combination using a non-negative linear least squares estimate.

Solution rejection:
An attenuation vector b * is calculated from each solution.These and the voxel's reconstructed attenuation estimate (b) are normalized by the corresponding attenuation for water (this is a pseudo-Hounsfield unit as outlined in [25]).The cosine angle between the normalized b * and b vectors are calculated for each solution.The solution with the smallest cosine angle is selected as the final solution for the respective voxel.

Algorithm Calibration
In its implementation within MARS systems, the M matrix in the MARS-MD algorithm is calibrated using scans of phantoms which contain the materials of interest at known concentrations and locations (Fig. 1).This is a two step process.First the effective linear attenuation for each material (for each concentration and energy bin) is estimated by taking the mean of respective regions in the reconstructed data.These values are then used to calculate the effective mass attenuation, which is used in the M matrix.If more than one concentration is used to calculate the mass attenuation for a given material then it is determined using the non-negative linear least squares estimate of Eqn.2.1 (where the mass attenuation is taken to be the unknown and concentration taken to be the known variables).The three phantoms that were scanned are known as the CaAu, GdI, and Multi-Contrast QA phantoms (Fig. 2).Each phantom is made from a 31 mm diameter PMMA, with nine 6 mm capillaries -eight around periphery and one in the centre.0.2 ml polypropylene Eppendrof tubes containing solutions of materials are inserted into each capillary.Each tube is sealed with paraffin wax.The Multi-Contrast phantom is used for MD calibration and contains vegetable oil; 2x water; 8 and 2 mg(Au)/ml gold chloride; 8 and 2 mg(Gd)/ml diluted MultiHance; 18 mg(I)/ml diluted Omnipaque 350; and 240 mg(Ca)/ml calcium chloride.The CaAu and GdI phantoms cover the same materials expect at a larger range of concentrations.The CaAu phantom contains water; 8, 4, 2, and 1 mg(Au)/ml gold chloride; and 240, 140, 70, and 35 mg(Ca)/ml calcium chloride.The GdI phantom contains water; 8, 4, 2, and 1 mg(Gd)/ml diluted MultiHance; and 18, 9, 4.5, and 2.25 mg(I)/ml diluted Omnipaque 350.

Multi-Contrast Phantom
The reconstruction of a single slice in the middle of the Multi-Contrast phantom is shown in Fig. 3.This was processed using three different material decomposition algorithms (Fig. 4), including the two MARS-MD variants and an unmodified non-negative linear least squares (LSQ).Each MD decomposed five energy bins into six material channels.The MD calibration was calculated from an average 11 slices (including the slice used for MD analysis) with material regions selected to cover the majority of each capillary in a similar way as shown in Fig. 1 (a).
The results from both MARS-MD variants demonstrate reasonable ability to separate the contents of each of the nine capillaries into their appropriate material channels.The exception to this being in the 2mg(Au)/ml and 2mg(Gd)/ml capillaries, where there is difficulty separating low concentration signal from the water background.Due to the enforced sparsity of at most two materials per voxel, there is less crossover of signal between material channels as compared to the non-negative LSQ results.A small amount of degradation of the decompositions are observed in the bottom left quadrant of the phantom, which result from minor beam hardening artefacts in the lower energy bins.
The most striking difference between the two MARS-MD variants is how the PMMA signal is split over different material channels.The PMMA has been mostly identified as a combination of gadolinium and water in the Segmentation variant.In this case, the higher signal density of the PMMA has enabled it to be grouped under the dense material category of the segmentation.This makes PMMA appear in the water channel, unlike the other two algorithms where it appears in the lipid channel.The lipid-water volume constraint used in the Segmentation variant also has a tendency of forcing PMMA signal into the water channel under situations where it processed in the soft tissue category.Despite the Angular Rejection variant not describing PMMA with a gadolinium component in this example, we have observed that it can happen when using different data acquisition parameters.

Concentration Phantoms
Two phantoms containing four concentrations each of gold, gadolinium, iodine, and calcium solutions were used to assess the quantitative results of the two MARS-MD variants.The concentrations in these phantoms were chosen to cover ranges of both good and poor detectability with the current MARS system.The decomposition of a single slice of these phantoms (excluding the lipid and water channels) is shown in Fig. 5. Line profiles presenting the MD estimated concentrations in selected capillaries is shown in Figs. 6 and 7.The two lowest concentrations of gadolinium, gold, and calcium have been excluded from the line profiles due to poor identification across the respective capillaries.Due to the way each algorithm has been constructed, the quantitative results for each should be almost equivalent.There are two main exceptions to this.Firstly, the gadolinium distribution is significantly different for reasons stated above.Secondly, the Segmentation variant has a tendency to identify slightly different material boundaries due to the restrictions on voxel materials from its initial segmentation procedure.This is most clearly demonstrated by the 18 mg(I)/ml capillary in Fig. 7 (b).
The lowest two concentrations of calcium (Fig 5) are mostly identified as gold.Misidentification such as this is not uncommon at low concentrations.It is also common observe low concentrations of gold misidentified as calcium when using different imaging parameters [26].
Most of the predicted concentrations shown in the line profile diagrams (Fig. 7) is a small underestimate when compared to the expected concentrations (horizontal dotted lines).Small deviations are expected in image space decomposition for a variety of reasons.For example -despite using narrow energy bins, beam hardening artefacts will still be present.In addition, the extent of the beam hardening will be different in each energy bin.It therefore is not surprising that we have observed that the size of this underestimation changes with modifications to energy bins and x-ray spectrum.

Pre-clinical Examples
Examples of previous publications where the Segmentation variant of the MARS-MD algorithm has been used include: soft tissue imaging using lamb meat [27]; arthritic cartilage imaging [28]; intrinsic biomarker and contrast identification in atherosclerosis [29]; and simultaneous imaging of multiple contrast agents [26].Each of these studies use different imaging parameters and provide independent assessments of MD results that are relevant to the given studies.The Angular Rejection variant is more recent and its applications are currently being investigated.

Discussion
Whereas development of MD techniques in the wider community has largely moved onto iterative joint-reconstruction methods, there is still benefit in further development of image based  techniques.Iterative joint-reconstruction techniques are typically very slow when compared to taking the image domain MD route.A more practical approach to the MD problem in spectral CT would be to construct an initial approximation of the decomposed volume using a fast filtered back projection (or similarly fast) reconstruction paired with an image space material decomposition, the results of which are then used for the initial guess in a joint-reconstruction decomposition.
The MARS-MD heuristic has a combinatorial approach to solving the MD problem.The typical issues associated with scalability have been minimized by only considering a small subset of material combinations.For the subsets used in the variants above the number of combinations per voxel scales linearly with the number of materials used in the decomposition (2N − 3 combinations for N ≥ 2 materials).This approach also provides an easy framework to reject solutions in lieu of others for unique reasons, such as unrealistic material concentrations.
One of the key features of the above MARS-MD algorithms is that they target solutions that are sparse in the material domain.Although this type of constraint has been used in other techniques [15,30], it is more common to find techniques which either use conservation constraints or promote spatial uniformity.Strict material sparsity as has been used here has the advantage that it minimizes the risk of fitting materials with K-edges in the imaging range to the noise.It can also turn an under-determined problem (more materials than energy bins) into an over-determined problem.There are drawbacks however.If a region contains more materials than the number allowed by the sparsity condition then it is not possible to correctly identify the contents.The MARS-MD algorithms described above, in this sense, do not facilitate either the identification of multiple high contrast materials within a single voxel nor the combination of high contrast materials and lipid.
One of the challenges of post-reconstruction MD is the question about how to calibrate them for a wide variety of applications and/or scan samples.Here we have calibrated directly from experimental phantom data.The validity of calibrations for each scan however are dependent on the degree of beam hardening within the sample.This in-turn depends on the size and composition of each scan sample -which will be different.Numerous practices have been followed in within the MARS research team's pre-clinical work to minimize this bias.The phantoms used for calibration are of a similar size and effective density to the samples generally scanned.Sometimes additional calibration vials have been included alongside initial samples to test the validity of the MD calibration and to correct it if necessary.
Lastly, the quality of results produced by any MD technique is highly dependent on the infor-mation measured from the spectral domain.What counts as good or sufficient spectral information can depend on things such as the decomposition bases that are used; what and how many energy bins are used; and even the degree of beam hardening caused by the sample.In other words, a set of optimal acquisition parameter for one application may not be optimal for other applications.We have given reference to a diverse range of applications using our methodology for this reason.This problem of spectral CT scan protocol optimality is currently an area of active research.

Conclusion
This paper has described two constrained image space material decomposition techniques that have been used in MARS spectral CT systems over the past several years.This description includes an outline of the heuristic that the MARS-MD algorithms have been constructed from; the commonly used Segmentation variant algorithm; and the more recent Angular Rejection variant algorithm.
A brief example of their application to spectral CT scans of generic multiple-contrast phantoms demonstrates reasonable material identification performance.

Figure 1 :
Figure 1: Example of a MARS-MD calibration using reconstructed images of a known phantom; (a) regions covering each material in the reconstructed images are averaged to estimate the effective linear attenuation for each material; (b) from this the effective mass attenuation for each material is estimated and used as the columns in the M matrix to calibrate the MARS-MD algorithm.

Figure 4 :
Figure 4: Material decomposition of Multi-Contrast Phantom using three different methodsthe Angular Rejection and Segmentation variants of MARS-MD, along with an unmodified nonnegative linear least squares decomposition.The MARS-MD variants show reasonable separation of the materials contained within each capillary.

Figure 5 :
Figure 5: Material decomposition of concentration phantoms using (a) Segmentation variant and (b) Angular Rejection variant of MARS-MD showing gold, gadolinium, calcium, and iodine channels fused with the highest energy bin.

Figure 6 :
Figure 6: Positions of line profiles displayed in Fig. 7 for GdI (left) and CaAu (right) phantoms.Profiles of low concentration capillaries with high levels of misidentification have been excluded.

Figure 7 :
Figure 7: Select line profiles of material decomposed concentration phantoms.Horizontal dotted lines indicate real concentrations of given materials.The profiles cover the following capillaries: (a) 8 mg(Gd)/ml and 4 mg(Gd)/ml; (b) 18 mg(I))/ml and 9 mg(I)/ml; (c) 4.5 mg(I)/ml and 2.25 mg(I)/ml; (d) 8 mg(Au)/ml and 4 mg(Au)/ml; and (e) 240 mg(Ca)/ml and 140 mg(Ca)/ml.In this example, most MD results are slightly less than expected yet still close to their desired concentrations.