BioMed Research International
Volume 2017 (2017), Article ID 2010512, 18 pages
https://doi.org/10.1155/2017/2010512

Research Article

Image Restoration for Fluorescence Planar Imaging with Diffusion Model

Xuanxuan Zhang, Yuzhu Gong, Yang Li, Xu Cao, and Shouping Zhu

Engineering Research Center of Molecular and Neuro Imaging of the Ministry of Education & School of Life Science and Technology, Xidian University, Xi’an, Shaanxi 710071, China

Correspondence should be addressed to Shouping Zhu; spzhu@xidian.edu.cn

Received 6 June 2017; Accepted 5 November 2017; Published 27 November 2017

Academic Editor: Toshiyuki Sawaguchi

Copyright © 2017 Xuanxuan Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Fluorescence planar imaging (FPI) is failure to capture high resolution images of deep fluorochromes due to photon diffusion. This paper presents an image restoration method to deal with this kind of blurring. The scheme of this method is conceived based on a reconstruction method in fluorescence molecular tomography (FMT) with diffusion model. A new unknown parameter is defined through introducing the first mean value theorem for definite integrals. System matrix converting this unknown parameter to the blurry image is constructed with the elements of depth conversion matrices related to a chosen plane named focal plane. Results of phantom and mouse experiments show that the proposed method is capable of reducing the blurring of FPI image caused by photon diffusion when the depth of focal plane is chosen within a proper interval around the true depth of fluorochrome. This method will be helpful to the estimation of the size of deep fluorochrome.

1. Introduction

Fluorescence imaging techniques have become indispensable tools for numerous biomedical applications attributing to the everlasting development of fluorescent probes [1]. With the help of various fluorescent probes and fluorescence reporter techniques [2, 3], fluorescence imaging techniques are capable of tracing biomedical processes at cellular and subcellular levels in vivo and noninvasively in wide applications such as gene expression, protein function, and cell therapy [4–8].

Up to the present, a number of fluorescence imaging techniques have been developed [9–15]. Microscopic fluorescence imaging techniques provide high spatial resolutions but suffer from small fields of vision. On the contrary, macroscopic fluorescence imaging techniques can capture whole-body images for small animals but with a limited spatial resolution. Fluorescence planar imaging (FPI) [16–20] is the most widely used macroscopic fluorescence imaging technique, which directly detects the fluorescence photons on the surface of an imaged small animal using camera. According to the locations of excitation light source and camera, FPI can be formed in two different modes [1, 21]: epi-illumination mode and transillumination mode. Epi-illumination mode places excitation source and camera at the same side of the imaged small animal, which collects fluorescent photons in the same direction of the reflected excitation lights; thus it is also called fluorescence reflectance imaging (FRI). The defect of this mode is the difficulty of the excitation of deep fluorochromes. As an alternative, transillumination mode places the imaged small animal between excitation light source and camera. This mode can easily excite the fluorochromes far away from camera but the images are more heavily contaminated by excitation lights than epi-illumination mode although the excitation lights are attenuated by filters.

Whichever mode is applied, FPI is incapable of imaging deep fluorochromes with high spatial resolution. It is well-known that the penetration depth of near-infrared light in tissues is several centimeters [22]. Nevertheless, due to the elastic scattering, near-infrared photons are diffused after several millimeters of propagation in tissues [23]. So the fluorescent images acquired with camera are blurred. The deeper the fluorochromes are, the more strongly the fluorescent photons are diffused and the more blurry the images are. This restricts the applications of FPI in many cases. For instance, when imaging deep tumors, it is difficult to estimate the sizes of tumors through FPI images because they are strongly blurred.

Image restoration techniques aim to eliminate or reduce the impact of image degradation such as blurring. The causations of blurring can be classified into three types [24]: medium-induced, optical, and mechanical. The blurring derived from photon diffusion belongs to the first and second types due to the elastic scattering in medium. Blurring can be described with linear or nonlinear models, which depends on the specific problem. The general linear model can be summarized as [24–26], where denotes noise, is the system matrix, and and are the blurry and expected images, respectively. The key of deblurring is the construction of , which is known as the point spread function (PSF) in many applications. Because the linear model is usually formed with a convolution like [24–26], deblurring is also called deconvolution. During the last two decades, the researches of deblurring in fluorescence imaging focused on microscopic fluorescence imaging techniques [26–33] which are known as techniques with almost no photon diffusion. In these investigations, researchers implemented deconvolution methods to deal with the blurring derived from imaging system, that is, the mechanical type of blurring through PSFs of imaging system.

In this paper, we aim to build a method to reduce the impact of the blurring derived from photon diffusion in FPI. This will be helpful to the estimation of the size of deep fluorochrome. The scheme of the proposed image restoration method is conceived based on a reconstruction scheme in fluorescence molecular tomography (FMT) [34–37], in which the diffusion model [38–40] is used to describe the photon propagation in tissues, the Born approximation [41–43] is applied to solve the diffusion equation, and the Kirchhoff approximation [44, 45] is implemented to obtain Green’s function. Different from the blurring in fluorescence microscopic imaging, the blurring in FPI is not caused by imaging system. Consequently, the construction method of system matrix in fluorescence microscopic imaging is not applicable to the deblurring in FPI. The primary contribution of this work is the construction of the system matrix for FPI. Through introducing the first mean value theorem for definite integrals, we define a new unknown parameter as the restoration target rather than the fluorescent yield. The new unknown parameter is a weighted average of the voxel values of fluorescent yield along detection direction. To construct the system matrix that converts this parameter to the blurry image, depth conversion matrix is defined, which consists of the weights of the voxels with different depths related to the same pixel of the expected image. Subsequently, the elements of depth conversion matrices related to a chosen plane named focal plane are selected to construct the system matrix according to a proportional relationship. Finally, the Levenberg-Marquardt method [46, 47] is applied to solve the system equation and acquire the restored image. Phantom and mouse experiments are carried out to validate the proposed method.

