Open Access
1 September 2005 Optical imaging of turbid media using independent component analysis: theory and simulation
Author Affiliations +
Abstract
A new imaging approach for 3-D localization and characterization of objects in a turbid medium using independent component analysis (ICA) from information theory is developed and demonstrated using simulated data. This approach uses a multisource and multidetector signal acquisition scheme. ICA of the perturbations in the spatial intensity distribution measured on the medium boundary sorts out the embedded objects. The locations and optical characteristics of the embedded objects are obtained from a Green's function analysis based on any appropriate model for light propagation in the background medium. This approach is shown to locate and characterize absorptive and scattering inhomogeneities within highly scattering medium to a high degree of accuracy. In particular, we show this approach can discriminate between absorptive and scattering inhomogeneities, and can locate and characterize complex inhomogeneities, which are both absorptive and scattering. The influence of noise and uncertainty in background absorption or scattering on the performance of this approach is investigated.

1.

Introduction

Optical probing of the interior of multiply scattering colloidal suspensions and biological materials has attracted much attention over the last decade. In particular, biomedical optical tomography and spectroscopy, which has the potential to provide functional information about brain activities and diagnostic information about tumors in breast and prostate, are being actively pursued.1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17 Simultaneous developments in experimental apparatus and techniques for object interrogation and signal acquisition,2, 4, 5, 18, 19 analytical models for light propagation,10, 20, 21, 22 and computer algorithms for image reconstruction6, 8 hold promise for realization of these potentials of optical tomography.

Multiple scattering of light in thick turbid media precludes direct imaging of embedded targets. One typically uses an inverse image reconstruction6, 8 (IIR) approach to reconstruct a map of the optical properties, such as absorption coefficient (μa) and scattering coefficient (μs) , of the medium by matching the measured light intensity distribution on the boundary of the turbid medium to that calculated by a forward model for the propagation of light in that medium. The commonly used forward models include the radiative transfer equation17, 23 (RTE), the diffusion approximation (DA) of the RTE (Refs. 6, 8), and random walk of photons.24, 25

The inversion problem is ill-posed and must be regularized to stabilize the inversion at a cost of reduced resolution.8, 26 Both iterative reconstruction and noniterative linearized inversion approaches have been used to solve the inversion problem in optical tomography, which is weakly nonlinear, with limited success. Reconstruction of images with adequate spatial resolution and accurate localization and characterization of the inhomogeneities remain a formidable task. The time required for data acquisition and image reconstruction is another important consideration.

In this paper, we present a novel algorithm based on the independent component analysis27, 28 (ICA) from information theory to locate absorptive and scattering inhomogeneities embedded in a thick turbid medium and demonstrate the efficacy using simulated data. ICA has been successfully applied in various other applications such as electroencephalogram and nuclear magnetic resonance spectroscopy.28, 29, 30, 31 We refer to this information theory-inspired approach as optical imaging using independent component analysis, abbreviated as, OPTICA. The novelty of OPTICA over other ICA applications is that OPTICA associates directly the independent components to the Green’s functions responsible for light propagation in the turbid medium from the inhomogeneities to the source and the detector, and therefore the retrieved independent components can be used to locate and characterize the inhomogeneities.

OPTICA uses a multisource illumination and multidetector signal acquisition scheme providing a variety of spatial and angular views essential for three-dimensional (3-D) object localization. Each object (or inhomogeneity) within the turbid medium alters the propagation of light through the medium. The influence of an object on the spatial distribution of the light intensity at the detector plane involves propagation of light from the source to the object, and from the object to the detector, and can be described in terms of two Green’sfunctions (propagators) describing light propagation fromsource to the object and that from the object to the detector, respectively.

The absorptive or scattering inhomogeneities illuminated by the incident wave are assumed to be virtual sources, and the perturbation of the spatial distribution of the light intensity on the medium boundary is taken to be some weighted mixture of signals arriving from these virtual sources. These virtual sources are statistically independent and can be recovered by ICA of the recorded data set. The number of leading independent components is same as the number of embedded objects. The location and characterization of inhomogeneities are obtained from the analysis of the retrieved virtual sources using an appropriate model of light propagation in the background medium.

The remainder of this paper is organized as follows. Section 2 presents a brief introduction to ICA and reviews the general theoretical framework for OPTICA. Section 3 presents the results from simulations for different configurations. Implications of these results are discussed in Sec. 4.

2.

Theoretical Formalism

2.1.

ICA

Blind source separation is a class of problem of general interest that consists of recovering unobserved signals or virtual sources from several observed mixtures. Typically the observations are the output of a set of sensors, where each sensor receives a different combination of the source signals. Prior knowledge about the mixture in such problems is usually not available. The lack of prior knowledge is compensated by a statistically strong but often physically plausible assumption of independence between the source signals. Over the last decade, ICA has been proposed as a solution to the blind source separation problem and has emerged as a new paradigm in signal processing and data analysis.27, 28, 31, 32

The simplest ICA model, an instantaneous linear mixture model,32 assumes the existence of n independent signals si(t) (i=1,2,,n) and the observation of at least as many mixtures xi(t) (i=1,2,,m) by mn sensors, these mixtures being linear and instantaneous, i.e., xi(t)=j=1naijsj(t) for each i at a sequence of time t . In a matrix notation,

Eq. 1

x(t)=As(t),ARm×n,
where A is the mixing matrix. The j ’th column of A gives the mixing vector for the j ’th virtual source. ICA can be formulated as the computation of an n×m separating matrix B whose output

Eq. 2

