# Diffuse Optical Imaging Using Decomposition Methods

Authors: Binlin Wu; M. Alrubaiee; W. Cai; M. Xu; S. K. Gayen
Journal: International Journal of Optics
 (2012)

---

## Abstract

Diffuse optical imaging (DOI) for detecting and locating targets in a highly scattering turbid medium is treated as a blind source separation (BSS) problem. Three matrix decomposition methods, independent component analysis (ICA), principal component analysis (PCA), and nonnegative matrix factorization (NMF) were used to study the DOI problem. The efficacy of resulting approaches was evaluated and compared using simulated and experimental data. Samples used in the experiments included Intralipid-10% or Intralipid-20% suspension in water as the medium with absorptive or scattering targets embedded.


---

## Body

## 1. Introduction

Diffuse optical imaging (DOI) for detection and retrieval of location information of targets in a highly scattering turbid medium may be treated as a blind source separation (BSS) problem [1, 2]. Various matrix decomposition methods, such as, independent component analysis (ICA) [3], principal component analysis (PCA) [4], and nonnegative matrix factorization (NMF) [5, 6] have been developed for solving the BSS problem and retrieving desired information.Xu et al. adapted ICA of information theory to develop optical tomography using independent component analysis (OPTICA) and demonstrated its application for diffuse imaging of absorptive, scattering, and fluorescent targets [7–11]. ICA assumes the signals from different targets to be independent of each other and optimizes a relevant measure of independence to obtain the ICs associated with different targets. The position coordinates of targets in three dimensions are determined from the individual components separately.PCA assumes that the PCs contributing to the signal areuncorrelatedand explain the most variance in the signal. PCA has been widely used in various applications, such as spectroscopy [12], face recognition [13], and neuroimaging [14]. NMF seeks to factorize a matrix into two nonnegative matrices (component signals and weights) and requires the contributions to signal and the weights of the components to be non-negative. It does not imply any relationship between the components. NMF has also been widely used in biological analysis [15] and spectral analysis [16].The objective of this study is to test and compare the efficacy of these algorithms when used to solve the DOI problem. Results are presented and compared using simulative data and experimental data using absorptive and scattering targets embedded in model scattering media. Our interest in solving the DOI problem derives from the need for a noninvasive modality for detecting, locating, and diagnosing breast tumors in early stages of growth.The remainder of the paper is organized as follows. In Section2, the formalisms of the three methods are introduced. Section 3 evaluates the resulting imaging approaches using simulated data. The approaches are further examined in Section 4 for experimental data acquired using absorptive and scattering targets embedded in model scattering media. Section 5 summarizes and discusses the results.


## 2. Formalism

### 2.1. Blind Source Separation Problem

Blind source separation (BSS), also known as blind signal separation, is a general problem in information theory that seeks to separate different individual signals from the measured signals, which are weighted mixtures of those individual signals. AssumingM individual signals, sj(t), j=1,…M, are linearly mixed instantaneously, the BSS problem is modeled as follows. The dimension of sj(t) is Ns, the number of sampling times. In this study, t will be replaced by spatial positions of the excitation light sources. A total of Nd detectors sense Nd different mixtures of sj(t). The mixture measured by the ith detector can be presented as xi(t)=∑j=1Maijsj(t), or X=AS, in a matrix notation, where A∈RNd×M is a mixing or weighting matrix, S∈RM×Ns, X∈RNd×Ns, and M<min(Ns,Nd). The objective of BSS is to retrieve the signals sj(t) and their weights, aij. ICA, PCA, and NMF are statistical analysis methods used to solve the BSS problem.


### 2.2. Diffuse Optical Imaging Problem

In DOI, one measures the signal at the sample boundary, which includes a weighted mixture of contributions from embedded targets. One uses the diffusion approximation [17–19] of the radiative transfer equation [20, 21] as the forward model to describe light propagation in a highly scattering turbid medium. The perturbation in the light intensity distribution measured on the boundary of the sample due to the presence of the targets (which are localized inhomogeneities in the optical properties within the sample volume) may be written, in the first-order Born approximation, as [22, 23](1)Δϕ(rd,rs)=-∫G(rd,r)δμa(r)cG(r,rs)d3r-∫δD(r)c∇rG(rd,r)⋅∇rG(r,rs)d3r,
where rs, rd, and r are the positions of a source of unit power, detector and target, respectively; G(r,rs) and G(rd,r) are the Green’s functions that describe light propagation from the source to the target and from the target to the detector, respectively; δμa and δD are the differences in absorption coefficient and diffusion coefficient between the targets and the background medium, respectively; and c is the light speed in the medium.A multisource illumination and multidetector signal acquisition scheme is used to acquire light transmitted through a scattering medium. For small absorptive targets, a perturbation data matrix is constructed using-Δϕ for all sources. The elements of the data matrix pertaining to absorptive targets represented by the first term in (1) may be written in a discrete form as(2)Xij=∑m=1MGd(ri,rm)τmGs(rm,rj)(i=1,2,…,Nd;j=1,2,…,Ns),
where ri, rj, and rm are the locations of the ith detector, jth source and mth target, respectively; Ns, Nd, and M are the numbers of sources, detectors, and targets, respectively; τm=δμa(rm)cδVmis the optical absorption strength of the mth target of volume δVm; Gs(rm,rj) and Gd(ri,rm) are the Green’s functions that describe light propagation from jth source to mth target and from mth target to ith detector, respectively. The number of targets is assumed to be less than that of sources and detectors, M<min(Nd,Ns).Themth target may be considered to be a virtual source of strength τmGs(rm,rj) excited by the real light source located at rj. The data matrix X={Xij} may be considered to be a set of combinations of light signals from all virtual sources mixed by a mixing matrix {Gd(ri,rm)}. Therefore, this problem can be treated as a BSS problem.As the second term in (1) suggests, each scattering target is represented by three colocated virtual sources of strength: τm∂pGs(rm,rj), where ∂p=∂/∂p, (p=x,y,z), and τm=δD(rm)cδVm, is the optical scattering strength of the mth target [8]. The mixing matrices become {∂pGd(ri,rm)}  for the three virtual sources generated by the mth target. The elements of the data matrix for scattering targets may be written as(3)Xij=∑m=1M∑p={x,y,z}∂pGd(ri,rm)τm∂pGs(rm,rj).Since one absorptive target is represented by one centrosymmetric virtual source, while three virtual sources (one centrosymmetric and two dumb-bell shaped) represent one scattering target [7, 8], the number and patterns of virtual sources may be used, in favorable situations, to identify the target as absorptive or scattering in nature. In this paper, only small targets are considered since all three algorithms are suited for small targets, and early detection, when the tumors are more amenable to treatment, is of practical interest.


### 2.3. DOI as a BSS Problem

The data matrix for the DOI problem may be written as(4)X=AS=∑m=1MAimSmj,
where A∈RNd×M, S∈RM×Ns, and X∈RNd×Ns. For absorptive targets,(5a)Aim=βmGd(ri,rm),Smj=αmGs(rm,rj),
while for scattering targets,(5b)Aim=βm∂pGd(ri,rm),Smj=αm∂pGs(rm,rj).{Smj}(j=1,2,…,Ns) and {Aim} (i=1,2,…Nd) are two-dimensional intensity distributions on the source and detector planes, respectively. Source and detector planes are the boundaries of the sample through which light enters and exits the sample volume, respectively. The scaling factors βm and αm are related to the target optical strength, τm=αmβm. The location of the target and the scaling factors can be retrieved using a least squares fitting via(6a)argminαm,βm,rm⁡{∑j[αm-1Smj-Gs(rm,rj)]2+∑i[βm-1Aim-Gd(ri,rm)]2},or(6b)argminαm,βm,rm{∑p{∑j[αm-1Smj-∂pGs(rm,rj)]2+∑i[βm-1Aim-∂pGd(ri,rm)]2}},for absorptive and scattering targets, respectively. However, when a scattering target is embedded deep in a turbid medium, only the τm∂zGs(rm,rj) virtual source remains significant. So, only p=z may be used for fitting in (6b) [8].


#### 2.3.1. ICA