2. Methods

The generation of fluorescence consists of two processes: excitation and emission. In the excitation process, photons from excitation source propagate to fluorochromes. Subsequently, fluorescent photons emitted from fluorochromes propagate to detectors in the emission process. Each process can be modeled by the diffusion equation with Robin-type boundary condition as follows [35–39]:

where denotes the position, is the photon density, is the source term, and and are the absorption and diffusion coefficients, respectively. The diffusion coefficient is defined as , where is the reduced scattering coefficient. Ω is the domain of the object and ∂Ω is the corresponding boundary. denotes the outward normal vector and is a coefficient related to the reflective index mismatch at boundary [40]. For the excitation process, the source term is determined by the location of excitation source and commonly approximated as an isotropy point source located one scattering length below the surface when a collimated source is used [39, 40]. As for the emission process, the source term is determined by the distribution of the photon density for excitation as well as the fluorescent yield of fluorochrome.

In this paper, the Born approximation [41–43] is used to solve the diffusion equations as follows:

where is Green’s function solution to the diffusion equation and is the fluorescent yield of fluorochrome. denotes the photon density for emission at the position of detector when a point source is located at position . denotes the photon density for excitation at position when the source is located at . If the excitation source is a point source, is also a Green’s function solution to the diffusion equation; otherwise, is the convolution of Green’s function and the distribution function of the source. is the fluorescent photon density for a pair of source and detector. The analytic formula of Green’s function solution to the diffusion equation can be achieved only for infinite space, semi-infinite space, and several simple geometries. To obtain Green’s function in geometries with arbitrary boundaries, the Kirchhoff approximation is implemented [44, 45].

Let us consider an imaging situation with transillumination mode as shown in Figure 1(a). A collimated source is used to excite fluorochrome and a planar detector is applied to capture fluorescent images. The imaged object is assumed to be a cube with a spherical fluorescent target located at the center. An illustration of the image restoration problem is shown in Figure 1(b). The fluorescent image acquired with the detector should be a blurry image due to the photon diffusion. The image we expect to achieve through the image restoration (hereinafter abbreviated as expected image) should be a projection along the detection direction, that is, -axis.

(a)

(b)

Figure 1: Schematic of the transillumination FPI and illustration of the image restoration problem. (a) Schematic of the transillumination FPI with a collimated source and planar detector. A sphere embedded in a cube is assumed to be the imaged fluorescent target. (b) Illustration of the image restoration problem.

Equation (2) can be written in a three-dimensional Cartesian coordinate system as follows:

Then we introduce the first mean value theorem for definite integrals:

where is a continuous function on , is an integrable function that does not change sign on , and is a value in . Through (3) and (4), the following equation can be formed:

In (5), a new parameter independent of -axis displaces the fluorescent yield . Subsequently, we discretize (5) with step lengths of , , and and obtain the following equation:

where is the value of the th pixel of the blurry image, denotes the corresponding pixel location, is a component of corresponding to the th pixel of the expected image , is the volume of voxel, and is the number of pixels. , , and are the numbers of voxels along -axis, -axis, and -axis, respectively. is a weight that converts the expected image to the blurry image and is a component of after the discretization of -axis where is the index of -coordinate.

An illustration of the linear model in (6) is shown in Figure 2, in which the image size is 8 × 8 and the pixel index of the blurry image is assumed to be 15. Equation (6) describes the relationship between the expected image and the blurry image. The pixel values of the expected image are weighted averages of the voxel values of fluorescent yield along -axis, which can be described with

Figure 2: Illustration of the linear model described by (6) with an image. The voxels in yellow correspond to the yellow elements in the map of and the expected image. The map of and the expected image relate to a pixel of the blurry image assumed to be the 15th pixel and colored in green.

According to (6), the weight varies with the pixel location of the blurry image . Consequently, based on (7) the pixel value of the expected image is not the same for different pixels of the blurry image; that is, the expected image is not unique. Figure 3 illustrates the nonuniqueness of . Figures 3(b) and 3(c) show the maps of for two different pixels ( and ) of the blurry image and Figure 3(d) gives profiles of along -axis for the two pixels. These figures indicate the differences between the weights corresponding to different pixels of the blurry image. Figures 3(e) and 3(f) are the images of calculated with the weights assuming that the fluorescent yield is known, which show that the images of are distinctly different for different pixels of the blurry image. In addition, the computational error results in the minute differences between the center four pixels. As a result of the nonuniqueness of , it is infeasible to construct a system equation that converts the expected image to the blurry image through a combination of (6) for all the pixels of the blurry image.

Figure 3: Illustration of the nonuniqueness of the expected image . (a) An 8 × 8 blurry image. (b) and (c) Maps of weights related to two different pixels ( and ) of the blurry image. (d) Weights for and as a function of index of -coordinate. (e) and (f) Images of calculated with related to (b) and (c).

In order to construct the system equation, firstly, we express (7) for all the pixels of the blurry image with the following matrix equation:

Because is not unique, we use a superscript on to denote the differences caused by different pixels of the blurry image in (8). The matrix on the left of (8) consists of the weights corresponding to all the pixels of the blurry image and a certain pixel of the expected image. The elements of each row of this matrix are arranged according to the -coordinate, that is, the depth. Each element determines the contribution of the fluorochrome at a certain depth to a certain pixel of the blurry image. Thus, we name this matrix depth conversion matrix.

To construct the system equation, we must build a set of equations that describe the relationship between the pixels of the blurry image and a stationary expected image. From (6), we know that the pixel values of the blurry image are the summations of . If we can displace the in (8) with the same , the relationship between the components and the pixel value of a stationary expected image will be formed; then the system equation can be constructed. We achieve this purpose through a proportional relationship derived from (8) as follows:

In (9), the pixel value of the expected image is fixed as but the weight is changed. The fluorescent yield is not known in image restoration. To obtain the weights, we manually choose a depth to approximate them as follows:

where the subscript denotes the index of -coordinate according to the chosen depth. The fluorescent signals are conceived to come from the plane at the chosen depth which is named focal plane.

Based on (6) and (10), we can construct the system equation as follows:

where is the system matrix that converts the expected image to the blurry image . Figure 4 is an illustration of the construction of system matrix, which shows the calculation process of the element of system matrix in the 3rd row and 63rd column. Firstly, the depth conversion matrix of the 63rd pixel of the expected image is calculated through (6); then the elements at the 5th column are selected according to a chosen depth shown as the red dotted line in the top center subfigure and the elements of the 1st row are summed to calculate ; finally, the elements of the system matrix are composed through (12). Although the system matrix is a square matrix, the inversion of is an ill-posed problem in practical application. The Levenberg-Marquardt method [46, 47] is implemented to solve (11) as follows:

where denotes the vector of the expected image for the th iteration, is the vector of the blurry image, is the regularization parameter, is the trace of , and is the identity matrix.

Figure 4: Illustration of the construction of system matrix. The calculation of the element of system matrix in the 3rd row and 63rd column is illustrated. Squares in blue and green denote the corresponding locations of the pixels of the expected image and the blurry image, respectively. The focal plane is shown with a red dotted line.

In general, the image restoration consists of three steps. Firstly, the depth conversion matrices for all the pixels of the expected image are calculated. Subsequently, the system matrix is assembled with the depth conversion matrices and a chosen focal plane. Finally, the system equation is solved with the Levenberg-Marquardt method to achieve the restored image.

3. Experimental Setup

3.1. Phantom Experiments

Phantom experiments were carried out to validate the proposed image restoration method. The phantom was a  cm³ cuboid tank made of perspex as shown in Figure 5(a). The cuboid tank was filled with diluted Intralipid-20% with a volume concentration of 5% and the height of the Intralipid-20% in the tank was 3 cm. A transparent glass capillary tube with a diameter of 0.3 cm was immersed in the tank. Holes were drilled on the wall of the tank with a thickness of 5 mm for the fastening of the tube. The depth of the tube was 2 cm. The distance between the center of the tube and the boundary along -axis was 1 cm. The distance between the bottom of the tube and the boundary along -axis was 1.1 cm. 20 μL Cy5.5 dye with five different concentrations of 4, 6, 8, and 10 μmol/L was successively filled into the tube as the fluorescent target. The fluorescent images of the phantom were acquired with a transillumination FPI system as shown in Figure 5(b). The phantom was placed on a lift table and an electron multiplying charge-coupled device (EMCCD) camera (iXon3 888, Andor Technologies, UK) coupled with a 50 mm f/0.95 lens (DO-5095, Navitar, USA) was placed above the phantom to capture images. A  nm bandpass filter (BrightLine, Semrock, USA) in front of the camera was used to capture fluorescent image. A 671 nm laser (CrystaLaser LC, Reno, USA) below the lift table was used to excite the fluorescent target within the phantom through a pinhole at the center of the lift table. Two daylight lamps placed on both sides of the camera were used to obtain white light images. The whole imaging system was covered with a black box to block environmental lights.

(a)

(b)

Figure 5: Phantom and imaging system. (a) Geometry of the phantom. (b) Schematic of the transillumination FPI.

3.2. Mouse Experiment

For the further validation of the proposed image restoration method, mouse experiment was implemented, which was conducted under the protocol approved by the Fourth Military Medical University Animal Care and Use Committee. The imaged object was a nude mouse. Before the experiment, the mouse was anesthetized with 10% sodium pentobarbital through intraperitoneal injection. A transparent glass capillary tube with a diameter of 0.3 cm filled with 20 μL Cy5.5 dye with a concentration of 1 μmol/L was planted into the abdomen of the mouse. Subsequently, the mouse with a supine position was fastened on a cardboard with black tapes. A slot is located at the center of the cardboard to let excitation light get through. The used FPI system was the same with the phantom experiments. After the acquisition of fluorescent images, X-ray computed tomography (X-CT) scan of the mouse was carried out and the filtered back-projection (FBP) method [48] was applied to accomplish the reconstruction of X-CT. Then the surface of a section of the mouse torso, which is required in the image restoration, was extracted through the segmentations of the X-CT images. The X-CT reconstruction results of the mouse and the extracted surface are shown in Figure 6.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 6: X-CT reconstruction results of the imaged mouse. (a)–(c) Three-dimensional volume rendering of the reconstruction results. (d)–(f) Surface of a section of the mouse torso extracted from X-CT images fused with volume rendering. (g)–(i) Representative slices of the reconstruction results. The first to third columns are the results for coronal, transversal, and sagittal positions, respectively.

4. Results

4.1. Phantom Experiments

The results of phantom experiments are shown in Figure 7. For all the image restoration results, the depth of focal plane (DFP) was set as 2 cm, the voxel size was set as 0.05 cm, and the optical coefficients and were set as 0.02 cm⁻¹ and 10 cm⁻¹, respectively. The regularization parameter was empirically chosen as 10⁻⁴ and 200 iterations were executed. Figure 7(a) is an illustration of the location of target, focal plane, and detected plane in -plane. Figures 7(b)–7(e) are the original fluorescent images for concentrations of 4, 6, 8, and 10 μmol/L normalized with the maximum of (e) while Figures 7(f)–7(i) are the corresponding restored images normalized with the maximum of (i). The magenta square in Figures 7(f)–7(i) shows the true location of target. Figure 7(j) provides profiles along the white dotted line in (e)–(i) as well as the true profile. Because the tube is a cylinder, the true profile is a curve with a formula of , where is the radius of the tube and is the distance away from the center of the tube. Figure 7(k) shows a linear fitting of the maximum of (f)–(i). The full widths at half maximum (FWHMs) of the profiles of the original and the restored images in Figure 7 are shown in Table 1.

