1. Introduction

Over the last few years, several optical tomography applications that use diffuse light to examine tissue in vivo have been developed. These techniques are able to retrieve the 3D distribution of a wide set of intrinsic and extrinsic optical contrasts, and have many preclinical and clinical applications, such as breast imaging [1–4] and 3D or quantitative imaging of fluorophore concentration in live animal samples, as in fluorescence diffuse optical tomography (fDOT). fDOT is also named fluorescence molecular tomography (FMT) when there are molecular imaging probes involved [5–7].

The fluorescence diffuse optical tomography (fDOT) forward problem can be expressed as a linear system of equations

Solving [Eq. (1)] is equivalent to finding the predicted data that are as close as possible to real data, in other words, solving the least square problem

One of the main issues of the fDOT inverse problem is that it is ill-posed in the Hadamard sense [8,9], yielding multiple non-unique and unstable solutions to the reconstruction problem. Therefore, small perturbations in acquired data may lead to arbitrarily large changes in the reconstructed images.

In order to obtain a meaningful solution the problem must be regularized. The most common method is the Tikhonov regularization, which consists in replacing the least squares problem [Eq. (2)] by the following functional minimization problem:

In many implementations, the regularization parameter, α, is manually selected, using a sequence of regularization parameters and selecting the parameter value that leads to the best results, as judged by the operator . However, this procedure is subjective and slow. To avoid this, several automatic methods for selecting regularization parameters have been suggested over the years, for instance: the Unbiased predictive risk estimator method (UPRE), the Discrepancy principle (DP) and strategies based on the calculation of the variance of the solution that require prior knowledge of the noise, or others such as Generalized cross-validation (GCV) and L-curve (LC) that do not need a-priori information.

More details of these methods can be found in a survey paper that describes some of these methods [10], and in the chapter 7 of the book [11]. As Vogel emphasized in his book, depending on the particular inverse problem some methods may perform better than others, and some methods may even fail.

Regarding fDOT, time-resolved (fluorescence) diffuse optical tomography, diffuse optical tomography and optical fluorescence tomography, the reported strategies for the selection of the regularization parameter have been the manual selection [12], selection by using some acceptable reconstruction variability (non-biased estimators: mean and standard deviation of the reconstruction) [13], choosing from the resolution-variance plot to correspond to equal noise variance [14], the L-curve [1,15,16] and a variant of the L-curve method [17].

L-curve method (LC) has been extensively analyzed [18,19] and applied in different areas. Recently, it has been found that L-curve returns good regularization parameter values for electrical impedance tomography of the brain [20] and for diffuse optical topography [21]. In both studies several methods were compared arriving to the conclusion that in practice, selection methods without a-priori information, such as GCV and L-curve, were more robust and there was not significant difference between GCV and L-curve in accuracy. The only relevant difference is that GCV is more computationally expensive for large systems than L-curve [22].

In terms of diffuse optical tomography, the authors of [15] emphasized that the L-curve analysis yielded overly-smooth solution in same cases.

Lately, the U-curve method has been proposed by [23,24] for the selection of the regularization parameter in inverse problem. The U-curve was tested on some numerical examples of the Fredholm integral equation of first kind [23,24] and super-resolution [25]. These works concluded that

We validated this method by confirming that Picard’s condition [26] is fulfilled and inspecting the noise level of the reconstructed images to ensure that the U-curve method yields a satisfactory regularized solution. Phantom and ex-vivo mice data were acquired with a custom-made fDOT system [27,28].

This paper is organized as follows: Section 2.1 describes our in-house parallel-plate non-contact fDOT setup. The forward problem is detailed in section 2.2. Section 2.3 briefly mentions the singular value decomposition (SVD), which provides a way of analyzing the ill-posedness of a problem. Section 2.4 introduces the L-curve method. Section 2.5 presents the selection of the regularization parameter by the U-curve method. The experiments are explained in section 2.6. Section 2.7 describes how we validate the parameter selection. In section 3, we summarize the main results obtained. We finalize by discussing the issues raised and draw conclusions in section 4.