OPTICA assumes that the virtual sources areindependent of each other [8]. So, they can be retrieved through an iterative process which seeks to maximize the independence among the components. In practice, the independent components are found by maximizing some measure of non-Gaussianity, such as kurtosis (the fourth-order cumulant), of the unmixed components. A Matlab program for ICA was adopted from http://sccn.ucsd.edu/eeglab/. The location of the target can be retrieved by fitting the independent component intensity distributions (ICIDs) to Green’s functions or derivatives of Green’s functions using (6a) and (6b).


#### 2.3.2. PCA

PCA assumes that the virtual sources are uncorrelated so that the correlation (covariance) between them is ideally zero and minimal in practice. The covariance matrix ofS, cov(S) should be diagonal. The general process of PCA is as follows. The data matrix X=AS+𝒩, where 𝒩 is random noise added to the data, A and S the same as defined in (4). When S is mean centered, elements of the mean-centered matrix S′ are defined as(7a)Smj′=Smj-1Ns∑j=1NsSmj.
Similarly,(7b)Xij′=Xij-1Ns∑j=1NsXij.PCA looks for a matrix P that decomposes X into virtual sources, S=PX. It also holds that S′=PX′, since P is just a rotation matrix which does not change the center of the data.(8)cov⁡(S)=S′S′T=(PX′)(PX′)T=PX′X′TPT=Λ,
where Λ=diag{λ1,λ2,…}. The eigenvalues λm are variances in the covariance matrix. Therefore, X′X′TPT=PTΛ, where PT is orthonormal. PCA is realized by eigenvalue decomposition (EVD) of the covariance matrix of X. The eigenvectors with leading eigenvalues (largest variances) are selected to be the PCs using the L-curve [24].SinceX=PTS≈AS,  A  is determined as a matrix including only PCs. S is calculated as S≈(ATA)-1ATX. Rows of S and columns of A represent principal component intensity distributions (PCIDs) on the source plane and detector plane, respectively and are proportional to the images of the virtual sources projected on the source and detector planes. The target positions are determined using (6a) and (6b).


#### 2.3.3. NMF

NMF is a group of multivariate analysis algorithms that factorize a matrixX into A, and  S:X=AS, A and S are nonnegative [6]. Unlike ICA and PCA, NMF does not imply any relationship between the retrieved components; instead, it just enforces non-negativity constraints on A and S. There are various algorithms developed to solve NMF, such as the multiplicative update method [5] and alternating least squares method [25, 26].In the multiplicative update implementation of NMF,A and S can be found by minimizing the square of Euclidean distance ∥X-AS∥2 as the cost function, where A≥0 and S≥0, using the multiplicative update rule:(9a)Aik⟵Aik(XST)ik(ASST)ik,(9b)Skj⟵Skj(ATX)kj(ATAS)kj.The alternating least squares implementation of NMF uses alternate least squares steps to estimateA (or S), and use that estimate to optimize S (or A),   repeating the alternative steps until the desired optimization is obtained. Nonnegativity is ensured by setting any negative element of A or S equal to 0.An NMF toolbox was obtained fromhttp://cogsys.imm.dtu.dk/toolbox/ to perform NMF computation. A built-in command nnmf  is also available in Matlab (R2011a).NMF algorithm requires that the non-negativity assumption must hold in the problem. In particular, for absorptive targets, whenX is constructed with -Δϕ, τm should be positive, that is, the targets should be more absorbing than the background. If the targets have weaker attenuation properties than the background, X needs to be constructed with +Δϕ instead. For scattering targets,  X should be treated similarly to keep its elements positive.When NMF is applied to a scattering target, only the centrosymmetric component shows up properly, since the other two components have dumb-bell shape which includes negative values [8]. So without any prior knowledge or some other experimental means to assess if the target is absorptive or scattering, NMF may not distinguish between the two possibilities.The decomposition methods can be applied with different sample geometries such as slab and cylindrical geometries, and different measurement domains such as time-resolved domain, frequency domain, and continuous wave (CW). In this paper, Green’s functions for slab geometry [23] with CW measurement were used for simulation and experiments.


## 2.1. Blind Source Separation Problem

Blind source separation (BSS), also known as blind signal separation, is a general problem in information theory that seeks to separate different individual signals from the measured signals, which are weighted mixtures of those individual signals. AssumingM individual signals, sj(t), j=1,…M, are linearly mixed instantaneously, the BSS problem is modeled as follows. The dimension of sj(t) is Ns, the number of sampling times. In this study, t will be replaced by spatial positions of the excitation light sources. A total of Nd detectors sense Nd different mixtures of sj(t). The mixture measured by the ith detector can be presented as xi(t)=∑j=1Maijsj(t), or X=AS, in a matrix notation, where A∈RNd×M is a mixing or weighting matrix, S∈RM×Ns, X∈RNd×Ns, and M<min(Ns,Nd). The objective of BSS is to retrieve the signals sj(t) and their weights, aij. ICA, PCA, and NMF are statistical analysis methods used to solve the BSS problem.


## 2.2. Diffuse Optical Imaging Problem

In DOI, one measures the signal at the sample boundary, which includes a weighted mixture of contributions from embedded targets. One uses the diffusion approximation [17–19] of the radiative transfer equation [20, 21] as the forward model to describe light propagation in a highly scattering turbid medium. The perturbation in the light intensity distribution measured on the boundary of the sample due to the presence of the targets (which are localized inhomogeneities in the optical properties within the sample volume) may be written, in the first-order Born approximation, as [22, 23](1)Δϕ(rd,rs)=-∫G(rd,r)δμa(r)cG(r,rs)d3r-∫δD(r)c∇rG(rd,r)⋅∇rG(r,rs)d3r,
where rs, rd, and r are the positions of a source of unit power, detector and target, respectively; G(r,rs) and G(rd,r) are the Green’s functions that describe light propagation from the source to the target and from the target to the detector, respectively; δμa and δD are the differences in absorption coefficient and diffusion coefficient between the targets and the background medium, respectively; and c is the light speed in the medium.A multisource illumination and multidetector signal acquisition scheme is used to acquire light transmitted through a scattering medium. For small absorptive targets, a perturbation data matrix is constructed using-Δϕ for all sources. The elements of the data matrix pertaining to absorptive targets represented by the first term in (1) may be written in a discrete form as(2)Xij=∑m=1MGd(ri,rm)τmGs(rm,rj)(i=1,2,…,Nd;j=1,2,…,Ns),
where ri, rj, and rm are the locations of the ith detector, jth source and mth target, respectively; Ns, Nd, and M are the numbers of sources, detectors, and targets, respectively; τm=δμa(rm)cδVmis the optical absorption strength of the mth target of volume δVm; Gs(rm,rj) and Gd(ri,rm) are the Green’s functions that describe light propagation from jth source to mth target and from mth target to ith detector, respectively. The number of targets is assumed to be less than that of sources and detectors, M<min(Nd,Ns).Themth target may be considered to be a virtual source of strength τmGs(rm,rj) excited by the real light source located at rj. The data matrix X={Xij} may be considered to be a set of combinations of light signals from all virtual sources mixed by a mixing matrix {Gd(ri,rm)}. Therefore, this problem can be treated as a BSS problem.As the second term in (1) suggests, each scattering target is represented by three colocated virtual sources of strength: τm∂pGs(rm,rj), where ∂p=∂/∂p, (p=x,y,z), and τm=δD(rm)cδVm, is the optical scattering strength of the mth target [8]. The mixing matrices become {∂pGd(ri,rm)}  for the three virtual sources generated by the mth target. The elements of the data matrix for scattering targets may be written as(3)Xij=∑m=1M∑p={x,y,z}∂pGd(ri,rm)τm∂pGs(rm,rj).Since one absorptive target is represented by one centrosymmetric virtual source, while three virtual sources (one centrosymmetric and two dumb-bell shaped) represent one scattering target [7, 8], the number and patterns of virtual sources may be used, in favorable situations, to identify the target as absorptive or scattering in nature. In this paper, only small targets are considered since all three algorithms are suited for small targets, and early detection, when the tumors are more amenable to treatment, is of practical interest.


