1. Introduction

Biomedical optical imaging is an imaging modality that uses light to probe structural and functional variations in tissue. One prominent example is the detection of a tumor in a human breast [1]. With recent developments, wherein a scanning beam and a CCD camera are employed to replace the illumination and detection fiber-optics of the conventional optical tomography system, the size of generated data sets is orders of magnitude larger than those acquired with a fiber-based system and the solutions demand fast inversion algorithms [2, 3]. Fast inversion algorithms rely on a linearized inversion procedure on the scattering field produced by the inhomogeneities in the turbid medium, i.e., the difference between the “dirty” image of the medium with embedded targets and the “clean” one due to the (homogeneous) host medium alone. As an alternative approach, Optical Imaging using Independent Component Analysis (OPTICA) first unmixes the signal from targets embedded in a turbid medium buried in the scattering field using independent component analysis (ICA). It then detects, locates, and characterizes targets based on the analysis of unmixed independent components [4–6]. OPTICA has been demonstrated to detect and locate small absorptive [6, 7], scattering [8], or fluorescent targets [5] in tissue phantoms and model breast samples in a transmission geometry.

The reflection (backscattering) geometry which is most amenable for probing targets hidden in the subsurface of a turbid medium poses more difficulties than the transmission geometry. Several backscattering geometry related issues are: (1) The estimation of the clean image of the homogeneous host medium from the measured data set obtained at an array of scanning positions can not be obtained simply by averaging the measured data set. The subsurface inhomogeneities affect the backscattering images appreciably. (2) The commonly used diffusion model in optical imaging breaks down, in general, for probing the subsurface and computational expensive methods have to be used to model backscattering light at a close source-detector separation [9]. (3) The unknown surface properties of the biological sample such as superficial structure and roughness impact the spatial distribution of the backscattering light. These issues are critical in analyzing images measured in a backscattering geometry and are addressed in this article. We first propose a new procedure to synthesize a “clean” image of a homogeneous host medium numerically from a set of raster-scanned “dirty” images of a medium with embedded targets recorded by a large frame CCD camera. The independent intensity distribution on the surface of the medium corresponding to individual targets is unmixed using ICA of the difference between the set of dirty images and the clean image. A novel analytical method, which marches the target to the surface until a match with the retrieved independent component is accomplished, is then introduced to locate the targets, automatically accounting for the unknown surface properties of the sample. With clean image synthesis and target numerical marching, three-dimensional localization of objects in biological tissue using ICA in a reflection geometry is demonstrated for the first time by imaging a small piece of cancerous prostate tissue embedded in a host consisting of normal prostate tissue.

2. Theory

2.1. Clean image synthesis

Fast inversion algorithms require as input the scattering field due to the sought-after targets. This scattering field is the difference between the measured intensity distribution of backscattering light and a clean image of the homogeneous host medium without any embedded targets. The image for a target-free host medium is not available in most applications and must be generated from the measured data set of the medium with targets embedded inside. Denote the recorded images I(ρ_d,ρ_s) where ρ_d covers the whole 2D array (for example, ρ_d enumerates all pixels on a CCD image) for a series of beam scanning positions at ρ_s on the sample surface (z_d = z_s = 0). One crude method to obtain the image of the host medium is simply the average of all the array images after shifting the scanning position to the origin. This results in, however, only a distorted version of the real background image. To illustrate this point, assume one thin slice absorptive object δμ_a of volume ΔV is located at r′ = (ρ′, z′) with the extension Δz ≪ 1 along the axial direction (δμ_a is constant within ΔV and 0 outside), denote the exact Green's function for light propagation in the host medium G ₀ and the intensity of the incident pencil beam I ₀, the measured array image is given by

We have set the speed of light to be unity for clarity. Scattering targets can be treated in a similar fashion.

The average of shifted array images (shifting the scanning position for each image ρ_s to the origin 0) gives

where the error term is given by

with $u(\rho \prime ,z\prime )={\Sigma }_{{\rho }_s}\delta {\mu }_a(\rho \prime -{\rho }_s,z\prime )$ the superposition of the absorber at all possible shifted positions, N_s is the total number of scanning positions, and A_s ≡ N_sa ² _s is the total area of the scanning grid with step size a_s, provided that the scanning area A_s is sufficiently large and the pitch a_s is fine enough such that the scanning area covers the whole range of the Green's function above the noise floor on the z = z′ plane and the value of the Green's function does not change appreciably within one scanning grid.