2. Materials and methods

2.1 fDOT experimental set up

In the non-contact parallel-plate configuration [29], the study sample is gently compressed between two parallel anti-reflective plates.

The excitation laser beam enters the sample perpendicularly through the first plate, and the transmitted light emerging through the opposite plate is recorded with a CCD camera.

Figure 1 shows our parallel-plate fDOT experimental setup.

 

Fig. 1 Experimental setup.

Download Full Size | PDF

In our prototype, the light emerging from a laser diode is focused onto the sample at the desired points (source locations) using two mirrors moved by galvanometers, thus making it possible to choose the number and spatial distribution of sources. The laser emitter wavelength is set at 675±5 nm, and the power delivered to the sample is controlled by means of a TTL signal that modulates the laser duty cycle. Typical power values are in the 1-mW range. All the components of the set-up are placed inside a light-shielded box. The acquisition process is controlled by in-house software hosted on a PC workstation.

All fluorescence images are taken in transmission and are recorded by placing a 10-nm bandwidth filter centered at 720 nm in front of the camera lens, while for the transmitted excitation images a 10-nm bandwidth filter centered at 675 nm is used. The filters are placed in front of the camera using a motorized filter wheel (see Fig. 1). The acquired images are divided by their respective laser power.

For each source, a variable number and distribution of detectors can be defined over the CCD sensor field of view (FOV), thus making it possible to retrieve the fluorescent and excitation photon density at the desired points on the sample surface.

2.2 Weight matrix formulation. Forward problem

Calculating the weight matrix requires a theoretical model (forward problem) that predicts photon propagation through the diffusive medium. The propagation of photons through biological media can be described by the radiative transfer equation [30]. Under highly scattering conditions such as those presented in this work, the diffuse approximation is derived from the radiative transfer equation [30], which provides simple analytical and numerical expressions for light transport and can be applied to a great variety of scattering systems, thus enabling their optical characterization [31]. Photon propagation was modeled as described in [32,33], taking into account the effects of the planar boundaries of the slab in the light transport model.

In order to accurately obtain the 3D distribution of fluorophore concentration, we used the normalized Born approximation previously introduced in [6,34]. This involves normalizing the measured fluorescence intensity to the corresponding measured intensity at the excitation wavelength for each source rs and detector rd, and assuming the contribution of each voxel in the mesh to be correctly described by the Born approximation as described in[35]

We can write Eq. (5) as Eq. (1) to simplify.

2.3 Tikhonov solution based on SVD

A common solution to the Tikhonov minimization problem Eq. (3) is given in terms of the singular value decomposition (SVD) of the weight matrix [36]. The SVD of a matrix generates a triplet of matrices [36]

The solution of the inverse problem of Eq. (1), in terms of SVD, can be written analytically as [37]

Since this problem is ill-posed, the singular values gradually decay to zero. The left and right singular vectors associated with the smaller singular values tend to have more sign changes, thus increasing the number of unstable solutions. The purpose of the regularization is to filter out the contribution of small singular values to the solution. In particular, the Tikhonov regularization [Eq. (3)] can be written in terms of the solution [Eq. (7)] as [37]

2.4 L-curve method

The L-curve is a log-log plot, for the different $\alpha >0$, of the image norm $\left(\Vert f_{\alpha }\Vert \right)$ against the corresponding residual norm $\left(\Vert W{f}_{\alpha }-d\Vert \right)$ [18,19].The L-curve is L-shaped. Applying the L-curve routine available in the Matlab Regularization Toolbox [38], we take the value of the regularization parameter, ${\alpha }_l$, which corresponds to the point of maximum curvature on the graph.

2.5 Selection of the noise threshold by the U-curve method (regularization parameter)

The U-curve is a plot of the sum of the inverse of the regularized solution norm $\left(\Vert f_{\alpha }\Vert \right)$ and the corresponding residual norm $\left(\Vert W{f}_{\alpha }-d\Vert \right)$, for $\alpha >0$ on a log-log scale

