Multi-scale Image Fusion Between Pre-operative Clinical CT and X-ray Microtomography of Lung Pathology

TL;DR

This paper presents a multi-scale non-rigid registration method to align clinical CT and high-resolution micro-CT images of lung tissue, enabling detailed analysis of lung pathology at multiple scales.

Contribution

It introduces a novel multi-scale registration approach that effectively aligns clinical CT with micro-CT images of lung tissue, facilitating multi-resolution image fusion for medical analysis.

Findings

Reduced surface distance between CT and micro-CT from 3.3mm to 2.3mm

Successful alignment enabling future fusion with microscopic images

Method handles resection-induced volume differences

Abstract

Computational anatomy allows the quantitative analysis of organs in medical images. However, most analysis is constrained to the millimeter scale because of the limited resolution of clinical computed tomography (CT). X-ray microtomography ($\mu$CT) on the other hand allows imaging of ex-vivo tissues at a resolution of tens of microns. In this work, we use clinical CT to image lung cancer patients before partial pneumonectomy (resection of pathological lung tissue). The resected specimen is prepared for $\mu$CT imaging at a voxel resolution of 50 $\mu$m (0.05 mm). This high-resolution image of the lung cancer tissue allows further insides into understanding of tumor growth and categorization. For making full use of this additional information, image fusion (registration) needs to be performed in order to re-align the $\mu$CT image with clinical CT. We developed a multi-scale non-rigid registration approach. After manual initialization using a few landmark points and rigid alignment, several levels of non-rigid registration between down-sampled (in the case of $\mu$CT) and up-sampled (in the case of clinical CT) representations of the image are performed. Any non-lung tissue is ignored during the computation of the similarity measure used to guide the registration during optimization. We are able to recover the volume differences introduced by the resection and preparation of the lung specimen. The average ($\pm$ std. dev.) minimum surface distance between $\mu$CT and clinical CT at the resected lung surface is reduced from 3.3 $\pm$ 2.9 (range: [0.1, 15.9]) to 2.3 mm $\pm$ 2.8 (range: [0.0, 15.3]) mm. The alignment of clinical CT with $\mu$CT will allow further registration with even finer resolutions of $\mu$CT (up to 10 $\mu$m resolution) and ultimately with histopathological microscopy images for further macro to micro image fusion that can aid medical image analysis.

Tables

Table 1. TABLE I: Average minimum distances (AvgDist) in mm between the lung surface extracted from clinical CT and μ 𝜇 \mu CT before and after non-rigid registration.

AvgDist	before registration				after registration
[mm]	mean	std	min.	max.	mean	std	min.	max.
case 1	3.6	3.1	0.0	17.2	2.6	2.8	0.0	15.9
case 2	2.9	2.7	0.1	14.6	2.0	2.7	0.0	14.7
mean	3.3	2.9	0.1	15.9	2.3	2.8	0.0	15.3

Equations

$$\mathbf{T}(\vec{x})=\sum _{i,j,k}{\beta }^3(\frac{x}{{\delta }_x}-i)\times {\beta }^3(\frac{y}{{\delta }_y}-j)\times {\beta }^3(\frac{y}{{\delta }_z}-k)\times {\vec{\varphi }}_{ijk},$$

$$\begin{gathered} \mathcal{O}\left(I_{\mathrm{p}},I_{\mathrm{s}}\left(\mathbf{T}\right);\{\vec{\varphi }\}\right)=\left(1-\alpha -\beta -\gamma \right)\times \mathcal{S} \\ -\alpha \times {\mathcal{C}}_{\mathrm{smooth}}(\mathbf{T})-\beta \times {\mathcal{C}}_{\mathrm{volpres}}(\mathbf{T})-\gamma \times {\mathcal{C}}_{\mathrm{inconsistency}}(\mathbf{T}) \end{gathered}$$

$$\mathcal{S}\equiv \mathbf{NMI}=\frac{H(R)+H(F(\mathbf{T}))}{H(R,F(\mathbf{T}))}$$

$$\begin{gathered} {\mathcal{C}}_{\mathrm{smooth}}=\frac{1}{N}\sum _{\vec{x}\forall \Omega }({\left|\frac{{\partial }^2\mathbf{T}\left(\vec{x}\right)}{\partial x^2}\right|}^2+{\left|\frac{{\partial }^2\mathbf{T}\left(\vec{x}\right)}{\partial y^2}\right|}^2+{\left|\frac{{\partial }^2\mathbf{T}\left(\vec{x}\right)}{\partial z^2}\right|}^2 \\ +2\times \left[{\left|\frac{{\partial }^2\mathbf{T}\left(\vec{x}\right)}{\partial xy}\right|}^2+{\left|\frac{{\partial }^2\mathbf{T}\left(\vec{x}\right)}{\partial yz}\right|}^2+{\left|\frac{{\partial }^2\mathbf{T}\left(\vec{x}\right)}{\partial xz}\right|}^2\right]). \end{gathered}$$

