3D reconstruction of light flux distribution on arbitrary surfaces from 2D multi-photographic images

Xueli Chen; Xinbo Gao; Duofang Chen; Xiaopeng Ma; Xiaohui Zhao; Man Shen; Xiangsi Li; Xiaochao Qu; Jimin Liang; Jorge Ripoll; Jie Tian

doi:10.1364/OE.18.019876

1. Introduction

As a promising molecular imaging technology, optical imaging (OI) has been attracting increased attention and been widely applied in many domains in recent years because of its significant advantages in temporal resolution, imaging contrast and sensibility, no ionizing radiation and cost-effectiveness [1–3]. In particular, optical tomography has become a valuable tool for the biomedical imaging field, since it is a noninvasive imaging technique that is capable of three-dimensionally recovering the location and concentration of the luminescent probe inside a small living animal from surface measurements of the transmitted light flux [4–6]. Among the various optical tomography modalities, fluorescence molecular tomography [6] and bioluminescence tomography [3,7] have developed rapidly and become research focuses over the past years because of their excellent performance, which can recover the spatial distribution of the concentration of targeted cells for tracking tumor cells, stem cells, immune cells and bacteria as well as for imaging gene expression [2,5–7].

In contrast to planar imaging, optical tomography can recover the three-dimensional (3D) spatial distribution information of luminescent probes inside the animals, specifically their location and emission powers, from measurements taken on the surface of the body. The measured surface signals are incorporated into a tomographic procedure written as a linear system of equations which are solved for the unknown, 3D spatial distribution of the internal luminescent probes [8–15]. In order to measure the light flux traversing the surface, the ideal experimental setup would be such that the detectors are not in contact with the body surface. Such non-contact detection system offers simple experimental procedures and enables multi-view tomographic projections. Some algorithms have been developed to model free-space propagation of diffuse photon flux from a body surface with arbitrary shape to a set of non-contact detectors [16–19]. Uniform frameworks have also been developed to recover the internal target directly using two-dimensional (2D) photographic images [20,21]. In addition, a simple method has been proposed to establish a relationship of diffuse photon flux between the non-contact detectors and the body surface in the case of small numerical aperture measurement [22], which is under the assumption that a one-to-one mapping exists between the non-contact detectors and the body surface. However, an effective and practical method is lacking that reconstructs photon flux on a body surface of arbitrary shape from the measurements on a set of non-contact detectors, which is a procedure necessary in practical applications of optical tomography.

In this paper, we propose a novel and effective method that reconstructs photon flux on a body surface of arbitrary shape accurately from multi-view 2D photographic images measured by a CCD camera. Taking an example from the principle of optical path reversibility and the ideology of reciprocity theorem between sources and detectors, this method assumes that each pixel unit of the photographic images captured by the CCD camera is viewed as a Lambertian source. As a result, an analytical expression which accounts for this type of free-space light propagation is addressed from radiosity theory based on the Lambertian source approximation [23] and integrated with the radiance theorem described in literature [19]. Particularly, the introduction of the virtual detector plane and the lens-dependent visibility factor contributes to construct a specific relationship between the non-contact detectors and the body surface, which facilitates to eliminate the influence of overlapping adjacent projection images. The performance of the proposed method is demonstrated with comparison experiments involving both regularly and mouse-shaped phantoms. Finally, the potential of this approach to improve the current status of optical tomography is further discussed.

2. Theory and method

2.1 Lambertian source theory

In order to model the angular distribution of diffuse light on exiting a diffusive medium, a common theory employed is to assume a Lambertian source distribution of light when traveling in free-space [24]. A Lambertian source is such that irradiates isotropically into the surrounding space. Following this approximation, the differential power $d P (r_{d})$ measured by a detector at position $r_{d}$ of differential area and orientation $d A_{d}$ emitted from a source $r_{s}$ of differential area and orientation $d A_{s}$ , in direction s within a solid angle $d Ω$ is [24,25]:

d P (r_{d}) = \frac{1}{π} J (r_{s}) \frac{\cos θ_{s} \cos θ_{d}}{{| r_{d} - r_{s} |}^{2}} d A_{d} d A_{s},

where

J (r_{s})

is the flux density at the source point

r_{s}

;

\cos θ_{s} = (d A_{s} • s) / d A_{s}

is the cosine dependence of the Lambert’s law;

\cos θ_{d} = - (d A_{d} • s) / d A_{d}

accounts for detector orientation with respect to the direction s that always points to

r_{d}

;

d A_{s}

and

d A_{d}

are differential areas of the source and detectors, respectively. Equation (1) distinctively depicts the transport characteristic of diffuse light in free-space emitted from a Lambertian source, as the case shown in Fig. 1 .

Fig. 1 Schematic diagram of the transport characteristic of diffuse light in free-space.

Download Full Size | PDF

2.2 Surface flux reconstruction algorithm

Reconstruction of 3D surface flux distribution is achieved by employing a novel and practical theoretic expression to depict the corresponding relationship between the 2D flux density map captured by a non-contact CCD camera and 3D flux distribution on the body surface, and by capitalizing on computed tomography techniques to accurately register the surface of arbitrary geometries. The theoretic expression employed here is derived under the assumptions that the surface of arbitrary geometries and the corresponding multi-view photographic images at the CCD camera have been obtained. The surface of arbitrary geometries is reconstructed using an X-ray computed tomography (CT) system employing the FDK algorithm, a reconstruction algorithm for 3D cone beam CT developed by Feldkamp, Davis and Kress [26]. In particular, multi-view 2D photographic images should be registered onto the coordinate system of the reconstructed surface before the reconstruction of surface flux distribution is performed. Moreover, taking examples from the principle of optical path reversibility and the reciprocity theorem, we assume that each pixel unit in the photographic images registered at the CCD camera acts as a Lambertian source irradiating the body’s surface. Thus, corresponding relationship between the 2D flux density map at the CCD camera and the 3D flux distribution on the body surface satisfies the transport characteristic of a Lambertian source. In addition, the concept of a virtual detector plane [16], which would be the focal plane of the objective of the CCD camera, is introduced as a transit station in the process of surface flux reconstruction to relate the CCD camera’s measurement and the reconstructed arbitrary surface.