## 2.3. DOI as a BSS Problem

The data matrix for the DOI problem may be written as(4)X=AS=∑m=1MAimSmj,
where A∈RNd×M, S∈RM×Ns, and X∈RNd×Ns. For absorptive targets,(5a)Aim=βmGd(ri,rm),Smj=αmGs(rm,rj),
while for scattering targets,(5b)Aim=βm∂pGd(ri,rm),Smj=αm∂pGs(rm,rj).{Smj}(j=1,2,…,Ns) and {Aim} (i=1,2,…Nd) are two-dimensional intensity distributions on the source and detector planes, respectively. Source and detector planes are the boundaries of the sample through which light enters and exits the sample volume, respectively. The scaling factors βm and αm are related to the target optical strength, τm=αmβm. The location of the target and the scaling factors can be retrieved using a least squares fitting via(6a)argminαm,βm,rm⁡{∑j[αm-1Smj-Gs(rm,rj)]2+∑i[βm-1Aim-Gd(ri,rm)]2},or(6b)argminαm,βm,rm{∑p{∑j[αm-1Smj-∂pGs(rm,rj)]2+∑i[βm-1Aim-∂pGd(ri,rm)]2}},for absorptive and scattering targets, respectively. However, when a scattering target is embedded deep in a turbid medium, only the τm∂zGs(rm,rj) virtual source remains significant. So, only p=z may be used for fitting in (6b) [8].


### 2.3.1. ICA

OPTICA assumes that the virtual sources areindependent of each other [8]. So, they can be retrieved through an iterative process which seeks to maximize the independence among the components. In practice, the independent components are found by maximizing some measure of non-Gaussianity, such as kurtosis (the fourth-order cumulant), of the unmixed components. A Matlab program for ICA was adopted from http://sccn.ucsd.edu/eeglab/. The location of the target can be retrieved by fitting the independent component intensity distributions (ICIDs) to Green’s functions or derivatives of Green’s functions using (6a) and (6b).


### 2.3.2. PCA

PCA assumes that the virtual sources are uncorrelated so that the correlation (covariance) between them is ideally zero and minimal in practice. The covariance matrix ofS, cov(S) should be diagonal. The general process of PCA is as follows. The data matrix X=AS+𝒩, where 𝒩 is random noise added to the data, A and S the same as defined in (4). When S is mean centered, elements of the mean-centered matrix S′ are defined as(7a)Smj′=Smj-1Ns∑j=1NsSmj.
Similarly,(7b)Xij′=Xij-1Ns∑j=1NsXij.PCA looks for a matrix P that decomposes X into virtual sources, S=PX. It also holds that S′=PX′, since P is just a rotation matrix which does not change the center of the data.(8)cov⁡(S)=S′S′T=(PX′)(PX′)T=PX′X′TPT=Λ,
where Λ=diag{λ1,λ2,…}. The eigenvalues λm are variances in the covariance matrix. Therefore, X′X′TPT=PTΛ, where PT is orthonormal. PCA is realized by eigenvalue decomposition (EVD) of the covariance matrix of X. The eigenvectors with leading eigenvalues (largest variances) are selected to be the PCs using the L-curve [24].SinceX=PTS≈AS,  A  is determined as a matrix including only PCs. S is calculated as S≈(ATA)-1ATX. Rows of S and columns of A represent principal component intensity distributions (PCIDs) on the source plane and detector plane, respectively and are proportional to the images of the virtual sources projected on the source and detector planes. The target positions are determined using (6a) and (6b).


### 2.3.3. NMF

NMF is a group of multivariate analysis algorithms that factorize a matrixX into A, and  S:X=AS, A and S are nonnegative [6]. Unlike ICA and PCA, NMF does not imply any relationship between the retrieved components; instead, it just enforces non-negativity constraints on A and S. There are various algorithms developed to solve NMF, such as the multiplicative update method [5] and alternating least squares method [25, 26].In the multiplicative update implementation of NMF,A and S can be found by minimizing the square of Euclidean distance ∥X-AS∥2 as the cost function, where A≥0 and S≥0, using the multiplicative update rule:(9a)Aik⟵Aik(XST)ik(ASST)ik,(9b)Skj⟵Skj(ATX)kj(ATAS)kj.The alternating least squares implementation of NMF uses alternate least squares steps to estimateA (or S), and use that estimate to optimize S (or A),   repeating the alternative steps until the desired optimization is obtained. Nonnegativity is ensured by setting any negative element of A or S equal to 0.An NMF toolbox was obtained fromhttp://cogsys.imm.dtu.dk/toolbox/ to perform NMF computation. A built-in command nnmf  is also available in Matlab (R2011a).NMF algorithm requires that the non-negativity assumption must hold in the problem. In particular, for absorptive targets, whenX is constructed with -Δϕ, τm should be positive, that is, the targets should be more absorbing than the background. If the targets have weaker attenuation properties than the background, X needs to be constructed with +Δϕ instead. For scattering targets,  X should be treated similarly to keep its elements positive.When NMF is applied to a scattering target, only the centrosymmetric component shows up properly, since the other two components have dumb-bell shape which includes negative values [8]. So without any prior knowledge or some other experimental means to assess if the target is absorptive or scattering, NMF may not distinguish between the two possibilities.The decomposition methods can be applied with different sample geometries such as slab and cylindrical geometries, and different measurement domains such as time-resolved domain, frequency domain, and continuous wave (CW). In this paper, Green’s functions for slab geometry [23] with CW measurement were used for simulation and experiments.


## 2.3.1. ICA

OPTICA assumes that the virtual sources areindependent of each other [8]. So, they can be retrieved through an iterative process which seeks to maximize the independence among the components. In practice, the independent components are found by maximizing some measure of non-Gaussianity, such as kurtosis (the fourth-order cumulant), of the unmixed components. A Matlab program for ICA was adopted from http://sccn.ucsd.edu/eeglab/. The location of the target can be retrieved by fitting the independent component intensity distributions (ICIDs) to Green’s functions or derivatives of Green’s functions using (6a) and (6b).


## 2.3.2. PCA

PCA assumes that the virtual sources are uncorrelated so that the correlation (covariance) between them is ideally zero and minimal in practice. The covariance matrix ofS, cov(S) should be diagonal. The general process of PCA is as follows. The data matrix X=AS+𝒩, where 𝒩 is random noise added to the data, A and S the same as defined in (4). When S is mean centered, elements of the mean-centered matrix S′ are defined as(7a)Smj′=Smj-1Ns∑j=1NsSmj.
Similarly,(7b)Xij′=Xij-1Ns∑j=1NsXij.PCA looks for a matrix P that decomposes X into virtual sources, S=PX. It also holds that S′=PX′, since P is just a rotation matrix which does not change the center of the data.(8)cov⁡(S)=S′S′T=(PX′)(PX′)T=PX′X′TPT=Λ,
where Λ=diag{λ1,λ2,…}. The eigenvalues λm are variances in the covariance matrix. Therefore, X′X′TPT=PTΛ, where PT is orthonormal. PCA is realized by eigenvalue decomposition (EVD) of the covariance matrix of X. The eigenvectors with leading eigenvalues (largest variances) are selected to be the PCs using the L-curve [24].SinceX=PTS≈AS,  A  is determined as a matrix including only PCs. S is calculated as S≈(ATA)-1ATX. Rows of S and columns of A represent principal component intensity distributions (PCIDs) on the source plane and detector plane, respectively and are proportional to the images of the virtual sources projected on the source and detector planes. The target positions are determined using (6a) and (6b).


## 2.3.3. NMF

