# Application of the topological gradient method to tomography

## Application du gradient topologique à la tomographie

D. Auroux* — L. Jaafar-Belaid** — B. Rjaibi**

**ABSTRACT.** A new method for parallel beam tomography is proposed. This method is based on the topological gradient approach. The use of the topological asymptotic analysis for detecting the main edges of the data allows us to filter the noise while inverting the Radon transform. Experimental results obtained on noisy data illustrate the efficiency of this promising approach in the case of Magnetic Resonance Imaging. We also study the sensitivity of the algorithm with respect to several regularization and weight parameters.

**RÉSUMÉ.** Une nouvelle méthode de reconstruction pour la tomographie par faisceaux parallèles est proposée. Cette méthode est basée sur l'approche du gradient topologique. En détectant les contours sur les données grâce à l'analyse asymptotique topologique, il est possible de filtrer le bruit dans le processus d'inversion de la transformée de Radon. Des résultats expérimentaux obtenus sur des données bruitées illustrent les possibilités de cette approche prometteuse dans le domaine de traitement d'images IRM. Nous étudions également la sensibilité de l'algorithme par rapport aux différents paramètres de régularisation et pondération.

**KEYWORDS:** Topological asymptotic expansion, topological gradient, tomographic reconstruction, tomography, Radon transform.

**MOTS-CLÉS:** Développement asymptotique topologique, gradient topologique, reconstruction tomographique, tomographie, transformation de Radon.

<footnote>* Laboratoire J. A. Dieudonné,
Université de Nice Sophia Antipolis,
Parc Valrose, 06108 Nice cedex 2,
France.
auroux@unice.fr</footnote>

<footnote>** ENIT-LAMSIN, Campus Universitaire,
BP 37, 1002 Le Belvédère, Tunis, Tunisie.
lamia.belaid@esstt.rnu.tn
badreddine.rjaibi@lamsin.rnu.tn</footnote>
# 1. Introduction

The tomographic problem can be seen as the reconstruction of an object from measurements which are line integrals of the processed object at some given orientations (or view angles). This problem has many applications including medical imaging, electron microscopy, radio astronomy, ... [11]. This reconstruction highly depends on the amount of available line integrals. Classically, when line integrals are available from many directions and when the measurements are nearly noise free, the filtered backprojection reconstruction gives good results. Unfortunately, these conditions are not possible in real life applications, and the interpretation of results will be degraded. Another drawback is the computational cost required by traditional methods for tomography. Currently, the fastest algorithms require $O(n^2 \log n)$ operations for reconstructing $n^2$ pixels.

The goal of this paper is to present a new method for the tomographic problem. This method is based on the topological gradient approach, that has been introduced for topological optimization purpose [1, 2, 12, 13, 16, 19]. The idea of topological optimization is to create a partition of a given domain (in our case, an image) into an optimal design and its complementary. A common way to consider the restoration problem is to identify the edges of the image, in order to preserve them in the restoration process. Then, the image is smoothed elsewhere. This technique was successfully used for several problems in image processing [4, 5, 6, 14]. The computational time was very interesting too. On the other hand, it has been shown in [14] that a classical way to restore a given image $f$ is exactly the Tikhonov regularisation applied to the inversion of a compact operator. For all these reasons, the topological gradient approach seems appropriate for the tomographic problem, which is well known to be an ill-posed problem. There is indeed no continuity of the reconstructed image with respect to the data.

The Radon transform represents a set of parallel line integral projections of a 2D function $f$ at different angles $\theta$, where $f$ is the processed image, defined on a bounded open and convex domain $\Omega$ of $\mathbb{R}^2$. These projections (or sinogram) are given by:

$$ \mathcal{R}(f)(\theta, r) = \iint f(x, y)\delta(r - x \cos \theta - y \sin \theta)dxdy, $$

where $r$ and $\theta$ are the polar coordinates and $\delta(\cdot)$ is the Dirac delta function. We refer to [8, 11, 15, 17, 20] for the historical development of tomography and some review of all the methods proposed in the literature for tomographic reconstruction.

