Optical-CT Dual-Modality Mapping Base on DRR Registration

.


Introduction
Fluorescence Molecular Tomography (FMT) has been used frequently for the noninvasive study of cancer biological behavior [1][2][3]. With specific imaging agents, FMT has powerful capability in detection and quantification of important biological processes by reconstructing in vivo distribution of agents [4,5]. Combining FMT with computerized tomography (CT), both functional signals and anatomical information can be obtained. However, when FMT/CT dual-modality imaging technology is used to calculate internal bioluminescent source's size and location, it is necessary to establish a mapping between two-dimensional (2D) optical images and three-dimensional (3D) CT data in order to reconstruct 3D energy distribution on surface of the imaging object [6,7].
The registration of the dual-modality images is the basis of mapping. The registration between the optical image and the CT volume data is a registration between the 2D and 3D images. Markers are often used as indicators in the registration of different modality images [8]. The registration of 2D and 3D can also be transformed into the registration of 3D and 3D through projecting reconstruction techniques. Ning [9] reconstructs 3D optical surfaces by acquiring multi-orientation 2D images and then complete the 3D registration between optical image and CT image. However, multi-orientation optical image's acquisition is complicated and variations of the fluorescent energy over time will reduce the reconstruction accuracy of the internal bioluminescent source. In addition, 3D image can also be projected into 2D image to achieve the registration between two modal 2D images. However, this approach is works well for two modal images with high similarity, such as X-ray images and CT [10,11].
Fluorescent image is used along with CT to map functional information on tissue structures due to its high sensitivity. However, 2D fluorescent image acquired from CCD is short of three-dimensional coordinate information. In this paper, we propose an optical-CT dual-modality image mapping algorithm based on Digitally Reconstructed Radiograph (DRR) registration. The algorithm can directly map the 2D optical image onto the surface of the 3D CT data to obtain the 3D light intensity distribution on the tissue surface.

DRR Virtual Image
The 3D CT data can be converted to the 2D DRR virtual images through the ray casting algorithm [12,13]. As shown in Fig. 1, the virtual point light source simulates a conventional X-ray light source and emits a number of virtual rays passing through the CT dataset to project on the DRR imaging plane, and the projected points of all rays on the imaging plane constitute a DRR virtual image. For a given CT volume data J, the DRR virtual image can be obtained from the following model [14]: ( ) is the intensity of the DRR virtual image at , is the 3D geometric transformation of the CT volume data, ( , ) is the ray casting model, is a set of parameter points on the light.

Registration of Optical Images and DRR Virtual Image
For the fluorescence imaging system (AOIS, Animal Optical Imaging System, Nanjing University of Aeronautics and Astronautics, China), the white-light image and the fluorescence image are located in the same coordinate system. The white-light image provides the external contour information of the imaging object, and the fluorescence image provides light intensity distribution. Therefore, registration of 2D fluorescent image and 3D CT image can be achieved by registration of 2D white-light image with external contour and 3D CT image. The white-light image and the DRR virtual image computed from CT are aligned iteratively until an optimal matching is achieved by adjusting the optimal transformation parameters and similarity measure. The flow chart based on DRR registration of these two modalities images is shown in Fig. 2.

Similarity Measure Registration Based on Mutual Information
In medical image registration, normalized mutual information (NMI) is widely used in multi-modality image registration because of its high registration accuracy and good robustness. Given two images and , the NMI between them is defined as [15]: where ( ) and ( ) represent the entropies of images and , respectively, and ( , ) represents the joint entropy of images and , which are defined as follows: where ( ) and ( ) respectively represent the gray-level distribution of images and B.
, ( , ) is the joint probability distribution of gray level.

