# Tomographic Reconstruction of Homogeneous 2D Geometric Models with Unknown Attenuation

Zenith Purisha, Samuli Siltanen

► To cite this version:

Zenith Purisha, Samuli Siltanen. Tomographic Reconstruction of Homogeneous 2D Geometric Models with Unknown Attenuation. 26th Conference on System Modeling and Optimization (CSMO), Sep 2013, Klagenfurt, Austria. pp.247-256, 10.1007/978-3-662-45504-3_24 . hal-01286431

HAL Id: hal-01286431

https://inria.hal.science/hal-01286431

Submitted on 10 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Tomographic reconstruction of homogeneous 2d
geometric models with unknown attenuation*

Zenith Purisha and Samuli Siltanen

Department of Mathematics and Statistics, University of Helsinki,
Finland
{zenith.purisha,samuli.siltanen}@helsinki.fi

**Abstract.** A new method is presented for tomographic reconstruction of objects with homogeneous attenuation. The method is based on parametric representation with Non-Uniform Rational B-Splines (NURBS) and statistical inversion with a Markov Chain Monte Carlo (MCMC) algorithm. The method recovers the approximate boundary curve shape and the attenuation value of two-dimensional homogeneous objects. The boundary can be represented by NURBS with few parameters, reducing the number of degrees of freedom. However, this leads to a nonlinear inverse problem, and therefore statistical inversion is used. One of the benefits of the approach is that the reconstruction is automatically in the form of the geometrical representation in industrial CAD format or CNC configuration. Computational results are presented with two different simulated homogeneous geometric models and sparsely sampled tomographic data. The new method outperforms the baseline method (filtered back-projection) in image quality but not in computational speed.

**Keywords:** tomography; homogeneous; CAD, NURBS; Bayesian inversion; MCMC

# 1 Introduction

Creating a virtual model of a given physical object is increasingly important in, for example, reverse engineering and game development. The details of reconstructing the model depend on the kind of measurements that are available about the object. For example laser scanning and digital photography are popular methods providing surface information. In this work we concentrate on sparsely sampled X-ray tomography measurements.

Consider a three-dimensional cylindrical object $\Omega \times \mathbb{R}$ with a simply connected base $\Omega \subset \mathbb{R}^2$. Furthermore, assume that we only know that the object is homogeneous: the X-ray attenuation coefficient has an unknown but constant value $c > 0$ inside the object.

<footnote>* This work was supported by the Academy of Finland through the Finnish Centre of Excellence in Inverse Problems Research 2012–2017, decision number 250215.</footnote>
We discuss the situation where we have X-ray projection data of a transversal slice of the object. In other words, we have access to a collection of line integrals of the function $f: \mathbb{R}^2 \to \mathbb{R}$ defined by

$$f(x, y) = \begin{cases} c & \text{for } (x, y) \in \Omega, \\ 0 & \text{for } (x, y) \in \mathbb{R}^2 \setminus \Omega. \end{cases} \qquad (1)$$

The angular sampling of the X-ray data can be very sparse, allowing for quick measurement process with low radiation dose. Our aim is to recover two things: the boundary $\partial\Omega \subset \mathbb{R}^2$ represented as a parameterized curve and the attenuation coefficient $c$.

We want our method to be practically useful in industrial environments. Computer numerical control (CNC) machines are widely used in modern production facilities, and they use computer-aided design (CAD) models. The proposed tomographic method represents the unknown boundary curve in Non-Uniform Rational Basis Spline (NURBS) form, which is the standard in CAD software. This direct connection with industrial standards is the main motivation behind the proposed method.

NURBS curves are represented by a relatively small number of parameters: a set of planar control points and a related knot vector. In this paper we fix the knot vector, so the information to be recovered consists merely of the control points and the attenuation parameter. The low dimensionality of this problem formulation offers computational advantages. However, there is a complication as well: the linear inverse problem of X-ray tomography becomes nonlinear in this parameterization. Therefore, we resort to the very general framework of Bayesian inversion [4,10].

In Bayesian inversion, limited measurement data is complemented by a priori information using the Bayes formula. This way the ill-posed inverse problem is recast in a well-posed form of exploring the posterior probability distribution. As explained in [5], in the case of X-ray tomography this involves a discrete attenuation model and a Monte Carlo Markov Chain (MCMC) method for sampling the posterior. Usually the large number of pixels in the reconstructed image leads to MCMC sampling in a very high-dimensional space (one dimension for each pixel). In our case the posterior distribution is defined in a relatively low-dimensional space: one dimension for the attenuation value plus two dimensions for each control point. This enables efficient MCMC sampling.

