Clean image synthesis and target numerical marching for optical imaging with backscattering light.

Scanning backscattering imaging and independent component analysis (ICA) are used to probe targets hidden in the subsurface of a turbid medium. A new correction procedure is proposed and used to synthesize a "clean" image of a homogeneous host medium numerically from a set of raster-scanned "dirty" backscattering images of the medium with embedded targets. The independent intensity distributions on the surface of the medium corresponding to individual targets are then unmixed using ICA of the difference between the set of dirty images and the clean image. The target positions are localized by a novel analytical method, which marches the target to the surface of the turbid medium until a match with the retrieved independent component is accomplished. The unknown surface property of the turbid medium is automatically accounted for by this method. Employing clean image synthesis and target numerical marching, three-dimensional (3D) localization of objects embedded inside a turbid medium using independent component analysis in a backscattering geometry is demonstrated for the first time, using as an example, imaging a small piece of cancerous prostate tissue embedded in a host consisting of normal prostate tissue.

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(ρ,ρ) where ρ covers the whole 2D array (for example, ρ enumerates all pixels on a CCD image) for a series of beam scanning positions at ρ on the sample surface (z = z = 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 δμ of volume ΔV is located at r′ = (ρ′, z′) with the extension Δz ≪ 1 along the axial direction (δμ is constant within ΔV and 0 outside), denote the exact Green's function for light propagation in the host medium G0 and the intensity of the incident pencil beam I0, 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 ρ to the origin 0) gives

where the error term is given by

with the superposition of the absorber at all possible shifted positions, N is the total number of scanning positions, and A ≡ N2 is the total area of the scanning grid with step size a, provided that the scanning area A is sufficiently large and the pitch a 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 δμΔV are located on z = z planes (j = 1,2,...,N), the averaging over the shifted images yields

where h = Σ h1(ρ,0;z) 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 ρ − ρ and is shared for the difference images ΔI(ρ,ρ) obtained at all scanning positions ρ. 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 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, and ΔI = ΔI(ρ,ρ)−h(ρ−ρ,0), respectively.

2.2. Numerical target marching

After the difference images ΔI(ρ,ρ) 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, G0(r, r) and G0(r, r), 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 to r by a two-step process: forward propagating light first bouncing back by the bottom semi-infinite medium (z > z), and then diffusing from the layer z = z to the top surface z = z < z. This approximation is justified as the probing beam still retains its forward direction at the shallow depth (|z − z| ~ a few l) 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 G0(ρ,0), is exactly the clean image, , 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 G0(ρ,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, G0(r, r), 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 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 is the Fourier transform of the jth independent intensity on the detector plane, G0(q,0) is the Fourier transform of the clean image, , and g(q, z, z) is the Fourier transform of g(ρ, z, z) representing the propagation of a plane wave of spatial modulation frequency q from the z = z plane to the surface z = z.

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 migrating to the z = z plane. In the center-moved diffusion model [14-16], g(ρ, z, z) = exp(−κr)/r where is the attenuation coefficient for incident beam of intensity modulation frequency ω, D0 = l/3 is the diffusion coefficient, l is the transport mean free path, is the absorption coefficient for the host medium, , and z* = z + l is the location of the effective source by displacing the incident point of the collimated beam one l along the incident direction. The Green's function g in the Fourier space is given by g(q, z, z) = 2π exp(−Q|z − z*|)/Q and Eq. (8) is simplified to:

where .

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 G0(r, r) from inverse Fourier transform of Eq. (9) to the retrieved jth independent intensity distribution on the detector plane. G0(r, r) 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 × 38mm2. 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.8mm2. One pixel of the image corresponds to 0.086×0.086mm2 on the sample surface.

Fig. 1.

Experiment setup.

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 and μ′ ≃ 0.53mm−1, 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 μ′ 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 and the correction ratio image 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 represents the propagator, G0(r,0), the intensity distribution of backscattering light from a collimated beam incident at the origin.

Fig. 2.

The clean host image and the correction ratio. The clean image is shown in a 10-base logarithm scale. The correction ratio ( ) is displayed on the right pane.

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.

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).

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 = 0 and r) of diffuse optical imaging, the integral ∫G0(r;ρ′, z′)G0(ρ′, z′;0)d2ρ′ in Eq. (3) equals approximately fG0(r,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

where f 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 G0(r,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 G0(r,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 G0(r,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.

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.