Improved Similarity Measure Registration
The registering accuracy based on the NMI is often suffered from image missing, noise, and excessive local extremes. This paper proposes an improved Hausdorff Normalized Mutual Information (HNMI) similarity measure for optical image and DRR virtual image registration, which combines image's mutual information and contour features.
After the average filtering and Maximum Between-Class Variance (OTSU) segmenting, the contour feature point sets of images A and B are obtained by sampling the contours according to uniform angle from centroid individually. Hausdorff distance between contour feature point sets A and B is defined as follows [16]: where ( ) = ∈ ‖ − ‖, ( ) ( ) represents the th distance after sorting, ℎ is in the range from 0.6 to 0.8, is the number of points in the set ′. So, for images A and B to be registered, the improved HNMI function is defined as: where 1 and 2 are constants used to balance the contributions of ( ′, ′) and ( , ). After doing a lot of tests, 1 is set 0.2 and 2 is set 1 in the following phantom and mouse experiments.
( , ) is the NMI of two images, and ( ′, ′) is the Hausdorff distance between point sets ′ and ′.

Multi-Resolution Registration Strategy Based on Gaussian Pyramid
In order to avoid the side effect of local optimization on HNMI and to improve the registration accuracy and operating efficiency, Gaussian Pyramids are introduced for multi-resolution image registration [17]: ( , ) represents the th layer image, and respectively are the row and column of the image, ( , ) is a window function.
The process of Gaussian Pyramid for multi-resolution registration is shown in Fig. 3. Optical image and DRR virtual image are respectively resampled to obtain three sets of image sequences with different resolutions. The rough registration is obtained from the image with the lowest resolution in the top layer, and the transformation parameters are used as the initial parameters of the next layer of image registration.

Figure 3: Multi-resolution registration based on Gaussian pyramid
With respect to the selection of the optimization algorithm, particle swarm optimization (PSO) has low convergence precision but strong global optimization ability, and powell optimization has fast convergence rate but is easily dropped in local minimization [18]. The contour spindle axis method is usually used for rough registration of different modal images [19]. Therefore, after roughly alignment using contour spindle axis method in the top layer, contour similarity distance measure and the PSO algorithm are chosed to align two modal images in the second layer. The contour similarity distance measure is defined as follows: where × is the resolution of image, ( , ) is the pixel value of the image.
In this end, optical and DRR images are precisely registered by combining the improved HNMI measure and the Powell algorithm.

Registration Evaluation Index
To evaluate the proposed registration method, some marked points were set on the surface of the imaging object. The mean error (ME) and root mean squared error (RMSE) are often used as evaluation indexes of the registration. ME and RMSE are defined as follows: At a marked point , ( ) is the registration coordinate and ( ) is the real coordinate. is the number of marked points.
The optimal registration parameters can be obtained after the registration between DRR virtual image and the optical image. Then, according to the coordinate correspondence relationship between DRR virtual image and 3D CT data, the fluorescence distribution can be reconstructed on the surface of CT volume data.

Free Space Flux Mapping Model
The fluorescent photons propagate along a straight line in free space, and the Lambert's cosine law is used to describe the transmission of photons from the surface of the imaging object to the CCD detectors [20].
The light flux reconstruction on the surface of the object is to simplify the imaging lens in the optical imaging system into the pin-hole imaging model [21]. As shown in Fig. 4, the energy mapping relationship between the diffused light on the surface of the imaging object and the CCD detectors is as follows: where is the light flux density at the pin-hole with a value approximately equal to the light flux density ( )at the CCD planar; ( ) is the light flux density at the bin ; and are the area of and respectively. is the number of sampling points.

Optical-CT Dual-Modality Image Mapping
Optical-CT image registration based on DRR establishes the mapping relationship between the 2D optical image and the 3D CT data. So, the light intensity distribution on the surface of CT is completed following these steps: Step_1 Acquire 2D optical image and 3D CT data; Step_2 Set the initial projection angle after the 3D CT data being down-sampled; Step_3 Obtain the DRR virtual image under the setting angle from CT, and perform zero-filling to make it consistent with the same resolution of the optical image; Step_4 Multi-resolution spatial registration is performed on the DRR virtual image and optical image to obtain the optimal transform parameter and similarity measure value; Step_5 The optical image with external contour and the DRR virtual image are aligned iteratively until the termination condition is achieved by adjusting the optimal transformation parameters and similarity measure; Step_6 The 3D light intensity distribution on the surface of CT is obtained according to the coordinate correspondence obtained by the Eq. (12).