NMF is a group of multivariate analysis algorithms that factorize a matrixX into A, and  S:X=AS, A and S are nonnegative [6]. Unlike ICA and PCA, NMF does not imply any relationship between the retrieved components; instead, it just enforces non-negativity constraints on A and S. There are various algorithms developed to solve NMF, such as the multiplicative update method [5] and alternating least squares method [25, 26].In the multiplicative update implementation of NMF,A and S can be found by minimizing the square of Euclidean distance ∥X-AS∥2 as the cost function, where A≥0 and S≥0, using the multiplicative update rule:(9a)Aik⟵Aik(XST)ik(ASST)ik,(9b)Skj⟵Skj(ATX)kj(ATAS)kj.The alternating least squares implementation of NMF uses alternate least squares steps to estimateA (or S), and use that estimate to optimize S (or A),   repeating the alternative steps until the desired optimization is obtained. Nonnegativity is ensured by setting any negative element of A or S equal to 0.An NMF toolbox was obtained fromhttp://cogsys.imm.dtu.dk/toolbox/ to perform NMF computation. A built-in command nnmf  is also available in Matlab (R2011a).NMF algorithm requires that the non-negativity assumption must hold in the problem. In particular, for absorptive targets, whenX is constructed with -Δϕ, τm should be positive, that is, the targets should be more absorbing than the background. If the targets have weaker attenuation properties than the background, X needs to be constructed with +Δϕ instead. For scattering targets,  X should be treated similarly to keep its elements positive.When NMF is applied to a scattering target, only the centrosymmetric component shows up properly, since the other two components have dumb-bell shape which includes negative values [8]. So without any prior knowledge or some other experimental means to assess if the target is absorptive or scattering, NMF may not distinguish between the two possibilities.The decomposition methods can be applied with different sample geometries such as slab and cylindrical geometries, and different measurement domains such as time-resolved domain, frequency domain, and continuous wave (CW). In this paper, Green’s functions for slab geometry [23] with CW measurement were used for simulation and experiments.


## 3. Simulation

The sample was considered to be a 40 mm thick uniform scattering slab with lateral dimension of80mm×80mm, as shown in Figure 1. Its absorption and diffusion coefficients were taken to be μa=0.003 mm−1 and D=1/3 mm (transport mean free path, lt=1 mm), respectively, which are similar to the average value of those parameters for human breast tissue. An absorptive and a scattering point targets were placed at (50, 60, 15) mm and (30, 30, 25) mm, respectively. The index of refraction n of the medium was taken to be 1.33. The speed of light is 2.998×108 m/s or 299.8 mm/ns in vacuum, and 225.4 mm/ns in the medium. The absorption coefficient of the absorptive target was set to be higher than the background with Δμa=0.001 mm−1, while the diffusion coefficient was taken to be the same as that of background. The diffusion coefficient of the scattering target was set to be lower than the background (higher scattering coefficient) with ΔD=-0.1 mm (lt=0.7 mm), while the absorption coefficient was taken to be the same as that of the background. The volumes of both targets are set to be 8 mm3. The optical strengths of the absorptive and scattering targets were ΔμacΔV=1.803 mm3/ns and ΔDcΔV=-180.3 mm5/ns, respectively. The incident CW beam step scanned the sample at 21×21 grid points covering an 80×80 mm2 area, with a step size of 4 mm. Light on the opposite side was recorded at 41×41 grid points covering the same area. Multiplicative Gaussian noise of 5% was added to the simulated data. The data matrix X was then obtained using (2) and (3) directly and analyzed using the three different algorithms.Figure 1
Light intensity distribution on the detector plane is recorded when a point source scans on the source plane.


### 3.1. ICA Analysis

One independent component for the absorptive target and three independent components for the scattering target were retrieved by ICA. The independent component intensity distributions (ICIDs) on the detector plane are shown in Figures2(a), 2(c), 2(d), and 2(e). Similar ICIDs were obtained on the source plane. Figure 2(g) shows the centrosymmetric ICID for the scattering target, and Figure 2(i) shows the ICID for the absorptive target, on the source plane.Figure 2
ICA-extracted two-dimensional intensity distribution on the detector plane of: (a) the centrosymmetric component; (c) and (d) dumb-bell shaped components of the scattering target; (e) the absorptive target. Similar intensity distribution on the source plane of: (g) the centrosymmetric component of the scattering target and (i) the absorptive target for comparison. Fits to the spatial intensity profile on the detector plane along the white dashed line (shown in figures) of the centrosymmetric component of the scattering target is shown in (b), and that of the absorptive target is shown in (f). Corresponding fits to spatial profiles on the source plane are displayed in (h) and (j), respectively.The components in either the detector plane or the source plane can, in principle, be used to extract position and optical strength of the target(s). However, in our experimental arrangement signal is collected by a1024×1024 pixels CCD camera, while the source plane is scanned in an x-y array of points, which is much smaller than the number of pixels in the CCD camera. Consequently, the resolution in the detector plane is much better, and the data set more robust than the source side. So, we used the images on the detector plane for retrieving target information using experimental data. While it would not matter in simulation, to be consistent with experimental situations, we employed detector plane images when using simulated data as well for all three algorithms. Table 1 lists the locations and strengths of the absorptive and scattering targets retrieved by fitting the spatial intensity profile of the centrosymmetric components on the detector plane to Green’s functions and derivatives of Green’s functions using (6a) and (6b), respectively, as shown in Figures 2(b) and 2(f). Figures 2(h) and 2(j) show the corresponding fits to the profiles on the source plane.Table 1
Positions and optical strengths retrieved using simulated data and ICA, PCA, and NMF algorithms.
TargetKnown position (mm)AlgorithmFitted position (mm)Error (mm)Known  strength*Fitted strength*Error (%)Sca.(30, 30, 25)ICA(29.9, 30.0, 25.1)(0.1, 0, 0.1)−180.3−179.90.22PCA(30.0, 30.0, 25.0)(0, 0, 0)−180.3−180.10.11NMF(30.0, 30.0, 25.0)(0, 0, 0)−180.3−178.51Abs.(50, 60, 15)ICA(50.1, 60.2, 15.0)(0.1, 0.2, 0)1.8031.8261.28PCA(50.1, 60.1, 14.9)(0.1, 0.1, 0.1)1.8031.8120.5NMF(50.1, 60.1, 15.0)(0.1, 0.1, 0)1.8031.8030*The unit for absorption strength of the target is mm3/ns and for scattering strength is mm5/ns.


### 3.2. PCA Analysis

Eigenvalue equation of the covariance matrix ofX was solved. The eigenvalues found by PCA were sorted in a descending order. Figure 3 shows a plot of leading 20 eigenvalues on a logarithmic scale.Figure 3
A logarithmic plot of the first 20 PCA eigenvalues.First four leading eigenvalues were selected for PCs. The corresponding PCIDs were calculated. The PCIDs on the detector plane are shown in Figure4. Similar images for PCIDs on the source plane were obtained. The scattering target has one centrosymmetric (Figure 4(a)) component and two dumb-bell shaped (Figures 4(c) and 4(d)) components, while the absorptive target has only one component (Figure 4(e)).PCA-extracted two-dimensional intensity distribution on the detector plane of: (a) the centrosymmetric component; and (c) and (d) dumb-bell shaped components of the scattering target; (e) the absorptive target. Green’s function fits to the spatial intensity profiles along the dashed line (shown in figures) of the (b) centrosymmetric component of the scattering target and (f) absorptive target, respectively, to retrieve the locations of the two targets.
(a)(b)(c)(d)(e)(f)Figures4(b) and 4(f) show fits to the spatial intensity profile of the centrosymmetric component of the scattering target and that of the absorptive target, respectively, to retrieve the locations of the two targets. The locations and optical strengths of the targets retrieved by PCA are also shown in Table 1.


### 3.3. NMF Analysis

The mixing matrix and virtual sources were retrieved from the data matrixX using NMF as explained in Section 2.3.3. As in the other two approaches, only one component is extracted for the absorptive target. Since NMF has a non-negativity constraint, only the centrosymmetric component for the scattering target is obtained. Nonnegative component intensity distributions (NCIDs) on detector planes are shown in Figure 5. Similar images for NCIDs on source plane were also obtained using the virtual sources in S. The results are also shown in Table 1.NMF-extracted two-dimensional intensity distribution on the detector plane of: (a) the centrosymmetric component of the scattering target; (c) the absorptive target. Fits to the corresponding spatial intensity profiles along the dashed line (shown in figures) are shown in (b) and (d), respectively.
(a)(b)(c)(d)


### 3.4. Results and Discussion

The positions and optical strengths of the targets retrieved by ICA, PCA, and NMF algorithms are shown in Table1, and compared to the known values. The retrieved results using all three algorithms from this simulated data are in excellent agreement with the known values.


