Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Tumor localization using diffuse optical tomography and linearly constrained minimum variance beamforming

Open Access Open Access

Abstract

We present a tumor localization method for diffuse optical tomography using linearly constrained minimum variance (LCMV) beam-forming. Beamforming is a spatial filtering technique where signals from certain directions can be enhanced while noise and interference from other directions are suppressed. In our method, we tessellate the domain into small voxels and regard each voxel as a possible position of abnormality (e.g., tumor). We then design a spatial filter based on the linearly constrained minimum variance criterion and apply it to each voxel in the domain. The abnormality is localized by observing the peak in the filter output signals. We test our method using simulated 3D examples. We assume a cubic transmission geometry and consider different cases where the abnormality is an absorber, a scatterer, and both. We also give examples showing the resolution of our method and its performance under different perturbation levels and noise levels. Simulation results show that LCMV beamforming can localize the abnormality well with good computational efficiency. It can be used alone for tumor localization and also as an effective preprocessing tool for improving the image reconstruction performances of other inverse methods.

©2007 Optical Society of America

1. Introduction

Diffuse optical tomography (DOT) is becoming a useful complement to the current non-invasive imaging modalities. It has a temporal resolution an order of magnitude faster than fMRI, with the advantages of being portable, low-cost, and non-inonizing. Furthermore, it offers unique physiological information about oxy-haemoglobin and deoxy-haemoglobin concentrations that are unavailable from other imaging methods [1]. For these reasons, DOT is being used in areas such as functional brain imaging [2, 3] and optical mammography [4, 5].

The forward model in DOT describes the photon propagation in tissue and is usually given by the diffusion approximation (DA) of the Boltzmann transport equation [6]. The inverse problem involves reconstructing the spatially varying absorption and scattering coefficients of human tissue. It is typically ill-posed and underdetermined because of the large number of unknown parameters and limited number of measurements. Different methods, such as (truncated) singular value decomposition [7], algebraic reconstructed techniques [8], and regularized minimum norm methods [9], have been proposed to solve the inverse problem. In order to improve the resolution of the reconstructed image, prior information has also been incorporated, leading to Bayesian-based estimation methods [10, 11].

Solving the inverse problem is usually computationally expensive, especially in 3D. To ameliorate this problem,a scanning method has been proposed in [12, 13], where data from different source-detector geometries are used to evaluate the depth inclusion, and lateral coordinates are determined by the position of maximal contrast. This method does not involve intensive computation but generally underestimates the difference in optical properties between abnormalities and surrounding tissue [4]. In this paper, we propose a beamforming technique as an alternative to overcome the computational issue associated with the inverse problem. Beamforming is essentially a spatial filter widely used in radar and sonar signal processing [14]. A beamformer operates on the output of an array of sensors in order to enhance the amplitude of a coherent wavefront relative to the background noise and directional interference. In medical imaging applications, beamforming has been used in MEG source imaging [15, 16], EEG source localization [17, 18], and ultrasonic image processing [19]. In all these applications, the beamformer provides a spatial scanning scheme where the source (e.g., a dipole current in the brain) is localized by observing the peak in the filter output signal. In the field of DOT, the location of a tumor in the breast or an activation spot in the brain is determined by observing the distributions of the absorption and scattering coefficients. In this case, we tessellate the whole domain into many small voxels and regard each voxel as a possible location for an abnormality where the optical properties are different from the background tissue. We then design a spatial filter based on the linearly constrained minimum variance beamforming, which scans the domain on a voxel basis. The filter output reflects the contribution of the optical parameters in each voxel to the measurement, and the peak in the filter’s output signal indicates the location of the abnormality. Note that at the current stage, we regard the inverse problem more as a source localization rather than an image reconstruction, since the focus here is to determine the locations of the abnormalities. Since our method scans the domain voxel by voxel instead of trying to recover all the voxels in the whole domain all at once, it is not restricted by the number of voxels and is fast to implement. Thus, it can be used alone to localize the tumor or as a preprocessing tool to provide prior information for other image reconstruction methods (e.g., to help choose the initial searching vector for some iterative inverse approaches [8]).

This paper is organized as follows: In Section 2, we describe the forward and measurement models in DOT, assuming a first-order perturbation approach and Rytov approximation to DA. In Section 3, we introduce the linearly constrained minimum variance beamforming and its application in DOT to solve the inverse problem. Numerical examples are given in Section 4, and Section 5 concludes the paper.

2. Forward and measurement models

In this section, we first briefly describe the forward model based on the Rytov approximation of the diffusion approximation, assuming both absorbing and scattering heterogeneities. See [7], [20, 21] for more detailed derivations and discussions on Rytov approximation. We then construct the measurement model by considering additive Gaussian noise.

2.1. Forward model

In DOT, the diffusion approximation (DA) of the Boltzmann transport equation is well accepted for modeling the photon propagation in tissue. In the frequency domain, DA is given by [22]