Experiments and Results
In order to verify the feasibility and effectiveness of the mapping algorithm proposed in this paper, we designed phantom and tumor mouse experiments. The optical images were acquired by the Animal Optical Imaging System (AOIS). White-light and fluorescence images in different orientations, -45º, 0º, and 45º, can be collected by rotating CCD, as shown in Fig. 5.

Phantom Experiment
A cloud-shaped homogeneous cylinder (Fig. 6) is used in the phantom experiment. The radius of each arcs is 10 mm and the height of the cylinder is 35 mm. The hollow cylinder with a 2.5 mm radius and 15 mm height is filled with Indocyanine Green (ICG). The volume data of the phantom with 15 marked points is 420 × 420 × 254, and each voxel is 0.15 × 0.15 × 0.15 mm. In order to verify the improved HNMI registration method, the registration centers are respectively set at -45º, 0º and 45º and the initial deviation angles of the DRR projection is set as -10º ~ 10º with step size of 1º. The optimal similarity measure values between the DRR virtual image and the optical image in different orientations are calculated, as shown in Fig. 7. Fig. 7(a) shows the variation of HNMI, and Fig. 7(b) shows the variation of NMI. Fig. 7 indicates that the HNMI has a good convergence, while the NMI exists multiple local extreme values. The 3D coordinates of these marked points on the optical images can be calculated after the registration of white-light image and DRR virtual image. The registration error is shown in Tab. 1. The ME is less than 1.87 mm and RMSE is less than 0.48 mm. The results is satisfied the requirements of experiments. To evaluate the performance of the mapping algorithm based on DRR optical-CT registration, the reconstruction results are also compared with the forward simulation results of COMSOL, as shown in Fig. 8. Fig. 8(a) is the surface light flux distribution of the COMSOL forward simulation, Fig. 8(b) is the surface light flux reconstruction based on the proposed mapping algorithm, in which there are 2011 nodes, 1228 triangular patches and 10331 tetrahedral elements. A total of 1014 surface sampling points are selected, the intensity distribution of these nodes is shown in Fig. 8(c). The orange curve is the result of the mapping algorithm and the blue curve is the COMSOL forward simulation result. The ME is 3.77% and the CF is 0.9854. Both light flux distributions have highly similarity and consistency.

Tumor Mouse Experiment
ICG was injected through the tail vein in the mouse experiment. Abdomen partial CT data of mouse was obtained from micro-CT system (Hiscan-2000, Suzhou Hiscan Information Technology Co., Ltd) with a scan parameter of 40 kVp, 200 μA. White-light images and fluorescent images were acquired from the AOIS system. The experimental scheme is the same with the phantom experiment.
The optimal similarity measure values between the DRR virtual images and the optical images in different orientations (-45º, 0º and 45º) are calculated, as shown in Fig. 9. Fig. 9(a) and Fig. 9(b) indicated the variation of HNMI and NMI respectively. The result shows that the HNMI is superior to the NMI in terms of convergency.  Based on the registration results of the optical-CT dual-modality images, 3D energy distribution of the biological surface is reconstructed in Fig. 11. Fig. 11(a) is the reconstruction result of the COMSOL forward simulation, and Fig. 11(b) is the reconstruction result based on the proposed algorithm, in which there are 2224 nodes, 8753 triangular patches and 10936 tetrahedral elements. A total of 970 sampling points on the surface of the mouse grid were selected. The distribution of the light intensity of these nodes is shown in Fig. 11(c). The orange curve is the result of the mapping algorithm and the blue curve is result of the COMSOL forward simulation. The ME is 3.14% and the CF is 0.9883.  To further verify the accuracy of the proposed registration algorithm, 3D bioluminescent source inside the mouse was reconstructed based on the surface light flux distribution. The fusion of the internal source and CT volume data is shown in Fig. 12.

Conclusion
This paper proposes an optical-CT dual-modality image mapping algorithm based on DRR registration. In the process of registration, the improved HNMI similarity measure is used to achieve spatial registration of optical image and DRR virtual image. The phantom and mouse experiments indicated that the proposed algorithm can align the different orientation 2D optical image to the 3D CT image and further reconstruct the 3D surface light flux distribution with high precision. This study is the basis for positioning and quantifing the 3D bioluminescent source inside the imaging object.