The a priori information we use is rather simple: we assume that we have an upper bound for the diameter of the two-dimensional shape (transversal slice) under measurement. Also, we assume that the curve does not have too small details (parts with very high curvature) and choose the number of control points to be as small as possible while still capable of representing the smallest details in the curve.

We demonstrate the novel NURBS-MCMC method using two simulated non-convex examples. See Figure 3 below. The reconstruction algorithm is found to recover the attenuation coefficient quite precisely and the boundary shape with reasonable accuracy from very sparsely sampled X-ray data (only 18 projection directions).
This paper is organized as follows. In section 2, we discuss the theory of NURBS curves. In section 3 we present the X-ray measurement model. Section 4 is devoted to the description of Bayesian inversion. In section 5, we present the reconstruction results, and in section 6 we conclude our findings.

## 2 NURBS description for parametric curve

We model an unknown object boundary $\partial\Omega$ by a continuous curve $S: [0, 1] \to \mathbb{R}^2$. In our computational problem, we construct $S$ using NURBS that are widely used as computationally fast and robust representations of curves.

The basic building blocks of NURBS are the following:

1. Control points $\mathbf{p}_1, \dots, \mathbf{p}_n$. These planar locations $\mathbf{p}_i \in \mathbb{R}^2$ are, roughly speaking, points of attraction for the NURBS curve, where $i = 1, 2, \dots, n$. Throughout the paper we denote by $n$ the number of control points.

2. Knots $t_1, t_2, \dots, t_K \in [0, 1]$, with ordered as follows:

$$0 = t_1 \le t_2 \le \dots \le t_K = 1,$$

where $K > n$. The knot are used to divide the interval $[0, 1]$ into suitable pieces. We collect the knots into a knot vector $[t_1 \ t_2 \ \dots \ t_K]$.

3. Basis function ($N_{i,p}(t)$) specifies how strongly the control point $\mathbf{p}_i$ attracts the NURBS curve. The first-order basis function is

$$N_{i,1}(t) = \begin{cases} 1 & \text{if } t_i \le t < t_{i+1}, \\ 0 & \text{otherwise.} \end{cases}$$

Higher-order basis functions are defined recursively as

$$N_{i,p}(t) = \frac{t - t_i}{t_{i+p} - t_i} N_{i,p-1}(t) + \frac{t_{i+p+1} - t}{t_{i+p+1} - t_{i+1,p-1}} N_{1+i,p-1}(t),$$

where $p$ is the order of the basis function and $K = n + p$.

The general form of NURBS curve can be written as

$$S(t) = \sum_{i=1}^{n} \mathbf{p}_i R_{i,p}(t). \qquad (2)$$

The Rational in NURBS comes from the rational function $R_{i,p}(t) = \frac{\omega_i N_{i,p}(t)}{\sum_{i=0}^{n} \omega_i N_{i,p}(t)}$, where the weights $\omega_i \ge 0$, for all $i$. In this preliminary result, we use the same weights for all control points.

## 3 Tomographic measurement model

Consider a continuous tomography model $f: \mathbb{R}^2 \to \mathbb{R}$ as in (1), where $f(x, y) \ge 0$ and $\text{supp}(f) \subset \Omega$ with bounded $\Omega \subset \mathbb{R}^2$.
<figure><img src="image_2.png" /><figcaption>Fig. 1. An X-ray travels along a homogeneous target slice. The shade of gray describes a constant coefficient attenuation c inside Ω.</figcaption></figure>

Consider an X-ray traveling through a two dimensional object along a straight line as shown in the Fig. 1. In this specimen, the slice of the target is in square defined by $0 \le x \le a$ and $0 \le y \le b$. Assume that an X-ray penetrates along the horizontal path $0 \le x \le a$ and $y = b_1$.

Let us consider that the X-ray has the initial intensity $I_0 = I(0)$ and the intensity becomes smaller, say $I_1 = I(1)$ after it passes the object. This situation can be modeled using $f(x, y)$, an attenuation coefficient function, as:

$$ \frac{dI(x)}{I(x)} = -f(x, b_1)dx, $$

where $I(x)$ is the intensity of the X-ray at the point $(x, b_1)$ while passing through the source to the detector.

In tomographic imaging, we want to collect information about $f$ using different angles. Let us consider the Radon transform, denoted by $\mathcal{R}$, as follows. Assume $\alpha \in \mathbb{R}$ as an angle measured in radians:

$$ \alpha = \begin{bmatrix} \cos \alpha \\ \sin \alpha \end{bmatrix} \in \mathbb{R}^2, $$

the unit vector with angle $\alpha$ with respect to the x-axis.

The radon function of the function $f$ depends on the angular parameter $\alpha$ and on a linear parameter $s \in \mathbb{R}$ as follows:

$$ \mathcal{R}f(s, \alpha) = \int_{\mathbf{x} \cdot \alpha = s} f(\mathbf{x})d\mathbf{x}^{\perp}, $$

where $d\mathbf{x}^{\perp}$ is the one dimensional Lebesgue measure along the line $\{\mathbf{x} \in \mathbb{R}^2 : \mathbf{x} \cdot \alpha = s\}$.

For computational reasons, we need a discrete model. In this case, we construct two discrete models: a pixel-based object model and a NURBS-based object model, a model where the boundary $\partial\Omega$ is expressed as a NURBS curve as shown in Fig. 2.

In the pixel-based model, the line integral is discretized using the standard pencil-beam model. We use the pixel-based Matlab routine `radon.m` for simulating parallel-beam tomographic data.
In the NURBS-based model, the line integral is discretized by moving to pixel-based model using an operator $\mathcal{B}$ defined by

$$ \mathcal{B}(\mathbf{p}, c) = \begin{cases} c, & \text{if the pixel center is inside the NURBS curve,} \\ 0, & \text{if the pixel center is outside the NURBS curve.} \end{cases} \qquad (3) $$

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Left: the NURBS-based object model where $\partial\Omega$ is a NURBS curve. The inside of the curve is set to be $c$ and the outside is set to 0. Right: the pixel image.</figcaption></figure>

Assuming that the knot vector is fixed, the degrees of freedom in our NURBS model are the control points $\mathbf{p}_1, \dots, \mathbf{p}_n$ together with the attenuation.

In the simulation, we measure two simple homogeneous shapes that have different attenuation. To avoid inverse crime [7], we produce the synthetic phantoms $\Omega_1$ and $\Omega_2$ without using NURBS. Those objects are set to be homogeneous inside with attenuation values 2 and 3,5, respectively, as shown in Fig. 3.

The objects are measured with the resolution 64 × 64 using parallel beam geometry as shown in Fig. 4. From the source, the X-ray penetrates through the objects and a sensor detects the projection images from different directions. Sparse full angle data, 0°, 10°, 20°, ..., 170°, are applied to obtain the projections and each direction consists of 95 lines.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Homogeneous phantoms</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Parallel beam X-ray measurement geometry. There are 6 different directions (0°, 30°, 60°, 90°, 120°, 150°) and 15 lines. Black dots show the locations of the X-ray source at different times of measurement. The thick line represents the detector measuring the intensity of the X-rays after passing through the target.</figcaption></figure>

## 4 Bayesian inversion for control points and attenuation value

This section presents the Bayesian approach to handle the inverse problem. This measurement data, **m**, is used to get information about other quantities. In this case, we encounter a nonlinear inverse problem, which need to be solved by recovering **B** that depends on **p** and **c**.

We model the problem as the following form:

$$ \mathbf{m} = \mathcal{R}(\mathcal{B}(\mathbf{p}, \mathbf{c})) + \varepsilon, \quad (4) $$

where $\varepsilon$ is the error of the measurement.

The Bayesian inversion approach is based on the relations between probability distributions to model the inadequacy of information in an inverse problem. Before performing the collection of measurement data, we construct a model for a priori knowledge. Since the control points are presented in polar coordinates, i.e. $\mathbf{p}_i = (r_i \sin \theta_i, r_i \cos \theta_i)$, we assume that the angle of each parameter is not less than $\theta_i^{\min}$ and not more than $\theta_i^{\max}$, and the distance of each parameter from the central point of the object is nonnegative and not more than $r_i^{\max}$. In this case, $\Omega_1$ and $\Omega_2$ have $r_i^{\max}$ values that equal to 15 and 30, respectively, and the maximum of the attenuation value $c_{\max}$ is 5 for both objects.
We formulate the prior condition as follows

$$ \pi(\mathbf{m} | (\mathbf{p}, c)) = \begin{cases} \exp(-\frac{1}{2\sigma_2^2} \|(\mathbf{p}, c) - (\tilde{\mathbf{p}}, \tilde{c})\|_2^2) & \text{for } 0 \le r_i \le r_i^{\text{max}} \text{ and } 0.1 < c < c_{\text{max}} \\ 0 & \text{and } \theta_i^{\text{min}} \le \theta_i \le \theta_i^{\text{max}}, \\ 0 & \text{otherwise,} \end{cases} \tag{5} $$