Figure 2 illustrates the reconstruction procedure of 3D surface flux distribution. Firstly, the non-contact detection system is modeled as a simple optical system which is comprised of two components: a simplified thin-lens and a CCD camera plane, using the simplification theory of an imaging objective [19], where the focal plane of the CCD camera plane is selected as the virtual detector plane [27]. Secondly, light flux at the CCD camera plane is mapped onto the corresponding virtual detector plane based on the radiance theorem which states that radiance is conserved through a loss-less optical system [19]. Thirdly, the corresponding relationship of points between the virtual detector plane and the body surface is established with the help of a visibility factor which is relevant to the simplified thin-lens. Merging with the visibility factor, the flux at the surface of the body is finally determined based on the Lambertian source theory presented in Eq. (1). In conclusion, given the multi-view 2D photographic images, the reconstructed surface of arbitrary geometries and the imaging information of non-contact detection system, the flux density of any point at the surface of the body can be calculated using the following equation:

J (r) = \frac{1}{π} \int_{Ω} ℜ B (r_{d}) T (r_{d}, r) d Ω,

where Ω is detection space constructed by the CCD camera;

J (r)

is flux density of point r at the surface of the body;

B (r_{d})

is gray-scale value of point

r_{d}

in the detection space Ω; ℜ is a conversion factor that converts the gray-scale value into the corresponding flux based on the intrinsic properties of the camera:

ℜ = Θ Q^{- 1} E_{0} t_{p}^{- 1},

where Θ is a specific value relevant to CCD camera and defined as a ratio of electron number corresponding to full well over maximum value of gray-scale in the photographic image; Q is quantum efficiency of the CCD camera;

t_{p}

is the exposure time for acquiring each image;

E_{0}

is the unit photon energy and can be calculated as follows:

E_{0} = \frac{h c}{λ},

where h is the Planck constant; c and λ are the velocity and wavelength of light registered at the CCD camera.

Fig. 2 Diagram for the procedure of the reconstruction of 3D surface flux distribution.

Download Full Size | PDF

In Eq. (2), a transfer function that accounts for the relationship of energy transformation from the plane of the CCD camera to the body surface has been introduced:

T (r_{d}, r) = τ (f) ξ (r_{d}, r) \frac{1}{{| r - r_{d} - \frac{t u^{2}}{f \cos θ} s (r_{d}, r_{d}^{*}) |}^{2}} \cos θ_{s} \cos θ_{d},

where

τ (f)

is the optical transmittance of the camera lens;

ξ (r_{d}, r)

is a visibility factor that discards the detection points not visible from the surface and depends on the camera lens; t is magnification ratio of the non-contact detection system; u denotes an object distance between the virtual detector plane and the simplified thin-lens; f is the focal length of the camera lens;

r_{d}^{*}

is a conjugated point of

r_{d}

at the virtual detector plane and satisfies the lens law with

r_{d}

;

s (r_{d}, r_{d}^{*})

is a unit vector which denotes the direction from

r_{d}

to

r_{d}^{*}

; θ is an included angle between the aforementioned unit direction vector and the optical axis;

θ_{s}

and

θ_{d}

are included angles between the unit direction vector from

r_{d}^{*}

to r and the normal vector of the surface or virtual detector point, respectively. The physical significance of the parameters is shown in Fig. 2. It should be noted that it is distinctly important to accurately determine the position and size of the virtual detector plane for achieving the most approximate reconstructed flux distribution, which can be determined according to the simplified theory of the camera lens presented in literature [19].

Algorithm 1 Implementation flowchart for light flux reconstruction on arbitrary surfaces

Require: Optical measurements, 3D surface of subject and information of optical system. 1: Preprocess the optical measurements and co-register with subject’s surface. 2: Discretize the subject’s surface and obtain the number of differential surface unit M . 3: Discretize the CCD camera plane and obtain the total number N . 4: for i = 1 to M do 5: Compute the coordinates of point

r^{i}

on the subject’s surface. 6: Let

k = 0

,

J (r^{i}) = 0

. 7: for j = 1 to N do 8: Compute the coordinates of point

r_{d}^{j}

at the CCD camera plane. 9: Compute the visibility factor

ξ (r_{d}^{j}, r^{i})

. 10: if

ξ (r_{d}^{j}, r^{i}) = 1

then 11: Compute

J^{k + 1} (r^{i})

using the discretized form of Eqs. (2)-(5). 12: Compute

J (r^{i})

by

J (r^{i}) = J (r^{i}) + J^{k + 1} (r^{i})

. 13: Update k by

k = k + 1

. 14: else 15: Break the current inner loop. 16: end if 17: end for 18: end for

Equation (2) is a generalized formula for depicting the reconstructed flux density of a surface point from multi-view 2D photographic images registered at the CCD camera. Integrating Eq. (2) over the whole body surface, we obtain the following analytical expression for accurately determining the reconstructed flux distribution on the whole body surface:

J = \frac{1}{π} \iint ℜ B (r_{d}) T (r_{d}, r) d Ω d S .

To solve the above integral equation, discretization of the whole body surface and of the CCD camera plane is needed. Thus, the analytical expression presented in Eq. (6) is rewritten as the following matrix form:

J = AB,

where the components of the matrix A and the vectors

J, B

are defined as:

{\begin{cases} J_{i} = J (r^{i}) \\ A_{i j} = \frac{1}{π} ℜ T (r_{d}^{j}, r^{i}) \\ B_{j} = B (r_{d}^{j}) \end{cases},

where

1 \leq i \leq M

and

1 \leq j \leq N

; M is the number of discrete elements on the body surface and N is the number of discrete elements at the CCD camera. Solving Eq. (7), we can obtain the 3D flux distribution on the whole body surface.

To sum up, the implementation flowchart for the proposed algorithm is decomposed into pseudo-code form and presented in Algorithm 1.

3. Experiments and results

To evaluate the performance of the proposed method, comparison experiments took place, including the measurements of regularly and mouse-shaped phantoms in order to recover the 3D distribution of the surface flux. Furthermore, an in vivo small animal imaging experiment was also performed to illustrate the potential of the proposed method in the application of optical tomography. For experimental measurements, multi-view 2D photographic images and anatomic structure of the subjects under study were acquired by a dual-modality ZKKS-Direct 3D molecular imaging system (jointly developed by Guangzhou Zhongke Kaisheng Medical Technology CO., Ltd and Xidian University), as shown in Fig. 3 . The OI system consists of a liquid-cooled back-illuminated CCD camera (Princeton Instruments/Acton 2048B) coupled to a Micro camera lens (Micro-Nikkor 55 mm, f/2.8) to acquire the multi-view 2D photographic images of light emitted from the surface of subjects under study. The micro CT system integrates an X-ray tube (Oxford Instruments, series 5000) and a flat panel sensor (Hamamatsu, C7921CA-02) to obtain high-quality 3D anatomic structure using an accelerated FDK reconstruction algorithm [28]. To avoid anatomical deformation of the subjects during the experiments, the two systems are orthogonally mounted around the imaging holder that is fixed onto a computer controlled rotation stage.

Fig. 3 Dual-modality OI/micro CT system. (1) CCD camera, (2) Camera lens, (3) X-ray tube, (4) X-ray detector, (5) Lifting tables, (6) Rotation stage, (7) Translation stages, (8) Mouse holder.

Download Full Size | PDF

3.1 Experiment verifications based on phantoms

Two groups of imaging experiments were carried out to validate the accuracy of the proposed algorithm in this subsection. In the experiments, two different kinds of surfaces were employed, namely a cubic or cylindrical phantom (the ‘regular’ surface from now on) and a mouse-shaped phantom (termed the ‘irregular’ surface from now on). All of the phantoms utilized in the experiments were made of nylon. To simulate the luminescent light source, a fluorescent solution extracted from a red luminescent light stick (Glowproducts) was injected into the phantom. The peak wavelength of the luminescence light generated by this luminescent solution is approximately 650 nm, whereas the optical properties of the phantom used were measured at 660 nm by a time-correlated single photon counting system [29]. The absorption coefficient was found to be $μ_{a} = 0.0138 m m^{- 1}$ and the reduced scattering coefficient was $μ_{s}^{'} = 0.91 m m^{- 1}$ . To validate the reconstructed flux distribution on the phantom surface, a software platform of Molecular Optical Simulation Environment (MOSE) was also employed to obtain the flux distribution on the phantom surface and its simulation results were compared with the reconstructed results of the proposed method. MOSE, which is based on the Monte Carlo method to simulate light propagation in turbid media, is currently under development and regularly being updated. This software can be freely accessed and downloaded on the web site: http://www.mosetm.net. Sufficient experiments have demonstrated that the light flux on the phantom surface simulated by MOSE can served as a standard to evaluate the accuracy of the reconstruction results [30].

3.1.1 Regular surface flux reconstruction

Two different phantoms were designed to evaluate the performance of the proposed method in the case of regular surface, including a cubic phantom and a cylindrical phantom. In order to better scale the discrepancy between the reconstructed and simulated flux, the mean error (ME) and the correlation factor (CF) are introduced:

\begin{matrix} \bar{e} = \sum_{i = 1}^{N} | P^{r e c (i)} - P^{s i m (i)} | / N, and \\ ρ = \sum_{i = 1}^{N} (P^{r e c (i)} - {\bar{P}}^{r e c}) (P^{s i m (i)} - {\bar{P}}^{s i m}) / ((N - 1) σ (P^{r e c}) σ (P^{s i m})), \end{matrix}

where

\bar{e}

and ρ are the value of ME and CF, respectively; the superscript rec represents the flux reconstructed by the proposed method and sim the simulated flux of MOSE;

\bar{P}

and σ are the mean value and standard deviation of the flux that can be either

P^{r e c}

or

P^{s i m}

; N is the total sample number on the phantom surface and i represents the ith sample. The CF indicates the degree of correlation between the reconstructed and simulated flux while the ME describes the discrepancy of them. Accordingly, the closer CF gets to unity, and the closer ME gets to zero, the better the performance of the proposed method is.

3.1.1.1 Cubic phantom experiment

In the cubic phantom experiment, a homogenous phantom of $45 m m$ in lateral dimension was utilized. One small hole of $2.95 m m$ radius and $25.2 m m$ depth from the top surface was drilled at the center of the phantom in order to place the flourophore, as shown in Fig. 4 . Luminescent solutions of $150 μ l$ volume were injected into the small hole, which forms a light source of $5.4 m m$ height. After implanting the light source, the rest of the small hole was filled with a solid piece of the same material so as to avoid light leakage. After that, the cubic phantom was mounted on a 360° rotation stage. By rotating the phantom a complete 360°, four views of 2D photographic images were acquired for the surface flux recovery while the anatomical structure of the cubic phantom was obtained by the aforementioned dual-modality OI/micro CT system. Since the light source is located at the center of the phantom, the four views of the photographic images yielded the same light flux distribution. Therefore, one of the images, 90° view in this experiment, was selected to reconstruct the flux distribution, as shown in Fig. 5(a) . In this figure, the phantom holder mounted on the rotation stage can be seen at the bottom of the image.

Fig. 4 Experimental setup for the cubic phantom based comparison experiment. (a) Picture of the cubic phantom used in the optical imaging experiment; (b) Numerical phantom used in the simulation.

Download Full Size | PDF

Fig. 5 Comparison results between the reconstructed flux distribution and the measured or simulated results. (a) 90° view photographic image captured by the CCD camera; (b) Simulated results of MOSE; (c) Reconstructed flux distribution of the proposed method.

Download Full Size | PDF

Reconstruction of flux distribution on the selected surface of the phantom is shown in Fig. 5 (c). Figure 5(b) depicts the corresponding simulated result of MOSE. Both the reconstruction and simulation were based on the same density of mesh which consisted of 12288 triangular elements and 35937 nodes. Similar tendency and good agreement were observed in this case. In particular, the reconstructed flux distribution performed well in accordance with the photographic image acquired by the OI system due to the flat characteristics of the phantom surface. As an example, we are also able to recover the position of phantom holder previously mentioned (see the bottom of Fig. 5(c)). However, there was a slight difference between the reconstructed and simulated flux distribution, which might be caused by the non-uniformity of the cubic phantom. In the simulation of MOSE, we made an assumption that the experimental setup was absolutely ideal, assuming a uniformity of the optical properties and the accurate position and size of the light source. A comparison of the curves of the results reconstructed by the proposed method and simulated by MOSE or measured by CCD camera is shown in Fig. 6 . Figure 6(a)-(c) present the comparisons between the reconstructed flux and the measured data and Fig. 6(d)-(f) show the comparisons between the reconstructed and simulated flux. In these figures, the solid lines represent the simulated results or the measured data while the asterisks show the reconstructed flux of the proposed method. From Fig. 6, we also find that the reconstructed flux was more consistent with the measured data than the simulated results. Although the reconstructed flux exhibits slightly less intensity around the peak of the curves in Fig. 6(e) and (f), the ME and CF still show the conformity between the reconstructed flux and the simulated results, as listed in Table 1 . Comparison of the results indicated that the proposed method performed very well in the reconstruction of flux distribution on the flat surface in the case of the cubic phantom.

Fig. 6 Compared curves between the reconstructed flux and the measured or simulated results at height 0, 5 and 10 mm from the center of phantom. (a)-(c) Comparisons of the reconstructed and measured flux, (d)-(f) Comparisons of the reconstructed and simulated flux; (a) and (d) 0 mm, (b) and (e) 5 mm, (c) and (f) 10 mm.

Download Full Size | PDF

Table 1. Error comparisons between the reconstructed and the measured or simulated results

View Table | View all tables in this article

3.1.1.2 Cylindrical phantom experiment

In the cylindrical phantom experiment, a homogeneous phantom of $27 m m$ in height and $15 m m$ in radius was used to conduct the comparison. Similarly to the case presented previously, one small hole of $1 m m$ radius and $16 m m$ depth was drilled in the phantom to locate the fluorescent solution, as shown in Fig. 7(a) . A luminescent light source of $4 m m$ height and $1 m m$ radius was embedded into the phantom with its center at (9, 0, 0), as shown in Fig. 7(c). After the cylinder phantom was mounted on the rotation stage and rotated a complete 360°, four equally spaced views 2D photographic images and the corresponding anatomical structure were acquired utilizing the aforementioned dual-modality OI/micro CT system. Figure 7(b) presents the diagram of four different perspectives and the corresponding measured four-view 2D photographic images are shown in Fig. 8(a)-(d) . Using these measurements, a 3D flux distribution on the phantom surface was reconstructed using the proposed method.

Fig. 7 Experimental setup for the cylindrical phantom based comparison experiment. (a) Picture of the cylindrical phantom used in the imaging experiment; (b) Four different views for capturing 2D photographic images; (c) Numerical phantom used in the simulation.

Download Full Size | PDF

Fig. 8 Four-view 2D photographic images measured at the CCD camera. (a), (b), (c) and (d) are front-, left-, back- and right-view, respectively.

Download Full Size | PDF

To examine how the reconstructed flux distribution was affected by the number of images used in the reconstruction, we reconstructed the flux distribution on the phantom surface using one, two, three and four-view images, respectively. Comparisons between the reconstructed and simulated flux distribution are shown in Fig. 9(a)-(l) . Similarly, all of the results presented in Fig. 9 were based on the same density of mesh that consisted of 259200 elements. For each comparison, results for three detector z positions: z_d = 0.0 mm, z_d = 5.0 mm and z_d = 10.0 mm are examined being the recovered flux normalized to their maximum value versus the complete 360° projections. The ME and CF between the reconstructed and simulated results are listed in Table 2 . Similar tendency and good agreement between them are observed in the cases examined. From Fig. 9, some interesting results are summarized. Firstly, we find that the reconstructed flux seems slightly sparse around the peak value of the curves, which are intrinsically caused by the inadequate discretization of phantom surface. If the discretization of the phantom surface is refined sufficiently, the reconstructed flux distribution would be smoother. Figure 10(a)-(c) present the fact that the smoothness of the results changed better with the increase of mesh density. Secondly, the reconstructed flux is slightly lower in intensity than the simulated results of MOSE around the valley of the curves. This is mainly caused by the insufficient measurements (projections) in the experiment. If more photographic images were acquired at different projections during the experiment and included in the reconstruction, the more accurate the reconstructed surface flux would be. Thirdly, the reconstructed flux gets to be closer to the simulated results as the image number increases, as shown in Fig. 9. In order to compare the effect of the number of projections used, Fig. 10(d) shows the logarithm of the reconstructed and the simulated flux at the central height of the phantom, which shows the improvement of the reconstructed flux more clearly as the number of projections used increases. A similar variation trend also can be seen in Table 2. By increasing the image number, the ME turns to be smaller and the CF gets to be closer to unity. Lastly, the errors are tolerable for practical applications in all the cases examined. From the error comparisons listed in Table 2, we conclude that one single photographic view is sufficient to reconstruct the surface flux distribution in this experiment, with the ME less than 0.03. As a result, an interesting and useful conclusion can be addressed that a practical reconstructed flux distribution can be achieved using one or some of acquired view images if the intensities of these images are significantly larger than others, such as the case presented in this experiment. As expected, the proposed method worked well in the reconstruction of flux distribution on the curved surface in the case of the cylindrical phantom, using different numbers of photographic images.

Fig. 9 Comparison results between the reconstructed and the simulated flux at height 0, 5 and 10 mm from the center of the phantom surface. (a)-(c) Reconstruction using one photographic image shown in Fig. 7(a), (d)-(f) Reconstruction using two adjacent views photographic images shown in Fig. 7(a) and (b), (g)-(i) Reconstruction using three photographic images shown in Fig. 7(a), (b) and (c), (j)-(l) Reconstruction using four views photographic images shown in Fig. 7(a)-(d); (a), (d), (g) and (j) 0 mm; (b), (e), (h) and (k) 5 mm; (c), (f), (i) and (l) 10 mm.

Download Full Size | PDF

Table 2. Error comparisons between the reconstructed and the simulated results

View Table | View all tables in this article

Fig. 10 Reconstructed flux at height 0 mm from the center of phantom surface. (a)-(c) Reconstruction using different density of surface meshes, including 48600, 194400 and 388800 elements, respectively; (d) Logarithmic curve for the reconstruction using various photographic views images compared with the simulated flux.

Download Full Size | PDF

3.1.2 Irregular surface flux reconstruction

To demonstrate the capability of the proposed method in the reconstruction of 3D surface flux distribution in the case of irregular surface, a mouse shaped phantom was designed and used to conduct the comparison experiment, as shown in Fig. 12(a) . Similar to the two cases previously described, a small hole of $1.5 m m$ in radius was drilled at the head of the mouse phantom to locate the light source where luminescent solutions of $20 μ l$ volume were injected. The dimensions of the light source were $3 m m$ in diameter and $2.8 m m$ in height with its center at (20.80, 19.52, 20.80) mm, dimensions which were obtained by the micro-CT system. After the mouse shape phantom was mounted on the rotation stage and rotated a complete 360°, the photographic images from the ventral, left lateral, dorsal and right lateral views of the mouse phantom and the corresponding CT images were acquired using our dual-modality OI/micro CT system. Using the software Amira^TM (Mercury Computer Systems, MA), the surface structure of the mouse phantom used in the reconstruction was extracted from the anatomical structure which was reconstructed from the CT images using the FDK algorithm. Figure 11(a)-(d) shows the four views of the photographic images normalized to their maximum intensity. In addition, four acquisition views of the mouse phantom were also acquired to be used for registering the photographic images onto the coordinate system of the extracted surface structure.

Fig. 12 Experimental setup and results for the mouse phantom based comparison experiment. (a) Picture of the mouse phantom used in the imaging experiment; (b) Simulated flux distribution of MOSE; (c) Reconstructed flux distribution of the proposed method; (d)-(f) Compared curves between the reconstructed and the simulated flux at height 13, 17 and 21 mm, respectively.

Download Full Size | PDF

Fig. 11 Four-view 2D photographic images measured at the CCD camera. (a) Ventral view; (b) Left lateral view; (c) Dorsal view; (d) Right lateral view.

Download Full Size | PDF

Using the extracted surface structure and four views that have been registered to the coordinate system of the surface structure, we reconstructed the surface flux distribution utilizing the method proposed in this work and obtained the simulated flux using the platform of MOSE, as shown in Fig. 12. Figure 12(a) describes the mouse shaped phantom adopted in the experiment; Fig. 12(b) shows the simulated flux distribution and the corresponding calculated results are shown in Fig. 12(c). Both the reconstructed and simulated fluxes were based on the same density of mesh that consisted of 50000 triangular elements and 25002 nodes. Moreover, results for three different detector z positions: z = 13 mm, z = 17 mm, and z = 21 mm were examined. Figure 12(d)-(f) shows plots of the light flux normalized to the maximum value versus the complete 360° detection positions. The solid lines show the simulated flux of MOSE and the asterisks represent the reconstructed flux of the proposed method. From Fig. 12, a similar tendency and comparable distribution between the simulated and reconstructed flux are observed, with the average ME and CF being about 0.0537 and 0.9034. However, a bigger discrepancy between them can also be found in this experiment. This might be caused by several reasons. Firstly, we employed a mouse shape phantom with arbitrary geometries and irregular surface in this experiment so that the surface structure was more complicated than that in the regular surface experiment. The complicated and irregular surface induced the inaccuracy of the registering results between the 2D photographic images and the surface structure. Poor registering results would lead to poor reconstruction of the surface flux distribution. Secondly, the non-uniformity of the phantom material performed greater influence on the reconstruction of the surface flux distribution in this experiment due to the irregularity and complexity of phantom surface. Lastly, the information loss in the acquisition of photographic images led to the reconstructed flux distribution smaller and more convergent than the simulated results, as shown in Fig. 12(b) and (c). If more photographic images were acquired and employed in the reconstruction, the reconstructed flux distribution would be significantly improved.

In order to further illustrate the effectiveness of the reconstructed surface flux distribution, we recovered the light source using both the reconstructed and the simulated flux distributions, respectively. Except the surface flux distribution, the other parameters for reconstruction were identical in both cases. In order to retrieve the position and intensity of the light source, the phantom domain was discretized to a tetrahedral mesh that consisted of 18341 tetrahedral elements and 3916 nodes and a potential permissible source region was set as ${(x, y, z) | 18 < x < 24, 10 < y < 25, 19 < z < 22}$ according to the flux distribution on the phantom surface. Employing the self-developed adaptive hp finite element method (hp-FEM) [31], the source positions were accurately recovered. Figure 13 presents the recovered results of source positions, where Fig. 13(a) gives the recovered light source using the simulated flux and the recovered result using the reconstructed flux is shown in Fig. 13(b). The red cylinders represent the actual sources determined by the CT images and the small domains beside the cylinders express the recovered sources. From Fig. 13, we find that the recovered sources approach the actual ones well in both cases. In the case of using the reconstructed surface flux distribution, the recovered source position is (19.66, 20.30, 20.93), which is 1.39 mm away from the actual source, i.e. in the order of one transport mean free path. For the simulated results, the recovered source position is (20.01 20.29 21.86), with a distance of 1.53 mm from the actual one, which is a bit worse than the result using the reconstructed surface flux. As a result, the reconstructed flux is effective and can be comparable to the simulated flux of MOSE in the light source recovery.

Fig. 13 Comparison results of the internal light source reconstruction. (a) Reconstruction using the simulated surface flux of MOSE; (b) Reconstruction using the reconstructed surface flux of the proposed method.

Download Full Size | PDF

The comparison results presented in this experiment evinced that the proposed method is capable of reconstructing the flux distribution on a complicated and irregular surface and overall exhibits excellent performance.

3.2 In vivo imaging experiment

To further study the potential of the proposed method in the application of optical tomography, an in vivo imaging experiment of a living mouse with an implanted light source was conducted. This experiment was designed to examine the ability to recover the internal light source in a heterogeneous mouse model using a reconstructed surface flux distribution. In the experiment, a living mouse was used to construct the heterogeneous mouse model. A luminescent catheter, which was filled with the fluorescent solution mentioned previously was used as the internal light source. The dimensions of the catheter were $1.4 m m$ in diameter and $4.5 m m$ in length. Prior to performing the experiment, the mouse was anesthetized and the catheter was implanted into the epigastric torso of the mouse. Photographic images of the light flux distribution and reference white-light images were both acquired from the ventral, left lateral, dorsal and right lateral of the living mouse by the aforementioned dual-modality OI/micro CT system. The exposure time was set as 1.5 min. Based on the micro-CT images acquired from the complete 360° projections, the surface and volumetric mesh utilized in the reconstruction were generated using Amira. They contained 4,000 triangle meshes, 4,975 nodes and 27,092 tetrahedral elements. Under the precondition that the photographs had been registered to the coordinate system of the triangle meshes, which was realized by registering 2D planar images to 3D space of CT images based on some labeled markers [32], the flux distribution on the mouse surface was reconstructed using the proposed method and presented in Fig. 14 . From the micro-CT images, the catheter source can be clearly identified, as shown in Fig. 15 . Using the reconstructed surface flux distribution, we recovered the internal catheter source based on the heterogeneous mouse model utilizing the hp-FEM. In order to construct the heterogeneous mouse model, several main organs, such as adipose, heart, lungs, liver and kidneys were segmented from micro-CT images using a combination method of threshold and interactive segmentation, and the corresponding optical properties were calculated according to the literature [33] and listed in Table 3 . The recovered position of the catheter source is shown in Fig. 15, where Fig. 15(a) presents the 3D view of the recovered result and the corresponding coronal, axial and sagittal views are shown in Fig. 15(b)-(d). The actual central position of the internal catheter source measured via CT reconstruction is at (29.28, 22.28, 26.64) mm and the corresponding recovered central position is (30.136, 21.7582, 26.3759) mm, with a distance error being 1.0367 mm. This study experiment showed that the proposed method provides a reconstructed surface flux distribution which enables accurately recovering the internal light source in an optical tomography application. In the present study, the light source intensity and the instrument factors were not calibrated; however, the possibility to converts them to absolute concentration units is currently pursued by integrating a calibration procedure into our non-contact detection system for which work is ongoing.

Fig. 14 Reconstructed surface flux of living mouse based in vivo imaging experiment.

Download Full Size | PDF

Fig. 15 Recovered internal catheter source of living mouse based in vivo imaging experiment using the reconstructed surface flux distribution. (a) 3D view of the recovered internal source result; (b)-(d) the corresponding coronal, axial and sagittal views, respectively. The position of the recovered source is marked with red circle.

Download Full Size | PDF

Table 3. Optical properties for the organs of mouse

View Table | View all tables in this article

4. Discussion and conclusion

We have developed a novel and practical method for 3D flux reconstruction of arbitrarily shaped surface from multi-view 2D photographic images. Reconstructed results for both regular and irregular surfaces have been illustrated. Comparisons with simulated results obtained from the platform of Molecular Optical Simulation Environment (MOSE) have demonstrated the accuracy and effectiveness of the proposed method. Furthermore, the potential of the method in its application on optical tomography is illustrated by a qualitative study of the source reconstruction experiment as well, in which the position of the internal implanted source is accurately recovered based on the reconstructed surface flux distribution using appropriate algorithms. However, there are several deficiencies required that need to be improved. Firstly, the surface of the subject under study has to be discretized adequately to obtain a smooth reconstructed flux distribution, which will lead to a high computational cost. Secondly, use of more projection views from the photographic images employed in the reconstruction will improve the quality of the reconstructed surface flux distribution. However, by including additional projections superposition between two adjacent images must be accounted for. Lastly, quantitative reconstruction (for example for bioluminescent sources) is still lacking in this approach both for the reconstruction of the surface flux and for the internal light source, since they are obviously directly related.

Overall, we believed that this approach will expose the complete projection imaging and the tomographic potential of optical tomography as applied to in vivo small animal imaging. These preliminary results show that the proposed method is of great potential and will prove crucial in the field of tomographic imaging if quantitation and high resolution is pursued. Further study will concentrate on the acceleration of the proposed method, the elimination of superposition problem caused by an increase in the number of images and the quantitative reconstruction of surface flux distribution.

Acknowledgments

This work is supported by the Program of the National Basic Research and Development Program of China (973) under Grant No.2006CB705700, the Cheung Kong Scholars and Innovative Research Team in University (PCSIRT) under Grant No.IRT0645, the Chair Professors of Cheung Kong Scholars Program of Ministry of Education of China, CAS Hundred Talents Program, the National Natural Science Foundation of China (NNSFC) under Grant Nos. 30873462, 30900334, 60771068, the Shaanxi Provincial Natural Science Foundation Research Project under Grant No. 2009JQ8018, and the Fundamental Research Funds for the Central Universities. J. Ripoll acknowledges support from EU FP7 grant FMT-XCT and the Bill and Melinda Gates Foundation TB Drug Accelerator Program.

References and links

1. R. Weissleder and V. Ntziachristos, “Shedding light onto live molecular targets,” Nat. Med. 9(1), 123–128 (2003). [CrossRef] [PubMed]

2. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005). [CrossRef] [PubMed]

3. J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27(5), 48–57 (2008). [CrossRef] [PubMed]

4. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50(4), R1–R43 (2005). [CrossRef] [PubMed]

5. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]