Inverting the Radon transform is quite easy if the data are unnoisy. There exists indeed explicit inversion formulas, based on the Fourier transform or the backprojection operator. But in the presence of noisy data, even for a very small level of noise, these formulas do not allow one to recover a good estimation of the initial image. Thus it is necessary to consider a way to filter the noise. A simple but efficient algorithm is the Filtered Back-Projection (FBP) reconstruction method. This is probably the most widely used technique for inverting the 2D Radon transform, as it is very fast and quite accurate. The idea is to filter some Fourier modes of the projections, and then to invert the Radon transform [17].

Another way to deal with noise consists in relaxing the problem: if we denote by $g$ the noisy data, the idea is to regularize the inversion problem $\mathcal{R}(f) = g$, and to solve it in a least-square sense:

$$ \min_{f} \|\mathcal{R}(f) - g\|^2 + \lambda\phi(f), $$

where $\phi$ is a regularization function on the reconstructed image $f$, and $\lambda$ is a regularization coefficient. Several regularization functions have been studied in order to preserve the
main edges of the image. For instance, one can consider the total variation norm $\phi(f) = \int_{\Omega} |\nabla f| dx$ (see e.g. [18]). Also, it is possible to consider different regularizations on homogeneous zones and image edges: $\phi(f) = \int_{\Omega} |\nabla f|^{\alpha} dx$, with $\alpha = 2$ in homogeneous zones and $\alpha = 1$ on the edges [9]. This method leads to a piecewise homogeneous image with sharp edges. However, the main issue is to find the edge set of the image. We use then the topological asymptotic analysis for this purpose, in order to have an efficient regularization term and to filter the noise, but also to preserve the main characteristics of the data, while inverting the operator $\mathcal{R}$.

This paper is organized as follows: Section 2 is devoted to a quick presentation of a standard filtered backprojection reconstruction method. The topological gradient approach and its application to tomography is developed in section 3. In section 4, we report the results of many numerical experiments and comparisons. We also discuss the complexity of our algorithm. Finally, some conclusions and perspectives are given in section 5.

## 2. A filtered backprojection reconstruction method

We describe briefly in this section a Fourier reconstruction algorithm based on the well known Fourier slice theorem [17]. We denote by $\mathcal{R}$ the Radon transform and $\mathcal{R}^{\#}$ the backprojection operator. Then, for a given function $f$,

$$ \mathcal{R}^{\#} \mathcal{R} f = f * h, \quad \text{where } h(x, y) = \frac{1}{\sqrt{x^2 + y^2}}. $$

One has then to deconvolute the backprojection, in order to estimate the solution $f$. This deconvolution is often performed with the Fourier transform (and its inverse). This method is quite easy to implement, but it does not produce good results in the presence of noisy data. A very simple idea is then to add a filter step in the backprojection method. More precisely, a filtered backprojection reconstruction method processes into the following steps:

1) Compute the 1D Fourier transform $\mathcal{F}_r$ of the projections in the variable $r$:

$$ \mathcal{F}_r \mathcal{R}(f)(r, \theta). $$

2) Use a ramp filter for filtering the projections:

$$ \tilde{f}(\rho, \theta) = \mathcal{F}_r^{-1} |r| \mathcal{F}_r \mathcal{R}(f)(r, \theta), $$

where $\mathcal{F}_r^{-1}$ denotes the inverse Fourier transform in the variable $r$.