When a total of N slices of strength δμ_jΔV_j are located on z = z_j planes (j = 1,2,...,N), the averaging over the shifted images yields

where h = Σ_j h ₁(ρ_d,0;z_j) and the difference image between the dirty images and Ī can be written as:

Here the error term h in Eqs. (4) and (5) only depends on ρ_d − ρ_s and is shared for the difference images ΔI(ρ_d,ρ_s) obtained at all scanning positions ρ_s. Hence a simple estimation of h is provided by averaging all shifted difference images which are minimally perturbed by the embedded objects, such as,

where B denotes the perimeter of the scanning grid which contains a total of N_B scanning positions. After h has been obtained, the clean image of the host medium and the difference between the dirty images and the clean one are given by, $\overline{I^c}=\overline{I}({\rho }_d,0)+h({\rho }_d,0)$ and ΔI^c = ΔI(ρ_d,ρ_s)−h(ρ_d−ρ_s,0), respectively.

2.2. Numerical target marching

After the difference images ΔI^c(ρ_d,ρ_s) have been obtained, ICA can be used to unmix the signal arising from individual targets and obtain the independent intensity distribution due to the jth target (j = 1,2,...,N). The jth independent component constitutes the projection of the Green's functions, G ₀(r_s, r_j) and G ₀(r_d, r_j), on the source and detector plane, respectively [4]. In the reflection geometry, the intensity of the backscattering light under illumination of a pencil incident beam is strongly peaked at the point of incidence. There is no analytical solution for the Green's function in this case and numerical simulations such as Monte Carlo techniques have been commonly used instead to generate the Green's function [9].

To overcome these difficulties, we approximate light propagation from r_j to r_d by a two-step process: forward propagating light first bouncing back by the bottom semi-infinite medium (z > z_j), and then diffusing from the layer z = z_j to the top surface z = z_d < z_j. This approximation is justified as the probing beam still retains its forward direction at the shallow depth (|z_j − z_s| ~ a few l_t) whereas the propagation direction of the bounced back light is randomized within the backward hemisphere, and light backscattered without interaction first with the bottom semi-infinite medium is negligible. The main advantage of this approximation is that the light propagator for the first process of light backscattering from the host medium, i.e., the Green's function G ₀(ρ_d,0), is exactly the clean image, $\overline{I^c}$ , which has been computed directly from the measured data set as outlined above without assumption of any particular light propagation model in the medium. Another advantage of using the clean image from experimentally measured data set as the Green's function G ₀(ρ_d,0) is that this light propagator automatically incorporates the blurring effect due to the unknown superficial structure and surface roughness of a biological sample.

The projection of the Green's functions for the jth target, G ₀(r_d, r_j), on the detector plane can now be rewritten as:

using the first Rayleigh-Sommerfeld integral [10, 11] where the integration is performed over the z = z_j plane and g is the Green's function inside an infinite homogeneous medium for a diffuse photon density wave [12]. Performing Fourier transform over the lateral coordinates, a simple relationship is obtained:

Here $G_{0_j}\left(q\right)$ is the Fourier transform of the jth independent intensity on the detector plane, G ₀(q,0) is the Fourier transform of the clean image, $\overline{I^c}$ , and g(q, z_d, z_j) is the Fourier transform of g(ρ, z_d, z_j) representing the propagation of a plane wave of spatial modulation frequency q from the z = z_j plane to the surface z = z_d.

To incorporate the detection condition where only photons escaping the medium in the normal direction is being detected, we first use the optical reciprocal property [13] and consider instead for a normally incident beam at r_d migrating to the z = z_j plane. In the center-moved diffusion model [14–16], g(ρ, z_d, z_j) = exp(−κr)/r where $\kappa =\sqrt{\frac{({\mu }_{a_0}-i\omega )}{D_0}}$ is the attenuation coefficient for incident beam of intensity modulation frequency ω, D0 = l_t/3 is the diffusion coefficient, l_t is the transport mean free path, ${\mu }_{a_0}$ is the absorption coefficient for the host medium, $r=\sqrt{{\rho }^2+{(z_j-z_d^{*})}^2}$ , and z ^* _d = z_d + l_t is the location of the effective source by displacing the incident point of the collimated beam one l_t along the incident direction. The Green's function g in the Fourier space is given by g(q, z_d, z_j) = 2π exp(−Q|z_j − z ^* _d|)/Q and Eq. (8) is simplified to:

where $Q\equiv \sqrt{q^2+{\kappa }^2}$ .

Equation (9) is the main expression for numerical target marching. The deeper the target is hidden inside the medium, the stronger damping of the high spatial frequency details in the independent intensity distribution corresponding to that target. The degree to which the independent intensity distribution flattens and loses its high spatial frequency details comparing to the clean image provides the basis to obtain the depth of the target.

The position of the jth target is obtained by fitting G ₀(r_d, r_j) from inverse Fourier transform of Eq. (9) to the retrieved jth independent intensity distribution on the detector plane. G ₀(r_d, r_j) needs to be deconvoluted with the incident beam profile before fitting to remove the effect of the beam profile on the clean image.

3. Experiment

Figure 1 displays the schematic diagram of the experimental setup. The sample is illuminated by a collimated diode laser (635 nm) in the direction close to the normal to the surface. Two galvanometric mirrors and controllers (General Scanning Inc. Lumonics-GSIL, Bedford, M.A., Model DSC W/HCI) are used to scan the illuminating beam along the x- and y-directions on the front surface of the sample. The two dimensional image of the backscattering light from the sample in the normal direction is recorded by a CCD camera (Photometrix CH350L, 1024 × 1024 pixels, 16 bit) after passing through a neutral density filter (NF) and lens (L). A custom software (GSIL, WinMCL Plus) is used to control the scanning speed, the two dimensional scanning grid, and to trigger the recording of the backscattering image at each scanning position. A CCD control board (Photometrix PCI-X 01-490-400 with PVCAM driver) and the corresponding imaging software (Digital Optics V++) are used to record the 2D backscattering images. The spot size of the scanning beam is ~ 500µm. The sample studied consists of a large piece of normal prostate tissue with a small piece of cancerous prostate tissue (4mm × 4mm × 1.5mm) embedded inside at a depth of z = 3.0 mm from the front surface. The thickness of the whole sample is 10 mm and with lateral dimensions of 45 × 38mm². The embedded and the host prostate tissue were verified to be cancerous and normal, respectively, by pathology. The scanning grid is 8 × 8 with a step size of 2.13 mm. The scanning window could be adjusted and/or enlarged to contain the target if the target is found to be not well contained inside the scanning window. The CCD camera captures the backscattering light within a front surface window of size 87.8 × 87.8mm². One pixel of the image corresponds to 0.086×0.086mm² on the sample surface.

The optical property of prostate tissue reported in the literature varies over a wide range [17–19]. The absorption and reduced scattering coefficients for normal prostate tissue (the host medium) in this study were ${\mu }_{a_0}\simeq 0.026{\mathrm{mm}}^{-1}$ and μ′_s ≃ 0.53mm⁻¹, respectively, from fitting the diffusion model to the clean image for the source-detector separation larger than 9mm. The cancerous prostate tissue has a much smaller μa (about 10 times smaller) than the normal prostate tissue whereas the cancerous prostate tissue has a slightly larger μ′_s than the normal prostate tissue at the probing wavelength 635 nm [20]. The piece of cancerous prostate tissue embedded inside the model prostate tissue sample behaves predominantly as an absorption inhomogeneity.

4. Results

The clean host image $\overline{I^c}$ and the correction ratio image $\frac{h}{\overline{I^c}}$ is shown in Fig. 2. The correction is computed using a total of 30 images measured on the perimeter of the scanning grid. The clean image is shown in a 10-base logarithm scale. The correction ratio is displayed on the right pane. The clean image $\overline{I^c}$ represents the propagator, G ₀(r_d,0), the intensity distribution of backscattering light from a collimated beam incident at the origin.