The regularized solution norm $\left(\Vert f_{\alpha }\Vert \right)$ and the corresponding residual norm $\left(\Vert W{f}_{\alpha }-d\Vert \right)$are calculated numerically.

As previously mentioned, the weight matrices that describe our system were decomposed by SVD following [Eq. (6)]. Tikhonov regularization was used to reconstruct the images, and the regularization parameter was selected using the U-curve method explained above.

In the experiments presented here, α took 200 values equally spaced from $\alpha ={\sigma }_r^{2/3}$ to $\alpha ={\sigma }_0^{2/3}$.

2.6 Experiments

We present two fDOT reconstructions corresponding to two experiments: one consisted of a slab phantom in which the tip of a capillary containing a known concentration of fluorophore was inserted into the center. For the second experiment we used a euthanized mouse in which the tip of a capillary containing a known concentration of fluorophore was inserted into the esophagus.

2.6.1 Phantom experiments

We made a slab-shaped agar-based phantom (8 x 6 x 1.5 cm) using intralipid and India ink to obtain an absorption coefficient of approximately μ_a =0.3 cm⁻¹ and a reduced scattering coefficient of μ_s=10 cm⁻¹ [39]. A capillary with its tip filled with 6 µl at 30µM of Alexa Fluor 700 (Invitrogen, Carlsbad, California, USA) was inserted into the phantom, with the tip positioned at the center of the slab.

We built a weight matrix of fDOT data collected with 20 x 20 x 10 mesh points for a 1.5 x 1.5 x 1.5-cm FOV. Equally spaced 10 x 10 sources and 12 x 12 detectors, respectively, were selected. The center of the mesh FOV was aligned with the center of the slab.

Figure 2 shows that the mesh FOV for the reconstruction is three-dimensional and the sources and detectors are equally spaced in an area of [1.5 x 1.5] on the front and back plates (Fig. 2(c)).

 

Fig. 2 (a) Geometrical configuration of the slab (black) and the mesh FOV (red). The capillary tip is represented by the black sphere. (b) Detail of the mesh FOV. The grey dots represent the positions of the voxels; the red dot represents the capillary tip. (c) Geometrical configuration of the sources and detectors for the slice corresponding to y=0.75; the rectangles represent the reconstructed image voxels.

Download Full Size | PDF

2.6.2 Ex-vivo mouse data

A euthanized mouse was imaged with a capillary inserted into the esophagus. The tip of the capillary (<1.5 mm thick) was filled with 6 µl at 30 µM of Alexa Fluor 700 (Invitrogen, Carlsbad, California, USA).

We constructed a weight matrix of fDOT data assembled with 20 x 20 x 10 mesh points for a 1.4 x 1.4 x 1.5-cm FOV centered on the chest of the mouse. Equally spaced 6 x 6 sources and 10 x 10 detectors were selected. The mouse was gently compressed between two anti-reflective plates (until 1.5 cm approximately) in order to conform its geometry as much as possible to that of a slab. We can assume the axis x of the FOV along the width of the mouse, the axis y along the length of the mouse and the axis z along the height of the mouse.

2.7 Validation of the regularization parameter obtained by the U-curve method

Once the weight matrix was decomposed by SVD, the images were reconstructed using Tikhonov regularization with α parameters in the 10⁻¹ to 10⁻⁶ range, which included the U-curve-based regularization parameter, ${\alpha }_u$.

The noise content and resolution of the images for each α value were then evaluated. To assess image resolution, we followed the procedure described in [15], assuming that the FWHM of the point spread function (PSF) of the capillary tip (that can be considered as a single isolation region) is directly related to the resolution performance of the system. We measured the noise present in the images as the standard deviation in a region of the image with no signal. The resolution versus noise graph is plotted in Figs. 6 and 11

 

Fig. 6 (a) Profiles taken in the x direction, corresponding to the line z=y=7.5 mm for each regularization parameter. FWHM improves when the regularization parameter decreases. b) Resolution vs. noise plot.

Download Full Size | PDF

 

