Image Restoration for Fluorescence Planar Imaging with Diffusion Model

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.


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][5][6][7][8].
Up to the present, a number of fluorescence imaging techniques have been developed [9][10][11][12][13][14][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][17][18][19][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 wellknown 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, 2 BioMed Research International 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 0 = + [24][25][26], where denotes noise, is the system matrix, and 0 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][25][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][27][28][29][30][31][32][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][35][36][37], in which the diffusion model [38][39][40] is used to describe the photon propagation in tissues, the Born approximation [41][42][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.

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][36][37][38][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 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][42][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. em ( , ) denotes the photon density for emission at the position of detector when a point source is located at position . Φ ex ( , ) denotes the photon density for excitation at position when the source is located at . If the excitation source is a point source, Φ ex ( , ) is also a Green's function solution to the diffusion equation; otherwise, Φ ex ( , ) 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. 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 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.
G ？Ｇ (r dm , x i , y 1 , z 1 ) Φ ？Ｒ (x i , y 1 , z 1 , r s ) 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 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 . 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.
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 thecoordinate, 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: (1) , ( , , ) (1) . (9) In (9), the pixel value of the expected image is fixed as (1) 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 (1) 63 ; finally, the elements 3,63 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.
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.

Phantom Experiments.
Phantom experiments were carried out to validate the proposed image restoration method. The phantom was a 3 × 3 × 3.5 cm 3 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 alongaxis 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 711 ± 25 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.

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.

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 −1 and 10 cm −1 , respectively. The regularization parameter was empirically chosen as 10 −4 and 200 iterations were executed. Figure 7 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 (ℎ) = √ 2 − ℎ 2 , 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.
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, = 0.02 cm −1 , = 10 cm −1 , voxel size was 0.05 cm, and = 10 −4 ). 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 −1 and the reduced scattering coefficients are set as 5, 10, 15, and 20 cm −1 . Table 3 lists the FWHMs of the profiles of restored images with different optical coefficients. Figure 10 shows the  Figure 11 shows the restored images and profiles with regularization parameters of 10 −2 , 10 −4 , 10 −6 , and 10 −8 . 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.

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 1.8 × 2.9 × 2 cm 3 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 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.   Figures 12(a3), (b2), (c2), (d2), (e2), and (f2) are shown in Figure 12(g). All the images are normalized with the maximums. 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 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)   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 (1) corresponds to the first pixel of the blurry image and the weight of each voxel ( ) , = em ( , , , )Φ ex ( , , , ) 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, Φ ex and em 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 onplane and -plane separately. V and V denote the locations of the true target and the virtual target, respectively. Because the restored image (1) corresponds to the first pixel of the blurry image, we only consider the situation = 1 . In order to achieve the same pixel value 1 , the weights of the true target and the virtual target should be the same; that is, the distance V should equal V and 1 V should equal 1 V . Therefore, the virtual target should be located at the intersection of the focal plane and the circles with and 1 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 V 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 (‖ − true ‖ 2 2 ) 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. 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 −1 and 10 cm −1 when the wavelength is between 632.8 nm and 751 nm [49]. Therefore, we consider 0.02 cm −1 and 10 cm −1 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.

Discussion
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 3 × 3 × 3 cm cube and the voxel size is 0.05 cm, the image size will be 61 × 61 and the number of voxels will be 61 × 61 × 61. 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 em and Φ ex . For the above case, it is required to calculate em and Φ ex 61 × 61 × 61 + 61 × 61 ×(61×61−1) 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 961 × 961 and 961 × 31, respectively, and em and Φ ex 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    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 epiillumination 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 Φ ex 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.

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.