# Image Reconstruction With Local Directional Total Variation Regularization Using Tomographic Incompleteness Maps

Matthieu Laurendeau, Laurent Desbat, Frédéric Jolivet, Guillaume Bernard, Sébastien Georges, Simon Rit

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

To cite this version:

Matthieu Laurendeau, Laurent Desbat, Frédéric Jolivet, Guillaume Bernard, Sébastien Georges, et al. Image Reconstruction With Local Directional Total Variation Regularization Using Tomographic Incompleteness Maps. The 8th International Conference on Image Formation in X-Ray Computed Tomography, Marc Kachelrieß, Aug 2024, Bamberg, Germany. pp.376-379. hal-04704936

HAL Id: hal-04704936

https://hal.science/hal-04704936v1

Submitted on 21 Sep 2024

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.

Public Domain
# Image Reconstruction With Local Directional Total Variation Regularization Using Tomographic Incompleteness Maps

Matthieu Laurendeau, Frédéric Jolivet, Laurent Desbat, Guillaume Bernard, Sébastien Gorges and Simon Rit

***Abstract***—Limited angle acquisition is a well-known challenge in computed tomography reconstruction because the lack of data generates severe geometric artifacts in reconstructed images. Recently, directional total variation (DTV) regularization has shown promising results for this kind of problem but it requires fine-tuning of a global directional hyperparameter in addition to the regularization weight. In this work, we propose a new regularization, called local directional total variation (LDTV), which is a DTV based on local directional weights determined from tomographic incompleteness maps. We evaluated LDTV and compared it to state-of-the-art algorithms on simulated two-dimensional acquisitions of the Forbild head phantom with a source trajectory made of two orthogonal arcs of 60° each. The reconstructed images show that the LDTV regularization performs better in this geometry for both noiseless and noisy data.

***Index Terms***—limited angle problem, tomographic incompleteness, directional total variation (DTV), local DTV (LDTV)

## I. INTRODUCTION

IN X-ray computed tomography (CT), some applications require limited space which implies a reduced acquisition angle. This situation could arise both in the medical field, with tomosynthesis for example, and in the industrial field where the aim is to minimize the size of the system and reduce its scanning time.

Several image reconstruction algorithms have been developed to address the limited angle problem. Iterative methods have demonstrated improved results compared to filtered-backprojection (FBP) algorithms [1]. They minimize a cost function made of a data fidelity term and a regularization term. The latter implements an image prior compensating the lack of data to obtain more plausible results.

M. Laurendeau is with Univ Lyon, INSA-Lyon, Université Claude Bernard Lyon 1, UJM-Saint Étienne, CNRS, Inserm, CREATIS UMR 5220, U1294, F-69373, Lyon, France, with Univ. Grenoble Alpes, CNRS, UMR 5525, VetAgro Sup, Grenoble INP, TIMC, 38000 Grenoble, France and with Thales AVS, Moirans, France (email: matthieu.laurendeau@insa-lyon.fr).

F. Jolivet is with Thales AVS, Moirans, France (email: frederic.jolivet@thalesgroup.com).

L. Desbat is with Univ. Grenoble Alpes, CNRS, UMR 5525, VetAgro Sup, Grenoble INP, TIMC, 38000 Grenoble, France (email: laurent.desbat@univ-grenoble-alpes.fr).

G. Bernard is with Thales AVS, Moirans, France (email: guillaume.bernard@thalesgroup.com).

S. Gorges is with Thales AVS, Moirans, France (email: sebastien.gorges@thalesgroup.com).

S. Rit is with Univ Lyon, INSA-Lyon, Université Claude Bernard Lyon 1, UJM-Saint Étienne, CNRS, Inserm, CREATIS UMR 5220, U1294, F-69373, Lyon, France (email: simon.rit@creatis.insa-lyon.fr).

Recently, a new regularization, directional total variation (DTV), was developed [2] and showed promising results for the limited angle problem. This regularization strengthens the image prior in the direction where data is missing. In [3], the authors studied orthogonal arcs with limited angular range for dual-energy CT application and showed similarly good results using DTV.