y(t)=Bx(t)=Cs(t),BRn×m,CBARn×n,
is an estimate of the vector s(t) of the source signals.

The basic principle of ICA can be understood in the following way. The central limit theorem in probability theory tells us that the distribution of independent random variables tends toward a Gaussian distribution under certain conditions. Thus, a sum of multiple independent random variables usually has a distribution that is closer to Gaussian than any of the original random variables. In Eq. 2 yi(t)=jCijsj(t) , as a summation of independent random variables sj(t) , is usually more Gaussian than sj(t) , while yi(t) becomes least Gaussian when it in fact equals one of the sj(t) . This heuristic argument shows that ICA can be intuitively regarded as a statistical approach to find the separating matrix B such that yi(t) is least Gaussian. This can be achieved by maximizing some measure of non-Gaussianity, such as maximizing kurtosis32, 33 (the fourth-order cumulate), of yi(t) .

ICA separates independent sources from linear instantaneous or convolutive mixtures of independent signals without relying on any specific knowledge of the sources except that they are independent. The sources are recovered by a maximization of a measure of independence (or, a minimization of a measure of dependence), such as non-Gaussianity and mutual information between the reconstructed sources.31, 32 The recovered virtual sources and mixing vectors from ICA are unique up to permutation and scaling.31, 32

2.2.

Optical Imaging Using ICA

The classical approach to propagation of multiply scattered light in turbid media, which assumes that phases are uncorrelated on scales larger than the scattering mean free path ls , leads to the RTE in which any interference effects are neglected.34 The RTE does not admit closed-form analytical solutions in bounded regions and its numerical solution is computationally expensive. The commonly used forward models in optical imaging of highly scattering media is6, 8 the DA to RTE.

The approach OPTICA can be applied to different models of light propagation in turbid media, such as the diffusion approximation,6, 8 the cumulant approximation,20, 22, 35 the random walk model,10, 24 and radiative transfer17, 34 when they are linearized. The diffusion approximation is valid when the inhomogeneities are deep within a highly scattering medium. We discuss only the formalism of OPTICA in the diffusion approximation in this paper.

In this diffusion approximation, the perturbation of the detected light intensities on the boundaries of the medium, the scattered wave field, due to absorptive and scattering objects (inhomogeneities) can be written as3, 13

Eq. 3

ϕsca(rd,rs)=G(rd,r)δμa(r)cG(r,rs)drδD(r)crG(rd,r)rG(r,rs)dr,
to the first order of the Born approximation36 when illuminated by a point source of unit power, where rs , r , and rd are the positions of the source, the inhomogeneity, and the detector, respectively; δμa=μa,objμa and δD=DobjD are the differences in absorption coefficient and diffusion coefficient, respectively, between the inhomogeneity and the background; c is the speed of light in the medium; and G(r,r) is the Green’s function describing light propagation from r to r inside the background turbid medium of absorption and diffusion coefficients μa and D .

Equation 3 is written in the frequency domain and does not explicitly show the modulation frequency ω of the incident wave for clarity. The following formalism applies to continuous wave, frequency-domain, and time-resolved measurements. The time-domain measurement is first Fourier transformed over time to obtain data over many different frequencies.

The Green’s function G for a slab geometry in DA is given by

4.

G(r,r)G(ρ,z,z)=14πDk=[exp(κrk+)rk+exp(κrk)rk]
rk±=[ρ2+(zz±2kd)]12,
for an incident amplitude-modulated wave of modulation frequency ω , where k=0,±1,±2,,ρ=[(xx)2+(yy)2]12 is the distance between the two points r=(x,y,z) and r=(x,y,z) projected onto the xy plane; κ=[(μaiωc)D]12 is chosen to have a nonnegative real part; and the extrapolated boundaries of the slab are located at z=0 and z=d=Lz+2ze , respectively, where Lz is the physical thickness of the slab and the extrapolation length ze should be determined from the boundary condition of the slab.37, 38, 39 Greens’ functions in Eq. 3 for other geometries can be obtained either analytically or numerically.40, 41

2.2.1.

Absorptive inhomogeneities

We first consider absorptive inhomogeneities. Under the assumption that absorptive inhomogeneities are localized, that is, the j ’th one is contained in volume Vj centered at rj(1jJ) , the scattered wave field in Eq. 3 can be rewritten as

Eq. 5

ϕsca(rd,rs)=j=1JG(rd,rj)qjG(rj,rs),
where qj=δμa(rj)cVj is the absorption strength of the j ’th inhomogeneity. The scattered wave is in a form of an instantaneous linear mixture of Eq. 1. One absorptive inhomogeneity is represented by one virtual source qjG(rj,rs) with a mixing vector G(rd,rj) .

As the virtual source qjG(rj,rs) at the j ’th inhomogeneity is independent of the virtual sources at other locations, ICA can be used with the observations obtained for the light source at nJ different positions to separate out both virtual sources sj(rs) and the mixing vectors27, 32 aj(rd) . The j ’th virtual source sj(rs) and the j ’th mixing vector aj(rd) provide the scaled projections of the Green’s function on the source and detector planes, G(rj,rs) and G(rd,rj) , respectively. We can write

6.

sj(rs)=αjG(rj,rs),
aj(rd)=βjG(rd,rj),
where αj and βj are scaling constants for the j ’th inhomogeneity.

Both the location and strength of the j ’th object can be computed by a simple fitting procedure using Eq. 6. We adopted a least-squares fitting procedure given by:

Eq. 7