Fig. 11 (a) Profiles taken in the x direction, corresponding to the line z=8 mm, y=0.7 mm for each regularization parameter. While the regularization parameter decreases, the FWHM improves. (b) Resolution vs. noise plot.

Download Full Size | PDF

Generally, the use of low regularization parameters means that the resolution of the reconstructions improves while the image noise increases [15]. It is not possible to define an optimal parameter common to all imaging applications, since for each case the user may ask for different noise and resolution tolerances.

An indication of the ill-posedness of the problem is the rate of decrease of the singular values, σ_i. The Discrete Picard’s condition (DPC) [26] provides us with an objective assessment of this fact.

The data vector d is said to satisfy the DPC if the data space coefficients $\left|u_i^{T}d\right|$, on average, decay to zero faster than the singular values σ_i [26].

To compute a satisfactory solution by means of the Tikhonov regularization, Picard’s condition has to be fulfilled [26], since it determines how well the regularized solution approximates the unknown, exact solution. In ill-posed problems, there may be a point, where the data become dominated by errors and the DPC fails. In these cases, a suitable regularization term should fulfill the DPC, therefore, it facilitates to locate before the point where the data becomes dominated by errors.

To confirm that Picard’s condition was fulfilled, we plotted on the same graph the $\left|u_i^{T}d\right|$ coefficients, their corresponding singular value, and the quotient of both by applying the Picard routine available in the Matlab Regularization Toolbox [38].

3. Results

3.1 Phantom

Figure 3 shows the U-curve plot on a log-log scale for the phantom experiment.

 

Fig. 3 U-curve plots on log-log scale. (Minimum which corresponds to ${\alpha }_u$ = 4.38*10⁻²).

Download Full Size | PDF

In this case, the U-curve shows a minimum which corresponds to a regularization parameter ${\alpha }_u$ = 4.38*10⁻².

Figure 4 shows the L-curve plot on a log-log scale provided using [38]:

 

Fig. 4 L-curve plots on log-log scale. (Minimum which corresponds to ${\alpha }_l$ = 5.65*10⁻²).

Download Full Size | PDF

The L-curve did not exhibit a neat corner. The failure to find a sharp corner was due to the high ill-posedness of this problem.

Figure 5 depicts the fDOT reconstruction with α parameters in the 10⁻¹ to 10⁻⁶ range and with ${\alpha }_u$ = 4.38*10⁻².

 

Fig. 5 Slices in the y-z plane of the reconstructions obtained for the α parameter in the 10⁻¹ to 10⁻⁶ range. The result for ${\alpha }_u$ =4.38*10⁻² (obtained from the U-curve) is showed at the bottom center. At the bottom right: slice indicating the phantom fluorescence concentration.

Download Full Size | PDF

Figure 6(a) shows profiles taken in the x direction (corresponding to z=7.5 mm and y=7.5 mm) and the FWHM (Full width half maximum) versus noise plot, outlining the behavior of the resolution and the image noise of the reconstructions depending on the regularization parameter. It has been shown in [15] that the general trend is for resolution to increase together with image noise while the regularization parameter decreases, and this trend can be seen clearly in Fig. 6. For α values of 10⁻¹, 4.38*10⁻² (${\alpha }_u$), and 10⁻², the FWHM of the profiles decreases. For α = 10⁻³ the trend is truncated because image noise starts to prevail, while for α<10⁻³ the reconstruction is noise only, and the object is no longer visible in the reconstructed images. Taking these considerations into account, we observe a heuristic range of ${\alpha }_u$ values that produces reconstructed images with a reasonable amount of noise and resolution. This range includes the optimum value obtained by the U-curve method, namely, ${10}^{-1}\le \alpha <{10}^{-3}$.

In Figs. 7(a) and 7(b), we plot the noisy SVD components of the solution and the right-hand side of the phantom study. One interesting aspect is the severe ill-posedness of the problem, which is seen by the fact that the singular values decay exponentially (Fig. 7(c)).

 