Table 1: Full widths at half maximum of the profiles of the original and restored images with different concentrations of Cy5.5.

Concentration (μmol/L)
4
6
8
10

FWHM (cm)

    Restored
0.21
0.21
0.22
0.22

    Original
1.86
1.85
1.86
1.87

Figure 7: Results of phantom experiments. (a) Illustration of the location of target, focal plane, and detected plane in -plane. (b)–(e) Original fluorescent images for concentrations of 4, 6, 8, and 10 μmol/L, respectively. (f)–(i) Restored images of (b)–(e) in which the magenta square denotes the true location of target. (j) Profiles along the white dotted line in (e)–(i). (k) Linear fitting of the maximum of (f)–(i). The parameters are given in the top yellow box.

The DFP, the optical coefficients, the voxel size, and the regularization parameter affect the results. We tested the effect of the four factors through restoring the fluorescent image for the concentration of 8 μmol/L. During the test of each factor, the tested factor was changed while the other three were set as above (i.e., DFP was 2 cm,  cm⁻¹,  cm⁻¹, voxel size was 0.05 cm, and ). The two optical coefficients were also tested separately through changing one and fixing another. The results are shown in Figures 8–11. Figure 8 provides the restored images and profiles when the DFPs are 3 cm, 2.5 cm, 2 cm, 1.5 cm, and 1 cm and Table 2 shows the corresponding deviations of the centers of restored targets from the true center. Figure 9 gives the restored images and profiles when the absorption coefficients are set as 0.01, 0.02, 0.03, and 0.04 cm⁻¹ and the reduced scattering coefficients are set as 5, 10, 15, and 20 cm⁻¹. Table 3 lists the FWHMs of the profiles of restored images with different optical coefficients. Figure 10 shows the restored images and profiles with voxel sizes of 0.05, 0.075, 0.1, and 0.15 cm and Table 4 provides the corresponding FWHMs as well as the computational time. To show the voxel size in the images, Figures 10(c)–10(f) are shown without interpolation. Computational time as a function of voxel size with a sampling interval of 0.01 cm is also provided as Figure 10(f). Figure 11 shows the restored images and profiles with regularization parameters of 10⁻², 10⁻⁴, 10⁻⁶, and 10⁻⁸. Table 5 lists the corresponding derivations and FWHMs. For Figure 11(f), the lower right bulk is considered as the restored target. All the images are normalized with the maximums.

Table 2: Deviations of the centers of restored targets from the true center with different DFPs.

DFP (cm)
3
2.5
2
1.5
1

Deviation (cm)
0.26
0.10
0.04
0.31
0.74

Table 3: Full widths at half maximum of the profiles of restored images with different optical coefficients.

(cm⁻¹)
0.01
0.02
0.03
0.04
0.02
0.02
0.2
0.02

(cm⁻¹)
10
10
10
10
5
10
15
20

FWHM (cm)
0.21
0.22
0.22
0.21
0.22
0.22
0.22
0.21

Table 4: Full widths at half maximum of the profiles of restored images and computational time with different voxel sizes.

Voxel size (cm)
0.05
0.075
0.1
0.15

FWHM (cm)
0.22
0.20
0.19
0.15

Computational time (s)
201.0
23.1
5.4
0.5

Table 5: Deviations of the centers of restored targets from the true center and FWHMs of the profiles of restored images with different regularization parameters.

λ
10⁻²
10⁻⁴
10⁻⁶
10⁻⁸

Deviation (cm)
0.13
0.04
0.30
0.51

FWHM (cm)
0.39
0.22
0.16
0.10

Figure 8: Effect of DFP. (a)–(e) Illustration of the location of target, focal plane, and detected plane in -plane. (f)–(j) Restored images corresponding to (a)–(e). (k) Original fluorescent image. (l) Profiles along the white dotted line in (f)–(j). The fixed parameters are given in the top yellow box.

Figure 9: Effect of optical coefficients. (a) Illustration of the location of target, focal plane, and detected plane in -plane. (b) Original fluorescent image. (c)–(f) Restored images when the absorption coefficients are 0.01, 0.02, 0.03, and 0.04 cm⁻¹, respectively. (g)–(j) Restored images when the reduced scattering coefficients are 5, 10, 15, and 20 cm⁻¹, respectively. (k) Profiles along the white dotted line in (c)–(f). (l) Profiles along the white dotted line in (g)–(j). The fixed parameters are given in the yellow box.

Figure 10: Effect of voxel size. (a) Original fluorescent image. (b)–(e) Restored images when the voxel sizes are 0.05, 0.075, 0.1, and 0.15 cm, respectively. (f) Computational time as a function of voxel size with a sampling interval of 0.01 cm. (g) Profiles along the white dotted line in (c)–(f). The fixed parameters are given in the top yellow box.

Figure 11: Effect of regularization parameter. (a) Original fluorescent image. (b) Illustration of the location of target, focal plane, and detected plane in -plane. (c)–(f) Restored images when the regularization parameter is 10⁻², 10⁻⁴, 10⁻⁶, and 10⁻⁸, respectively. (g) Profiles along the white dotted line in (c)–(f). The fixed parameters are given in the top yellow box.

4.2. Mouse Experiment

The results of mouse experiment are shown in Figure 12. A part of the torso of the mouse with a size of about  cm³ was used to model the light propagation. For fine images, a smaller voxel size, 0.03 cm, than the values in the phantom experiments was chosen. The optical coefficients and were set as 0.3 cm⁻¹ and 10 cm⁻¹, respectively. The regularization parameter was empirically chosen as 10⁻³ and the image restorations were terminated after 200 iterations. The white light image, original fluorescent image, and a fused image of them are given in Figures 12(a1)–12(a3). Five different DFPs (1.1 cm, 0.9 cm, 0.67 cm, 0.4 cm, and 0.2 cm) were tested and the corresponding results are shown as Figures 12(b1)–12(f3). The results are shown with the restored image fused with the white light image as well as a coronal X-CT projection image. Profiles along the white dotted line in Figures 12(a3), (b2), (c2), (d2), (e2), and (f2) are shown in Figure 12(g). All the images are normalized with the maximums.

