2.1.

Determination of the Point-Spread Function

Many techniques have been employed to account for the instrumental response of the imaging system¹¹ and to obtain the PSF of a camera system.^{12, 13, 14} The most direct means of measuring the PSF is to take images of a point source. The CCD camera system used in this study was a TE/CCD 512 EFT (Princeton Instruments, Trenton, New Jersey) with a $16\text{-}\text{bit}$ analog-to-digital converter and an array of $512\times 512\text{pixels}$ corresponding to an area of the image of approximately $37\times 37{\mathrm{mm}}^2$ . The collection optics consisted of a commercial camera zoom lens (Nikkor, $50\text{-}\mathrm{mm}$ focal length, Nikon, Japan). The aperture of the camera was adjusted to a medium level. The camera chip was cooled by a Peltier cooling to $-35°\mathrm{C}$ to minimize dark noise. Additionally, we used water cooling for the camera.

Figure 1 shows the experimental setup. The beam of a HeNe laser emitting at $633\mathrm{nm}$ (LHRR-0100, Laser 2000, Germany) was coupled into an optical single-mode fiber with a core diameter of $\sim 10\mu \mathrm{m}$ . The fiber served as a point source so that the light of the laser appears in the center of a CCD pixel. The exposure time could be varied by means of a shutter. To minimize stray light, the experiment was performed in a dark room. Unwanted light strayed and reflected by parts of the experimental setup was carefully blocked with a black cloth. The end of the fiber was positioned in the image plane of the camera system. With this setup, the strongest signal generated by the laser on the CCD sensor occupied just one pixel. Before each measurement, a dark image was recorded and subtracted.

Fig. 1

Experimental setup for the determination of the PSF. The dashed line illustrates a beam scattered by the objective of the camera system, causing the PSF not to be a delta function.

To simplify the characterization of the PSF, we assumed that it is radially symmetric across the sensor. The exposure was divided into concentric rings, and pixel values in each ring were summed and divided by the number of pixels in the ring. Pixel values from different areas but with the same radial distance from the source were averaged. Radial averaging helps to increase the signal-to-noise ratio: As the signal decreases for larger distances to the illumination point, the number of pixels averaged increases because of the increasing radius of the ring. This leads to an improvement of the signal-to-noise ratio.

2.2.

Experimental Setup for Phantom Measurements

To experimentally examine the effect of the PSF on our CCD camera system and for the proper assessment of the experimental setup, measurements on tissue phantoms were performed. The phantom should mimic biological tissue and allow easy tuning of the optical properties. Several publications concerning the production of phantoms have been reviewed.^{15, 16} Studies of Lipovenös (Fresenius Kabi, Germany), an intravenous nutrient containing 20% soy oil, or similar fat solutions described measurements of its scattering properties and found that the absorbance can be neglected in many cases.^{17, 18} As an absorber, we used trypan blue (Biochrom, Germany), an organic dye, due to its negligible scattering. We determined the absorption properties of undiluted trypan blue by measuring the collimated transmission. Both substances were mixed with distilled water. Changing the concentrations of the added substances, the optical properties of the obtained solutions could be varied. We made up the tissue phantoms in a cylindrical black bin with a height of $\sim 18\mathrm{cm}$ and a diameter of $\sim 10\mathrm{cm}$ , so the phantoms could be regarded to be a semi-infinite half-space. The index of refraction of the phantom was assumed to be $n=1.33$ , because the phantom consists almost entirely of water. An analytic expression for the reflectance derived by means of the diffusion model was used for the determination of the optical properties.¹⁹

The experimental setup is provided in Fig. 2 . Again the same HeNe laser emitting at $633\mathrm{nm}$ served as a light source. The laser beam was deflected onto the specimen by a mirror at an angle of incidence of $10\mathrm{deg}\text{to}15\mathrm{deg}$ to avoid the detection of specularly reflected light. We confirmed by Monte Carlo simulations that this small angle had a negligible influence on the spatially resolved reflectance. The diffusely reflected light was detected by the same CCD camera mounted perpendicular to the specimen surface as described in Sec. 2.1. The distance between the front of the camera and the sample was fixed at $288\mathrm{mm}$ . Neutral density filters were used to adjust the incident laser power to avoid saturation of the detector system. Again, we assumed that the reflectance is radially symmetric, and thus, we divided the exposure into concentric rings. Pixel values in each ring were summed and divided by the number of pixels in the ring to improve the signal-to-noise ratio.

Fig. 2

Experimental setup for phantom measurements.

When directing the beam into the middle of the objective of the camera and fixing the distance between the laser and the camera to $288\mathrm{mm}$ , as mentioned in Sec. 2.1, we were able to determine the diameter of the laser beam. Referring to an intensity of $(1∕e)\approx 37\%$ of the maximum value, we received a beam diameter of $\sim 0.15\mathrm{mm}$ .

2.3.

Data Analysis

