# Optimisation of reconstruction for the registration of CT liver perfusion sequences

Blandine Romain, Veronique Letort, Olivier Lucidarme, Florence

d'Alché-Buc, Laurence Rouet

► To cite this version:

Blandine Romain, Veronique Letort, Olivier Lucidarme, Florence d'Alché-Buc, Laurence Rouet. Optimisation of reconstruction for the registration of CT liver perfusion sequences. SPIE Medical Imaging 2012, Feb 2012, San Diego, United States. 10.1117/12.910929. hal-00744964

HAL Id: hal-00744964

https://hal.science/hal-00744964

Submitted on 30 Aug 2022

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
# Optimisation of reconstruction for the registration of CT liver perfusion sequences

B. Romain<sup>abc</sup>, V. Letort<sup>b</sup>, O. Lucidarme<sup>e</sup>, F. d'Alché-Buc<sup>cd</sup>, L. Rouet<sup>a</sup>

<sup>a</sup> Medisys Research Lab, Philips Healthcare, Suresnes, France

<sup>b</sup> Ecole Centrale Paris, MAS, Chatenay-Malabry, France

<sup>c</sup> IBISC EA 4526, Université d'Evry, Evry, France

<sup>d</sup> INRIA (TAO/AMIB), LRI CNRS 8623, Orsay, France

<sup>e</sup> AP-HP, Hospital La Pitié-Salpétrière, Paris, France

## ABSTRACT

**Objective.** CT abdominal perfusion is frequently used to evaluate tumor evolution when patients are undergoing anti-tangiogenic therapy. Parameters depending on longer-term dynamics of the diffusion of the contrast medium (e. g. permeability) could help assessing the treatment efficacy. To this end, dynamic image sequences are obtained while patients breath freely. Prior to any analysis, one needs to compensate the respiratory motion. The goal of our study is to optimize the CT reconstruction parameters (filter of reconstruction, thickness of image volumes) for our registration method. We also aim at proposing relevant criteria allowing to quantify the registration quality. **Methods.** Registration is computed in 4 steps: z-global rigid registration, local refinements with multiresolution blockmatching, regularization and warping. Two new criteria are defined to evaluate the quality of registration: one for spatial evaluation and the other for temporal evaluation. **Results.** The two measures decrease after registration (58% and 10% average decrease for the best reconstruction parameters for the spatial and temporal criteria respectively) which is consistent with visual inspection of the images. They are therefore used to determine the best combination of reconstruction parameters.

**Keywords:** spatio-temporal registration, perfusion, free-breathing, validation criteria, blockmatching, multiresolution

## 1. INTRODUCTION

Liver cancer can be treated by tumor resection surgery or, for inoperable patients, by tumor anti-angiogenesis drug treatments. To increase the survival rate of inoperable patients, one main challenge is to develop tools for clinicians allowing the earliest possible evaluation of tumor response to a given drug, in order to consider a potential change of treatment in case of negative response. Anatomical measurements, such as tumor size, are currently the main criteria to perform the assessment of disease evolution (RECIST 1.1).¹ However, anatomical changes may be insensitive to, or provide markedly delayed indications of, the response to treatment. Functional imaging of the microcirculation, based on the acquisition of time sequences combined with contrast injection, provides additional in vivo information. Depending on the model used for perfusion, the acquisition duration may vary from 30 seconds up to 5 minutes. Parameters depending on longer-term dynamics of the diffusion of the contrast medium (e.g. permeability) can be derived to assess treatment efficiency. Providing time dependent parameters such as contrast intake curves is challenging in the context of liver, due to the motion of the organ during acquisitions. Two main approaches are possible to deal with respiratory motion, namely breath hold (during one or successive periods)²⁻⁴ and free breathing. Free breathing is selected in the current study, since a less stressed patient usually provides sequences with more regular motion. So, a step allowing to compensate the respiratory motion is necessary.