Fig. 7 Plots of the SVD components of the phantom and Picard’s plot. (a) Noisy SVD components of the solution. (b) Noisy SVD components of the right-hand side. (c) Decay of the singular values. The regularization parameter provided by U-curve method (${\alpha }_u$ =4.38*10⁻²) is plotted as a horizontal dashed blue line. (d) Picard’s plot with the maximum and minimum parameter of the heuristic acceptable range plotted as two horizontal dashed black lines (10⁻¹ and 10⁻³).

Download Full Size | PDF

Figure 7(d) illustrates Picard’s plot showing the maximum and minimum parameter of the heuristically acceptable range plotted as two horizontal dashed lines (10⁻¹ and 10⁻³).

Picard’s plot makes it possible to compare the SVD coefficients of the right-hand side with the singular values and their quotient. Singular values below 10⁻³, on average, decay to zero faster than those of the respective $\left|u_i^{T}d\right|$coefficients.

The blue line is the decay of the singular values ${\sigma }_i$, the green crosses correspond to $\left|u_i^{T}d\right|$, and the red circles represent the quotient $\frac{\left|u_i^{T}d\right|}{{\sigma }_i}$.

We can observe that the singular values above the heuristic acceptable range of parameters 10⁻¹ and 10⁻³, particularly the singular values above the U-curve cut-off (4.38*10⁻²), fulfill Picard’s condition.

3.2 Ex-vivo experiments

Figure 8 shows the U-curve plot on a logarithmic scale with a minimum which corresponds to the regularization value ${\alpha }_u$ =5.72*10⁻².

 

Fig. 8 U-curve plot on the log-log scale. (Minimum which corresponds to ${\alpha }_u$ = 5.72*10⁻²).

Download Full Size | PDF

The curve is not really U-shaped, showing that fewer singular values remain than in the phantom experiment.

Figure 9 depicts the fDOT reconstruction with α parameters in the 10⁻¹ to 10⁻⁶ range and with ${\alpha }_u$ =5.72*10⁻².

 

Fig. 9 Slices in the y-z plane of the 3D render of the reconstructions obtained for the α parameter in the 10^-1 to 10^-6 range. The result for ${\alpha }_u$= 5.72*10^-2 (obtained from the U-curve) is showed at the bottom center. Dye concentration and volume are indicated at the bottom right.

Download Full Size | PDF

Figure 10 shows a coronal view of a 3D render of the Tikhonov reconstruction in terms of SVD for the regularization parameter obtained using the U-curve method. The reconstruction is merged with the white light image of the mouse.

 

Fig. 10 Coronal view of a 3D render of the reconstruction for ${\alpha }_u$ = 5.72*10⁻² (obtained from U-curve). The white light image is shown behind as a reference image. Dye concentration and volume are indicated on the image.

Download Full Size | PDF

Figure 11(a) shows profiles taken in the x direction, corresponding to z=8 mm and y=0.7 mm, and the FWHM (resolution) vs. noise plot, which outlines the behavior of the resolution and image noise of the reconstructions depending on the regularization parameter. Again, we observe a range of α values that produces reconstructed images with a reasonable amount of noise and resolution. The U-curve-based value fall within this range, which is ${10}^{-1}\le \alpha \le {10}^{-2}$.

In this case, we can say that the U-curve gives minimum noise while retaining the best resolution possible.

Similar to Fig. 7(d), Fig. 12 shows that the singular values above the heuristic acceptable range of parameters, particularly the singular values above the U-curve cut-off, fulfill Picard’s condition.

 

Fig. 12 Picard’s plot with the maximum and minimum parameter of the heuristic acceptable range plotted as two horizontal dashed black lines (10⁻¹ and 10⁻²).

Download Full Size | PDF

The blue line is the decay of the singular values ${\sigma }_i$, the green crosses correspond to $\left|u_i^{T}d\right|$, the red circles represent the quotient $\frac{\left|u_i^{T}d\right|}{{\sigma }_i}$, and the two horizontal dashed lines represent the heuristic acceptable range (10⁻¹ and 10⁻²).