minrj,αj,βj{rs[αj1sj(rs)G(rj,rs)]2 +rd[βj1aj(rd)G(rd,rj)]2}.
The fitting yields the location rj and the two scaling constants αj and βj for the j ’th inhomogeneity whose absorption strength is then given by qj=αjβj .

2.2.2.

Scattering inhomogeneities

For scattering inhomogeneities, a similar analysis shows the scattered wave can be written as

Eq. 8

ϕsca(rd,rs)=j=1Jgz(rj,rd)qjgz(rj,rs)+j=1Jρdjcosθdg(rj,rd)qjρsjcosθsg(rj,rs)+j=1Jρdjsinθdg(rj,rd)qjρsjsinθsg(rj,rs),
where qj=δD(rj)cVj is the diffusion strength of the j ’th scattering inhomogeneity of volume Vj(j=1,2,,J) ; ρdj=[(xdxj)2+(ydyj)2]12 , ρsj=[(xsxj)2+(ysyj)2]12 ; θd and θs are the azimuthal angles of rdrj and rsrj , respectively, and the two auxiliary functions are given by

Eq. 9

g(r,r)=14πDk=+[(κrk++1)exp(κrk+)(rk+)3(κrk+1)exp(κrk)(rk)3],
and

Eq. 10

gz(r,r)=14πDk=+{(zz+2kd)(κrk++1)exp(κrk+)(rk+)3(z+z2kd)(κrk+1)exp(κrk)(rk)3}.

The scattered wave from one scattering inhomogeneity is thus a mixture of contributions from (3J) virtual sources:

Eq. 11

qjgz(rj,rs),qjρsjcosθsg(rj,rs),qjρsjsinθsg(rj,rs),
with mixing vectors

Eq. 12

gz(rj,rd),ρdjcosθdg(rj,rd),ρdjsinθdg(rj,rd),
where 1jJ , respectively. Both the location and strength of the j ’th scattering object are computed by fitting the retrieved virtual sources and mixing vectors to Eqs. 11, 12 using a least-squares procedure, respectively.

There are, in general, three virtual sources of specific patterns (one centrosymmetric and two dumbbell-shaped) associated with one scattering inhomogeneity, whereas only one centrosymmetric virtual source is associated with one absorptive inhomogeneity. This difference may serve as the basis to discriminate absorptive and scattering inhomogeneities.

The only assumption made in OPTICA is that virtual sources are mutually independent. The number of inhomogeneities within the medium is determined by the number of the independent components presented in the multisource multidetector data set. No specific light propagation model is assumed in this step. The analysis of retrieved independent components from ICA then localizes and characterizes the absorptive and scattering inhomogeneities inside the turbid medium using an appropriate model of the light propagator. Extra independent components may appear, depending on the level of noise in the data. These components can be discarded and only the leading independent components must be analyzed to detect and characterize inhomogeneities of interest.

3.

Results

Simulations were performed for a 50-mm -thick slab, as shown schematically in Fig. 1 . The absorption and diffusion coefficients of the uniform slab are μa=1300mm1 and D=13mm , respectively, close to that of human breast tissue.42 The incident cw beam scans a set of 21×21 grid points covering an area of 90×90mm2 . The spacing between two consecutive grid points is 4.5mm . This light intensity on the other side of the slab is recorded by a CCD camera on 42×42 grid points covering the same area.

Fig. 1

Light intensity on one side of the slab is measured when a point source scans on the other side. Two inhomogeneities are placed at (50, 60, 20) and (30,30,30)mm inside the slab.

051705_1_028505jbo1.jpg

In the simulations presented in the following subsections, we fix the ratio of strength of absorption to that of diffusion to be 0.01, which produce perturbations of comparable magnitude on the light intensities measured on the detector plane from the absorption and scattering inhomogeneities. As the scattered wave is linear with respect to the absorption and diffusion strengths, we also set the strength of either absorption or diffusion to be unity in simulations for convenience.

3.1.

Absorptive Inhomogeneities

Two absorptive inhomogeneities, each of a unity absorption strength, are placed at positions (50, 60, 20) and (30,30,30)mm , respectively. Gaussian noise of 5% was added to the simulated light intensity change on the detector plane. OPTICA operates on a noisy scattering wave ϕsca(rd,rs)[1+n(rd,rs)] , where n(rd,rs) is a Gaussian random variable of a standard deviation 0.05.

ICA of the perturbations in the spatial intensity distributions provided corresponding independent intensity distributions on the source and detector planes. ICA-generated independent intensity distributions on the source and detector planes are shown in Fig. 2 for the two absorptive inhomogeneities. Locations of the absorptive objects are obtained from fitting these independent component intensity distributions to those of the diffusion approximation in a slab Eq. 4 by the least-squares procedure of Eq. 7. The first object is found at (50.0,60.0,20.0)mm and the second one at (30.0,30.0,30.1)mm . The coordinates of both objects agree to within 0.1mm of their known locations. The strengths of the two objects are q1=1.00 and q2=0.99 , respectively, with an error not greater than 1% of the true values.

Fig. 2

Normalized independent spatial intensity distributions on the input (or source) plane (the first row), the exit (or detector) plane (the second row), and the least-squares fitting using Eq. 7 (the third row). The left column is for the first absorptive inhomogeneity at (50,60,20)mm and the right column is for the second absorptive inhomogeneity at (30,30,30)mm . On the third row, the horizontal profile of intensity distributions on the source plane (diamonds) and on the detector plane (circles) are displayed, and solid lines show the respective Green’s function fit used for obtaining locations and strengths of objects. The noise level is 5%.

051705_1_028505jbo2.jpg

3.2.

Discrimination between Absorptive and Scattering Inhomogeneities