Since the PSF is not just a delta function, effects of the PSF on the intensity values far away from the brightest pixel cannot be ignored without appreciable radiometric consequences. In a real optical system, images obtained are actually the convolution of the pixel intensity values with the normalized PSF of the system. We carried out a two-dimensional (2-D) convolution algorithm based on fast Fourier transform (FFT) and inverse FFT techniques. The Fourier transform of the convoluted raw image of the camera, ${\overline{I}}_{\mathit{\text{raw}}}$ , equals the product of the Fourier transforms of the normalized PSF, ${\overline{I}}_{\mathit{PSF}}$ , and the true image, ${\overline{I}}_{\mathit{\text{true}}}$ ,

The absorption and the reduced scattering coefficients were obtained by (1) fitting the experimental data to an analytical solution of the diffusion equation for a homogenous infinite turbid medium according to Kienle and Patterson,¹⁹ and (2) fitting the 2-D convolution of the preceding solution of the diffusion equation and of the PSF to the experimental data. We used a convolution algorithm because a deconvolution process tends to increase the noise level in the image and so enhances the contrast of both the signal and the noise. The extrapolated boundary condition was regarded. The fitting range contained only distances larger than $1∕{\mu }_s^{\prime }$ because the diffusion equation is not valid for distances shorter than one transport mean free path.¹⁹ It included three orders of magnitude from this starting distance. For large ${\mu }_s^{\prime }$ , the fitting range included only distances larger than $0.8\mathrm{mm}$ because of the finite beam diameter used in the experiments. If the optical parameters of the phantoms were roughly known, the starting parameters were chosen to vary at about 30 to 40% from the expected values. If the optical parameters were completely unknown, we started with ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1.0{\mathrm{mm}}^{-1}$ . As the obtained optical parameters depend on the fitting range and the fitting range—as described earlier—depends on the starting value of ${\mu }_s^{\prime }$ , the fitting algorithm sometimes had to be accomplished twice or triply. This was especially the case if the starting parameter of ${\mu }_s^{\prime }$ differed strongly from the real value.

The best fit was reached with a Levenberg-Marquardt algorithm by varying ${\mu }_a$ and ${\mu }_s^{\prime }$ , a multiplicative constant due to the unknown conversion of pixel counts to reflectance values, and an additive constant because of the drifting of the camera background. In the fitting procedure, we assumed a Poisson-like distribution of noise, so the weight of every value was calculated by its inverse intensity value measured by the camera. The optical properties could be found most accurately in this way.

3.1.

Point-Spread Function of the CCD Camera System

The first measurements were performed without using a neutral density filter. The exposure time was set to obtain maximum counts in the central bright spot without saturating the pixel. We used $t=0.01\mathrm{s}$ in this case. The bright spot had digital counts on the order of $3\cdot {10}^4$ . Pixels immediately surrounding the bright spot consisted of intensities of ${10}^3$ digital counts or lower. The challenge was that outside of the few bright pixels, the signal level was low, which led to a low signal-to-noise ratio. It was, therefore, necessary to raise digital counts in that area. Thus, another measurement had to be accomplished with a longer exposure time $(t=1\mathrm{s})$ , which resulted in saturation of the central pixels. This led to an improvement of the signal-to-noise ratio for the pixels far away from the point source but also increased the intensity signal of the whole curve. For further interpretation, we took only the part of the curve with intensity values below saturation into account.

Because of the uncertain effective exposure time, i.e., the unknown opening time of the shutter, it is inevitable to determine the exact factor between the opening times of both measurements. This can be achieved by comparing two nonsaturated measurements with different shutter opening times. As a result, we received the PSF with a high signal-to-noise ratio even for large distances to the illumination point.

We determined the PSF in this way with the point source illuminating several different positions of the CCD chip. Small variations between the single PSFs were found only for large radial distances, i.e., distances between the image of the point source and the other pixels larger than $20\mathrm{mm}$ , resulting from the drifting of the camera’s dark noise. Despite the fact that we took an associating dark image immediately before the real measurement, it is not possible to determine the low intensity values exactly. Figure 3 shows the average of five different PSFs, which have been determined for five different positions of the CCD chip. As can be seen, relatively large error bars appear for distances far away from the illumination point, whereas for small distances, the errors are smaller.

Fig. 3

Mean PSF measured for five different positions on the CCD chip. Error bars show standard deviation. Errors bars are not centered because of logarithmic scaling.

Therefore, careful characterization of the camera’s PSF is essential because the accuracy of the PSF will affect the accuracy of the received optical coefficients. Each optical system will have a characteristic PSF that must be measured precisely. The PSF tends to spread energy everywhere and lowers the contrast. Images and boundaries of bright objects will be seriously affected by the PSF. This can make per-pixel interpretation erroneous because pixels that map places on the phantom surface where no light is remitted can be increased by a bias. The larger the extension of the source is, the worse is the bias. One way to remove the artificial radiance variation is to deconvolute the image with the system’s PSF and return the intensity values back to the true values.

3.2.

Optical Parameters of Phantom Measurements