Figure 12: Results of mouse experiment. (a1)–(a3) White light image, original fluorescent image, and a fused image of them, respectively. (b1), (c1), (d1), (e1), and (f1) Illustrations of the locations of focal planes on a part of a sagittal X-CT projection image. (b2), (c2), (d2), (e2), and (f2) Restored images fused with white light images when the DFPs are set as 1.1 cm, 0.9 cm, 0.67 cm, 0.4 cm, and 0.2 cm, respectively. (b3), (c3), (d3), (e3), and (f3) Corresponding restored images fused with a coronal X-CT projection image. (g) Profiles along the white dotted line in (a3), (b2), (c2), (d2), (e2), and (f2). All the scale bars denote 0.5 cm. The fixed parameters are shown in the yellow box.

5. Discussion

Figure 7 and Table 1 show that the proposed method is capable of restoring the blurry images caused by photon diffusion when the DFP equals the depth of target. The profiles in Figure 7(j) and Table 1 demonstrate that the size of the restored fluorescent target in terms of FWHM (~0.22 cm) is close to the size of the true target (0.3 cm) while the FWHMs of the original fluorescent images are around 1.86 cm. Figure 7(k) indicates the pixel values of the restored images are approximately proportional to the concentration of Cy5.5 dye. The quality of image restoration relies on the choice of DFP as shown in Figure 8 and Table 2. When the DFP is deeper than the depth of target, the target can also be well restored but the restored location deviates from the true location as shown in Figures 8(f) and 8(g) as well as the deviations listed in Table 2. The deeper the DFP is, the farther the restored location deviates from the true location. For example, the deviation for a DFP of 3 cm is 0.26 cm but 0.10 cm for a DFP of 2.5 cm. On the other hand, when the DFP is shallower than the depth of target, the restored target tends to be spread around the image as shown in Figures 8(i) and 8(j). As referred in the derivation of (10), the fluorescent signals are conceived to come from the focal plane. Consequently, if the focal plane is not located at the true depth of target, the image restoration will result in a virtual target in the focal plane and the location of this virtual target varies with DFP. From (10) and (6), we know that the restored image corresponds to the first pixel of the blurry image and the weight of each voxel is a function of the location of the pixel of the blurry image , the location of source , and the location of voxel . If we ignore the effect of boundary condition, and are functions of the distance between and and the distance between and [45]. We use Figure 13 to illustrate the effect of DFP, where the fluorescent target is assumed to be a point. In Figure 13, the problem is illustrated on -plane and -plane separately. and denote the locations of the true target and the virtual target, respectively. Because the restored image corresponds to the first pixel of the blurry image, we only consider the situation . In order to achieve the same pixel value , the weights of the true target and the virtual target should be the same; that is, the distance should equal and should equal . Therefore, the virtual target should be located at the intersection of the focal plane and the circles with and as their centers as shown in the left column in Figure 13. However, the intersection does not exist for most situations as shown in the right column in Figure 13. In these situations, the voxel within the circle intersecting the focal plane and with a weight closest to the weight at could be considered as the virtual target. The above analysis ignores the effect of boundary condition and the target is a point. The actual situation is much more complicated due to the boundary of the imaged object as well as the volume of the fluorescent target. We analyze the location of the virtual target through finding an equivalent of the true target within the focal plane. However, the virtual target may not be well restored because the weight of the virtual target may be quite different from the one of the true target. Figure 14 provides the error of system matrix () as a function of DFP for the phantom experiments, where the depth of target is 2 cm and the system matrix with a DFP of 2 cm is assumed as the true system matrix. It can be observed in Figure 14 that the error of system matrix with a shallower focal plane rises much faster than the one with a deeper focal plane. This explains why the target cannot be well restored when the DFP is shallower than the depth of target in Figures 8(i) and 8(j). The results of mouse experiment in Figure 12 and Table 6 are consistent with the results of phantom experiments in Figure 8. Because of the tiny size of the mouse, the locations of the restored target vary slightly when the DFP is deeper than the true depth of target in Figure 12. In general, Figures 7, 8, and 12 as well as Tables 1, 2, and 6 indicate the proposed method is capable of restoring the blurry images caused by photon diffusion when the DFP is set properly. The target can be well revealed when the DFP is within an interval around the true depth of target rather than an exact value.

Table 6: Full widths at half maximum of the profiles of restored images in mouse experiment with different DPFs.

DFP (cm)
0.2
0.4
0.67
0.9
1.1
Original

FWHM (cm)
0.40
0.35
0.22
0.19
0.18
0.88

Figure 13: Illustration of the effect of DFP with an 8 × 8 image and a point target. The left column shows the case when the intersection of the focal plane and the circles with and as their centers exists. The right column shows the case when the intersection does not exist. The middle column shows the blurry image and the expected image. The pixel in light blue denotes the true location of target while the pixel in dark blue denotes the location of virtual target in the focal plane.

Figure 14: Error of system matrix as a function of DFP for the phantom experiments. The depth of target is 2 cm.

Generally, the proposed method is more appropriate for the cases in which the depths of targets can be estimated because the location of the restored target varies with DFP. This restricts the applications of the proposed method. On the contrary, we might take use of the effect of DFP to estimate the depth of target. It is known that commonly the location of the maximum of the blurry image should correspond to the center of the fluorescent target. Therefore, we could restore the blurry image with a set of DFPs ranging from 0 to the thickness of the imaged object and compare the center of the restored images with the center of the blurry image to estimate the depth of the target.