Methods for registration of images acquired on free-breathing patients are being developed for contrast-enhanced MRI of pulmonary⁵ or myocardial,⁶ but their adaptation to other modalities and organs is not straightforward. The current study focuses on tumor characterization of liver or kidney using dynamic contrast-enhanced CT (DCE-CT). As opposed to MRI, three main constraints have to be taken into account. First, the image field of view is limited to about 6 cm in z-direction, perpendicular to axial plane. This limit has an important consequence for motion compensation, since the main motion of the liver is also in z-direction, and its maximal amplitude is
about 2 to 3 cm. It implies that in most cases, only truncated motions of liver can be acquired within the image field of view. Second, the use of time sequences of CT requires to limit as much as possible the radiation dose for each acquisition. As a consequence, images are acquired at sparse times. Third, images are noisy also for dose consideration. However, noise levels can be different depending on the type of CT reconstruction that is used.

The goal of the present study is to evaluate the impact of the CT reconstruction parameters on the quality of the registrations. For this purpose, sequences have been reconstructed with various parameters. Two new criteria, a spatial one (Area Ratio) and a smoothness one (Smoothness Criterion), are proposed to assess the registration quality.

## 2. MATERIALS AND METHODS

### 2.1 Data acquisition

8 data acquisitions are performed with a 256 slices CT (Brilliance iCT 256, Philips Healthcare, The Netherlands) on 6 patients. All patients gave informed consent. The dynamic CT protocol acquires 48 images every 2.5 seconds (80 kV, 80 mAs, rotation time of 0.33 seconds, dose of iodine 80 cc). With 8 cm detector coverage in z axis, the CT scanner allows for 5.5 cm effective z coverage in a single rotation with a 3D axial cone beam correction. This coverage is large enough to keep most of the lesion in the field of view (350 mm). In the axial plane, data are reconstructed with a pixel size of 0.68 mm × 0.68 mm (Figure 1).

<figure><img src="image_2.png" /><figcaption>Figure 1. Image data at time 0, start of injection, (on the left) and 15 seconds after injection (on the right). On the right figure: kidney is enhanced in white square (right) and aorta with hepatic artery are enhanced in white rectangle (middle).</figcaption></figure>

In order to evaluate the impact of noise on the registration quality, two reconstruction parameters are evaluated: slice thickness and filter of reconstruction.

Regarding slice thickness, use of a multi-slice scanner offers the possibility to vary slice thickness during reconstructions. We evaluate slice thicknesses of 2.5 mm and 1.25 mm. Slice increment, another reconstruction parameter, is fixed to half of the slice thickness value in order to enable imaging of small structures.

During CT reconstruction, smoothing filters may be applied in the sinogram domain. The amount of smoothing can be selected during reconstruction. Two levels of filtering are evaluated: smooth (filter 1) and sharp (filter 2). The level of noise increases when moving from filter 1 to filter 2. Examples of noise in the liver are presented in table 1.
<table><caption>Table 1. Example of noise inside a region of interest in the liver depending on the reconstruction parameters</caption><thead><tr><td>Thickness</td><td>Filter 1</td><td>Filter 2</td></tr></thead><tbody><tr><td>2.5 mm</td><td>5.51</td><td>8.07</td></tr><tr><td>1.25 mm</td><td>11.24</td><td>18.02</td></tr></tbody></table>

## 2.2 Registration

Each reconstruction provides a dynamic sequence $DS = (I_1, ..., I_T)$, where $I_t$ ($t \in (1, ..., T)$) is a 3D image acquired at time $t$ and $T$ the number of acquisition times. We note $I_r$, the reference volume for registration.

### 2.2.1 Masking

To improve robustness and computation time of the registration algorithm (see 2.2.2), a preprocessing step is applied, that consists in computing masked images (see figure 2). Indeed, bones (namely spine and ribs) have different motion from the rest of the abdomen and can introduce perturbations in the registration of the targeted organs. Removing these regions from the images to register helps improving the robustness of the algorithm. To decrease the computation time, a mask is also applied on the background which does not provide information and may also introduce errors in the registration.

