Optical-CT Dual-Modality Mapping Base on DRR Registration

Qingyang Zang,Dongsheng Li,Chunxiao Chen,＊ and Jianfei Li

1Department of Biomedical Engineering,Nanjing University of Aeronautics & Astronautics,Nanjing,211106,China.

Abstract:Optical-CT dual-modality imaging requires the mapping between 2D fluorescence images and 3D body surface light flux.In this paper,we proposed an optical-CT dual-modality image mapping algorithm based on the Digitally Reconstructed Radiography (DRR) registration.In the process of registration,a series of DRR images were computed from CT data using the ray casting algorithm.Then,the improved HMNI similarity strategy based on Hausdorff distance was used to complete the registration of the white-light optical images and DRR virtual images.According to the corresponding relationship obtained by the image registration and the Lambert’s cosine law based on the pin-hole imaging model,the 3D light intensity distribution on the surface of the object could be solved.The feasibility and effectiveness of the mapping algorithm are verified by the irregular phantom and mouse experiments.

Keyword:DRR; FMT; dual-modality registration; mapping

1 Introduction

Fluorescence Molecular Tomography (FMT) has been used frequently for the noninvasive study of cancer biological behavior [1-3].With specific imaging agents,FMT has powerful capability in detection and quantification of important biological processes by reconstructingin vivodistribution 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.

2 Theory and Method

2.1 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]:

where HT(p) is the intensity of the DRR virtual image at p,T is the 3D geometric transformation of the CT volume data,L(p,r) is the ray casting model,r is a set of parameter points on the light.

Figure1:DRR virtual image converted from CT

2.2 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.

Figure2:Flow chart based on DRR registration

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.

2.2.1 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 A and B,the NMI between them is defined as [15]:

where H(A) and H(B) represent the entropies of images A and B,respectively,and H(A,B)represents the joint entropy of images A and B,which are defined as follows:

where PA(a) and PB(b) respectively represent the gray-level distribution of images A and B.PA,B(a,b) is the joint probability distribution of gray level.

2.2.2 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 imagesAandBare 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]:

So,for imagesAandBto be registered,the improved HNMI function is defined as:

where w1 and w2 are constants used to balance the contributions of H(A′,B′) and NMI(A,B).After doing a lot of tests,w1 is set 0.2 and w2 is set 1 in the following phantom and mouse experiments.NMI(A,B) is the NMI of two images,and H(A′,B′) is the Hausdorff distance between point sets A′andB′.

2.2.3 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]:

where Gk(i,j) represents the kth layer image,i and j respectively are the row and column of the image,W(m,n) 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.

Figure3: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 m×n is the resolution of image,I(i,j) 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.

2.3 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 i,Pcal(i)is the registration coordinate and Pexp(i)is the real coordinate.N 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.

2.4 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 dPRis the light flux density at the pin-hole with a value approximately equal to the light flux density dP(rd)at the CCD planar; dP(r) is the light flux density at the bin r; Srand SRare the area of r and dR respectively.

Figure4:Light flux reconstruction of 3D surface based on pin-hole imaging model

ME and correlation factor (CF) are introduced as evaluation indexes for the light flux reconstruction.CF is defined as follows:

where Pcaland Pstdrepresent the light flux value at the sampling point by the mapping algorithm and the forward simulation method respectively.P is the mean,and σ(·) is the standard deviation of the surface light flux.N is the number of sampling points.

2.5 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_1Acquire 2D optical image and 3D CT data;

Step_2Set the initial projection angle after the 3D CT data being down-sampled;

Step_3Obtain 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_4Multi-resolution spatial registration is performed on the DRR virtual image and optical image to obtain the optimal transform parameter and similarity measure value;

Step_5The 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_6The 3D light intensity distribution on the surface of CT is obtained according to the coordinate correspondence obtained by the Eq.(12).

3 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.

Figure5:Image’s acquisition in different angles

3.1 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.

Figure6:Phantom experiment.(a) Phantom model; (b) Marked points on phantom surface

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 θzof 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.

Figure7:Phantom experiment in different orientations with initial deviation angle.(a) Variation of HNMI; (b)Variation of NMI

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.

Table1:Registration error (mm)

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.

Figure8:Surface light flux distribution of phantom.(a) Result of COMSOL forward simulation; (b)Result of proposed mapping algorithm; (c) Light flux distribution of each point on the surface

3.2 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.

Figure9:Mouse experiment in different orientations with initial deviation angle.(a) Variation of HNMI;(b) Variation of NMI

Fig.10 shows the optical-CT dual-modality image registration results of a mouse.Fig.10(a) is the white-light image,Fig.10(b) is the fusion image of fluorescence and white-light image,Fig.10(c) is the DRR virtual image computed from CT,and Fig.10(d) is the registration result.

Figure10:Optical-CT dual-modality images registration.(a) White-light image; (b) Optical fusion image;(c) DRR virtual image; (d) Registration result

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.

Figure11:Surface light flux reconstruction results of the mouse.(a) The result of COMSOL forward simulation; (b) The result based on the proposed algorithm; (c) Light flux comparison at points on the surface of the mouse

Figure12:The fusion of two modal images.(a) CT images; (b) Fusion images

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.

4 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.