.D(r)Urω+cμa(r)Urω+jωUrω=cq0rω,

where U is the fluence (W/m2), μ a the absorption coefficient (m-1), D the diffusion coefficient (m2/s), c the speed of light in the medium (m/s), and q 0 the source (W/m3). The term D is defined as D = c/3(μ a + μ ś), where μ ś = (1 - g)μ s is the reduced scattering coefficient (m-1) with μ s representing the scattering coefficient (m-1), and g the average of the cosine of the scattering angle. Since μ aμ ś for most biological tissues [6], we assume D = 1/3μ ś in our calculations.

In this paper, we consider heterogeneous scattering and absorption coefficients in the domain. Using a first-order perturbation approach, we can write

μa(r)=μa0+δμa(r),
D(r)=D0+δD(r),

where μ a0 and D 0 denote the homogenous background parts, and δμ a(r) and δD(r) the heterogenous parts. Based on this expansion, the DA can be solved using the Rytov approximation [7], where the photon fluence at position r due to a source at rs is given as

Urrs=U0rrseϕscrrs,

and we assume that (∇ϕ sc)2δμ a/D (i.e., the scattered field is slowly varying). In Eq. 4, U 0 represents the homogenous photon density corresponding to μ a0 and D 0, and ϕ sc the diffuse Rytov phase caused by δμ a(r) and δD(r). The task is then to find the relationship between δμ a(r)(δD) and ϕ sc(r,r s). After we discretize the whole domain into N equivolume voxels and assume N s source fibers and N d detectors, the forward model can be expressed as [7]

ϕscc=AcΔ,

where the superscript “c” denotes complex values, Δ ∈ R2N the change of optical parameters in each voxel, ϕ c sc ∈ CNsd the Rytov phase, Ac ∈ CNsd×N the weighting matrix, and N sd = N s × N d the number of source-detector pairs. More specifically,

ϕscc=[ϕscrs1rd1,...,ϕscrs1rdk,...,ϕscrsNsrdNd]T,
Δ=[δμa(r1),...,δμa(rj),...,δμa(rN)Δa,δD(r1),...,δD(rj),...,δD(rN)Δs]T,
Ac=[A111caA11NcaA111csA11NcsAik1caAikNcaAik1csAikNcsANsNd1caANsNdNcaANsNd1csANsNdNcs]=[AcaAcs].

The elements

Aikjca=U0rsirjG(rjrdk)h3D0U0rsirdk,
Aikjcs=U0rsirj.G(rjrdk)h3D0U0rsirdk

represent the effect of the i th source on the j th absorbing (scattering) voxel measured at the k th detector, where G(∙) denotes the Green function and h 3 the volume of each voxel. Note that (i) Eq. (5) can be used for different simple geometries such as infinite, semi-infinite, or slab, by modifying the elements in A c using the method of image sources [23] and extrapolated zero boundary conditions [24, 25]; and (ii) Eq. (6) and (8) can be easily extended to consider multiple frequencies.

Let ϕ sc = [R ϕ cs cT I ϕ cs cT]T, A = [RA cT IA cT]T, where R(∙) takes the real part of the matrix element and I(∙) the imaginary part; then (5) can be rewritten as

ϕsc=AΔ,

with ϕ ∈ RNtot, A ∈ RNtot×N, and N tot = 2×N sd. We will use this expression as our forward model.

2.2. Measurement model

We consider additive Gaussian noise and formulate the measurement model as

y=ϕsc+e,

where y ∈ RNtot denotes the measurement vector and e ∈ RNtot the noise vector. In DOT the optical signals are most often corrupted by shot noise, which follows a Poisson distribution. However, when the number of photons is very large, a Poisson distribution can be well approximated by a Gaussian distribution [26]. Since the noise variance is expected to be proportional to the number of photons at the detector, we choose the noise covariance matrix as

=σ2[0000],

where

0=[ϕ11c00ϕNdNsc],

“|∙|” denotes the module of a complex number, and σ represents the relaxing coefficient. In (13) and (14), we assumed that the noise is spatially uncorrelated for simplicity.

3. Tumor localization using LCMV beamforming

The inverse problem in DOT is to reconstruct the distributions of δμ a or μ ś from the measurement y. It is usually ill-posed, and different regularization techniques [9] and Bayesian approaches [10, 11] have been proposed to ameliorate this problem. However, in most cases, solving the inverse problem is associated with some Newton-type optimization procedure (e.g., the Levenberg-Marquardt algorithm) or regularized iterative approaches that induce the forward solver at each iteration and require knowledge of the gradient. Therefore, the whole procedure is computationally expensive, especially in 3D.