The optical coefficients and are prerequisites for the image restoration. In the phantom experiments, we used the diluted Intralipid-20% with a volume concentration of 5% as the diffusion medium, the absorption coefficient and reduced scattering coefficient of which are, respectively, around 0.02 cm⁻¹ and 10 cm⁻¹ when the wavelength is between 632.8 nm and 751 nm [49]. Therefore, we consider 0.02 cm⁻¹ and 10 cm⁻¹ as the true optical coefficients of the diffusion medium and tested the effect of the optical coefficients through setting them as 50%, 100%, 150%, and 200% of the true values. The results in Figure 9 and Table 3 indicate that the effect of optical coefficients is slight (FWHMs are around 0.22). Based on this conclusion, we used homogeneous optical coefficients in the mouse experiment rather than a heterogeneous model.

The results of the image restoration greatly depend on the voxel size as shown in Figure 10 and Table 4. In general, a small voxel size would result in a fine image but consumes more time on computation and more memories on the storage of system matrix and depth conversion matrices. For example, when the imaged object is a  cm cube and the voxel size is 0.05 cm, the image size will be and the number of voxels will be . It results in the fact that the size of system matrix is 3,721 × 3,721 and the size of each depth conversion matrix is 3,721 × 61. When the data are saved with double precision, the system matrix and each depth conversion matrix will occupy about 105.6 and 1.7 megabytes of memories, respectively. For the storage of all the depth conversion matrices and the system matrix, about 6.4 gigabytes is required. To reduce the requirement of the memories, we can save only the elements of depth conversion matrices required during the calculation of system matrix. Moreover, the calculations of system matrix focus on the calculations of weights , that is, the calculations of and . For the above case, it is required to calculate and 61 × 61 × 61 + 61 × 61 × times, that is, 14,069,101 times at least for a single DFP. In parallel, if the voxel size is 0.1 cm, the sizes of system matrix and each depth conversion matrix will be and , respectively, and and will be calculated 952,351 times for a single DFP. Through Figure 10(f) and Table 4, it can be found that the computational time exponentially increases with the decrease of the voxel size. For example, the computational time is only 0.5 seconds for a voxel size of 0.15 cm but 201 seconds for a voxel size of 0.05 cm; that is to say, when the voxel size reduces 3-fold, the computation time increases about 400-fold. This limits the choice of voxel size. In general, the voxel size should be as small as possible when the computational time and memories are not limited.

Figure 11 and Table 5 show the effect of regularization parameter. The function of the regularization parameter is controlling the effect of regularization. The regularization aims to suppress the impact of ill-posedness and noise through smoothing the solution vector. A large regularization parameter will result in oversmoothness like Figure 11(c) with a FWHM of 0.39 cm; on the contrary, a small one may lead to image distortion because of noise like Figures 11(e) and 11(f) with deviations of 0.30 cm and 0.51 cm. The choice of the regularization parameter depends on the noise level of data. Due to the fact that noise level is unable to be achieved in some applications, the regularization parameter is usually determined empirically [50, 51].

A defect of the proposed method is the requirement of the geometry of the imaged object, which is not necessary in FPI. The geometry is a prerequisite for modeling the photon propagation in the imaged object. Thus, an additional imaging modality such as X-CT and magnetic resonance imaging (MRI) is indispensable. As an alternative, the imaged object can be immersed into a tank with regular geometry filled with matching fluids.

Although the proposed method in this paper is based on a transillumination FPI system, it may be feasible for epi-illumination mode in theory as long as the location of source is available. If the source is still a point source, the methods for the two modes have no difference except the locations of source. When a planar source is implemented, the distribution of the photon density for excitation is not a Green’s function solution to the diffusion equation anymore but a convolution of Green’s function and the distribution function of the source. In general, transillumination mode is more appropriate because the aim of the proposed method is the imaging of deep fluorescent targets while transillumination mode can easily excite deep fluorochromes.

The proposed method is conceived based on the reconstruction scheme in FMT. It is well known that FMT reconstructs a three-dimensional distribution of fluorescent yield with multiple projection images from different angles, while the proposed method can be considered as the reconstruction of a two-dimensional image from a single blurry image. It is infeasible for the reconstruction method in FMT to directly reconstruct the distribution of fluorescent yield from a single image because the mismatch between the quantities of data and unknown parameters will result in great degradation of reconstruction results [52]. Through the definition of a new unknown parameter in (6), the proposed method balances the quantities of data and unknown parameters to make the reconstruction possible. In addition, although FMT can provide three-dimensional distribution of fluorescent yield, it cannot displace FPI because of the stability, simplicity, convenience, and fast data acquisition of FPI. Therefore, there are still quantities of applications accomplished with FPI especially those requiring high data acquisition speed such as pharmacokinetics [16, 17].

The results of phantom and mouse experiments prove that the proposed image restoration method is capable of revealing fluorescent targets beneath diffusion medium. It is well known that the near-infrared light is highly scattered within the tissues of living bodies and it is difficult to capture the shapes of deep fluorescent targets beneath tissues with cameras. This limits the applications of FPI. For example, size of tumor is an important indicator in oncology. Fluorescence planar imaging is a usual imaging technique for the research of tumor on animal model through fluorescent molecular probes. However, the estimation of size of tumor through FPI images is feasible for only subcutaneous tumor while the size of in situ tumor is difficult to be estimated. The proposed method provides the ability of the estimation of size of deep fluorescent targets like in situ tumors.

6. Conclusion