The mask $M_t$ to be applied on the volume $I_t$ can be defined as $R \setminus S$ where $R$ and $S$ are two nested ellipsoidal cylinders:

* the internal cylinder corresponds to the spine and is called $S$

* the external cylinder corresponds to the interior of the abdomen and is delimited by the ribs, that mark the frontiers between the abdomen and the background; it is called $R$

The contours of these two zones are approximated by two 2D ellipses calculated within an arbitrary axial plan and extended to the 3D volumes. To compute the equations of these ellipses, we first need to extract the centers and some points belonging to the contour of each of these two areas.

1. Contour points extraction:

First, a threshold of 1600 HU is applied on the volume to extract a set $B$ of voxels of bones (spine contour and ribs). The contour of the spine $S$ is defined as the largest 3D connected component. These points are projected in an arbitrary axial plane, thus defining the set of points of the contour of $C_S = \{(x_i, y_i), i = 1...N_S\}$ with $N_S$ the number of contour points of $S$.

Second, we consider the remaining extracted pixels of $B$: they are considered as belonging to ribs and, once projected in an arbitrary axial plane, define the set of points $C_R = \{(x_i, y_i), i = 1...N_R\}$ of the contour of $R$ with $N_R$ the number of contour points.

2. Centers computation:

For $S$, the center $(x_c, y_c)_S$ is defined as the barycenter of the contour points $C_S$.

As regards the center of $R$, a threshold of -500 HU is applied to consider regions of air in the volume. The largest connected component is the background. The center $(x_c, y_c)_R$ is assumed to be the 2D barycenter of pixels that not belong to the background.

3. Ellipse fitting:

Computing the equation of an ellipse, once its center $(x_c, y_c)$ and $N$ points of its contour $(x_i, y_i)_{i=1...N}$ are known, comes down to minimizing the distance from these points to the ellipse, as follows:
$$ (\alpha_{opt}, a_{opt}, b_{opt}) = \arg\min_{a,b,\alpha} \sum_{i=1}^{N} \left| \frac{((x_i - x_c)\cos\alpha + (y_i - y_c)\sin\alpha)^2}{a^2} + \frac{((x_i - x_c)\sin\alpha - (y_i - y_c)\cos\alpha)^2}{b^2} - 1 \right|^2 \quad (1) $$

where a and b are respectively the major and minor axis of the ellipse, and α is the angle between the two axis.

Let us note the vectors $X_c = (x_i - x_c)_{i=1...N}$ and $Y_c = (y_i - y_c)_{i=1...N}$. We can define the matrix of points as $M = (X_c^2 Y_c^2 X_c Y_c)$ and the vector of parameters $P = (P_1 P_2 P_3)$

where $P_1 = (\frac{\cos^2\alpha}{a^2} + \frac{\sin^2\alpha}{b^2})$, $P_2 = (\frac{\sin^2\alpha}{a^2} + \frac{\cos^2\alpha}{b^2})$ and $P_3 = -(2\cos\alpha\sin\alpha(\frac{1}{a^2} - \frac{1}{b^2}))$.

So, the problem of minimization can be written under the following form:

$$ P_{opt} = \arg\min_{P} \|M.P - 1\|_2 \quad (2) $$

A singular value decomposition allows solving this equation (2).

With the two ellipses on an axial plane, we propagate these ellipses on each plane of volumes. Finally, the 3D mask consists of an ellipsoidal cylinder with extrusion of an internal ellipsoidal cylinder (figure 2).

<figure><img src="image_3.png" /><figcaption>Figure 2. Mask example: white zone is the valid region, excluding spine, ribs and background</figcaption></figure>

After computing a mask for each image volume $I_t$, we define the domain to consider for the registration $\Omega_{r,t} = M_r \cap M_t$ as the intersection of the reference image mask $I_r$ and the current image mask $I_t$.

### 2.2.2 Spatio-temporal registration