In this section we propose solving the inverse problem from a spatial filter point of view. Our method is simple and fast to implement. The key ideas are described as follows. For simplicity, we consider localizing only the abnormality in μ a. Namely, we let A c = A ca and Δ = Δa. From (5), we can see that ϕ sc is actually a superposition of the Rytov phases resulting from many voxels. More specifically, the l th column of matrix A reflects the contribution of the l th voxel in Δ on ϕ sc. Since Δ contains only the change of μ a compared with the background, then, assuming a localized abnormality, it is clear that ideally only a limited number of elements in Δ are nonzero, whose locations correspond to the abnormality’s position. Therefore, the inverse solution is essentially to find the elements in Δ that provide the greatest contribution to the measurements.

Suppose ϕ sc is composed of the fluence due to L(non-zero valued) voxels; then (12) can be equivalently expressed as

y=l=1LHind{l}δμa(rind{l})+e,

where ind{} denotes the set of indices for which Δ is nonzero and H ind{l} the ind{l} column of A. Based on this idea, we can think of solving the inverse problem as finding a spatial filter W ∈ RNtot for each voxel location such that the filter output reflects the contribution of that voxel to the output. Mathematically, we are now looking for a vector W so that by scanning the whole domain voxel by voxel, the filter output

ϕ̂i=WiTy,i=1,2,...,Ntot

would reflect the effect of the i th element of Δ. For the ideal narrowband spatial filter, it should satisfy

WiTHi={1,iind{}0,otherwise

In this case, it is easy to derive that when noise is absent, the filter output is the absorption coefficients at the abnormality location.

However, it is generally impossible to have completely zero output in the stop band. Therefore, we need to have a criterion such that the filter is optimal in that sense. The LCMV method provides such a guideline. The idea is to find a W such that the variance of the filter output ϕ̂ is minimized while satisfying the linear constraint (17). Physically, that is equivalent to minimizing the contribution to the filter output due to the signals in the stop band. Since the variance of a random signal is a measure of the strength, and for a vector case, it is defined as the sum of the variance of each component, the LCMV beamforming is mathematically expressed so as to find a filter Wi, i = 1,2,...,N tot, for each voxel in the medium such that

minWitrC(ϕ̂)subjecttoWiTHi=1,

where “tr(∙)” denotes the trace of a matrix and C(∙) the covariance matrix. Using the method of Lagrange multipliers, we obtain the solution to W as

WiT=[HiTC1(y)Hi]1HiTC1(y).

A more detailed derivation can be found in the Appendix and [17].

Two comments should be noted. First, unlike most other image reconstruction approaches, Eq. (19) requires estimating the covariance matrix C(y). In practice, C(y) is generally unknown and needs to be estimated from the measurements. Assuming we have M temporal observations, y 1,y 2,...,yM, each of them following Eq. (12), one common estimator is the unbiased sample covariance matrix

Ĉ(y)=1M1m=1M(ymy-)(ymy-)T,

where ym represents the m th temporal sample, and y¯ = ∑M m=1 ym/M the sample mean. The validity of this estimation is based on two conditions: (i) we assume that the underlying μ a distribution is constant for all measurements ym, and (ii) the number of observations M must be larger than the dimension of the measurement vector ym, i.e., N tot. The second condition is for Ĉ to be nonsingular and typically, M is chosen to be a few times larger than N tot. When M is smaller than or equal to N tot, the sample covariance matrix (20) would not be a good approximation of the true covariance matrix. In this case, a novel shrinkage covariance estimator can be applied, which is well conditioned and always positive definite [27]. Since how to estimate the covariance matrix is not the main idea of beamforming (although it is an important factor), we just use (20) as our estimator for simplicity.

 figure: Fig. 1.

Fig. 1. Simulation setup. The domain is composed of a 8×8×6 cm3 cube, and two possible spherical abnormalities with radius 1 cm are inserted. Twenty-five sources (circles) are placed on the bottom surface (z = 0 cm) and twenty-five detectors (black dots) on the top (z = 6 cm). The range along the x and y axes is [-4,4] cm.

Download Full Size | PDF

Our second comment is that unlike other inverse methods that usually reconstruct the whole vector Δ at one time, our method regards each voxel in the domain as a possible abnormality and attempts to estimate each element’s contribution to the measurement by performing a voxel-based scanning. This property favors computation, since the number of voxels N will not impose a restriction for solving the inverse problem as in most methods. More specifically, looking at Eq. (19), we can see that the major computational load comes from calculating C -1(y), which has a computation complexity O(N 3 tot). Compared with typical Newton methods, which have either an O(N 3) complexity if using direct matrix inversion or an O(N 2) complexity for each iteration if some iterative approaches (e.g., steepest descent) are employed, our method offers faster implementation.

4. Numerical examples

In this section, we provide numerical examples to validate the applicability of LCMV beam-forming for localizing abnormalities in human tissue. We also analyze its localization performance in terms of resolution, robustness to the absorption perturbation level, and robustness to the noise level. We consider a 3D transmission geometry as shown in Fig. 1, where the origin is at the center of the bottom surface (z = 0 cm), and the sides are of length 8 cm, 8 cm, and 6 cm along the x, y, and z directions, respectively. We place 25 sources on the bottom surface with 1.75 cm interdistance and 25 detectors on the top surface (z = 6 cm) with 1.5 cm interdistance. We choose a source modulation frequency of 200MHz and wavelength λ = 750nm

4.1. Tumor localization results

We first demonstrate the localization results using the LCMV beamforming. The background optical parameters are selected to be μ a0 = 0.05 cm-1, μ s0́ = 9.5 cm-1, and c = 22 cm/ns [28]. We consider the following four cases for illustration:

 figure: Fig. 2.

Fig. 2. Original (a) and reconstructed (b) δμ a distributions in cm-1, assuming only one absorbing abnormality. Small images show the cross-section layers at different z values with a 0.5 cm distance.

Download Full Size | PDF

 figure: Fig. 3.

Fig. 3. Original (a) and reconstructed (b) δμ ś distributions in cm-1, assuming only one scattering abnormality. Small images show the cross-section layers at different z values with a 0.5 cm distance.

Download Full Size | PDF

  • (A) There is only one spherical absorbing abnormality centered at [-1.5,1.25,2.9] cm with radius R = 1 cm and δμ a(r) = 0.2 cm-1 (see 2(a));
  • (B) There is only one spherical scattering abnormality centered at [2, -1.25,1.75] cm with R = 1 cm and δμ ś(r) = -4 cm-1 (see Fig. 3(a));
  • (C) There is only one spherical abnormality that is both absorptive and scattering. It is centered at [-1.5,1.25,2.9] cm with δμ a(r) = 0.2 cm-1 and δμ ś(r) = 2 cm-1 (see Fig. 4 (a));
  • (D) Both abnormalities specified in (A) and (B) are shown in the domain (see Fig. 5(a)).

For each case, we set up the forward model (5) using the method of image sources and assuming extrapolated zero boundary conditions. We divide the domain into small voxels with size 4×4×5 mm3, leading to a total number of 4,800 voxels. We add a Gaussian noise with σ = 0.01 to obtain the measurements.

The reconstructed images using the LCMV beamforming method are shown in Fig. 2(b), Fig. 3(b), Fig. 4(b), and Fig. 5(b), respectively. We make the following observations:

First, when there is only one abnormality in the medium, whether it be absorptive, scattering, or both, LCMV beamforming can recover its location fairly well. A peak is clearly shown at the abnormality’ position; see Fig. 2(b), Fig. 3(b) and Fig. 4(b). However, note that the filter output does not recover the real optical values but rather an index of where the problem is.

Second, when both the absorber and scatterer are present in the medium, the abnormalities can also be recovered. However, both parameters are somehow interrelated to each other, as we can see two peaks in the reconstructed images (Fig. 5(b)). In particular, since the absolute change of μ ś (δμ ś ≈ 4–5 cm-1) is much bigger than the change of μ aμ a ≈ 0.1–0.2 cm-1), which is also the situation in reality, we can observe a clear peak shown in the recovered absorption images at the scatterer’s location. Its intensity is far above the recovered absorption coefficient so that the latter can hardly be seen. One way to overcome this problem and increase the resolution of the recovered δμ a is to apply a filter that removes the effect of δμ ś. This filter can be readily designed based on the information obtained from the recovered δμ ś. The filtered result is shown in Fig. 6, where we can see the absorber’s location more clearly.

 figure: Fig. 4.