## 3.1. ICA Analysis

One independent component for the absorptive target and three independent components for the scattering target were retrieved by ICA. The independent component intensity distributions (ICIDs) on the detector plane are shown in Figures2(a), 2(c), 2(d), and 2(e). Similar ICIDs were obtained on the source plane. Figure 2(g) shows the centrosymmetric ICID for the scattering target, and Figure 2(i) shows the ICID for the absorptive target, on the source plane.Figure 2
ICA-extracted two-dimensional intensity distribution on the detector plane of: (a) the centrosymmetric component; (c) and (d) dumb-bell shaped components of the scattering target; (e) the absorptive target. Similar intensity distribution on the source plane of: (g) the centrosymmetric component of the scattering target and (i) the absorptive target for comparison. Fits to the spatial intensity profile on the detector plane along the white dashed line (shown in figures) of the centrosymmetric component of the scattering target is shown in (b), and that of the absorptive target is shown in (f). Corresponding fits to spatial profiles on the source plane are displayed in (h) and (j), respectively.The components in either the detector plane or the source plane can, in principle, be used to extract position and optical strength of the target(s). However, in our experimental arrangement signal is collected by a1024×1024 pixels CCD camera, while the source plane is scanned in an x-y array of points, which is much smaller than the number of pixels in the CCD camera. Consequently, the resolution in the detector plane is much better, and the data set more robust than the source side. So, we used the images on the detector plane for retrieving target information using experimental data. While it would not matter in simulation, to be consistent with experimental situations, we employed detector plane images when using simulated data as well for all three algorithms. Table 1 lists the locations and strengths of the absorptive and scattering targets retrieved by fitting the spatial intensity profile of the centrosymmetric components on the detector plane to Green’s functions and derivatives of Green’s functions using (6a) and (6b), respectively, as shown in Figures 2(b) and 2(f). Figures 2(h) and 2(j) show the corresponding fits to the profiles on the source plane.Table 1
Positions and optical strengths retrieved using simulated data and ICA, PCA, and NMF algorithms.
TargetKnown position (mm)AlgorithmFitted position (mm)Error (mm)Known  strength*Fitted strength*Error (%)Sca.(30, 30, 25)ICA(29.9, 30.0, 25.1)(0.1, 0, 0.1)−180.3−179.90.22PCA(30.0, 30.0, 25.0)(0, 0, 0)−180.3−180.10.11NMF(30.0, 30.0, 25.0)(0, 0, 0)−180.3−178.51Abs.(50, 60, 15)ICA(50.1, 60.2, 15.0)(0.1, 0.2, 0)1.8031.8261.28PCA(50.1, 60.1, 14.9)(0.1, 0.1, 0.1)1.8031.8120.5NMF(50.1, 60.1, 15.0)(0.1, 0.1, 0)1.8031.8030*The unit for absorption strength of the target is mm3/ns and for scattering strength is mm5/ns.


## 3.2. PCA Analysis

Eigenvalue equation of the covariance matrix ofX was solved. The eigenvalues found by PCA were sorted in a descending order. Figure 3 shows a plot of leading 20 eigenvalues on a logarithmic scale.Figure 3
A logarithmic plot of the first 20 PCA eigenvalues.First four leading eigenvalues were selected for PCs. The corresponding PCIDs were calculated. The PCIDs on the detector plane are shown in Figure4. Similar images for PCIDs on the source plane were obtained. The scattering target has one centrosymmetric (Figure 4(a)) component and two dumb-bell shaped (Figures 4(c) and 4(d)) components, while the absorptive target has only one component (Figure 4(e)).PCA-extracted two-dimensional intensity distribution on the detector plane of: (a) the centrosymmetric component; and (c) and (d) dumb-bell shaped components of the scattering target; (e) the absorptive target. Green’s function fits to the spatial intensity profiles along the dashed line (shown in figures) of the (b) centrosymmetric component of the scattering target and (f) absorptive target, respectively, to retrieve the locations of the two targets.
(a)(b)(c)(d)(e)(f)Figures4(b) and 4(f) show fits to the spatial intensity profile of the centrosymmetric component of the scattering target and that of the absorptive target, respectively, to retrieve the locations of the two targets. The locations and optical strengths of the targets retrieved by PCA are also shown in Table 1.


## 3.3. NMF Analysis

The mixing matrix and virtual sources were retrieved from the data matrixX using NMF as explained in Section 2.3.3. As in the other two approaches, only one component is extracted for the absorptive target. Since NMF has a non-negativity constraint, only the centrosymmetric component for the scattering target is obtained. Nonnegative component intensity distributions (NCIDs) on detector planes are shown in Figure 5. Similar images for NCIDs on source plane were also obtained using the virtual sources in S. The results are also shown in Table 1.NMF-extracted two-dimensional intensity distribution on the detector plane of: (a) the centrosymmetric component of the scattering target; (c) the absorptive target. Fits to the corresponding spatial intensity profiles along the dashed line (shown in figures) are shown in (b) and (d), respectively.
(a)(b)(c)(d)


## 3.4. Results and Discussion

The positions and optical strengths of the targets retrieved by ICA, PCA, and NMF algorithms are shown in Table1, and compared to the known values. The retrieved results using all three algorithms from this simulated data are in excellent agreement with the known values.


## 4. Experiments

### 4.1. Experimental Materials and Methods

In this Section, the algorithms are evaluated using experimental data for absorptive and scattering targets embedded in model scattering media whose absorption and scattering properties are adjusted to mimic the average values of those parameters for human breast tissues. Two different experiments were carried out with two different samples. The first sample used a250mm×250mm×50mm transparent plastic container filled with Intralipid 10% suspension in water as the background medium. The concentration of Intralipid 10% was adjusted to provide [27, 28] an absorption coefficient of μa~0.003 mm−1, and a transport mean-free path lt~1.43 mm at 785 nm. The second sample used a similar container with dimension of 250mm×250mm×60mm filled with Intralipid 20% suspension in water. The concentration of Intralipid 20% was adjusted to provide [27, 28] μa~0.003 mm−1, and lt~1 mm at 785 nm. These optical parameters of the medium were selected to be similar to the average values of those parameters for human breast tissue. The thickness of the samples was also comparable to that of a typical compressed female human breast.In the first experiment, two absorptive targets were embedded in the medium. The targets were ~10-mm diameter glass spheres filled Indocyanine green (ICG) dye dissolved in Intralipid-20% suspension in water to obtain an absorption coefficientμa=1.15 mm−1 at 785 nm, and to match the background scattering coefficient of 1.97 mm−1. The targets were placed at (57.2, 18.1,   20.0) mm and (19.9, 48.1, 25.0) mm, respectively.In the second experiment, two scattering targets were embedded, which were also ~10 mm diameter glass spheres, filled with Intralipid-20% suspension in water. The transport mean free path,lt was adjusted to be 0.25 mm, with scattering coefficient  μs≈11 mm−1, and absorption coefficient μa same as that of the background medium. The targets were placed in the middle plane (z=30 mm) in the container with a lateral distance of 40 mm from each other (center to center).The experimental setup is schematically shown in Figure6. A 10 mW 785 nm diode laser beam was used to illuminate the first sample, while a 100 mW 785 nm diode laser beam was used for the second sample. The input surface (source plane) of the samples was scanned across the laser beam in an x-y array of grid points to realize the multi-source interrogation of the samples. The transmitted light from the exit surface (detector plane) was recorded by a 1024pixel×1024pixel (pixel size = 24 μm) CCD camera (Photometrics CH350) equipped with a 60 mm focal-length camera lens. Each pixel of the CCD camera can be considered to be a detector implementing the multidetector signal acquisition arrangement. A set of 16 bit 1024pixel×1024pixel images were acquired. The two samples were scanned in an array of 11×12 and 11×15 grid points, respectively, with a step size of 5 mm in both cases. The processes of scanning and data acquisition were controlled by a personal computer. At all scan positions, raw transillumination images of the samples were recorded by the computer for further analysis.Figure 6
A schematic diagram of the experimental arrangement used for imaging objects embedded in a turbid medium. The inset at the bottom shows the 2D array in the input plane that was scanned across the incident laser beam; the inset to the right shows a typical raw image recorded by the CCD. (CCD: charge coupled device, and PC: personal computer).