6. V. Ntziachristos, C. H. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. 8(7), 757–761 (2002). [CrossRef] [PubMed]

7. G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. McLennan, and M. Henry, “In vivo mouse studies with bioluminescence tomography,” Opt. Express 14(17), 7801–7809 (2006), http://www.opticsinfobase.org/as/viewmedia.cfm?id=97670&seq=0. [CrossRef] [PubMed]

8. A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004), http://www.opticsinfobase.org/ol/ViewMedia.cfm?id=81637&seq=0. [CrossRef] [PubMed]

9. X. Song, D. Wang, N. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007), http://www.opticsinfobase.org/VJBO/viewmedia.cfm?uri=oe-15-26-18300&seq=0. [CrossRef] [PubMed]

10. W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. V. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13(18), 6756–6771 (2005), http://www.opticsinfobase.org/jdt/viewmedia.cfm?id=85344&seq=0. [CrossRef] [PubMed]

11. Y. J. Lv, J. Tian, W. X. Cong, G. Wang, J. Luo, W. Yang, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express 14(18), 8211–8223 (2006), http://www.opticsinfobase.org/jot/viewmedia.cfm?id=97939&seq=0. [CrossRef] [PubMed]

12. J. C. Feng, K. B. Jia, G. R. Yan, S. P. Zhu, C. H. Qin, Y. J. Lv, and J. Tian, “An optimal permissible source region strategy for multispectral bioluminescence tomography,” Opt. Express 16(20), 15640–15654 (2008), http://www.opticsinfobase.org/oe/viewmedia.cfm?uri=oe-16-20-15640&seq=0. [CrossRef] [PubMed]