Fig. 4. Original (a) and reconstructed (b) δμ a and δμ ś distributions for one abnormality that is both absorptive and scattering. Small images show the cross-section layers of each parameter at different z values at 0.5 cm intervals.

Download Full Size | PDF

These examples show that our method can effectively localize the abnormalities in tissue. The biggest advantage of our method is its simplicity of implementation. After estimating the sample covariance matrix, the inverse solution itself takes only 112 seconds (using Matlab 7 on a PC with Pentium 4 2.6GHz CPU and 1G of RAM), which is efficient for 3D cases.

This method does not recover the exact values of δμ a and δμ ś, but it can be useful as an effective preprocessing tool. For example, the abnormality can first be localized preliminarily using LCMV beamforming, and then the a priori information obtained can be used to initiate starting points for more accurate (and also more time-consuming) methods. We will illustrate this application in the next subsection.

 figure: Fig. 5.

Fig. 5. Original (a) and reconstructed (b) δμ a and δμ ś distributions, assuming two abnormalities, one of which is absorbing and the other scattering. Small images show the cross-section layers at different z values at 0.5 cm intervals.

Download Full Size | PDF

 figure: Fig. 6.

Fig. 6. Reconstructed δμ a and δμ ś distributions after filtering out the effect of δμ ś in the recovered δμ a images. Small images show the cross-section layers of each parameter at different z values at 0.5 cm intervals.

Download Full Size | PDF

 figure: Fig. 7.

