Tomography: mathematical aspects and applications Paolo Facchi1,2, Marilena Ligab`o3, Sergio Solimini3 1DipartimentodiFisicaandMECENAS,Universit`adiBari,I-70126Bari,Italy 5 2INFN,SezionediBari,I-70126Bari,Italy 1 3DipartimentodiMeccanica,MatematicaeManagement,PolitecnicodiBari, I-70125Bari,Italy 0 2 n Abstract. In this article we present a review of the Radon transform and the instability of the tomographic reconstruction process. We show some new a J mathematical results in tomography obtained by a variational formulation of the reconstruction problem based on the minimization of a Mumford-Shah type 9 functional. Finally,weexhibitaphysicalinterpretationofthisnewtechniqueand 2 discusssomepossiblegeneralizations. ] h p PACSnumbers: 03.65.Wj,42.30.Wb,02.30.Uu - h t a m Keywords: Radon transform; integral geometry; image segmentation; calculus of variations. [ 1 v 5 8 3 7 0 . 1 0 5 1 : v i X r a Tomography: mathematical aspects and applications 2 Tomogram Tomogram Figure 1. Distributionandtwotomograms. 1. Introduction Figure 2. Reconstructionproblemsfromdiversefieldsmaybe The primary goal of tomography is to determine the unitedwithintheframeworkoftomography. internal structure of an object without cutting it, namelyusingdataobtainedbymethodsthatleavethe objectunderinvestigationundamaged. Thesedatacan its square root is defined by Fourier transform (see be obtained by exploiting the interaction between the Theorem 1). We now observe that the formula above object and various kinds of probes including X-rays, has built in a remarkable duality: first one integrates electrons, and many others. After its interaction with overthesetofpointsinaline,thenoneintegratesover the object under investigation, the probe is detected the set of lines passing through a given point. This to produce what we call a projected distribution or formulacanbeextendedtotheN-dimensionalcaseby tomogram, see Fig. 1. computingtheintegralsofthefunctionf onallpossible Tomography is a rapidly evolving field for its hyperplanes. This suggests to consider the transform broad impact on issues of fundamental nature and for f f defined as follows. If f is a function on RN (cid:55)→ R its important applications such as the development of then f is the function defined on the space of all R diagnostic tools relevant to disparate fields, such as possible (N 1)-dimensional planes in RN such that, − lunedì 3 novembre 14 engineering, biomedical and archaeometry. Moreover, given a hyperplane λ, the value of f(λ) is given by R tomography can be a powerful tool for many the integral of f along λ. The function f is called R reconstruction problems coming from many areas of Radon transform of f. research, such as imaging, quantum information and There exist several important generalizations of computation, cryptography, lithography, metrology the Radon transform by John [2], Gel’fand [3], and many others, see Fig. 2. Helgason [4] and Strichartz [5]. More recent analysis From the mathematical point of view the has been boosted by Margarita and Volodya Man’ko reconstruction problem can be formulated as follows: and has focused on symplectic transforms [6], on the one wants to recover an unknown function through deep relationship with classical systems and classical the knowledge of an appropriate family of integral dynamics [7, 8, 9], on the formalism of star product transforms. It was proved by J. Radon [1] that a quantization[10,11,12],andonthestudyofmarginals smooth function f(x,y) on R2 can be determined along curves that are not straight lines [13, 14]. explicitly by means of its integrals over the lines in In quantum mechanics the Radon transform of R2. Let f(X,θ) denote the integral of f along the the Wigner function [15, 16, 17] was considered in the R line xcosθ+ysinθ =X (tomogram). Then tomographic approach to the study of quantum states (cid:90) 2π dθ [18, 19] and experimentally realized with different 1 f(x,y)=( ∆)2 f(xcosθ+ysinθ,θ) , (1) particles and in diverse situations. For a review − R 4π 0 on the modern mathematical aspects of classical and where ∆ = ∂2 + ∂2 is the Laplacian on R2, and quantumtomographysee[20]. Goodreviewsonrecent ∂x2 ∂y2 Tomography: mathematical aspects and applications 3 tomographic applications can be found in [21] and contains a short introduction to the Radon transform, in[22],whereparticularemphasisisgivenonmaximum its dual map and the inversion formula. Section 3 is likelihood methods, that enable one to extract the devotedtoabriefdiscussionontheill-posednessofthe maximum reliable information from the available data tomographic reconstruction and to the introduction of can be found. regularization methods. In Section 4 a MS functional As explained above, from the mathematical point is applied to tomography as a regularization method. of view, the internal structure of the object is In particular, in Subsection 4.1 the piecewise constant described by an unknown function f (density), that is model and known results are discussed together connected via an operator to some measured quantity with a short list of some interesting open problems. g (tomograms). The tomographic reconstruction Finally, in Section 5 we present an electrostatic problem can be stated as follows: for given data interpretation of the regularization method based on g, the task is to find f from the operator equation the MS functional, which motivates us to introduce an f = g. There are many problems related to the improved regularization method, based on the Blake- R implementation of effective tomographic techniques Zisserman functional [29], as a relaxed version of the due to the instability of the reconstruction process. previous one. There are two principal reasons of this instability. The first one is the ill-posedness of the reconstruction 2. The Radon transform: definition and problem: in order to obtain a satisfactory estimate inversion formula of the unknown function it is necessary an extremely preciseknowledgeofitstomograms,whichisingeneral Consider a body in the plane R2, and consider a beam physically unattainable [23]. The second reason is the of particles (neutrons, electrons, X-rays, etc.) emitted discrete and possibly imperfect nature of data that by a source. Assume that the initial intensity of the allowstoobtainonlyanapproximationoftheunknown beam is I . When the particles pass through the body 0 function. theyareabsorbedorscatteredandtheintensityofthe Thefirstquestioniswhetherapartialinformation beam traversing a length ∆s decreases by an amount still determines the function uniquely. A negative proportional to the density of the body µ, namely answer is given by a theorem of Smith, Solomon ∆I/I = ∆sµ(s), (2) and Wagner [24], that states: “A function f with − compact support in the plane is uniquely determined so that by any infinite set, but by no finite set of its (cid:18) (cid:90) s (cid:19) I(s)=I exp µ(r)dr . (3) tomograms”. Therefore, it is clear that one has to 0 − 0 abandon the request of uniqueness in the applications Adetectorplacedattheexitofthebodymeasuresthe of tomography. Thus, due to the ill-posedness of final intensity I(s) and then from reconstruction problem and to the loss of uniqueness in the inversion process, a regularization method has I(s) (cid:90) s ln = µ(r)dr (4) to be introduced to stabilize the inversion. − I0 0 A powerful approach is the introduction of one can record the value of the density µ integrated a Mumford-Shah (MS) functional, first introduced on a line. If another ray with a different direction is in a different context for image denoising and considered, with the same procedure one obtains the segmentation [25]. The main motivation is that, in value of the integral of the density on that line. many practical applications, one is not only interested Themathematicalmodeloftheabovesetupisthe in the reconstruction of the density distribution f, following: Given a smooth function f(x) on the plane, but also in the extraction of some specific features or x R2, and a line λ, consider its tomogram, given by patterns of the image. An example is the problem of ∈ (cid:90) the determination of the boundaries of inner organs. f(λ)= f(x)dm(x), (5) R By minimizing the MS functional, one can find not λ only (an approximation of) the function but also its where dm is the Euclidean measure on the line λ. In sharp contours. Very recently a MS functional for this way, we have defined an operator that maps a R applicationstotomographyhasbeenintroducedinthe smooth function f on the plane R2 into a function f R literature [26, 27, 28]. Some preliminary results in this on P2, the manifold of the lines in R2. context are already available but there are also many We ask the following question: If we know the interestingopenproblemsandpromisingresultsinthis family of tomograms ( f(λ))λ∈P2, can we reconstruct R direction, as we will try to explain in the second part the density function f? The answer is affirmative and of this article. in the following we will see how to obtain this result. The article is organized as follows. Section 2 Let us generalize the above definitions to the case ofanN-dimensionalspace. Letf beafunctiondefined Tomography: mathematical aspects and applications 4 N n R R potential(f]()[(fx))(x) n R I R x � aanN x y y � | � | charge distributionf(y) X Figure 4. I(Rf)(x) is the potential at x generated by the chargedistributionf. onRN. Inanalogywith (RN)wedefine (R SN−1) S S × lunedì 27 ottobre 14 as the space of C∞ functions g on R SN−1 which for giovedì 30 ottobre 14any integers m 0, any multiindex ×α NN, and any ≥ ∈ differential operator D on SN−1 satisfy (cid:12)(cid:12) ∂α(Dg) (cid:12)(cid:12) Figure 3. Parametrizationofthehyperplaneλusingitssigned sup (cid:12)X m (X,ξ)(cid:12)<+ . (9) distanceX fromtheoriginandaunitvectorξ perpendicularto X∈R,ξ∈SN−1(cid:12)| | ∂Xα (cid:12) ∞ λ. The space (PN) is then defined as the set of g S ∈ (R SN−1) satisfying g( X, ξ)=g(X,ξ). S × − − on RN, integrable on each hyperplane in RN and let Now we want to obtain an inversion formula, PN be the manifold of all hyperplanes in RN. The namely we want to prove that one can recover a lunedì 27 ottobre 14 Radon transform of f is defined by Eq. (5), where dm function f on RN from the knowledge of its Radon istheEuclideanmeasureonthehyperplaneλ. Thuswe transform. In order to get this result we need a have an operator , the Radon transform, that maps preliminary lemma, whose proof can be found in [20], R a function f on RN into a function f on PN, namely which suggests an interesting physical interpretation. R f f. Itsdualtransform,alsocalledbackprojection Lemma1. Letf (RN)andV(x)=1/x,x RN, op(cid:55)→eraRtor, g g associates to a function g on PN the x=0 . Then ∈S | | ∈ function Ig(cid:55)→onIRN given by ((cid:54) f)(x)=aNf V, (10) (cid:90) I R ∗ g(x)= g(λ)dµ(λ), (6) where aN depends only on the dimension N, and I x∈λ denotes the convolution product, ∗ (cid:90) where dµ is the unique probability measure on the (f g)(x)= f(y)g(x y)dy. (11) compact set {λ ∈ PN|x ∈ λ} which is invariant under ∗ RN − the group of rotations around x. A physical interpretation of Lemma 1 is the Let us consider the following covering of PN following: if f is a charge distribution, then the R SN−1 PN, (X,ξ) λ, (7) potential at the point x generated by that charge is × → (cid:55)→ exactly ( f(x)),seeFig.4. Notice,however,thatthe where SN−1 is the unit sphere in RN. Thus, the I R potential of a point charge scales always as the inverse equation of the hyperplane λ is distance independently of the dimension N, and thus λ= x RN X ξ x=0 , (8) it is Coulomb only for N = 3. The only dependence { ∈ | − · } on N is in the strength of the elementary charge a . with a b denoting the Euclidean inner product of N · This fact is crucial: indeed, the associated Poisson a,b RN. See Fig. 3. Observe that the pairs ∈ equationinvolvesanN-dependent(fractional)powerof (X,ξ),( X, ξ) R SN−1 are mapped into the − − ∈ × the Laplacian, which appears in the inversion formula same hyperplane λ PN. Therefore (7) is a ∈ for the Radon transform. double covering of PN. Thus PN has a canonical manifold structure with respect to which this covering Theorem 1. Let f (RN). Then ∈S mapping is differentiable. We identify continuous 1 N−1 (differentiable) functions g on PN with continuous f(x)= 2(2π)N−1(−∆) 2 I(Rf)(x) (12) (differentiable) functions g on R SN−1 satisfying where ( ∆)α, with α > 0, is a pseudodifferential g(X,ξ)=g( X, ξ). × operator−whose action is − − We will momentarily work in the Schwartz space (cid:90) dk (RN) of complex-valued rapidly decreasing functions (( ∆)αf)(x)= k 2αfˆ(k)eik·x , (13) S − RN | | (2π)N Tomography: mathematical aspects and applications 5 where fˆis the Fourier transform of f, reason, regularization methods have to be introduced (cid:90) to stabilize the inversion in the presence of data noise. fˆ(k)= f(x)e−ik·x dx. (14) We discuss ill-posed problems only in the RN framework of linear problems in Hilbert spaces [30]. The proof of Theorem 1 can be found in [20]. Let H,K be Hilbert spaces and let A be a linear Equation (12) says that, modulo the final action of bounded operator from H into K . The problem N−1 (−∆) 2 , the function f can be recovered from its given g K , find f H such that Af =g (17) Radon transform f by the application of the dual ∈ ∈ mapping : first oRne integrates over the set of points is called well-posed by Hadamard (1932) if it is inahyperIplaneandthenoneintegratesoverthesetof uniquely solvable for each g K and if the solution ∈ hyperplanes passing through a given point. Explicitly depends continuously on g. Otherwise, (17) is called we get ill-posed. Thismeansthatforanill-posedproblemthe (cid:90) operatorA−1 eitherdoesnotexist,orisnotdefinedon 1 f(x)= ( ∆)N2−1 f(ξ x,ξ)dξ, (15) all of K , or is not continuous. The practical difficulty 2(2π)N−1 − SN−1R · with an ill-posed problem is that even if it is solvable, which has the following remarkable interpretation. thesolutionofAf =gneednotbeclosetothesolution Note that if one fixes a direction ξ SN−1, then of Af =g(cid:15) if g(cid:15) is close to g. ∈ the function f(ξ x,ξ) is constant on each plane In general A−1 is not a continuous operator. R · perpendiculartoξ,i.e.itisa(generalized)planewave. To restore continuity we introduce the notion of a Therefore,Eq.(15)givesarepresentationoff interms regularizationofA−1. Thisisafamily(T ) oflinear γ γ>0 of a continuous superposition of plane waves. A well- continuous operators T : K H which are defined γ known analogous decomposition is given by Fourier on all K and for which → transform. When N = 2, one recovers the inversion lim T g =A−1g (18) γ formula (1) originally found by Radon [1]. γ→0 on the domain of A−1. Obviously T + as γ (cid:107) (cid:107) → ∞ 3. Instability of the inversion formula with γ 0 if A−1 is not bounded. With the help of a → noisy data regularization we can solve (17) approximately in the following sense. Let g(cid:15) K be an approximation to ∈ We have defined the Radon transform of any function g such that g g(cid:15) (cid:15) . Let γ((cid:15)) be such that, as f (RN) as f. The following theorem [4] contains (cid:15) 0, (cid:107) − (cid:107) ≤ ∈S R → the characterization of the range of the Radon linear γ((cid:15)) 0, T (cid:15) 0. (19) operator and the extension of to the space of → (cid:107) γ((cid:15))(cid:107) → square intRegrable functions L2(RN)R. Then, as (cid:15)→0, T g(cid:15) A−1g T (g(cid:15) g) Theorem 2. The Radon transform is a linear one- (cid:107) γ((cid:15)) − (cid:107)≤(cid:107) γ((cid:15)) − (cid:107) to-one mapping of S(RN) onto SHR(PN), where the +(cid:107)(Tγ((cid:15))−A−1)g(cid:107) space H(PN) is defined as follows: g H(PN) if T (cid:15)+ (T A−1)g 0. and onSly if g (PN) and for any integ∈erSk N the ≤(cid:107) γ((cid:15))(cid:107) (cid:107) γ((cid:15))− (cid:107)→ ∈ S ∈ (20) integral (cid:90) Hence, Tγ((cid:15))g(cid:15) is close to A−1g if g(cid:15) is close to g. g(X,ξ)Xk dX (16) The number γ is called a regularization parameter. R Determining a good regularization parameter is one of is a homogeneous polynomial of degree k in ξ1,...,ξN. the crucial points in the application of regularization Moreover, the Radon operator can be extended to a methods. R continuous operator from L2(RN) and L2(PN). There are several methods for constructing In medical imaging, computerized tomography a regularization as the truncated singular value is a widely used technique for the determination of decomposition, the method of Tikhonov-Phillips or the density f of a sample from measurements of some iterative methods [30]. In the following section the attenuation of X-ray beams sent through the we present a regularization method based on the material along different angles and offsets. The minimization of a Mumford-Shah type functional. measured data g are connected to the density f via the Radon Transform . To compute the density 4. Mumford-Shah functional for the R distributionf theequationg = f hastobeinverted. simultaneous segmentation and reconstruction R Unfortunately it is a well known fact that is not of a function R continuously invertible on L2(PN) [4], and this imply that the problem of inversion is ill-posed. For this In many practical applications one is not only interested in the reconstruction of the density Tomography: mathematical aspects and applications 6 distributionf butalsointheextractionofsomespecific featureswithin theimage whichrepresentsthe density distribution of the sample. For example, the planning of surgery might require the determination of the boundaries of inner organs like liver or lung or the separationofcancerousandhealthytissue. Segmenting a digital image means finding its homogeneous regions and its edges, or boundaries. Of course, the Figure 5. Left: image of an eye (g). Center: contours of the homogeneous regions are supposed to correspond to image in the Mumford-Shah model (edges Γ). Right: piecewise meaningful parts of objects in the real world, and smoothfunctionapproximatingtheimage(cartoonu)[32]. the edges to their apparent contours. The Mumford- Shah variational model is one of the principal models is small, a lot of edges are allowed and we get a fine of image segmentation. It defines the segmentation segmentation. As β increases, the segmentation gets problem as a joint smoothing/edge detection problem: coarser. Formoredetailsonthemodelseetheoriginal given an image g(x), one seeks simultaneously a paper [25], and the book [31]. “piecewise smoothed image” u(x) with a set Γ of The minimization of the MS functional in (21) is abrupt discontinuities, the “edges” of g. The original performed among the admissible pairs (Γ,u) such that Mumford-Shah functional [25], is the following: Γ is closed and u C1(D Γ). It is worth noticing (cid:90) ∈ \ MS(Γ,u)= u g 2 +α u(x)2 dx that in this model there are two unknowns: a scalar (cid:107) − (cid:107)L2(D) |∇ | function u and the set Γ of its discontinuities. For this D\Γ +β N−1(D Γ), (21) reason this category of problems is often called “Free H ∩ Discontinuities Problem”. Existence of minimizers of where the MS functional in (21) was proven by De Giorgi, D RN is an open set (screen); Carriero, Leaci in [33] in the framework of bounded • ⊂ variation functions without Cantor part (space SBV) Γ RN is a closed set (set of edges); • ⊂ introducedbyAmbrosioandDeGiorgiin[34]. Further u:D R (cartoon); regularity properties for optimal segmentation in the • → u denotes the distributional gradient of u; Mumford-Shah model were shown in [35, 36, 37, 38]. • ∇ g L2(D) is the datum (digital image); Here we present a variation of the MS functional, • ∈ adapted to the inversion problem of the Radon α,β >0 are parameters (tuning parameters); • transform. Moreprecisely,weconsideraregularization N−1 denotes the (N 1)-dimensional Hausdorff method that quantifies the edge sets together with • H − measure. images, i.e. a procedure that gives simultaneously The squared L2 distance in (21) plays the role a reconstruction and a segmentation of f (assumed of a fidelity term: it imposes that the cartoon u to be supported in D RN) directly from the ⊂ approximate the image g. The second term in the measured tomograms g, based on the minimization of functional imposes that the cartoon u be piecewise the Mumford-Shah type functional smooth outside the edge set Γ. In other word this J (Γ,f)= f g 2 term favors sharp contours rather than zones where a MS (cid:107)R − (cid:107)L2(PN) (cid:90) thin layer of gray is used to pass smoothly from white +α f(x)2 dx+β N−1(Γ). (22) to black or viceversa. Finally the third term in the D\Γ|∇ | H functional imposes that the contour Γ be “small” and The only difference between the functionals MS and as smooth as possible. What is expected from the J is the first term, i.e. the fidelity term, that MS minimization of this functional is a sketchy, cartoon- ensures that the reconstruction for f is close enough like version of the given image together with its to a solution of the equation f = g, whereas the contours. See Fig. 5. other terms play exactly the saRme role explained for TheminimizationoftheMS functionalrepresents thefunctionalMS. Asexplainedabove, inadditionto a compromise between accuracy and segmentation. the reconstruction of the density f, we are interested The compromise depends on the tuning parameters in the reconstruction of its singularity set Γ, i.e. the α and β which have different roles. The parameter setofpointswherethesolutionf isdiscontinuous. The α determines how much the cartoon u can vary, if α maindifferencewithrespecttothestandardMumford- is small some variations of u are allowed, while as α Shah functional (21) is that we have to translate the increases u tends to be a piecewise constant function. information about the set of sharp discontinuities of g The parameter β represents a scale parameter of the (and hence on the space of the Radon transform) into functional and measure the amount of contours: if β information about the strong discontinuities of f. Tomography: mathematical aspects and applications 7 4.1. The piecewise constant model Ω(apartitionofthedomainDwithatmostmdistinct regions satisfying the non degeneracy condition (23)). Here we will review the results obtained by Ramlau So the problem is to find f˜ PC (D) such that and Ring [26] concerning the minimization of (22) ∈ m m restricted to piecewise constant functions f, and then (cid:88) f˜= f˜χ PC (D), (26) consider some interesting open problems. For medical k Ω˜k ∈ m k=1 applications, it is often a good approximation to where restrict the reconstruction to densities f that are constant with respect to a partition of the body, as (Ω˜,f˜)=arg min J (Γ,f). (27) β the tissues of inner organs, bones, or muscles have (Ω,f) approximately constant density. It is clear that f˜ will depend on the regularization We introduce the space PCm(D) as the space parameter β and on the error level (cid:15). of piecewise constant functions that attain at most Now we can state the results concerning the m different function values, where D is an open and functional J in (25). There are several technical β bounded subset of RN. In other words, each f details necessary for the precise statement and proof ∈ PCm(D) is a linear combination of m characteristic of the theorems, for which we refer to the original functions χΩk of sets (Ωk)k=1,...,m which satisfy paper[26]. Herewewillgiveasimplifiedversionofthe theorems with the purpose of explain the main goal, m (cid:88) χ =χ a.e. without too many technical details. The first result is Ωk D about the existence of minimizers of the functional J k=1 β in (25). We assume that the Ω ’s are open relatively to D and k we set Γ = ∂Ω for the boundary of Ω with respect Theorem 3. For all g(cid:15) L2(PN) there exists a to the tokpologykrelative to the open dkomain D. In minimizer (Ω˜(cid:15)β,f˜β(cid:15)) of the fu∈nctional Jβ in (25), with this situation the edge set will be given by the union β >0. of the boundaries of Ω ’s. For technical reasons it is k The second result regards the stable dependence necessary to assume a nondegeneracy condition on the of the minimizers of the functional J in (25) on the β admissible partitions of D: error level (cid:15). Ω=(Ωk)k=1,...,m is admissible if LN(Ωk)≥δ, (23) Theorem 4. Let (g(cid:15)n)n∈N be a sequence of functions forsomeδ >0, forallk =1,...,m, where N denotes in L2(PN) and let g(cid:15) L2(PN). For all n N, the Lebesgue measure on RN. L let (Ω˜(cid:15)n,f˜(cid:15)n) denote the∈minimizers of the funct∈ional β β It turns out to be convenient to split the Jβ with initial data g(cid:15)n. If g(cid:15)n g(cid:15) in L2(PN), as information encoded in a typical function, n + , then there exists a subse→quence of (Ω˜(cid:15)n,f˜(cid:15)n) (cid:88)m su→ch th∞at β β f =k=1fkχΩk ∈PCm(D), (24) (Ω˜β(cid:15)nj,f˜β(cid:15)nj)→(Ω˜(cid:15)β,f˜β(cid:15)) into a “geometrical” part described by the m-tuple of as j + , and (Ω˜(cid:15),f˜(cid:15)) is a minimizer of J with → ∞ β β β pairwise disjoint sets Ω = (Ω ) which cover initial data g(cid:15). Moreover, the limit of each convergent k k=1,...,m D up to a set of measure zero and a “functional” subsequence of (Ω˜(cid:15)n,f˜(cid:15)n) is a minimizer of J with β β β part given by the m-tuple of values f = (f ) . initial data g(cid:15). k k=1,...,m We also use the notation Γ = (Γ ) , for the k k=1,...,m FinallythelastTheoremisaregularizationresult. boundaries Γ =∂Ω of Ω . k k k As usual when dealing with inverse problems, Theorem 5. Let f∗ PCm(D) be given, ∈ we have to assume that the data g are not exactly m ktonmowognr,ambustg(cid:15)thoaftaw(heypaorethoetnilcyal)gievxeanctndoaistyasmeteagsuwriethd f∗ =(cid:88)fk∗χΩ∗k, k=1 (cid:107)g(cid:15)−Ifgw(cid:107)Le2(rPeNst)r≤ict(cid:15).the functional (22) to functions in and let g∗ = f∗. Assume we have noisy data R PCm(D) we obtain that the second term (involving g(cid:15) ∈ L2(PN) with (cid:107)g(cid:15) − g∗(cid:107)L2(PN) ≤ (cid:15). Let us the derivatives of f) disappears, therefore it remains choosetheparameterβ =β((cid:15))satisfyingtheconditions to minimize the functional β((cid:15)) 0 and (cid:15)2/β((cid:15)) 0 as (cid:15) 0. For any m sequen→ce (cid:15) 0, let (Ω˜n→,f˜n) denote→the minimizers Jβ(Ω,f)=(cid:107)Rf −g(cid:15)(cid:107)2L2(PN)+β(cid:88)HN−1(Γk), (25) of the funnct→ional Jβ((cid:15)n) with initial data g(cid:15)n and k=1 regularization parameter β = β((cid:15)n). Then there exists overPC (D),withrespecttothefunctionalvariablef a convergent subsequence of (Ω˜n,f˜n). Moreover, for m (avectorofmcomponents)andthegeometricvariable Tomography: mathematical aspects and applications 8 every convergent subsequence with limit (Ω˜,f˜) the Proof. We know that f = 1 ∆ ( f), therefore −2(2π)2 I R function (cid:90) (cid:90) f˜=(cid:88)m f˜kχΩ˜k ∈PCm(D) R×S2|Rf(X,ξ)|2 dX dξ = R3f(x)I(Rf)(x)dx k=1 (cid:90) 1 = [ ∆ ( f)](x) ( f)(x)dx is a solution of the equation Rf = g∗ with a minimal R3 2(2π)2 − I R I R perimeter. Moreoveriff∗ istheuniquesolutionofthis (cid:90) 1 equation then the whole sequence converges = ( f)(x)2 dx. 2(2π)2 R3|∇I R | f˜n f∗, in L2(D), → when n + . Lemma 3. For all real valued f (R3) define ∈S → ∞ 1 Finally, let us list some open problems in this E = ( f) (31) −2(2π)2∇I R context: and Is the nondegeneracy condition (23) necessary? • 1 Can one find an a priori optimal value for the ϕ= ( f). (32) • 2(2π)2I R number m of different values? Then Is it possible to give an a priori estimate on the • L∞-norm of the solution (maximum principle)? f = E = ∆ϕ. (33) ∇· − And finally, it would be very important for Proof. • applications to prove the existence of minimizers 1 of the functional J not restricted to piecewise MS E = ( f)= ∆ϕ constant functions f. ∇· − 2(2π)2∇·∇I R − 1 We observe that all these problems are quite natural, = ( ∆) ( f)=f, 2(2π)2 − I R and have been completely solved in the case of the standard Mumford-Shah functional MS in (21), see where we used the inversion formula (28). e.g. [31, 39]. Now we consider a measured tomogram g : P3 → 5. Electrostatic interpretation of JMS R and let us assume that g =Rf0 for some f0 :R3 → R. By Lemma 2-3 it follows immediately that the In this section we restrict our attention to the fidelityterm(cid:107)Rf−g(cid:107)2L2(P3)canberewrittenasfollows: 3-dimensional case. We propose an electrostatic f g 2 =2(2π)2 E E 2 interpretation of the regularization method based on (cid:107)R − (cid:107)L2(P3) (cid:107) − g(cid:107)L2(R3) the functional JMS discussed in the previous section. =2(2π)2(cid:107)∇ϕ−∇ϕg(cid:107)2L2(R3), (34) The intent is to give a physical explanation of the where fidelityterm f g 2 inthefunctional(22),that provide the (cid:107)inRtui−tion(cid:107)Lf2o(rP3)an improved regularization E =−2(21π)2∇I(Rf), Eg =−2(21π)2∇I(g) (35) method. ForN =3,theinversionformula(12)andthe are the corresponding electric fields, while electrostaticidentity(10)particularize,respectively,as follows: for all f (R3) one gets 1 1 1 ∈S ϕ= 2(2π)2I(Rf), ϕg = 2(2π)2I(g) (36) f = ∆ ( f) (28) −2(2π)2 I R are the corresponding potentials. With respect to the standard Mumford-Shah and functional MS in (21), the new fidelity term in the I(Rf)=a3f ∗V, (29) functional JMS in (22) controls the distance between where V(x) = 1/x and a is a constant. We present the Radon transform of f and the tomographic data 3 two preliminary L|em| mas. g. The relevant difference with respect to the original functional is that the function f and its Radon Lemma 2. For all real valued f (R3) one has transform f are defined in different spaces. Let ∈S R 1 us try to interpret the fidelity term f g 2 f 2 = ( f) 2. (30) (cid:107)R − (cid:107)L2(P3) (cid:107)R (cid:107) 2(2π)2(cid:107)∇I R (cid:107) from a physical point of view. A key ingredient for this goal is the electrostatics formulation of the Radon transform. This formulation can be summarized as follows: if we consider, in dimension 3, a function f, Tomography: mathematical aspects and applications 9 we can think at it as a charge distribution density; if among admissible triplets (Γ ,Γ ,v), where 0 1 we apply to f first the Radon operator and then its R D RN is an open set; adjoint weobtain,uptoaconstant,theelectrostatic • ⊂ potentiaIl generated by the charge distribution f. This Γ0,Γ1 RN are closed sets; • ⊂ formulation can be stated in any dimension N: the Γ0 isthesetofdiscontinuitiesofv (jumpset),and • difference with general potential theory in dimension Γ1 thethesetofdiscontinuitiesof v (creaseset); ∇ N is that, in tomography, the potential produced by a v : D R, v C2(D (Γ0 Γ1)) C(D Γ0) is point charge always scales like 1/x, which is the case • a scala→r functi∈on; \ ∪ ∩ \ | | ofelectrostaticpotentialonlyindimension3. Fromthe ∆v denotes the distributional Laplacian of v; electrostatic formulation of the Radon transform we • v L2(D) is the datum (grey intensity levels of can prove that the fidelity term in the functional J 0 MS • ∈ the given image); actuallyimposesthattheelectricfieldproducedbythe charge distribution f must be close to the “measured α,β,γ >0 are parameters; • electric field”. Therefore we conclude that the term N−1 denotes the (N 1)-dimensional Hausdorff (cid:107)Rf −g(cid:107)2L2(P3) is a fidelity term in this weaker sense. • Hmeasure. − Using this property based on the electrostatic The Blake-Zisserman functional allows a more precise interpretation of the tomographic reconstruction, we segmentation than the Mumford-Shah functional in can try to minimize some appropriate functionals in the sense that also the curvature of the edges of the the new variables E (electric field) or ϕ (electric original picture is approximated. On the other hand, potential) and then compute the corresponding f minimizers may not always exist, depending on the (charge density). We manipulate the functional J MS values of the parameters β,γ and on the summability as follows: assumption on v . We refer to [29, 40, 41, 42, 43, 44] 0 JMS(Γ,f) for motivation and analysis of variational approach to (cid:90) image segmentation and digital image processing. In = f g 2 +α f(x)2 dx+β 2(Γ) (cid:107)R − (cid:107)L2(P3) R3\Γ|∇ | H particular see [40, 41, 45] for existence of minimizer (cid:90) results and [42] for a counterexample to existence =2(2π)2 E Eg 2+α ( E)(x)2 dx and [46, 47] for results concerning the regularity of (cid:107) − (cid:107) R3\Γ|∇ ∇· | minimizers. +β 2(Γ) Equation (37) implies that the functional J H (cid:90) MS can be rewritten in terms of the vector field E and =2(2π)2 E E 2+α ∆E 2 dx+β 2(Γ) (cid:107) − g(cid:107) R3\Γ| | H of the discontinuities set of ∇ · E, i.e. the set of creases of E, using the terminology of the Blake- =F(Γ,E) (37) Zisserman model. The fact that in the functional F where F is a new functional depending on a vector the discontinuities set of E is not present depends on function E and on a set Γ, and we used the fact that the fact that we are assuming that the charge density E =0, since E is conservative. f inthefunctionalJ donotconcentrateonsurfaces ∇∧ MS We observe that the functional or on lines. If we admit concentrated charge layers we (cid:90) F(Γ,E)=2(2π)2 E E 2+α ∆E 2 dx can consider the Blake-Zisserman model for the vector g (cid:107) − (cid:107) R3\Γ| | function E as a relaxed version of the Mumford-Shah +β 2(Γ) (38) model for the charge f. In other words we propose to H investigatetheconnectionsbetweenminimizersofJ MS is a second order functional for a vector field E and minimizers of the higher order functional J : BZ in which appears the measure of the set Γ that is J (Γ ,Γ ,E)= E E 2 the set of discontinuities of f and thus is the set BZ 0 1 (cid:107) − g(cid:107)L2(R3) (cid:90) of discontinuities of E. In the functional F ∇ · +α ∆E(x)2 dx we recognize some similarities with a famous second- R3\(Γ0∪Γ1)| | order free-discontinuity problem: the Blake-Zisserman +β 2(Γ ) 0 model. Thismodelisbasedontheminimizationofthe H +γ 2(Γ Γ ), (40) Blake-Zisserman functional 1 0 H \ BZ(Γ ,Γ ,v)= v v 2 with the additional constraint E = 0. The 0 1 (cid:107) − 0(cid:107)L2(D) main advantage of this approach i∇s t∧hat the functional (cid:90) +α ∆v(x)2 dx JBZ is a purely differential functional, while the D\(Γ0∪Γ1)| | functional JMS is an integro-differential one. We +β N−1(Γ D) expect that some results about the Blake-Zisserman 0 +γ HN−1((Γ ∩ Γ ) D), (39) modelthatcouldberephrasedintotomographicterms 1 0 H \ ∩ Tomography: mathematical aspects and applications 10 wouldprovideimmediatelynewresultsintomography. [3] Gel’fand I M and Shilov G E 1966 Generalized Functions: Conversely all the peculiar tomographic features as Properties and Operations vol5,(AcademicPress). [4] HelgasonS1980TheRadonTransform vol5,(Boston,MA: the intrinsic vector nature of the variable E, the fact Birkh¨auser). that its support cannot be bounded and the extra- [5] Strichartz R S 1982 American Mathematical Monthly 89, constraint E =0,motivatenewresearchdirections 377. in the stu∇dy∧of free-discontinuities problems. For [6] Mancini S, Man’ko V I and Tombesi P 1995 Quantum Semiclass. Opt.7,615. example, an interesting result in this context would [7] Man’ko O V and Man’ko V I 1997 J. Russ. Laser Res. 18, be the determination of a good hypothesis on the 407. datum E that ensure that the charge density f do [8] AsoreyM,FacchiP,Man’koVI,MarmoG,PascazioSand g not concentrate. SudarshanECG2007Phys. Rev. A76,012117. [9] FacchiP,Ligabo`MandPascazioS2010JournalofModern We conclude this section with some comments: Optics 57,239. We proved that the measured data g are actually [10] Man’ko O V, Man’ko V I and Marmo G 2000 Phys. Scr. • 62,446. the measured electric field produced by the [11] Man’ko O V, Man’ko V I and Marmo G 2002 J. Phys. A: unknown charge density, so the term f Math. Gen.35,699. (cid:107)R − g 2 in the functional is a fidelity term in a [12] Man’ko V I, Marmo G, Simoni A and Ventriglia F 2006 (cid:107)L2(P3) Open Sys. & Information Dyn.13,239. weak sense. [13] Man’ko M A, Man’ko V I and Mendes R V 2001 J. Phys. The problem of the reconstruction of the charge A34,8321. • can be rephrased into a reconstruction problem [14] AsoreyM,FacchiP,Man’koVI,MarmoG,PascazioSand SudarshanECG2008Phys. Rev. A77,042115. for the electric field. The electric field is an [15] WignerEP1932Phys. Rev.40,749. irrotational vector field, so the new minimization [16] MoyalJ1949Proc. Camb. Phil. Soc.45,99. problem is actually a constrained minimization. [17] HillaryM,O’ConnellRF,ScullyMOandWignerE1984 In order to avoid this constraint one could Phys. Rep.106,121. [18] BertrandJandBertrandP1987Found. Phys.17,397. reformulate the reconstruction problem in terms [19] VogelKandRiskenH1989Phys. Rev. A40,2847. of the electric potential ϕ (E = ϕ) obtaining [20] FacchiPandLigab`oM2010AIP Conf. Proc.1260,3. −∇ a third-order functional in which the fidelity term [21] Deans S R 2007 The Radon transform and some of its is applications,(NewYork: DoverPublications). [22] Paris M G A and Rˇeh´aˇcek J 2004 Quantum State ϕ ϕ 2 , (41) Estimation, Lecture Notes in Physics vol. 649, (Berlin: (cid:107)∇ −∇ g(cid:107)L2(R3) Springer). where the potentials are given by (36). [23] Asorey M, Facchi P, Florio G, Man’ko V I, Marmo G, Pascazio S, Sudarshan E C G 2011 Phys. Lett. A 375, All this considerations hold true in dimension 3. • 861. In a generic dimension n 3 the situation [24] Smith K T, Solmon D C and Wagner S L 1977 Bull. Am. ≥ is quite different because the inversion formula Math. Soc.83,1227. for the Radon transform involves a (possibly [25] Mumford D and Shah J 1989 Commun. Pure Appl. Math. 42,577. fractional) power of the Laplacian. In this case [26] RamlauRandRingW2010InverseProblems 26,115001. the electrostatic description of tomography given [27] RamlauRandRingW2007J. Comput. Phys.221,539. in this section fails. In order to restore it, [28] RondiLandSantosaF2001ESAIMControlOptim.Calc. it is necessary to consider another Radon-type Var.6,517. [29] Blake A and Zisserman A 1987 Visual reconstruction, transformwhichinvolvesintegralsoff overlinear (Cambridge: TheMITPress). manifolds with codimension d such that (n [30] Natterer F 1986 The mathematics of computerized tomog- − d)/2=1, i.e. d=n 2, see e.g. [4, 20]. raphy vol32,(SIAMClassicinAppliedMathematics). − [31] Morel J M and Solimini S 1995 Variational methods in image segmentation vol14,(Boston: Birkh¨auser). Acknowledgments [32] Yuille A 2012 http://en.wikipedia.org/wiki/File:Mumford shah-example.png.(Theoriginalimagehasbeenpartially WethankG.Devillanova,G.FlorioandF.Maddalena modifiedhere.) [33] De Giorgi E, Carriero M and Leaci A 1989 Arch. Ration. for for helpful discussions. Mech. Anal.108,195. This work was supported by “Fondazione Cassa [34] DeGiorgiEandAmbrosioL1992Atti Accad. Naz. Lincei di Risparmio di Puglia” and by the Italian National 8 168,89. Group of Mathematical Physics (GNFM-INdAM). [35] Dal Maso G, Morel J M and Solimini S 1988 Acta Matematica 82,199. [36] Ambrosio L, Fusco N and Pallara D 2000 Functions of References Bounded Variation and Free Discontinuity Problems, (Oxford: OxfordUniversityPress). [1] RadonJ1917Bre. Sachsische Akad. Wiss.69262. [37] Lops F A, Maddalena F and Solimini S 2001 Annales de [2] John F 1955 Plane waves and spherical means: Applied l’institut Henri Poincar´e (C) Analyse non lin´eaire 18, to Partial Differential Equations, (New York: Wiley 639. Interscience). [38] BonnetAandDavidG2001Ast´erique 274,pp.vi+259. [39] FuscoN2003Milan J. Math.71,95.