Although DTV could improve visual image quality, it requires fine-tuning of the directional parameter, which could be time-consuming and sensitive to the acquisition geometry and the scanned object. Moreover, the regularization applies the same directional weight for all pixels of the image, which may not accurately describe spatial variability of data incompleteness due to limited angle acquisitions.

In this work, we present the local directional total variation (LDTV) regularization based on the tomographic incompleteness map previously developed in [4] which uses the incompleteness metric of [5]. This regularization locally adjusts both the strength and the directional weight parameters by taking into account the geometry of the scanner. We simulate noiseless and noisy projections of the Forbild head phantom [6] with a limited angle geometry made of double orthogonal arcs (DOA), compute the corresponding tomographic incompleteness map and compare several regularizations, including state-of-the-art DTV [1] and our LDTV method.

## II. MATERIALS AND METHODS

### A. Geometry and tomographic incompleteness

Similarly to [3], we used a two-dimensional (2D) fan-beam geometry with limited angle DOA. An illustration of the geometry is provided in Fig. 1. We used 120 source positions, evenly distributed every degree on two orthogonal arcs of 60° each. The source-to-object distance (SOD) and the source-to-detector distance (SDD) were 510 mm and 1020 mm, respectively. In addition to having an acquisition limited to an angular range of 120°, this geometry was selected because the direction of missing data is spatially variable with both horizontal and vertical directions.

Monochromatic fan-beam projections were simulated from a reference image of the Forbild head phantom [6] in the region $\Omega \subset \mathbb{R}^2$ with 256 x 256 pixels and 1 mm spacing (Fig. 1). The linear detector was large enough (1200 mm, 1 mm spacing) to avoid lateral truncation of the projections.

We computed the tomographic incompleteness map of the scanner geometry, as described in [4], in the source trajectory
<figure><img src="image_2.png" /><figcaption>Fig. 1. Double orthogonal arcs (DOA) geometry. The two 60° blue dashed arcs represents the 120 sources distributed every degree. The image of the Forbild head phantom is placed at the center of the arcs.</figcaption></figure>

plane. The tomographic incompleteness criterion [5], noted $I(\mathbf{p}, \mathbf{n}) \in \mathbb{R}^+$, measures how much a 3D point fulfills Tuy's condition, i.e., for a continuous source trajectory $\Gamma$ and without truncations, a point $\mathbf{p}$ can be reconstructed if all planes passing through $\mathbf{p}$ intersect $\Gamma$. One advantage of the incompleteness criterion $I$ is that it allows evaluation of the tomographic incompleteness for a discrete source trajectory. The value of $I(\mathbf{p}, \mathbf{n})$ is defined as the tangent of the minimal angle between the plane $\Pi_{\mathbf{p},n}$ passing through $\mathbf{p}$ with co-direction $\mathbf{n}$ and all the X-ray lines defined by the source points and $\mathbf{p}$ [5]. When $I(\mathbf{p}, \mathbf{n})$ reaches 0, $\Pi_{\mathbf{p},n}$ intersects the source trajectory and larger values of $I(\mathbf{p}, \mathbf{n})$ indicate a larger angle between $\Pi_{\mathbf{p},n}$ and the closest source position. As in [4], the tomographic incompleteness map calculates for all pixels $\mathbf{p}$ in the image the most incomplete co-direction:

$$ \mathbf{n}_{\infty}(\mathbf{p}) = \arg\max_{\mathbf{n}} I(\mathbf{p}, \mathbf{n}) \quad \forall \mathbf{n} \in S_{1/2}^2, \quad (1) $$

with $S_{1/2}^2$ the unit hemisphere, and the corresponding incompleteness value:

$$ I_{\infty}(\mathbf{p}) = I(\mathbf{p}, \mathbf{n}_{\infty}) = \max_{\mathbf{n}} I(\mathbf{p}, \mathbf{n}) \quad \forall \mathbf{n} \in S_{1/2}^2. \quad (2) $$