Fig. 7. Reconstructed μ a distributions using SIRT with (a) and without (b) the prior information provided by the LCMV solutions. The true abnormality location is shown in Fig. 2(a). Small images show the cross-section layers of each parameter at different z values at 0.5 cm intervals.

Download Full Size | PDF

4.2. Application of the LCMV beamforming as a preprocessing tool

In this subsection, we illustrate the application of the LCMV beamforming as a preprocessing tool for other reconstruction methods. For simplicity, we consider reconstructing the absorption coefficient only. The original μ a distribution is given in Fig. 2(a).

In Fig. 7(a), we show the reconstruction results using the simultaneous iterative reconstruction technique (SIRT), where the zero vector is used as the initial searching value (i.e., no prior information is applied). As we can see, although a peak can be seen at the abnormality’s location, the contrast is not very pronounced. In addition, the reconstructed δμ a value is not correct as well, as it is very biased toward the initial searching vector. We then perform the reconstruction by initializing the Δ vector based on the information obtained from the LCVM beamforming results (Fig. 2(b)). More specifically, we set the δμ a value to 0.2 at the positions where the LCMV filter output is more than half the peak value. The inverse results are shown in Fig. 7(b). We can see that the reconstruction result is greatly improved with a more localized abnormality.

4.3. Performance analysis

In this subsection, we analyze the localization performance of LCMV beamforming under different tumor conditions. For simplicity, we consider only absorptive abnormalities. We first study the effects of the absorption perturbation level and noise level on the localization results. Then we present preliminary results on the resolution of the beamforming method.

4.3.1. Effect of the absorption perturbation level

We assume one absorbing abnormality whose shape and position were specified in the case (A) of Subsection 4.1; see Fig. 2(a). The reconstructed result for the perturbation level δμ a = 0.2 cm-1 was shown in Fig. 2(b). According to [7], the perturbation in the tumor’s absorption value ranges from 0.02 cm-1 to 0.3 cm-1 for most biological tissues. Therefore, we reduce δμ a to 0.1 cm-1 and 0.05 cm-1 and applied our LCMV beamforming method. The results are shown in Fig. 8. Comparing Fig. 2(b) and Fig. 8, we observe that as the perturbation level decreases, the filter output signal becomes smaller as well. However, the abnormality can be reconstructed correctly with a pronounced contrast in all three cases, indicating the robustness of our method to small absorption perturbations.

 figure: Fig. 8.

Fig. 8. Reconstructed δμ a distributions in cm-1, assuming only one absorbing abnormality as shown in Fig. 2(a). (a) For true δμ a = 0.05 cm-1; (b) For true δμ a = 0.1 cm-1. Small images show the cross-section layers at different z values with a 0.5 cm distance.

Download Full Size | PDF

 figure: Fig. 9.

Fig. 9. Reconstructed δμ a distributions in cm-1 with different noise levels, assuming only one absorbing abnormality as shown in Fig. 2(a). (a) For σ = 0.1; (b) For σ = 1. Small images show the cross-section layers at different z values with a 0.5 cm distance.

Download Full Size | PDF

4.3.2. Effect of noise level

We assume the same setup as in Subsection 4.3.1 with δμ a = 0.2 cm-1. We increase the noise level to σ = 0.1 and σ = 1, and the localization results are shown in Fig. 9. Comparing with Fig. 2(b) where σ = 0.01, we observe that as the noise level increases, the size of the recovered abnormality increases as well (i.e., the resolution worsens). Nevertheless, the center of the abnormality can still be recovered correctly and the tumor is detected with clear contrast.

4.3.3. Resolution of the LCMV beamforming in DOT

We present two examples to illustrate the ability of LCMV beamforming to resolve two closely spaced abnormalities in DOT. We assume two spherical absorptive abnormalities with radius R = 0.75 cm and δμ a = 0.2 cm-1. In the first case, the centers of these two spheres are at [-1.5, 0.8, 2] cm and [1.5, -0.8, 4] cm, with a distance of 4 cm between two centers and 2.5 cm between the two closest points on the spheres. In the second case, the centers are moved to [-1, 0.5, 2.5] cm and [1, -0.5, 3.5] cm, and the distance becomes 2.45 cm between centers and 1 cm between the closet points. The localization results are shown in Fig. 10. We see that the abnormalities can be recovered with limited resolution when the two abnormalities are more than 2 cm apart (Fig. 10(a)). The contrast of the reconstructed image is not as good as for the single tumor case, but the two abnormalities can still be separated and the centers be recovered at the right places. When the two tumors are less than 1.5 cm apart, our LCMV beamforming can not differentiate them. Only one big sphere is shown in the reconstructed image (Fig. 10(b)), with no indication of whether this is from one big tumor or two small separated ones.

 figure: Fig. 10.