Figure 3 displays the line profile of the clean image, the independent component originating from the cancerous prostate target, and the fitting of the Green's function to the independent component along the vertical direction. The full width at half maximum (FWHM) is 0.56mm for the clean host image and FWHM expands to 3.98mm for the independent intensity distribution on the detector plane arising from the target. The depth of the target is obtained by numerically marching the host image through the medium via Eq. (9) after deconvolution with the beam profile in the Fourier space until it matches the independent intensity on the detector plane. The depth from data fit is 3.01 mm, in excellent agreement with the known value. For comparison, the same analysis without applying the correction procedure and using the uncorrected difference images (ΔI) yields the depth for the target at 2.75 mm, performing much worse.

5. Discussion and conclusion

The proposed procedure to obtain a clean image of a target-free host medium from raster-scanned dirty images of the medium with embedded targets is applicable in both transmission and reflection geometries. The magnitude of the correction h directly links to the strength of the targets embedded in the host medium. For example, in the transmission geometry (the object sits on the plane z = z′ between r_s = 0 and r_d) of diffuse optical imaging, the integral ∫G ₀(r_d;ρ′, z′)G ₀(ρ′, z′;0)d ² ρ′ in Eq. (3) equals approximately fG ₀(r_d,0) where the scaling factor f depends on the optical property and the geometry of the host medium and the depth of the object z′. The value of f is well approximated by the corresponding ratio of the Green's function in the Fourier space with the lateral spatial frequency set to 0. The Green's function in the Fourier space under the diffusion approximation are well known [21]. The total strength of embedded objects in the transmission geometry is proportional to the error term via

 

Fig. 1. Experiment setup.

Download Full Size | PDF

 

Fig. 2. The clean host image and the correction ratio. The clean image is shown in a 10-base logarithm scale. The correction ratio ( $\left(\frac{h}{\overline{I^c}}\right)$ ) is displayed on the right pane.

Download Full Size | PDF

where f_j is the scaling factor for the jth target.

One major challenge of applying backscattering light to probe the subsurface (at the depth up to a few transport mean free paths) of a turbid medium is that the diffusion approximation to radiative transfer is inadequate in this case and computationally much expensive methods need to be used [3, 9]. The novelty of numerical target marching is that it alleviates this difficulty using the Green's function G ₀(r_d,0) of light reflection from a semi-infinite medium, i.e., the clean image of the host medium, obtained directly from experimentally measured array images. As the result, no light propagation model is assumed in obtaining this Green's function G ₀(r_d,0), and the blurring effect due to the unknown superficial structure and surface roughness of a biological sample on light propagation is automatically incorporated in G ₀(r_d,0). This new approach, which builds the proper Green's function for optical imaging from experimentally measured clean image, should be applicable in modeling light propagation in general when the commonly used diffusion model is inadequate or incorrect.

 

Fig. 3. The line profile of the clean host image (Left), the independent component originating from the cancerous prostate target (Middle), and the fitting of the Green's function to the independent component along the vertical direction (Right).

Download Full Size | PDF

The geometry of choice to detect prostate cancer is to use backscattering light through rectum. As the rectal wall is relatively thin (~2.5–3.0mm) [22], backscattering light can penetrate the rectal wall and examine the peripheral zone where most prostate cancer initiates. The developed methods are, in particular, useful for optical imaging of prostate cancer through rectum in the backscattering geometry. We have demonstrated in our prior work that optical imaging using independent component analysis to be an effective method to locate and characterize absorption, scattering, or fluorescent targets as small as several millimeters in the transmission geometry [5–8]. As the FWHM of the Green's function for light backscattering from subsurface in the reflection geometry is much smaller than that of the Green's function for light transmission through a slab in the transmission geometry, we expect even smaller size of inhomogeneity can be detected in the backscattering geometry provided a sufficient signal to noise ratio.

In conclusion, clean image synthesis of a host medium and target numerical marching have been presented for optical imaging. These techniques are, in particular, useful for imaging in a reflection (backscattering) geometry. Employing clean image synthesis and target numerical marching, three-dimensional localization of objects embedded inside a turbid medium using independent component analysis in a backscattering geometry has been demonstrated for the first time with an example imaging a small piece of cancerous prostate tissue embedded in a host consisting of normal prostate tissue.