Opuscula Math. 39, no. 4 (2019), 453–482 https://doi.org/10.7494/OpMath.2019.39.4.453 Opuscula Mathematica APPLICATIONS OF PDES INPAINTING TO MAGNETIC PARTICLE IMAGING AND CORNEAL TOPOGRAPHY Andrea Andrisani, Rosa Maria Mininni, Francesca Mazzia, Giuseppina Settanni, Alessandro Iurino, Sabina Tangaro, Andrea Tateo, and Roberto Bellotti Communicated by Zdzisław Jackiewicz Abstract. In this work we propose a novel application of Partial Differential Equations (PDEs) inpainting techniques to two medical contexts. The first one concerning recovering of concentration maps for superparamagnetic nanoparticles, used as tracers in the framework of Magnetic Particle Imaging. The analysis is carried out by two set of simulations, with and without adding a source of noise, to show that the inpainted images preserve the main properties of the original ones. The second medical application is related to recovering data of corneal elevation maps in ophthalmology. A new procedure consisting in applying the PDEs inpainting techniques to the radial curvature image is proposed. The images of the anterior corneal surface are properly recovered to obtain an approximation error of the required precision. We compare inpainting methods based on second, third and fourth-order PDEs with standard approximation and interpolation techniques. Keywords: PDEs inpainting, medical imaging, magnetic particle imaging, radial curvature image, anterior surface of a cornea. Mathematics Subject Classification: 35Q68, 68U10, 82D80, 92C55, 94A08. 1. INTRODUCTION Medical and biomedical experts often draw their conclusions by the analysis of tomographic or 2D-3D images acquired by various techniques and methodologies. Sometimes, these images may require reconstruction of empty areas, as they were missing data for different reasons. In the framework of image restoration, this task is often achieved by means of some mathematical techniques and algorithms, globally collected under the name of inpainting. Inpainting techniques are mainly used in restoration of old paintings by museum artists, in recovering scratches in a camera picture, in altering scenes c Wydawnictwa AGH, Krakow 2019 453 454 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. in photographs, or in restoration of motion pictures (the reader can refer to [5, 13] and references therein). In this work we consider PDEs-based digital inpainting methods which are usually applied on images with large-scale missing domains, mimicking the basic techniques of professional art restorators [3, 50]. They were first introduced in the pioneering works of Bertalmio et al. [2, 3], where damaged areas of an image are restored by means of transport equations, which smoothly propagates broken isophote lines inside the inpainting region. We review and implement six PDEs inpainting methods: three of second order, one of third order and two of fourth order. These methods are compared to two numerical techniques available in the literature: the automated smoothing procedure, based on a penalized least squares method, implemented in the function inpaintn [22]; the spline interpolation algorithm specified in the MATLAB R griddata function. We applied and tested the six PDEs inpainting techniques in the following two medical contexts: 1. Magnetic Particle Imaging (MPI), to recover concentration maps of superparamagnetic nano-particles. MPI is a tracer-based method for recovering medical images. The spatial distribution c(x) of a tracer material injected inside the volume under investigation is computed applying external magnetic fields to the field of view, and then numerically processing and analyzing the signal s(t) induced by the material. Similarly to other imaging methods such as PET or SPECT, MPI assumes superparamagnetic nanoparticles as the tracer material. Consequently, the MPI reconstruction problem consists in deriving the nanoparticles concentration c from the induced signal s. The pioneering work in MPI dates back to 2005, when Gleich and Weizenecker introduced this new medical imaging technique in [24]. We will show that inpainting based on PDEs is a useful technique to solve in a novel way the MPI reconstruction problem. The nanoparticles distribution c(x) applies only to a subset of image points. Concentrations in the remaining image points have to be determined by some interpolation or approximation algorithms. Then the inpainting approach can be used as an approximant to reconstruct final MPI images. In this work we test the PDEs inpainting techniques on standard “phantom” objects presenting different geometric features, which are typically used for testing MPI devices in preclinical investigations (see for example [24, 32]). The analysis carried out by two set of simulations in the MPI framework, with and without adding a source of noise, shows that the inpainted images preserve the main properties of the original ones. 2. Ophthalmology, to reconstruct elevation maps of the anterior surface of human corneas from real data given in output by the corneal topographer PrecisioHDTM of the iVis Technologies [30]. In the treatment of refractive problems, 3D mapping of the anterior surface of a human cornea, generated by a corneal topographer, is specifically designed for diagnosis and custom vision surgery. Sometimes it occurs that some areas of the detected domain are missed. Usually, these areas occur near the boundary of the scanned region, because of acquisition problems related to the laser technology or of random measuring errors when, for instance, the patient closes his eyes or moves his pupils. This motivates our idea herien proposed of Applications of PDEs inpainting to magnetic particle imaging. . . 455 applying inpainting techniques based on PDEs for recovering missed data and filling the empty areas of the scanned domain. The inpainting techniques applied to the elevation data did not give satisfactory results as expected. Thus we propose a new approach based on radial curvature that improves the results on the reconstruction of the elevation maps. Some recent applications of inpainting techniques to medical images recovering and reconstruction can be found in the literature [7]. Nevertheless, the use of PDEs digital inpainting in the medical frameworks mentioned above constitutes a novelty. The article is organized as follow. After a brief review of PDEs inpainting techniques in Section 2, their application to the medical problems mentioned above is discussed in Sections 3 and 4. Some tests and simulations to validate the use of these PDEs inpainting techniques have been also performed and compared to other approximation methods. Finally, Section 5 concludes. The achieved results are of comparable accuracy with those derived by standard and well tested numerical techniques usually applied in medical field, even better in the presence of noise or when, in the MPI reconstruction problem, Brownian relaxation for the superparamgnetic material is considered. 2. BRIEF REVIEW OF INPAINTING MODELS BASED ON PDES As observed in the Introduction, the inpainting technique has been proposed and numerically implemented in several ways. With regard to digital inpainting methods based on PDEs, in the pioneering work of Bertalmio et al. [3] the PDE governing the algorithm was given by the following transport equation: ∇T u · ∇∆u = 0, (2.1) where the gradient ∇T u = (−∂y u, ∂x u) indicates the isophote directions of an image u. The Laplacian ∆u takes the role of a smoothness estimator. This idea was further developed in [2], where the smooth propagation along the isophote lines was mathematically reformulated by means of the PDE analogous to the 2D stationary Euler equation of the classical fluid dynamics. One of the novelties in Bertalmio’s et al. approach was the use of higher-order PDEs for inpainting, against previous works concerning digital image processing that mainly treat second-order PDEs [1, 14, 37, 47, 48]. Caselles et al. [14] singled out a set of elliptic equations to the approximation problem from data on curves and/or points, following the axiomatic approach previously proposed in [1] for the multiscale image processing. The authors recovered well known PDEs, such as the Laplace’s and the p-Laplace’s equation, and anisotropic diffusions such as Total Variation (TV) and, in particular, Absolute Minimal Lipschitz Extension (AMLE) equation, which admits a unique viscous solution with bounded gradient in most applicative problems. Shen and Chan [53] considered constrained TV equation for their variational approach to disocclusion, derived by the denoising model of Rudin and Osher [48] and Rudin et al. [47], with solutions in the space of bounded variation (BV) functions. (Disocclusion is 456 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. the ability to reconstruct partially occluded objects in a visual perception; it is often associated with inpainting, aiming at similar tasks). Generally anisotropic diffusions as TV equation present some advantage with respect to the Laplace’s equation, being their solutions able to preserve edges in disocclusion and denoising/edge-detection problems as well [46,48,53]. Nevertheless, as a second-order model, TV equation suffers of artificial corners – i.e. rump effect – and limited connectivity between distant edges in the reconstructed image. In this respect TV inpainting model significantly deviates from typical human attitudes in visual perception, as pointed out by the Gestalt school of psychology (Kanizsa [31]). To overcome these drawbacks, Shen and Chan [52] and Chan et al. [15] also considered higher-order inpainting/disocclusion PDEs models: the curvature driven diffusion (CDD) and the Euler’s elastica models, respectively. Other contributions to inpainting arise from segmentation models (see [43] for an overview of such topics), aiming at a segmented reconstruction of inpainting domains. They are based on the minimization of some functionals F . Denoted by Ω ⊂ R2 the image domain, D ⊂ Ω the inpainting sub-domain and f the known image on Ω \ D, the general expression of F acting on an “object function” u and a set of edges Γ, is the following Z   F (u, Γ) = µ r u(x), Du(x), . . . , D(m) u(x) dx + α |Γ| Ω\Γ + λ0 Z (2.2) |u(x) − f (x)|2 dx, Ω\D where D(·) denotes the order of derivative of u. The parameter λ0 is a constant usually called fidelity parameter, µ and α are constants, |Γ| is the one dimension Hausdorff measure of the set Γ, while r(·) is a local function of the derivatives of u (u usually lies 2 in the space of BV functions on Ω). For example, r(·) can be defined as |∇u| [21, 53] in the inpainting model derived by the Mumford-Shah segmentation model [44], or 2 |∆u| [13] in the inpainting model based on the Blake-Zisserman image model [6]. Notice that the main difference occurring from segmentation models to inpainting models lies in the last integral on the right-hand side of (2.2) – the fidelity term – which is now defined onto the complementary set of the inpainting sub-domain D instead of the whole image domain Ω. To overcome the limitations related to low-order models, Esedoglu and Shen [21] suggested to add a curvature term to |Γ| on the right-hand side of (2.2). The inpainting Mumford-Shah-Euler (MSE) image model reads as Z Z
2 α + βk 2 (s) ds FM SE (u, Γ) = µ |∇u| (x)dx + Γ Ω\Γ + λ0 Z |u(x) − f (x)|2 dx, Ω\D where the parameter β is a constant and k(s) is the curvature function. Applications of PDEs inpainting to magnetic particle imaging. . . 457 The gradient descent for a MSE approximating functional – resulting from a De Giorgi conjecture [19] about the Γ-convergence of elliptic functionals – reads [21]
∂t u = ∇ · z 2 ∇u + λ0 χ(Ω \ D) (f − u) , ′   ′′ β W (z)  2 ∂t z = − |∇u| z − α + 2 W (z) 2ε∆z − 2ε 4ε ′  W (z)  + 4β∆ 2ε∆z − . 4ε (2.3) (2.4) The boundary conditions are given by ∂u/∂ν = ∂z/∂ν = ∂g/∂ν = 0, with g = 2ε∆z − W ′ (z)/(4ε), and ν being the normal vector of the frontier, say ∂Ω, of Ω. Here χ(Ω \ D) denotes the indicator function on the set Ω \ D, and z is an “edge signature function”, i.e. a function that nearly equals 1, except for a small tubular neighborhood of Γ of size ε where it is close to zero. Finally, W is a suitable function of z. Note that if we get rid of the first two terms on the right-hand side of (2.4), we get the Cahn–Hilliard equation. The idea of Bertozzi et al. [5] was to propose a shape inpainting model (the BEGCH model) – i.e. working only for binary images – based just on the Cahn–Hilliard equation term plus the fidelity term, that is    ′ ∂ u = −∆ ε∆u − W (u) + λ χ(Ω \ D) (f − u) in Ω 0 t ε (2.5) ∇u · ν = ∇ (∆u) · ν = 0 on ∂Ω. The function W is usually given by the double-well potential, namely W (u) = (u − 1)2 (u + 1)2 or W (u) = u2 (u − 1)2 . Nevertheless, a combination of logarithmic terms instead of polynomials can define a different function W, as proposed in [16]. The BEGCH model presents some advantages over the Γ-convergence approximation of the MSE inpainting model. In particular, it results to be up to an order of O(N log(N )), where N is the number of pixel interested by the inpainting process, regardless of the complexity of the image [5]. Some generalizations of the BEGCH model for grey-values or color images, keeping the simplicity and computational cost of the original one, can be found in [7, 12, 17]. A further generalization of the classical Cahn–Hilliard equation with dynamic boundary conditions has been proposed in [40], which could be useful for restoration of motion pictures. In [12] Burger et al. derived their TV-H−1 inpainting equation model by considering a total variation functional with finite values when applied to functions |u| ≤ 1 (observe that the energy functional associated to the Cahn–Hilliard equation Γ-converges, in L1 topology as ε → 0, to a total variation functional with finite values limited to indicator functions u = 0 or 1 [41,42]). A different approach was followed by Bosch and Stoll [7] for grey-values images and, independently, by Cherfils et al. [17] for color images. They both considered a vector-valued BEGCH equation, with each component describing the distribution of a given grey-value intensity or color of the restored image. With respect to (2.5), an interaction term among the vector components is added to avoid ambiguity in assigning grey-values or colors to the image pixels. 458 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. 2.1. NUMERICAL IMPLEMENTATION OF SOME INPAINTING MODELS In this paper we have applied the following second-, third- and forth-order PDEs based inpainting techniques to the medical contexts described in Sections 3 and 4. a) Harmonic Equation (HE) [50]: ( ∆u − λ0 χ(Ω \ D) (f − u) = 0 in Ω, u=f on ∂Ω. (2.6) This equation has been numerically solved by the gradient descent method. To this end we applied the function heat_equation [49] in the MATLAB R software package. b) Free Noise Harmonic Equation (FHE): ( ∆u = 0 in D, u=f in Ω \ D. (2.7) This equation has been numerically solved by using finite difference method for partial derivatives and the resulting linear system is solved for the unknown elements. To this end we applied the function inpaint_nans (with input parameter 0) [20] in the MATLAB R software package. c) TV equation:   ( ∇u + λ0 χ(Ω \ D) (f − u) = 0 in Ω, ∇ · |∇u| u=f on ∂Ω. (2.8) This equation admits at least a solution in the space of the BV functions in Ω when f is a BV function as well [15]. This equation has been numerically solved by the gradient descent method. To this end we applied the function bvl2_inpainting_convs [49] in the MATLAB R software package. d) CDD equation:   ( ∇ · G(x,|k|)∇u + λ0 χ(Ω \ D) (f − u) = 0 in Ω, |∇u| u=f on ∂Ω, where k denote the curvature of u, while G is defined as follows ( 1 if x ∈ Ω \ D, G(x, s) = g(s) if x ∈ D. (2.9) (2.10) The function g(s) denotes the intensity of the diffusion coefficient in (2.9) driven by the curvature of u. Following Brito and Chen [10], we set g(s) = a+bsp , with model 459 Applications of PDEs inpainting to magnetic particle imaging. . . parameters a << b, such that a numerical solution of (2.9) can be determined by a faster fixed point method with respect to the gradient descent method proposed in [52]. e) Euler’s elastica equation (EE): ( ∇ · V + λ0 χ(Ω \ D)(f − u) = 0 in Ω, ∇u · ν = ∇ (φ′ (k)|∇u|) · ν = 0 on ∂Ω with V = φ(k) (2.11) ∇u ∇T u − ∇ (φ′ (k)|∇u|) · ν, |∇u| |∇u| and φ(k) = a + bk p . A numerical solution of (2.11) can be determined by the algorithm proposed by Tai et al. [54], a faster and efficient method based on a constrained minimization problem related to an augmented Lagrangian functional. Other algorithms to efficiently solve (2.11) can be found in [9, 11]. f) T V − H −1 equation: ∆p + λ0 χ(Ω \ D) (f − u) = 0, where ( |Du (Ω)| , T V (u) = +∞, p ∈ ∂T V (u), (2.12) if |u(x)| ≤ 1 a.e. in Ω, otherwise. (2.13) Elements p ∈ ∂T V (u) are given by [12, Theorem 1.5],[57, Proposition 4.1]:    ∇u  a.e. on {|u| < 1}, −∇ ·     |∇u|  ∇u (2.14) p = −∇ · |∇u| + p̃, p̃ ≤ 0 a.e. on {u = −1},      ∇u  −∇ · |∇u| + p̃, p̃ ≥ 0 a.e. on {u = +1}, provided that ∇u/|∇u| · ν = 0 on ∂Ω plus some other conditions on the Cantor part of the distributional derivative Du of u. Here the convention ∇u/|∇u| = 0 if ∇u = 0 is used. Further, p̃ ∈ ∂χ1 (u) with ( 1 if |u| ≤ 1 a.e. in Ω, χ1 (u) = +∞ otherwise. With a suitable choice of p in (2.14), TV-H−1 equation (2.12) now reads i  ( h ∇u ∆ ∇ · |∇u| − λ0 χ(Ω \ D) (f − u) = 0 in Ω, ∇u |∇u| ·ν =0 on ∂Ω. (2.15) It is worth noting that the first term of (2.12), or (2.15), can be seen as the gradient flow for the functional TV in (2.13) with respect to the H −1 norm. That justifies 460 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. the name to equation (2.12) given by Burger et al. [12]. Furthermore, equation (2.15) differs from the above mentioned TV inpainting model [53] by the Laplacian operator behind the first term on the right-hand side of (2.8). The existence of a solution u to the equation (2.12) in the space of BV functions on Ω is proved in [12, Theorem 1.4], provided f ∈ L2 (Ω). Equation (2.15) has been numerically solved by the gradient descent method. To this end we applied the function bvnegh_inpainting_convs [49] in the MATLAB R software package. g) Free Noise biharmonic equation (FBHE): ( ∆(∆u) = 0 in D, u=f in Ω \ D. (2.16) This equation has been numerically solved by using finite difference method for partial derivatives and the resulting linear system is solved for the unknown elements. It is implemented in the function inpaint_nans (with input parameter 3) [20] in the MATLAB R software package. About the TV, CDD and TV-H−1 models, we considered a regularized version – as done in [12, 52, 53] respectively – by replacing |∇u| in (2.8), (2.9) and (2.15) with p (∇u)2 + δ 2 , being δ a small number known as regularity parameter. To numerically solve the TV (2.8) and TV-H−1 (2.15) PDEs, a semi-implicit convexity splitting method (see [5] and references therein) has provided an unconditionally gradient stable time-discretization scheme [51], allowing larger time steps for the related gradient flow equations compared to the HE (2.6) (see Table 1). Finally, we observe that inpainting PDEs don’t always admit a unique solution, as pointed out by Chan et al. [15] for the TV and EE equations. The numerical solutions obtained by the iterative methods mentioned above may actually depend on the initial guess function u0 . In our tests we have set u0 = 0 and u0 = f respectively inside and outside the inpainting domain D. 3. MAGNETIC PARTICLE IMAGING MPI assumes superparamagnetic nanoparticles as the tracer material. The spatial distribution c(x) of a tracer material injected inside the volume under investigation is computed applying external magnetic fields to the field of view, and then numerically processing and analyzing the signal s(t) induced by the material. The reaction of the magnetic nanoparticles to the external excitation fields can be modeled by a non-linear response, which is filtered, analysed and used for reconstructing the final image. A combination of a static selection field and a time changing magnetic field (drive field) is applied to the tracer material; their superposition saturates all the particles inside the volume except those in a small area, called field free point (FFP), which rapidly moves around the scanned region along paths as the Cartesian trajectories of Figure 1. 461 0.004 0.004 0.002 0.002 y−axis (m) y−axis (m) Applications of PDEs inpainting to magnetic particle imaging. . . 0 −0.002 −0.004 0 −0.002 −0.002 0 0.002 0.004 −0.004 x−axis (m) −0.002 0 0.002 0.004 x−axis (m) Fig. 1. An example of FFP Cartesian trajectories, along the x-(left) and y-(right) direction The induced signal s(t), usually measured by means of some receive coils, just depends on the tracer concentration c around the FFP instantaneous position xF F P . Consequently, the MPI reconstruction problem consists in deriving the nanoparticles concentration c from the induced signal s. Basically, there exist two main lines of research in the literature for solving the MPI reconstruction problem, as described in [27]: the measurement based (MB) approach [24, 34], and the model based approach including x-space (XS) [25, 26]. The MB approach is based on the inverse problem expressing the non-linear relation between the concentration c and the frequency response of the signal. The Fourier transform of the induced signal, sf t , is linearly related to c by the following linear system sf t = S c. (3.1) Here S is a matrix called system function; it depends on the physical configuration of the apparatus and on the nanoparticles properties. By measuring sf t and inverting the system (3.1), the concentration c can be recovered by computing the solution in the least squares sense or with the help of regularization techniques, such as the singular value decomposition (SVD), the Tikhonov’s method or more advantageous iterative methods as the Kaczmarz method, the CGLS and LSQR approach [32, 36]. The XS approach, proposed by Goodwill et al. in [25, 26], differs from MB because the reconstruction problem is not formulated as an inverse problem. It is based instead on the assumption that the nanoparticles distribution c – once convolved with an opportune point spread function – is proportional to the induced signal s generated on a receive coil over the rate of change of the nanoparticles saturation. Set û the axis of the receive coil, it results [26] s(t) ∝ (c ∗ [û · H v̂F F P (t)]) (xF F P (t)) . |vF F P (t)| (3.2) Here vF F P = |vF P P | v̂F F P is the FFP velocity vector and H is a space variable function with values in Rd × Rd . H indicates the area where the paramagnetic nanoparticles are not yet saturated. It is usually modeled by the derivative of the Langevin function (see, for instance, [26, 32]). 462 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. From (3.2) one deduces that the point spread function – the term between the square brackets – actually depends on v̂F F P (t), which varies with time. However, by adopting Cartesian trajectories for the FFP ([32]), v̂F F P (t) is mainly directed along a fixed direction. Thus one can assume v̂F F P constant, so simplifying the deconvolution process. Observe also that the inversion points on the Cartesian trajectories are usually neglected, as they can generate ambiguity in (3.2) being both s and vF F P close to zero. On the one hand, XS is a simpler approach than those methods based on the inversion of the system function. On the other hand, XS works only for uniformly spatial drive fields. In recent years, a different Hybrid XS (HXS) approach, combining both the MB and XS, was also proposed in [55, 56]. An interesting overview of all methods techniques adopted in MPI, together with their preclinical investigations, can be found in [33]. 3.1. APPLICATION OF PDES INPAINTING TO MPI In this section inpainting based on PDEs is showed to be a useful technique to solve in a novel way the MPI reconstruction problem by the application of the XS approach. Formula (3.2), describing the nanoparticles distribution, applies only to those points belonging to the FFP trajectory. Concentrations in the remaining image points have to be determined by some interpolation or approximation algorithms. Then the inpainting approach can be used as an approximant to reconstruct final MPI images by the XS method. In this preliminary study we test the PDEs inpainting techniques mentioned in Section 2.1 on standard “phantom” objects presenting different geometric features, which are typically used for testing MPI devices in preclinical investigations (see for example [24, 32]). We consider two phantom objects, a “bar” and the letter “P”, whose binary images are shown in Figure 2. Both these images are of the same size, say 8 × 8 mm2 , and are sampled at 100 × 100 pixels. The phantom “bar” is of about 6 mm length and 0.7 mm thick, while the object “P” is composed by seven distinct spheres of 0.8 mm length diameter. The distance between two spheres is equal to 1.2 mm, thus close to the theoretical intrinsic resolution of the simulated apparatus as described hereafter. Fig. 2. Images used for the simulation study. Left: Phantom “P”, Right: Phantom “Bar”. Size: 8 × 8 mm2 , resolution: 100 × 100 pixels To carry out MPI simulations, we considered iron-oxide nanoparticles (Resovistr ) with core diameter D = 30 nm, and a gradient for the selection field equal to 5 T/m, Applications of PDEs inpainting to magnetic particle imaging. . . 463 in order to achieve a theoretical resolution of ≈ 0.5 mm at a temperature T = 310 K [25]. Drive fields along x- and y- directions oscillated at 25 kHz and 25/24 kHz to implement Cartesian trajectories for the FFP (Figure 1). The frequency sampling was set to 2.5 MHz. For the sake of simplicity nanoparticles was assumed to have a pure Langevin-type behaviour [32]. To avoid anisotropy in the resolution of 2D reconstructed images by XS (see [26]), the simulated scanning process for each phantom object was performed in both x- and y-directions of the Cartesian trajectories (Step 1). The obtained partially filled images have now to be approximated. In the present paper we performed this approximation by processing the images with the following six different PDEs inpainting procedures: the second-order HE (2.6), FHE (2.7) and TV (2.8) equations, the third-order CDD (2.9), and the forth-order EE (2.11) and TV-H−1 (2.12) equations. The parameter values used in the corresponding numerical algorithms are listed in Table 1 (Step 2). Table 1. Parameter values used in the inpainting PDEs algorithms. For a full description about the meaning of the parameters we refer to [5, 10, 12, 15, 50–53] HE FHE TV CDD EE TV-H−1 Time step τ 0.01 − 100 − − 100 Mesh size hx hy 1 1 1 1 1 1 1 1 1 1 1 1 Fidelity parameter λ0 100 − 100 100 100 100 Regularity parameter δ − − 0.01 0.01 − 0.01 Convexity splitting coefficients C1 C2 − − − − 100 100 − − − − 100 100 Model parameters a b p − − − − − − − − − 0.01 1 2 1 1 2 − − − The resulting images ux and uy was firstly processed by a Wiener’s deconvolution (Step 3) and successively pmerged pixel by pixel into the final image u by a geometrical mean, that is u(i, j) = ux (i, j)uy (i, j), as done in [56] (Step 4). At this last step a further processing by a high-pass filter was performed. Figure 3 illustrates these steps as regard to the HE inpainting method applied to the phantom object “P”. Observe that in our simulation study the XS approach was applied by neglecting any problem related to the “fundamental frequency” loss (the reader can refer to [45] and references therein for a discussion about this issue). Moreover, sources of “noise” were not considered in this first set of simulations. In Figure 4 the performances of the six PDEs inpainting methods mentioned above are compared to two numerical techniques available in the literature: 1) the automated smoothing procedure, based on a penalized least squares method, implemented in the function inpaintn [23] created by Garcia [22]; 2) the spline interpolation algorithm specified in the MATLAB R griddata function (with the v4 option). A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. 1 0.004 0.8 y−axis (m) 0.002 0.6 0 0.4 −0.002 −0.004 0.2 −0.002 0 x−axis (m) 0.002 0.004 0 1 0.004 0.8 0.002 y−axis (m) 464 0.6 0 0.4 −0.002 0.2 −0.004 −0.002 0 0.002 0.004 0 x−axis (m) Fig. 3. HE inpainting procedure applied to the phantom “P” object. Top row: scanned XS images by the Cartesian trajectories, in false colour; the inpainting masks are given by the black points (Step 1). Second row: corresponding inpainting recovering (Step 2). Third row: Wiener deconvolution (Step 3). Bottom row: True phantom object (left) and final merged image u (right), further processed by a high-pass filter (Step 4) 465 Applications of PDEs inpainting to magnetic particle imaging. . . HE FHE TV CDD r=0.857 PSNR=16.75 SSIM=0.504 CPU=2.665 r=0.868 PSNR=17.08 SSIM=0.513 CPU=0.262 r=0.812 PSNR=15.89 SSIM=0.490 CPU=11.41 r=0.827 PSNR=16.08 SSIM=0.492 CPU=0.409 inpaintn griddata (v4) EE −1 TV−H r=0.835 PSNR=16.37 SSIM=0.493 CPU=77.06 r=0.857 PSNR=16.41 SSIM=0.513 CPU=6.538 r=0.868 PSNR=17.01 SSIM=0.508 CPU=2.952 r=0.870 PSNR=16.87 SSIM=0.520 CPU=2.293 HE FHE TV CDD r=0.850 PSNR=18.39 SSIM=0.572 CPU=2.101 r=0.869 PSNR=18.80 SSIM=0.597 CPU=0.081 r=0.790 PSNR=17.34 SSIM=0.492 CPU=10.33 r=0.838 PSNR=18.24 SSIM=0.518 CPU=0.316 inpaintn griddata (v4) r=0.871 PSNR=18.83 SSIM=0.597 CPU=2.952 r=0.877 PSNR=19.08 SSIM=0.597 CPU=3.673 EE r=0.848 PSNR=18.60 SSIM=0.561 CPU=77.10 −1 TV−H r=0.854 PSNR=18.75 SSIM=0.576 CPU=2.398 Fig. 4. Comparison among the true phantom objects and the final reconstructed images by eight different inpainting methods. For each image the correlation coefficient r, the peak signal-to-noise ratio P SN R, the structural similarity SSIM and the CP U are reported 466 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. By analogy with the analysis performed in [7], for each method we computed three parameters measuring the quality of the reconstruction with respect to the original image: the correlation coefficient r, the peak signal-to-noise ratio PSNR (the higher the better, see for example [7]), the mean structural similarity SSIM introduced by Wang et al. [58] (values are between −1 and 1, those closer to 1 indicate a better quality of the reconstructed image). The computational time CPU in seconds was also calculated. All these parameters are reported in Figure 4 for each reconstructed image. Note that in order to evaluate the CPU parameter, the following stopping criterion was applied to the six PDEs iterative methods: kun − un−1 k1 < η, kun−1 k1 (3.3) where k·k1 denotes the L1 norm and un the n-th iterative solution. The threshold value η was set equal to 0.5 · 10−5 . Note that the criterion (3.3) could not be applied to inpaintn, because this threshold value was never reached in our tests. Thus we stop the iterative procedure after 1000 iterations. All computations have been executed in MATLAB R R2010a on a 64-bit computer with a single CPU type Intel R CoreTM i3-4005U at 1.70 GHz, dual core, and with 12 GB main memory available. Figure 4 shows that the FHE, the inpaintn and the griddata methods provide the best recovered MPI images for each phantom, even if with a faster computational time for FHE. The high computational time for EE is essentially due to the complexity of the model, if compared to the unconditionally stable TV-H−1 . Further tests involving the AMLE model, the original transport equation (2.1), the vector-valued BEGCH, and its modification by using a logarithmic term as a double well potential, did not improve the results showed in Figure 4. Finally, we wish to point out some issues relating the different results obtained by the second-order HE, the third-order CDD and the forth-order TV-H−1 method. Consider the recovered MPI image ux for the object “P”. The image ux reconstructed by HE, compared with the corresponding ones obtained by CDD and TV-H−1 , highlights a (more evident) trace of the FFP trajectories (see Figure 5). HE CDD TV−H−1 THEOR. Fig. 5. Phantom object “P”: a zoomed section of the recovered MPI image ux inpainted by the HE, CDD and TV-H−1 procedures, compared with the theoretical image obtained by eq. (3.2) 467 Applications of PDEs inpainting to magnetic particle imaging. . . HE FHE TV CDD r=0.801 PSNR=16.15 SSIM=0.385 CPU=2.645 r=0.780 PSNR=15.86 SSIM=0.371 CPU=0.261 r=0.766 PSNR=15.29 SSIM=0.389 CPU=15.89 r=0.753 PSNR=15.43 SSIM=0.358 CPU=2.326 inpaintn griddata (v4) r=0.778 PSNR=15.80 SSIM=0.364 CPU=3.013 r=0.764 PSNR=15.63 SSIM=0.358 CPU=3.579 EE r=0.801 PSNR=16.10 SSIM=0.378 CPU=81.30 −1 TV−H r=0.801 PSNR=16.06 SSIM=0.383 CPU=6.856 HE FHE TV CDD r=0.797 PSNR=17.48 SSIM=0.421 CPU=2.585 r=0.754 PSNR=16.82 SSIM=0.413 CPU=0.085 r=0.725 PSNR=16.44 SSIM=0.413 CPU=6.202 r=0.760 PSNR=16.72 SSIM=0.380 CPU=1.657 EE TV−H−1 inpaintn griddata (v4) r=0.816 PSNR=17.80 SSIM=0.465 CPU=79.45 r=0.820 PSNR=17.97 SSIM=0.470 CPU=7.471 r=0.781 PSNR=17.20 SSIM=0.398 CPU=3.042 r=0.784 PSNR=17.26 SSIM=0.357 CPU=3.606 Fig. 6. Comparison among the final reconstructed images by eight different inpainting methods for the second set of simulations. For each image the correlation coefficient r, the peak signal-to-noise ratio P SN R, the structural similarity SSIM and the CP U are reported 468 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. On the one hand, this indicates that the second-order method is less powerful in connecting distant points of the image. On the other hand, a similar or better performance of the HE method with respect to CDD and TV-H−1 can be observed in terms of the final image quality. Indeed, a further comparison with the theoretical image obtained by (3.2) shows that CDD performs an inhomogeneous completion, while disconnected areas of the original image are “excessively glued” together when the TV-H−1 method is applied. In a second set of simulations (Fig. 6) we considered a more realistic situation by assuming that: a) superparamagnetic nanoparticles have lognormal size distribution with diameter d = 30 ± 3.2 nm; b) receive coils signals are corrupted by a source of white noise with signal-to-noise ratio (SNR) equal to 20 (in line with what reported in [24], Figure 4); c) superparamagnetic nanoparticles are subjected to Brownian relaxation ([32,59] and references therein), with relaxation time equal to 5 · 10−7 s. Relaxation introduces a delay on scans as well as a significant asymmetric blurring of the image in the scanning direction ([18]). Besides it further reduces the SNR. In contrast to the results reported in Figure 4, Figure 6 shows that TV-H−1 and HE methods provide the best recovered MPI images for each phantom object. Also in this case the higher computational time of the EE method penalizes its good performance. This result is expected because, contrary to FHE, inpaintn and griddata, these inpainting methods work also on the original data to denoise the image. 4. APPLICATION OF PDES INPAINTING IN CORNEAL TOPOGRAPHY In this Section we present the use of PDEs inpainting in the reconstructed elevation map of the cornea anterior surface. We do not consider the problem of data acquisition, but we use an approximation of real data of the anterior surface of the cornea. This approximated data are generated on a grid with 360 × 101 uniform points in polar coordinates, namely (θ, ρ) ∈ [−π, π] × [0, 5], which corresponds to the reconstructed elevation map of the cornea anterior surface from data given in output by the corneal topographer PrecisioHDTM of the iVis Technologies [29, 30]. The measuring scale of a corneal topographer is about the order of millimeters. Hence, the efficiency of an approximating algorithm applied to scanned data is measured by a reconstruction error within the order of 1 µm. It is also important to reconstruct the radial derivatives and the radial curvature of the original surface with enough precision, since the curvature of a cornea is responsible of aberrations that can impact the imaging performance of the eye. Specifically, about the 80% of the dioptric power resides in the cornea, a region of the human eye in hemispherical shape of maximum diameter 11-12 mm. Sometimes it occurs that some areas of the detected domain are missed. Usually, these areas occur near the boundary of the scanned region, as shown in Figure 7, which plots the elevation map of a real human eye. That can be caused by some Applications of PDEs inpainting to magnetic particle imaging. . . 469 acquisition problems due to the laser technology, or by random measuring errors when, for instance, the patient closes his eyes or moves his pupils. It is worth noting that, even if the surface seems to be regular, its radial curvature does not present the same regularity, as shown in Figure 8. Fig. 7. Real Eye Elevation Map with Missing Data Left: Contour plot. Right: 3D plot Fig. 8. Real Eye Radial Curvature Map with Missing Data. Left: Contour plot. Right: 3D plot The main goal is then to reconstruct the original anterior surface of the cornea by recovering the elevation and the radial curvature with an accuracy up to 1 µm and 10 µm, respectively. What mentioned above motivates our idea proposed in this work of applying inpainting techniques based on PDEs for recovering missed data and filling the empty areas of the scanned domain. To this end, we first considered the elevation maps of cornea anterior surfaces, acquired by a corneal topographer, as images to be inpainted. As in the classical inpainting setting, areas with missing data are called inpainting region. From now on, we call this approach “elevation inpainting”, because the elevation maps of the corneal anterior surface are considered as grey-values of the image to be inpainted. Since the quality of the reconstructed surface by the elevation inpainting was not satisfactory, we implemented a new procedure, which we call “radialinpainting”. The procedure consists first in applying the PDEs inpainting 470 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. techniques to the radial curvature image and then in properly recovering the image of the anterior corneal surfaces. This approach considerably reduces the approximation error. We applied the following PDEs inpainting techniques, described in Section 2.1, for both the elevation and radial inpainting procedures: the second-order HE (2.6), FHE (2.7), TV (2.8) methods; the third-order CDD (2.9) method; the forth-order TV-H−1 (2.12) and FBHE (2.16) methods. Further, these inpainting techniques have been compared with the automated smoothing procedure implemented in the function inpaintn [22, 23]. The corresponding parameter values are listed in Table 1. 4.1. ELEVATION INPAINTING We started by considering the elevation maps in polar coordinates, (θ, ρ) ∈ [−π, π] × × [0, 5]. We applied digital inpainting techniques to these images in order to reconstruct the empty areas with missing data. Observe that the initial value of the angle θ is set equal to the first value, say θ0 ≥ −π, corresponding to a complete information in the image. To evaluate the performance of the inpainting method we computed the radial curvature values of the inpainted surface. As a first numerical test, we considered a real human eye surface without missing data, and we applied an inpainting mask to the original surface to artificially remove some data points. Thus we were able to compute the approximation error with respect to the original surface. Figure 9 compares the original surface and the corresponding surface with missing data regions to be inpainted, with the surfaces reconstructed by some of the inpainting algorithms mentioned above. The picture highlights the failure of the elevation inpainting when the HE is applied. A similar failure results from the application of the other inpainting methods not included in the picture (TV,TV-H−1 , inpaintn). Better performances are given by the FBE and FBHE algorithms, even if the absolute curvature error is high (see Figure 10). This result has a theoretically explanation: the missing data regions are placed on the border of the surface, but all the applied inpainting methods are suitable for internal inpainting regions, because they need boundary conditions. Only the FHE, FBHE allow incomplete information on the border, so justifying a better behaviour. Figure 10 reports, for all the compared algorithms, the maximum absolute error of the radial curvature and elevation, computed in the region (θ, ρ) ∈ [−π, π] × [0, ρ∗ ] for different values of ρ∗ , ρθ ≤ ρ∗ ≤ 5, where ρθ denotes the higher radius with complete information. This is important in order to recover a region larger than the original one with missing data, but without considering the border where the error is higher. The maximum absolute errors in the radial curvature are very high for all the algorithms, higher than 10 µm when the radius ρ∗ > ρθ . The test results lead us to conclude that this inpainting approach is unable to reconstruct elevations and radial curvatures with the desired precision. 471 Applications of PDEs inpainting to magnetic particle imaging. . . Fig. 9. Example 1. Comparison among the true eye and the final elevation map images reconstructed by the four best performing inpainting methods. Top line (from left to right): original elevation, original elevation with missing data, HE method. Bottom line (from left to right): FHE, FBHE, CDD methods Max Absolute Error in the radial curvature on know data, 10 = 4.45 3 HE TV FHE FBHE TVH1 CDD Inpaintn 10 2 10 1 10 0 10 -1 10 -2 10 -3 4.3 4.4 4.5 4.6 4.7 Max Absolute Error in the elevation on know data, 4.8 4.9 5 = 4.45 10 2 HE TV FHE FBHE TVH1 CDD Inpaintn 10 0 10 -2 10 -4 10 -6 10 -8 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5 Fig. 10. Example 1. Maximum Absolute error (in mm) in the curvature (top line) and in the elevation (bottom line) between the real known area of the eye and the reconstructed area stretching to different radius ρθ ≤ ρ∗ ≤ 5, with ρθ = 4.45 472 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. 4.2. RADIAL CURVATURE INPAINTING As we pointed out in Section 4, the main goal of a corneal reconstruction problem consists in a final approximation with maximum absolute error not larger than 1 µm. At the same time, the radial curvature of the approximating surface should not differ too much from that of the original surface. Unfortunately, as shown in the previous section, the inpainting techniques applied to the elevation data do not give satisfactory results as expected, especially for the radial curvature approximation. Thus we propose a radial curvature based approach improving the results on the reconstruction of the surface as well as of the curvature. Our idea consists in applying digital inpainting techniques to the radial curvature values instead of to the elevation data of the corneal surface. Given a continuous function f (θ, ρ) approximating the elevation data, the radial curvature can be computed, for each angle θ, using the following formula c(θ, ρ) = p f ′′ (θ, ρ) (1 + f ′ (θ, ρ)2 )3 , 0 ≤ ρ ≤ ρθ , (4.1) where ρθ ≤ 5 denotes the higher radius with complete information corresponding to the fixed θ value. Since we know only the elevation data in the grid mesh points, we need to find a tool to approximate the first and second derivatives in (4.1). To this end, we use the object oriented c++ library QIBSH++ that implements Hermite Quasi Interpolant for 1D and 2D functions based on B-spline [4, 8, 28, 38]. Applying inpainting techniques to the curvature allows to overcome the main issue, explained in Section 4.1, when we deal directly with the surface. Indeed, the gradient of the curvature on the boundary is usually very small (for the sphere it is exactly zero) and so the boundary conditions imposed by the applied inpainting algorithms are good approximations of the unknown ones. The procedure here proposed consists in three steps. As a first step, the available elevation data of the surface are used to compute the curvature, along each radial direction. After approximating the partial derivatives of the surface by a generalized backward difference formula (GBDF), as described in [8, 35, 39], we compute a quasi-interpolant Hermite spline and the related curvature. We choose a GBDF approximation of order 4 and a spline quasi-interpolants of degree 3. In a second step, the radial curvature, c(θ, ρ), is recovered using inpainting techniques. Then, the final recovered surface is achieved by radially solving a second-order Ordinary Differential Equation (ODE) to compute the elevation of the real eye starting from the inpainted curvature values. The function f to be recovered is solution of the following differential equation: p f ′′ (θ, ρ) = c(θ, ρ) (1 + f ′ (θ, ρ)2 )3 with fixed initial radius 0 ≪ ρ0 < ρθ , and initial values f (θ, ρ0 ), f ′ (θ, ρ0 ). The second order ODE is solved using the MATLAB R function ODE45. Applications of PDEs inpainting to magnetic particle imaging. . . 473 The last step of the algorithm relies in a final regularization computed by a tensor product Hermite Quasi Interpolant based on B-spline, available in the QIBSH++ library. The quasi-interpolation method, as described in [8, 38], requires function and first derivative values for computing the approximation. Moreover, the choice of a quasi-interpolant is convenient for two main reasons. Firstly, as often it occurs in real applications, interpolating data are affected by errors; hence, to relax the interpolation conditions, smoothing or quasi-interpolation are commonly preferred to interpolation, which may be too stringent. Secondly, the use of a quasi-interpolant has a computational cost less than a global interpolant, this allowing to reduce time and computational effort. The above described numerical approach was applied to the same eye surface described in Section 4.1. As shown in Figure 11, the empty areas have been reconstructed by the curvature inpainting, no oscillation appears in the final inpainted surface. Moreover, by comparing the error results in Figure 13 and in Figure 10, it is evident the better performance of the radial curvature based approach. The original surface and curvature are recovered in the inpainting domains with a maximum absolute error of about 10−3 mm and 10−2 mm, respectively. The best recovered radial curvature in Figure 12 is obtained by the application of FHE and FBHE algorithms. Other real corneal surfaces with different shapes and irregularities, showing larger inpainting domains were also tested. We show the maximum absolute errors in the reconstruction of two additional real eyes illustrated in Figures 14, 16. Fig. 11. Example 1. Radial Curvature Inpainting. Comparison among the true eye and the final elevation map images reconstructed by six different methods. Top line (from left to right): original elevation, original elevation with missing data, HE, FHE methods. Bottom line (from left to right): FBHE, TV-H−1 , CDD, inpaintn methods 474 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. Fig. 12. Example 1. Radial Curvature Inpainting. Comparison among the curvature of the true eye and the curvature of the final elevation map images reconstructed by six different methods. Top line (from left to right): original curvature, original curvature with missing data, HE, FHE methods. Bottom line (from left to right): FBHE, TV-H−1 , CDD, inpaintn methods Max Absolute Error in the radial curvature on know data, 10 = 4.45 -1 HE TV FHE FBHE TVH1 CDD Inpaintn 10 -2 10 -3 10 -4 10 -5 10 -6 4.3 4.4 4.5 4.6 4.7 Max Absolute Error in the elevation on know data, 4.8 4.9 5 = 4.45 10 -2 HE TV FHE FBHE TVH1 CDD Inpaintn 10 -3 10 -4 10 -5 10 -6 10 -7 10 -8 10 -9 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5 Fig. 13. Example 1. Radial Curvature Inpainting. Maximum Absolute error (in mm) in the curvature (top line) and in the elevation (bottom line) between the real known area of the eye and the reconstructed area stretching to different radius ρθ ≤ ρ∗ ≤ 5, with ρθ = 4.45 475 Applications of PDEs inpainting to magnetic particle imaging. . . Fig. 14. Example 2. Top Left: Real Eye Elevation Contour Plot. Top Right: Real Eye Elevation Map with Missing Data. Bottom Left: Real Eye Radial Curvature Contour Plot. Bottom Right: Real Eye Radial Curvature Map with Missing Data Max Absolute Error in the radial curvature on know data, = 2.25 10 -1 10 -2 HE TV FHE FBHE TVH1 CDD Inpaintn 10 -3 10 -4 10 -5 10 -6 1.5 2 2.5 3 3.5 Max Absolute Error in the elevation on know data, 4 4.5 5 4 4.5 5 = 2.25 10 -3 10 -4 10 -5 HE TV FHE FBHE TVH1 CDD Inpaintn 10 -6 10 -7 10 -8 10 -9 1.5 2 2.5 3 3.5 Fig. 15. Example 2. Radial Curvature Inpainting. Maximum Absolute error (in mm) in the curvature (top line) and in the elevation (bottom line) between the real known area of the eye and the reconstructed area stretching to different radius ρθ ≤ ρ∗ ≤ 5, with ρθ = 2.25 476 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. Fig. 16. Example 3. Top Left: Real Eye Elevation Contour plot. Top Right: Real Eye Elevation Map with Missing Data. Bottom Left: Real Eye Radial Curvature Contour plot. Bottom Right: Real Eye Radial Curvature Map with Missing Data Max Absolute Error in the radial curvature on know data, = 2.6 10 -1 HE TV FHE FBHE TVH1 CDD Inpaintn 10 -2 10 -3 10 -4 10 -5 10 -6 2 2.5 3 3.5 Max Absolute Error in the elevation on know data, 4 4.5 5 4 4.5 5 = 2.6 10 -3 10 -4 10 -5 10 -6 HE TV FHE FBHE TVH1 CDD Inpaintn 10 -7 10 -8 2 2.5 3 3.5 Fig. 17. Example 3. Radial Curvature Inpainting. Maximum Absolute error (in mm) in the curvature (top line) and in the elevation (bottom line) between the real known area of the eye and the reconstructed area stretching to different radius ρθ ≤ ρ∗ ≤ 5, with ρθ = 2.6 Applications of PDEs inpainting to magnetic particle imaging. . . 477 In both cases illustrated in Figures 14 and 16 we used an inpainting region larger then the original one, so that to have some information about the elevation data inside the inpainting domain. Thus we computed the elevation and radial curvature errors in this region. The results listed in Figures 15 and 17, highlight, once again, the good performance of the curvature based approach, which, in this case, is able to recover the eye surface with a maximal error no larger than 0.1 µm. 5. CONCLUSION In this work we proposed a novel application of inpainting techniques based on PDEs to two different medical contexts: the first one concerning visual reconstruction of superparamagnetic nanoparticles distributions in MPI clinical images; the second one is related to recovering data of corneal elevation maps in ophthalmology. The analysis carried out by two set of simulations in the MPI framework, with and without adding a source of noise, has shown that the inpainted images preserve the main properties of the original ones. Further, the achieved results are of comparable accuracy with those derived by standard and well tested numerical techniques usually applied in medical field, even better in the presence of noise and when Brownian relaxation is considered. As for the ophthalmological framework, we have implemented a new procedure consisting in applying first the PDEs inpainting techniques to the radial curvature image. Then, the images of the anterior corneal surfaces are properly recovered to improve their quality by reducing the approximation error to the order of precision required. However, we believe that this work is a first attempt to apply PDEs inpainting techniques to the medical contexts described above and that there is still much to improve. Indeed, further developments of this study will investigate on searching the best inpainting model to improve the overall quality of both the intrinsic MPI images and the corneal surface reconstruction, e.g. with a lower gridding definition. In particular, for MPI a lower frequency sampling to minimize the storage data could be used, while for corneal reconstruction the application of PDEs inpainting from scattered volumetric data given in output by the corneal topographer could be also carried out. Acknowledgements The authors thankfully acknowledge the financial support by the National Research Project AMIDERHA “Minimally invasive advanced systems for diagnostics and radiotherapy”, leading proponent MEDIS, PON R&C 2007-2013, PON02-00576-3329762. REFERENCES [1] L. Alvarez, F. Guichard, P.L. Lions, J.M. Morel, Axioms and fundamental equations of image processing, Arch. Ration. Mech. An. 123 (1993), 199–257. [2] M. Bertalmio, A.L. Bertozzi, G. Sapiro, Navier–Stokes, fluid dynamics, and image and video inpainting, Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, vol. 1, 2001, 1–355. 478 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. [3] M. Bertalmio, G. Sapiro, V. Caselles, C. Ballester, Image Inpainting, Proceedings of the 27th annual conference on Computer graphics and interactive techniques, ACM Press/Addison-Wesley Publishing Co. 2000, 417–424. [4] E. Bertolazzi, F. Mazzia, The object oriented C++ library QIBSH++ for Hermite quasi interpolation, 2019. [5] A. Bertozzi, S. Esedoglu, A. Gillette, Analysis of a two-scale Cahn-Hilliard model for binary image inpainting, Multiscale Model. Sim. 6 (2007), 913–936. [6] A. Blake, A. Zisserman, Visual reconstruction, Artificial Intelligence Series, The MIT Press, 2003. [7] J. Bosch, M. Stoll, A fractional inpainting model based on the vector-valued Cahn–Hilliard equation, SIAM J. Imaging Sci. 8 (2015), 2352–2382. [8] C. Bracco, C. Giannelli, F. Mazzia, A. Sestini, Bivariate hierarchical Hermite spline quasi-interpolation, BIT Num. Math. 56 (2016), 1165–1188. [9] K. Bredies, T. Pock, B. Wirth, A convex, lower semicontinuous approximation of Euler’s elastica energy, SIAM J. Math. Anal. 47 (2015), 566–613. [10] C. Brito-Loeza, K. Chen, Multigrid method for a modified curvature driven diffusion model for image inpainting, J. Comput. Math. 26 (2008), 856–875. [11] C. Brito-Loeza, K. Chen, Fast numerical algorithms for Euler’s elastica inpainting model, International Journal of Modern Mathematics 5 (2010), 157–182. [12] M. Burger, L. He, C.B. Schönlieb, Cahn-Hilliard inpainting and a generalization for grayvalue images, SIAM J. Imaging Sci. 2 (2009), 1129–1167. [13] M. Carriero, A. Leaci, F. Tomarelli, Free gradient discontinuity and image inpainting, J. Math. Sci. 181 (2012), 805–819. [14] V. Caselles, J.M. Morel, C. Sbert, An axiomatic approach to image interpolation, IEEE T. Image Process. 7 (1998), 376–386. [15] T.F. Chan, S.H. Kang, J. Shen, Euler’s elastica and curvature-based inpainting, SIAM J. Appl. Math. 63 (2006), 564–592. [16] L. Cherfils, H. Fakih, A. Miranville, On the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation with logarithmic nonlinear terms, SIAM J. Imaging Sci. 8 (2015), 1123–1140. [17] L. Cherfils, H. Fakih, A. Miranville, A Cahn–Hilliard system with a fidelity term for color image inpainting, J. Math. Imaging Vis. 54 (2016), 117–131. [18] L.R. Croft, P. Goodwill, M. Ferguson, K. Krishnan, S. Conolly, Relaxation in x-space magnetic particle imaging, T.M. Buzug, J. Borgert (eds.), Magnetic Particle Imaging, vol. SPPHY 140, Springer-Verlag, Berlin Heidelberg, 2012, pp. 149–153. [19] E. De Giorgi, Area-minimizing oriented boudaries, Sem. di Mat. de Scuola Norm. Sup. Pisa, 1960–1961, 1–56 [in Italian]; Engl. transl. in: L. Ambrosio, G. Dal Maso, M. Forti, M. Miranda, S. Spagnolo (eds.), Ennio De Giorgi: Selected Papers, Springer, Berlin, 2006, pp. 231–263. [20] J. D’Errico, Interpolates (& extrapolates) nan elements in a 2d array, https://www. mathworks.com/matlabcentral/fileexchange/4551, 2012 (accessed January 23, 2017). Applications of PDEs inpainting to magnetic particle imaging. . . 479 [21] S. Esedoglu, J. Shen, Digital inpainting based on the Mumford–Shah–Euler image model, Eur. J. Appl. Math. 13 (2002), 353–370. [22] D. Garcia, Robust smoothing of gridded data in one and higher dimensions with missing values, Comput. Stat. Data An. 54 (2010), 1167–1178. [23] D. Garcia, Inpaint over missing data in 1-D, 2-D, 3-D,... ND arrays, http://it. mathworks.com/matlabcentral/fileexchange/27994-inpaint-over-missing-datain-1-d--2-d--3-d--n-d-arrays, 2017 (accessed November 20, 2017). [24] B. Gleich, J. Weizenecker, Tomographic imaging using the nonlinear response of magnetic particles, Nature 435 (2005), 1214–1217. [25] P.W. Goodwill, S.M. Conolly, The X-space formulation of the magnetic particle imaging process: 1-D signal, resolution, bandwidth, SNR, SAR, and magnetostimulation, IEEE T. Med. Imaging 29 (2010), 1851–1859. [26] P.W. Goodwill, S.M. Conolly, Multidimensional X-space magnetic particle imaging, IEEE T. Med. Imaging 30 (2011), 1581–1590. [27] M. Grüttner, T. Knopp, J. Franke, M. Heidenreich, J. Rahmer, A. Halkola, C. Kaethner, J. Borgert, T.M. Buzug, On the formulation of the image reconstruction problem in magnetic particle imaging, Biomed. Tech. 58 (2013), 583–591. [28] A. Iurino, BS hermite quasi-interpolation methods for curves and surfaces, Ph.D. thesis, Dipartimento di Matematica, Università degli Studi di Bari Aldo Moro, Bari, May 2014. [29] IVIS Technologies, Corneal tomography and topography, http://www.ivistechnologies. it/presentations.php, 2017 (accessed June 5, 2017). [30] IVIS Technologies, Ivis technologies, http://www.ivistechnologies.it, 2017 (accessed June 5, 2017). [31] G. Kanizsa, Organization in Vision, Praeger Pub Text, 1979. [32] T. Knopp, T.M. Buzug, Magnetic particle imaging: an introduction to imaging principles and scanner instrumentation, Springer Science & Business Media, 2012. [33] T. Knopp, N. Gdaniec, M. Möddel, Magnetic particle imaging: From proof of principle to preclinical applications, Phys. Med. Biol. 62 (2017), R124–R178. [34] T. Knopp, A. Weber, Sparse reconstruction of the magnetic particle imaging system matrix, IEEE T. Med. Imaging 32 (2013), 1473–1480. [35] R. Lamour, F. Mazzia, Computation of consistent initial values for properly stated index 3 DAEs, BIT Num. Math. 49 (2009), 161–175. [36] J. Lampe, J. Bassoy, J. Rahmer, J. Weizenecker, H. Voss, B. Gleich, J. Borgert, Fast reconstruction in magnetic particle imaging, Phys. Med. Biol. 57 (2012), 1113–1134. [37] S. Masnou, J.M. Morel, Level lines based disocclusion, Proceedings 1998 International Conference on Image Processing, vol. 2, IEEE, 1998, pp. 259–263. [38] F. Mazzia, A. Sestini, The BS class of Hermite spline quasi-interpolants on nonuniform knot distributions, Bit Numer. Math. 49 (2009), 611–628. [39] F. Mazzia, A. Sestini, Quadrature formulas descending from BS Hermite spline quasi-interpolation, J. Comput. Appl. Math. 236 (2012), 4105–4118. 480 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. [40] R.M. Mininni, A. Miranville, S. Romanelli, Higher-order Cahn-Hilliard equations with dynamic boundary conditions, J. Math. Anal. Appl. 449 (2017), 1321–1339. [41] L. Modica, S. Mortola, γ-convergence limit for a family of elliptic functionals, Boll. Un. Mat. Ital. A (5) 14 (1977), 526–529 [in Italian]. [42] L. Modica, S. Mortola, An example of γ-convergence, Boll. Un. Mat. Ital. B (5) 14 (1977), 285–299 [in Italian]; Engl. transl. in: J. Jost, X. Li-Jost, Calculus of Variations, Cambridge Studies in Advanced Mathematics, Cambridge University Press, 1998, pp. 248–256. [43] J.M. Morel, S. Solimini, Variational methods in image segmentation: with seven image processing experiments, Progress in nonlinear differential equations and their applications, Birkhäuser, 1995. [44] D. Mumford, J. Shah, Optimal approximations by piecewise smooth functions and associated variational problems, Commun. Pur. Appl. Math. 42 (1989), 577–685. [45] R. Orendorff, D. Hensley, J. Konkle, P. Goodwill, S. Conolly, Computationally efficient image reconstruction via optimization for X-space MPI, 5th International Workshop on Magnetic Particle Imaging (IWMPI), IEEE, 2015. [46] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE T. Pattern. Anal. 12 (1990), 629–639. [47] L.I. Rudin, S. Osher, Total variation based image restoration with free local constrains, IEEE Comput. Soc. Press 1st International Conference on Image Processing, 1994, 31–35. [48] L.I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992), 259–268. [49] C.B. Schönlieb, Higher-order total variation inpainting, https://it.mathworks.com/ matlabcentral/fileexchange/34356-higher-order-total-variation-inpainting, 2012 (accessed January 23, 2017). [50] C.B. Schönlieb, Partial differential equation methods for image inpainting, Cambridge Monographs on Applied and computational Mathematics, Cambridge University Press, 2015. [51] C.B. Schönlieb, A. Bertozzi, Unconditionally stable schemes for higher order inpainting, Commun. Math. Sci. 9 (2011), 413–457. [52] J. Shen, T.F. Chan, Non-texture inpainting by curvature-driven diffusions, J. Vis. Commun. Image R. 12 (2001), 436–449. [53] J. Shen, T.F. Chan, Mathematical models for local nontexture inpaintings, SIAM J. Appl. Math. 62 (2002), 1019–1043. [54] X.C Tai, J. Hahn, J. Chung, A fast algorithm for Euler’s elastica model using augmented Lagrangian method, SIAM J. Imaging Sci. 4 (2011), 313–344. [55] A. Tateo, A. Andrisani, A. Iurino, G. Settanni, P.F. Stifanelli, P. Larizza, F. Mazzia, R.M. Mininni, S. Tangaro, R. Bellotti, A hybrid approach for FFP velocity gridding in MPI reconstruction, 5th International Workshop on Magnetic Particle Imaging (IWMPI), 2015. Applications of PDEs inpainting to magnetic particle imaging. . . 481 [56] A. Tateo, A. Iurino, G. Settanni, A. Andrisani, P.F. Stifanelli, P. Larizza, F. Mazzia, R.M. Mininni, S. Tangaro, R. Bellotti, Hybrid X-space: a new approach for MPI reconstruction, Phys. Med. Biol. 61 (2016), 4061–4077. [57] L. Vese, A study in the bv space of a denoising-deblurring variational problem, Appl. Math. Opt. 44 (2001), 131–161. [58] Z. Wang, A.C. Bovik, H.R. Sheikh, E.P. Simoncelli, Image quality assessment: From error visibility to structural similarity, IEEE T. Image Process. 13 (2004), 600–612. [59] Y. Wu, Z. Yao, G. Kafka, D. Farrell, M. Griswold, R. Brown, The effect of relaxation of magnetic particle imaging, Magnetic Nanoparticles. Proceedings of the First International Workshop on Magnetic Particle Imaging, T.M. Buzug et al. (eds.), World Scientific Publishing, 2010, pp. 113–119. Andrea Andrisani andrea.andrisani@uniba.it Università degli Studi di Bari Aldo Moro Dipartimento di Matematica 70125 Bari, Italy Rosa Maria Mininni rosamaria.mininni@uniba.it Università degli Studi di Bari Aldo Moro Dipartimento di Matematica 70125 Bari, Italy Francesca Mazzia francesca.mazzia@uniba.it Università degli Studi di Bari Aldo Moro Dipartimento di Informatica 70125 Bari, Italy Giuseppina Settanni giuseppina.settanni@uniba.it Università degli Studi di Bari Aldo Moro Dipartimento di Matematica 70125 Bari, Italy 482 A. Andrisani, R.M. Mininni, F. Mazzia, G. Settanni, A. Iurino, S. Tangaro, et al. Alessandro Iurino aiurino@hotmail.it Università degli Studi di Bari Aldo Moro Dipartimento Interateneo di Fisica “M. Merlin” 70125 Bari, Italy Sabina Tangaro sonia.tangaro@ba.infn.it Istituto Nazionale di Fisica Nucleare sezione di Bari, 70125 Bari, Italy Andrea Tateo andrea.tateo@ba.infn.it Università degli Studi di Bari Aldo Moro Dipartimento Interateneo di Fisica “M. Merlin” 70125 Bari, Italy Roberto Bellotti roberto.bellotti@uniba.it Università degli Studi di Bari Aldo Moro Dipartimento Interateneo di Fisica “M. Merlin” 70125 Bari, Italy Istituto Nazionale di Fisica Nucleare sezione di Bari, 70125 Bari, Italy Received: January 15, 2019. Accepted: January 29, 2019.