Fig. 10. Reconstructed δμ a distributions in cm-1, assuming two spherical absorbing abnormalities with R = 0.75 cm and δμ a = 0.2 cm-1. The centers of the two spheres are at (a) [- 1.5, 0.8, 2] cm and [1.5, -0.8, 4] cm; and (b) [-1, 0.5, 2.5] cm and [1, -0.5, 3.5] cm. Small images show the cross-section layers at different z values with a 0.5 cm distance.

Download Full Size | PDF

These examples show that the current LCMV beamforming algorithm has limited abilities to separate two closely spaced tumors using DOT. In fact, this limitation of beamformer is not restricted to DOT. It also exists in areas such as EEG/MEG source localization, radar and sonar signal processing, and is related to factors such as the signal-to-noise ratio, the number of source-detector pairs, and the filter design. Different approaches have been proposed to improve this issue. For example, Sekihara et al. [16] has proposed a multi-resolution scheme for MEG source imaging. An adaptive beamforming was also proposed in [29] for acoustic imaging.

5. Conclusion

We proposed a tumor localization scheme for diffuse optical tomography by using linearly constrained minimum variance beamforming. Our method is based on a voxel-by-voxel scanning approach, and thus it is fast to implement. It can be used alone to localize abnormailities that are absorbing, scattering, or both. It can also be used as an effective preprocessing method to improve the localization performances of other reconstruction approaches.

We analyzed the performance of our method under different tumor conditions. We found that the LCMV beamforming is fairly robust to small absorption perturbation levels. As the noise level increases, our beamformer can detect the location of the abnormality center correctly, whereas the resolution becomes worse. If there are two closely spaced abnormalities in the domain, LCMV beamforming fails to separate them when they are less than 1.5 cm apart. This last issue can be further improved by using the multiresolution approach [16] or adaptive beamforming [29]. Our future work includes implementing our method to real optical data and finding schemes to improve its resolution performance.

Appendix

We give a brief derivation in this Appendix for Eq. (18) using the method of Lagrange multipliers. We consider a more general case where Wi and Hi can be matrices instead of a row (column) vector as required by (18). For this reason, we omit the subscript i in the following derivation.

Let 2D be a matrix of Lagrange multipliers; then finding the matrix W satisfying (18) is equivalent to finding a W that minimizes the function

LWD=tr{WTCW+2(WTHI)D}.

Since trM = trMT for any square matrix M, we have

LWD=tr{WTCW+(WTHI)D+DT(HTWI)}.

After some mathematical manipulation,

LWD=tr{(WT+DTHTC1)C(W+C1HL)LLTLTHTC1HL}.

In (23), since only the first term in the bracket contains W and the covariance matrix C is positive definite, the minimum of L(W,D) is obtained by setting the first term to zero; i.e.,

W=C1HD.

Then the Lagrange multiplier matrix D is obtained by substituting (24) into the constraint WTH = I, and we have

DTHTC1H=I,

or

DT=(HTC1H)1.

Finally, substituting (26) back into (24), we obtain

WT=[HTC1H]1HTC1.

Acknowledgments

The authors thank Mr. Patricio S. La Rosa and Dr. Carlos M. Muravchik for their helpful comments.

References and links

1. A. H. Barnett, J. P. Culver, A. G. Sorensen, A. M. Dale, and D. A. Boas, “Bayesian estimation of optical properties of the human head via 3D structural MRI,” in Photon Migration and Diffuse-Light Imaging, D. A. Boas, ed. Proc. SPIE 5138,139–147 (2003). [CrossRef]  

2. G. Strangman, D. Boas, and J. Sutton, “Non-invasive neuroimaging using near-infrared light,” Biol. Psychiatry 52,679–693 (2002). [CrossRef]   [PubMed]  

3. A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. 20,435–442 (1997). [CrossRef]   [PubMed]  

4. D. Grosenick, T. Moesta, H. Wabnitz, J. Mucke, C. Stroszcynski, R. Macdonald, P. Schlag, and H. Rinnerberg, “Time-domain optical mammography: Initial clinial results on detection and characterization of breast tumors,” Appl. Opt. 42,3170–3186 (2003). [CrossRef]   [PubMed]  

5. X. Intes, J. Ripoll, Y. Chen, S. Nioka, A. Yodh, and B. Chance, “In vivo continuous-wave optical breast imaging enhanced with Indocyanine Green,” Med. Phys. 30,1039–1047 (2003). [CrossRef]   [PubMed]  

6. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15,R41–R93 (1999). [CrossRef]  

7. M. A. O’Leary, “Imaging with diffuse photon density waves,” Ph.D. thesis, University of Pennsylvania (1996).

8. R. J. Gaudette, D. H. Brooks, C. A. DiMarzio, M. E. Kilmer, E. L. Miller, T. Gaudette, and D. A. Boas, “A comparison study of linear reconstruction techniques for diffuse optical tomographic imaging of absorption coefficient,” Phys. Med. Biol. 45,1051–1070 (2000). [CrossRef]   [PubMed]  