13. H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,” Opt. Lett. 31(3), 365–367 (2006). [CrossRef] [PubMed]

14. H. Dehghani, S. C. Davis, and B. W. Pogue, “Spectrally resolved bioluminescence tomography using the reciprocity approach,” Med. Phys. 35(11), 4863–4871 (2008). [CrossRef] [PubMed]

15. A. D. Klose, B. J. Beattie, H. Dehghani, L. Vider, C. Le, V. Ponomarev, and R. Blasberg, “In vivo bioluminescence tomography with a blocking-off finite-difference SP3 method and MRI/CT coregistration,” Med. Phys. 37(1), 329–338 (2010). [CrossRef] [PubMed]

16. J. Ripoll, R. B. Schulz, and V. Ntziachristos, “Free-space propagation of diffuse light: theory and experiments,” Phys. Rev. Lett. 91(10), 103901 (2003). [CrossRef] [PubMed]

17. J. Ripoll and V. Ntziachristos, “Imaging scattering media from a distance: theory and applications of noncontact optical tomography,” Mod. Phys. Lett. B 18(28 & 29), 1403–1431 (2004). [CrossRef]

18. R. B. Schulz, J. Peter, W. Semmler, C. D’Andrea, G. Valentini, and R. Cubeddu, “Quantifiability and image quality in noncontact fluorescence tomography,” Proc. SPIE 5859, 58590Z (2005). [CrossRef]