In the second example, one absorptive object of absorption strength of 0.01 is placed at (50,60,20)mm and one scattering object of diffusion strength of negative unity (corresponding to an increase in scattering for the inhomogeneity) is placed at (30,30,30)mm , respectively. We added 5% Gaussian noise to the simulated light intensity change on the detector plane.

Figure 3 shows the ICA-generated independent intensity distributions on the source and detector planes and the least-squares fitting. The first column corresponds to the absorptive inhomogeneity. The second through fourth columns correspond to the scattering object, which produces one pair of centrosymmetric and two pairs of dumbbell-shaped virtual sources and mixing vectors. The absorptive inhomogeneity is found to be at (50.2,60.3,20.2)mm with a strength q1=0.0101 . The scattering object produces three pairs (one centrosymmetric and two dumbbell-shaped) of virtual sources and mixing vectors centering around the position (x,y)=(30,30)mm (see the second through fourth columns in Fig. 3). The dumbbell-shaped virtual source or mixing vector comprises one bright part and its antisymmetric dark counterpart. The resolved position and strength of the scattering object are found to be (30.0,30.0,30.0)mm and q2=0.99 , (32.1,32.4,30.2)mm and q2=0.96 , and (31.3,30.2,27.1)mm and q2=1.05 , respectively, through fitting to the individual pair. For the scattering object, the best result is obtained from the fitting to the first pair of centrosymmetric virtual source and mixing vector from the scattering object. Taking the position and strength of the scattering object to be that from fitting the centrosymmetric virtual source and mixing vector, the error of the resolved positions of both objects is within 0.3mm of their known locations. The error of the resolved strengths of both objects is approximately 1% of the true values.

Fig. 3

Normalized independent spatial intensity distribution on the source plane, i.e., virtual sources (the first row) and on the detector plane, i.e., mixing vectors (the second row), and the least-squares fitting (the third row). The first column is for the first absorptive inhomogeneity at (50,60,20)mm and the second through fourth columns are for the second scattering inhomogeneity at (30,30,30)mm and represent the centrosymmetric and two dumbbell-shaped pairs of virtual sources and mixing vectors, respectively. The dumbbell comprises one bright part and its antisymmetric dark counterpart. In the third row, the profile of intensity distributions on the source plane (diamonds) and on the detector plane (circles) are displayed, and solid lines show the respective Green’s function fit used for obtaining locations and strengths of objects. The X coordinate is the horizontal coordinate. The ρ coordinate is the coordinate on the symmetrical line passing through the dumbbell axis. The small dark circular region appearing near the right-upper corner of the normalized independent spatial intensity distribution on the first row and in the fourth column is an artifact.

051705_1_028505jbo3.jpg

3.3.

Colocated Absorptive and Scattering Inhomogeneities

For one complex inhomogeneity that is both absorptive and scattering, the two pairs of dumbbell-shaped virtual sources and mixing vectors produced by its scattering abnormality can be used to obtain its scattering strength. By subtracting the scattering contribution off the measured scattered wave, our procedure can be applied again to the cleaned data and we proceed to obtain its absorption strength. The third example considers a complex inhomogeneity at (30,30,20)mm with strengths of absorption q1=0.01 and diffusion q2=1 (corresponding to a decrease in scattering), respectively. We added 5% Gaussian noise to the simulated light intensity change on the detector plane.

Figure 4 shows the ICA-generated independent intensity distributions on the source and detector planes and the least-squares fitting. The first and second columns correspond to the pairs of dumbbell-shaped virtual sources and mixing vectors produced by its scattering component. The position and strength of this diffusive component is obtained to be (32.7,33.0,20.5)mm and q2=0.95 , and (31.7,30.1,20.4)mm and q2=0.96 by fitting the two individual dumbbell-shaped pairs, respectively. The position and strength of the diffusive component is found to be (30.9,30.9,20.4)mm and q2=0.95 if both dumbbell-shaped virtual sources and mixing vectors are used in fitting. The third column of Fig. 4 corresponds to its absorptive component obtained by first removing the scattering contribution from the measured scattered wave. The depth and strength of the absorption component is found to be (30.8,30.7,32.7)mm and q1=0.0091 .

Fig. 4

Normalized independent spatial intensity distribution on the source plane (the first row) and on the detector plane (the second row), and the least-squares fitting (the third row) for one inhomogeneity located at (30,30,20)mm with strengths of absorption q1=0.01 and diffusion q2=1 . The first and second columns correspond to the pairs of dumbbell-shaped virtual sources and mixing vectors produced by its scattering component. The third column corresponds to its absorptive component obtained by first removing the scattering contribution. In the third row, the profile of intensity distributions on the source plane (diamonds) and on the detector plane (circles) are displayed, and solid lines show the respective Green’s function fit used to obtain locations and strengths of objects. The X coordinate is the horizontal coordinate and the ρ coordinate is the coordinate on the symmetrical line passing through the dumbbell axis.

051705_1_028505jbo4.jpg

The error in positioning the scattering component is less than 1mm and the error of the resolved strength of the scattering strength is 5% . The errors in positioning and the resolved strength of the absorptive component equal 3mm and 10% , respectively, which are larger than those of the scattering component because the error is amplified when the estimated scattering component is used in subtraction off its contribution to the scattered wave in our procedure.

3.4.

Effect of Noise

To demonstrate the effect of noise on the performance of OPTICA, different levels of Gaussian noise were added to the simulated light intensity change on the detector plane.