9. B. Pogue, T. McBride, J. Prewitt, U. Osterberg, and K. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38,2950–2961 (1999). [CrossRef]  

10. M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Phys. Med. Biol. 50,2837–2858 (2005). [CrossRef]   [PubMed]  

11. J. C. Ye, C. A. Bouman, K. J. Webb, and R. P. Millane, “Nonlinear multigrid algorithms for Bayesian optical diffuse tomography,” IEEE Trans. Image Process. 10,909–922 (2001). [CrossRef]  

12. D. Grosenick, H. Wabnitz, 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). [CrossRef]  

13. D. Grosenick, H. Wabnitz, K. T. Moesta, J. Mucke, M. Möller, C. Stroszczynski, J. Stöbel, B. Wassermann, P. M. Schlag, and H. Rinneberg, “Concentration and oxygen saturation of haemoglobin of 50 breast tumours determined by time-domain optical mammography,” Phys. Med. Biol. 49,1165–1181, (2004). [CrossRef]   [PubMed]  

14. B. D. Van Veen and K. M. Buckley, “Beamforming: A versatile approach to spatial filtering,” IEEE ASSP. Magazine 5,4–24 (1988). [CrossRef]  

15. K. Sekihara and S. S. Nagarajan, “Neuromagnetic source reconstruction and inverse modeling,” in Proceedings of IEEE EMBS Asian-Pacific Conference on Biomedical Engineering (Institute of Electrical and Electronics Engineers, Keihanna, Japan, 2003), pp.20–22.

16. K. Sekihara, S. S. Nagarajan, D. Poeppel, A. Marantz, and Y. Miyashita, “Reconstructing spatio-temporal activities of neural sources using an MEG vector beamforming technique,” IEEE Trans. Biomed. Eng. 48,760–771 (2001). [CrossRef]   [PubMed]  

17. B. D. Van Veen, W. van Drongelen, M. Yuchtman, and A. Suzuki, “Localization of brain electrial activity via linearly constrained minimum variance spatial filtering,” IEEE Trans. Biomed. Eng. 44,867–880 (1997). [CrossRef]   [PubMed]  

18. M. E. Spencer, R. M. Leahy, J. C. Mosher, and P. S. Lewis, “Adaptive filters for monitoring localized brain activity from surface potential time series,” in Conference record of the twenty-sixth Asilomar conference on Signals, Systems and Computers (Institute of Electrical and Electronics Engineers, Pacific Grove, CA, 1992), pp.156–161.

19. J.-F. Synnevag, A. Austeng, and S. Holm, “Minimum variance adaptive beamforming applied to medical ultrasound imaging,” in Proceedings of IEEE Ultrasonics Symposium, (Institute of Electrical and Electronics Engineers, Rotterdam, Netherlands, 2005), pp.1199–1202.

20. A. C. Kak and M. Slaney, Principles of computerized tomographic imaging (IEEE Press, 1988).

21. D. A. Boas, “A fundamental limitation of linearized algorithms for diffuse optical tomography,” Opt. Express 1,404–413 (1997). [CrossRef]   [PubMed]  

22. V. Ntziachristos, B. Chance, and A. G. Yodh, “Differential diffuse optical tomography,” Opt. Express 5,230–242 (1999). [CrossRef]   [PubMed]  

23. M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties,” J. Appl. Opt. 28,2331–2336 (1989). [CrossRef]  

24. R. Aronson, “Boundary conditions for diffusion of light,” J. Opt. Soc. Am. A 12,2532–2539 (1995). [CrossRef]  

25. R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11,2727–2741 (1994). [CrossRef]  

26. S. O. Rice, “Mathematical analysis of random noise,” Bell Syst. Tech. J. 23,282–332 (1944).

27. J. Schäfer and K. Strimmer, “A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics,” Statist. Appl. Genet. Mol. Biol. 4, Article 32 (2005).

28. M. J. Holboke, B. J. Tromberg, X. Li, N. Shah, J. Fishkin, D. Kidney, J. Butler, B. Chance, and A.G. Yodh, “Three-dimensional diffuse optical mammography with ultrasound localization in a human subject,” J. Biomed. Opt. 5,237–247 (2000). [CrossRef]   [PubMed]  

29. M. Papazoglou and J. L. Krolik, “High resolution adaptive beamforming for three-dimensional acoustic imaging of zooplankton,” J. Acoust. Soc. Am. 100,3621–3630 (1996). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (10)