19. X. Chen, X. Gao, X. Qu, J. Liang, L. Wang, D. Yang, A. Garofalakis, J. Ripoll, and J. Tian, “A study of photon propagation in free-space based on hybrid radiosity-radiance theorem,” Opt. Express 17(18), 16266–16280 (2009), http://www.opticsinfobase.org/VJBO/viewmedia.cfm?uri=oe-17-18-16266&seq=0. [CrossRef] [PubMed]

20. R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Noncontact optical tomography of turbid media,” Opt. Lett. 28(18), 1701–1703 (2003). [CrossRef] [PubMed]

21. Q. Z. Zhang, L. Yin, Y. Y. Tan, Z. Yuan, and H. B. Jiang, “Quantitative bioluminescence tomography guided by diffuse optical tomography,” Opt. Express 16(3), 1481–1486 (2008), http://www.opticsinfobase.org/josab/viewmedia.cfm?id=149907&seq=0. [CrossRef] [PubMed]

22. H. Meyer, A. Garofalakis, G. Zacharakis, S. Psycharakis, C. Mamalaki, D. Kioussis, E. N. Economou, V. Ntziachristos, and J. Ripoll, “Noncontact optical imaging in mice with full angular coverage and automatic surface extraction,” Appl. Opt. 46(17), 3617–3627 (2007). [CrossRef] [PubMed]