$${\mathcal{C}}_{\mathrm{volpres}}=\frac{1}{N}\sum _{\vec{x}\forall \Omega }{\left[\log \left(\det \left(\mathrm{Jac}\left(\mathbf{T}\left(\vec{x}\right)\right)\right)\right)\right]}^2$$

$$\begin{gathered} {\mathcal{C}}_{\mathrm{inconsistency}}=\sum _{\vec{x}\forall \Omega }{\Vert {\mathbf{T}}_{\mathrm{forward}}\left({\mathbf{T}}_{\mathrm{backward}}\left(\vec{x}\right)\right)\Vert }^2 \\ +\sum _{\vec{x}\forall \Omega }{\Vert {\mathbf{T}}_{\mathrm{backward}}\left({\mathbf{T}}_{\mathrm{forward}}\left(\vec{x}\right)\right)\Vert }^2 \end{gathered}$$

Full text

Keywords: image fusion, non-rigid registration, computed tomography (CT), X-ray microtomography ($\mu$CT)

I INTRODUCTION

Clinical computed tomography (CT) is used for the diagnostic imaging of the living human (in-vivo imaging). As a result, most computational analysis is constrained to the millimeter scale because of the limited resolution of clinical CT. At this millimeter scale, pulmonary blood vessels and lung lobes can be observed. However, finer detailed anatomy is not observable. X-ray microtomograph ($\mu$CT) on the other hand allows imaging of ex-vivo tissues at a resolution of tens of microns. At this $\mu$m-scale, the alveoli and bronchiole regions can be clearly observed [1].

The prospective study of pre-operative imaging together with the high-resolution image analysis of resected tissue after surgery may provide us with unique opportunities to verify and potentially improve imaging protocols for cancer diagnostics [2, 3, 4].

In this work, we use clinical CT to image lung cancer patients before partial pneumonectomy (resection of pathological lung tissue). The resected specimen is then prepared for $\mu$CT imaging allowing the imaging of $\mu$m-scale anatomy. This study investigates the use of non-rigid intensity-based registration in order to establish a scale-seamless registration between clinical CT and $\mu$CT, with the ultimate aim of allowing a seamless navigation between anatomical scales inside the human body [1]. Related work is the registration of in-vivo and ex-vivo MRI of surgically resected specimens by [5, 2].

II METHOD

The resected specimen is prepared for $\mu$CT imaging at a voxel resolution of 50 $\mu$m (0.05 mm). This high-resolution image of the lung cancer tissue allows further insides into understanding of tumor growth and categorization. For making full use of this additional information, image fusion (registration) needs to be performed in order to re-align the $\mu$CT image with clinical CT. We developed a multi-scale non-rigid registration approach.

II-A Establishing scale-seamless registration between $\mu$CT and clinical CT

We use a non-rigid registration method to align the $\mu$CT specimen with the pre-operative clinical CT. This method is driven by the intensity similarity between the images. Any non-lung tissue is ignored during computation of the similarity measure $\mathcal{S}$ used to guide the registration during optimization. This is achieved via simple thresholding and morphological operations to extract the lung region in clinical CT and $\mu$CT as pre-processing step.

A coarse-to-fine approach is proposed in order to capture first the largest deformations and then the smaller differences between both images. This is achieved with a four-level multi-resolution pyramid as illustrated in Fig. 1. In order to compensate for the large differences in resolution between the $\mu$CT (in the case of $\mu$CT) and clinical CT images, we use a down-sample and up-sample (in the case of clinical CT) pyramid and perform the registration optimizations at each level in a coarse-to-fine fashion.

A good initialization for the registration algorithm is important. We use a few (3$\sim$5) manual point in order to provide affine alignment of the specimen with the pre-operative CT image. After affine initialization, we establish non-rigid alignment using the B-spline registration method, also known as 3D free form deformation based registration of [6] with the implementation provided by [7].

A 3D cubic B-Splines deformation model uses a lattice of control points $\{\vec{\phi}\}$. The spacing between each control point is uniform and denoted as $\delta_{x}$, $\delta_{y}$, and $\delta_{z}$ along the $x$-, $y$-, and $z$-axis respectively. For each voxel $\vec{x}$ in the domain $\Omega$ of the target image. the deformation $\mathbf{T}(\vec{x})$ can be computed as:

[TABLE]

where $\beta^{3}$ represents the cubic B-Spline function.

The images are aligned by finding the transformation which maximizes the following objective function:

[TABLE]