### 4.2. Analysis and Results

A region of interest (ROI) was cropped out from each image. Then, every5×5 pixels in each cropped image were binned to one pixel to enhance signal-to-noise ratio. A background image was generated by calculating an average image for all scan positions to approximate the transillumination image without target(s) embedded.This averaging method for generating background image is suitable for small targets used in our experiments, as the ratio of the volume of the sample to that of the target was quite high (~500 : 1). Forin vivo imaging of tumors in early stages of growth, the breast-to-tumor volume ratio will be similarly high, and the averaging method will be applicable. Alternative approaches for generating a background image include using image of (a) a phantom that has the same average optical properties as the sample [29]; (b) the healthy contralateral breast for breast imaging [30]; (c) the sample obtained using light of wavelength for which the target(s) and the background have identical optical properties [31]. Still another approach is to compute the background using an appropriate forward model [32]. A more detailed discussion of this important issue appears in one of our earlier publications [33].The background image was also cropped and binned corresponding to the ROI for each scan position. Perturbation in the light intensity distribution,Δϕ due to targets in each image was found by subtracting the background image from the image. The data matrix X was then constructed using the light intensity perturbations at all scan positions. ICA, PCA, and NMF decomposition algorithms were performed on the data matrix separately. Results are shown and discussed below.


#### 4.2.1. Absorptive Targets

The images on the detector plane obtained using the ICA, PCA, and NMF algorithms are shown in Figures7, 8, and 9, respectively. Similar images on the source plane were also obtained using all three algorithms. The right side of each figure shows the corresponding spatial intensity profile. Locations of the targets are extracted from fits to these spatial intensity profiles, as described in Section 2.3 using (6a) and (6b). The results are presented in Table 2. In Figure 7, images on the source plane are shown in (e) and (g), and Green’s function fits to their spatial profiles are shown in (f) and (h) for comparison.Table 2
Known positions versus retrieved positions of the absorptive targets using ICA, PCA, and NMF algorithms.
TargetKnown position  (mm)AlgorithmFitted position (mm)Error (mm)1(57.2, 18.1, 20)ICA(57.4, 18.2, 21.5)(0.2, 0.1, 1.5)PCA(57.4, 18.2, 20.6)(0.2, 0.1, 0.6)NMF(57.4, 18.2, 19.5)(0.2, 0.1, 0.5)2(19.9, 48.1, 25)ICA(18.2, 46.7, 24.7)(1.7, 1.4, 0.3)PCA(18.2, 47.6, 25.9)(1.7, 0.5, 0.9)NMF(18.2, 47.6, 23.3)(1.7, 0.5, 1.7)Figure 7
ICA-generated ICIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed lines are shown in (b) and (d). ICIDs on the source plane are shown in (e) and (g); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (f) and (h).PCIDs on the detector plane are shown in (a) and (c); and corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)NCIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)It follows from the comparison of the results in Table2 that the positions retrieved by all three algorithms are in good agreement with the known positions. The errors in the retrieved locations (x,y,z) of the two targets were within 1.7 mm. The PCIDs were not totally separated. Some “residue” was observed in one PCID from the other. ICA and NMF separated two components from this dataset more clearly.


#### 4.2.2. Scattering Targets

The “images” corresponding to the centrosymmetric components of the virtual sources (targets) on the detector plane obtained using the ICA, PCA, and NMF algorithms are shown in Figures10, 11, and 12, respectively. Similar images on the source plane were also obtained. The right side of each figure shows the corresponding spatial intensity profile. Locations of the targets are extracted from fits to these spatial intensity profiles, as described in Section 2.3 using (6a) and (6b). The results are presented in Table 3.Table 3
Known positions versus retrieved positions of the scattering targets using ICA, PCA, and NMF algorithms.
TargetKnown position  (mm)AlgorithmFitted position  (mm)Error (mm)1(13.0, 28.0, 30.0)ICA(12.6, 28.7, 29.1)(0.4, 0.7, 0.9)PCA(12.6, 28.7, 28.6)(0.4, 0.7, 1.4)NMF(12.0, 28.5, 33.0)(1.0, 0.5, 3.0)2(53.3, 28.5, 30.0)ICA(51.0, 31.8, 26.8)(2.3, 3.3, 3.2)PCA(50.9, 31.8, 26.7)(2.4, 3.3, 3.3)NMF(53.3, 28.0, 30.3)(0.0, 0.5, 0.3)ICA-generated ICIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)PCIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)NCIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)Both targets were detected by all three algorithms. The target locations retrieved by three algorithms are shown in Table3 and compared with known locations. Overall, all three algorithms detect and locate the scattering and the absorptive targets with good accuracy, the maximum deviation of any one coordinate from the known value being ~3 mm. Since the maximum difference between the known and retrieved position coordinates was larger for the scattering targets, we calculated the squared correlation coefficient γ to assess the fitting quality. NMF retrieves the position coordinates better (within 0.5 mm) for the scattering Target 2 than done by ICA and PCA (deviation from known values being between 2-3 mm). NMF retrieved the position coordinates for Target 1 with 3.0 mm error in z direction, which is not as good as that done by ICA and PCA. But γ is 0.783 and 0.778 in the fittings for ICA and PCA, respectively, as compared to 0.993 for NMF, indicating that the quality of the fitting is better for NMF. The quality of fitting is presumably affected by the efficacy of decomposition. The decomposed NCIDs by NMF were more “clean” than those decomposed by ICA and PCA. We ascribe the observed higher errors in ICA and PCA estimates of the position coordinates of the scattering Target 2 than the NMF estimates to the interference from the other virtual source (corresponding to Target 1) in ICA (Figure 10(c)) and PCA (Figure 11(c)) images. It is commonly believed that errors in locating a scattering target are higher than that for locating an absorptive target, and the results of this study conform to that notion.


## 4.1. Experimental Materials and Methods

In this Section, the algorithms are evaluated using experimental data for absorptive and scattering targets embedded in model scattering media whose absorption and scattering properties are adjusted to mimic the average values of those parameters for human breast tissues. Two different experiments were carried out with two different samples. The first sample used a250mm×250mm×50mm transparent plastic container filled with Intralipid 10% suspension in water as the background medium. The concentration of Intralipid 10% was adjusted to provide [27, 28] an absorption coefficient of μa~0.003 mm−1, and a transport mean-free path lt~1.43 mm at 785 nm. The second sample used a similar container with dimension of 250mm×250mm×60mm filled with Intralipid 20% suspension in water. The concentration of Intralipid 20% was adjusted to provide [27, 28] μa~0.003 mm−1, and lt~1 mm at 785 nm. These optical parameters of the medium were selected to be similar to the average values of those parameters for human breast tissue. The thickness of the samples was also comparable to that of a typical compressed female human breast.In the first experiment, two absorptive targets were embedded in the medium. The targets were ~10-mm diameter glass spheres filled Indocyanine green (ICG) dye dissolved in Intralipid-20% suspension in water to obtain an absorption coefficientμa=1.15 mm−1 at 785 nm, and to match the background scattering coefficient of 1.97 mm−1. The targets were placed at (57.2, 18.1,   20.0) mm and (19.9, 48.1, 25.0) mm, respectively.In the second experiment, two scattering targets were embedded, which were also ~10 mm diameter glass spheres, filled with Intralipid-20% suspension in water. The transport mean free path,lt was adjusted to be 0.25 mm, with scattering coefficient  μs≈11 mm−1, and absorption coefficient μa same as that of the background medium. The targets were placed in the middle plane (z=30 mm) in the container with a lateral distance of 40 mm from each other (center to center).The experimental setup is schematically shown in Figure6. A 10 mW 785 nm diode laser beam was used to illuminate the first sample, while a 100 mW 785 nm diode laser beam was used for the second sample. The input surface (source plane) of the samples was scanned across the laser beam in an x-y array of grid points to realize the multi-source interrogation of the samples. The transmitted light from the exit surface (detector plane) was recorded by a 1024pixel×1024pixel (pixel size = 24 μm) CCD camera (Photometrics CH350) equipped with a 60 mm focal-length camera lens. Each pixel of the CCD camera can be considered to be a detector implementing the multidetector signal acquisition arrangement. A set of 16 bit 1024pixel×1024pixel images were acquired. The two samples were scanned in an array of 11×12 and 11×15 grid points, respectively, with a step size of 5 mm in both cases. The processes of scanning and data acquisition were controlled by a personal computer. At all scan positions, raw transillumination images of the samples were recorded by the computer for further analysis.Figure 6
A schematic diagram of the experimental arrangement used for imaging objects embedded in a turbid medium. The inset at the bottom shows the 2D array in the input plane that was scanned across the incident laser beam; the inset to the right shows a typical raw image recorded by the CCD. (CCD: charge coupled device, and PC: personal computer).