Figure 5 shows the case presented in Fig. 3 of Sec. 3(B) with, instead, 10 and 20% Gaussian noise added to the scattered wave. The resolved absorptive inhomogeneity is at (50.2,60.3,20.1)mm with strength 0.0101 at 10% noise, and at (50.1,60.3,20.5)mm with strength 0.0102 at 20% noise. The resolved position and strength of the scattering object are found to be (30.0,30.1,30.0)mm and q2=0.98 , (32.1,32.4,30.4)mm and q2=0.95 , and (31.4,30.1,27.5)mm and q2=1.00 , respectively, through fitting to the pair of centrosymmetric and two pairs of dumbbell-shaped virtual sources and mixing vectors [see the second to fourth columns of Fig. 5(a)], respectively, at 10% noise. The resolved values become (28.9,27.0,32.9)mm and q2=0.59 from fitting the pair of centrosymmetric virtual source and mixing vector [see the second column of Fig. 5(b)], and (30.3,32.3,26.6)mm and q2=1.33 from fitting the first pair of dumbbell-shaped virtual source and mixing vector [see the third column of Fig. 5(b)], respectively, at 20% noise. The dumbbell-shaped virtual source in the source plane, of the second pair of dumbbell-shaped virtual source and mixing vector, is deformed and the fitting is not shown [see the fourth column of Fig. 5(b)]. The deformation of dumbbell appears first in the source plane with the increase of noise as the grid spacing on the source plane is larger than that in the detector plane in the simulation.

Fig. 5

Same as Fig. 3 with (a) 10% and (b) 20% Gaussian noise. The dumbbell-shaped virtual source on the source plane in the fourth column of (b) is deformed and the least-squares fitting is not shown.

051705_1_028505jbo5.jpg

The error in localization and characterization of scattering inhomogeneities increases rapidly with the increase of noise, from 0.1mm in positioning and 2% in strength at 10% noise to 3mm in positioning and 50% in strength at 20% noise. On the other hand, the effect of noise on localization and characterization of absorptive inhomogeneities is much smaller, the errors at both noise levels are less than 0.5mm in positioning and 2% in strength. The results in Sec. 3(B) and this section are summarized in Table 1 .

Table 1

Comparison of known and OPTICA determined positions and strengths of absorption (Abs) and scattering (Sca) inhomogeneities under different levels of Gaussian noise.

NoiseTargetKnown Position (x,y,z) (mm)KnownStrengthResolved Position (x,y,z) (mm)ResolvedStrengthError inPosition (mm)Error inStrength (%)
5%Abs(50, 60, 20)0.01(50.2, 60.3, 20.2)0.0101 0.3 1
Sca(30, 30, 30) 1 (30.0, 30.0, 30.0) 0.99 0.1 1
10%Abs(50, 60, 20)0.01(50.2, 60.3, 20.1)0.0101 0.3 1
Sca(30, 30, 30) 1 (30.0, 30.1, 30.0) 0.98 0.1 2
20%Abs(50, 60, 20)0.01(50.1, 60.3, 20.5)0.0102 0.5 2
Sca(30, 30, 30) 1 (28.9, 27.0, 32.9) 0.59 3 50

3.5.

Effect of Uncertainty in Background

In the examples discussed here, we have assumed the light intensities change measured on the detector plane is obtained with an exact knowledge about the background. To examine the effect of uncertainty in background optical property on the performance of OPTICA, we model the error in the estimation of the background absorption or diffusion coefficients as a uniform Gaussian random field f(r) . The Gaussian noise addressed in Sec. 3(D) is set to be zero here. OPTICA operates on a “dirty” scattering wave ϕsca(rd,rs)+δϕsca(rd,rs) , where δϕsca(rd,rs) is the change in the scattered wave from that of a uniform background of absorption μa (or diffusion D ) to that of a background of absorption μa+f(r) [or diffusion D+f(r) ]. The magnitude of the background uncertainty is represented by the signal-to-noise ratio (SNR) defined by

Eq. 13

SNR(dB)=10log10rdrsϕsca(rd,rs)2rdrsδϕsca(rd,rs)2.

Figures 6(a), 6(b), 6(c) show the case presented in Fig. 3 of Sec. 3(B) with 40, 34, and 10dB SNR due to background absorption uncertainty, respectively. The resolved absorptive inhomogeneity is at (50.2,60.3,20.1)mm with strength 0.0101 at 40dB SNR, and at (50.2,60.3,20.1)mm with strength 0.0100 at 34dB SNR, and at (50.6,60.3,20.3)mm with strength 0.0090 at 10dB SNR.

Fig. 6

Same as Fig. 3 with (a) 40, (b) 34, and (c) 10dB background absorption uncertainty.

051705_1_028505jbo6.jpg

The resolved position and strength of the scattering object are found to be (30.1,30.1,30.0)mm and q2=0.99 , (32.1,32.9,30.0)mm and q2=0.95 , and (31.4,30.0,27.5)mm and q2=1.01 , respectively, through fitting to the pair of centrosymmetric and two pairs of dumbbell-shaped virtual sources and mixing vectors [see the second to fourth columns of Fig. 5(a)], respectively, at 40dB SNR. The resolved values become (31.6,31.7,25.3)mm and q2=0.52 and (31.7,29.6,31.7)mm and q2=0.78 at 34dB and 10dB SNRs, respectively, from fitting the pair of centrosymmetric virtual source and mixing vector [see the second columns of Figs. 5(b) and 5(c)]. The dumbbell-shaped virtual sources and mixing vectors, especially the dumbbell-shaped virtual sources on the source plane are deformed and the fitting are not shown [see the third and fourth columns of Figs. 5(b) and 5(c)]. The results for the influence of background absorption uncertainty on the performance are summarized in Table 2 .