23. M. Firbank, S. R. Arridge, M. Schweiger, and D. T. Delpy, “An investigation of light transport through scattering bodies with non-scattering regions,” Phys. Med. Biol. 41(4), 767–783 (1996). [CrossRef] [PubMed]

24. A. Ishimaru, Wave propagation and scattering in random media (Academic, New York, 1978).

25. S. Chandrasekhar, Radiative Transfer (Dover, New York, 1960).

26. L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A 1(6), 612–619 (1984). [CrossRef]

27. M. Aggarwal and N. Ahuja, “A pupil-centric model of image formation,” Int. J. Comput. Vis. 48(3), 195–214 (2002). [CrossRef]

28. G. Yan, J. Tian, S. Zhu, Y. Dai, and C. Qin, “Fast cone-beam CT image reconstruction using GPU hardware,” J. XRay Sci. Technol. 16, 225–234 (2008).

29. D. Qin, H. Zhao, Y. Tanikawa, and F. Gao, “Experimental determination of optical properties in turbid medium by TCSPC technique,” Proc. SPIE 6434, 64342E (2007). [CrossRef]

30. H. Li, J. Tian, F. Zhu, W. Cong, L. Wang, E. Hoffman, and G. Wang, “A mouse optical simulation environment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo method,” Acad. Radiol. 11(9), 1029–1038 (2004). [CrossRef] [PubMed]