Due to the contrast medium injection, intensity of a given tissue may vary between different images of the sequence. As a consequence, the choice of the reference volume on which to register the other images is difficult. To avoid the risk of error accumulation by using multiple references, we use a single reference volume $I_r$ roughly corresponding to the enhancement phase (1.5 minutes after injection).

Our registration approach involves four main steps:

1. **Global z-translation**: since the main motion is in the z-direction, the first step consists in a global evaluation of the z-translation. We select half of reference image volume $I_r$, and we search the sub volume of $I_t$ that best matches according to a similarity measure. So, the search range is 27.5 mm (half of volume thickness). To increase the robustness of the algorithm to contrast enhancement, we choose difference entropy (DE) as similarity measure.⁷ Difference between $I_r$ and $I_t$, noted $I_{r-t}$, is used to evaluate DE. The histogram $H(I_{r-t})$ of $I_{r-t}$ normalized by number of image voxels $N$ with $B$ bins is calculated and the difference entropy is defined by:
$$DE = - \sum_{i=1}^{B} H(I_{r-t})_i * \log(H(I_{r-t})_i) \qquad (3)$$

2. **Multi-resolution blockmatching:** initialized by the z-translation found in the first step, a 3D registration method⁸ is computed. The following parameters are selected: three multi-resolution levels in (x,y), one in z; (2.8 mm, 2.8 mm, 10 mm) search range in (x, y, z) directions; block sizes of (22 mm, 22 mm, 7 mm). Finally, we obtain, for each block, motion vector fields corresponding to translations in x, y and z dimensions. The multiresolution approach allow to regularize a little the vector motion fields and to have a better computation time. The blockmatching algorithm is used in our method to refine local displacements, due to independent movements of organs and also to little elastic motion in the liver.

3. **Regularization:** a regularization step is necessary to smooth the motion vector fields. A Gaussian filter ($\sigma = 1.2$) is applied in 3D.

4. **Warping:** warping consists in reconstructing new image sequences ($J_1, ..., J_T$) with less respiratory motion. A trilinear interpolation is used.

## 2.3 Evaluation criteria

Since the registration has spatial and temporal dimensions, we propose two criteria to validate our method.

* **Spatial evaluation criterion:**

A criterion, named Area Ratio (AR), is based on Sum of Square Differences (SSD) between $J_t$ and $I_r$ (see Figures 3 and 4).

$$AR = \int_{0}^{T} \frac{\sum_{(x,y,z) \in \Omega_{r,t}} |J_t(x,y,z) - I_r(x,y,z)|^2}{\sum_{(x,y,z) \in \Omega_{r,t}} |I_t(x,y,z) - I_r(x,y,z)|^2} dt \qquad (4)$$

* **Temporal evaluation criterion:**

A region of interest within the left kidney is defined. Kidney is selected since its presence in a given axial slice is highly sensitive to motion. Inside this region, mean intensities $\overline{J}_t$ and $\overline{I}_t$ are calculated after and before registration, for each acquisition time. Since the size of this region is large (> 1000 pixels), $\overline{J}_t$ and $\overline{I}_t$ can be considered noise free. The quality of the registration is assessed by the smoothness of the curve after registration. Note that there are some intrinsic low-frequency variations, often sigmoid-shaped, due to the diffusion of the contrast medium. Thus, we propose a new evaluation criterion based on curve smoothness: the two curves are independently filtered by a median time-filter with a window span of 3 points and the sum of the absolute difference between the smoothed curve and the original curve is computed. The ratio of the smoothness after registration over before registration is used as a second temporal criterion, named Smoothness Criterion (SC) (equation 5).

$$SC = \frac{\sum_{t=1}^{T} |\overline{J}_t - \text{Median}(\overline{J}_t, 3)|}{\sum_{t=1}^{T} |\overline{I}_t - \text{Median}(\overline{I}_t, 3)|} \qquad (5)$$
### 3. RESULTS

Differences between the image $I_t$ and the reference image $I_r$ are lower after registration than before, as illustrated in Figures 3 and 4.