where $\mathbf{p} = \{\mathbf{p}_i\}$ and $(\tilde{\mathbf{p}}, \tilde{c})$ is a priori information of the position of control points and the attenuation value. After examining the measurement setting and prior information, we can model the conditional probability of $\mathbf{m}$, which is called the likelihood function

$$ \pi(\mathbf{m} | (\mathbf{p}, c)) = \exp\left(-\frac{1}{2\sigma_1^2} \|\mathcal{R}(\mathcal{B}(\mathbf{p}, c)) - \mathbf{m}\|_2^2\right). \tag{6} $$

By given observed data, $\mathbf{m}$, the conditional probability $\pi(\mathbf{p}, c | \mathbf{m})$ of $\mathbf{p}$ and $c$ can be expressed as follows

$$ \pi(\mathbf{p}, c | \mathbf{m}) = \frac{\pi(\mathbf{p}, c)\pi(\mathbf{m} | (\mathbf{p}, c))}{\pi(\mathbf{m})}, \tag{7} $$

which is called the posterior distribution. To solve our inverse problem, we need to explore this distribution.

As a common method to represent statistical estimates, we apply the conditional mean (CM) of the unknown $\mathbf{p}$ and $c$. Since CM is defined as

$$ (\mathbf{p}^{\text{CM}}, c^{\text{CM}}) = \int_{\mathbb{R}^N} (\mathbf{p}, c)\pi(\mathbf{p}, c | \mathbf{m}) d(\mathbf{p}, c), $$

finding the estimate leads to the integration problems. Typically, the integration is over a high-dimensional space. To unfold this issue, a Markov chain Monte Carlo (MCMC) technique is recommended to generate a sample from the posterior distribution. For a general introduction to Bayesian inversion and properties of MCMC computation see [2,3,4].

By applying the CM estimate to the samples $\{\mathbf{p}_{1,l}, \mathbf{p}_{2,l}, \dots, \mathbf{p}_{n,l}, c_l\}$, we get

$$ \mathbf{p}_i^{\text{CM}} \approx \frac{1}{N} \sum_{l=1}^{N} \mathbf{p}_{i,l}^{\text{CM}} \text{ and } c_i^{\text{CM}} \approx \frac{1}{N} \sum_{l=1}^{N} c_l^{\text{CM}}, $$

where $\mathbf{p}^{\text{CM}} = \{\mathbf{p}_i^{\text{CM}}\}$, $i = 1, 2, \dots, n$ and $N$ is the number of evaluations. For the NURBS curve reconstruction with $N$ evaluations, it is written as $S_N^{\text{CM}}$.

# 5 Computational results

In this section, numerical examples are presented. We use Metropolis-Hastings as sampling algorithm to generate control points and attenuation with 1 000 000 iterations (applied also to the Radon transform and its adjoint). In each
<figure><img src="image_6.png" /><figcaption>Fig. 5. The thin black line is the target curve. The thick black line is the reconstruction of the NURBS curve, $S_N^{CM}$. Top: reconstructions for $\Omega_1$, bottom: reconstructions for $\Omega_2$. The black circle markers are the control points, $\mathbf{p}^{CM}$.</figcaption></figure>

iteration, the weights are set to be equal while the order and knot vector of NURBS curve are set to be fixed. The order is set to be 3 because it is widely used in practical application and to avoid heavy calculation times. As a default knot vector in CAD, the open uniform knot vector is chosen for $\Omega_1$ and $\Omega_2$, $[0 \ 0 \ 0 \ \frac{1}{11} \ \frac{2}{11} \ \frac{3}{11} \ \frac{4}{11} \ \frac{5}{11} \ \frac{6}{11} \ \frac{7}{11} \ \frac{8}{11} \ \frac{9}{11} \ \frac{10}{11} \ 1 \ 1 \ 1]$ and $[0 \ 0 \ 0 \ \frac{1}{7} \ \frac{2}{7} \ \frac{3}{7} \ \frac{4}{7} \ \frac{5}{7} \ \frac{6}{7} \ 1 \ 1 \ 1]$, respectively.

The NURBS curves as in the rightmost Fig. 5 are achieved. By using the mapping as in (3), both final shape reconstructions are presented in the middle of Fig. 6 and Fig. 7. The error in the shape reconstructions is given as follows. Denote $O$ as the image of the original 2D object and $O^{rec}$ as the image of the reconstruction. Set $O\setminus O^{rec}$ for points that belong to the original object but not to the reconstruction and $O^{rec}\setminus O$ for points that belong to the reconstruction but not to the original object. The relative error in the reconstruction is written as