Fig. 1.
Fig. 1. Simulation setup. The domain is composed of a 8×8×6 cm3 cube, and two possible spherical abnormalities with radius 1 cm are inserted. Twenty-five sources (circles) are placed on the bottom surface (z = 0 cm) and twenty-five detectors (black dots) on the top (z = 6 cm). The range along the x and y axes is [-4,4] cm.
Fig. 2.
Fig. 2. Original (a) and reconstructed (b) δμ a distributions in cm-1, assuming only one absorbing abnormality. Small images show the cross-section layers at different z values with a 0.5 cm distance.
Fig. 3.
Fig. 3. Original (a) and reconstructed (b) δμ ś distributions in cm-1, assuming only one scattering abnormality. Small images show the cross-section layers at different z values with a 0.5 cm distance.
Fig. 4.
Fig. 4. Original (a) and reconstructed (b) δμ a and δμ ś distributions for one abnormality that is both absorptive and scattering. Small images show the cross-section layers of each parameter at different z values at 0.5 cm intervals.
Fig. 5.
Fig. 5. Original (a) and reconstructed (b) δμ a and δμ ś distributions, assuming two abnormalities, one of which is absorbing and the other scattering. Small images show the cross-section layers at different z values at 0.5 cm intervals.
Fig. 6.
Fig. 6. Reconstructed δμ a and δμ ś distributions after filtering out the effect of δμ ś in the recovered δμ a images. Small images show the cross-section layers of each parameter at different z values at 0.5 cm intervals.
Fig. 7.
Fig. 7. Reconstructed μ a distributions using SIRT with (a) and without (b) the prior information provided by the LCMV solutions. The true abnormality location is shown in Fig. 2(a). Small images show the cross-section layers of each parameter at different z values at 0.5 cm intervals.
Fig. 8.
Fig. 8. Reconstructed δμ a distributions in cm-1, assuming only one absorbing abnormality as shown in Fig. 2(a). (a) For true δμ a = 0.05 cm-1; (b) For true δμ a = 0.1 cm-1. Small images show the cross-section layers at different z values with a 0.5 cm distance.
Fig. 9.
Fig. 9. Reconstructed δμ a distributions in cm-1 with different noise levels, assuming only one absorbing abnormality as shown in Fig. 2(a). (a) For σ = 0.1; (b) For σ = 1. Small images show the cross-section layers at different z values with a 0.5 cm distance.
Fig. 10.
Fig. 10. Reconstructed δμ a distributions in cm-1, assuming two spherical absorbing abnormalities with R = 0.75 cm and δμ a = 0.2 cm-1. The centers of the two spheres are at (a) [- 1.5, 0.8, 2] cm and [1.5, -0.8, 4] cm; and (b) [-1, 0.5, 2.5] cm and [1, -0.5, 3.5] cm. Small images show the cross-section layers at different z values with a 0.5 cm distance.

Equations (27)

Equations on this page are rendered with MathJax. Learn more.

. D ( r ) U r ω + c μ a ( r ) U r ω + jωU r ω = c q 0 r ω ,
μ a ( r ) = μ a 0 + δ μ a ( r ) ,
D ( r ) = D 0 + δD ( r ) ,
U r r s = U 0 r r s e ϕ sc r r s ,
ϕ sc c = A c Δ ,
ϕ sc c = [ ϕ sc r s 1 r d 1 , . . . , ϕ sc r s 1 r d k , . . . , ϕ sc r s N s r d N d ] T ,
Δ = [ δ μ a ( r 1 ) , . . . , δ μ a ( r j ) , . . . , δ μ a ( r N ) Δ a , δ D ( r 1 ) , . . . , δD ( r j ) , . . . , δ D ( r N ) Δ s ] T ,
A c = [ A 111 ca A 11 N ca A 111 cs A 11 N cs A ik 1 ca A ikN ca A ik 1 cs A ikN cs A N s N d 1 ca A N s N d N ca A N s N d 1 cs A N s N d N cs ] = [ A ca A cs ] .
A ikj ca = U 0 r s i r j G ( r j r d k ) h 3 D 0 U 0 r s i r d k ,
A ikj cs = U 0 r s i r j . G ( r j r d k ) h 3 D 0 U 0 r s i r d k
ϕ sc = A Δ ,
y = ϕ sc + e ,
= σ 2 [ 0 0 0 0 ] ,
0 = [ ϕ 11 c 0 0 ϕ N d N s c ] ,
y = l = 1 L H ind { l } δ μ a ( r ind { l } ) + e ,
ϕ ̂ i = W i T y , i = 1,2 , . . . , N tot
W i T H i = { 1 , i ind { } 0 , otherwise
min W i tr C ( ϕ ̂ ) subject to W i T H i = 1 ,
W i T = [ H i T C 1 ( y ) H i ] 1 H i T C 1 ( y ) .
C ̂ ( y ) = 1 M 1 m = 1 M ( y m y - ) ( y m y - ) T ,
L W D = tr { W T CW + 2 ( W T H I ) D } .
L W D = tr { W T CW + ( W T H I ) D + D T ( H T W I ) } .
L W D = tr { ( W T + D T H T C 1 ) C ( W + C 1 HL ) L L T L T H T C 1 HL } .
W = C 1 HD .
D T H T C 1 H = I ,
D T = ( H T C 1 H ) 1 .
W T = [ H T C 1 H ] 1 H T C 1 .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.