<figure><img src="image_4.png" /><figcaption>Figure 3. Norm L2 of differences between current image $I_t$ and reference image $I_r$ with respect to time before registration (dash line) and after registration (solid line)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Images of difference: on the top, before registration; on the bottom, after registration</figcaption></figure>

In Figure 5, the corresponding curves of mean intensities in the left kidney are smoother after registration than before (see also Table 2).

<figure><img src="image_6.png" /><figcaption>Figure 5. Mean intensities inside a region of interest (left kidney) on an axial slice as a function of time of acquisition, before registration (dash line) and after registration(solid line) for CT5.</figcaption></figure>
Table 2 allows the identification of the best set of parameters for the two criteria, for each examination (boldface). Low values of these criteria correspond to images of very good quality in term of visual inspection. The two criteria lead to similar conclusions. With criterion SC, the best results of registration are obtained with thickness 2.5 mm (Th) and filter 1 (F) for all examinations. It ranges from 0.18 to 0.71, with a mean of 0.42 which represent a 58% improvement of the smoothness criterion from the original to the registered sequence. With criterion AR, the best parameter set also included a thickness of 2.5 mm for all the examinations, but no clear difference can be made between the two filters.

<table><thead><tr><th colspan="2">Exams</th><th colspan="2">CT1</th><th colspan="2">CT2</th><th colspan="2">CT3</th><th colspan="2">CT4</th></tr><tr><th>Th</th><th>F</th><th>SC</th><th>AR</th><th>SC</th><th>AR</th><th>SC</th><th>AR</th><th>SC</th><th>AR</th></tr></thead><tbody><tr><td rowspan="2">2.5 mm</td><td>1</td><td><strong>0.46</strong></td><td><strong>0.82</strong></td><td><strong>0.56</strong></td><td><strong>0.95</strong></td><td><strong>0.71</strong></td><td><strong>0.93</strong></td><td><strong>0.19</strong></td><td><strong>0.96</strong></td></tr><tr><td>2</td><td>0.61</td><td>0.86</td><td>0.84</td><td><strong>0.95</strong></td><td>1.16</td><td>0.97</td><td>0.64</td><td>0.96</td></tr><tr><td rowspan="2">1.25 mm</td><td>1</td><td>1.14</td><td>0.96</td><td>1.13</td><td>0.99</td><td>1.32</td><td>0.99</td><td>0.92</td><td>0.99</td></tr><tr><td>2</td><td>0.99</td><td>0.97</td><td>1.18</td><td>0.99</td><td>1.34</td><td>1.00</td><td>0.88</td><td>0.99</td></tr></tbody></table>

<table><caption>Table 2. Evaluation criteria on the 8 examinations: ratio of SC (resp. AR) after registration to before registration. Th: thickness; F: filter.</caption><thead><tr><th colspan="2">Exams</th><th colspan="2">CT5</th><th colspan="2">CT6</th><th colspan="2">CT7</th><th colspan="2">CT8</th></tr><tr><th>Th</th><th>F</th><th>SC</th><th>AR</th><th>SC</th><th>AR</th><th>SC</th><th>AR</th><th>SC</th><th>AR</th></tr></thead><tbody><tr><td rowspan="2">2.5 mm</td><td>1</td><td><strong>0.69</strong></td><td><strong>0.80</strong></td><td><strong>0.18</strong></td><td><strong>0.86</strong></td><td>0.17</td><td><strong>0.90</strong></td><td>0.32</td><td><strong>0.98</strong></td></tr><tr><td>2</td><td>0.73</td><td>0.89</td><td><strong>0.18</strong></td><td>0.90</td><td><strong>0.16</strong></td><td>0.91</td><td><strong>0.19</strong></td><td>1.06</td></tr><tr><td rowspan="2">1.25 mm</td><td>1</td><td>0.99</td><td>0.95</td><td>0.81</td><td>0.95</td><td>0.91</td><td>0.96</td><td>1.04</td><td>0.99</td></tr><tr><td>2</td><td>1.00</td><td>0.98</td><td>0.91</td><td>0.98</td><td>0.92</td><td>0.97</td><td>1.03</td><td>0.99</td></tr></tbody></table>

