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

Abstract: We present a tumor localization method for diffuse optical tomography using linearly constrained minimum variance (LCMV) beamforming. 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.


Introduction
Diffuse optical tomography (DOT) is becoming a useful complement to the current noninvasive 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.

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.

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] where U is the fluence (W/m 2 ), μ a the absorption coefficient (m −1 ), D the diffusion coefficient (m 2 /s), c the speed of light in the medium (m/s), and q 0 the source (W/m 3 ).The term D is defined as D = c/3(μ a + μ s ), where μ s = (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 μ s for most biological tissues [6], we assume D = 1/3μ s 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 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 r s is given as 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] φ where the superscript "c" denotes complex values, Δ ∈ R 2N the change of optical parameters in each voxel, φ c sc ∈ C N sd the Rytov phase, A c ∈ C N sd ×N the weighting matrix, and N sd = N s × N d the number of source-detector pairs.More specifically, φ c sc = [φ sc (r s1 , r d1 ),... ,φ sc (r s1 , r dk ),... ,φ sc (r sN s , r dN d )] T , (6) Δ = [δ μ a (r 1 ),... ,δ μ a (r j ),... ,δ μ a (r N ) Δ a , δ D(r 1 ),... ,δ D(r j ),... ,δ D(r N ) The elements 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φ cT sc I φ cT sc ] T , A = [RA cT I A 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Δ, (11) with φ ∈ R N tot , A ∈ R N tot ×N , and N tot = 2 × N sd .We will use this expression as our forward model.

Measurement model
We consider additive Gaussian noise and formulate the measurement model as where y ∈ R N tot denotes the measurement vector and e ∈ R N tot 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 where "| • |" 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.

Tumor localization using LCMV beamforming
The inverse problem in DOT is to reconstruct the distributions of δ μ a or δ μ s 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 where ind{} denotes the set of indices for which Δ is nonzero and H ind{l} the ind{l} th column of A. Based on this idea, we can think of solving the inverse problem as finding a spatial filter W ∈ R N tot 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 would reflect the effect of the i th element of Δ.For the ideal narrowband spatial filter, it should satisfy 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 W i , i = 1, 2,...,N tot , for each voxel in the medium such that min 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 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 ,... ,y M , each of them following Eq.( 12), one common estimator is the unbiased sample covariance matrix where y m represents the m th temporal sample, and ȳ = ∑ M m=1 y m /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 y m , and (ii) the number of observations M must be larger than the dimension of the measurement vector y m , 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.
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 voxelbased 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.

Numerical examples
In this section, we provide numerical examples to validate the applicability of LCMV beamforming 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.

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:  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 mm 3 , 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 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 μ s (δ μ s ≈ 4-5 cm −1 ) is much bigger than the change of μ a (δ μ a ≈ 0.1-0.2cm −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 δ μ s .This filter can be readily designed based on the information obtained from the recovered δ μ s .The filtered result is shown in Fig. 6, where we can see the absorber's location more clearly.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 δ μ s , 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.Fig. 6.Reconstructed δ μ a and δ μ s distributions after filtering out the effect of δ μ s in the recovered δ μ a images.Small images show the cross-section layers of each parameter at different z values at 0.5 cm intervals.

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.

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.

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

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]  beamforming [29].Our future work includes implementing our method to real optical data and finding schemes to improve its resolution performance.

Fig. 1 .
Fig. 1.Simulation setup.The domain is composed of a 8×8×6 cm 3 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) δ μ s 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 δ μ s 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 δ μ s 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. 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. 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.