Numerically, for this 2D geometry, we only evaluated 720 homogeneously distributed co-directions of $S_{1/2}^2$ in the source trajectory plane, i.e. on the unit half circle. Fig. 2 illustrates the 2D incompleteness map for our DOA geometry with colored arrows, each giving two pieces of information: the direction encodes the direction $\mathbf{n}_{\infty}$ normal to the most incomplete line and the color encodes the incompleteness value $I_{\infty}$. In Fig. 2, the map is roughly separated in two regions: the upper half with horizontal co-directions and tomographic incompleteness values $I_{\infty}(\mathbf{p})$ in the range [0.27, 0.34], and the lower half with vertical co-directions and tomographic incompleteness values of about 0.40. This is expected as the hole between the two arcs of the source trajectory creates an incomplete direction defined by the midpoint of the arcs $(0, -510)$ for pixels in the upper half, while for pixels in the lower half, the incompleteness is larger in the horizontal direction. In summary, the studied geometry implies a significant tomographic incompleteness $I_{\infty}(\mathbf{p}) > 0.26 \forall \mathbf{p} \in \Omega$, the upper half is mainly incomplete in the vertical direction while the lower half is even more incomplete but in the horizontal direction.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Tomographic incompleteness map of the DOA geometry. The source points are represented in blue. The arrows show the most incomplete co-direction $\mathbf{n}_{\infty}$ and their color the corresponding tomographic incompleteness value $I_{\infty}$.</figcaption></figure>

The aim of the proposed algorithm is to account for this incompleteness map in the regularization of an iterative reconstruction algorithm adapted to the limited angle problem.

## B. Local Directional Total Variation

Data collected by the X-ray tomography system using our DOA geometry was modelled by the following discrete linear system:

$$ \mathbf{y} = \mathbf{A}\mathbf{x} + \epsilon \quad (3) $$

where $\mathbf{x} \in \mathbb{R}^{N=N_i \times N_j}$ is the sought attenuation map, $\mathbf{A} \in \mathbb{R}^{M \times N}$ the discrete X-ray transform of the DOA geometry, $\mathbf{y} \in \mathbb{R}^M$ the projection data and $\epsilon \in \mathbb{R}^M$ the associated noise. In a limited angle framework as the DOA geometry, reconstructing $\mathbf{x}$ with FBP leads to poor results due to the lack of data. In recent years, iterative reconstruction algorithms with an image prior have given the best results by constraining the problem to a space of plausible solutions [1]. To estimate $\mathbf{x}$, one can solve the following optimization problem:

$$ \mathbf{x}^* = \arg\min_{\mathbf{x}} \frac{1}{2} \|\mathbf{y} - \mathbf{A}\mathbf{x}\|_2^2 + \mathrm{R}(\mathbf{x}) \quad (4) $$

where the first term, the data fidelity term, ensures that the final solution $\mathbf{x}^*$ matches the projection data $\mathbf{y}$ and $\mathrm{R}$ is the regularization term embedding the image prior.

Total variation (TV) regularization, which promotes piece-wise constant solutions, has proven efficient for CT reconstruction, in particular for the limited angle problem [1]. Recently, Zhang et al. [2] developed DTV regularization based on the anisotropic version of TV, which weights differently the regularization on the horizontal and the vertical axes of the image. It can be written as:

$$ R_{\mathrm{DTV}}(\mathbf{x}) = \lambda \sum_{i=0}^{N_i-1} \sum_{j=0}^{N_j-1} \left| \beta(\mathbf{x}_{i,j} - \mathbf{x}_{i+1,j}) \right| + \left| \sqrt{1 - \beta^2}(\mathbf{x}_{i,j} - \mathbf{x}_{i,j+1}) \right| \quad (5) $$

where $\lambda \in \mathbb{R}^+$ controls the strength of the regularization with respect to the data fidelity term and $\beta \in [0, 1]$ is the directional weight. DTV is well adapted to the limited-angle problem where the lack of data is in the same direction for
every pixel of the reconstructed image. The parameter $\beta$ is tuned to regularize more this specific direction. However, it is not obvious to guess the best $\beta$ because it depends on the geometry and the scanned object. Moreover, DTV defines only one direction parameter while in some geometries, e.g. the DOA geometry, the incompleteness direction is not the same for every pixel in the image, as explained above.

We propose to use the tomographic incompleteness map to locally adjust the regularization. We formulate our LDTV regularization as follows:

$$R_{\text{LDTV}}(\mathbf{x}) = \sum_{i=0}^{N_i-1} \sum_{j=0}^{N_j-1} \tilde{I}_{\infty}(\mathbf{x}_{i,j}) (|\mathbf{n}_{\infty 2}(\mathbf{x}_{i,j}) (\mathbf{x}_{i,j} - \mathbf{x}_{i+1,j})| + |\mathbf{n}_{\infty 1}(\mathbf{x}_{i,j}) (\mathbf{x}_{i,j} - \mathbf{x}_{i,j+1})|) \quad (6)$$

with $\mathbf{n}_{\infty 1}$ and $\mathbf{n}_{\infty 2}$ the first and second components of the co-direction $\mathbf{n}_{\infty}$ and

$$\tilde{I}_{\infty}(\mathbf{x}) = \lambda_{\min} + \frac{\lambda_{\max} - \lambda_{\min}}{\max I_{\infty}} I_{\infty}(\mathbf{x}) \quad (7)$$

the tomographic incompleteness value linearly adjusted between the hyperparameters $\lambda_{\min} \in \mathbb{R}^+$ and $\lambda_{\max} \in \mathbb{R}^+$. The regularization is stronger when the incompleteness level is higher and the horizontal and vertical directional weights are defined by $\mathbf{n}_{\infty 2}$ and $\mathbf{n}_{\infty 1}$ respectively as they design the direction normal to the most incomplete line.

### C. Numerical experiments

The cost function of Eq. 4 is convex but non-smooth when dealing with TV-like regularizations. Aiming at a fair comparison between several regularization terms, we have used the same FISTA-TV optimization algorithm with non-negative constraints for all iterative methods [7].

We generated noiseless projections of the Forbild head phantom and corrupted them with pre-log Poisson noise with $10^6$ photons per pixel in air. We reconstructed images using the noiseless data with four algorithms: conventional FBP with a Hamming window and three iterative reconstruction algorithms with anisotropic TV, DTV and the proposed LDTV regularization terms. The two last have also been used to reconstruct images from noisy data. For the TV, DTV and LDTV, we performed 1200 iterations of FISTA with 60 denoising iterations which was enough to achieve the convergence according to the relative error with the ground truth (GT).

Image quality was quantified with the peak signal-to-noise ratio (PSNR). For the regularization terms, the hyperparameters were automatically tuned by minimizing the normalized root mean squared error (NRMSE) with the GT image using 600 iterations of the Nelder-Mead downhill simplex algorithm.

## III. RESULTS

The reconstructed images with noiseless projections are shown in Fig. 3. FBP reconstruction was severely impacted by the lack of data, and neither the inner left ear nor the profile were close to the GT. The high frequency noise in this result demonstrated that limiting the solution to be more piece-wise

constant was relevant. TV reconstruction, with $\lambda = 2.9 \times 10^{-4}$, really increased the visual quality of the image with a PSNR score of 28 dB. The image was more piece-wise constant, as expected, but there were residual artifacts at the border of the head, and it had difficulty to reconstruct high frequency information in the two ears. The best hyperparameters for DTV regularization were $\lambda = 2.8 \times 10^{-4}$, $\beta = 0.97$ and the final reconstruction had a PSNR score of 30 dB, notably better than the reconstruction with TV regularization. One can see the effect of the horizontal directional regularization on the left ear compared to the TV result. Also, the small disks in the right ear were closer to the GT. Finally, the reconstruction with LDTV regularization outperformed the other reconstruction algorithms with a high PSNR value of 36 dB confirmed by the visual impression of a higher image quality. The best hyperparameters for LDTV were $\lambda_{\min} = 3.3 \times 10^{-5}$ and $\lambda_{\max} = 2.6 \times 10^{-4}$. The borders of the head were better reconstructed, as one can notice above the left and right inner ears or at the bottom of the phantom. Moreover, the left inner ear showed 7 squares as the GT and the right ear was more accurate. Lastly, the profiles show that LDTV was able to reconstruct the two low contrast spheres at the top of the phantom unlike the other reconstruction algorithms.

DTV and LDTV images reconstructed from noisy data are shown in Fig. 4. DTV had a PSNR of 27 dB with $\lambda = 6.4 \times 10^{-4}$ and $\beta = 0.99$. Visually, noise clearly degraded the result. The reconstructed image with LDTV regularization had a PSNR of 30 dB. The optimal hyperparameters were $\lambda_{\min} = 2.7 \times 10^{-4}$ and $\lambda_{\max} = 9.3 \times 10^{-4}$ and the visual image quality was improved compared to DTV. The right ear and the border of the head were closer to the GT. However, the left ear and the low contrast profile were visually degraded compared to the result without noise.