In conclusion, an image restoration method is proposed to deal with the blurring caused by photon diffusion in FPI. The method is conceived based on the reconstruction scheme in FMT. The primary contribution of this work is the construction of system matrix, which is achieved through the definition of a new unknown parameter and depth conversion matrices with a chosen focal plane. The new unknown parameter is defined through the first mean value theorem for definite integrals and represents a weighted average of the fluorescent yields along the detection direction. Results of phantom and mouse experiments indicate that the proposed image restoration method is capable of reducing the blurring of fluorescent image caused by photon diffusion when the depth of focal plane is chosen within a proper range around the true depth of fluorochrome. The effect of optical coefficients is slight. The quality of the restored image greatly depends on the voxel size but it is limited by the computational time and memories. The regularization parameter also influences the results that a large regularization parameter would result in oversmoothness while a small one might lead to image distortion. This method will be helpful to the estimation of the size of deep fluorochrome.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the Program of National Natural Science Foundation of China under Grants nos. 81227901, 61405149, 81230033, and 61471279, the Program of the National Key Research and Development Program of China under Grant no. 2016YFC0103802, and the Fundamental Research Fund for the Central Universities under Grant no. XJS17049.

References

V. Ntziachristos, “Fluorescence molecular imaging,” Annual Review of Biomedical Engineering, vol. 8, pp. 1–33, 2006.

R. Weissleder and V. Ntziachristos, “Shedding light onto live molecular targets,” Nature Medicine, vol. 9, no. 1, pp. 123–128, 2003.

T. F. Massoud and S. S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes & Development, vol. 17, no. 5, pp. 545–580, 2003.

R. Y. Tsien, “The green fluorescent protein,” Annual Review of Biochemistry, vol. 67, pp. 509–544, 1998.

M. Funovics, R. Weissleder, and C.-H. Tung, “Protease sensors for bioimaging,” Analytical and Bioanalytical Chemistry, vol. 377, no. 6, pp. 956–963, 2003.

R. Y. Tsien, “Building and breeding molecules to spy on cells and tumors,” FEBS Letters, vol. 579, no. 4, pp. 927–932, 2005.

C. Chamberlain and K. M. Hahn, “Watching proteins in the wild: Fluorescence methods to study protein dynamics in living cells,” Traffic, vol. 1, no. 10, pp. 755–762, 2000.

Y. Li, W. Liu, and F. Liu, “Primed 3D injectable microniches enabling low-dosage cell therapy for critical limb ischemia,” Proceedings of the National Acadamy of Sciences of the United States of America, vol. 111, no. 37, pp. 13511–13516, 2014.

G. Alexandrakis, E. B. Brown, R. T. Tong et al., “Two-photon fluorescence correlation microscopy reveals the two-phase nature of transport in tumors,” Nature Medicine, vol. 10, no. 2, pp. 203–207, 2004.

X. Liu, X. Guo, F. Liu et al., “Imaging of indocyanine green perfusion in mouse liver with fluorescence diffuse optical tomography,” IEEE Transactions on Biomedical Engineering, vol. 58, no. 8, pp. 2139–2143, 2011.

A. Ruggiero, J. P. Holland, J. S. Lewis, and J. Grimm, “Cerenkov luminescence imaging of medical isotopes,” Journal of Nuclear Medicine, vol. 51, no. 7, pp. 1123–1130, 2010.

Z. Hu, Y. Qu, K. Wang et al., “In vivo nanoparticle-mediated radiopharmaceutial-excited fluorescence molecular imaging,” Nature Communications, vol. 6, p. 7560, 2015.

X. T. Ding, K. Wang, B. Jie, Y. L. Luo, Z. H. Hu, and J. Tian, “Probability method for Cerenkov luminescence tomography based on conformance error minimization,” Biomedical Optics Express, vol. 5, no. 7, pp. 2091–2112, 2014.

Y. Lv, J. Tian, W. Cong et al., “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Optics Express, vol. 14, no. 18, pp. 8211–8223, 2006.

A. Zaheer, R. E. Lenkinski, A. Mahmood, A. G. Jones, L. C. Cantley, and J. V. Frangioni, “In vivo near-infrared fluorescence imaging of osteoblastic activity,” Nature Biotechnology, vol. 19, no. 12, pp. 1148–1154, 2001.

W. Cai, H. Guang, C. Cai, and J. Luo, “Effects of temperature on multiparametric evaluation of hindlimb ischemia with dynamic fluorescence imaging,” Journal of Biophotonics, vol. 10, no. 6-7, pp. 811–820, 2017.

H. Guang, C. Cai, S. Zuo, W. Cai, J. Zhang, and J. Luo, “Multiparametric evaluation of hindlimb ischemia using time-series indocyanine green fluorescence imaging,” Journal of Biophotonics, vol. 10, no. 3, pp. 456–464, 2017.

M. Choi, K. Choi, S.-W. Ryu, J. Lee, and C. Choi, “Dynamic fluorescence imaging for multiparametric measurement of tumor vasculature,” Journal of Biomedical Optics, vol. 16, no. 4, Article ID 046008, 2011.

V. Ntziachristos, G. Turner, J. Dunham et al., “Planar fluorescence imaging using normalized data,” Journal of Biomedical Optics, vol. 10, no. 6, Article ID 064007, 2005.

W. Buchalla, Á. M. Lennon, M. H. Van Der Veen, and G. K. Stookey, “Optimal camera and illumination angulations for detection of interproximal caries using quantitative light-induced fluorescence,” Caries Research, vol. 36, no. 5, pp. 320–326, 2002.

F. Leblond, S. C. Davis, P. A. Valdés, and B. W. Pogue, “Pre-clinical whole-body fluorescence imaging: review of instruments, methods and applications,” Journal of Photochemistry and Photobiology B: Biology, vol. 98, no. 1, pp. 77–94, 2010.

F. F. Jobsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science, vol. 198, no. 4323, pp. 1264–1267, 1977.

V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nature Methods, vol. 7, no. 8, pp. 603–614, 2010.

T. F. Chan and J. Shen, Theory and Computation of Variational Image Deblurring, World Scientific, Singapore, 2006.