## 4.2. Analysis and Results

A region of interest (ROI) was cropped out from each image. Then, every5×5 pixels in each cropped image were binned to one pixel to enhance signal-to-noise ratio. A background image was generated by calculating an average image for all scan positions to approximate the transillumination image without target(s) embedded.This averaging method for generating background image is suitable for small targets used in our experiments, as the ratio of the volume of the sample to that of the target was quite high (~500 : 1). Forin vivo imaging of tumors in early stages of growth, the breast-to-tumor volume ratio will be similarly high, and the averaging method will be applicable. Alternative approaches for generating a background image include using image of (a) a phantom that has the same average optical properties as the sample [29]; (b) the healthy contralateral breast for breast imaging [30]; (c) the sample obtained using light of wavelength for which the target(s) and the background have identical optical properties [31]. Still another approach is to compute the background using an appropriate forward model [32]. A more detailed discussion of this important issue appears in one of our earlier publications [33].The background image was also cropped and binned corresponding to the ROI for each scan position. Perturbation in the light intensity distribution,Δϕ due to targets in each image was found by subtracting the background image from the image. The data matrix X was then constructed using the light intensity perturbations at all scan positions. ICA, PCA, and NMF decomposition algorithms were performed on the data matrix separately. Results are shown and discussed below.


### 4.2.1. Absorptive Targets

The images on the detector plane obtained using the ICA, PCA, and NMF algorithms are shown in Figures7, 8, and 9, respectively. Similar images on the source plane were also obtained using all three algorithms. The right side of each figure shows the corresponding spatial intensity profile. Locations of the targets are extracted from fits to these spatial intensity profiles, as described in Section 2.3 using (6a) and (6b). The results are presented in Table 2. In Figure 7, images on the source plane are shown in (e) and (g), and Green’s function fits to their spatial profiles are shown in (f) and (h) for comparison.Table 2
Known positions versus retrieved positions of the absorptive targets using ICA, PCA, and NMF algorithms.
TargetKnown position  (mm)AlgorithmFitted position (mm)Error (mm)1(57.2, 18.1, 20)ICA(57.4, 18.2, 21.5)(0.2, 0.1, 1.5)PCA(57.4, 18.2, 20.6)(0.2, 0.1, 0.6)NMF(57.4, 18.2, 19.5)(0.2, 0.1, 0.5)2(19.9, 48.1, 25)ICA(18.2, 46.7, 24.7)(1.7, 1.4, 0.3)PCA(18.2, 47.6, 25.9)(1.7, 0.5, 0.9)NMF(18.2, 47.6, 23.3)(1.7, 0.5, 1.7)Figure 7
ICA-generated ICIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed lines are shown in (b) and (d). ICIDs on the source plane are shown in (e) and (g); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (f) and (h).PCIDs on the detector plane are shown in (a) and (c); and corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)NCIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)It follows from the comparison of the results in Table2 that the positions retrieved by all three algorithms are in good agreement with the known positions. The errors in the retrieved locations (x,y,z) of the two targets were within 1.7 mm. The PCIDs were not totally separated. Some “residue” was observed in one PCID from the other. ICA and NMF separated two components from this dataset more clearly.


### 4.2.2. Scattering Targets

The “images” corresponding to the centrosymmetric components of the virtual sources (targets) on the detector plane obtained using the ICA, PCA, and NMF algorithms are shown in Figures10, 11, and 12, respectively. Similar images on the source plane were also obtained. The right side of each figure shows the corresponding spatial intensity profile. Locations of the targets are extracted from fits to these spatial intensity profiles, as described in Section 2.3 using (6a) and (6b). The results are presented in Table 3.Table 3
Known positions versus retrieved positions of the scattering targets using ICA, PCA, and NMF algorithms.
TargetKnown position  (mm)AlgorithmFitted position  (mm)Error (mm)1(13.0, 28.0, 30.0)ICA(12.6, 28.7, 29.1)(0.4, 0.7, 0.9)PCA(12.6, 28.7, 28.6)(0.4, 0.7, 1.4)NMF(12.0, 28.5, 33.0)(1.0, 0.5, 3.0)2(53.3, 28.5, 30.0)ICA(51.0, 31.8, 26.8)(2.3, 3.3, 3.2)PCA(50.9, 31.8, 26.7)(2.4, 3.3, 3.3)NMF(53.3, 28.0, 30.3)(0.0, 0.5, 0.3)ICA-generated ICIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)PCIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)NCIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)Both targets were detected by all three algorithms. The target locations retrieved by three algorithms are shown in Table3 and compared with known locations. Overall, all three algorithms detect and locate the scattering and the absorptive targets with good accuracy, the maximum deviation of any one coordinate from the known value being ~3 mm. Since the maximum difference between the known and retrieved position coordinates was larger for the scattering targets, we calculated the squared correlation coefficient γ to assess the fitting quality. NMF retrieves the position coordinates better (within 0.5 mm) for the scattering Target 2 than done by ICA and PCA (deviation from known values being between 2-3 mm). NMF retrieved the position coordinates for Target 1 with 3.0 mm error in z direction, which is not as good as that done by ICA and PCA. But γ is 0.783 and 0.778 in the fittings for ICA and PCA, respectively, as compared to 0.993 for NMF, indicating that the quality of the fitting is better for NMF. The quality of fitting is presumably affected by the efficacy of decomposition. The decomposed NCIDs by NMF were more “clean” than those decomposed by ICA and PCA. We ascribe the observed higher errors in ICA and PCA estimates of the position coordinates of the scattering Target 2 than the NMF estimates to the interference from the other virtual source (corresponding to Target 1) in ICA (Figure 10(c)) and PCA (Figure 11(c)) images. It is commonly believed that errors in locating a scattering target are higher than that for locating an absorptive target, and the results of this study conform to that notion.


## 4.2.1. Absorptive Targets

The images on the detector plane obtained using the ICA, PCA, and NMF algorithms are shown in Figures7, 8, and 9, respectively. Similar images on the source plane were also obtained using all three algorithms. The right side of each figure shows the corresponding spatial intensity profile. Locations of the targets are extracted from fits to these spatial intensity profiles, as described in Section 2.3 using (6a) and (6b). The results are presented in Table 2. In Figure 7, images on the source plane are shown in (e) and (g), and Green’s function fits to their spatial profiles are shown in (f) and (h) for comparison.Table 2
Known positions versus retrieved positions of the absorptive targets using ICA, PCA, and NMF algorithms.
TargetKnown position  (mm)AlgorithmFitted position (mm)Error (mm)1(57.2, 18.1, 20)ICA(57.4, 18.2, 21.5)(0.2, 0.1, 1.5)PCA(57.4, 18.2, 20.6)(0.2, 0.1, 0.6)NMF(57.4, 18.2, 19.5)(0.2, 0.1, 0.5)2(19.9, 48.1, 25)ICA(18.2, 46.7, 24.7)(1.7, 1.4, 0.3)PCA(18.2, 47.6, 25.9)(1.7, 0.5, 0.9)NMF(18.2, 47.6, 23.3)(1.7, 0.5, 1.7)Figure 7
ICA-generated ICIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed lines are shown in (b) and (d). ICIDs on the source plane are shown in (e) and (g); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (f) and (h).PCIDs on the detector plane are shown in (a) and (c); and corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)NCIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)It follows from the comparison of the results in Table2 that the positions retrieved by all three algorithms are in good agreement with the known positions. The errors in the retrieved locations (x,y,z) of the two targets were within 1.7 mm. The PCIDs were not totally separated. Some “residue” was observed in one PCID from the other. ICA and NMF separated two components from this dataset more clearly.