## IV. DISCUSSION

The LDTV method uses the incompleteness map (Fig. 2) and therefore applies a stronger regularization in the lower half of the image than the upper half, in the horizontal directions and vertical directions respectively. DTV uses a single parameter and cannot locally adjust the regularization direction. The optimal directional hyperparameter of DTV was horizontal, which is consistent with the most incomplete direction of the incompleteness map and supports its use in LDTV.

For this work, we did not study the dependence of the regularization hyperparameters to the object. We have used an optimization algorithm to find the hyperparameters that minimize the NRMSE with the GT. The optimal hyperparameters will not only depend on the scanner geometry but also on the scanned object. The LDTV regularization adjusts the directional weights using the tomographic incompleteness, but the optimal hyperparameters $\lambda_{\min}$ and $\lambda_{\max}$ will probably depend on the characteristics of the scanned object.

The LDTV method gave satisfactory results in this 2D framework. The tomographic incompleteness maps deliver three-dimensional (3D) knowledge on the CT geometry. Computing LDTV regularization for a 3D images was beyond the scope of this work.
<figure><img src="image_4.png" /><figcaption>Fig. 3. GT and results of the four reconstruction algorithms using noiseless data with a gray window of [0.9, 1.2]. The left column shows the full reconstructed image with the corresponding PSNR value at the top-right. In the right column, the two zoom images represent the inner left and right ears, and the orange profile is taken along the line passing through the two low contrast spheres at the top of the phantom.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Reconstructed images using DTV and LDTV regularization from noisy projections. The gray window is [0.9, 1.2]. As in Fig. 3, the full reconstructed images with their corresponding PSNR values are in the left column beside the zooms on the ears and the plot through the two low contrast spheres.</figcaption></figure>

V. CONCLUSION

This article presented the LDTV, a new regularization using local tomographic incompleteness direction and level. To illustrate its effect, we have simulated and reconstructed noiseless and noisy data from the Forbild head phantom in a limited angle geometry. We have compared our method to state-of-the-art reconstruction algorithms. The results showed that LDTV regularization is a good solution to account for the acquisition geometry in the regularization term. LDTV may therefore be used to improve image quality of CT images reconstructed from projections acquired with a limited angle geometry.

REFERENCES

[1] E. Sidky, C.-M. Kao, and X. Pan, "Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT," *Department of Radiology, University of Chicago, 5841 S. Maryland Ave., Chicago, IL, 60637, USA*, vol. 14, pp. 119-139, Apr. 2006.

[2] Z. Zhang, B. Chen, D. Xia, E. Y. Sidky, and X. Pan, "Directional-TV algorithm for image reconstruction from limited-angular-range data," *Medical Image Analysis*, vol. 70, p. 102030, may 2021.

[3] B. Chen, Z. Zhang, D. Xia, E. Y. Sidky, and X. Pan, "Dual-energy CT imaging over non-overlapping, orthogonal arcs of limited-angular ranges," *Journal of X-Ray Science and Technology*, vol. 29, no. 6, pp. 975-985, 10 2021.

[4] M. Laurendeau, L. Desbat, G. Bernard, F. Jolivet, S. Gorges, F. Morin, and S. Rit, "Three-dimensional maps of the tomographic incompleteness of cone-beam ct scanner geometries," in *Fully 3D Reconstruction In Radiology and Nuclear Medicine*, 2023.

[5] R. Clackdoyle and F. Noo, "Quantification of tomographic incompleteness in cone-beam reconstruction," *IEEE Transactions on Radiation and Plasma Medical Sciences*, vol. 4, no. 1, pp. 63-80, Jan. 2020.

[6] Forbild. 2003. [Online]. Available: www.dkfz.de/en/roentgenbildgebung/ct/CT_Phantoms/Head_Phantom/Head_Phantom.html

[7] A. Beck and M. Teboulle, "Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems," *IEEE Transactions on Image Processing*, vol. 18, no. 11, pp. 2419-2434, nov 2009.