M. Rostami, O. Michailovich, and Z. Wang, “Image deblurring using derivative compressed sensing for optical imaging application,” IEEE Transactions on Image Processing, vol. 21, no. 7, pp. 3139–3149, 2012.

J. Markham and J.-A. Conchello, “Fast maximum-likelihood image-restoration algorithms for three-dimensional fluorescence microscopy,” Journal of the Optical Society of America A: Optics and Image Science, and Vision, vol. 18, no. 5, pp. 1062–1071, 2001.

P. J. Verveer, M. J. Gemkow, and T. M. Jovin, “A comparison of image restoration approaches applied to three- dimensional confocal and wide-field fluorescence microscopy,” Journal of Microscopy, vol. 193, no. 1, pp. 50–61, 1999.

H. Yoo, I. Song, and D.-G. Gweon, “Measurement and restoration of the point spread function of fluorescence confocal microscopy,” Journal of Microscopy, vol. 221, no. 3, pp. 172–176, 2006.

A. Squire and P. I. H. Bastiaens, “Three dimensional image restoration in fluorescence lifetime imaging microscopy,” Journal of Microscopy, vol. 193, no. 1, pp. 36–49, 1999.

D. Sud and M.-A. Mycek, “Image restoration for fluorescence lifetime imaging microscopy (FLIM),” Optics Express, vol. 16, no. 23, pp. 19192–19200, 2008.

J. Kim, S. An, S. Ahn, and B. Kim, “Depth-variant deconvolution of 3D widefield fluorescence microscopy using the penalized maximum likelihood estimation method,” Optics Express, vol. 21, no. 23, pp. 27668–27681, 2013.

D. U. Campos-Delgado, O. G. Navarro, E. R. Arce-Santana, A. J. Walsh, M. C. Skala, and J. A. Jo, “Deconvolution of fluorescence lifetime imaging microscopy by a library of exponentials,” Optics Express, vol. 23, no. 18, pp. 23748–23767, 2015.

N. Patwary and C. Preza, “Image restoration for three-dimensional fluorescence microscopy using an orthonormal basis for efficient representation of depthvariant point-spread functions,” Biomedical Optics Express, vol. 6, no. 10, article no. A020, pp. 3826–3841, 2015.

V. Ntziachristos, C.-H. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nature Medicine, vol. 8, no. 7, pp. 757–760, 2002.

H. Yi, X. Zhang, J. Peng et al., “Reconstruction for Limited-Projection Fluorescence Molecular Tomography Based on a Double-Mesh Strategy,” BioMed Research International, vol. 2016, Article ID 5682851, pp. 1–11, 2016.

H. Zhang, G. Geng, X. Wang, X. Qu, Y. Hou, and X. He, “Fast and Robust Reconstruction for Fluorescence Molecular Tomography via L1-2 Regularization,” BioMed Research International, vol. 2016, Article ID 5065217, pp. 1–9, 2016.

X. Zhang, F. Liu, S. Zuo et al., “Reconstruction of fluorophore concentration variation in dynamic fluorescence molecular tomography,” IEEE Transactions on Biomedical Engineering, vol. 62, no. 1, pp. 138–144, 2015.

S. R. Arridge and J. C. Schotland, “Optical tomography: Forward and inverse problems,” Inverse Problems, vol. 25, no. 12, Article ID 123010, 2009.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems, vol. 15, no. 2, pp. R41–R49, 1999.

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Medical Physics, vol. 22, no. 11, pp. 1779–1792, 1995.

F. Liu, X. Liu, D. Wang, B. Zhang, and J. Bai, “A Parallel excitation based fluorescence molecular tomography system for whole-body simultaneous imaging of small animals,” Annals of Biomedical Engineering, vol. 38, no. 11, pp. 3440–3448, 2010.

M. A. O'Leary, D. A. Boas, X. D. Li, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Optics Letters, vol. 21, no. 2, pp. 158–160, 1996.

V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Optics Letters, vol. 26, no. 12, pp. 893–895, 2001.

J. Ripoll, M. Nieto-Vesperinas, R. Weissleder, and V. Ntziachristos, “Fast analytical approximation for arbitrary geometries in diffuse optical tomography,” Optics Letters, vol. 27, no. 7, pp. 527–529, 2002.

J. Ripoll, V. Ntziachristos, R. Carminati, and M. Nieto-Vesperinas, “Kirchhoff approximation for diffusive waves,” Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 64, no. 5, Article ID 051917, 2001.

K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Quarterly of Applied Mathematics, vol. 2, pp. 164–168, 1944.

D. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM Journal on Applied Mathematics, vol. 11, no. 2, pp. 431–441, 1963.

I. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” Journal of the Optical Society of America A: Optics, Image Science & Vision, vol. 1, no. 6, pp. 612–619, 1984.

P. D. Ninni, F. Martelli, and G. Zaccanti, “Intralipid: towards a diffusive reference standard for optical tissue phantoms,” Physics in Medicine and Biology, vol. 56, no. 2, pp. N21–N28, 2011.

B. Brooksby, S. Jiang, H. Dehghani et al., “Combining near-infrared tomography and magnetic resonance imaging to study in vivo breast tissue: Implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” Journal of Biomedical Optics, vol. 10, no. 5, Article ID 051504, 2005.

G. Zhang, F. Liu, H. Pu, W. He, J. Luo, and J. Bai, “A direct method with structural priors for imaging pharmacokinetic parameters in dynamic fluorescence molecular tomography,” IEEE Transactions on Biomedical Engineering, vol. 61, no. 3, pp. 986–990, 2014.

X. Cao, B. Zhang, F. Liu, X. Wang, and J. Bai, “Reconstruction for limited-projection fluorescence molecular tomography based on projected restarted conjugate gradient normal residual,” Optics Letters, vol. 36, no. 23, pp. 4515–4517, 2011.