Multi-modality imaging with structure-promoting regularisers
MMulti-modality imagingwith structure-promoting regularisers
Matthias J. Ehrhardt
Institute for Mathematical Innovation, University of Bath, UK [email protected]
Abstract
Imaging with multiple modalities or multiple channels is becoming increasingly important forour modern society. A key tool for understanding and early diagnosis of cancer and dementiais PET-MR, a combined positron emission tomography and magnetic resonance imaging scannerwhich can simultaneously acquire functional and anatomical data. Similarly in remote sensing, whilehyperspectral sensors may allow to characterise and distinguish materials, digital cameras offer highspatial resolution to delineate objects. In both of these examples, the imaging modalities can beconsidered individually or jointly. In this chapter we discuss mathematical approaches which allowto combine information from several imaging modalities so that multi-modality imaging can be morethan just the sum of its components.
Many tasks in almost all scientific fields can be posed as an inverse problem of the form Ku = f (1)where K is a mathematical model that connects an unknown quantity of interest u to measured data f . The task is to recover u from data f under the model K . In practice this task is difficult becauseof measurement errors in the data f and inaccuracies in the model K . Moreover, in many cases themodel (1) lacks information we have at hand about the unknown quantity u such as its regularity. Inthis chapter we are interested in the situation when have a-priori knowledge about the ”structure” of u from a second measurement v which we want to exploit in the inversion. Throughout this chapterwe will refer to v as the side information . Intuitively, this is the case when u and v describe differentproperties of the same geometry (in medicine: anatomy). We will be more precise in Section 2 where wediscuss mathematical models for structural similarity. The two notions we will discuss in detail that theedges of the two images u and v having similar 1) locations [1–5] and 2) directions [3, 5–16]. Real-worldexamples for these mathematical models are numerous as we will see in the next section. Historically the first application where information from several modalities was combined was positronemission tomography (PET) and magnetic resonance imaging (MRI) in the early 1990s [17]. Sharinginformation between two different imaging modalities is motivated by the fact that all images will behighly influenced by the same underlying anatomy, see Figure 1. Since single-photon emission computedtomography (SPECT) imaging is both mathematically and physically similar to PET imaging, mostof the proposed models can be directly translated and often models are proposed for both modalitiessimultaneously, see e.g. [18–21]. Over the years there always has been research in this direction, seee.g. [18–20, 22–34], which was intensified with the advent of the first simultaneous PET-MR scanner in2011 [35], see e.g. [4, 5, 11–13, 36–43].The same motivation applies to other medical imaging techniques, for example multi-contrast MRI,see e.g. [10, 44–48]. In multi-contrast MRI multiple acquisition sequences are used to acquire data of thesame patient, see Figure 2 for a T - and a T -weighted image with share anatomy. Other special casesare the combination of anatomical MRI (e.g. T -weighted) and magnetic particle imaging [14], functionalMRI (fMRI) and anatomical MRI [49] as well as anatomical ( H) and fluorinated gas ( F) MRI [50]. A1 a r X i v : . [ ee ss . I V ] J u l igure 1: PET-MR and PET-CT.
A low resolution functional PET image (left) is to be reconstructedwith the help of an anatomical MRI (middle) or CT image (right). As is evident from the images, allthree images share many edges due to the same underlying anatomy. Note that the high soft tissuecontrast in MRI makes it favourable over CT for this application. Images curtesy of P. Markiewicz andJ. Schott.Figure 2:
Multi-contrast MRI.
The same MRI scanner can produce different images depending onthe acquisition sequence such as T -weighted (left) and T -weighted images (right). Images courtesy ofN. Burgos.related imaging task is quantitative MRI (such as Magnetic Resonance Fingerprinting [51]) [52–55] whereone aims to reconstruct quantitative maps of tissue parameters (e.g. T , T , proton density, off-resonancefrequency), but regularisers coupling these maps have not been used to date. The idea to couple channelshas also been used for parallel MRI [56].Starting from the 1990s, mathematical models were developed that make use of the expected corre-lations between colour channels of RGB images [58–60], see Figure 3. Research in this field is still veryactive today, see e.g. [2, 8, 61–64].In remote sensing observations are often available from multiple sensors either mounted on a planeor on a satellite. For example a hyperspectral camera with low spatial resolution and a digital camerawith higher spatial resolution may be used simultaneously, see Figure 4. This situation naturally invitesfor the fusion of information, see [15, 65–71] and references therein. In some situations the response ofthe cameras to certain wavelengths is (assumed to be) known such that the data can be fused makinguse of this knowledge. This is commonly referred to as, see e.g. pansharpening [68–70]. It is importantto note that this assumption is sometimes not fulfilled and many of the aforementioned algorithms areflexible enough to fuse data in this more general situation.Dual and spectral computed tomography (CT) is becoming increasingly popular in (bio-) medicalimaging and material sciences due to its ability to distinguish different materials which would not beFigure 3: Color imaging.
The color image (left) is composed of three color channels (right) all of whichshow similar edges due to the same scenery. Images courtesy of M. Ehrhardt.2igure 4:
Hyperspectral imaging + photography.
A nowadays common scenario is that multiplecameras are mounted on a plane or satellite for remote sensing. While one camera carries spectralinformation (right), the other has high spatial resolution (left). Images courtesy of D. Coomes.Figure 5:
Spectral CT.
Standard (white-beam) CT on the left and three channels (28, 34 and 39 keV)of spectral CT on the right of an iodine-stained lizard head reconstructed by CIL [57]. The spectralchannels clearly show a large increase in intensity from 28 to 34 keV, thereby revealing the presence,location and concentration of iodine. Images courtesy of J. Jorgensen and R. Warr.possible using a single energy, see Figure 5. Since the energy channels have a very different signal-to-noise ratio, coupling them within the reconstruction allows to transfer information from high signal tolow signal channels [9, 72–74].In geophysics, the coupling between modalities has been used to model similarity between electricalresistivity and seismic velocity [6, 7], estimating conductivity from multi-frequency data [75], invertinggravity and seismic tomography [75] and controlled-source electromagnetic resistivity inversion [76]. Foran overview and more details on examples in geophysics see in [3, 77] and references therein.Ideas from multi-modality imaging have recently also been used for art restoration. When a canvas ispainted on both sides, an x-ray image shows the superposition of both paintings. The x-ray informationcan then be separated using photos of both sides of the canvas [78].Other examples that were considered in the literature are combining anatomical information andelectrical impedance tomography [16, 79], CT and MRI [80], photo-acoustic and optical coherence to-mography [81], x-ray fluorescence and transmission tomography [82] and various channels in multi-modalelectron tomography [83]. The combination of various imaging modalities into one system may eventuallylead to what is sometimes referred to as omni-tomography [84].Image reconstruction with side information is mathematically similar to multi-modal image registra-tion and thus it is not surprising that both fields share a lot of mathematical models, see e.g. [85–88].
Inverse problems of the form (1) can be solved using variational regularisation, i.e. framed as theoptimisation problem u α ∈ arg min u D ( Au, f ) + α R ( u ) . (2)Here the data fidelity D : Y × Y → R ∞ := R ∪ {∞} measures how close the estimated data Au fits theacquired data f . The regulariser (also referred to as the prior ) R : X → R ∞ defines which propertiesof the image u we favour and which we do not. The trade-off between data fitting and regularisationcan be chosen using the regularisation parameter α >
0. Problems of this form have been extensivelystudied, see for instance [89–93] and references therein.Three popular regularisers for imaging are the squared H -semi norm (H ), the total variation (TV)[94, 95] and the total generalised variation (TGV) [96–98]. It is common to model images as functions3 : Ω ⊂ R d → R . If u is smooth enough, then these regularisers are defined asH ( u ) = (cid:90) Ω |∇ u ( x ) | d x (3)TV( u ) = (cid:90) Ω |∇ u ( x ) | d x (4)TGV( u ) = inf ζ (cid:90) Ω |∇ u ( x ) − ζ ( x ) | + β | Eζ ( x ) | d x . (5)Here ∇ u : Ω → R d , [ ∇ u ] i = ∂ i u denotes the gradient of u , Eζ : Ω → R d × d , [ Eζ ] i,j = ( ∂ i ζ j + ∂ j ζ i ) / ζ : Ω → R d , see [98] for more details, and | · | denotesthe Euclidean/Frobenius norm. For TV and TGV it is of interest to develop other formulations whichare well-defined even when u is not smooth. For simplicity, we do not go into more detail in this directionbut refer the interested reader to the literature, e.g. [95, 96].All three regularisers promote solutions with different smoothness properties. H promotes smoothsolutions with small gradients everywhere, whereas TV promotes solutions which have sparse gradients,i.e. the images are piecewise constant and appear cartoon-like. The latter also leads to the staircaseartefact which can be overcome by TGV which promotes piecewise linear solutions. None of theseregularisers are able to encode additional information on the location or direction of edges. The contributions in this chapter are threefold.
Overview over existing methods
We provide an overview on existing mathematical models forstructural similarity which are related to the shared location or direction of edges. We then discussvarious regularisers which promote similarity in this sense.
Higher order models
Existing methods focus on incorporating additional information into regularis-ers modelling first-order smoothness. We extend existing methodology to second-order smoothness usingthe total generalised variation framework.
Extensive numerical comparison
We highlight the properties of the discussed regularisers and thedependence on various parameters using two inverse problems: tomography and super-resolution.
One can think of the setting (1) with extra information v as a special case when multiple measurements K i u i = f i i = 1 , . . . , m (6)are taken. If m = 2 and one inverse problem is considerably less ill-posed, then this can be solvedfirst to guide the inversion of the other. Some of the described models can be extended to the moregeneral case (e.g. an arbitrary number of modalities) or the joint recovery of both/all unknowns, see e.g.[3–9, 12, 38, 41, 56, 58, 63, 75–77, 82, 83, 99], but it is out of the scope of this chapter to provide anoverview on those. For an overview up to 2015, see [100]. A few recent contributions are summarized in[101].Model (6) may include several special cases i) multiple measurements of the same unknown, i.e. u i = u and ii) measurements correspond to different states of the same unknown, e.g. in dynamicimaging u i = u ( · , t i ). The former case is covered by the standard literature when concatenating themeasurements and the systems models, i.e. ( Ku ) i := K i u and f = ( f , . . . , f m ). The latter has beenwidely studied in the literature, too, see e.g. [102–104] and references therein. Both of these are in generalunrelated to multi-modality imaging. 4 .4.2 Other models for similarity The earliest contributions to structure-promoting regularisers for multi-modality imaging were madein the early 1990s by Leahy and Yan [17] who used a segmentation of an anatomical MRI image toenhance PET reconstruction. This is achieved by carefully handcrafting a regulariser which can encodethis information. In this chapter we will use the same strategy but in a continuous setting which isindependent of the discretisation and will not rely on a segmentation of the side information v . Theseideas were subsequently refined in various directions [18–20, 22–25, 27, 28, 33, 34, 44] of which Bowshersprior [23] remains most popular today.Other models that can combine information of multiple modalities are based on coupled diffusion[1, 61, 99], level-sets [30], information theoretic priors (joint entropy, mutual information) [21, 26, 29, 37],Bregman distances [49, 65, 66, 105, 106], Bregman iterations [64, 107], the structure tensor [108], jointdictionary learning [47, 78, 109], common edge weighting [41] and deep learning [48]. Most of thesemethods are very different to what will be described in this chapter. There are some similarities betweenthe methods of this chapter and methods which are based on the Bregman distance of the total variation[49, 64–66, 105–107] but a detailed treatment is outside the scope of this section. In this section we define mathematical models where we aim to capture the similarities as shown inFigures 1 to 5. We start by explicitly stating two definitions which capture structural similarity whichhave been used implicitly in the literature. The first is based on the location of edges or the edge set[1–5, 41, 56, 64] and the second is based on direction of edges or the shape of an object [3, 5–9, 12].The latter is essentially the same as Definition 5.1.6 in [100] except for the degenerate case when either ∇ u ( x ) = 0 or ∇ v ( x ) = 0. Definition 1 (Structural similarity with edge sets) . Two differentiable images u, v : Ω → R are said tobe structurally similar in the sense of edge sets if E u = E v (7) where E u = { x ∈ Ω | ∇ u ( x ) (cid:54) = 0 } . We also write u e ∼ v to denote that u and v are structurally similar inthe sense of edge sets. Definition 2 (Structural similarity with parallel level sets) . Two differentiable images u, v : Ω → R aresaid to be structurally similar in the sense of parallel level sets if u e ∼ v and for all x ∈ E u there is ∇ u ( x ) (cid:107) ∇ v ( x ) . (8) We also write u d ∼ v to denote that u and v are structurally similar in the sense of parallel level sets. Remark 1.
For smooth images u and v , their gradients are perpendicular to their level sets, i.e. u − ( s ) = { x ∈ Ω | u ( x ) = s } . Thus parallel gradients is equivalent to parallel level sets which explains the naming.The notion that the structure of an image is contained in its level sets dates back to [110]. Remark 2.
By definition, similarity with parallel level sets (Definition 2) is stronger than the definitionthat only involves edge sets (Definition 1). An example of two images u and v which have the sameedge set but do not have parallel level sets is the following. u, v : Ω ⊂ R → R , u ( x ) = x , v ( x ) = x .Clearly they have the same edge set since E u = E v = Ω , but they do not have parallel level sets since ∇ u ( x ) = [1 , but ∇ v ( x ) = [0 , . Remark 3.
Two images u and v have parallel level sets if and only if u e ∼ v and for all x ∈ E u thereexists α ∈ R such that ∇ u ( x ) = α ∇ v ( x ) . (9) Examples of images which have parallel level sets include:1.
Function value transformations . Let f : R → R be smooth and strictly monotonic, i.e. f (cid:48) > or f (cid:48) < . Then v := f ◦ u d ∼ u . This is readily to be seen from the fact that ∇ v ( x ) = f (cid:48) ( u ( x )) ∇ u ( x ) (cid:54) =0 if and only if ∇ u ( x ) (cid:54) = 0 . . Local function value transformations . Let f i : R → R be smooth and strictly monotonic and u = (cid:80) i u i where u i are smooth functions whose gradients have mutually disjoint support. Then v := (cid:80) i f i ◦ u i d ∼ u . Remark 4.
It has been argued in the literature that many multi-modality images z : Ω → R m essentiallydecompose as z i ( x ) = τ i ( x ) ρ ( x ) (10) where ρ ( x ) describes its structure and τ is a material property, see e.g. [63, 111]. Since the materialdoes not change arbitrarily, it is natural to assume that τ i is slowly varying or even piecewise constant.In the latter case, if x is such that ∇ τ i ( x ) = 0 , then we have ∇ z i ( x ) = τ i ( x ) ∇ ρ ( x ) , (11) in particular if τ i , τ j (cid:54) = 0 , then z i d ∼ z j . This property is also related to the material decomposition inspectral CT, see e.g. [112–114]. Measuring the degree of similarity with respect to the previous two definitions of structural similarityis not easy and we will now discuss a couple of ideas from the literature. Here and for the rest ofthis chapter, we will make frequent use of the vector-valued representation of a set of images z : Ω → R , z ( x ) := [ u ( x ) , v ( x )]. We denote by J its Jacobian, i.e. J : Ω → R d × , J i,j = ∂ i z j .With the definition of the Jacobian we see that u e ∼ v if and only if (cid:90) Ω | J ( x ) | d x = (cid:90) Ω |∇ u ( x ) | d x = (cid:90) Ω |∇ v ( x ) | d x (12)where | x | := 1 if x (cid:54) = 0 and 0 else.Similarly, by definition u d ∼ v if and only if u e ∼ v and (a) rank J ( x ) = 1 for all x ∈ E u . (a) isequivalent to (b) a vanishing determinant, i.e. det J (cid:62) ( x ) J ( x ) = 0. Simple calculations, see e.g. [100],show that det J (cid:62) ( x ) J ( x ) = |∇ u ( x ) | |∇ v ( x ) | − (cid:104)∇ u ( x ) , ∇ v ( x ) (cid:105) , (13)where we use the notation (cid:104) x, y (cid:105) = x (cid:62) y for the inner product of two column vectors x and y . In orderto get further equivalent statements we turn to the singular values of the Jacobian which are given by σ / ( x ) = 12 (cid:20) | J ( x ) | ± (cid:113) | J ( x ) | − det J (cid:62) ( x ) J ( x ) (cid:21) (14)with | J ( x ) | = |∇ u ( x ) | + |∇ v ( x ) | , see e.g. [100]. Since σ ( x ) ≥ σ ( x ) ≥ σ ( x ) = 0 or (d) the vector of singular vectors σ ( x ) = [ σ ( x ) , σ ( x )] is 1-sparse. Many of the abstract models from the previous section to measure the degree of similarity with respectto the previous two definitions of structural similarity are computationally challenging as they relateto non-convex constraints. In this section we will define convex structure-promoting regularisers whichmake them computationally tractable.
We first look at isotropic models which only depend on gradient magnitudes rather than directions, thuspromote structural similarity in the sense of edge sets, Definition 1.6 mage η = − η = − η = − Figure 6: Influence of the parameter η on estimation of edge location . The images on the right showthe scalar field w : Ω → [0 ,
1] which locally weights the influence of the regulariser, see (18). Here ”black”denotes 0 and ”white” denotes 1.First, based on (12) if we approximate | J ( x ) | by | J ( x ) | , thenJTV( u ) = (cid:90) Ω | J ( x ) | d x = (cid:90) Ω (cid:112) |∇ u ( x ) | + |∇ v ( x ) | d x (15) ≤ (cid:90) Ω |∇ u ( x ) | + |∇ v ( x ) | d x = TV( u ) + TV( v ) (16)with equality if and only if E u ∩ E v = ∅ . This regulariser is called joint total variation in some commu-nities, see e.g. [3, 5, 11, 56] and vectorial total variation in others, see e.g. [2]. Remark 5.
Note that
JTV has the favourable property that if ∇ v = 0 , then JTV( u ) = TV( u ) , so thatit reduces to a well defined regularisation in u in this degenerate case. Note that this property also holdslocally. Remark 6.
We would also like to note that there is a connection between
JTV and the singular valuesof J . Let σ , σ : Ω → [0 , ∞ ) be the two singular values of J , then we have JTV( u ) = (cid:90) Ω (cid:113) σ ( x ) + σ ( x ) d x . (17)Another strategy to favour edges at similar locations while reducing to a well-defined regulariser inthe degenerate case is to introduce local weighting. Let w : Ω → [0 ,
1] be an edge indicator function for v such that w ( x ) = 1 when ∇ v ( x ) = 0 and a small value whenever |∇ v ( x ) | is large. For example, choose w ( x ) = η (cid:112) η + |∇ v ( x ) | (18)which is illustrated in Figure 6. The figure shows that with a medium η the weight w in (18) shows themain structures of the images so that these can be promoted in the other image. If η is too small, thenalso unwanted structures are captured in w such as a smooth background variation. If η is too large,then the structures start to disappear.For regularisers which are based on the image gradient ∇ u , the weighting w can be used to favouredges at certain locations by replacing ∇ by w ∇ . For instance, for H (3), TV (4) and TGV (5) thisstrategy results in wH ( u ) = (cid:90) Ω | w ( x ) ∇ u ( x ) | d x = (cid:90) Ω w ( x ) |∇ u ( x ) | d x (19)wTV( u ) = (cid:90) Ω | w ( x ) ∇ u ( x ) | d x = (cid:90) Ω w ( x ) |∇ u ( x ) | d x (20)wTGV( u ) = inf ζ (cid:90) Ω | w ( x ) ∇ u ( x ) − ζ ( x ) | + β | Eζ ( x ) | d x (21)7 mage η = − η = − η = − Figure 7: Influence of the parameter η on estimation of edge location and direction . The images onthe right show the vector field ξ : Ω → R d which locally defines the influence of the regulariser, see e.g.(22). Here ”black” denotes that the magnitude of ξ , i.e. | ξ ( x ) | , is 0 and a bright colour denotes that | ξ ( x ) | is 1. The colours show the direction of the vector field ξ modulo its sign.which we will refer to as weighted squared H -semi norm , weighted total variation and weighted totalgeneralised variation . wTV was used in [1, 10, 115]. A variant of wTV has been considered for singlemodality imaging in [116, 117] and extended to a variant of wTGV [118]. Remark 7.
The parameter η in w , see (18) , should be chosen in relation to |∇ v ( x ) | . A common strategyis to normalise the side information first such that sup x ∈ Ω |∇ v ( x ) | = 1 . Then desirable values of η areusually within the range [0 . , . The same idea which resulted in isotropically ”weighted” variants of common regularisers can be usedanisotropically, i.e. by making the local weights vary with direction. Let us denote the anisotropicweighting by D : Ω → R d × d . Similar to the isotropic variant, one would like the weight to become theidentity matrix, i.e. D ( x ) = I , when ∇ v ( x ) = 0. In order to promote parallel level sets it is desirablethat D ( x ) ∇ u ( x ) should be small if ∇ u ( x ) (cid:107) ∇ v ( x ) and D ( x ) ∇ u ( x ) = ∇ u ( x ) if ∇ u ( x ) ⊥ ∇ v ( x ). Forexample D ( x ) = I − γξ ( x ) ξ (cid:62) ( x ) , ξ ( x ) = ∇ v ( x ) (cid:112) η + |∇ v ( x ) | (22)for γ ∈ (0 ,
1] (usually close to 1) and η > ∇ v ( x ) = 0 then ξ = 0 such that D ( x ) = I . Moreover, if ∇ u ( x ) (cid:107) ∇ v ( x ), then there exists an α such that ∇ u ( x ) = α ∇ v ( x )and D ( x ) ∇ u ( x ) = (cid:20) I − γη + |∇ v ( x ) | ∇ v ( x ) ∇ v (cid:62) ( x ) (cid:21) ∇ u ( x ) (23)= (cid:20) − γ |∇ v ( x ) | η + |∇ v ( x ) | (cid:21) ∇ u ( x ) . (24)The scalar weighting factor converges to 1 − γ for |∇ v ( x ) | → ∞ . Finally, if ∇ u ( x ) (cid:107) ∇ v ( x ), then clearly D ( x ) ∇ u ( x ) = ∇ u ( x ).The example of the matrix-field D : Ω → R d × d in (22) is determined by the vector-field ξ : Ω → R d which we visualise in Figure 7. The colours show the direction of the vector-field modulo its sign (since ξ ( x ) ξ (cid:62) ( x ) is invariant to a change of sign) and the brightness indicate its magnitude | ξ ( x ) | . Note thatimages appear as colour versions of Figure 6 which shows the isotropic weighting w .8able 1: Examples of first-order structure-promoting regularisers, see (32).regulariser definition B ( x ) y m φ ( x )H (3) y d | x | wH (19) w ( x ) y d | x | dH (25) D ( x ) y d | x | TV (4) y d | x | wTV (20) w ( x ) y d | x | dTV (26) D ( x ) y d | x | JTV (16) [ y, ξ ( x )] d × | x | TNV (31) [ y, ξ ( x )] d × | x | ∗ Using a matrix-field in common regularisers lead to their ”directional” variantdH ( u ) = (cid:90) Ω | D ( x ) ∇ u ( x ) | d x (25)dTV( u ) = (cid:90) Ω | D ( x ) ∇ u ( x ) | d x (26)dTGV( u ) = inf ζ (cid:90) Ω | D ( x ) ∇ u ( x ) − ζ ( x ) | + β | Eζ ( x ) | d x . (27) Remark 8.
There is a connection between the particular choice of the matrix-field D in (22) and theJacobian J . | D ( x ) ∇ u ( x ) | = |∇ u ( x ) − γη + |∇ v ( x ) | (cid:104)∇ u ( x ) , ∇ v ( x ) (cid:105)∇ v ( x ) | (28)= |∇ u ( x ) | − γη + γ (2 − γ ) |∇ v ( x ) | ( η + |∇ v ( x ) | ) (cid:104)∇ u ( x ) , ∇ v ( x ) (cid:105) . (29) For η = 0 , γ = 1 and |∇ v ( x ) | = 1 , then with (13) we have | D ( x ) ∇ u ( x ) | = |∇ u ( x ) | |∇ v ( x ) | − (cid:104)∇ u ( x ) , ∇ v ( x ) (cid:105) = det J (cid:62) ( x ) J ( x ) . (30) Thus, dH corresponds to penalising the determinant. This regulariser is widely used for joint recon-struction in geophysics under the name cross-gradient function since it is also the cross product of ∇ u ( x ) and ∇ v ( x ) , see e.g. [6, 7, 76, 77]. Similarly the dTV used for instance in medical imaging[10, 11, 13, 14, 16, 50] and remote sensing [15] can be seen as penalising the square root of the determi-nant. Another strategy to promote parallel level sets is via nuclear norm of the Jacobian which is definedas | J ( x ) | ∗ = (cid:80) min( d, i =1 σ i ( x ) where σ i ( x ) denotes the i th singular value of J ( x ). Using the nuclear normpromotes sparse vectors of singular values σ ( x ) = [ σ ( x ) , σ ( x )] and thereby parallel level sets. As aregulariser TNV( u ) = (cid:90) Ω | J ( x ) | ∗ d x (31)this strategy became known as total nuclear variation , see [9, 12, 63, 73].All first-order regularisers of this section can be readily summarised in the following standard form J ( u ) = (cid:90) Ω φ [ B ( x ) ∇ u ( x )] d x (32)where B ( x ) : R d → R m is a an affine transformation and φ : R m → R . For details how B and φ canbe chosen for specific regularisers to fit this framework, see Table 1. It is useful for Jacobian-basedregularisers to use the reweighted Jacobian [ ∇ u ( x ) , ξ ( x )] with ξ ( x ) = η ∇ v ( x ) instead. Note that the solution to variational regularisation (2) with either first- (32) or second-order structuralregularisation (5), (21), (27) can be cast into the general non-smooth composite optimisation formmin x F ( Ax ) + G ( x ) (33)9able 2: Mapping the variational regularisation models into the composite optimisation framework (33).In all cases we choose A x = Ku , F ( y ) = D ( y , b ) and G ( x ) = ı ≥ ( u ).regulariser definition x A x A x F ( y ) F ( y )H (3) u ∇ u - α (cid:107) y (cid:107) -wH (19) u w ∇ u - α (cid:107) y (cid:107) -dH (19) u D ∇ u - α (cid:107) y (cid:107) -TV (4) u ∇ u - α (cid:107) y (cid:107) , -wTV (20) u w ∇ u - α (cid:107) y (cid:107) , -dTV (26) u D ∇ u - α (cid:107) y (cid:107) , -JTV (16) u [ ∇ u,
0] - α (cid:107) y − [0 , ξ ] (cid:107) , -TNV (31) u [ ∇ u,
0] - α (cid:107) y − [0 , ξ ] (cid:107) ∗ , -TGV (5) ( u, ζ ) ∇ u − ζ Eζ α (cid:107) y (cid:107) , αβ (cid:107) y (cid:107) , wTGV (21) ( u, ζ ) w ∇ u − ζ Eζ α (cid:107) y (cid:107) , αβ (cid:107) y (cid:107) , dTGV (27) ( u, ζ ) D ∇ u − ζ Eζ α (cid:107) y (cid:107) , αβ (cid:107) y (cid:107) , with F ( y ) = (cid:80) ni =1 F i ( y i ) and Ax = [ A x, . . . , A n x ], see Table 2. We denote by (cid:107) · (cid:107) , , (cid:107) · (cid:107) and (cid:107) · (cid:107) ∗ , discretisations of z (cid:55)→ (cid:90) Ω | z ( x ) | d x , z (cid:55)→ (cid:90) Ω | z ( x ) | d x and z (cid:55)→ (cid:90) Ω | z ( x ) | ∗ d x . (34) A popular algorithm to solve (33) and therefore (2) is the primal-dual hybrid gradient (PDHG) [119, 120],see Algorithm 1. It consists of two simple steps only involving basic linear algebra and the evaluationof the operator A and its adjoint A ∗ . Moreover, it involves the computation of the proximal operator of τ G and the convex conjugate of σ F ∗ where τ and σ are scalar step sizes. The proximal operator of afunctional H is defined as prox H ( z ) := arg min x (cid:26) (cid:107) x − z (cid:107) + H ( x ) (cid:27) . (35)The proximal operator can be computed in closed-form for (cid:107) · (cid:107) , and (cid:107) · (cid:107) . It also also be computedin closed-form for (cid:107) · (cid:107) ∗ , if either the number channels or the dimension of the domain are strictly lessthan 5, i.e. m, d <
5, see [63] for more details. Note also that the proximal operator of α F ( · − ξ ) canbe readily computed based on the proximal operator of F . More details on proximal operators, convexconjugates and examples can be found for example in [121–124].For some applications (e.g. x-ray tomography) a preconditioned [42, 125] or randomised [42, 126]variant can be useful but we will not consider these here for simplicity. Algorithm 1
Primal-dual hybrid gradient (PDHG) to solve (33). Default values given in brackets.
Input: iterates x (= 0), y (= 0), step size parameter ρ (= 1) Initialize: extrapolation x = x , step sizes σ = ρ/ (cid:107) A (cid:107) , τ = 0 . / ( ρ (cid:107) A (cid:107) ) for k = 1 , . . . do x + = prox τ G ( x − τ A ∗ y ) y + = prox σ F ∗ ( y + σA (2 x + − x )) Since the operator norms (cid:107) A i (cid:107) , i = 1 , . . . n can vary significantly, it is often advisably to ”prewhiten” theproblem by recasting it as min x ˜ F ( ˜ Ax ) + G ( x ) . (36)with ˜ F ( y ) := (cid:80) ni =1 F i ( (cid:107) A i (cid:107) · y i ) and ˜ A i x := A i x/ (cid:107) A i (cid:107) . Then trivially (cid:107) ˜ A i (cid:107) = 1 , i = 1 , . . . , n so that alloperator norms are equal. Note that the proximal operator of σ ˜ F is simple to compute if the proximal10 -ray ground truth side information data super-resolution ground truth side information data Figure 8:
Test cases for numerical experiments . Top: x-ray reconstruction from sparse views andfailed detectors, bottom: super-resolution by a factor of 5 and Gaussian noise.operators of σ F i , i = 1 , . . . , n are simple to compute, since[prox σ ˜ F ( y )] i = λ − i [prox σλ i F i ( λ i y i )] , (37)for any λ i >
0, see for instance [92, Lemma 6.136].
This section describes numerical experiments to compare first- and second-order structure-promotingregularisers.
Software
The numerical computations are carried out in Python using ODL (version 1.0.0.dev0) [127]and ASTRA [128, 129] for computing line integrals in the tomography example. The source code whichreproduces all experiments in this chapter can be found at https://github.com/mehrhardt/Multi-Modality-Imaging-with-Structural-Priors . Data
We consider two test cases with different characteristics, both of which are visualised in Figure8. The first test case, later referred to as x-ray , is parallel beam x-ray reconstruction from only 15 viewswhere additionally some detectors are broken. The latter is modelled by salt-and-pepper noise where 5% of all detectors are corrupted. We aim to recover an image with domain [ − , discretised with 200 pixels. The simulated x-ray camera has 100 detectors and a width of 3 in the same dimensions as theimage domain. Therefore, the challenges are 1) sparse views, 2) small number of detectors and 3) brokendetectors.The second test case, which we refer to as super-resolution , considers the task of super-resolution.Also here we aim to recover an image with domain [ − , discretised with 200 pixels. The forwardoperator is integrating over 5 pixels, thus mapping images of size 200 to images of size 40 . In addition,Gaussian noise of mean zero and standard deviation of 0.01 is added. Algorithmic parameters
We chose the default value ρ = 1 for balancing the step sizes in PDHG andran the algorithm for 3,000 iterations without choosing a specific stopping criterion.11 dge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
759 29.50 .
896 23.80 . wH ( α = ) .
887 32.90 .
962 26.40 . w T V ( α = ) .
632 34.30 .
960 29.40 . w T GV ( α = , β = e − ) Figure 9: Effect of edge weighting on locally weighted models for test case x-ray : increasing edgeparameter η from left to right. All other parameters where tuned to maximize the PSNR and visualimage quality. The multiplicative scaling of an unconstrained optimisation problem is arbitrary, nevertheless we reportthe absolute values here for completeness. For simplicity, all regularisation parameters are shown asmultiples of 1e −
4. The figures at the bottom right of each image are PSNR and SSIM. x-ray
Effect of edge weighting
All structure-promoting regularisers described in Section 3 have in commonthat they rely to some extend on the size of edges in the side information, i.e. |∇ v ( x ) | . For JTV andTNV the actual values of |∇ v ( x ) | matter so that a parameter η is needed to correct for this. For allother regularisers a parameter η is needed to decide which edges to trust and which not. The effect ofthis edge weighting parameter η on all described regularisers is illustrated in Figures 9, 10 and 11. Thelocally weighted regularisers (i.e. wH , wTV and wTGV) and the directional regularisers (i.e. dH , dTVand dTGV) have in common that if η is too small, then small artefacts around the edges appear. Thiseffect is more pronounced in locally weighted regularisers. If η is too large, then the structure-promotingeffect becomes too small. For joint total variation and total nuclear variation similar effects exist withreverse relationship to η . Effect of regularisation
The effect of the regularisation parameter α on the solution is illustratedin Figures 12, 13 and 14. All regularisers show the same behaviour if α is too small or too large. Ifthe regularisation parameter is chosen too small then artefacts from inverting an ill-posed operator areintroduced and if it is chosen too large then all regularisers oversmooth the solution. Note that allstructure-promoting regularisers have an increased robustness in areas of shared structures.12 dge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) edge weightinglow: η =1e − η =1e − η =1e − .
862 30.90 .
914 24.20 . d H ( α = ) .
959 37.30 .
976 27.20 . d T V ( α = ) .
938 39.90 .
978 30.80 . d T GV ( α = , β = e − ) Figure 10: Effect of edge weighting on directional models for test case x-ray : increasing edge parameter η from left to right. All other parameters where tuned to maximize the PSNR and visual image quality( γ = 1). edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =2e − . η =2e + . η =2e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =2e − . η =2e + . η =2e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =2e − . η =2e + . η =2e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =2e − . η =2e + . η =2e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =2e − . η =2e + . η =2e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =2e − . η =2e + . η =2e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =2e − . η =2e + . η =2e + . T NV ( α = ) Figure 11: Effect of edge weighting on joint total variation and total nuclear variation for test case x-ray :increasing edge parameter η from left to right. All other parameters where tuned to maximise the PSNRand visual image quality. 13 egularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) regularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) regularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) regularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) regularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) regularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) regularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) regularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) regularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) regularizationlow medium high α =1 19.60 . α =5 24.70 . α =50 21.90 . H α =5 27.20 . α =20 29.50 . α =100 26.00 . wH ( η = e − ) α =5 27.90 . α =20 30.90 . α =100 28.40 . d H ( η = e − , γ = ) Figure 12: H -semi norm based structure-promoting regularisers for test case x-ray : increasing theregularisation parameter α from left to right. All other parameters where tuned to maximise the PSNRand visual image quality. All regularisers in this figure reduce to the H -semi norm in areas where theside information is flat. 14 egularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) regularizationlow medium high α =50 24.00 . α =150 28.30 . α =300 25.60 . T V α =50 19.70 . α =300 32.90 . α =600 27.30 . w T V ( η = e − ) α =50 24.00 . α =300 37.30 . α =600 30.80 . d T V ( η = e − , γ = ) α =50 24.70 . α =200 34.00 . α =600 29.30 . J T V ( η = ) α =50 27.60 . α =200 36.30 . α =600 29.10 . T NV ( η = ) Figure 13:
Total variation based structure-promoting regularisers for test case x-ray : increasing theregularisation parameter α from left to right. All other parameters where tuned to maximise the PSNRand visual image quality. All regularisers in this figure reduce to the total variation in areas where theside information is flat. 15 egularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) regularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) regularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) regularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) regularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) regularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) regularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) regularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) regularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) regularizationlow medium high α =50 20.70 . α =150 28.50 . α =300 25.20 . T GV α =50 11.90 . α =200 34.30 . α =600 26.90 . w T GV ( η = e − ) α =50 14.20 . α =200 39.90 . α =600 29.80 . d T GV ( η = e − ) Figure 14:
Total generalised variation based structure-promoting regularisers for test case x-ray :increasing the regularisation parameter α from left to right. All other parameters where tuned tomaximise the PSNR and visual image quality ( β = 5e − .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . H .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . .
751 TV 28.30 .
893 TGV 28.50 . .
896 wTV 32.90 .
962 wTGV 34.30 . .
914 dTV 37.30 .
976 dTGV 39.90 . .
958 TNV 36.30 . Figure 15: Comparison of structure-promoting regularisers for test case x-ray . All parameters wheretuned to maximise the PSNR and visual image quality.
Comparison of regularisers
All eleven regularisers are compared in Figure 15. It can be seen that thestructure-promoting regularisers perform much better in terms of PSNR and SSIM as their non-structure-promoting counterparts. Moreover, one can observe an interesting effect that the structure-promotingregularisers also perform visually better in regions where the structure is not shared, e.g. the outer ringof circles. This effect is most dominant for dTGV where the circle at the top left is clearly visible, whileit is difficult to spot for many of the other regularisers. super-resolution
Effect of edge weighting
Figures 16, 17 and 18 show the effect of the edge weighting parameter η .One can make similar observations as in Figures 9, 10 and 11 for the test case x-ray . In addition, onecan observe from the close-ups that if η is too small (or too large for JTV and TNV), then ghostingartefacts may appear. Note that these are present for TNV even for a moderate choice of η . Comparison of regularisers
All regularisers are compared in Figure 19 for the test case super-resolution .It can be noted from all images that introducing structural information allows to resolve some of theinner circles which have been merged for regularisers which are not structure-promoting. Moreover, alltotal generalised variation based regularisers do not perform much better than the total variation based17 dge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
877 33.60 .
885 27.50 . wH ( α = . ) .
932 32.10 .
949 29.30 . w T V ( α = ) .
953 32.60 .
968 29.20 . w T GV ( α = , β = e − ) Figure 16: Effect of edge weighting on locally weighted models for test case super-resolution : increas-ing edge parameter η from left to right. All other parameters where tuned to maximise the PSNR andvisual image quality. 18 dge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) edge weightinglow: η =5e − η =5e − η =5e − .
895 35.00 .
900 27.70 . d H ( α = . ) .
954 35.30 .
957 30.50 . d T V ( α = ) .
975 35.80 .
975 30.40 . d T GV ( α = , β = e − ) Figure 17: Effect of edge weighting on directional models for test case super-resolution : increasingedge parameter η from left to right. All other parameters where tuned to maximise the PSNR and visualimage quality ( γ = 0 . edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =1e − . η =1e + . η =1e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =1e − . η =1e + . η =1e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =1e − . η =1e + . η =1e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =1e − . η =1e + . η =1e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =1e − . η =1e + . η =1e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =1e − . η =1e + . η =1e + . T NV ( α = ) edge weightinglow medium high η =2e − . η =2e + . η =2e + . J T V ( α = ) η =1e − . η =1e + . η =1e + . T NV ( α = ) Figure 18: Effect of edge weighting on joint total variation and total nuclear variation for test case super-resolution : increasing edge parameter η from left to right. All other parameters where tunedto maximise the PSNR and visual image quality. 19 .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . H .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . .
804 TV 28.00 .
915 TGV 28.10 . .
885 wTV 32.10 .
949 wTGV 32.60 . .
900 dTV 35.30 .
957 dTGV 35.80 . .
969 TNV 35.30 . Figure 19: Comparison of structure-promoting regularisers for test case super-resolution . All param-eters where tuned to maximise the PSNR and visual image quality.regularisers. The directional regularisers as well as JTV and TNV perform best in terms of PSNR forthis example.
The median computing times for the numerical experiments are reported in Table 3. The computingtime of PDHG is mainly influenced by the dimensions of the models, the proximal operator and theforward model. As can be seen from the table, H and TV are roughly the same fast. TGV whichuses a second primal variable in the space of the image gradient is significantly slower with about twicethe computational cost. In all three cases, introducing isotropic weights (i.e. wH , wTV and wTGV)increases the cost by about 6 seconds, and anisotropic weights (i.e. dH , dTV and dTGV) by about 12seconds. JTV is more costly than dTV but not as costly as TGV. TNV is by far the most costly of allalgorithms due to the need to compute singular value decompositions of 2 × x-ray super-resolution x-ray super-resolution H , TV or TGV based regulariser is desirable depends on each individual application. For eachof them, there is a clear trend that one achieves better image quality by introducing more information, i.e.first isotropic information and then anisotropic information, each of which increases their computationalcost. However, the increase in computational cost is so small that for most applications the directionalvariant is likely to be favoured. This chapter introduced fundamental mathematical concepts on the structure of images and how struc-tural similarity between images can be measured. The fundamental building blocks are the similaritybased on edge sets and parallel level sets. These notions lead to several classes of structure-promotingregularisers all of which are convex and thereby lead to tractable optimisation problems when used invariational regularisation for linear inverse problems with convex data fits. While some of the regularis-ers are smooth and others are non-smooth, the resulting optimisation problem for all of them can beefficiently computed by PDHG. The effectiveness of these regularisers for the promotion of structure hasbeen observed in many applications and was also illustrated in this chapter on two simulation studies.
The mathematical framework for structure-promoting regularisers is by now well established and fairlymature. Open problems reside in practical problems in the translation of these techniques to applicationswhich will also motivate further mathematical research.
Misregistration
The biggest open problem is misregistration. All of the described regularisers assumethat both images are perfectly aligned. Even in scanners which have two imaging modalities in thesame system such as PET-MR, this assumptions is never perfectly fulfilled. This issue has not beenaddressed much in the literature. In [40], the authors proposed an alternating approach between imagereconstruction and image registration with some success. In [15], the problem was formulated as a blinddeconvolution problem so that translations can be compensated with a shifted kernel. The resultingoptimisation algorithm is related to alternating gradient descent, thus alternates between incrementalupdates of each variable. This approach seems fairly efficient and robust but is limited to compensatingtranslations and fails to find large deformations. The latter was approached heuristically in [71] but agenerally applicable strategy is still to be found.
Extensions beyond two modalities
It is natural to consider the case that more than one imageis available as side information. For instance, in some remote sensing applications a colour photographwith high spatial resolution is available. Similarly in PET-MR, images of more than one MR sequence21ight be available. This setting has also been considered in [39] for a purely discrete model. Some ofthe regularisers to promote structural similarity in this chapter naturally extend to multiple images asside information, but this has not yet been properly investigated.
Applications
As we illustrated in Section 1, there are many applications where structure-promotingregularisers were already used or are on the horizon. The list of potential target applications growssteadily with more and more multi-modality scanners being introduced. Next to the misregistrationmentioned before, the biggest hurdle in real applications is the interpretation of images that were createdby fusing information from several modalities. A common question is: ”Which edges can I trust?” sinceoften the reconstruction from multi-modality data would be performed on a finer resolution than forthe single-modality case. For example, for PET-MR the reconstruction of PET data with an alreadyreconstructed MR image as side information can be performed on the native MRI resolution. The answermight be that such an image should not be interpreted as a PET image, but in fact as a synergistic PET-MR image.
Joint reconstruction
Throughout this chapter the focus was on improving the reconstruction of oneimage with the aid of another modality used as side information. Since the other image is rarely acquireddirectly, it is natural to aim to reconstruct both images simultaneously rather than sequentially. Whileconceptually appealing this strategy leads to many more complications than the approach discussedin this chapter which is sometimes referred to as one-sided reconstruction. While the mathematicalframework for one-sided reconstruction is quite mature, the framework for joint reconstruction is despitea lot of research effort in the last 10 years still in its infancy. Fundamental problems like computationallytractable and efficient coupling of modalities is still unsolved. The appealing strategy of making use ofthe solid mathematical foundations of one-sided reconstruction for joint reconstruction in a mathematicalsound and computationally tractable way is still not possible to date.
Acknowledgements
The author acknowledges support from the EPSRC grant EP/S026045/1 and the Faraday InstitutionEP/T007745/1. Moreover, the author is grateful to all his collaborators which indirectly contributed tothis chapter over the last couple of years.
References [1] S. R. Arridge, V. Kolehmainen, and M. J. Schweiger, “Reconstruction and Regularisation in OpticalTomography,” in
Mathematical Methods in Biomedical Imaging and Intensity-Modulated RadiationTherapy (IMRT) (A. Censor, Y. and Jiang, M. and Louis, ed.), Scuola Normale Superiore, 2008.[2] X. Bresson and T. F. Chan, “Fast Dual Minimization of the Vectorial Total Variation Norm andApplications to Color Image Processing,”
Inverse Problems and Imaging , vol. 2, pp. 455–484, 112008.[3] E. Haber and M. Holtzman-Gazit, “Model Fusion and Joint Inversion,”
Surveys in Geophysics ,pp. 675–695, 8 2013.[4] F. Knoll, T. Koesters, R. Otazo, F. Boada, and D. K. Sodickson, “Simultaneous MR-PET Recon-struction using Multi Sensor Compressed Sensing and Joint Sparsity,” in
International Society forMagnetic Resonance in Medicine , vol. 22, 2014.[5] M. J. Ehrhardt, K. Thielemans, L. Pizarro, D. Atkinson, S. Ourselin, B. F. Hutton, and S. R.Arridge, “Joint Reconstruction of PET-MRI by exploiting Structural Similarity,”
Inverse Problems ,vol. 31, p. 015001, 1 2015.[6] L. A. Gallardo and M. A. Meju, “Characterization of Heterogeneous Near-Surface Materials byJoint 2D Inversion of DC Resistivity and Seismic Data,”
Geophysical Research Letters , vol. 30,no. 13, p. 1658, 2003.[7] L. A. Gallardo and M. A. Meju, “Joint Two-Dimensional DC Resistivity and Seismic Travel TimeInversion with Cross-Gradients Constraints,”
Journal of Geophysical Research , vol. 109, no. B3,pp. 1–11, 2004. 228] M. J. Ehrhardt and S. R. Arridge, “Vector-Valued Image Processing by Parallel Level Sets,”
IEEETransactions on Image Processing , vol. 23, no. 1, pp. 9–18, 2014.[9] D. Rigie and P. La Riviere, “Joint Reconstruction of Multi-Channel, Spectral CT Data via Con-strained Total Nuclear Variation Minimization,”
Physics in Medicine and Biology , vol. 60, pp. 1741–1762, 2015.[10] M. J. Ehrhardt and M. M. Betcke, “Multi-Contrast MRI Reconstruction with Structure-GuidedTotal Variation,”
SIAM Journal on Imaging Sciences , vol. 9, no. 3, pp. 1084–1106, 2016.[11] M. J. Ehrhardt, P. J. Markiewicz, M. Liljeroth, A. Barnes, V. Kolehmainen, J. Duncan, L. Pizarro,D. Atkinson, B. F. Hutton, S. Ourselin, K. Thielemans, and S. R. Arridge, “PET Reconstructionwith an Anatomical MRI Prior using Parallel Level Sets,”
IEEE Transactions on Medical Imaging ,vol. 35, pp. 2189–2199, 8 2016.[12] F. Knoll, M. Holler, T. Koesters, R. Otazo, K. Bredies, and D. K. Sodickson, “Joint MR-PETReconstruction using a Multi-Channel Image Regularizer,”
IEEE Transactions on Medical Imaging ,vol. 36, no. 1, 2016.[13] G. Schramm, M. Holler, A. Rezaei, K. Vunckx, F. Knoll, K. Bredies, F. Boada, and J. Nuyts,“Evaluation of Parallel Level Sets and Bowsher’s Method as Segmentation-Free Anatomical Priorsfor Time-of-Flight PET Reconstruction,”
IEEE Transactions on Medical Imaging , vol. 62, no. 2,pp. 590 – 603, 2017.[14] C. Bathke, T. Kluth, and P. Maass, “Improved Image Reconstruction in Magnetic Particle Imagingusing Structural a priori Information,”
International Journal of Magnetic Particle Imaging , vol. 3,no. 1, 2017.[15] L. Bungert, D. A. Coomes, M. J. Ehrhardt, J. Rasch, R. Reisenhofer, and C.-B. Sch¨onlieb, “BlindImage Fusion for Hyperspectral Imaging with the Directional Total Variation,”
Inverse Problems ,vol. 34, no. 4, p. 044003, 2018.[16] V. Kolehmainen, M. J. Ehrhardt, and S. R. Arridge, “Incorporating Structural Prior Informationand Sparsity into EIT using Parallel Level Sets,”
Inverse Problems and Imaging , vol. 13, no. 2,pp. 285–307, 2019.[17] R. M. Leahy and X. Yan, “Incorporation of Anatomical MR Data for Improved Functional Imagingwith PET,” in
Information Processing in Medical Imaging , pp. 105–120, Springer, 1991.[18] J. E. Bowsher, V. E. Johnson, T. G. Turkington, R. J. Jaszczak, C. E. Floyd, and R. E. Coleman,“Bayesian Reconstruction and Use of Anatomical A Priori Information for Emission Tomography,”
IEEE Transactions on Medical Imaging , vol. 15, no. 5, pp. 673—-686, 1996.[19] A. Rangarajan, I.-T. Hsiao, and G. Gindi, “A Bayesian Joint Mixture Framework for the Inte-gration of Anatomical Information in Functional Image Reconstruction,”
Journal of MathematicalImaging and Vision , vol. 12, no. 3, pp. 199–217, 2000.[20] C. Chan, R. Fulton, D. D. Feng, W. Cai, and S. Meikle, “An Anatomically based Regionally Adap-tive Prior for MAP Reconstruction in Emission Tomography,” in
IEEE Nuclear Science Symposiumand Medical Imaging Conference , pp. 4137–4141, 2007.[21] J. Nuyts, “The Use of Mutual Information and Joint Entropy for Anatomical Priors in EmissionTomography,” in
IEEE Nuclear Science Symposium and Medical Imaging Conference , pp. 4149–4154, IEEE, 2007.[22] C. Comtat, P. E. Kinahan, J. A. Fessler, T. Beyer, D. W. Townsend, M. Defrise, and C. J. Michel,“Clinically Feasible Reconstruction of 3D Whole-Body PET/CT Data using Blurred AnatomicalLabels,”
Physics in Medicine and Biology , vol. 47, pp. 1–20, 1 2002.[23] J. E. Bowsher, H. Yuan, L. W. Hedlund, T. G. Turkington, G. Akabani, A. Badea, W. C. Kurylo,C. T. Wheeler, G. P. Cofer, M. W. Dewhirst, and G. A. Johnson, “Utilizing MRI Information toEstimate F18-FDG Distributions in Rat Flank Tumors,” in
IEEE Nuclear Science Symposium andMedical Imaging Conference , pp. 2488–2492, 2004.2324] K. Baete, J. Nuyts, W. Van Paesschen, P. Suetens, and P. Dupont, “Anatomical-based FDG-PETReconstruction for the Detection of Hypo-Metabolic Regions in Epilepsy,”
IEEE Transactions onMedical Imaging , vol. 23, no. 4, pp. 510–9, 2004.[25] C. Chan, R. Fulton, D. D. Feng, and S. Meikle, “Regularized Image Reconstruction with anAnatomically Adaptive Prior for Positron Emission Tomography,”
Physics in Medicine and Biol-ogy , vol. 54, no. 24, pp. 7379–400, 2009.[26] J. Tang and A. Rahmim, “Bayesian PET Image Reconstruction Incorporating Anato-FunctionalJoint Entropy,”
Physics in Medicine and Biology , vol. 54, no. 23, pp. 7063–75, 2009.[27] A. Bousse, S. Pedemonte, D. Kazantsev, S. Ourselin, S. R. Arridge, and B. F. Hutton, “WeightedMRI-based Bowsher Priors for SPECT Brain Image Reconstruction,” in
IEEE Nuclear ScienceSymposium and Medical Imaging Conference , pp. 3519–3522, 2010.[28] S. Pedemonte, A. Bousse, B. F. Hutton, S. R. Arridge, and S. Ourselin, “Probabilistic GraphicalModel of SPECT/MRI,” in
Machine Learning in Medical Imaging , pp. 167–174, 2011.[29] S. Somayajula, C. Panagiotou, A. Rangarajan, Q. Li, S. R. Arridge, and R. M. Leahy, “PET ImageReconstruction using Information Theoretic Anatomical Priors,”
IEEE Transactions on MedicalImaging , vol. 30, no. 3, pp. 537–549, 2011.[30] J. Cheng-Liao and J. Qi, “PET Image Reconstruction with Anatomical Edge Guided Level SetPrior,”
Physics in Medicine and Biology , vol. 56, pp. 6899–6918, 2011.[31] K. Vunckx, A. Atre, K. Baete, A. Reilhac, C. M. Deroose, K. Van Laere, and J. Nuyts, “Evaluationof Three MRI-based Anatomical Priors for Quantitative PET Brain Imaging,”
IEEE Transactionson Medical Imaging , vol. 31, no. 3, pp. 599–612, 2012.[32] D. Kazantsev, S. R. Arridge, S. Pedemonte, A. Bousse, K. Erlandsson, B. F. Hutton, andS. Ourselin, “An Anatomically Driven Anisotropic Diffusion Filtering Method for 3D SPECTReconstruction,”
Physics in Medicine and Biology , vol. 57, no. 12, pp. 3793–3810, 2012.[33] A. Bousse, S. Pedemonte, B. A. Thomas, K. Erlandsson, S. Ourselin, S. R. Arridge, and B. F.Hutton, “Markov Random Field and Gaussian Mixture for Segmented MRI-based Partial VolumeCorrection in PET,”
Physics in Medicine and Biology , vol. 57, pp. 6681–705, 10 2012.[34] B. Bai, Q. Li, and R. M. Leahy, “Magnetic Resonance-guided Positron Emission TomographyImage Reconstruction,”
Seminars in Nuclear Medicine , vol. 43, pp. 30–44, 2013.[35] G. Delso, S. Furst, B. Jakoby, R. Ladebeck, C. Ganter, S. G. Nekolla, M. Schwaiger, S. I. Ziegler,S. F¨urst, B. Jakoby, R. Ladebeck, C. Ganter, S. G. Nekolla, M. Schwaiger, and S. I. Ziegler, “Per-formance measurements of the Siemens mMR integrated whole-body PET/MR scanner,”
Journalof Nuclear Medicine , vol. 52, pp. 1914–22, dec 2011.[36] M. J. Ehrhardt, K. Thielemans, L. Pizarro, P. J. Markiewicz, D. Atkinson, S. Ourselin, B. F.Hutton, and S. R. Arridge, “Joint Reconstruction of PET-MRI by Parallel Level Sets,” in
IEEENuclear Science Symposium and Medical Imaging Conference , 2014.[37] J. Tang and A. Rahmim, “Anatomy Assisted PET Image Reconstruction Incorporating Multi-Resolution Joint Entropy,”
Physics in Medicine and Biology , vol. 60, no. 1, pp. 31–48, 2015.[38] A. Mehranian, M. Belzunce, C. Prieto, A. Hammers, and A. J. Reader, “Synergistic PET andSENSE MR Image Reconstruction using Joint Sparsity Regularization,”
IEEE Transactions onMedical Imaging , vol. 37, no. 1, pp. 20 – 34, 2018.[39] A. Mehranian, M. A. Belzunce, F. Niccolini, M. Politis, C. Prieto, F. Turkheimer, A. Hammers,and A. J. Reader, “PET Image Reconstruction using Multi-Parametric Anato-Functional,”
PhysMed Biol , 2017.[40] Y.-j. Tsai, S. Member, A. Bousse, S. Ahn, W. Charles, S. Arridge, B. F. Hutton, S. Member,and K. Thielemans, “Algorithms for Solving Misalignment Issues in Penalized PET / CT Recon-struction Using Anatomical Priors,” in
IEEE Nuclear Science Symposium and Medical ImagingConference Proceedings (NSS/MIC) , IEEE, 2018.2441] Y. Zhang and X. Zhang, “PET-MRI joint reconstruction with common edge weighted total varia-tion regularization,”
Inverse Problems , vol. 34, no. 6, p. 065006, 2018.[42] M. J. Ehrhardt, P. J. Markiewicz, and C.-B. Sch¨onlieb, “Faster PET reconstruction with non-smooth priors by randomization and preconditioning,”
Physics in Medicine & Biology , vol. 64,no. 22, p. 225019, 2019.[43] D. Deidda, N. A. Karakatsanis, P. M. Robson, Y. J. Tsai, N. Efthimiou, K. Thielemans, Z. A.Fayad, R. G. Aykroyd, and C. Tsoumpas, “Hybrid PET-MR list-mode kernelized expectationmaximization reconstruction,”
Inverse Problems , vol. 35, no. 4, 2019.[44] B. Bilgic, V. K. Goyal, and E. Adalsteinsson, “Multi-Contrast Reconstruction with Bayesian Com-pressed Sensing,”
Magnetic Resonance in Medicine , vol. 66, pp. 1601–15, 12 2011.[45] J. Huang, C. Chen, and L. Axel, “Fast Multi-contrast MRI Reconstruction,”
Magnetic ResonanceImaging , vol. 32, no. 10, pp. 1344–52, 2014.[46] D. K. Sodickson, L. Feng, F. Knoll, M. Cloos, N. Ben-Eliezer, L. Axel, H. Chandarana, K. T.Block, and R. Otazo, “The Rapid Imaging Renaissance: Sparser Samples, Denser Dimensions, andGlimmerings of a Grand Unified Tomography,” in
Proceedings of SPIE , vol. 9417, pp. 94170G1–14,2015.[47] P. Song, L. Weizman, J. F. Mota, Y. C. Eldar, and M. R. Rodrigues, “Coupled Dictionary Learningfor Multi-Contrast MRI Reconstruction,” in
International Conference on Image Processing , no. 2,pp. 2880–2884, 2018.[48] L. Xiang, Y. Chen, W. Chang, Y. Zhan, W. Lin, Q. Wang, and D. Shen, “Deep-Learning-BasedMulti-Modal Fusion for Fast MR Reconstruction,”
IEEE Transactions on Biomedical Engineering ,vol. 66, no. 7, pp. 2105–2114, 2019.[49] J. Rasch, V. Kolehmainen, R. Nivajarvi, M. Kettunen, O. Gr¨ohn, M. Burger, and E. M. Brinkmann,“Dynamic MRI reconstruction from undersampled data with an anatomical prescan,”
Inverse Prob-lems , vol. 34, no. 7, 2018.[50] A. J. Obert, M. Gutberlet, A. L. Kern, T. F. Kaireit, R. Grimm, F. Wacker, and J. Vogel-Claussen,“1H-guided reconstruction of 19F gas MRI in COPD patients,”
Magnetic Resonance in Medicine ,no. January, pp. 1–11, 2020.[51] D. Ma, V. Gulani, N. Seiberlich, K. Liu, J. L. Sunshine, J. L. Duerk, and M. A. Griswold, “MagneticResonance Fingerprinting,”
Nature , vol. 495, pp. 187–92, mar 2013.[52] M. Davies, G. Puy, P. Vandergheynst, and Y. Wiaux, “A Compressed Sensing Framework forMagnetic Resonance Fingerprinting,”
SIAM Journal on Imaging Sciences , vol. 7, no. 4, pp. 2623–2656, 2013.[53] S. Tang, C. Fernandez-Granda, S. Lannuzel, B. Bernstein, R. Lattanzi, M. Cloos, F. Knoll, andJ. Asslander, “Multicompartment magnetic resonance fingerprinting,”
Inverse Problems , vol. 34,no. 9, 2018.[54] G. Dong, M. Hinterm¨uller, and K. Papafitsoros, “Quantitative magnetic resonance imaging: Fromfingerprinting to integrated physics-based models,”
SIAM Journal on Imaging Sciences , vol. 12,no. 2, pp. 927–971, 2019.[55] M. Golbabaee, Z. Chen, Y. Wiaux, and M. Davies, “CoverBLIP: accelerated and scalable iterativematched-filtering for magnetic resonance fingerprint reconstruction,”
Inverse Problems , vol. 36,no. 1, p. 015003, 2020.[56] C. Chen, Y. Li, and J. Huang, “Calibrationless Parallel MRI with Joint Total Variation Regular-ization,” in
Medical Image Computing and Computer-Assisted Intervention , pp. 106–114, 2013.[57] E. Ametova, G. Fardell, J. S. Jørgensen, W. R. B. Lionheart, E. Papoutsellis, E. Pasca, D. Sykes,M. Turner, R. Warr, and P. J. Withers, “Core Imaging Library (CIL),” 2019.[58] G. Sapiro and D. L. Ringach, “Anisotropic Diffusion of Multivalued Images with Applications toColor Filtering,”
IEEE Transactions on Image Processing , vol. 5, pp. 1582–6, 1 1996.2559] P. Blomgren and T. F. Chan, “Color TV: Total Variation Methods for Restoration of Vector-ValuedImages.,”
IEEE Transactions on Image Processing , vol. 7, no. 3, pp. 304–309, 1998.[60] N. Sochen, R. Kimmel, and R. Malladi, “A General Framework for Low Level Vision,”
IEEETransactions on Image Processing , vol. 7, no. 3, pp. 310–318, 1998.[61] D. Tschumperl´e and R. Deriche, “Vector-Valued Image Regularization with PDEs: A CommonFramework for Different Applications.,”
IEEE Transactions on Pattern Analysis and MachineIntelligence , vol. 27, pp. 506–17, 4 2005.[62] B. Goldluecke, E. Strekalovskiy, and D. Cremers, “The Natural Vectorial Total Variation WhichArises from Geometric Measure Theory,”
SIAM Journal on Imaging Sciences , vol. 5, pp. 537–563,1 2012.[63] K. M. Holt, “Total Nuclear Variation and Jacobian Extensions of Total Variation for Vector Fields,”
IEEE Transactions on Image Processing , vol. 23, no. 9, pp. 3975–3989, 2014.[64] M. M¨oller, E.-M. Brinkmann, M. Burger, and T. Seybold, “Color Bregman TV,”
SIAM Journalon Imaging Sciences , vol. 7, pp. 2771–2806, 12 2014.[65] C. Ballester, V. Caselles, L. Igual, J. Verdera, and B. Roug´e, “A Variational Model for P+XSImage Fusion,”
International Journal of Computer Vision , vol. 69, pp. 43–58, 4 2006.[66] M. M¨oller, T. Wittman, A. L. Bertozzi, and M. Burger, “A Variational Approach for SharpeningHigh Dimensional Images,”
SIAM Journal on Imaging Sciences , vol. 5, pp. 150–178, 1 2012.[67] F. Fang, F. Li, C. Shen, and G. Zhang, “A Variational Approach for Pan-Sharpening,”
IEEETransactions on Image Processing , vol. 22, no. 7, pp. 2822–2834, 2013.[68] L. Loncan, L. B. De Almeida, J. M. Bioucas-Dias, X. Briottet, J. Chanussot, N. Dobigeon, S. Fabre,W. Liao, G. A. Licciardi, M. Simoes, J. Y. Tourneret, M. A. Veganzones, G. Vivone, Q. Wei, andN. Yokoya, “Hyperspectral Pansharpening: A Review,”
IEEE Geoscience and Remote SensingMagazine , vol. 3, no. 3, pp. 27–46, 2015.[69] N. Yokoya, C. Grohnfeldt, and J. Chanussot, “Hyperspectral and multispectral data fusion: Acomparative review of the recent literature,”
IEEE Geoscience and Remote Sensing Magazine ,vol. 5, no. 2, pp. 29–56, 2017.[70] J. Duran, A. Buades, B. Coll, C. Sbert, and G. Blanchet, “A survey of pansharpening methodswith a new band-decoupled variational model,”
ISPRS Journal of Photogrammetry and RemoteSensing , vol. 125, pp. 78–105, 2017.[71] L. Bungert, M. J. Ehrhardt, and R. Reisenhofer, “Robust Blind Image Fusion for MisalignedHyperspectral Imaging Data,” in
Proceedings in Applied Mathematics & Mechanics , vol. 18,p. e201800033, 2018.[72] R. Foygel Barber, E. Y. Sidky, T. Gilat Schmidt, and X. Pan, “An algorithm for constrained one-step inversion of spectral CT data,”
Physics in Medicine and Biology , vol. 61, no. 10, pp. 3784–3818,2016.[73] D. S. Rigie, A. A. Sanchez, and P. J. La Rivi´ere, “Assessment of vectorial total variation penaltieson realistic dual-energy CT data,”
Physics in Medicine and Biology , vol. 62, no. 8, pp. 3284–3298,2017.[74] D. Kazantsev, J. S. Jørgensen, M. S. Andersen, W. R. Lionheart, P. D. Lee, and P. J. With-ers, “Joint image reconstruction method with correlative multi-channel prior for x-ray spectralcomputed tomography,”
Inverse Problems , vol. 34, no. 6, 2018.[75] E. Haber and D. W. Oldenburg, “Joint Inversion: A Structural Approach,”
Inverse Problems ,vol. 13, pp. 63–77, 1997.[76] M. A. Meju, R. L. Mackie, F. Miorelli, A. S. Saleh, and R. V. Miller, “Structurally-tailored 3Danisotropic CSEM resistivity inversion with cross-gradients criterion and simultaneous model cali-bration,”
Geophysics , vol. 84, no. 6, pp. 1–62, 2019.2677] L. A. Gallardo and M. A. Meju, “Structure-Coupled Multiphysics Imaging in Geophysical Sci-ences,”
Reviews of Geophysics , vol. 49, pp. 1–19, 2011.[78] N. Deligiannis, J. F. Mota, B. Cornelis, M. R. Rodrigues, and I. Daubechies, “Multi-modal dictio-nary learning for image separation with application in art investigation,”
IEEE Transactions onImage Processing , vol. 26, no. 2, pp. 751–764, 2017.[79] J. P. Kaipio, V. Kolehmainen, M. Vauhkonen, and E. Somersalo, “Inverse Problems with StructuralPrior Information,”
Inverse Problems , vol. 15, no. 3, pp. 713–729, 1999.[80] Y. Xi, J. Zhao, J. Bennett, M. Stacy, A. Sinusas, and G. Wang, “Simultaneous CT-MRI recon-struction for constrained imaging geometries using structural coupling and compressive sensing,”
IEEE Transactions on Biomedical Engineering , 2015.[81] P. Elbau, L. Mindrinos, and O. Scherzer, “Quantitative reconstructions in multi-modal photoa-coustic and optical coherence tomography imaging,”
Inverse Problems , vol. 34, no. 1, 2018.[82] Z. W. Di, S. Leyffer, and S. M. Wild, “Optimization-Based Approach for Joint X-Ray Fluorescenceand Transmission Tomographic Inversion,”
SIAM Journal on Imaging Sciences , vol. 9, no. 1, pp. 1–23, 2016.[83] R. Huber, G. Haberfehlner, M. Holler, and K. Bredies, “Total Generalized Variation regularizationfor multi-modal electron tomography,”
Nanoscale , pp. 1–38, 2019.[84] G. Wang, J. Zhang, H. Gao, V. Weir, H. Yu, W. Cong, X. Xu, H. Shen, J. Bennett, M. Furth,Y. Wang, and M. Vannier, “Towards Omni-Tomography - Grand Fusion of Multiple Modalities forSimultaneous Interior Tomography,”
PloS one , vol. 7, p. e39700, jan 2012.[85] W. M. Wells III, P. Viola, H. Atsumi, S. Nakajima, and R. Kikinis, “Multi-Modal Volume Regis-tration by Maximization of Mutual Information,”
Medical Image Analysis , vol. 1, no. 1, pp. 35–51,1996.[86] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, and P. Suetens, “Multimodality ImageRegistration by Maximization of Mutual Information,”
IEEE Transactions on Medical Imaging ,vol. 16, pp. 187–98, apr 1997.[87] J. P. W. Pluim, J. B. A. Maintz, and M. A. Viergever, “Image Registration by Maximization ofCombined Mutual Information and Gradient Information,”
IEEE Transactions on Medical Imag-ing , vol. 19, no. 8, pp. 809–14, 2000.[88] E. Haber and J. Modersitzki, “Intensity Gradient based Registration and Fusion of Multi-ModalImages,” in
Medical Image Computing and Computer-Assisted Intervention , vol. 46, pp. 726–733,Berlin Heidelberg: Springer-Verlag, 2006.[89] H. W. Engl, M. Hanke, and A. Neubauer,
Regularization of Inverse Problems . Mathematics andIts Applications, Springer, 1996.[90] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen,
Variational Methods inImaging , vol. 167. 2008.[91] K. Ito and B. Jin,
Inverse Problems - Tikhonov Theory and Algorithms . World Scientific Publishing,2014.[92] K. Bredies and D. A. Lorenz,
Mathematical Image Processing . Birkh¨auser Basel, 1 ed., 2018.[93] M. Benning and M. Burger, “Modern regularization methods for inverse problems,”
Acta Numerica ,vol. 27, pp. 1–111, 2018.[94] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear Total Variation based Noise Removal Algorithms,”
Physica D: Nonlinear Phenomena , vol. 60, no. 1, pp. 259–268, 1992.[95] M. Burger and S. Osher, “A Guide to the TV Zoo,” in
Level Set and PDE Based ReconstructionMethods in Imaging , vol. 2090 of
Lecture Notes in Mathematics , pp. 1–70, Springer, 2013.2796] K. Bredies, K. Kunisch, and T. Pock, “Total Generalized Variation,”
SIAM Journal on ImagingSciences , vol. 3, no. 3, pp. 492–526, 2010.[97] K. Bredies and M. Holler, “Regularization of Linear Inverse Problems with Total GeneralizedVariation,”
Journal of Inverse and Ill-Posed Problems , vol. 22, no. 6, pp. 871–913, 2014.[98] K. Bredies and M. Holler, “A TGV-Based Framework for Variational Image Decompression, Zoom-ing, and Reconstruction. Part II: Numerics,”
SIAM Journal on Imaging Sciences , vol. 8, no. 4,pp. 2851–2886, 2015.[99] S. R. Arridge and A. Simmons, “Multi-Spectral Probabilistic Diffusion Using Bayesian Classifi-cation,” in
Scale-Space Theories in Computer Vision (B. M. ter Haar Romeny, L. Florack, J. J.Koenderink, and M. A. Viergever, eds.), pp. 224–235, Berlin: Springer, 1997.[100] M. J. Ehrhardt,
Joint Reconstruction for Multi-Modality Imaging with Common Structure . Phdthesis, University College London, 2015.[101] S. R. Arridge, M. Burger, and M. J. Ehrhardt, “Preface to special issue on joint reconstructionand multi-modality / multi-spectral imaging,”
Inverse Problems , vol. 36, p. 020302, 2020.[102] U. Schmitt and A. K. Louis, “Efficient algorithms for the regularization of dynamic inverse prob-lems: I. Theory,”
Inverse Problems , vol. 18, no. 3, pp. 645–658, 2002.[103] U. Schmitt, A. K. Louis, C. Wolters, and M. Vaukhonen, “Efficient algorithms for the regularizationof dynamic inverse problems: II. Applications,”
Inverse Problems , vol. 18, no. 3, pp. 659–676, 2002.[104] T. Schuster, B. Hahn, and M. Burger, “Dynamic inverse problems: Modelling - Regularization -numerics,”
Inverse Problems , vol. 34, no. 4, 2018.[105] V. Estellers, J. Thiran, and X. Bresson, “Enhanced Compressed Sensing Recovery With Level SetNormals,”
IEEE Transactions on Image Processing , vol. 22, pp. 2611–2626, 7 2013.[106] D. Kazantsev, W. R. B. Lionheart, P. J. Withers, and P. D. Lee, “Multimodal Image Reconstruc-tion using Supplementary Structural Information in Total Variation Regularization,”
Sensing andImaging , vol. 15, p. 97, 1 2014.[107] J. Rasch, E.-M. Brinkmann, and M. Burger, “Joint Reconstruction via Coupled Bregman Iterationswith Applications to PET-MR Imaging,”
Inverse Problems , vol. 34, no. 1, p. 014001, 2018.[108] V. Estellers, S. Soatto, and X. Bresson, “Adaptive Regularization With the Structure Tensor,”
IEEE Transactions on Image Processing , vol. 24, no. 6, pp. 1777–1790, 2015.[109] P. Song, X. Deng, J. F. C. Mota, N. Deligiannis, P.-L. Dragotti, and M. Rodrigues, “MultimodalImage Super-resolution via Joint Sparse Representations induced by Coupled Dictionaries,”
IEEETransactions on Computational Imaging , pp. 1–1, 2019.[110] V. Caselles, B. Coll, and J.-M. Morel, “Geometry and Color in Natural Images,”
Journal of Math-ematical Imaging and Vision , vol. 16, no. Section 2, pp. 89–105, 2002.[111] R. Kimmel, R. Malladi, and N. Sochen, “Images as Embedded Maps and Minimal Surfaces:Movies, Color, Texture, and Volumetric Medical Images,”
International Journal of Computer Vi-sion , vol. 39, no. 2, pp. 111–129, 2000.[112] J. A. Fessler, I. Elbakri, P. Sukovic, and N. H. Clinthorne, “Maximum-likelihood dual-energytomographic image reconstruction,” in
SPIE: Medical Imaging , vol. 4684, pp. 1–25, 2002.[113] B. Heismann, B. Schmidt, and T. Flohr,
Spectral Computed Tomography . SPIE Press, 2012.[114] Y. Long and J. A. Fessler, “Multi-material decomposition using statistical image reconstructionfor spectral CT,”
IEEE Transactions on Medical Imaging , vol. 33, no. 8, pp. 1614–1626, 2014.[115] F. Lenzen and J. Berger, “Solution-Driven Adaptive Total Variation Regularization,” in
SSVM ,pp. 203–215, 2015. 28116] M. Hinterm¨uller and M. M. Rincon-Camacho, “Expected Absolute Value Estimators for a Spa-tially Adapted Regularization Parameter Choice Rule in L1-TV-based Image Restoration,”
InverseProblems , vol. 26, no. 8, p. 085005, 2010.[117] Y. Dong, M. Hinterm¨uller, and M. M. Rincon-Camacho, “Automated Regularization ParameterSelection in Multi-Scale Total Variation Models for Image Restoration,”
Journal of MathematicalImaging and Vision , vol. 40, no. 1, pp. 82–104, 2011.[118] K. Bredies, Y. Dong, and M. Hinterm¨uller, “Spatially Dependent Regularization Parameter Se-lection in Total Generalized Variation Models for Image Restoration,”
International Journal ofComputer Mathematics , pp. 1–15, 2012.[119] E. Esser, X. Zhang, and T. F. Chan, “A General Framework for a Class of First Order Primal-Dual Algorithms for Convex Optimization in Imaging Science,”
SIAM Journal on Imaging Sciences ,vol. 3, no. 4, pp. 1015–1046, 2010.[120] A. Chambolle and T. Pock, “A First-Order Primal-Dual Algorithm for Convex Problems withApplications to Imaging,”
Journal of Mathematical Imaging and Vision , vol. 40, no. 1, pp. 120–145, 2011.[121] H. H. Bauschke and P. L. Combettes,
Convex Analysis and Monotone Operator Theory in HilbertSpaces . 2011.[122] P. L. Combettes and J. C. Pesquet, “Proximal splitting methods in signal processing,”
SpringerOptimization and Its Applications , vol. 49, pp. 185–212, 2011.[123] N. Parikh and S. P. Boyd, “Proximal Algorithms,”
Foundations and Trends in Optimization , vol. 1,no. 3, pp. 123–231, 2014.[124] A. Chambolle and T. Pock, “An Introduction to Continuous Optimization for Imaging,”
ActaNumerica , vol. 25, pp. 161–319, 2016.[125] T. Pock and A. Chambolle, “Diagonal Preconditioning for First Order Primal-Dual Algorithms inConvex Optimization,” in
Proceedings of the IEEE International Conference on Computer Vision ,pp. 1762–1769, 2011.[126] A. Chambolle, M. J. Ehrhardt, P. Richt´arik, and C.-B. Sch¨onlieb, “Stochastic Primal-Dual Hy-brid Gradient Algorithm with Arbitrary Sampling and Imaging Applications,”
SIAM Journal onOptimization , vol. 28, no. 4, pp. 2783–2808, 2018.[127] J. Adler, H. Kohr, and O. ¨Oktem, “Operator Discretization Library (ODL),” 1 2017.[128] W. van Aarle, W. J. Palenstijn, J. De Beenhouwer, T. Altantzis, S. Bals, K. J. Batenburg, andJ. Sijbers, “The ASTRA Toolbox: A platform for advanced algorithm development in electrontomography,”
Ultramicroscopy , vol. 157, pp. 35–47, 2015.[129] W. van Aarle, W. J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. De Been-houwer, K. Joost Batenburg, and J. Sijbers, “Fast and flexible X-ray tomography using the ASTRAtoolbox,”