which combines a similarity measure, $\mathcal{S}$, and three penalty constraint terms, $\mathcal{C}_{\mathrm{smooth}}$, $\mathcal{C}_{\mathrm{volpres}}$, and $\mathcal{C}_{\mathrm{inconsistency}}$. Each term is weighted against each other by user-defined weights $\alpha$, $\beta$, and $\gamma$.

The similarity measure used between the reference ($R$) and floating image ($F$) is normalized mutual information (NMI):

[TABLE]

where $H(R)$ and $H(F(\mathbf{T}))$ the two marginal entropies, and $H(R,F(\mathbf{T})$) is the joint entropy. Its computation requires a joint histogram which is filled by using a Parzen Window (PW) approach [7, 8].

The three constraint terms are used to encourage realistic deformations. The bending energy describes the smoothness of the deformation and is defined as:

[TABLE]

The volume-preserving penalty term discourages large expansions/contractions, and is defined as:

[TABLE]

In addition we prevent the occurrence of folding in the transformation using a folding correction scheme [9]. For each transformed voxel that would cause a negative Jacobian determinant, its influence on its neighborhood control points is computed. The control point positions is then changed until the determinant value is positive.

In an ideal case, the transformations from $F$ to $R$ (forward) and $R$ to $F$ (backward) are the inverse of each other, e.g. $\mathbf{T}_{\mathrm{forward}}=\mathbf{T}_{\mathrm{backward}}^{-1}$ and $\mathbf{T}_{\mathrm{backward}}=\mathbf{T}_{\mathrm{forward}}^{-1}$ [10]. Hence, we include a penalty term that encourages inverse consistency of both transformations. We follow the approach of [11] using compositions of $\mathbf{T}_{\mathrm{forward}}$ and $\mathbf{T}_{\mathrm{backward}}$ and add

[TABLE]

The following parameters were found empirically by visual examination of the registration results. We use a four-level multi-resolution pyramid with a maximum of 500 iterations per level. Both the image and B-spline control point grid resolutions are doubled with increasing resolution levels. The final control point spacing between voxels is 5. The objective function weights are set to $\alpha=10^{-4}$, $\beta=10^{-12}$, and $\gamma=0.1$. These parameters were found to recover the majority of the deformation between the two images, while preventing unrealistic deformations from occurring. We used the open-source software111http://sourceforge.net/projects/niftyreg by [7, 12] for this study.

III RESULTS

We collected pre-operative clinical CTs from two patients before partial pneumonectomy. The clinically acquired CT images have dimensions of [512, 512, 435$\sim$554], and voxel spacings of [0.625, 0.625, 0.6] mm. The $\mu$CT images used for registration had dimensions of [278$\sim$512, 278$\sim$512, 346$\sim$538] with isotropic voxel spacings of $0.111\sim 0.127$ mm. Note that we downsampled the original $\mu$CT images (circa 50$\mu$m resolution) by a factor of two for the registration experiments. The tube voltage for $\mu$CT was 90kVp, and tube current was 110 $\mu$A.

3$\sim$5 corresponding points where chosen by an expert clinician (SN) in order to provide an initial affine alignment of the specimen with the pre-operative CT image. We then crop the clinical CT image to the extent of the aligned $\mu$CT image for subsequent non-rigid registration.

Figure 2 shows the alignment of $\mu$CT lung specimen with pre-operative clinical CT before and after non-rigid registration for one case. A qualitatively better alignment of $\mu$CT with the target clinical CT regions can be observed.

In order to quantitatively evaluate the registration performance, we measure the average minimum surface distance (AvgDist) between the lung surface extracted from clinical CT and $\mu$CT before and after non-rigid registration. Figure 3 shows the extracted surfaces used for measurement before and after non-rigid alignment. The AvgDist measures are given in Table I. A reduction from 3.3 $\pm$ 2.9 (range: [0.1, 15.9]) to 2.3 mm $\pm$ 2.8 (range: [0.0, 15.3]) mm on average can be observed. This is a significant improvement with p $<$ 0.001 (Wilcoxon signed rank test).

IV CONCLUSIONS

We presented a method for non-rigid alignment between images of pre-operative clinical CT and x-ray microtomography ($\mu$CT) of lung pathology. After manual initialization using a few landmark points and affine alignment, several levels of non-rigid registration between down-sampled (in the case of $\mu$CT) and up-sampled (in the case of clinical CT) representations of the image are performed. This allows us to recover the volume differences introduced by the resection and preparation of the lung specimen. The proposed multi-scale image fusion approach will allow further registration with even finer resolutions of $\mu$CT (up to 10 $\mu$m resolution) [13] and ultimately with histopathological microscopy images for further macro to micro image fusion that can aid medical image analysis.

ACKNOWLEDGEMENTS

This paper was supported by MEXT KAKENHI (26108006 and 15H01116).