Medical scanner is the most used application of tomographic reconstruction. It allows to explore the interior of a human body. In the same way, industrial tomography explores the interior of an object and is often used for non-destructive testing.
We are interested here in a very specific application of tomographic reconstruction for a physical experiment described later. The goal of this experiment is to study the behavior of a material under a shock. We obtain during the deformation of the object an X-ray radiography by high speed image capture. We suppose this object is radially symmetric, so that one radiograph is enough to reconstruct the 3D object.
Several authors have proposed techniques (standard in medical tomography) for tomographic reconstruction when enough projections (from different points of view) are available: this allows to get an analytic formula for the solution (see for instance  or ). These methods cannot be used directly when only a few number of projections is known. Some alternative methods have been proposed in order to partially reconstruct the densities (see for instance ). We are interesting here in single view tomographic reconstruction for radially symmetric object (see for instance  for a more complete presentation of the subject). As any tomographic reconstruction, this problem leads to an ill-posed inverse problem. As we only have one radiograph, data are not very redundant and the ill-posed character is even more accurate.
We present here a tomographic method adapted to this specific problem, originally developed in 
, and based on a curve evolution approach. The main idea is to add some a priori knowledge on the object we are studying in order to improve the reconstruction. The object may then be described by a small set of characters (in this case, they will be curves) which are estimated by the minimization of an energy functional. This work is very close to another work by Feng and al. The main difference is the purpose of the work: whereas they are seeking recovering textures, we are looking for accurate edges. It is also close to the results of Bruandet and al . However, the present work handles very noisy data and highly unstable inverse problems, and shows how this method is powerful despite these perturbations. Further, we take here into account the effects of blur (which may be non-linear) and try to deconvolve the image during the reconstruction.
Let us mention at this point that our framework if completely different from the usual tomographic point of view, and usual techniques (such as filtered back-projection) are not adapted to our case. Indeed, usually, as the X-rays are supposed to be parallel (this is also the case here), the “horizontal” slices of the object are supposed to be independent and studied separately. Usual regularization techniques deal with one slice and regularize this particular slice. Here, because of the radial symmetry, the slices are composed of concentric annulus and do not need any regularization. The goal of this work is to add some consistency between the slices in order to improve the reconstruction.
The paper is organized as follows. First we present the physical experiment whose data are extracted and explain what are the motivations of the work. Next, we introduce the projection operator. In Section 4, we present a continuous model with the suitable functional framework and prove existence result. Section 5 is devoted to formal computation of the energy derivative in order to state some optimality conditions. In Section 6, a front propagation point of view is adopted and the level set method leads to a non local Hamilton-Jacobi equation. In the last section, we present some numerical results and give hints for numerical schemes improvement.
This work is part of some physical experiments whose goal is the study of the behavior of shocked material. The present experiment consists in making a hull of well known material implode using surrounding explosives. The whole initial physical setup (the hull, the explosives …) are radially symmetric. A reasonable assumption is to suppose that during the implosion, everything remains radially symmetric.
Physicists are looking for the shape of the interior at some fixed time of interest. At that time, the interior may be composed of several holes which also may be very irregular. Figure 3 is a synthetic object that contains all the standard difficulties that may appear. These difficulties are characterized by:
Several disconnected holes.
A small hole located on the symmetry axis (which is the area where the details are difficult to recover).
Smaller and smaller details on the boundary of the top hole in order to determine a lower bound detection.
To achieve this goal, a X-rays radiograph is obtained. In order to extract the desired informations, a tomographic reconstruction must be performed. Let us note here that, as the object is radially symmetric, a single radiography is enough to compute the reconstruction.
A radiography measures the attenuation of X-rays through the object. A point on the radiography will be determined by its coordinates in a Cartesian coordinates system where the -axis will be the projection of the symmetry axis. If is the intensity of the incident X-rays flux, the measured flux at a point is given by
where the integral operates along the ray that reaches the point of the detector, is the infinitesimal element of length along the ray and is the local attenuation coefficient. For simplicity, we will consider that this coefficient is proportional to the material density. To deal with linear operators, we take the Neperian logarithm of this attenuation and will call the transformation
the projection operator.
Through the rest of the paper, in order to simplify the expression of the projection operator, we will suppose that the X-ray source is far enough away from the object so that we may consider that the rays are parallel, and orthogonal to the symmetry axis. As a consequence, the horizontal slices of the object may be considered separately to perform the projection.
As the studied object is radially symmetric, we will work in a system of cylinder coordinates where the -axis is the symmetry axis. The object is then described by the density at the point , which is given by a function which depends only on by symmetry. In the text, the notation will always refer to the density of the object. A typical function is given in Figure 2. It represents an object composed of concentric shells of homogeneous materials (called the “exterior” in what follows) surrounding a ball (called the “interior”) of another homogeneous material that contains some empty holes. This figure may be viewed as a slice of the object by a plane that contains the symmetry axis. To recover the 3D-object, it suffices to perform a rotation of this image around the axis. For instance, the two round white holes in the center are in fact the slice of a torus. As the looked-after characteristic of the object is the shape of the holes, we will focus only on the interior of the object (see Figure 3). We here handle only binary objects composed of one homogeneous material (in black) and some holes (in white).
3. A variational approach
3.1. The projection operator
We first explicit the projection operator and its adjoint.
In the case of a radially symmetric object, the projection operator, denoted by , is given, for every function with compact support, by
Proof - Consider a 3D-object which is described by a function (Cartesian coordinates system). The projection operator is
In the case of radially symmetric object, we parametrize the object by a function with cylinder coordinates. Therefore
Then, we have, for
We perform the following change of variable
For , we have
Using the change of variable and the fact that is even, by symmetry, we get
Operator may be defined by density on measurable functions such that all the partial applications belong to . Then, all the functions belong to the space of bounded mean oscillations functions :
where . For more details, one can refer to .
In the sequel, we will need to handle functions that are defined on (instead of on ). We thus define the operator for function with compact support by
although this has no more physical meaning. Here, the function is defined by
We shall also need the back-projection that is the adjoint operator of ; it can be computed in a similar way.
The adjoint operator (in ) of the projection operator is given, for every function with compact support by :
Proof - The adjoint operator of is the unique operator such that, for every and in with compact support,
Using (3.1) and Fubini’s theorem, we get
So we obtain the following expression for the back projection :
3.2. Toward a continuous model
Thanks to the symmetry, this operator characterizes the Radon transform of the object and so is invertible; one radiograph is enough to reconstruct the object. The inverse operator is given, for an almost everywhere differentiable function with compact support, and for every , by
Because of the derivative term, the operator is not continuous. Consequently, a small variation on the measure leads to significant errors on the reconstruction. As our radiographs are strongly perturbed, applying to our data leads to a poor reconstruction. Due to the experimental setup they are also two main perturbations:
A blur, due to the detector response and the X-ray source spot size.
Others perturbations such as scattered field, motion blur… also exist
but are neglected in this study.
We denote by the effect of blurs. We will consider the following simplified case where is supposed to be linear
where is the usual convolution operation, is the projected image and is a positive symmetric kernel.
A more realistic case stands when the convolution operates on the intensity : then is of the form
where is the multiplicative coefficient between the density and the attenuation coefficient. Some specific experiments have been carried out to measure the blur effect. Consequently, we will suppose that, in both cases, the kernel is known. The linear blur is not realistic but is treated here to make the computations simpler for the presentation.
The noise is supposed for simplicity to be an additive Gaussian white noise of mean 0, denoted by. Consequently, the projection of the object will be
The comparison between the theoretical projection and the perturbed one is shown on Figure 4. The reconstruction using the inverse operator applied to is given by Figure 5. The purpose of the experiment is to separate the material from the empty holes and consequently to precisely determine the frontier between the two areas, which is difficult to perform on the reconstruction of Figure 5.
It is clear from Figure 5 that the use of the inverse operator is not suitable. In order to improve the reconstruction, we must add some a priori knowledge on the object to be reconstructed. Indeed the object that we reconstruct must satisfy some physical property.
We chose to stress on two points:
The center of the object is composed of one homogeneous known material’s density with some holes inside.
There cannot be any material inside a hole.
In a previous work , J.M. Lagrange reconstructed the exterior of the object. In this reconstruction, the density of the material at the center of the object is known, only the holes are not reconstructed. In other words, we can reconstruct an object without holes and we can compute (as and are known) the theoretical projection of this reconstruction. We then act as if the blurred projection was linear and subtract the projection of the non-holed object to the data. In what follows, we will call experimental data this subtracted image which corresponds to the blurred projection of a “fictive” object of density 0 with some holes of known “density” . Consequently, the space of admitted objects will be the set of functions that take values in . This space of functions will be denoted in the sequel.
The second hypothesis is more difficult to take into account. We chose in this work to tackle the problem via an energy minimization method where the energy functional is composed of two terms: the first one is a matching term , the second one is a penalization term which tries to handle the second assumption. The matching term will be a -norm which is justified by the Gaussian white noise.
In the case where is not linear, the exact method to remove the exterior is to operate the blur function on the addition of the known exterior and the center. For the sake of simplicity, we will not use this method and will consider that the errors are negligible when subtracting the projections.
Let us first describe more precisely the set . This time, the functions will be defined on , with values still in , with compact support. Therefore, such a function will be characterized by the knowledge of the curves that limit the two areas where is equal to and to . Indeed, as the support of the function is bounded, these curves are disjoint Jordan curves and the density of the inside is whereas the density of the outside is 0. Consequently, the energy that we will consider will be a function of where is a set of disjoint Jordan curves. For mathematical reasons, we must add an extra-assumption : the curves are
so that the normal vector of the curves is well-defined (as an orthogonal vector to the tangent one).
In this continuous framework, the matching term is just the usual -norm between and the data (where is given by (3.1)). So, the first term is
For the penalization term, we choose
where denotes the length of the curves . Let us remark that this penalization term may be also viewed as the total variation (up to a multiplicative constant) of the function
because of the binarity. Eventually, the total energy functional is
which is an adaptation of the well-known Mumford-Shah energy functional introduced in . The “optimal” value of may depend on the data.
3.3. A continuous model in BV space
The previous analysis gives the mains ideas for the modelization. Now, we make it precise using an appropriate functional framework. Let be a bounded open subset with Lipschitz boundary. We shall consider bounded variation functions. Recall that the space of such functions is
where denotes the space of functions with compact support in . The space endowed with the norm
is a Banach space.
If its derivative in (distributions) is a bounded Radon measure denoted and is the total variation of on . Let us recall useful properties of BV-functions ():
Let be an open subset with Lipschitz boundary.
1. If , we get the following decomposition for :
where is the absolutely continuous part of with respect of the Lebesgue measure and is the singular part.
2. The map from to is lower semi-continuous (lsc) for the topology.
3. with compact embedding.
4. with compact embedding.
We precise hereafter an important continuity property of the projection operator .
The projection operator is continuous from to for every and .
Proof - It is a direct consequence of Hölder inequality. Let be and . Its support is included in which is included in some where and only depends on . It is clear that is defined everywhere on and, for every
where . Therefore
The computations in the case are similare and lead to the same inequality (with some additional absolute values). As
here and in the sequel denotes a generic constant depending on and . So
As is bounded, this yields
As we consider the length of curves, the most suitable functional space to set a variational formulation of the reconstruction problem is . Therefore, we consider the following minimization problem
stands for the - norm, and
The operator is given by (3.3). Without loss of generality, we may assume (for simplicity) that
At last, “”, is the binarity constraint. We have mentionned that the image takes its values in where . With the change of variable , we may assume that the image values belong to .
A similar problem has been studied in  with smoother projection operator and convex constraints. This is not our case. The pointwise constraint “ ” is a very hard constraint. The constraint set is not convex and its interior is empty for most usual topologies.
Now we may give the main result of this section :
Problem admits at least a solution.
Proof - Let be a minimizing sequence. It satisfies ; so
Therefore the sequence is - bounded. As is bounded as well, the sequence is bounded in . Thus it converges (extracting a subsequence) to some for the weak-star topology.
Estimate (3.8) implies the weak convergence of to in for every . Thanks to the continuity property of proposition 3.4, we assert that weakly converges to in . We get
with the lower semi-continuity of the norm.
Moreover is compactly embedded in . This yields that strongly converges to in . As is lsc with respect to - topology, we get
As the pointwise constraint is obviously satisfied, is a solution to .
4. Computation of the energy derivative
Now we look for optimality conditions. Unfortunately we cannot compute easily the derivative of the energy in the framework. Indeed we need regular curves and we do not know if the minimizer provides a curve with the required regularity. Moreover, the set of constraints is not convex and it is not easy to compute the Gâteaux- derivative (no admissible test functions).
So we have few hope to get classical optimality conditions and we rather compute minimizing sequences. We focus on particular ones that are given via the gradient descent method inspired by . Formally, we look for a family of curves such that
so that decreases as . Let us compute the energy variation when we operate a small deformation on the curves . In other word, we will compute the Gâteau derivative of the energy for a small deformation :
We will first focus on local deformations . Let be a point of . We consider a local reference system which center is and axis are given by the tangent and normal vectors at and we denote by the new generic coordinates in this reference system. With an abuse of notation, we still denote . We apply the implicit functions theorem to parametrize our curve: there exist a neighborhood of and a function such that, for every ,
Eventually, we get a neighborhood of , a neighborhood of and a function such that
The local parametrization is oriented along the outward normal to the curve at point (see figure 6). More precisely, we define the local coordinate system where is the usual tangent vector,
is the direct orthonormal vector; we set the curve orientation so thatis the outward normal. The function if then defined on by
This parametrization is described on figure 6. We then consider a local (limited to ) deformation along the normal vector. This is equivalent to handling a function whose support is included in . The new curve obtained after the deformation is then parametrized by
This defines a new function :
We will also set . This deformation is described on figure 6.
The Gâteau derivative for the energy has already been computed in  and is
where curv denotes the curvature of the curve and is the parametrization of .
It remains to compute the derivative for the matching term. First we estimate : a simple computation shows that
Now we compute where (resp. ) is the curve associated to the function (resp. ):
To simplify the notations, we denote by so that we need to compute
As is zero out of the neighbourhood , we have
In the case , we have,
As the function is continuous (and thus bounded on ), we may pass to the limit by dominated convergence and get
In the case , we have
and we obtain the same limit as in the nonnegative case.
Finally, the energy derivative is
If we set , we get
As formula (4.12) may be written
where denotes the outward pointing normal unit vector of the curve , denotes the usual scalar product in and is a positive coefficient that depends on the curvilinear abscissa .
The latter expression is linear and continuous in , this formula is also true for a non-local deformation (which can be achieved by summing local deformations).
5. Front propagation and level set method
The goal of the present section is to consider a family of curves that will converge toward a local minimum of the functional energy. From equation (4.13), it is clear that if the curves evolve according to the differential equation
the total energy will decrease.
To implement a numerical scheme that discretizes equation (5.14), it is easier to use a level set method (see  for a complete exposition of the level set method). Indeed, equation (5.14) may present some instabilities, in particular when two curves collide during the evolution or when a curve must disappear. All these evolutions are handled easily via the level set method.
The level set method consists in viewing the curves as the 0-level set of a smooth real function defined on . The function that we are seeking is then just given by the formula
We must then write an evolution PDE for the functions that corresponds to the curves . Let be a point of the curve and let us follow that point during the evolution. We know that this point evolves according to equation 5.14
We can re-write this equation in terms of the function recognizing that
where stands for the gradient of with respect to , denotes the euclidean norm. The evolution equation becomes
Then, as the point remains on the curve , it satisfies . By differentiating this expression, we obtain
which leads to the following evolution equation for :
The above equation is an Hamilton-Jacobi equation which involves a non local term (through and ). Such equations are difficult to handle especially when it is not monotone (which is the case here). In particular, existence and/or uniqueness of solutions (even in the viscosity sense) are not clear. The approximation process is not easy as well and the numerical realization remains a challenge though this equation is a scalar equation which is easier to discretize than the vectorial one. Here we used the discrete scheme described in  to get numerical results.
6. Results and discussion
6.1. Explicit scheme
We briefly present the numerical scheme. We used an explicit scheme in time and the spatial discretization has been performed following  . We set
so that the equation (5.15) is
We set with et . The explicit Euler scheme gives :
The curvature term is computed as
where stands for the partial derivative with respect to . The discrete approximation of the gradient is standard:
are defined in the same way. A usual approximation for is given by :
The non local term is exactly computed.
6.2. Numerical results
The previous scheme has been implemented on a 3.6 GHz PC. A classical reinitialization process has been used each 500 iterations. The test image size was 256 256 pixels. The other parameters of the computation were set to
and the blur kernel is a Gaussian kernel of standard deviation 5 pixels.
The computed image is quite satisfying (see figure 7.) However, we note a bad reconstruction along the symmetry axis due to the problem geometry and a lack of information. Moreover we have to improve the algorithme behavior. Indeed, we observe numerical instability (in spite of the regularization process) that leads to a very small time step choice. Therefore the computational time is quite long (about 2.5 hours). In addition, classical stopping criteria are not useful here : the expected solution corresponds to a “flat” level of function and the difference between two consecutive iterates means no sense. An estimate of the cost function decrease is not appropriate as well (we observe oscillations). We decided to stop after a large enough number of iterations (here 20 000).
In spite of all these disadvantages, this method is satisfactory considering the low signal to noise ratio of the radiograph. These good results can be explained by the strong assumptions that we add (in particular the binary hypothesis) which are verified by our synthetic object. Anyway, the method has beeen successfully tested on “real” images as well, that is imgaes of objets with the same kind of properties (“almost” binary) but we cannot report them here (confidential data). A semi-implicit version of the algorithm is actually tested to improve stability.
-  I. Abraham, R. Abraham Technical Report CEA (2001).
-  L. Ambrosio, N. Fusco et D. Pallara,Functions of bounded variation and free discontinuity problems, Oxford mathematical monographs, Oxford University Press, 2000.
J.P. Bruandet, F. Peyrin, J.M. Dinten, M. Barlaud 3D tomographic reconstruction of binary images from cone-beam
projections: a fast level-set approach
2002 IEEE International Symposium on Biomedical Imaging, p. 677-80, (2002)
-  E. Casas, K. Kunisch and C. Pola, Regularization by Functions of Bounded Variation and Applications to Image Enhancement, Applied Mathematics and Optimization, 40:229Ð257 (1999)
J.-M. Dinten Tomographie à partir d’un nombre limité de projections
: régularisation par champs markovien
PHD thesis, Université d’Orsay Paris-Sud (1990)
N. J. Dusaussoy Image reconstruction from
SPIE’s international symposium on optics, imaging and instrumentation. San Diego (1994)
H. Feng, W. Karl, D. Castanon A curve
evolution approach to object-based tomographic reconstruction
IEEE Trans. on Image Proc., 12, 44-57, (2003)
K. Hanson Tomographic reconstruction of
axially symmetric objects from a single radiograph
High Speed Photography, 491, (1984)
G. Herman Image reconstruction from projections: the fundamentals of
Academic Press (1980)
J.M. Lagrange Reconstruction tomographique à partir d’un petit nombre de vues
PHD thesis, ENS Cachan (1998)
D. Mumford - J. Shah Optimal approximations by
piecewise smooth functions and associated variational problem
Comm. Pure and Appl. Math. 42, 577-685 (1989)
J.A. Sethian Theory, algorithm and
applications of level set method for propagating interfaces
Iserles, A. (ed.), Acta Numerica Vol. 5, 1996. Cambridge: Cambridge University Press. 309-395 (1996)
-  E. Stein Singular integrals and differentiability properties of functions, Princeton University Press