Table 2

Comparison of known and OPTICA determined positions and strengths of absorption (Abs) and scattering (Sca) inhomogeneities under different levels of background absorption uncertainty.

SNR (dB)TargetKnown Position (x,y,z) (mm)KnownStrengthResolved Position (x,y,z) (mm)ResolvedStrengthError inPosition (mm)Error inStrength (%)
40Abs(50, 60, 20)0.01(50.2, 60.3, 20.1)0.0101 0.3 1
Sca(30, 30, 30) 1 (30.1, 30.1, 30.0) 0.99 0.1 1
34Abs(50, 60, 20)0.01(50.2, 60.3, 20.1)0.0100 0.3 1
Sca(30, 30, 30) 1 (31.6, 31.7, 25.3) 0.52 5 50
10Abs(50, 60, 20)0.01(50.6, 60.3, 20.3)0.0090 0.3 10
Sca(30, 30, 30) 1 (31.7, 29.6, 31.7) 0.78 5 50

Figures 7(a) and 7(b) show the case presented in Fig. 3 of Sec. 3(B) with 34 and 10dB SNR due to background scattering uncertainty. The resolved absorptive inhomogeneity is at (50.1,60.3,20.1)mm with strength 0.0100 at 34dB SNR and at (49.9,60.5,20.1)mm with strength 0.0099 at 10dB SNR.

Fig. 7

Same as Fig. 3 with (a) 34 and (b) 10dB background scattering uncertainty.

051705_1_028505jbo7.jpg

The resolved position and strength of the scattering object are found to be (30.0,30.1,30.0)mm and q2=0.99 , (32.2,33.0,30.0)mm and q2=0.96 , and (32.3,29.3,27.1)mm and q2=1.08 , respectively, through fitting to the pair of centrosymmetric and two pairs of dumbbell-shaped virtual sources and mixing vectors [see the second to fourth columns of Fig. 5(a)], respectively, at 34dB SNR. The resolved position and strength of the scattering object are found to be (31.7,31.1,32.5)mm and q2=0.75 , (30.9,31.4,27.5)mm and q2=1.08 , respectively, through fitting to the pair of centrosymmetric and the first pair of dumbbell-shaped virtual sources and mixing vectors [see the second to fourth columns of Fig. 5(b)] at 10dB SNR. The results for the influence of background scattering uncertainty on the performance are summarized in Table 3 .

Table 3

Comparison of known and OPTICA determined positions and strengths of absorption (Abs) and scattering (Sca) inhomogeneities under different levels of background scattering uncertainty.

SNR (dB)TargetKnown Position (x,y,z) (mm)KnownStrengthResolved Position (x,y,z) (mm)ResolvedStrengthError inPosition (mm)Error inStrength (%)
34Abs(50, 60, 20)0.01(50.1, 60.3, 20.1)0.0100 0.3 1
Sca(30, 30, 30) 1 (30.0, 30.1, 30.0) 0.99 0.1 1
10Abs(50, 60, 20)0.01(49.9, 60.5, 20.1)0.0099 0.5 1
Sca(30, 30, 30) 1 (31.7, 31.1, 32.5) 0.75 2.5 25

The uncertainty in the background absorption or diffusion coefficients affects the performance of OPTICA in a similar fashion as the noise does discussed in Sec. 3(D). The error in localization and characterization of scattering inhomogeneities increases rapidly while the error in localization and characterization of absorptive inhomogeneities only increases mildly with the increase of the uncertainty in the background optical property. The uncertainty in background scattering has a less adverse effect on the performance of OPTICA than that in background absorption.

4.

Discussion

The simulational study of OPTICA presented in this paper demonstrates its potential in optical imaging of objects in turbid media. It was shown to be able to locate and characterize absorptive and scattering inhomogeneities within highly scattering medium. In particular, OPTICA can discriminate between absorptive and scattering inhomogeneities and locate and characterize complex inhomogeneities, which is both absorptive and scattering. The accuracy of localization and characterization of inhomogeneities is high. In the cases investigated for concentrated inhomogeneities within a tissue-emulating slab of thickness of 50mm , the errors in resolved object locations are not greater than 1mm and the errors in the resolved optical strengths are 2% under favorable noise levels and reliable background estimations.

Noise at higher levels and/or larger uncertainty in the optical property of the background medium makes it difficult to discriminate between absorptive and scattering inhomogeneities. In such a situation, other corroborative evidence, such as multiwavelength measurements, are required to determine the nature of inhomogeneities. Noise at higher levels and/or larger uncertainties in the optical property of the background medium also introduces larger errors in localization and characterization of scattering inhomogeneities. The accuracy of localization and characterization of absorptive inhomogeneities is only mildly affected by the amount of noise and/or uncertainty in the range investigated.

OPTICA unmixes inhomogeneities based on the mutual statistical independence between them and does not rely on a Gaussian distribution of light intensity on the surface of the embedding medium. OPTICA has several salient features. First, OPTICA provides the independent components due to the inhomogeneities with minimal processing of the data and does not have to resort to any specific light propagation model for obtaining this information. Specific light propagation models are needed only in the later stage to determine location and optical strength. Second, different geometries, or even an arbitrary shaped boundary, can be used with OPTICA. Although we used the slab geometry in the work reported in this paper, the approach does not depend on any specific geometry. Third, the approach is fast and is amenable to near real-time detection and localization of objects in a turbid medium, which is a key consideration for in vivo medical imaging.