$$ \frac{(\text{area}(O^{\text{rec}}) + \text{area}(O^{\text{rec}}\setminus O))}{\text{area} O} \times 100\% \qquad (8) $$

By applying (8), the relative errors of $\Omega_1$ and $\Omega_2$ reconstructions using NURBS-MCMC are 15% and 8.3%, respectively. Recovered chains of attenuation values of $\Omega_1$ and $\Omega_2$ after burn-in period have relative errors 9.26% and 1.47%, respectively.

The rightmost images in Fig. 6 and Fig. 7 show recovered shapes using filtered back projection (FBP). The reconstruction uses the resolution 64 × 64. To assess the error in the reconstructed attenuation value, a representative rectangular
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Left: Original $\Omega_1$. Center: NURBS-MCMC reconstruction. Right: FBP reconstruction. Both are using error 0.1%.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Left: Original $\Omega_2$. Center: NURBS-MCMC reconstruction. Right: FBP reconstruction. Both are using error 0.1%.</figcaption></figure>

region of interest is picked from the inside the reconstruction. The mean and the standard deviation of the recovered attenuation values are computed as we can see in Table 1, while Table 2 shows computation times for both methods.

## Discussion and conclusions

Reconstruction using the NURBS-MCMC method in nonlinear inverse problem can recover measurement data successfully. Homogeneous objects $\Omega_1$ and $\Omega_2$ are recovered by only $2n + 1$ parameters: 25 and 17, respectively. Those recovered data are geometrical representations which are automatically set to CAD or CNC configuration. In the middle of Fig. 6 and Fig. 7, the vector graphic form is converted to be 512 $\times$ 512.

<table><caption>**Table 1.** Mean and standard deviation of FBP reconstruction.</caption><thead><tr><td></td><td>mean</td><td>standard deviation</td></tr></thead><tbody><tr><td>Ω<sub>1</sub></td><td>1.9795</td><td>0.11</td></tr><tr><td>Ω<sub>2</sub></td><td>3.49</td><td>0.14</td></tr></tbody></table>
<table><caption>**Table 2.** Computation time (in seconds) for all reconstruction methods.</caption><thead><tr><td>FBP</td><td>NURBS-MCMC</td></tr></thead><tbody><tr><td>1</td><td>18 000</td></tr></tbody></table>

In filtered backprojection, the reconstruction is represented by pixel images and consequently doing a segmentation to represent the shape is nontrivial.

Nevertheless, the slowness of computation is a shortcoming of the proposed method as we can see in Table 2, but by implementing parallel computing, the problem can be handled.

## References

1. Bertrand, C., et al.: A probabilistic solution to the MEG inverse problem via MCMC methods: the reversible jump and parallel tempering algorithms. Biomedical Engineering, IEEE Transactions on 48.5, 533-542 (2001)

2. Gamerman, D., and Lopes, H.F.,: Markov Chain Monte Carlo In Practice. Chapman and Hall/CRC (1996)

3. Gilks, W. R., et al.; Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference. Chapman and Hall/CRC (2006)

4. Kaipio, J., Somersalo, E.: Statistical and Computational Inverse Problems. Springer (2004)

5. Kolehmainen, V. and Siltanen, S., et al.: Statistical inversion for medical x-ray tomography with few radiographs: II. Application to dental radiology. Phys. med. Biol. 48,1465-1490 (2003)

6. Marzouk, Y. M. and Habib N. N.: Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems. Journal of Computational Physics 228.6, 1862-1902 (2009)

7. Mueller J. and Siltanen S.: Linear and Nonlinear Inverse Problems with Practical Applications. SIAM, Computational Science and Engineering (2012).

8. Renken, F., and Subbaraya, G.: Nurbs-based solutions to inverse problems in droplet shape prediction, Computer methods in applied mechanics and engineering, 190, pp.1391-1406 (2000)

9. Rogers, D. F.: An Introduction to NURBS: with historical perspective, vol. 1 of Academic Press, Morgan Kaufmann (2001)

10. Tarantola, A.: Inverse problem theory and methods for model parameter estimation. Society for Industrial Mathematics (2005)

11. Wang, J., and Zabaras, N.: Hierarchical Bayesian models for inverse problems in heat conduction. Inverse Problems 21.1, 183 (2004)