## 4.2.2. Scattering Targets

The “images” corresponding to the centrosymmetric components of the virtual sources (targets) on the detector plane obtained using the ICA, PCA, and NMF algorithms are shown in Figures10, 11, and 12, respectively. Similar images on the source plane were also obtained. The right side of each figure shows the corresponding spatial intensity profile. Locations of the targets are extracted from fits to these spatial intensity profiles, as described in Section 2.3 using (6a) and (6b). The results are presented in Table 3.Table 3
Known positions versus retrieved positions of the scattering targets using ICA, PCA, and NMF algorithms.
TargetKnown position  (mm)AlgorithmFitted position  (mm)Error (mm)1(13.0, 28.0, 30.0)ICA(12.6, 28.7, 29.1)(0.4, 0.7, 0.9)PCA(12.6, 28.7, 28.6)(0.4, 0.7, 1.4)NMF(12.0, 28.5, 33.0)(1.0, 0.5, 3.0)2(53.3, 28.5, 30.0)ICA(51.0, 31.8, 26.8)(2.3, 3.3, 3.2)PCA(50.9, 31.8, 26.7)(2.4, 3.3, 3.3)NMF(53.3, 28.0, 30.3)(0.0, 0.5, 0.3)ICA-generated ICIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)PCIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)NCIDs on the detector plane are shown in (a) and (c); corresponding Green’s function fits to the horizontal spatial profiles through the dashed line are shown in (b) and (d).
(a)(b)(c)(d)Both targets were detected by all three algorithms. The target locations retrieved by three algorithms are shown in Table3 and compared with known locations. Overall, all three algorithms detect and locate the scattering and the absorptive targets with good accuracy, the maximum deviation of any one coordinate from the known value being ~3 mm. Since the maximum difference between the known and retrieved position coordinates was larger for the scattering targets, we calculated the squared correlation coefficient γ to assess the fitting quality. NMF retrieves the position coordinates better (within 0.5 mm) for the scattering Target 2 than done by ICA and PCA (deviation from known values being between 2-3 mm). NMF retrieved the position coordinates for Target 1 with 3.0 mm error in z direction, which is not as good as that done by ICA and PCA. But γ is 0.783 and 0.778 in the fittings for ICA and PCA, respectively, as compared to 0.993 for NMF, indicating that the quality of the fitting is better for NMF. The quality of fitting is presumably affected by the efficacy of decomposition. The decomposed NCIDs by NMF were more “clean” than those decomposed by ICA and PCA. We ascribe the observed higher errors in ICA and PCA estimates of the position coordinates of the scattering Target 2 than the NMF estimates to the interference from the other virtual source (corresponding to Target 1) in ICA (Figure 10(c)) and PCA (Figure 11(c)) images. It is commonly believed that errors in locating a scattering target are higher than that for locating an absorptive target, and the results of this study conform to that notion.


## 5. Summary and Discussion

Diffusive optical imaging was modeled as a BSS problem. ICA, PCA, and NMF were used to decompose the data matrix and locate the targets embedded in a highly scattering turbid medium. Only the components corresponding to the targets were extracted from a large dataset for target detection and localization.It may be instructive to compare the objectives, scope, and computational complexity of these decomposition methods with model-based reconstruction methods. Decomposition methods obtain the 3D locations of targets (the number of targets is generally small). Based on the retrieved locations, the methods may then be further extended to retrieve size and optical property information of the targets [9]. The common practice of model-based inverse reconstruction methods is to discretize the sample volume into N×N×N voxels and estimate absorption and/or scattering coefficient in each voxel iteratively. Voxels with significantly different optical properties than the surrounding are regions of interest and may be identified as targets. While estimating the optical properties, the forward model is solved repeatedly to calculate the intensity of the multiply scattered light on the sample boundary. The difference between the intensity of the multiply scattered light predicted by the forward model and the experimental measurements is minimized by seeking an optimal set of the optical properties of every voxel in the sample volume. The number of variables is thus, on the order of  N3. To determine location(s) of target(s) in three dimensions, the decomposition methods process the data matrix to retrieve the main components (A and S). Here, A and S are two-dimensional matrices with the number of unknowns on the order of N2. The number of unknowns is, hence, reduced N times in the decomposition methods compared to the model-based approaches, which leads to a substantial saving in the computational time when N is large. No repeated solution of the forward model is involved in decomposition methods. Consequently, decomposition methods are considerably faster.A comparison of the computational complexity of these two types of approaches may shed further light on their relative computation economy. For a model-based iterative reconstruction method, an equation of the formb=Wx is solved to find the targets, where W is a weight matrix of sizeNdNs×Nv,Nd,Ns, and Nv are the numbers of detectors, sources and voxels, respectively, b is an NdNs×1 vector describing the perturbation in the detected light intensity due to the presence of targets, and x is the perturbation in the optical properties from the background values with dimension of Nv×1. The computational complexity is typically O(NdNsNv2) for a single iteration. For the decomposition approach, b is written as a 2D matrix X with dimension Nd×Ns. To decompose matrix X, the computational complexity per iteration is typically of order O(NdNk) for ICA [34], and O(NdNsNk) for NMF [16], where Nk is the number of components that relates to the number of targets and is usually a small number. For PCA using SVD, the complexity is O(Ns2Nk) [34]. The computational complexity of the intrinsic iterative process involved in the matrix decomposition algorithms is much lower than that in the model-based inverse reconstruction methods.All three matrix decomposition methods presented in this manuscript can potentially be used inin-vivo real-time breast cancer imaging. The three algorithms have different assumptions, which may lead to different favored conditions. In this study, the algorithms were evaluated using simulative and experimental data using model scattering media and absorptive and scattering targets. The (x,y,z) positions of the targets were retrieved with good accuracy. The decomposition provided by ICA is “cleaner” than that of the PCA. PCA did not clearly separate the two absorptive targets used in the first experiment. NMF decomposition seems to provide residue-free “cleaner” images than the other two methods in this study. However, since NMF is based on nonnegativity assumption, the results might deteriorate when such a non-negativity assumption does not hold well. While continuous wave measurements were used in the work presented in this paper, the approaches could be used with frequency domain and time domain measurements as well.The work presented here focuses on detecting and locating small targets, which derive impetus from the need to detect tumors in early stages of growth when those are more amenable to treatment. All three methods are applicable for extended targets as well and are expected to provide the “center of optical strength” as the location of the target.All three approaches are applicable for both scattering and absorbing targets and may be used in clinical setting. The contrast between a tumor and surrounding normal tissue can be due to differences in absorption, scattering, or both absorption and scattering properties and may depend significantly on the wavelength of light used. However,a priori knowledge of the optical characteristics (absorptive or scattering) is not crucial. As has been shown in (2) and (3), the expression for elements of the data matrix for absorptive targets involves Green’s Functions G, while that for scattering targets involves ∂G/∂z≈-κG,where κ=μa/D in CW [9]. This relationship with G provides basis for detection and localization of target(s), whether contrast is due to absorption, scattering, or both. We are using transillumination geometry, which is one of the approaches used by other researchers, and adequate signal for in vivo breast imaging is obtained [29, 35–38].In this paper, we presented results when the approaches were used to detect and obtain three-dimensional location information of the targets. We have demonstrated, while developing OPTICA  that a backprojection formalism can be further implemented to get a cross-section image of the target [11], or the retrieved target locations can be fed into other DOI methods as a priori information to get three-dimensional tomographic images. Since the approaches are suited for small targets, this hold promise for detecting and locating breast tumors in early stages of growth, which is crucially important for effective treatment. Further work involving ex vivo (model) and in vivo imaging of cancerous breast will be needed to establish the full potential of these approaches.



---

*Source: 185435-2012-04-05.xml*