2.1.

Experimental Setup

Most of the elements of the optical setup (illustrated on Fig. 1) used for this study are classical and are based on a common setup for reflectance molecular imaging. The light source (noted 1 on Fig. 1) is a 690-nm fibered laser (Intense HPD model 7404) which illuminates a tissue-like liquid phantom. A $690/10\mathrm{nm}$ clean up filter was used (noted 2 on Fig. 1). Fluorescence images are acquired with a CCD camera (PCO Pixelfly VGA, $640\times 480\text{pixels}$ images, noted 6 on Fig. 1) for each position of the object which rests on a motorized translation stage (noted 4 on Fig. 1). A cylindrical lens (noted 3 on Fig. 1) is used to focus the laser at 45 deg on the phantom along a line of length 4 cm and width 1 mm which sets the translation steps at 1 mm to fully illuminate the phantom.

Fig. 1

Optical setup used during the study [(a) schematic and (b) platform].

For this study, the excitation line is static and the object observed is translated. While this is not a problem for this proof of principle where only phantoms were studied, a better implementation of the setup would keep the object observed static and we would move the excitation line. This would indeed be more practical, in particular for in vivo imaging.

The optical setup was optimized to minimize the amount of excitation bleed-through by using proper filters. A fluorescence filter (Semrock Razoredge 808-nm long pass filter, noted 5 on Fig. 1) is in front of the camera to stop all excitation photons so as to detect a fluorescence signal. The laser power is 15 mW.

For all experiments presented here, the stack of fluorescence images was first acquired and then the stack of excitation images was acquired. For the acquisition of excitation images, the fluorescence filter is replaced manually with a neutral density filter. The fact that we switch filters manually does not cause any significant distortion between fluorescence and excitation images for these experiments, where only still optical phantoms are observed. To image living subjects, we should, however, be able to acquire fluorescence and excitation images at the same time to avoid distortions caused by movements of the subject.

For the contrast enhancement experiments, 80 images were acquired for each depth considered, with integration times ranging from 30 ms at 1 mm to 3 s at 1 cm. The acquisition was automated so that 1 s passed between each acquisition to ensure that the liquid phantom was completely still and did not move because of the translation. This leads to total scanning times ranging from 82.4 s at 1 mm to 320 s at 1 cm. No significant photobleaching was observed for any of the experiments.

For the resolution enhancement experiments, we used $100\text{-}\mu \mathrm{m}$ translation steps to be able to use thinner detection stripes. 400 images were acquired for each depth considered, with integration times ranging from 30 ms at 1 mm to 250 ms at 4 mm. This leads to total scanning times ranging from 412 s at 1 mm to 500 s at 4 mm. No significant photobleaching was observed for any of the experiments.

All the liquid phantoms used in this study were made with the same recipe to obtain tissue-like optical properties with an absorption coefficient ${\mu }_{\mathrm{a}}=0.05{\mathrm{cm}}^{-1}$ and a reduced scattering coefficient ${\mu }_{\mathrm{s}}^{\prime }=10{\mathrm{cm}}^{-1}$. The fluorescent inclusions are 1-mm diameter tubes containing $3\mu \mathrm{M}$ of indocyanine green encapsulated in lipid nanoparticles (LNP-ICG²¹) diluted in the same preparation as the phantoms to match the background optical properties. We added different amounts of LNP-ICG to the phantom preparation to obtain different ratios of fluorescence to background fluorescence. The LNP-ICG absorption and fluorescence emission spectra are presented on Fig. 2(d).

Fig. 2

(a) Fluorescence target used for the resolution study; (b), (c) phantoms used for the absorption correction study; and (d) indocyanine green encapsulated in lipid nanoparticles (LNP-ICG) fluorescence excitation and emission spectra.

Four different types of phantoms were used for this study:

2.2.

Image Processing Methods

2.2.2.

Method proposed