Figure 6 show the results of registration using a checkerboard. After registration, frontiers of organs are better aligned.

<figure><img src="image_7.png" /><figcaption>Figure 6. Checkerboards between image of reference (yellow) and current image (gray). At the left: before registration; at the right: after registration</figcaption></figure>

## 4. CONCLUSION

We propose two new criteria for quantifying the registration quality and allowing to optimize the parameters of reconstruction.
Results show that a thickness reconstruction parameter of 2.5 mm provides the best registration. This recon-
struction corresponds to the reconstruction with the lowest level of noise. Moreover, as noise is less affected by
the filter parameter than thickness parameter, the filter parameter has less impact on the registration quality.
Regarding the computation time, our method is relatively quick. To register one image volume of (512*512*39)
pixel size, it takes about 4 seconds (Intel(R) Xeon(R) CPU, E5430, 2.66 GHz).

ACKNOWLEDGMENTS

We thank Rafael Wiemker for his technical assistance, Philippe Coulon for his implication, Raphael Prevost and Vincent Auvray for their discussions.

REFERENCES

[1] Eisenhauer, E. A., Therasse, P., Bogaerts, J., Schwartz, L. H., Sargent, D., Ford, R., Dancey, J., Arbuck, S., Gwyther, S., Mooney, M., Rubinstein, L., Shankar, L., Dodd, L., Kaplan, R., Lacombe, D., and Verweij, J., "New response evaluation criteria in solid tumours: revised RECIST guideline (version 1.1).," *Eur J Cancer* **45**, 228-247 (Jan 2009).

[2] Sahani, D. V., Holalkere, N.-S., Mueller, P. R., and Zhu, A. X., "Advanced hepatocellular carcinoma: CT perfusion of liver and tumor tissue–initial experience.," *Radiology* **243**, 736–743 (Jun 2007).

[3] Koh, T. S., Thng, C. H., Hartono, S., Lee, P. S., Choo, S. P., Poon, D. Y. H., Toh, H. C., and Bisdas, S., "Dynamic contrast-enhanced CT imaging of hepatocellular carcinoma in cirrhosis: feasibility of a prolonged dual-phase imaging protocol with tracer kinetics modeling.," *Eur Radiol* **19**, 1184–1196 (May 2009).

[4] Ng, C. S., Chandler, A. G., Wei, W., Herron, D. H., Anderson, E. F., Kurzrock, R., and Charnsangavej, C., "Reproducibility of CT perfusion parameters in liver tumors and normal liver.," *Radiology* **260**, 762–770 (Sep 2011).

[5] Tokuda, J., Mamata, H., Gill, R. R., Hata, N., Kikinis, R., Padera, Jr, R. F., Lenkinski, R. E., Sugarbaker, D. J., and Hatabu, H., "Impact of nonrigid motion correction technique on pixel-wise pharmacokinetic analysis of free-breathing pulmonary dynamic contrast-enhanced MR imaging.," *J Magn Reson Imaging* **33**, 968–973 (Apr 2011).

[6] Wollny, G., Ledesma-Carbayo, M. J., Kellman, P., and Santos, A., "A new similarity measure for non-rigid breathing motion compensation of myocardial perfusion MRI.," *Conf Proc IEEE Eng Med Biol Soc* **2008**, 3389–3392 (2008).

[7] Buzug, T. M., Weese, J., Fassnacht, C., and Lorenz, C., "Image registration: convex weighting functions for histogram-based similarity measures," in [*CVRMed*], 203–212 (1997).

[8] Tzovaras, D., Strintzis, M. G., and Sahinoglou, H., "Evaluation of multiresolution block matching techniques for motion and disparity estimation," *Signal Processing: Image Communication* **6**(1), 59 – 67 (1994).