As is well known, the diffusion approximation to RTE, which is widely used in inverse image reconstruction, does not apply when the separation between any two of the source, the inhomogeneity and the detector is small, or when there are clear regions in the medium. A special treatment is also required when the medium has aligned microstructures, such as myofibrils, axons, or collagen fibers in tissues.43 The fact that a prior assumption of a specific model of light propagation in the medium is not assumed in the identification of independent components by ICA and is required only in a Green’s function analysis of the retrieved independent component is desirable, especially in situations that demand a more complex model than the conventional DA. Performing the fitting procedure for each identified independent component is much simpler and more transparent than matching the measured light intensity to a forward model iteratively. The quality of reconstruction of OPTICA is expected to be higher than the conventional approach when only an imperfect forward model is available.

OPTICA is most suited to detect small objects. Given its ability to identify low-contrast small objects, the approach is expected to be especially useful for detection of tumors at their early stages of development.

Acknowledgments

This work is supported in part by U.S. Army Medical Research and Materials Command, Office of Naval Research (ONR), New York State Office of Science, Technology, and Academic Research (NYSTAR), and City University of New York (CUNY) organized research programs. M. X. thanks the support by the Department of Army (Grant No. DAMD17-02-1-0516). M. Alrubaiee thanks the National Science Foundation (NSF) for an Advance Placement Fellowship. We acknowledge Dr. W. Cai for helpful discussions.

References

1. 

Medical Optical Tomography: Functional Imaging and Monitoring, (1993) Google Scholar

2. 

A. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today, 48 38 –40 (1995). 0031-9228 Google Scholar

3. 

M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography,” Opt. Lett., 20 426 –428 (1995). 0146-9592 Google Scholar

4. 

S. K. Gayen and R. R. Alfano, “Emerging optical biomedical imaging techniques,” Opt. Photonics News, 7 17 –22 (1996). 1047-6938 Google Scholar

5. 

J. C. Hebden, S. R. Arridge, and D. T. Delpy, “Optical imaging in medicine: I. Experimental techniques,” Phys. Med. Biol., 42 825 –840 (1997). https://doi.org/10.1088/0031-9155/42/5/007 0031-9155 Google Scholar

6. 

S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol., 42 841 –853 (1997). https://doi.org/10.1088/0031-9155/42/5/008 0031-9155 Google Scholar

7. 

W. Cai, S. K. Gayen, M. Xu, M. Zevallos, M. Alrubaiee, M. Lax, and R. R. Alfano, “Optical tomographic image reconstruction from ultrafast time-sliced transmission measurements,” Appl. Opt., 38 4237 –4246 (1999). 0003-6935 Google Scholar

8. 

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl., 15 R41 –R93 (1999). https://doi.org/10.1088/0266-5611/15/2/022 0266-5611 Google Scholar

9. 

D. Grosenick, H. Wabnitz, H. H. Rinneberg, K. T. Moesta, and P. M. Schlag, “Development of a time-domain optical mammograph and first in vivo applications,” Appl. Opt., 38 2927 –2943 (1999). 0003-6935 Google Scholar

10. 

V. Chernomordik, D. Hattery, A. H. Gandjbakhche, A. Pifferi, P. Taroni, A. Torricelli, G. Valentini, and R. Cubeddu, “Quantification by random walk of the optical parameters of nonlocalized abnormalities embedded within tissuelike phantoms,” Opt. Lett., 25 951 –953 (2000). 0146-9592 Google Scholar

11. 

V. A. Markel and J. C. Schotland, “Inverse scattering for the diffusion equation with general boundary conditions,” Phys. Rev. E, 64 035601 (2001). https://doi.org/10.1103/PhysRevE.64.035601 1063-651X Google Scholar

12. 

A. H. Hielscher and S. Bartel, “Use of penalty terms in gradient-based iterative reconstruction schemes for optical tomography,” J. Biomed. Opt., 6 183 –192 (2001). https://doi.org/10.1117/1.1352753 1083-3668 Google Scholar

13. 

M. Xu, M. Lax, and R. R. Alfano, “Time-resolved Fourier optical diffuse tomography,” J. Opt. Soc. Am. A, 18 1535 –1542 (2001). 0740-3232 Google Scholar

14. 

B. A. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Near-infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE J. Sel. Top. Quantum Electron., 9 199 –209 (2003). https://doi.org/10.1109/JSTQE.2003.813304 1077-260X Google Scholar

15. 

H. Dehghani, B. W. Pogue, S. P. Poplack, and K. D. Paulsen, “Multiwavelength three-dimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results,” Appl. Opt., 42 135 –145 (2003). 0003-6935 Google Scholar

16. 

Appl. Opt., 42 2869 –3329 (2003) (0003-6935) Google Scholar

17. 

W. Cai, M. Xu, and R. R. Alfano, “Three dimensional radiative transfer tomography for turbid media,” IEEE J. Sel. Top. Quantum Electron., 9 189 –198 (2003). 1077-260X Google Scholar

18. 

L. Wang, P. P. Ho, C. Liu, G. Zhang, and R. R. Alfano, “Ballistic 2-D imaging through scattering walls using an ultrafast optical Kerr gate,” Science, 253 769 –771 (1991). 0036-8075 Google Scholar

19. 

R. R. Alfano, X. Liang, L. Wang, and P. Ho, “Time-resolved imaging of translucent droplets in highly scattering media,” Science, 264 1913 –1914 (1994). 0036-8075 Google Scholar

20. 

W. Cai, M. Lax, and R. R. Alfano, “Analytical solution of the elastic Boltzmann transport equation in an infinite uniform medium using cumulant expansion,” J. Phys. Chem. B, 104 3996 –4000 (2000). https://doi.org/10.1021/jp994447+ 1089-5647 Google Scholar