The processing method we will present relies on the assumption that there is a relationship between the excitation and the autofluorescence signals defined as following:²²

Eq. (1)

$$\frac{A(r)}{E(r)}=\beta r,$$

where $E$ is the excitation, $A$ is the autofluorescence, $r$ is the distance, and $\beta$ is a proportionality coefficient which depends on the optical properties of the medium and the background fluorescence yield. The geometry and intensity profiles are depicted in Fig. 3.

Fig. 3

Schematic of the geometry studied: a homogeneous fluorescent medium is considered and we look at the excitation (red curves) and autofluorescence (blue curves) intensity profiles with respect to the distance $r$ from the excitation line.

If the assumption is true, knowing the excitation profile could give us an insight on the autofluorescence profile. Simulations and experimental validations were performed to test this relationship, with simulations’ parameters for the optical properties chosen to match the experimental conditions.

Two types of simulations were performed to check the validity of our hypothesis: first with the Monte-Carlo method (with ${10}^7$ photons launched) and also with the NIRFAST software^23,24 which is based on the diffusion approximation.

In both simulations, the medium is a homogeneous fluorescent one in a slab geometry with tissue-like optical properties (${\mu }_{\mathrm{a}}=0.05{\mathrm{cm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=10{\mathrm{cm}}^{-1}$, $g=0.9$).

To experimentally validate our hypothesis, measurements were acquired on a homogeneous autofluorescent phantom. Fluorescence and excitation images were acquired. As with the simulations, we studied the profiles ratios $A/E$ for several increasing levels of autofluorescence.

The results obtained both with simulations and experimentally are presented in Fig. 4.

Fig. 4

$A/E$ intensity ratios obtained for linearly increasing levels of autofluorescence with the: (a) Monte-Carlo method (20 levels of autofluorescence considered), (b) with NIRFAST (10 levels of autofluorescence considered), and (c) experimentally (10 levels of autofluorescence considered).

$A/E$ ratio profiles are presented for increasing levels of autofluorescence. We see the linear relationship of $A/E$ with $r$ as expected in all cases. For experimental data and a distance superior to 30 mm, the noise due to lack of photons becomes dominant. Subsequent experimental analyses will, therefore, only be done for distances to the source inferior to 20 mm.

We will now explain how we take advantage of this relationship between the excitation and autofluorescence signals to reduce the effects of autofluorescence.

Let us consider that for one excitation position $i$ and one column of the images $j$:

$$\begin{gathered} I_{\mathrm{tot}}(r,i,j)\approx F(r,i,j)+\alpha (i,j)E(r,i,j)+A(r,i,j), \\ {I}_{\mathrm{tot}}(r,i,j)\approx F(r,i,j)+\alpha (i,j)E(r,i,j)+\beta (i,j)E(r,i,j)r, \end{gathered}$$

where $I_{\mathrm{tot}}$ is the total signal detected on the camera, $F$ is the fluorescence signal of interest, $A$ is the autofluorescence, $E$ is the excitation leaks, and $\alpha$ and $\beta$ are the coefficients depending on the excitation position and the column considered for the intensity profiles.

By fitting $I_{\mathrm{tot}}$ with $\alpha E+\beta Er$, we can find the fluorescence signal of interest $F$ by subtraction.

The fitting method is the following: for each excitation position, we acquire both a fluorescence and an excitation image. For each position, we then fit the fluorescence intensity profile in each column of the image (480 columns with our detector) by using the corresponding excitation profile with the method of least squares. We then obtain $\alpha$ and $\beta$ parameters specific to each position. These parameters stay the same for excitation positions, where only the autofluorescence is excited, and have strong variations for positions surrounding the fluorescence target.

We studied two possibilities: we can either choose to use the parameters specific to each position (referred as “local parameters”) or to use the mean parameters (referred to as “global parameters”). As we will see with the following simulations made with NIRFAST, each possibility has its own advantages and drawbacks.

In the first case studied, we used the same homogeneous fluorescent medium as in the previous simulation and we added a fluorescent inclusion located 2 mm below the surface. The simulated phantom and the different intensity profiles of interest are depicted in Fig. 5.

Fig. 5

Intensity profiles comparison: (a) blue line, raw fluorescence $F$; dotted blue line, autofluorescence; red line, local parameters fit of the fluorescence $M_1$; green line global parameters fit of the fluorescence $M_2$; (b) blue line $F$; red line, $(F-M_1)$; green line $(F-M_2)$.

As expected, the intensity profile obtained when fitting with the local parameters (noted $M_1$, red curve on the left of Fig. 5) is close to the total fluorescence profile and is sensitive to the fluorescent inclusion, while the one obtained when using the global parameters (noted $M_2$, green curve on the left of Fig. 5) is closer to the autofluorescence profile and is not sensitive to the fluorescent inclusion.

After subtraction, both methods completely nullify the autofluorescence contribution. The global parameters’ method leads to a better dynamic around the fluorescence inclusion, but suffers from boundary effects for distances larger than 40 mm from the inclusion in this case.

For the second case studied, we used the same medium as in the previous simulation and we added five fluorescent inclusions located 2 mm below the surface. The simulated phantom and the different intensity profiles of interest are depicted in Fig. 6.

Fig. 6

Intensity profiles comparison: (a) blue line, raw fluorescence $F$; dotted blue line, autofluorescence; red line, local parameters fit of the fluorescence $M_1$; green line, global parameters fit of the fluorescence $M_2$; (b) blue line, $F$; red line, $(F-M_1)$; green line, $(F-M_2)$.

In the previous case, the parameters $\alpha$ and $\beta$ were only varying for the excitation positions surrounding the fluorescent inclusion. This is why taking their mean value allowed to get close to the value of excitation positions where the background signal was the only signal detected. In this case, the parameters $\alpha$ and $\beta$ also vary only for excitation positions surrounding the fluorescent inclusions. However, there are five inclusions in this case, while the excited field has the same size. As the number of excitation positions having parameters which depend on the fluorescent inclusions is higher, taking their mean value leads to an overestimation of the autofluorescence. This is why after subtraction, we obtain a negative level of signal with the global parameters.

Figure 7 presents how our model behaves experimentally for five different positions of excitation around a fluorescent inclusion. On the first line of this figure, we can see the fluorescence images and the two columns that will be used for the example: column 1 is located over the fluorescent inclusion and column 2 is located over an area where only the autofluorescence signal is detected. On the second line of the figure, the fluorescence intensity profiles and their fits are plotted for column 1 at the five positions considered. On the third line, the fluorescence intensity profiles and their fits are plotted for column 2 at the five positions considered.

Fig. 7

Example of fits obtained for five different positions of the excitation line around a fluorescent inclusion; (a) fluorescence images with the two columns used in this example, column 1 is above the fluorescent inclusion and column 2 is above the autofluorescence only; (b) fluorescence intensity profile and its fit for column 1 for the five positions; and (c) fluorescence intensity profile and its fit for column 2 for the five positions.

On this last line, where only the autofluorescence is considered, we see that the intensity profile and its fit overlap nearly perfectly. When subtracting its fit from the fluorescence intensity profile, the contribution from the autofluorescence is then completely suppressed.

On the other hand, on the second line where some of the fluorescence of interest is detected, some differences appear between the fluorescence intensity profile and its fit. These differences become more and more important as the excitation line gets closer to the fluorescent inclusion. In this case, when the fit is subtracted from the fluorescence intensity profile, the autofluorescence contribution is suppressed but the fluorescence of interest contribution remains.

3.1.

Contrast Enhancement

To compare the different detection schemes and quantify the improvements, we introduce the contrast $C_{T,N}$ defined as

The regions of interest were the same for both WF-FRI and our method. The neutral region of interest was chosen in the top edge of the images so that it was the farthest possible from the target region of interest while still being illuminated the same way, to ensure that the inhomogeneities of the illumination could not bias the results.

We will first present how our method improves the contrast compared to the classical WF-FRI in the case of a phantom with a single fluorescent inclusion at different increasing depths. Two different levels of background fluorescence were considered: a realistic one and a stronger one to test the limits of the method. The realistic level of background fluorescence has been chosen using Ref. 25. We have used a concentrations’ ratio of about 80 between the fluorescent inclusion and the background medium (this is based on the “ICG equivalent” signal of skin, which is the median one compared to the different organs presented in Ref. 25). For the stronger level of background fluorescence, we used a ratio of fluorescence to background fluorescence five times lower. We then had an “ICG equivalent” autofluorescence signal 1.5 times higher than the worst case described in Ref. 25. This extreme case allows us to test the limits of our method and could be compared to a case right after the fluorescence dye is injected in the patient, leading to a strong homogeneous nonspecific fluorescence signal which would be added to the existing background fluorescence signal.

In Fig. 8(a), we compare the contrasts observed at the ten depths considered with the WF-FRI and both fitting methods for the realistic level of background fluorescence.

Fig. 8

Comparison between the contrasts obtained with WF-FRI (blue line), the local parameters fitting method (red line) and the global parameters fitting method (green line) (a) for a realistic background fluorescence level and (b) a stronger background fluorescence level used to test the limits of the method.

Contrast is enhanced with both the local and global parameters fitting methods. In this simple case where there is a single fluorescent inclusion, the global parameters method offers better performance as mentioned before: with the global parameters method, the parameters used are the means on every excitation position. By doing this, the “fast” fluorescence variations corresponding to fluorescence inclusions or absorption heterogeneities are lost, but we are closer to the “slow” fluorescence variations corresponding to autofluorescence. This is why after subtracting the fit to the fluorescence image we obtain a very good contrast as we subtract only the autofluorescence contribution. In this case, the gain obtained with the local parameters varies between 1.2 at 1 mm and 3.1 at 10 mm while it varies between 1.4 and 4.4 with the global parameters.

In Fig. 8(b), we present the results in the same case as before, but with a stronger background fluorescence signal.

Even in the case with a stronger background signal, our fitting method enhances the contrast compared to the WF-FRI with both local and global parameters. However, the gain obtained is smaller than the one obtained in the case with less background signal. With the local parameters fitting method, the gain is low, with a maximum of 1.5 at 10 mm. With the global parameters fitting method, the gain varies between 1.3 at 1 mm and 2.3 at 10 mm.

3.2.

Resolution Enhancement

The second set of results was obtained with the fluorescent resolution target described in Sec. 2.1. The concentration of background fluorescence was set to have a realistic fluorescence to background ratio (as described before), and the target was immersed at three depths between 1 mm and 3 mm. The images obtained are presented in Fig. 9.

Fig. 9

Images of the resolution target at three different depths obtained with (a) WF-FRI, (b) the local parameters fitting method, and (c) the global parameters fitting method.

At 1 mm, it is possible to resolve the target with WF-FRI (a) but there is already some crosstalk between the different groups of inclusions leading to an overestimation of the signal produced by the largest inclusions.

The local parameters fitting method (b) offers the best performance. It enhances the resolution as it increases the peak to valley ratio between the fluorescent inclusions and the background, and the background fluorescence surrounding the inclusions is totally suppressed. The global parameters’ fitting method (right column) also suppresses the background fluorescence surrounding the inclusions, but the peak to valley ratio between the fluorescent inclusions and the background is slightly lower than with the local parameters’ fitting method.

At 2 mm, it becomes difficult to properly resolve the target with WF-FRI. Signals coming from the three largest groups of inclusions start to overlap and form one large fluorescent signal, making it impossible to distinguish the single inclusions in some groups.

The local parameters’ fitting method still offers the best performance as the resolution between the groups of inclusions is still good, but it becomes harder to see the inclusions in the group of inclusions with the smallest diameter. The global parameters’ fitting method gives results comparable to the ones obtained with the local parameters fitting method, but, as seen at 1 mm, the peak to valley ratio between the fluorescent inclusions and the background is slightly lower.

At 3 mm, the whole target only emits one large fluorescent signal with WF-FRI and the groups of inclusions are completely unresolved. With both fitting methods, the target is unresolved and we can only see one large signal due to the three largest groups of inclusions and the two smallest signals due to the two smallest groups of inclusions.

To summarize these results, the local parameters’ fitting method is better suited to enhance the resolution then the global parameters’ fitting method. As mentioned before, with the local parameters, the parameters used are the ones specific to each position. Contrary to the mean parameters, we do not lose all the “fast” fluorescence variations which is why the subtraction of the fit leads to a better resolution and a correction of some effects of the absorption heterogeneities.

However, the main interest of these fitting methods is the reduction of the background fluorescence signal.

3.3.

Effect of the Absorption

We wanted to test the robustness of our method in conditions where absorption heterogeneities are present. Two simple yet interesting cases were studied: a first one where an absorption heterogeneity is located right above a fluorescent inclusion, giving the impression that two fluorescent inclusions are present rather than one, and a second one where two fluorescent inclusions with the same fluorescence but different absorptions give different signals. The aim is not to fully correct the absorption effects but rather to test the limits of our method. In this section, we compare the results obtained with our fitting methods to the WF-FRI as previously, but also to the fluorescence with the excitation ratio method presented in Ref. 26

The first phantom used for this study is the one described in Fig. 2(b). A fluorescent inclusion is located 2 mm below the surface of the phantom and an absorption inclusion is right above it. Four different absorption coefficients were used for the absorption inclusion: $0.05{\mathrm{mm}}^{-1}$ (matching with the surrounding absorption of the phantom), 0.1, 0.2, $0.5{\mathrm{mm}}^{-1}$. In Fig. 10, we show the images obtained with the WF-FRI, both fitting methods, and the fluorescence by excitation ratio method for the control experiment without heterogeneity and for the experiment with the highest absorption coefficient considered.

Fig. 10

Images obtained with the different methods for two configurations: (a) control experiment, the heterogeneity has the same absorption coefficient as the surrounding phantom; and (b) the heterogeneity has an absorption coefficient (${\mu }_{\mathrm{a}}=0.5{\mathrm{cm}}^{-1}$) 10 times higher than that of the surrounding phantom.

We see in the WF-FRI images the influence of the absorption heterogeneity: due to the decrease of the fluorescence signal caused by the heterogeneity, it seems that there are two fluorescence sources rather than a single one. The same behavior is observed for both fitting methods. However, the fluorescence by excitation ratio method corrects this effect as the signal detected is the same for the control experiment and the heterogeneity experiment.

These remarks are best visualized in Fig. 11 where we have plotted the intensity profiles taken perpendicularly to the capillary for the different methods. The profiles for the four different absorption coefficients are shown.

Fig. 11

Fluorescence intensity profiles obtained with the different methods for four absorption coefficients for the heterogeneity: $0.05{\mathrm{mm}}^{-1}$ (blue line), $0.1{\mathrm{mm}}^{-1}$ (red line), $0.2{\mathrm{mm}}^{-1}$ (green line), $0.5{\mathrm{mm}}^{-1}$ (black line). (a) WF-FRI, (b) local parameters fitting method, (c) global parameters fitting method, and (d) fluorescene by excitation ratio.

With WF-FRI, we see that the drop in intensity becomes more and more visible as the absorption coefficient of the heterogeneity increases. As said previously, for the highest absorption considered, two different peaks begin to appear.

With the global parameters’ fitting method, no improvement is observed, thus the behavior is comparable with the one observed with WF-FRI.

With the local parameters’ fitting method, we see that the two peaks caused by the heterogeneity are still visible and are not corrected. However, we see that the intensity decrease due to the absorption is partly corrected: with WF-FRI, the signal decreases at 83% of the control experiment with the ${\mu }_{\mathrm{a}}=0.1{\mathrm{mm}}^{-1}$ heterogeneity, 71% with the ${\mu }_{\mathrm{a}}=0.2{\mathrm{mm}}^{-1}$ heterogeneity, and 57% with the ${\mu }_{\mathrm{a}}=0.5{\mathrm{mm}}^{-1}$ heterogeneity. With the local parameters’ fitting method, the signal stays at 100% of the control experiment with the ${\mu }_{\mathrm{a}}=0.1{\mathrm{mm}}^{-1}$ heterogeneity, decreases at 75% with the ${\mu }_{\mathrm{a}}=0.2{\mathrm{mm}}^{-1}$ heterogeneity, and 75% with the ${\mu }_{\mathrm{a}}=0.5{\mathrm{mm}}^{-1}$ heterogeneity.

With the fluorescence by excitation ratio method, we see that the two peaks caused by the heterogeneity are corrected. We also see that the intensity decrease due to the absorption is partly corrected: the signal decreases at 88% of the control experiment with the ${\mu }_{\mathrm{a}}=0.1{\mathrm{mm}}^{-1}$ heterogeneity, 75% with the ${\mu }_{\mathrm{a}}=0.2{\mathrm{mm}}^{-1}$ heterogeneity, and 60% with the ${\mu }_{\mathrm{a}}=0.5{\mathrm{mm}}^{-1}$ heterogeneity.

To review these first results, the local parameters’ fitting method seems to be suited to correct some of the effects of absorption heterogeneities, namely the decrease of the intensity detected. The results are slightly better than the fluorescence by excitation ratio, but this method also manages to correct the shape changes of the signal detected which are caused by the heterogeneities: there is only one flat peak compared to the local parameters’ fitting method which gives two peaks.

The second phantom used for this study is the one described in Fig. 2(c). Two fluorescent capillaries are located at the same depth, the distance between them is 1 cm. Two experiments were performed: one where both fluorescent inclusions have the same absorption coefficient as the background, and another one where one of the inclusions has an absorption coefficient ten times higher (the top capillary on Fig. 12). The images obtained with the WF-FRI, both fitting methods, and the fluorescence by excitation ratio method are shown in Fig. 12.

Fig. 12

Images obtained with the different methods for two configurations: (a) both capillaries have the same absorption coefficient; and (b) the top capillary has an absorption coefficient 10 times higher.

As expected, the higher absorption of the top capillary leads to a decrease of its signal intensity: it is twice as low as the one of the capillary having the same absorption as the background. None of the techniques seem to be able to correct the effects of absorption.

To have a quantitative idea of the absorption correction of the different techniques, we have plotted in Fig. 13 the intensity profiles taken perpendicularly to the capillaries for the different methods.

Fig. 13

Fluorescence intensity profiles obtained with the different methods for two configurations: both capillaries have the same absorption coefficient (blue line), the left capillary has an absorption coefficient 10 times higher: (a) WF-FRI, (b) local parameters fitting method, (c) global parameters fitting method, and (d) fluorescene by excitation ratio.

First, we can notice that the local parameters’ fitting method enhances the resolution compared with other methods, as was shown in the previous section. It also appears that no method totally corrects the effects of absorption. The only method that increases the signal of the higher absorption capillary is the local parameters’ fitting method: the signal is at 63% of its original level, contrary to the other methods where it is at 50%.

As in the previous case, the local parameters’ fitting method is suited to correct some of the effects of absorption heterogeneities, even in this case where the fluorescence by excitation ratio has no beneficial effect.