Joint bi-modal image reconstruction of DOT and XCT with an extended Mumford-Shah functional
Di He, Ming Jiang, Alfred K. Louis, Peter Maass, Thomas Page
aa r X i v : . [ ee ss . I V ] O c t Joint bi-modal image reconstruction of DOT andXCT with an extended Mumford-Shah functional
Di He ‡ , Ming Jiang , , Alfred K. Louis , Peter Maass ,Thomas Page § School of Applied Science, Beijing Information Science & TechnologyUniversity, Beijing 100192, China. LMAM, School of Mathematical Sciences, Peking University, Beijing 100871,China. Beijing International Center for Mathematical Research, and CooperativeMedianet Innovation Center, Peking University, Beijing 100871, China. Institute of Applied Mathematics, Saarland University, 66041 Saarbr¨ucken,Germany. Center for Industrial Mathematics, University of Bremen, 28359 Bremen,Germany. Automated Driving Lab China, BMW China Services LTD.E-mail: [email protected], [email protected], [email protected],[email protected], [email protected].
Abstract.
Feature similarity measures are indispensable for joint imagereconstruction in multi-modality medical imaging, which enable joint multi-modalimage reconstruction (JmmIR) by communication of feature information from onemodality to another, and vice versa. In this work, we establish an image similaritymeasure in terms of image edges from Tversky’s theory of feature similarity inpsychology. For joint bi-modal image reconstruction (JbmIR), it is found thatthis image similarity measure is an extended Mumford-Shah functional with apriori edge information proposed previously from the perspective of regularizationapproach. This image similarity measure consists of Hausdorff measures of thecommon and different parts of image edges from both modalities. By construction,it posits that two images are more similar if they have more common edges andfewer unique/distinctive features, and will not force the nonexistent structuresto be reconstructed when applied to JbmIR. With the Γ-approximation of theJbmIR functional, an alternating minimization method is proposed for the JbmIRof diffuse optical tomography and x-ray computed tomography. The performanceof the proposed method is evaluated by three numerical phantoms. It is foundthat the proposed method improves the reconstructed image quality by more than10% compared to single modality image reconstruction (SmIR) in terms of thestructural similarity index measure (SSIM)
Keywords : Single-modal image reconstruction (SmIR), joint multi-modal imagereconstruction (JmmIR), joint bi-modal image reconstruction (JbmIR), imagesimilarity measures, extended Mumford-Shah functional, diffuse optical tomography(DOT), x-ray computed tomography (XCT). ‡ Part of the work was conducted while Di He was a PH.D student of Peking University. § Part of the work was conducted while Thomas Page was a PH.D student of University Bremanand visiting student of Peking University. oint bi-modal image reconstruction
1. Introduction
Biomedical imaging aims at visualizing structural or functional information nec-essary for medical research and clinical diagnosis. Each imaging modality canprovide images of one particular physical, physiological, or biological distribution.Multi-modality medical imaging technique is to combine multiple imaging modal-ities into one hybrid imaging system, such as PET/MRI [Judenhofer et al., 2008,Catana, 2017], PET/XCT [Cherry, 2009, Bockisch et al., 2009, Delbeke et al., 2009,Even-Sapir et al., 2009, Kaufmann and Di Carli, 2009], SPECT/XCT [Cherry, 2009,Bockisch et al., 2009,Delbeke et al., 2009,Even-Sapir et al., 2009,Kaufmann and Di Carli, 2009],XCT/MRI [Wang et al., 2015], DOT/MRI [Panagiotou et al., 2009] and DOT/XCT[Yuan et al., 2010,Boas, 2014,Deng et al., 2015,Baikejiang et al., 2017], BLT/DOT/XCT[Yang et al., 2015], and an omni-tomography system of multiple modalities[Wang et al., 2012]. k Hybrid multi-modality imaging systems can visualize anatomi-cal and functional structures simultaneously, improve the overall imaging performance,especially avoiding the tempo-spatial artifacts because of scanning with different de-vices at different time and positions [Townsend, 2008]. It offers significant diagnosticadvantages that cannot be achieved by a single modality [Steiner, 2007]. Moreover, ifhybrid multi-modality imaging systems are perfectly calibrated, image registration isnot necessary. Multi-modality medical imaging within one hybrid imaging system iscalled hardware fusion in [Townsend, 2008].For single-modal image reconstruction (SmIR), theories as well as algorithms havealready been well developed [Natterer et al., 2002,Scherzer, 2010,Censor et al., 2008].Image reconstructions for hybrid multi-modality imaging systems can be conductedin several ways. The first approach is that all the image reconstructions areconducted separately by multiple SmIRs without sharing information across imagingmodalities, which is a simple application of multiple SmIRs for each modality. Thesecond approach is that all the image reconstructions are conducted sequentially byusing reconstructed images to guide the next image reconstruction. The sequentialreconstructions are usually started from high resolution imaging modalities such asXCT or MRI. The structural information of reconstructed images is used to guidethe next reconstructions by applying structural similarity. This approach is based onthe observation that images of the same object possess similar structural information,especially at key structures or strong image edges, though with different contrastsfrom different modalities, because they are images of the same anatomical structure.This approach is called model fusion in [Haber and Gazit, 2013] or software fusion in [Townsend, 2008]. The third approach is that all the image reconstructions arejointly conducted by sharing reciprocally structural information from one modality toanother, and vice versa. This is called joint inversion in [Haber and Gazit, 2013] andjoint multi-modal image reconstruction (JmmIR) in this paper. For the joint bi-modalimage reconstruction such as DOT/XCT in this paper, it is called joint bi-modal imagereconstruction (JbmIR).JmmIR can be achieved by using iteratively model fusion with an appropriateimage similarity measure of structural information in certain feature space for imagesof the underlying modalities [Haber and Gazit, 2013]. Given such an image similarity k The abbreviations are as follows in the order of appearance: positron emission tomography (PET),magnetic resonance imaging (MRI), x-ray computed tomography (XCT), single-photon emissioncomputed tomography (SPECT), diffuse optical tomography (DOT), bioluminescence tomography(BLT). oint bi-modal image reconstruction a priori edge informationproposed previously for model fusion from the perspective of regularization approachin [Page, 2015]. This image similarity measure consists of Hausdorff measures ofthe common and different parts of image edges from DOT and XCT, and can beapproximated with edge indicator functions by the Γ-convergence theory. With theΓ-approximation of the JbmIR functional, an alternating minimization method isproposed for the JbmIR of DOT and XCT. The performance of the proposed methodis evaluated with three numerical phantoms. It is found that the proposed methodimproves the image quality by more than 10% in terms of the structural similarityindex measure (SSIM) for both XCT and DOT.The structure of the paper is as follows. We first review the previous work onmodel fusion and JmmIR in §
2. Then in § §
4. In § § §
7, we perform numerical experiments with threephantoms with single modal image reconstruction with the Mumford-Shah (SmIR-MS) regularization with our JbmIR algorithm (JbmIR-MS). We discuss the results,relevant issues and future work in §
8. We conclude this paper in §
2. Related work
Since JmmIR can be achieved by using iteratively model fusion, we review the previousmethods for model fusion and discuss the challenges for JmmIR in this section. Thefollowing classification of methods is based on the major characteristics of methods,and is used to address the challenges for JmmIR, but is not intended to be a strictclassification. There are certainly overlaps among the classifications, e.g., derivativesor its discretization, finite differences, are used explicitly or implicitly in almost allthe methods. oint bi-modal image reconstruction For model fusion with MRI or CT, image edges from MRI or CT are used toreplace the line process in the Gibbs distribution of [Geman and Geman, 1984]. Thisapproach is first used for the reconstruction of PET or SPECT [Leahy and Yan, 1991,Fessler et al., 1992, Gindi et al., 1993]. If the edge information from MRI is imperfectand if the edge process is improperly weighted, using such information “blindly”can even lead to artifacts [Fessler et al., 1992,Zhang et al., 1994,Bowsher et al., 1996].In [Bowsher et al., 1996], a simultaneous segmentation and reconstruction methodfor emission computed tomography (ECT) image by incorporating high-resolutionanatomical information such as CT or MRI is proposed and evaluated with SPECTand MRI. Higher prior probabilities are assigned to ECT segmentations in whicheach ECT region stays within a single anatomical region in this incorporation[Bowsher et al., 1996]. In [Sastry and Carson, 1997] , a different method for theincorporation of anatomical information into PET image reconstruction is proposed:segmentations of MRI images are used to assign tissue composition to PET imagepixels, while PET images are modeled as the sum of activities of each tissue type,weighted by the assigned tissue composition.If segmentation from a high resolution imaging modality is available, apriori structural information is used for DOT, as the weight for an anisotropicregularization to preserve the edges in the DOT image under reconstruction in[Douiri et al., 2007]. The idea of [Douiri et al., 2007] is later applied for PET in[Chan et al., 2009] and for SPECT in [Kazantsev et al., 2012, Dewaraja et al., 2010].For the model fusion of DOT with XCT, a similarity measure with aLaplacian-type smoothing operator within each region is proposed to average theDOT image within a region, respectively, while allowing discontinuity betweendifferent regions. [Brooksby et al., 2005,Brooksby et al., 2006,Yalavarthy et al., 2007,Baikejiang et al., 2017]. In [Fang et al., 2010, Boas, 2014, Deng et al., 2015], anempirical compositional relationship between the x-ray intensities for adipose andfibroglandular tissue is further proposed with the segmentation information fromXCT, and used to guide the reconstruction of DOT in [Fang et al., 2010, Boas, 2014,Deng et al., 2015]. ¶ This approach requires both accurate anatomical image segmentation and regis-tration [Fessler et al., 1992,Gindi et al., 1993,Ouyang et al., 1994,Sastry and Carson, 1997].
To resolve the difficulty in accurate anatomical image segmentation, informationtheoretic regularization using the mutual information or joint entropy tomeasure image similarity of image intensities are studied in [Nuyts, 2007,Tang and Rahmim, 2009, Panagiotou et al., 2009, Somayajula et al., 2011]. MRI isused as the information theoretic anatomical prior for the reconstruction of PETin [Somayajula et al., 2011] and DOT in [Panagiotou et al., 2009]. It is found thatthe joint entropy is more robust than the mutual information to differences betweenthe image under reconstruction and the anatomical image provided, especiallywhen the anatomical image has structures that are not present in the image ¶ To be precise, it is the x-ray tomosynthesis, a tomographic technique of limited views, that is usedin [Fang et al., 2010,Boas, 2014,Deng et al., 2015], instead of the conventional XCT with full-scannedcompleted data. oint bi-modal image reconstruction
Arandom spatial reordering of corresponding pairs of voxels in the two images wouldproduce identical measures ” [Somayajula et al., 2011]. Therefore, It is necessaryto incorporate features that capture the local variation in image intensities. In[Somayajula et al., 2011], scale-space features, such as a Gaussian-blurred imageand its Laplacian, are added into the similarity measure to help the joint entropyreconstruction. This method is related to the following derivative-based approach.Another disadvantage of information theoretic regularization is that it isprone to different image contrast ranges from different modalities and may leadto undesired reconstruction results [Panagiotou et al., 2009, Somayajula et al., 2011,Weizman et al., 2016].
As mentioned in the introduction, unlike model fusion, the structural informationin a JmmIR process is not from an image already reconstructed, but is updatedprogressively during the joint reconstruction process. Most of the model fusionmethods in § § L -normof second order derivatives of image differences is used as the similarity measurefor the joint inversion of gravity and seismic tomography data. A generalizedcross-gradient procedure is developed for joint multiple parameter inversion in[Gallardo, 2007, Gallardo and Meju, 2011], where cross-products of imaging gradientsare utilized to measure the parallelism of gradients, a.k.a., the structural similarity ofphysical parameters.The total variation is one of the most widely used regularizations in imagereconstruction problems [Rudin et al., 1992]. For JmmIR, the joint total variation bysumming up the total variations of images is used to solve the joint inversion problemin [Haber and Gazit, 2013]. However, such a joint total variation encourages the jointsparsity but not the similarity of image gradients. In [Rigie and La Riviere, 2015],the total nuclear variation (TNV) is proposed for encouraging common edge locationsand alignment of their gradient vectors, and is applied to the multi-spectral CTreconstruction. In [Knoll et al., 2017], the second order total generalized variation(TGV) with nuclear-norm is applied for the joint image reconstruction of MRI andPET. It is found that the TGV nuclear norm is robust regarding unwanted transfer oint bi-modal image reconstruction u ) and MRI ( u ), anotherimage similarity measure for image gradients is proposed to measure the parallelismof image gradients as follows, called parallel level sets (PLS),PLS β ( u , u ) = Z Ω ϕ (cid:2) ψ ( k∇ u k β · k∇ u k β ) − ψ (cid:0) |h∇ u , ∇ u i| β (cid:1)(cid:3) (1)where Ω ⊂ R N is a domain that contains the object to be imaged, for N = 2 or 3, k z k β = ( k z k + β ) for some β > β , h· , ·i is theEuclidean inner product, φ and ψ are strictly increasing functions. Please note thatthis image similarity measure is symmetric in u and u . Please note that asymmetricPLS can be formulated along the same framework [Ehrhardt, 2015].However, the above derivative-based methods depend on the magnitude ofgradients. For JmmIR, images tend to have different ranges of contrasts anddifferent heights and signs of edges detected with derivatives [Rasch et al., 2017].Normalized PLS can be formulated to resolve this issue within the frameworkin [Ehrhardt, 2015]. Another approach to resolve this issue is using the infimalconvolution of generalized Bregman distances to formulate weighted similarity measurefor image gradients without the edge orientation constraint [Rasch et al., 2017], whereasymmetric weightings are also possible. This approach has also the advantage forjoint reconstruction to avoid the artificial transfer of non-shared structures betweenthe images [Rasch et al., 2017]. All work reviewed above are valuable contributions because they demonstrate theperformance and advantages of JbmIR over SmIR in a number of multi-modal imagingapplications. All the methods in this subsection use image gradients or derivatives asthe image structural representation. Similarity measures are proposed to measurethe parallelism or correlation of image gradients or derivatives, respectively. In spiteof their demonstrated success, none of them is working directly on image edge. Forimage reconstruction, feature reconstruction [Louis, 2011] provides the possibility fordirectly addressing structural information on features like image edges, but has notyet been used for joint multi-modal image reconstruction. It is our motivation todevelop an approach directly with image edge together with an appropriate similaritymeasure, which is expect to be more robust to data noise than the derivative-basedmethods.As discussed in [Haber and Gazit, 2013], the challenges in JmmIR are as follows.The first fundamental challenge is because of the difference of contrast ranges ordifferent scales of different imaging modalities [Rasch et al., 2017], or different physicalmeaning (and units) [Haber and Gazit, 2013]. Because structural similarity is the onlyclue for JmmIR, the second one is an appropriate representation of image structuralinformation, which should support to capture the local variation in image intensities,be independent of image contrasts, and be updated progressively during the jointreconstruction process. Consequently the third one is an appropriate similaritymeasure for the representation of image structural information. There could befeatures in one modality but nonexistent in another. Hence, an appropriate similaritymeasure should encourage the reconstruction of features only if there is sufficient oint bi-modal image reconstruction
3. Similarity measures for image edges
Many similarity measures explain similarity (or dissimilarity) as a distance d ( · , · ) ina metric space X of perceptual stimulus of objects, with d satisfying the followingmetric axioms [Tversky, 1977],(i) Minimality: d ( K , K ) ≥ d ( K , K ) = 0.(ii) Symmetry: d ( K , K ) = d ( K , K ).(iii) Triangle inequality: d ( K , K ) + d ( K , K ) ≥ d ( K , K ),where K , K and K ∈ X . This is called the geometric model ofsimilarity [Blough, 2001, Rorissa, 2007]. However, all three metric axioms areunnecessarily stringent and can be inconsistent with a number of psychologicalexperiments [Tversky, 1977, Mumford, 1991, Rorissa, 2007, Skopal and Bustos, 2011,Biasotti et al., 2016]. It is realized that “Similarity isn’t a metric anyway” in[Mumford, 1991] after knowing the work of [Tversky, 1977].In his prominent paper [Tversky, 1977], Tversky proposed his famous featurecontrast model for the similarity of binary features. Let ∆ = { a, b, c, · · · } be thedomain of objects under study. Assume that each object in ∆ is represented by a setof binary features. Unlike the geometric model of similarity which represent stimulusas points in a metric space, an object a is characterized by the set A of features that theobject a possesses. Features are binary in the sense that a given feature of an objecteither is or is not in its set of features A . An example of the set of features is theimage edge set of an image. Tversky proposed another set of axioms about similaritymeasures of binary features, which include the axioms of matching, independence,solvability, invariance, and proved mathematically that feature similarity measuresmust be of the following form [Tversky, 1977] + S ( a, b ) = p ( A ∩ B ) − γ p ( A \ B ) − γ p ( B \ A ) , (2)where A and B denote the sets of binary features associated with the objects a and b , respectively, γ and γ are nonnegative constants. p is an additive functionsuch that p ( A ∩ B ) = p ( A ) + p ( B ) whenever A ∩ B = ∅ . This representation iscalled the feature contrast model and has been validated with multiple psychologicalexperimental data sets [Tversky, 1977, Tversky and Gati, 1982] and applied in image + Interested readers are referred the paper [Tversky, 1977] for details. Although the axiomsin [Tversky, 1977] are translated in some papers, we found that the original presentation in[Tversky, 1977] is still the best. Please note that some translations are incomplete due to thecomplicated formulations of axioms, especially for the axioms of solvability and invariance. It isfor the same reason that we give up the translation of Tversky’s axioms in this paper. oint bi-modal image reconstruction A to B , the biggerthe similarity measure S ( A, B ). This formula (2) implies that a reasonable featuresimilarity measure is a linear combination of the measures of their common and theirdistinct features.Similarity measures thus obtained increase with addition of common featuresand/or deletion of distinctive features (i.e., features that belong to one object but notto the other) [Tversky, 1977]. They posit that two stimuli are more similar if they havemore common features and fewer unique/distinctive features [Rorissa, 2007]. Hence,similarity measures will not force the nonexistent features to be reconstructed whenapplied to JmmIR as discussed in § γ = γ . It should address thatasymmetric similarity measures are in general necessary for JmmIR. Because ofdifferent image contrasts among modalities, similarity measures should be asymmetricat least from a heuristic perspective. Some of the approaches in § §
2. It is a segmentation-freemethod with a smoothing Markov prior with weights determined by the MRI imageintensities at nearby voxels. It can be used in an asymmetrical way by using onlypartial derivatives in certain directions [Bowsher et al., 2004, Schramm et al., 2018].It is reported in [Schramm et al., 2018] that the PLS is slightly inferior compared tothe asymmetrical Bowsher method.For 2D images, a natural choice for the measure p of image edge sets is the 1-dimensional Hausdorff measure H , or the length of image edges. Therefore, we obtain, S ( u , u ) = H ( K ∩ K ) − γ H ( K \ K ) − γ H ( K \ K ) , (3)where K and K denote the image edge sets of images u and u , respectively.This similarity measure of image edges is closely related to the Mumford-Shahregularization functional [Mumford and Shah, 1989]. Higher-dimensional Hausdorffmeasures can also be used in (3) for higher-dimensional images [Jiang et al., 2014].
4. JbmIR based on image edges
In this section, we apply the similarity measure for image edges of the previous sectionto JbmIR and establish a general approach for JbmIR based on image edges. Withoutloss of generality, we restrict to the bi-modal case and 2D case in this paper to avoidnotational complexity, although the mathematical formulation can be extended togeneral multi-modal case.For multi-modality imaging systems, the image reconstruction problem is toestimate images u i from measurement data g i , for i = 1 , A i ( u i ) = g i , (4)where A i is the forward operator of underlying imaging modality. TheMumford-Shah functional was first proposed as a variational approach for oint bi-modal image reconstruction u i , K i ) = Z Θ i | A i ( u i ) − g i | + α i Z Ω \ K |∇ u i | + β i H ( K i ) , (5)where Θ i is the domain of the measurement, K i is the edge set of the image u i , α i and β i are positive regularization parameters, for i = 1 , E ( u , u , K , K ) = MS( u , K ) + MS( u , K ) − τ S ( u , u ) , (6)where τ > − ” is toensure the similarity of edge sets of images u and u during the joint reconstructionprocess because of the convention that the more similar edges, the bigger theirsimilarity. Let Φ i ( u i , K i ) = Z Θ i | A i ( u i ) − g i | + α i Z Ω \ K i |∇ u i | . (7)Then we find, E ( u , u , K , K )= Φ ( u , K ) + Φ ( u , K ) + β H ( K ) + β H ( K ) − τ H ( K ∩ K ) + τ γ H ( K \ K ) + τ γ H ( K \ K ) . (8)In the above equation, the terms involving only the edges K and K sum up to β H ( K ) + β H ( K ) − τ H ( K ∩ K ) + τ γ H ( K \ K ) + τ γ H ( K \ K ) (9)= β H ( K ∩ K ) + β H ( K \ K ) + β H ( K ∩ K ) + β H ( K \ K ) (10) − τ H ( K ∩ K ) + τ γ H ( K \ K ) + τ γ H ( K \ K ) (11)=( β + β − τ ) H ( K ∩ K ) + ( β + τ γ ) H ( K \ K ) + ( β + τ γ ) H ( K \ K ) . (12)By re-parametrization with the reuse of the same notations for the regularizationparameters β , β , γ , and γ to avoid notational complexity, this sum can be written oint bi-modal image reconstruction ∗ Ψ( K , K )= β [ H ( K ∩ K ) + γ H ( K \ K )] + β [ H ( K ∩ K ) + γ H ( K \ K )] . (16)In summary, the joint reconstruction functional can be written as E ( u , u , K , K )= X i =1 ,j = i Z Θ i | A i ( u i ) − g i | + α i Z Ω \ K i |∇ u i | + β i [ H ( K i ∩ K j ) + γ i H ( K i \ K j )] , (17)where α i , β i and γ i ( i = 1 ,
2) are positive regularization parameters.
5. Variational approximation for JbmIR
In [Page, 2015], an extended Mumford-Shah functional with a priori edge informationwas established for image reconstruction from the perspective of regularization theory.For an image reconstruction problem A ( u ) = g with similar notations as in (4) and (5) ,if partial information of image edge can be obtained, then the following reconstructionfunctional is proposed in [Page, 2015],MS K ( u, K )= Z Θ | A ( u ) − g | + α Z Ω \ K |∇ u | + β (cid:2) H ( K \ K ) + γ H ( K ∩ K ) (cid:3) , (18)where K is the partial information of image edge. Theoretical results and numericalexperiments were reported in [Page, 2015]. The advantages of the extended Mumford-Shah functional were demonstrated by numerical experiments with the model fusionof DOT from XCT.By ignoring the accurate values of parameters, we have the joint reconstructionfunctional in (17) E ( u , u , K , K ) = MS K ( u , K ) + MS K ( u , K ) . (19)However, because none of K and K is fixed, it is difficult to extend the theoreticalresults in [Page, 2015] on existence, regularization property, Γ-approximation to thejoint reconstruction functional. JbmIR based on image edge as formulated in (17)is to find the minimizer pair ( u , u , K , K ) given the measurement g and g and ∗ Let β ′ , β ′ , γ ′ , and γ ′ be positive parameters such that β + β − τ = β ′ + β ′ , (13) β + τγ = β ′ γ ′ , (14) β + τγ = β ′ γ . ′ (15)Here we assume that β + β > τ . This is a reasonable assumption because all the edge regularizationsin (12) reduce to the Mumford-Shah regularization on image edge with the regularization parameter β + β − τ , in the ideal case when the image edge sets of u and u perfectly match, i.e., K = K .The equation (16) holds with β , β , γ , and γ being replaced by β ′ , β ′ , γ ′ , and γ ′ , respectively. oint bi-modal image reconstruction K ( u, K ) in (18) is proposed,MS K ( u, v )= Z Θ i | A ( u ) − g | + α Z Ω v |∇ u | + β (cid:20)Z Ω (cid:18) ε |∇ v | + (1 − v ) ε (cid:19) (1 + γ ( v − v ) ) (cid:21) , (20)where 0 ≤ v ≤ ≤ v ≤ K and K , respectively, i.e., v ( x ) = 1 if x / ∈ K and v ( x ) = 0 if x ∈ K , with thesame definition for v and K . ε controls the “width” of edges. The Γ-convergence ofMS K ( u, v ) is only proved for γ = 0 in [Page, 2015]. Although the Γ-approximation in(20) is proposed from an heuristic perspective, its performance is demonstrated withnumerical experiments in [Page, 2015]. Based on the relation in (19), we propose thefollowing Γ-approximation to the joint reconstruction functional E in (17) as follows, F ( u , v , u , v )= X i =1 ,j = i Z Θ i | A i ( u i ) − g i | + α i Z Ω v i |∇ u i | + β i (cid:20)Z Ω (cid:18) ε i |∇ v i | + (1 − v i ) ε i (cid:19) (1 + γ i ( v j − v i ) ) (cid:21) , (21)where 0 ≤ v ≤ ≤ v ≤ K and K , respectively. ε i controls the “width” of edges. As ε → ε → v and v are expected toconverge to the characteristic functions of the edge sets K and K , respectively, inthe sense of L . The current formulation of the Γ-approximation (21) is based on theheuristics in [Page, 2015, § γ = γ = 0 might be possible by using the samemethod in [Page, 2015] where a lengthy proof for single modal image reconstructionwith the extended Mumford-Shah function with a priori edge information is providedfor γ = 0. In this work we concentrate in studying numerically how the proposedmethod can help to improve the quality of reconstructed images.
6. JbmIR of XCT and DOT
In this section we give a brief introduction of XCT and DOT, and apply the JbmIRreconstruction method in the previous section to the JbmIR of XCT and DOT.
In conventional X-ray tomography, the image contrast is due to the x-ray absorptionwhen X-ray beams pass through an object. The interaction of X-ray and the object oint bi-modal image reconstruction I ( x ) be the intensity of an X-ray and f ( x ) the X-ray attenuation coefficientat position x . Then along a straight line L , dIdl = − f dl, (22)where dl is the differential element along L . Let I be the initial intensity of the X-rayand I the intensity after passing through the object. As we assume that the beamtravels in the straight line L , from (22) it followsln I I = − Z L f dl. (23)The measurement I and the initial intensity I provide line integrals of the X-rayattenuation coefficient f along lines, if the measurement is conducted in multipledirections and positions. The operator that maps a function into the set of its lineintegrals is the Radon Transform, which is defined as the following in the case of 2DXCT, R ( f )( θ, s ) = Z L f dl, (24)for every line L crossing the domain Ω ⊂ R containing underlying objects, where θ =(cos ϕ, sin ϕ ) is the normal direction of the line L = { ( x, y ) ∈ R : x cos ϕ + y sin ϕ = s } and s ∈ R . Let Θ = S × R (where S is the unit circle), and u = f , and g = ln I I .Then the XCT is to estimate the image u from its measurement data g , R ( u ) = g . (25) The steady-state diffuse optical tomography (DOT) is a functional imaging modalitywhich uses near infrared light to illuminate biological tissue and reconstruct theinside distribution of optical properties using the light intensity measured on thesurface of the tissue. Light at near-infrared wavelength will be absorbed andscattered while passing through biological tissue. The optical properties have directbiological relevance since the absorption of light is due to the existence of oxy-/deoxy-haemoglobin. Please refer to [Arridge, 1999, Arridge and Schotland, 2009] for moredetails.The governing equation used for diffuse light u in DOT is the diffuseapproximation of the radiative transport equation (RTE), − div( D ∇ u ) + µ a u = 0 , in Ω (26)where D is the diffusion coefficient and µ a is the absorption coefficient. The incominglight q can be modeled by Robin boundary condition in the diffuse approximation, u + 2 D −→ n · ∇ u = q, on ∂ Ω (27)where −→ n is the outer normal on the boundary ∂ Ω of Ω. The measurement g is thenegative Neumann boundary values of the solution of equation (26) g = − D −→ n · ∇ u, on ∂ Ω . (28) oint bi-modal image reconstruction q , diffusion coefficient D and absorption coefficient µ a , the forwardoperator F in DOT maps ( µ a , D ) to the measurement data g in (28). The DOT isto estimate ( µ a , D ) from the measurement data g , F ( µ a , D ) = g . (29)In this work we consider a special case of DOT to demonstrate the performance ofthe proposed JbmIR method. We are interested in recovering the absorption coefficient µ a , and assume the diffusion coefficient D to be known in the following. This is aninteresting but non-linear and ill-posed inverse problem [Arridge and Schotland, 2009].For a known diffusion coefficient D , we introduce the new forward operator G as G ( µ a ) = F ( µ a , D ) (30)To use the same notations in the previous sections, we write u to represent theabsorption image µ a , and g for the measurement G ( µ a ). Then the DOT with knowdiffuse coefficient is to estimate u from the measurement data g , G ( u ) = g , (31)where the measurement is on the boundary Θ = ∂ Ω of Ω.
With the development of hybrid multi-modal imaging techniques, there are a numberof hybrid optical tomography systems reported in the literature together withreconstruction methods by model fusion from XCT or JbmIR with XCT, such asDOT/XCT [Yuan et al., 2010, Boas, 2014, Deng et al., 2015, Baikejiang et al., 2017],BLT/DOT/XCT [Yang et al., 2015]. To reduce the radiation dose and accelerate thedata acquisition, it is the x-ray tomosynthesis technique of limited views that is usedin the latter work [Baikejiang et al., 2017, Boas, 2014, Deng et al., 2015] rather thanthe conventional full-scan XCT.By applying the general JbmIR formulation in (21), we obtain the following jointreconstruction function for XCT and DOT, F ( u , v , u , v )= Z Θ | R ( u ) − g | + α i Z Ω v |∇ u | + β (cid:20)Z Ω (cid:18) ε |∇ v | + (1 − v ) ε (cid:19) (1 + γ ( v − v ) ) (cid:21) + Z Θ | G ( u ) − g | + α Z Ω v |∇ u | + β (cid:20)Z Ω (cid:18) ε |∇ v | + (1 − v ) ε (cid:19) (1 + γ ( v − v ) ) (cid:21) , (32)where 0 ≤ v ≤ ≤ v ≤ u and DOT image u , respectively. Please note that there are noweighting parameters in front of the two fidelity terms, but this does not mean bothfidelity terms contribute equally to the joint reconstruction. Because of the alternatingminimization algorithm discussed in the following, they never meet each other andtheir contributions are independent of each other. The difference of their contribution oint bi-modal image reconstruction β , β , γ and γ . During the trainingof the regularization parameters, the two fidelity terms also play roles independent ofeach other. The difference of their contribution, even if there are different weights forboth, will merged into the regularization parameters.The JbmIR of XCT and DOT in this work is to reconstruct the images u and u as well as the edge indicators v and v simultaneously in an iterative process. This isachieved by minimizing the objective functional (32) with an alternating minimizationalgorithm in the order of u , v , u and v . Each minimization for u , v , u and v is a gradient descending with the backtracking line search method by the Armijo-Goldstein stopping condition [Armijo, 1966, Wikipedia contributors, 2017], with theprevious reconstruction as initial value and the others being frozen at the latestreconstructed values. Hence, the proposed alternating minimization algorithm forthe JbmIR of XCT and DOT is an instance of stochastic gradient descending with theaforementioned order, and is described in Algorithm 1. Notice that we use arg minin the algorithm description to indicate our minimizing purpose and in each step weuse gradient descent method, therefore the results we get from the iteration are theapproximations of minimizers. Algorithm 1
Alternating minimization for JbmIR of XCT and DOT for i = 0 , , , . . . , n do u i +11 = arg min u F ( u , v i , u i , v i ) (33)with u i as the initial value; v i +11 = arg min v F ( u i +11 , v , u i , v i ) (34)with v i as the initial value; u i +12 = arg min u F ( u i +11 , v i +11 , u , v i ) (35)with u i as the initial value; v i +12 = arg min v F ( u i +11 , v i +11 , u i +12 , v ) (36)with v i as the initial value; end for Although there are many terms in the joint reconstruction functional F in (32),there is not so many computing terms at each minimization step in (33), (34) (35) and(36) because some terms can be treated as constants. For example, at the minimizationstep (33), only the first two terms at the first line on the right-hand side of (32) areinvolved in this minimization step, those terms not involving u can be treated asconstants and ignored. Especially the time consuming computation of the forwardprocess G for the DOT is not needed. We have carefully analyzed the dependentterms at each minimization step to reduce the computing load. Implementation andexperimental details are reported in the next section. oint bi-modal image reconstruction
7. Numerical experiments
In discrete case, the gradients are discretized with forward difference and thedivergence with backward difference, so that the discretized gradient corresponds tothe gradient of the discretized functional.For XCT, the initial image value is set as zero and the initial edge value is setas one. For DOT, zero and one initials cannot lead to good reconstruction results.We use gradient descent method for only the fidelity term to get a blurry image asthe initial image value, and use the Mumford-Shah functional segmenting this initialimage to get the initial edge value. The SmIR-MS and JbmIR-MS both use this sameinitial value setting.We also need to set up the stopping criterion of our algorithm. This is achieved bysetting up the iteration numbers involved. After a number of numerical experimentswith phantoms of similar structures as in this work, we find that after 80 wholeiterations with 10 steps for each variables in every iteration the SSIM values, whichwe use to evaluate the quality of the reconstructed images, of both XCT and DOTstabilize at relatively higher values, and in our examples over 0.85 for XCT and over0.69 for DOT.We evaluate the quality of the reconstructed image u rec via the structuralsimilarity index measure (SSIM) [Wang et al., 2004] using the corresponding trueimage u true as a reference image:SSIM( u rec , u true ) = (2 u true u rec + C )(2 σ u true u rec + C )( u true + u rec + C )( σ u true + σ u rec + C ) (37)where u true , u rec are averages, σ u true , σ u rec are variances and σ u true u rec is thecovariance.SSIM value is proposed to measure the similarity between two images consideringthree comparison measurements: luminance, contrast and structure. We use thedefault parameter settings that C = (0 . L ) , C = (0 . L ) , where L is the dynamicrange of the ground truth. The higher SSIM value is the more similar the two imagesare. SSIM=1 only happens for identical images. The computations are done using Matlab. For XCT, the Radon Transformation andthe adjoint operator written by Lut Justen from the Software-Documentation of theCenter for Industrial Mathematics, University of Bremen are implemented. For DOT,the TOAST++ package [Schweiger and Arridge, 2014] from Martin Schweiger andSimon Arridge is implemented for the forward operator and the Jacobian matrix ofthe forward operator.For XCT, the projection data are obtained from 30 views uniformly distributed in[0 , o ]. Because our phantoms are of radius 25 mm and image sizes are of 100 × n to our noise free data d . Since the fidelityterm R Θ i | A i ( u i ) − g i | shows that A i ( u i ) should be a good approximation of the data oint bi-modal image reconstruction g i in L norm, we scale the noise n to the data with regard to the L norm. Therelative noise level is defined as η = k n kk d k , where k · k is the Euclidean norm. In ournumerical experiments 5% relative Gaussian noise is added to the XCT data and 2%relative Gaussian noise is added to the DOT data. −20 −10 0 10 20−25−20−15−10−50510152025 Figure 1: Finite element mesh for DOT problem with 16 source positions (red points)and 16 detector positions (blue points).
The choice of the regularization parameters is non-trivial and should be done task-dependently [Schramm et al., 2018].In our proposed JbmIR-MS method, the “width” control parameters ε and ε areset as 1 × − , which are the same as in [Page, 2015]. There are three regularizationparameters α i , β i and γ i ( i = 1 ,
2) needed to be chosen properly for each modality. α i is the regularization parameter of the smoothing penalty term. The larger it is, thesmoother the reconstructed image will be when there is not any edge. β i weights howstrongly the edge penalizes. The larger it is, the fewer edges will be reconstructed. γ i controls the effect of the reference edge information. The effect of the reference edgeinformation increases as the value of γ i increases.For JbmIR-MS, from the XCT with a priori edge information numericalexperiments done by Thomas we can get that when the order of magnitude of α is around 10 , the order of magnitude of β is around 10 − or 10 − and γ is nomore than 10 we could be able to get XCT results with clear structures. We firstfix the XCT parameters as α = 9 × , β = 2 × − , γ = 10 to choose DOTparameters. For joint DOT reconstruction we first roughly partition the parameterspace and compute the SSIM value to find rough-optimal parameters. Then accordingto the reconstruction results and the parameter effect discussed above, we fine-tunethe parameters around this rough-optimal choice. For example, if the reconstructedimage is too smooth in non-edge area, we will decrease the value of α . If thereare too few edges, we will slightly decrease the value of β . If the influence of thereference edge information is too strong that we reconstruct wrong structures thatshould only appear in one modality and not in the other, we will slightly decrease the oint bi-modal image reconstruction γ . Here α is chosen from A = [0 . , . , . . . , × ], β is chosen from B = [1 × − , × − , . . . , . , γ is chosen from C = [1 , , . . . , , α = 1 × , β = 1 × − and γ = 5 produce a higher SSIM value of the reconstructed DOT image againstits true image, among all possible combinations from A × B × C . Then we fine-tunethis rough-optimal choice according to the discussion above and reach the parametersettings we use in our numerical experiments. After choosing the DOT parameters, wethen fine-tune the joint XCT parameters around α = 9 × , β = 2 × − , γ = 10.The final parameter settings for every experiment are listed in the description of thenumerical results.We compare our method with SmIR using Mumford-Shah functional(SmIR-MS).For SmIR-MS, the “width” control parameter ε is also set as 1 × − . There aretwo regularization parameters α and β for each modality. The reconstruction result issensitive to both of them. The larger α is the smoother the reconstructed image willbe when there is not any edge. The larger β is the fewer edges will be reconstructed.We first perform a search of the parameter space in a rough scale with the SSIMas the criterion. For both single modality reconstructions α and β are chosen from A = [0 . , . , . . . , × ] and B = [1 × − , × − , . . . , . ,
1] to determine theproper order of magnitudes first. For DOT reconstruction we find that α = 1 × and β = 1 × − produce a higher SSIM value of the reconstructed DOT image against itstrue image, among all possible combinations from A × B . For XCT reconstruction theparameters are α = 1 × and β = 1 × − . We then fine-tune this rough-optimalchoice according to the discussion above and our final parameter settings for everyexperiment are listed in the description of the numerical results. We design three pairs of phantoms of size 100 × u and the absorption coefficients u are listedin Table 1, Table 2 and Table 3. True u True u Figure 2: Phantom 1. Left: the true distribution of the X-ray attenuation coefficient u . Right: the true distribution of the absorption coefficient u . oint bi-modal image reconstruction u ( mm − ) 1 2 — 2.5DOT u ( mm − ) 0.01 0.03 0.03 — True u True u Figure 3: Phantom 2. Left: the true distribution of the X-ray attenuation coefficient u . Right: the true distribution of the absorption coefficient u .Table 2: Phantom 2Background Four circles small circleXCT u ( mm − ) 1 2 2.5DOT u ( mm − ) 0.01 0.03 — True u True u Figure 4: Phantom 3. Left: the true distribution of the X-ray attenuation coefficient u . Right: the true distribution of the absorption coefficient u .Table 3: Phantom 3Background Two circle and one ellipse Top circle small circleXCT u ( mm − ) 1 2 — 2.5DOT u ( mm − ) 0.01 0.03 0.03 — oint bi-modal image reconstruction To compare the results, we put the reconstructed images of JbmIR-MS and SmIR-MS together with the true distribution, and we adjust the display window all thesame as the ground truth, that is [0 , .
5] for XCT images and [0 , .
03] for DOTimages. Since Mumford-Shah functional provides edge information, we also comparethe reconstructed edge information using JbmIR-MS and SmIR-MS.
For JbmIR-MS we choose α = 8 . × , β = 1 . × − , γ =9 . , α = 1 × , β = 6 × − , γ = 5. For SmIR-MS XCT, we choose α = 8 . × , β = 8 × − . For SmIR-MS DOT, we choose α = 1 × , β = 5 × − . Figure 5 shows the reconstruction results. The three images in the first row from leftto right are the true distribution of the X-ray attenuation coefficient u , JbmIR-MS u and SmIR-MS u . The second row are the true distribution of the absorptioncoefficient u , JbmIR-MS u , and SmIR-MS u . According to the SSIM values,the image quality of reconstructed u is improved using JbmIR-MS (SSIM from0.79 to 0.87). The image quality of reconstructed u is improved using JbmIR-MS(SSIM= 0 .
69) then SmIR-MS (SSIM= 0 . True u SSIM=0.69 00.010.020.03 SmIR-MS u SSIM=0.61 00.010.020.03True u SSIM=0.87 00.511.522.5 SmIR-MS u SSIM=0.79 00.511.522.5
Figure 5: Phantom 1 reconstruction results. (top) The XCT phantom andreconstruction results, (bottom) The DOT phantom and reconstruction results. ForSmIR-MS u we have SSIM= 0 .
79, and for the JbmIR-MS u we have SSIM= 0 . u , we have SSIM= 0 .
61, for JbmIR-MS u we have SSIM= 0 . u from bi-modalityand single-modality. The second row are reconstructed edge images for absorptioncoefficient u from JbmIR-MS and SmIR-MS. We can find that JbmIR-MS provides oint bi-modal image reconstruction JbmIR-MS u edges SmIR-MS u edgesJbmIR-MS u edges SmIR-MS u edges Figure 6: Reconstructed edge images.
For JbmIR-MS we choose α = 8 . × , β = 1 . × − , γ =9 . , α = 1 × , β = 7 × − , γ = 5. For SmIR-MS XCT, we set α = 8 . × , β =8 × − . For SmIR-MS DOT, we set α = 1 × , β = 5 × − . The reconstructedresults are shown in Figure 7. The image quality of reconstructed u is improved usingJbmIR-MS (SSIM from 0.75 to 0.86). The image quality of reconstructed u is higherusing JbmIR-MS (SSIM= 0 .
70) then SmIR-MS (SSIM= 0 . For JbmIR-MS we choose α = 8 . × , β = 1 . × − , γ =9 . , α = 1 × , β = 6 × − , γ = 5. For SmIR-MS XCT, we choose α = 8 . × , β = 8 × − . For SmIR-MS DOT, we choose α = 1 × , β = 5 × − . The reconstruction results are shown in Figure 9. According to the SSIM values, theimage quality of reconstructed u is improved using JbmIR-MS (SSIM from 0.76 to0.87). The image quality of reconstructed u is higher using JbmIR-MS (SSIM= 0 . . oint bi-modal image reconstruction True u SSIM=0.70 00.010.020.03 SmIR-MS u SSIM=0.61 00.010.020.03True u SSIM=0.86 00.511.522.5 SmIR-MS u SSIM=0.75 00.511.522.5
Figure 7: Phantom 2 reconstruction results. (top) The XCT phantom andreconstruction results, (bottom) The DOT phantom and reconstruction results. ForSmIR-MS u we have SSIM= 0 .
75, and for JbmIR-MS u we have SSIM= 0 .
86, ForSmIR-MS u we have SSIM= 0 .
61, for JbmIR-MS u we have SSIM= 0 . JbmIR-MS u edges SmIR-MS u edgesJbmIR-MS u edges SmIR-MS u edges Figure 8: Reconstructed edge images oint bi-modal image reconstruction True u SSIM=0.69 00.010.020.03 SmIR-MS u SSIM=0.60 00.010.020.03True u SSIM=0.87 00.511.522.5 SmIR-MS u SSIM=0.76 00.511.522.5
Figure 9: Phantom 3 reconstruction results. (top) The XCT phantom andreconstruction results, (bottom) The DOT phantom and reconstruction results. ForSmIR-MS u we have SSIM= 0 .
76, and for JbmIR-MS u we have SSIM= 0 .
87, ForSmIR-MS u we have SSIM= 0 .
60, for JbmIR-MS u we have SSIM= 0 . JbmIR-MS u edges SmIR-MS u edgesJbmIR-MS u edges SmIR-MS u edges Figure 10: Reconstructed edge images oint bi-modal image reconstruction
8. Discussion
In this work, we establish an image similarity measure in terms of image edges fromTversky’s theory of feature similarity in psychology to solve the DOT and XCT jointreconstruction problem. This similarity consists of Hausdorff measures of the commonand different parts of image edges from both modalities. When applied to JbmIR it willimprove the reconstructed image quality and will not force the nonexistent structuresto be reconstructed.To implement the joint reconstruction, a Γ − approximation(21) is proposed forthe joint reconstruction functional. The current formulation of the Γ-approximationis based on the heuristics in [Page, 2015, § γ = γ = 0 might be possibleby using the same method in [Page, 2015] where a lengthy proof for single modalimage reconstruction with the extended Mumford-Shah function with a priori edgeinformation is provided for γ = 0. However, it is too tedious to reproduce such a proofin the current work and does not contribute any mathematical insight into this problemwith such a stringent condition in our opinion. Instead, we concentrate in studyingnumerically how the proposed method can help improve the quality of reconstructedimages in this work. In the Γ-approximation functional there isn’t any weightingparameters in front of the two fidelity terms, that is because they are independentof each other. We reconstruct the two modalities alternatively, and the informationcommunication is achieved by edge similarity, so we just choose proper regularizationparameters for penalty term to control the effect of the other modality and do notweight the fidelity term.Since circle is the simplest shape, we first design phantoms with circle internalstructure. For each example, the structure of the two phantoms are designed similarto each other. XCT should be able to reconstruct finer structure so we add a smallcircle to the XCT phantoms. DOT cannot distinguish such fine structure so we onlydesign big circle inner structure. Considering that DOT can distinguish different softtissue better than XCT, we add one more big circle in DOT phantom for Example1 and Example 3. A reasonable similarity measure should maintain the distinctivestructure in each modality and not force the nonexistent features in one modality tobe reconstructed. In Example 2, the main parts are set exactly the same as each otherto see if we can improve image quality more when the phantoms from two modalitiesshare more common structure. And in Example 3, we change the shape of the middlestructure to be an ellipse to see the performance of different methods for differentshape of the inner structure.Our JbmIR-MS is compared with SmIR-MS. From results shown in Figure 5 forExample 1, we can see that the quality of images reconstructed jointly is higher thanthat of images reconstructed under the SmIR-MS. The SSIM values increase 10 . .
78% for DOT and XCT respectively. Looking at the reconstructed images,JbmIR-MS u keeps the structure of the small circle and no extra top circle structureis reconstructed because of the information provided by u . As for JbmIR-MS u thestructure is more clear than SmIR-MS result especially the parts that are consistentwith u . It is more obvious in the reconstructed edge information shown in Figure 6.The results for Example 2 are shown in Figure 7. We can see that the imagequality improvement is greater than that in Example 1. The SSIM values increase14 .
85% and 14 .
44% for DOT and XCT respectively. We think that the reason thatthe enhancement in Example 2 is greater than that in Example 1 is because the oint bi-modal image reconstruction a priori structural information. Compare to the edge results in Example 1, the edgeof the top circle for DOT image is much sharper because the XCT phantom in thisexample can provide the edge information of the top circle.In Example 3 we can see that for different shapes of the internal structure,ellipse in our case, JbmIR-MS method also provides better reconstructed images,shown in Figure 9. The SSIM values increase 13 .
82% and 15 .
91% for DOT and XCTrespectively. And from the result figure we can clearly see the shape of the middleellipse in our joint reconstruction result for DOT problem.We find that the visual improvement of XCT reconstruction is not impressive,because the phantoms are easy for XCT. However, looking at their line profiles, theimprovement is easy to capture as shown in Figures 11(a), 11(b) and 11(c). It can beseen that the XCT results using JbmIR-MS are closer to the true distributions thanthat using SmIR-MS.From the numerical results we can see that minimizing the Mumford-Shahfunctional enables us to reconstruct images and obtain edges at the same time.Since the images from two modalities share similar structure, our JbmIR-MS methodadding the edge similarity measure to the Mumford-Shah functional implementsthe information communication between DOT and XCT modalities which helps toenhance the reconstructed image quality for both modalities. With proper choiceof regularization parameters, our similarity measure is able to help enhancing theimage quality of the common features, to remain distinctive features for each modalityand not force to reconstruct the nonexistent features. The smoothing term togetherwith the usage of separate edge variables to represent structure information enable usto get piecewise smooth reconstructed images with clear and sharp edge with noisedata. However there are a number of parameters, including regularization parameters,needed to be chosen carefully and the current parameter choosing is a trial-errorapproach. For 3D situation, the DOT/XCT bi-modal imaging system has alreadybeen set up. We are now preparing the optical phantoms [He et al., 2018] for realphysical experiment. And we are going to test the JbmIR-MS method with realisticdata acquired under our bi-modal imaging system.
9. Conclusion
In this paper, we have proposed an asymmetric edge similarity measure followingthe feature contrast model in [Tversky, 1977] to solve the DOT and XCT jointreconstruction problem. This edge similarity measure is able to implementthe information communication between the two modalities, thus improves thereconstructed image quality for both modalities. Meanwhile, it remains distinctivestructure for each modality and does not force the nonexistent edges to bereconstructed when some of the structures only exist in one modality.Representing the edge sets by edge indicator functions, a Γ-approximation (21) isproposed to make it easier to implement the joint reconstruction. To demonstrate theperformance of our JbmIR-MS algorithm, we have designed two dimensional numericalphantoms and calculate the minimizers of the Γ-approximation. To calculate theminimizers, we use gradient descent method with the negative gradient as the descentdirection and the step size chosen by backtracking line search with the Armijo-Goldstein stopping condition. oint bi-modal image reconstruction -25 -20 -15 -10 -5 0 5 10 15 20 250123 TrueJbmIR-MSSmIR-MS -25 -20 -15 -10 -5 0 5 10 15 20 250123
TrueJbmIR-MSSmIR-MS (a) XCT profiles for Phantom 1. -25 -20 -15 -10 -5 0 5 10 15 20 250123
TrueJbmIR-MSSmIR-MS -25 -20 -15 -10 -5 0 5 10 15 20 250123
TrueJbmIR-MSSmIR-MS (b) XCT profiles for Phantom 2. -25 -20 -15 -10 -5 0 5 10 15 20 250123
TrueJbmIR-MSSmIR-MS -25 -20 -15 -10 -5 0 5 10 15 20 250123
TrueJbmIR-MSSmIR-MS (c) XCT profiles for Phantom 3.
Figure 11: XCT profiles for the three phantoms and reconstructed images. Blackvertical lines on images on the left-side show the positions of lines to draw the lineprofiles on the right-side. For the line profiles, the blue line shows absorption valuesof the truth image, red circles shows absorption values from the SmIR-MS of XCT,and green stars shows absorption values from the JbmIR-MS of XCT. It can be seenthat the results from the JbmIR-MS are closer to the truth distributions than theSmIR-MS. oint bi-modal image reconstruction
Acknowledgments
This work was partially supported by the Sino-German Center (GZ 1025), NationalScience Foundation of China (61520106004, 61421062), and National Basic ResearchProgram of China (973 Program) (2015CB351803).
References [Armijo, 1966] Armijo, L. (1966). Minimization of functions having Lipschitz continuous first partialderivatives.
Pacific Journal of Mathematics , 16(1):1–3.[Arridge, 1999] Arridge, S. R. (1999). Optical tomography in medical imaging.
Inverse problems ,15(2):R41.[Arridge and Schotland, 2009] Arridge, S. R. and Schotland, J. C. (2009). Optical tomography:forward and inverse problems.
Inverse Problems , 25(12):123010.[Baikejiang et al., 2017] Baikejiang, R., Zhang, W., and Li, C. Q. (2017). Diffuse optical tomographyfor breast cancer imaging guided by computed tomography: A feasibility study.
Journal of X-Ray Science and Technology , 25(3):341–355.[Biasotti et al., 2016] Biasotti, S., Cerri, A., Bronstein, A., and Bronstein, M. (2016). Recent trends,applications, and perspectives in 3d shape similarity assessment.
Computer Graphics Forum ,35(6):87–119.[Blough, 2001] Blough, D. S. (2001). The perception of similarity. In Cook, R. B., editor,
AvianVisual Cognition . https://pigeon.psy.tufts.edu/avc/dblough/default.htm.[Boas, 2014] Boas, D. (2014). Digital breast tomosynthesis guided tomographic optical breastimaging (TOBI). Clinical trials 2014 - 2019, https://clinicaltrials.gov/ct2/show/NCT02033486.[Bockisch et al., 2009] Bockisch, A., Freudenberg, L. S., Schmidt, D., and Kuwert, T. (2009). Hybridimaging by SPECT/CT and PET/CT: Proven outcomes in cancer imaging.
Seminars inNuclear Medicine , 39(4):276–289.[Bowsher et al., 1996] Bowsher, J. E., Johnson, V. E., Turkington, T. G., Jaszczak, R. J., Floyd, C.,and Coleman, R. E. (1996). Bayesian reconstruction and use of anatomical a priori informationfor emission tomography.
IEEE Transactions on Medical Imaging , 15(5):673–686.[Bowsher et al., 2004] Bowsher, J. E., Yuan, H., Hedlund, L. W., Turkington, T. G., Akabani, G.,Badea, A., Kurylo, W. C., Wheeler, C. T., Cofer, G. P., Dewhirst, M. W., and Johnson, G. A.(2004). Utilizing MRI information to estimate F18-FDG distributions in rat flank tumors. InSeibert, J. A., editor, ,IEEE Nuclear Science Symposium - Conference Record, pages 2488–2492.[Brooksby et al., 2005] Brooksby, B., Jiang, S. D., Dehghani, H., Pogue, B. W., Paulsen, K. D.,Weaver, J., Kogel, C., and Poplack, S. P. (2005). Combining near-infrared tomographyresonance imaging to study in vivo and magnetic breast tissue: implementation of a laplacian-type regularization to incorporate magnetic resonance structure.
Journal of Biomedical Optics ,10(5):10.[Brooksby et al., 2006] Brooksby, B., Pogue, B. W., Jiang, S. D., Dehghani, H., Srinivasan, S., Kogel,C., Tosteson, T. D., Weaver, J., Poplack, S. P., and Paulsen, K. D. (2006). Imaging breastadipose and fibroglandular tissue molecular signatures by using hybrid mri-guided near-infraredspectral tomography.
Proceedings of the National Academy of Sciences of the United States ofAmerica , 103(23):8828–8833.[Canny, 1986] Canny, J. (1986). A computational approach to edge detection.
IEEE Transactionson Pattern Analysis and Machine Intelligence , 8(6):679 – 698.[Catana, 2017] Catana, C. (2017). Principles of simultaneous PET/MR imaging.
MagneticResonance Imaging Clinics of North America , 25(2):231–+. oint bi-modal image reconstruction [Censor et al., 2008] Censor, Y., Jiang, M., and Louis, A. K. (2008). Mathematical Methods inBiomedical Imaging and Intensity-modulated Radiation Therapy (IMRT) . Scuola NormaleSuperiore Edizioni della Normale.[Chan et al., 2009] Chan, C., Fulton, R., Feng, D. D., and Meikle, S. (2009). Regularized imagereconstruction with an anatomically adaptive prior for positron emission tomography.
Physicsin Medicine and Biology , 54(24):7379–7400.[Cherry, 2009] Cherry, S. R. (2009). Multimodality imaging: Beyond PET/CT and SPECT/CT.
Seminars in Nuclear Medicine , 39(5):348–353.[Delbeke et al., 2009] Delbeke, D., Schoder, H., Martin, W. H., and Wahl, R. L. (2009). Hybridimaging (SPECT/CT and PET/CT): Improving therapeutic decisions.
Seminars in NuclearMedicine , 39(5):308–340.[Deng et al., 2015] Deng, B., Brooks, D. H., Boas, D. A., Lundqvist, M., and Fang, Q. Q. (2015).Characterization of structural-prior guided optical tomography using realistic breast modelsderived from dual-energy x-ray mammography.
Biomedical Optics Express , 6(7):2366–2379.[Dewaraja et al., 2010] Dewaraja, Y. K., Koral, K. F., and Fessler, J. A. (2010). Regularizedreconstruction in quantitative SPECT using CT side information from hybrid imaging.
Physicsin Medicine and Biology , 55(9):2523–2539.[Douiri et al., 2007] Douiri, A., Schweiger, M., Riley, J., and Arridge, S. R. (2007). Anisotropicdiffusion regularization methods for diffuse optical tomography using edge prior information.
Measurement Science and Technology , 18(1):87–95.[Ehrhardt, 2015] Ehrhardt, M. J. (2015).
Joint Reconstruction for Multi-Modality Imaging withCommon Structure . PhD thesis, University College London.[Ehrhardt et al., 2015] Ehrhardt, M. J., Thielemans, K., Pizarro, L., Atkinson, D., Ourselin, S.,Hutton, B. F., and Arridge, S. R. (2015). Joint reconstruction of PET-MRI by exploitingstructural similarity.
Inverse Problems , 31(1):015001.[Even-Sapir et al., 2009] Even-Sapir, E., Keidar, Z., and Bar-Shalom, R. (2009). Hybrid imaging(SPECT/CT and PET/CT)-improving the diagnostic accuracy of functional/metabolic andanatomic imaging.
Seminars in Nuclear Medicine , 39(4):264–275.[Fang et al., 2010] Fang, Q. Q., Moore, R. H., Kopans, D. B., and Boas, D. A. (2010). Compositional-prior-guided image reconstruction algorithm for multi-modality imaging.
Biomedical OpticsExpress , 1(1):223–235.[Fessler et al., 1992] Fessler, J. A., Clinthorne, N. H., and Rogers, W. L. (1992). Regularized emissionimage-reconstruction using imperfect side information.
IEEE Transactions on Nuclear Science ,39(5):1464–1471.[Gallardo, 2007] Gallardo, L. A. (2007). Multiple cross-gradient joint inversion for geospectralimaging.
Geophysical Research Letters , 34(19).[Gallardo and Meju, 2011] Gallardo, L. A. and Meju, M. A. (2011). Structure-coupled multiphysicsimaging in geophysical sciences.
Reviews of Geophysics , 49(1).[Geman and Geman, 1984] Geman, S. and Geman, D. (1984). Stochastic relaxation, gibbsdistributions, and the bayesian restoration of images.
IEEE Transactions on Pattern Analysisand Machine Intelligence , 6:721–741.[Gindi et al., 1993] Gindi, G., Lee, M., Rangarajan, A., and Zubal, I. G. (1993). Bayesianreconstruction of functional images using anatomical information as priors.
IEEE Transactionson Medical Imaging , 12(4):670–680.[Haber and Gazit, 2013] Haber, E. and Gazit, M. H. (2013). Model fusion and joint inversion.
Surveys in Geophysics , 34(5):675–695.[Haber and Oldenburg, 1997] Haber, E. and Oldenburg, D. (1997). Joint inversion: a structuralapproach.
Inverse Problems , 13(1):63.[He et al., 2018] He, D., He, J., and Mao, H. (2018). Phantom preparation and optical propertydetermination.
Sensing and Imaging , 19(1):10.[Jiang et al., 2014] Jiang, M., Maass, P., and Page, T. (2014). Regularizing properties of theMumford–Shah functional for imaging applications.
Inverse Problems , 30(3):035007.[Judenhofer et al., 2008] Judenhofer, M. S., Wehrl, H. F., Newport, D. F., Catana, C., Siegel, S. B.,Becker, M., Thielscher, A., Kneilling, M., Lichy, M. P., Eichner, M., Klingel, K., Reischl, G.,Widmaier, S., Rocken, M., Nutt, R. E., Machulla, H. J., Uludag, K., Cherry, S. R., Claussen,C. D., and Pichler, B. J. (2008). Simultaneous PET-MRI: a new approach for functional andmorphological imaging.
Nature Medicine , 14(4):459–465.[Kaufmann and Di Carli, 2009] Kaufmann, P. A. and Di Carli, M. F. (2009). Hybrid SPECT/CTand PET/CT imaging: The next step in noninvasive cardiac imaging.
Seminars in NuclearMedicine , 39(5):341–347.[Kazantsev et al., 2012] Kazantsev, D., Simon, R. A., Pedemonte, S., Bousse, A., Erlandsson, K., oint bi-modal image reconstruction Hutton, B. F., and Ourselin, S. (2012). An anatomically driven anisotropic diffusion filteringmethod for 3D SPECT reconstruction.
Physics in Medicine & Biology , 57(12):3793.[Klann et al., 2011] Klann, E., Ramlau, R., and Ring, W. (2011). A Mumford-Shah level-setapproach for the inversion and segmentation of SPECT/CT data.
Inverse Problems & Imaging ,5(1):137–166.[Knoll et al., 2017] Knoll, F., Holler, M., Koesters, T., Otazo, R., Bredies, K., and Sodickson,D. K. (2017). Joint MR-PET reconstruction using a multi-channel image regularizer.
IEEEtransactions on medical imaging , 36(1):1–16.[Leahy and Yan, 1991] Leahy, R. and Yan, X. (1991). Incorporation of anatomical MR data forimproved functional imaging with PET. In
Information Processing in Medical imaging , pages105–120. Springer.[Louis, 1989] Louis, A. K. (1989).
Inverse und schlecht gestellte Probleme . Teubner.[Louis, 2011] Louis, A. K. (2011). Feature reconstruction in inverse problems.
Inverse Problems ,27(6):065010.[Marr, 1982] Marr, D. (1982).
Vision . Freeman.[Marr and Hildreth, 1980] Marr, D. and Hildreth, E. (1980). Theory of edge detection.
Proceedingsof the Royal Society of London. Series B , 207:187 – 217.[Mumford, 1991] Mumford, D. (1991). Mathematical theories of shape: do they model perception?
Proc.SPIE , 1570:1570 – 1570 – 9.[Mumford and Shah, 1989] Mumford, D. and Shah, J. (1989). Optimal approximations by piecewisesmooth functions and associated variational problems.
Communications on Pure and AppliedMathematics , 42(5):577–685.[Natterer, 2001] Natterer, F. (2001).
The Mathematics of Computerized Tomography . SIAM.[Natterer et al., 2002] Natterer, F., Wubbeling, F., and Wang, G. (2002). Mathematical methods inimage reconstruction.
Medical Physics , 29(1):107–108.[Nuyts, 2007] Nuyts, J. (2007). The use of mutual information and joint entropy for anatomicalpriors in emission tomography. In , IEEE NuclearScience Symposium, pages 4149–4154. IEEE Press.[Ouyang et al., 1994] Ouyang, X., Wong, W. H., Johnson, V. E., Hu, X., and Chen, C.-T. (1994).Incorporation of correlated structural images in PET image reconstruction.
IEEE Transactionson Medical Imaging , 13(4):627–640.[Page, 2015] Page, T. S. (2015).
Image reconstruction by Mumford-Shah regularization with a prioriedge information . PhD thesis, Universit¨at Bremen.[Panagiotou et al., 2009] Panagiotou, C., Somayajula, S., Gibson, A. P., Schweiger, M., Leahy, R. M.,and Arridge, S. R. (2009). Information theoretic regularization in diffuse optical tomography.
JOSA A , 26(5):1277–1290.[Perona and Malik, 1990] Perona, P. and Malik, J. (1990). Scale-space and edge detection usinganisotropic diffusion.
IEEE Transactions on Pattern Analysis and Machine Intelligence ,12(7):629 – 639.[Ramlau and Ring, 2007] Ramlau, R. and Ring, W. (2007). A Mumford–Shah level-set approach forthe inversion and segmentation of X-ray tomography data.
Journal of Computational Physics ,221(2):539–557.[Rasch et al., 2017] Rasch, J., Brinkmann, E.-M., and Burger, M. (2017). Joint reconstructionvia coupled bregman iterations with applications to PET-MR imaging.
Inverse Problems ,34(1):014001.[Rigie and La Riviere, 2015] Rigie, D. S. and La Riviere, P. J. (2015). Joint reconstruction of multi-channel, spectral ct data via constrained total nuclear variation minimization.
Physics inMedicine and Biology , 60(5):1741–1762.[Rondi, 2008] Rondi, L. (2008). Reconstruction in the inverse crack problem by variational methods.
European Journal of Applied Mathematics , 19(6):635–660.[Rondi and Santosa, 2001] Rondi, L. and Santosa, F. (2001). Enhanced electrical impedancetomography via the Mumford–Shah functional.
ESAIM: Control, Optimisation and Calculusof Variations , 6:517–538.[Rorissa, 2007] Rorissa, A. (2007). Relationships between perceived features and similarity of images:A test of tversky’s contrast model.
Journal of the American Society for Information Science& Technology , 58(10):1401–1418.[Rudin et al., 1992] Rudin, L. I., Osher, S., and Fatemi, E. (1992). Nonlinear total variation basednoise removal algorithms.
Physica D , 60:259 – 268.[Santini and Jain, 1997] Santini, S. and Jain, R. (1997). Similarity is a geometer.
Multimedia Toolsand Applications , 5(3):277–306.[Santini and Jain, 1999] Santini, S. and Jain, R. (1999). Similarity measures.
IEEE Transactions oint bi-modal image reconstruction on Pattern Analysis and Machine Intelligence , 21(9):871–883.[Sastry and Carson, 1997] Sastry, S. and Carson, R. E. (1997). Multimodality bayesian algorithmfor image reconstruction in positron emission tomography: a tissue composition model. IEEETransactions on Medical imaging , 16(6):750–761.[Scherzer, 2010] Scherzer, O. (2010).
Handbook of Mathematical Methods in Imaging . SpringerScience & Business Media.[Schramm et al., 2018] Schramm, G., Holler, M., Rezaei, A., Vunckx, K., Knoll, F., Bredies, K.,Boada, F., and Nuyts, J. (2018). Evaluation of parallel level sets and bowsher’s method assegmentation- free anatomical priors for time-of-flight pet reconstruction.
IEEE Transactionson Medical Imaging , 37(2):590–603.[Schweiger and Arridge, 2014] Schweiger, M. and Arridge, S. R. (2014). The Toast++ softwaresuite for forward and inverse modeling in optical tomography.
Journal of Biomedical Optics ,19(4):040801.[Skopal and Bustos, 2011] Skopal, T. and Bustos, B. (2011). On nonmetric similarity searchproblems in complex domains.
ACM Computing Surveys , 43(4).[Somayajula et al., 2011] Somayajula, S., Panagiotou, C., Rangarajan, A., Li, Q. Z., Arridge, S. R.,and Leahy, R. M. (2011). PET image reconstruction using information theoretic anatomicalpriors.
IEEE Transactions on Medical Imaging , 30(3):537–549.[Steiner, 2007] Steiner, G. (2007). Application and data fusion of different sensor modalities intomographic imaging.
Elektrotechnik und Informationstechnik , 124(7):232–239.[Tang and Rahmim, 2009] Tang, J. and Rahmim, A. (2009). Bayesian pet image reconstructionincorporating anato-functional joint entropy.
Physics in Medicine and Biology , 54(23):7063–7075.[Townsend, 2008] Townsend, D. W. (2008). Multimodality imaging of structure and function.
Physics in Medicine and Biology , 53(4):R1–R39.[Tversky, 1977] Tversky, A. (1977). Features of similarity.
Psychological review , 84(4):327.[Tversky and Gati, 1982] Tversky, A. and Gati, I. (1982). Similarity, separability, and the triangleinequality.
Psychological Review , 89(2):123.[Wang et al., 2015] Wang, G., Kalra, M., Murugan, V., Xi, Y., Gjesteby, L., Getzin, M., Yang, Q. S.,Cong, W. X., and Vannier, M. (2015). Vision 20/20: Simultaneous CT-MRI - next chapter ofmultimodality imaging.
Medical Physics , 42(10):5879–5889.[Wang et al., 2012] Wang, G., Zhang, J., Gao, H., Weir, V., Yu, H. Y., Cong, W. X., Xu, X. C., Shen,H. O., Bennett, J., Furth, M., Wang, Y., and Vannier, M. (2012). Towards omni-tomography— grand fusion of multiple modalities for simultaneous interior tomography.
Plos One , 7(6).[Wang et al., 2004] Wang, Z., Bovik, A. C., Sheikh, H. R., and Simoncelli, E. P. (2004). Imagequality assessment: from error visibility to structural similarity.
IEEE transactions on imageprocessing , 13(4):600–612.[Weizman et al., 2016] Weizman, L., Eldar, Y. C., and Ben Bashat, D. (2016). Reference-basedMRI.
Medical Physics , 43(10):5357–5369.[Wikipedia contributors, 2017] Wikipedia contributors (2017). Backtracking line search —Wikipedia, the free encyclopedia. [Online; accessed 20-May-2018].[Yalavarthy et al., 2007] Yalavarthy, P. K., Pogue, B. W., Dehghani, H., Carpenter, C. M., Jiang,S. D., and Paulsen, K. D. (2007). Structural information within regularization matrices improvesnear infrared diffuse optical tomography.
Optics Express , 15(13):8043–8058.[Yang et al., 2015] Yang, Y. D., Wang, K. K. H., Eslami, S., Iordachita, I., Patterson, M. S., andWong, J. W. (2015). Systematic calibration of an integrated x-ray and optical tomographysystem for preclinical radiation research.
Medical Physics , 42(4):1710–1720.[Yuan et al., 2010] Yuan, Z., Zhang, Q., Sobel, E. S., and Jiang, H. (2010). High-resolution X-ray guided three-dimensional diffuse optical tomography of joint tissues in hand osteoarthritis:Morphological and functional assessments.
Medical physics , 37(8):4343–4354.[Zhang et al., 1994] Zhang, Y., Fessler, J. A., Clinthorne, N. H., and Rogers, W. L. (1994).