We can observe that, as in the case of Fig. 7(d), the singular values above the heuristic acceptable range of parameters 10⁻¹ and 10⁻², particularly the singular values above the U-curve cut-off (5.72*10⁻²), fulfill Picard’s condition. The singular values decay faster than the singular values of Figs. 7(c)–7(d) (for phantom experiment), therefore the problem for the ex-vivo experiment is more ill-posed than the problem for the phantom experiment.

4. Discussion and conclusion

In this paper a U-curve-based method is utilized for the first time to select the regularization parameter in Tikhonov regularization reconstruction of fDOT. Suitable selection of the regularization parameters, in terms of Picard’s condition, image resolution and image noise, has been provided by the U-curve. Results are shown both on phantom and mouse data.

Choosing the correct regularization parameter is crucial for the reconstruction of DOT and fDOT data. Singular Value Analysis has been used to optimize experimental setups in optical tomography [12,40–43], and the Tikhonov regularization has recently been used to balance the effect of anatomical and optical information in fDOT problems based on a priori information. Therefore, an automatic method that enables us to choose the regularization parameter is paramount.

To our knowledge, the L-curve method is the only automatic strategy which does not need a priori knowledge of the noise and that has been successfully applied to fDOT [1,16,17].

Recently, it has been found in other fields that in practice methods without a-priori information, such as L-curve and GCV were robust [20,21].

We wish to recall that the aim of this paper was not to perform a comparison of these methods. Nevertheless, we would like to point out that the L-curve method presents several theoretical limitations and may fail to find a good regularization parameter when the solutions are very smooth [44], and examples of inverse problems where the L-curve does not converge have been found [45].

In the diffuse optical tomography case, the authors of [15] emphasized that the L-curve analysis yielded overly-smooth solution in same cases.

The GCV method, on the other hand, is more computationally expensive for large systems than the L-curve [22].

Only two experiments are presented in this study. However, when we used the U-curve method in other experiments with different aims [40,41], we always obtained satisfactory reconstructions, both for mice and for phantoms.

It can be seen that the L-curve calculated for the phantom experiment did not exhibit a neat corner (Fig. 4). However, the U-curve for the same experiment had a pronounced minimum (Fig. 3). Furthermore, this minimum was found in the interval given in the section 2.5 (thus not being necessary to calculate the U-curve for the α parameters out of this interval). Thus, this is a case that shows the robustness and the computation effectiveness of the U-curve method.

Besides, Figs. 7(d) and 12 (Picard’s plots) show how the U-curve regularization parameter satisfies Picard’s condition, thus assuring a satisfactory regularized solution. In Figs. 5 and 9 (reconstruction obtained for different α parameters) and Figs. 6 and 11 (profiles and FWHM versus noise plot), we can choose values for the regularization parameter that are lower than the value at which the reconstructed image started to be noisy. Figures 7(d) and 12 demonstrated that, for these lower values, Picard’s condition is satisfied, as the U-curve parameter is in this range.

Simple observation of the Picard’s plot can reveal a valid regularization parameter; however, the choice could prove to be more subjective, as it would depend on the visual observation of the user. An automatic selection of the threshold parameter may be simpler and more objective in most cases.

Furthermore, in agreement with [15], we can see clearly in Fig. 6(b) (FWHM versus noise plot) that the general trend is for resolution to increase together with image noise while the regularization parameter decreases.

We want to point out that for the ex-vivo experiment, the resolution is better (Fig. 11(b) versus Fig. 6(b)), due to the fact that the tip of the capillary is closer to the surface than in the phantom experiment. As the resolution of DOT systems are depth-dependent, resolution improves the closer the object is to either side of the slab [46].

Some advantages of our work are the following: First, we believe that the U-curve method may constitute a good alternative in the cases above mentioned, where the L-curve yields unsatisfactory results. Second, according to [23–25], we have found that an interval exists where the U-curve has a minimum (the optimal regularization parameter exists) and this fact can greatly increase the computational efficiency in selecting the regularization parameter.

We expect the automatic U-curve method of selection of regularization parameter to yield robust and useful results that can be applied to the reconstruction of fDOT images and studies of image performance by singular value analysis.