21. 

W. Cai, M. Lax, and R. R. Alfano, “Analytical solution of the polarized photon transport equation in an infinite uniform medium using cumulant expansion,” Phys. Rev. E, 63 016606 (2000). https://doi.org/10.1103/PhysRevE.63.016606 1063-651X Google Scholar

22. 

M. Xu, W. Cai, M. Lax, and R. R. Alfano, “Photon migration in turbid media using a cumulant approximation to radiative transfer,” Phys. Rev. E, 65 066609 (2002). https://doi.org/10.1103/PhysRevE.65.066609 1063-651X Google Scholar

23. 

A. D. Klose and A. Hielscher, “Fluorescence tomography with simulated data based on the equation of radiative transfer,” Opt. Lett., 28 1019 –1021 (2003). 0146-9592 Google Scholar

24. 

A. H. Gandjbakhche, G. H. Weiss, R. F. Bonner, and R. Nossal, “Photon path-length distributions for transmission through optically turbid slabs,” Phys. Rev. E, 48 810 –818 (1993). https://doi.org/10.1103/PhysRevE.48.810 1063-651X Google Scholar

25. 

A. H. Gandjbakhche, V. Chernomordik, J. C. Hebden, and R. Nossal, “Time-dependent contrast functions for quantitative imaging in time-resolved transillumination experiments,” Appl. Opt., 37 1973 –1981 (1998). 0003-6935 Google Scholar

26. 

Ill-Posed Problems in the Natural Sciences, MIR, Moscow (1987). Google Scholar

27. 

P. Comon, “Independent component analysis—a new concept,” Signal Process., 36 287 –314 (1994). https://doi.org/10.1016/0165-1684(94)90029-9 0165-1684 Google Scholar

28. 

A. J. Bell, “Information theory, independent component analysis, and applications,” Unsupervised Adaptive Filtering, 237 –264 Wiley, New York (2000). Google Scholar

29. 

D. Nuzillard and J.-M. Nuzillard, “Application of blind source separation to 1-D and 2-D nuclear magnetic resonance spectroscopy,” IEEE Signal Process. Lett., 5 209 –211 (1998). 1070-9908 Google Scholar

30. 

R. Vigário, J. Särelä, V. Jousmäki, M. Hämäläinen, and E. Oja, “Independent component approach to the analysis of EEG and MEG recordings,” IEEE Trans. Biomed. Eng., 47 589 –593 (2000). https://doi.org/10.1109/10.841330 0018-9294 Google Scholar

31. 

A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component Analysis, Wiley, New York (2001). Google Scholar

32. 

J.-F. Cardoso, “Blind signal separation: statistical principles,” Proc. IEEE, 86 2009 –2025 (1998). https://doi.org/10.1109/5.720250 0018-9219 Google Scholar

33. 

A. Hyvärinen and E. Oja, “Independent component analysis: algorithms and applications,” Neural Networks, 13 411 –430 (2000). https://doi.org/10.1016/S0893-6080(00)00026-5 0893-6080 Google Scholar

34. 

S. Chandrasekhar, Radiative Transfer, Dover, New York (1960). Google Scholar

35. 

M. Xu, W. Cai, M. Lax, and R. R. Alfano, “A photon transport forward model for imaging in turbid media,” Opt. Lett., 26 1066 –1068 (2001). 0146-9592 Google Scholar

36. 

P. M. Morse and H. Feshbach, Method of Theoretical Physics, McGraw-Hill, New York (1953). Google Scholar

37. 

M. Lax, V. Nayaramamurti, and R. C. Fulton, “Classical diffusion photon transport in a slab,” Laser Optics of Condensed Matter, 229 –237 Plenum, New York (1987). Google Scholar

38. 

J. X. Zhu, D. J. Pine, and D. A. Weitz, “Internal reflection of diffusive light in random media,” Phys. Rev. A, 44 3948 –3959 (1991). https://doi.org/10.1103/PhysRevA.44.3948 1050-2947 Google Scholar

39. 

R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromber, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A, 11 2727 –2741 (1994). 0740-3232 Google Scholar

40. 

S. R. Arridge, “Photon-measurement density functions. Part I: Analytic forms,” Appl. Opt., 34 7395 –7409 (1995). 0003-6935 Google Scholar

41. 

S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part 2: Finite-element-method calculations,” Appl. Opt., 34 8026 –8037 (1995). 0003-6935 Google Scholar

42. 

H. Heusmann, J. Kölzer, and G. Mitic, “Characterization of female breasts in vivo by time resolved and spectroscopic measurements in near infrared spectroscopy,” J. Biomed. Opt., 1 425 –434 (1996). https://doi.org/10.1117/1.429785 1083-3668 Google Scholar

43. 

J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E, 68 031908 (2003). https://doi.org/10.1103/PhysRevE.68.031908 1063-651X Google Scholar
©(2005) Society of Photo-Optical Instrumentation Engineers (SPIE)
M. Xu, Mohammad Alrubaiee, Swapan Kumar Gayen, and Robert R. Alfano "Optical imaging of turbid media using independent component analysis: theory and simulation," Journal of Biomedical Optics 10(5), 051705 (1 September 2005). https://doi.org/10.1117/1.2101568
Published: 1 September 2005
Lens.org Logo
CITATIONS
Cited by 13 scholarly publications.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Scattering

Absorption

Independent component analysis

Sensors

Light scattering

Diffusion

Optical imaging

Back to Top