31. R. Han, J. Liang, X. Qu, Y. Hou, N. Ren, J. Mao, and J. Tian, “A source reconstruction algorithm based on adaptive hp-FEM for bioluminescence tomography,” Opt. Express 17(17), 14481–14494 (2009), http://www.opticsinfobase.org/VJBO/viewmedia.cfm?uri=oe-17-17-14481&seq=0. [CrossRef] [PubMed]

32. B. J. Beattie, A. D. Klose, C. H. Le, V. A. Longo, K. Dobrenkov, J. Vider, J. A. Koutcher, and R. G. Blasberg, “Registration of planar bioluminescence to magnetic resonance and x-ray computed tomography images as a platform for the development of bioluminescence tomography reconstruction algorithms,” J. Biomed. Opt. 14(2), 024045 (2009). [CrossRef] [PubMed]

33. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225–4241 (2005). [CrossRef] [PubMed]

Compared objects	measurements	MOSE
Mean Error	0.0179	0.0214
Correlation Factor	0.9974	0.9957

	One view	Two view	Three view	Four view
Mean Error	0.0249	0.0206	0.0160	0.0156
Correlation Factor	0.9804	0.9858	0.9882	0.9881

	adipose	heart	lung	liver	kidney
$μ_{a} [m m^{- 1}]$	0.0050	0.0786	0.2630	0.4708	0.0881
$μ_{s}^{'} [m m^{- 1}]$	1.2273	1.0066	2.2091	0.7000	2.3585

Compared objects	measurements	MOSE
Mean Error	0.0179	0.0214
Correlation Factor	0.9974	0.9957

	One view	Two view	Three view	Four view
Mean Error	0.0249	0.0206	0.0160	0.0156
Correlation Factor	0.9804	0.9858	0.9882	0.9881

3D reconstruction of light flux distribution on arbitrary surfaces from 2D multi-photographic images

Abstract

1. Introduction

2. Theory and method

2.1 Lambertian source theory

2.2 Surface flux reconstruction algorithm

3. Experiments and results

3.1 Experiment verifications based on phantoms

3.1.1 Regular surface flux reconstruction

3.1.1.1 Cubic phantom experiment

3.1.1.2 Cylindrical phantom experiment

3.1.2 Irregular surface flux reconstruction

3.2 In vivo imaging experiment

4. Discussion and conclusion

Acknowledgments

References and links

Cited By

Figures (15)

Tables (3)

Equations (9)

Optics Express