3) Apply the continuous backprojection (or adjoint) operator $\mathcal{R}^{\#}$:

$$ f(x, y) = \mathcal{R}^{\#} \tilde{f}(\rho, \theta) = \int_{0}^{\pi} \tilde{f}(x\cos\theta + y\sin\theta)d\theta. $$

The FBP reconstruction algorithm provides a reconstructed image of good quality in the case of sufficient noise-free data. However, these conditions are rarely available in real life applications, inducing lower visual quality of the processed image.
### 3. Application of the topological gradient method to computed tomography

In this section, we use the topological gradient as a tool for the reconstruction problem in tomography.

#### 3.1. Regularized problem

We first recall that a standard approach for regularizing the ill-posed problem of tomographic imaging consists in the following optimization problem:

$$ \min_{f} \|Af - g\|^2 + \lambda\phi(f), $$

where A denotes a system matrix defining the discrete Radon transform, g is the measured projection, φ is a regularization functional and λ represents a parameter which controls the tradeoff between a good fit to the data and the smoothness of the solution. Inspired by [9, 10], in which the authors introduced the semi-norm of the total variation (TV) which is particularly efficient in recovering piecewise smooth functions without smoothing sharp discontinuities, we propose to consider the following minimization problem:

$$ \min_{f} \int_{\Omega} |Af - g|^2 dx + c_0 \int_{\Omega} |\nabla f|^{\alpha} dx, $$

where $c_0$ is a positive constant and $\alpha$ is equal to 1 on the edges and to 2 elsewhere.

The main issue is then to identify the edge set in an efficient and correct way.

#### 3.2. Topological gradient approach

For this purpose, the topological asymptotic analysis has showed its efficiency [1, 5, 14]. We now detail the application of this method to the current problem.

The idea is to detect the edges of the image in order to preserve them during the reconstruction process. The edges can be modelled by cracks in the standard thermal diffusion approach. These cracks are supposed to be highly insulating and they do not allow the temperature (or image intensity) to jump across the edges.

For a small $\rho \ge 0$, let $\Omega_\rho = \Omega \setminus \sigma_\rho$ be the perturbed domain by the insertion of a crack $\sigma_\rho = x_0 + \rho\sigma(n)$, where $x_0 \in \Omega$, $\sigma(n)$ is a straight crack, and $n$ is a unit vector normal to the crack. For a given function g, we consider the following problem:

$$ \begin{cases} -\mathrm{div}(c_0 \nabla f_\rho) + A^\# A f_\rho = A^\# g & \text{in } \Omega_\rho, \\ \partial_n f_\rho = 0 & \text{on } \partial\Omega_\rho, \end{cases} \quad [1] $$

where $n$ denotes the outward unit normal to $\partial\Omega_\rho$ and $A^\#$ is the dual discrete Radon transform. Equation 1 is the optimality system of the following penalized least square criterion:

$$ \int_{\Omega_\rho} c_0 |\nabla f_\rho|^2 dx + \int_{\Omega_\rho} |Af_\rho - g|^2 dx. $$

The minimization of this criterion allows us to invert the discrete Radon transform. To preserve edges as much as possible, we look for the leading term of

$$ j(\rho) = J(f_\rho) = \int_{\Omega_\rho} |\nabla f_\rho|^2 dx. $$
By considering the solution $v$ to the adjoint problem

$$ \begin{cases} -\mathrm{div}(c_0 \nabla v) + A^{\#}Av = -\partial_f J(f) & \text{in } \Omega, \\ \partial_n v = 0 & \text{on } \partial\Omega, \end{cases} \quad [2] $$

we obtain the following topological asymptotic expansion in the case of a thin crack with boundary condition $\partial_n u = 0$ on $\partial\sigma(n)$ [1, 6, 14]:

$$ j(\rho) - j(0) = \rho^2 G(x_0, n) + o(\rho^2), $$

with

$$ G(x_0, n) = -\pi c_0 (\nabla f(x_0).n)(\nabla v(x_0).n) - \pi |\nabla f(x_0).n|^2. $$

The topological gradient reads

$$ G(x, n) = \langle M(x)n, n \rangle $$

where $M(x)$ is the symmetric matrix defined by

$$ M(x) = -\pi c_0 \frac{\nabla f(x) \nabla v(x)^T + \nabla v(x) \nabla f(x)^T}{2} - \pi \nabla f(x) \nabla f(x)^T. \quad [3] $$

For a given $x$, $G(x, n)$ takes its minimal value when $n$ is the eigenvector associated to the lowest eigenvalue $\lambda_{min}$ of $M$. This value will be considered as the topological gradient. The edge set is then defined as follows:

$$ \sigma = \{x \in \Omega, \lambda_{min} < \alpha_0 < 0\}, $$

where $\alpha_0$ is a negative threshold. This threshold can either be tuned on a set of test images, or be defined such that a given percentage of the pixels are identified as edges. It can also be replaced by an automatic identification of the edges by using valley lines of the topological gradient [7].

The main advantage of this approach is to detect the edge set $\sigma$ in a very fast and efficient way. One has indeed to solve only two problems: the direct unperturbed problem (equation 1 with $\rho = 0$), and the adjoint problem (equation 2). Then, the matrix $M(x)$ (given by equation 3) can be assembled very quickly, and it directly provides the edge set.

We have then to calculate the reconstruction solution by solving the following perturbed problem:

$$ \begin{cases} -\mathrm{div}(c_1(x_0)\nabla f) + A^{\#}Af = A^{\#}g & \text{in } \Omega, \\ \partial_n f = 0 & \text{on } \partial\Omega, \end{cases} \quad [4] $$

with

$$ c_1(x_0) = \begin{cases} \frac{c_0}{|\nabla f(x_0)|}, & x_0 \in \sigma = \{x \in \Omega, \lambda_{min} < \alpha_0 < 0\}, \\ c_0 & \text{elsewhere, where } c_0 \text{ is a constant.} \end{cases} \quad [5] $$

In the homogeneous parts of the image (i.e. not on the edges), we use a constant regularization coefficient, in order to smooth the image. And in the edge set, we adapt the coefficient in order to regularize the problem with a TV (total variation) norm.
# 4. Numerical experiments

We now report the results of several numerical experiments, in order to study the efficiency of our approach.

## 4.1. Shepp-Logan phantom

We have tested our method on the well-known Shepp-Logan head phantom, for which projections data have been computed using a discrete Radon transform.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Reconstruction of the Shepp-Logan head phantom using the topological gradient method: (a) original image, (b) noised sinogram (SNR=24.5), (c) reconstructed image (PSNR=26.18).</figcaption></figure>

Figure 1 (a) shows the original image. A Gaussian noise is added to the data with a signal to noise ratio *SNR* = 24.5. The sinogram is represented in figure 1(b). Figure 1(c) shows the reconstruction result (Peak Signal to Noise Ratio: *PSNR* = 26.18). The reconstructed image is visually excellent, as it is possible to see all different ellipses. The threshold coefficient *α*₀ used for defining the edges is set to *α*₀ = -0.025 in this experiment.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Reconstruction of the Shepp-Logan head phantom using the FBP and standard TV approaches: (a) reconstructed image using FBP (PSNR=14.59), (b) FBP reconstruction followed by a Hamming filter (PSNR=15.91), (c) reconstructed image using the TV method (PSNR=22.43).</figcaption></figure>

In order to visualize the performance of our approach, we propose to give numerical results for standard FBP and TV approaches. The first method is shown in figure 2(a), and the corresponding PSNR is 14.59. In order to attenuate the distorsion observed in the reconstructed image, we apply then a Hamming filter on the image, but this technique over-smoothes the edges as seen in figure 2(b), and the visual quality of the result is then
highly reduced, even if the PSNR is slightly improved. A total variation-based reconstruction for computed tomography based on [20] has also been tested and the result is shown in figure 2(c), and PSNR = 22.43 in this case. These results confirm the fact that the image should be smoothed outside the edges, in order to recover a piecewise smooth image, while keeping the TV norm for preserving the edges.

## 4.2. Sensitivity to data noise

We now study the sensitivity of our method to the level of noise of the data. We added more and more noise to the sinogram, corresponding to a signal to noise ratio going from 24.5 to 15.5.

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Reconstruction of the Shepp-Logan head phantom using the topological gradient method. Reconstructed images for several noise levels: (a) SNR=22.5; (b) SNR=20; (c) SNR=15.5.</figcaption></figure>

One can see on figures 1(c) and 3(a)-(c) the evolution of the reconstructed image when the level of noise is increased. For a relatively high level of noise (SNR=15.5 on the sinogram), the reconstructed image is still of good quality, as most of the ellipses are visible.

Figure 4 shows the quality of the reconstructed image as a function of the level of noise in the sinogram. The noise level is represented by its signal to noise ratio (SNR). The quality of the reconstructed image is quantified by both its SNR and PSNR. One can see that the reconstructed image loses only 5 dB while nearly 10 dB of noise is added to the data. This shows that the topological gradient based method that we introduced in this paper is not very sensitive to the presence of noise on the data.

## 4.3. Comparison with other methods

We now compare more quantitatively the results obtained by several methods. We have also adapted the topological gradient method for edge detection to a full TV model, by considering the following regularization function instead of equation 5:

$$c_1(x_0) = \begin{cases} \frac{\varepsilon}{|\nabla f(x_0)|}, & x_0 \in \sigma, \\ \frac{c_0}{|\nabla f(x_0)|}, & \text{elsewhere,} \end{cases} \quad [6]$$

where $\varepsilon > 0$ is a very small coefficient compared with $c_0$. This is then equivalent to considering the standard total variation norm, with a very small regularization coefficient on the edge set (in order to preserve the edges), compared to the smooth parts of the image. In the following, we will consider the topological gradient methods based on these two
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Evolution of the PSNR (and SNR) of the reconstructed image by the topological gradient method, as a function of the level of noise (SNR) in the sinogram. different models, either TV or $L^1/L^2$ norms, where the regularization function is given by equations 6 and 5 respectively. We will also show the results given by the filtered backprojection method (with or without a Hamming filter).</figcaption></figure>

<table><caption>**Table 1.** Comparison between several reconstruction methods, quantified by several ratios and similarity indexes, for a noisy sinogram with SNR=24.5.</caption><thead><tr><td>Method</td><td>PSNR</td><td>SNR</td><td>SSIM</td><td>MSE</td></tr></thead><tbody><tr><td>FBP</td><td>14.59</td><td>33.02</td><td>0.41</td><td>0.0342</td></tr><tr><td>FBP + Hamming</td><td>15.91</td><td>31.71</td><td>0.52</td><td>0.0273</td></tr><tr><td>Topological gradient (TV)</td><td>22.43</td><td>33.75</td><td>0.82</td><td>0.0042</td></tr><tr><td>Topological gradient ($L^1/L^2$)</td><td>26.18</td><td>34.01</td><td>0.94</td><td>0.0023</td></tr></tbody></table>

Table 1 shows the PSNR (Peak Signal to Noise Ratio, in dB), SNR (Signal to Noise Ratio, in dB), SSIM (Structural Similarity) and MSE (Mean Squared Error) of the reconstructed images for several methods: FBP (Filtered BackProjection), FBP+Hamming (with a Hamming filter on the FBP image), and the topological gradient approach using either the TV (Total Variation) norm, or both TV and $L^2$ norms (see equation 5). The sinogram has a noise level such that the SNR is 24.5.

We can see in this table that the best results for all indicators are obtained with the hybrid $L^1/L^2$ norm within the topological gradient approach. The peak SNR is higher of more than 3 decibels, the structural similarity index is larger than 94% and the mean squared error is smaller than 1%. Note that both topological gradient methods provide much better results than the filtered backprojection, which is not surprising.
<table><caption>**Table 2.** Comparison between several reconstruction methods, quantified by several ratios and similarity indexes, for a noisy sinogram with SNR=20.</caption><thead><tr><th>Method</th><th>PSNR</th><th>SNR</th><th>SSIM</th><th>MSE</th></tr></thead><tbody><tr><td>FBP</td><td>14.32</td><td>32.94</td><td>0.30</td><td>0.0370</td></tr><tr><td>FBP + Hamming</td><td>15.56</td><td>31.72</td><td>0.46</td><td>0.0277</td></tr><tr><td>Topological gradient (TV)</td><td>19.05</td><td>32.37</td><td>0.69</td><td>0.0124</td></tr><tr><td>Topological gradient (L<sup>1</sup>/L<sup>2</sup>)</td><td>24.81</td><td>33.92</td><td>0.85</td><td>0.0033</td></tr></tbody></table>

Table 2 shows the same results as table 1, but for a higher level of noise: the SNR of the sinogram is now 20. This table confirms the relative insensitivity of the topological gradient method with respect to noise, as shown in the previous subsection. The PSNR has indeed been decreased of less than 2 dB, the SSIM index is close to 1, and the MSE is still smaller than 1%. The most noticeable point is that the method is much less degraded when the level of noise is increased if one uses the hybrid $L^1/L^2$ norm. The interest of the topological gradient is clearly to identify the edge set, allowing us to use at best the two different regularizations on the solution.

## 4.4. Sensitivity to the regularization coefficient

We now study the evolution of the solution when we tune the regularization coefficient $c_0$ in the topological gradient method (see e.g. equation 5).

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Evolution of the PSNR (in dB) of the reconstructed image by the topological gradient method, as a function of the regularization coefficient $c_0$.</figcaption></figure>
Figure 5 shows the evolution of the PSNR of the reconstructed image as a function of $c_0$. We can see that there is a plateau, between $c_0 = 5$ and 30 nearly, for which the solution is of best quality. For a quite larger value of the regularization coefficient (up to 60) or a smaller one (around 1), the PSNR of the solution is slightly reduced (by 1 to 2 decibels). When the level of noise is increased, the shape of the PSNR curve is exactly the same, but the plateau is translated to larger values of $c_0$. This is not surprising, as in all optimization problems, the regularization coefficient is directly related to the level of noise on the data. But as the plateau is relatively wide, it is quite easy to define an a priori regularization coefficient that will lead to a good reconstructed image.

## 4.5. Numerical results on other images

We finally present some results on other images.

Figure 6 shows a comparison between the previous approaches (filtered backprojection; FBP + Hamming filter; topological gradient + TV; topological gradient + $L^1/L^2$) on a brain image. We can see that the last two algorithms produce very similar results, from both visual point of view and PSNR/SSIM values. The images obtained with the filtered backprojection methods are not satisfactory, as they are either noisy or blurred. In this example, considering two different regularizations (one on the edge set and one in the smooth parts of the image) does not really improve the results. This is possibly due to the fact that the edges are not as sharp as in the previous synthetic example.

Figure 7 shows the same comparison, for another image. The conclusions are almost the same: the topological gradient methods give much better results than the filtered backprojection methods, and there is a slight improvement if one considers a mixed $L^1/L^2$ norm for the reconstruction, rather than a TV norm.

From these experiments, we can conclude that one should consider the $L^1/L^2$ norm for solving the reconstruction problem, as defined in equation 5, as this sometimes improves a lot the quality of the reconstructed image.

## 4.6. Computation cost

We now briefly present the computation cost of the topological gradient method for inverting the Radon transform. Note that the two approaches that we considered (TV and $L^1/L^2$ norms) have a similar cost.

As previously explained, one has to solve two problems in order to define the topological gradient: one direct and one adjoint PDE problem. These problems are linear, and quite easy to solve, provided the operator $A^#A$ is assembled (see e.g. equation 1). Then, once the topological gradient defined, and the edge set identified, one has to solve the perturbed problem (see equation 4). But then, as the edge set is relatively small, we solve this perturbed problem with a conjugate gradient method, preconditioned by the operator of the unperturbed PDE problem. This makes the preconditioned operator close to the identity.

We also use a discrete cosine transform (or, in an equivalent way, a Fourier transform) to solve the linear systems. The topological gradient algorithm requires then $O(n^2 \log(n))$ operations for an image with $n \times n$ pixels. If we compare with other iterative methods, the topological gradient method is both efficient and fast. We refer to [5, 14] for more details about the complexity of such topological gradient algorithms.
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Reconstruction of a brain image: (a) original image; (b) noisy sinogram (SNR=24); reconstructed images by: (c) FBP (PSNR=24.19, SSIM=0.31); (d) FBP + Hamming filter (PSNR=26.71, SSIM=0.61); (e) topological gradient + TV norm (PSNR=33.47, SSIM=0.85); (f) topological gradient + $L^1/L^2$ norms (PSNR=33.38, SSIM=0.848).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Reconstruction of a hand image: (a) original image; (b) noisy sinogram (SNR=24); reconstructed images by: (c) FBP (PSNR=14.34, SSIM=0.39); (d) FBP + Hamming filter (PSNR=16.77, SSIM=0.56); (e) topological gradient + TV norm (PSNR=24.32, SSIM=0.73); (f) topological gradient + $L^1/L^2$ norms (PSNR=24.53, SSIM=0.73).</figcaption></figure>
## 5. Conclusion

We have presented in this work a new method for the reconstruction problem from 2D tomographic data. It is based on the topological asymptotic analysis for identifying the edge set of the image. Then, the edges are preserved and the other parts of the image are smoothed during the inversion process.

To make this method relevant for real life applications, we have considered a spectral approach, based on the discrete cosine transform, and a preconditioned conjugate gradient method for solving the equations. In this case, the topological gradient algorithm requires $O(n^2 \log n)$ operations, where $n$ is the size of the image.

The numerical results show reconstructed images of good visual quality compared to the standard FBP and TV methods. Our approach also produces images of higher PSNR.

We have shown the relative insensitivity of this approach to the noise level, as the topological gradient algorithm appeared to be much less sensitive than a TV method for instance. We have also shown that it is possible to find an a priori estimation of the regularization coefficient, as the algorithm is not sensitive to variations around the optimal coefficient.

This kind of method is currently under investigation in the 3D case.

## 6. References

[1] AMSTUTZ S., HORCHANI I., MASMOUDI M., "Crack detection by the topological gradient method", *Control and Cybernetics*, vol. 34, num. 1, pp. 119-138, 2005.

[2] AMSTUTZ S., MASMOUDI M., SAMET B., "The topological asymptotic for the Helmholtz equation", *SIAM J. Control Optim.*, vol. 42, num. 5, pp. 1523-1544, 2003.

[3] AUBERT G., KORNPROBST P., "Mathematical Problems in Image Processing", *Applied Mathematical Sciences, Springer Verlag*, vol. 147, 2001.

[4] AUROUX D., "From restoration by topological gradient to medical image segmentation via an asymptotic expansion", *Math. Comput. Model.*, vol. 49(11-12), pp. 2191-2205, 2009.

[5] AUROUX D., JAAFAR BELAID L., MASMOUDI M., "Image restoration and classification by topological asymptotic expansion", *Variational Formulations in Mechanics: Theory and Applications*, CIMNE, Barcelona, Spain, pp. 23-42, 2007.

[6] AUROUX D., JAAFAR BELAID L., MASMOUDI M., "A topological asymptotic analysis for the regularized gray-level image classification problem", *Math. Model. Numer. Anal.*, vol. 41, pp. 605-625, 2007.

[7] AUROUX D., MASMOUDI M., "Image processing by topological asymptotic expansion", *J. Math. Imaging Vis.*, vol. 33, pp. 122-134, 2009.

[8] BASU S., BRESLER Y., "$O(N^2 \log_2 N)$ filtered backprojection reconstruction algorithm for tomography", *IEEE Trans. Image Proc.*, vol. 9, num. 10, pp. 1760-1773, 2000.

[9] BLOMGREN P., CHAN T., MULET P., WONG C. K., "Total variation image restoration: numerical methods and extensions", *Proc. Int. Conf. on Image Processing*, pp. 384-387, 1997.

[10] CHAN T., MARQUINA A., MULET P., "High order total variation based image restoration", *SIAM J. Sci. Comput.*, vol. 22, num. 2, pp. 503-516, 2000.

[11] DEAN S. R., "The Radon transform and some of its applications", *Wiley, New York*, 1983.

[12] GARREAU S., GUILLAUME P., MASMOUDI M., "The topological asymptotic for PDE systems: the elasticity case", *SIAM J. Control Optim.*, vol. 39, pp. 17-49, 2001.
[13] HASSINE M., MASMOUDI M., "The topological asymptotic expansion for the Quasi-Stokes problem", *ESAIM: Control Optim. Calc. Var.*, vol. 10, pp. 478-504, 2004.

[14] JAAFAR BELAID L., JAOUA M., MASMOUDI M., SIALA L., "Application of the topological gradient to image restoration and edge detection", *Engineer. Anal. Bound. Elements*, vol. 32, pp. 891-899, 2008.

[15] LEWITT R. M., "Reconstruction algorithms: Transform methods", *Proc. IEEE*, vol. 71, num. 3, pp. 390-408, 1983.

[16] MASMOUDI M., "The Topological Asymptotic, Computanional Methods for Control Appli-cations", *Int. Series GAKUTO*, R. Glowinski, H. Kawarada and J. Périaux (Eds.), Tokyo, Japan, vol. 16, pp. 53-72, 2001.

[17] NATTERER F., "The Mathematics of Computerized Tomographic Imaging", *Wiley, New York*, 1986.

[18] RUDIN, L., OSHER S., FATEMI E., "Nonlinear Total Variation based noise removal algo-rithms", *Physica D*, vol. 60, pp. 259-268, 1992.

[19] SOKOLOWSKI J., ZOCHOWSKI A., "Topological derivatives of shape functionals for elastic-ity systems", *Int. Ser. Numer. Math.*, vol. 139, pp. 231-244, 2002.

[20] ZHANG X. Q., FROMENT J., "Total variation based Fourier reconstruction and regularization for computer tomography", *IEEE Nuclear Science Symp. Conf. Record*, 2005.