As described earlier, we assumed that in liquid tissue phantoms, Lipovenös is the only scatterer and trypan blue is the only absorber; however, for very small concentrations of trypan blue, the absorption of water, which is ${\mu }_a\approx 0.0003{\mathrm{mm}}^{-1}$ at $633\mathrm{nm}$ , has to be regarded. We investigated phantoms with optical properties between ${\mu }_a\approx 0.003\dots 0.1{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }\times =0.55\dots 2.0{\mathrm{mm}}^{-1}$ , which are typical for biological tissues. Table 1 shows the results for the fit of the experimental data with the solution of the diffusion equation for an infinite homogenous turbid medium both when using and without using a convolution algorithm.

Table 1

Expected and calculated optical parameters of tissue phantoms without and with a convolution algorithm. The absorption coefficient of water, which is ≈0.0003mm−1 at 633nm , has been considered

In general, we found that the reduced scattering coefficient can be determined more accurately than the absorption coefficient for the range stated earlier. For both small absorption and reduced scattering coefficients, we found good agreement between the expected values and the optical parameters received with both fitting algorithms (rms error in ${\mu }_{a}6\%$ , rms error in ${\mu }_s^{\prime }$ ).

On the other hand, as can be seen in Table 1, for large optical parameters, the fit without the convolution algorithm leads to systematically decreased absorption coefficients and systematically increased reduced scattering coefficients. For ${\mu }_a=0.1003{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=2.00{\mathrm{mm}}^{-1}$ , the errors were 38% in ${\mu }_a$ and 31% in ${\mu }_s^{\prime }$ . Figure 4 shows a comparison of the measured reflectance and the intensities for both fitting algorithms. As can be seen, especially the part for distances larger than $8\mathrm{mm}$ was not fitted well without using the convolution algorithm. The additive constant was overestimated in this case, which led to a constant intensity for distances larger than $11\mathrm{mm}$ . So the bending of the curve cannot be described accurately in that way. However, as the measured intensity decreases slightly for larger distances, the systematically falsified optical parameters are explainable.

Fig. 4

Measured reflectance and fit to measured data with and without using a convolution algorithm for ${\mu }_a\approx 0.10{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }\approx 2.0{\mathrm{mm}}^{-1}$ . The fitting range starts at $0.8\mathrm{mm}$ and includes three orders of magnitude.

Only when considering the specific influence of the objective by convoluting the data with the PSF could a significant improvement be achieved, and systematic discrepancies could be eliminated. The measured data could be fitted quite exactly (see Fig. 4) The convoluted curve now agrees with the experimental data both for small and for larger distances to the central beam. The slight decrease of the intensity for distances larger than $10\mathrm{mm}$ is recorded quite well.

Again four fitting parameters had to be used. Despite the consideration of light that is scattered by the objective of the camera with the help of the convolution algorithm, we found that the use of an additive constant is essential for correct results of the optical parameters, although the additive constant takes over a fairly small value in this case. The dark noise of the CCD camera is not stable temporally despite cooling down the camera to $-35°\mathrm{C}$ . Although before all measurements, a camera background was determined and these background data were subtracted from the immediately following measurements, the drift of the dark noise led to nonsatisfying results. Thus, the additive fitting parameter accounts for the instabilities of the dark noise. In conclusion, the fitting algorithm with both a variable additive and a variable multiplicative constant allows accurate determination of the optical properties for tissue phantoms in the range of typical biological values.

It should be added that if a camera system with a stable background exists, it would be possible to accomplish the fitting algorithm with only three fitting parameters, neglecting the variable additive constant. If so, the convolution of the measured intensity values with the exact PSF would be even more important because the fitting algorithm is less flexible with only three fitting parameters.

With our approach, we found good agreement for the whole range of optical parameters with rms errors in ${\mu }_a$ of 8% and rms errors in ${\mu }_s^{\prime }$ of 4% (see Table 1). Relatively large errors appearing, e.g., for ${\mu }_a=0.0978{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.55{\mathrm{mm}}^{-1}$ are probably due to diffusion theory, which requires ${\mu }_a⪡{\mu }_s^{\prime }$ . This assumption is not fulfilled here.

We also accomplished determination of the optical properties of tissue phantoms with slightly varying intensities of the PSF for larger distances to the illumination point but within the standard deviation (see Fig. 3). The calculated optical parameters differ only marginally with an rms error in ${\mu }_a$ of 11% and an rms error in ${\mu }_s^{\prime }$ of 5%.

To check the consistency of the method, used we accomplished a blind study with 12 tissue phantoms and their optical parameters extending the whole range of Table 1. When running the fitting algorithm a few times (this had to be done because of the completely unknown reduced scattering coefficient, the invalidity of the diffusion equation for distances smaller than $1∕{\mu }_s^{\prime }$ , and the dependence of the beginning of the fitting range on the value chosen for ${\mu }_s^{\prime }$ at the beginning; see Sec. 2.3), the results were similar with errors in ${\mu }_a$ of 10% and errors in ${\mu }_s^{\prime }$ of 5%.