1 Introduction

In recent years, the problem of forming bio-medical images from sparse observations of diffuse photon density wave data (DPDW) has received considerable attention, [1, 4, 5, 3]. A common goal of such processing is the detection, localization, and characterization of regions in the body, which we refer to as anomalies, that have optical properties (space varying absorption coefficient, µ_a (r), and scattering coefficient, µ_s (r)) that are different from those of the surrounding tissue. By extracting such information about anomalies from the DPDW data, one can potentially draw medically useful conclusions related to the flow of blood in an area, the presence of a tumor, etc.

A typical approach to these problems is the “image-then-detect” method, where the first step is to use measured DPDW data to form an image in two dimensions, or volumetric rendering in 3-D, of a subcutanous region of interest. This step is followed by image postprocessing via standard techniques to extract and analyze anomalous areas. Thus, a major difficulty with these methods is the need to first form a pixel by pixel reconstruction of the underlying region, as this step requires the solution to an inverse scattering problem. Since the inverse scattering problem itself is highly ill-posed in the sense that small changes in the data can have dramatic, adverse affects on the estimated solution, no matter which model (Born, Rytov, exact, etc.) is used to approximate the forward scattering problem, the computed solution will be hopelessly contaminated by noise artifacts. To lessen the sensitivity to noise, some type of regularization technique must be employed - any number of techniques is possible [7].

A further difficulty with standard approaches is that a small amount of measured DPDW data is used to recover the preliminary image comprised of a large number of pixels. Our approach to the imaging problem is fundamentally different in that we parameterize the problem directly in terms of a small number of parameters that describe the shape, location and contrast of the anomaly as well as the value of the background. As a result, we are able to localize anomalies directly without postprocessing. This formulation requires the minimization of a functional with respect to only a small number (relative to the total number of pixels reconstructed) of parameters, ensuring the efficiency of the method.

The idea of reducing the number of unknowns by appropriate parameterization is not new. The authors began to explore this approach for DPDW data in [9]. Another such method that has been applied to DPDW problems is presented in [2]. In that work, the authors use the diffusion equation approximation with piecewise constant diffusion and absorption to model the forward problem. They assume that the finite values of the diffusion and absorption are known a priori, and they seek information about the geometry of the smooth region boundaries where the functions are discontinuous. In other words, they try to solve for the few, unknown parameters that describe the boundaries of each of their curves. Our method is significantly different than the method in [2] for several reasons. We aim to reconstruct maps of the absorption within only two regions, the anomaly and the background. However, we do not need to assume the values of the absorption in different regions a priori: we obtain such information as we simultaneously reconstruct boundaries of the two regions. In the purest form, our model does not restrict the absorption perturbation within the two regions to be constant, although we make that simplification for the results presented here. The method of parameterization of the boundaries in our work and the work in [2] is different. Unlike the method in [2], our method incorporates regularization that helps combat the ill-conditioning and the underdetermined nature of the problem. Finally, the results we present here are for a different geometry than the one used in [2].

The remainder of this paper is organized as follows. In §2, the mathematical model for the diffuse photon density wave system of interest is presented. We describe our algorithm in §3 and provide numerical examples in §4. Finally, we discuss conclusions and future work in §5.

2 Models

2.1 Forward DPDW Model

In this paper, we consider information extraction algorithms based on the first-Born approximation to the frequency domain, integral equation [9, equation 1] formulation of the diffusion equation describing the flux of the optical modulation envelope. Specifically, for a point source located at position r_s , the model we have is

where V is the region where the anomaly is assumed to exist (see Fig. 1), G(r, r′) is the Green’s function for the problem at hand, and g represents the space varying perturbation in the absorption coefficient about a known deterministic background used in the specification of G. The number D is the diffusion coefficient and v is the propagation velocity in the medium. Our Green’s function is constructed using the extrapolated boundary condition of [6]. We assume that the transmitters can be modeled as point sources on the air-tissue interface. Unlike the work in [9], we consider here a transmission geometry, so the scattered diffuse photon wavefields are measured by an array of N_r point receivers located below the bottom air-tissue interface (refer to Figure 1).

In this work it is assumed that we have a total of N_s point sources where the scattered fields generated by each source are observed over an array of N_r receivers. We are using a 2-dimensional model. Thus, for the ith source, at a given modulation frequency, we have a data vector, y_i , of length 2N_r (for DC, the vector length is N_r ) comprised of the in-phase and quadrature components of the measured scattered wavefield. Using a method of moments procedure to discretize (1) with pulse basis functions to represent g(r) [8], we have the discrete equations

where G_i is the discretization of the Born kernel associated with the ith source, g is the vector of pulse basis expansion coefficients for g(r), and n i is the noise vector. Finally, all y _i are collected for the various modulation frequencies into a single vector, ${\mathbf{y}}^T=\left[{\mathbf{y}}_1^T\dots {\mathbf{y}}_{N_s}^T\right]$ to arrive at the discretized data model

with G and n obtained by “stacking” all the G _i and n _i for each frequency. In our experiments, we take data at 0 and 200 MHz, so y is of length 3N_rN_s .

 

Fig. 1. Source/receiver configuration in the transmission geometry.

Download Full Size | PDF

2.2 A Model for g(r)

The model we use here to describe g(r), the unknown perturbation in the absorption coefficient, was originally described by the authors in [9]. Here, g(r) is written as

where S(r) is an indicator function that is one over the (unknown) support of the object and zero elsewhere. The vectors a₁, a₂ hold the expansion coefficients that effectively determine the weight of each function. The 1×p vectors B_i i=1, 2 have as their components functions b_i,j (r), j=1, …, p which are the expansion functions. In other words, we model the variation of g over the anomaly as a linear combination of certain functions and the variation of g on the background as a different linear combination of functions. The particular choice of the b_i,j depends on the application. In this paper, for instance, we assume that there is a homogeneous anomaly of contrast a _1,1 against a homogeneous background of value a _2,1: thus, B ₁=b _1,1(r)=b _2,1(r)=B ₂=1. That is, we assume g is piecewise constant. Use of higher order polynomials, trigonometric functions etc. provides greater flexibility in capturing true, underlying inhomogeneities, but adds a small amount of computational complexity to the problem. In this paper, the objective is to determine the structure of S and the a _i given a data vector y related to g via (1) assuming the B_i are known.

In order to determine the support of the anomaly, and hence the S(r), we define the contour of the anomaly as a B-spline curve c(s):

for a given set of x̂_i , ŷ_i expansion coefficients, or control points, and a set of periodic, quadratic functions, ${\beta }_{k_i}\left(s\right)$. Each ${\beta }_{k_i}\left(s\right)$ has support over a subinterval [k_i , k_{i +3}]∊[0, L], and the sequence {k ₀,…, k_{K -1}} is called the knot sequence. Since the basis was taken to be periodic, [x̂ ₀, ŷ ₀]=[x̂ _K-2, ŷ _K-2] and [x̂ ₁, ŷ ₁]=[x̂ _K-1, ŷ _K-1]. According to our model, if (x, y) is a point outside the curve, S(r)=0, otherwise, S(r)=1.

From the point of view of discrete implementation of our model, S will be a diagonal matrix which functions as the discrete equivalent to S(r) in (4). If the center of the Ith pixel in the discretization is inside the closed curve c(s) then we set S _ll =1, otherwise we set S _ll =0. Define the m×p matrices B₁, B₂ to have columns given by B ₁(x_j ,y_i ), B ₂(x_j , y_i ) (with (x_j , y_i ) denoting the pixel center coordinates and m denoting the total number of pixels) ordering lexicographically first by increasing i then by increasing j. With these definitions the discrete version of (4) is

3 Algorithm Description

Our algorithm seeks a good approximation to g(r) by successively improving approximations to a₁, a₂, and S. Let c*(s) denote the true contour of the true anomaly, and let c ₀(s) denote an initial guess¹ to c*(s). We note that is highly likely that the number and locations of distinct control points of c*(s) and c ₀(s) will be different.

The c ₀(s) defines an initial guess at S, which we use to estimate a₁, a₂ as the minimizers of

Note that if the number of columns in B₁ and B₂ are small (as is the case for problems of interest here) compared to 3N_sN_r , the height of GQ, this problem is quickly and easily solved². The total cost associated with any curve c(s) (hence, with our estimate of the solution g(r)) therefore is

where the first term enforces fidelity to the data while the second denotes the regularizer. The regularizer enforces a priori information we have about the shape of the anomaly and serves also to stabilize the least squares solution. In this work, we chose Ω(c) to be

where (x̂_i , ŷ_i ) denotes the ith pair of control points for the current estimate, c(s), of c*(s). The idea is that large gaps between adjacent control points are penalized more, thereby discouraging the algorithm from choosing non-physical anomalies. For example, since there are neither sources nor receivers on the sides of the domain we are trying to reconstruct, without regularization we are more likely to reconstruct tall, narrow anomalies in presence of noisy data. Clearly a value of λ that is too large or too small will cause the reconstruction to be under or over dependent on the data, respectively. We do not address the selection of a near optimal parameter here, but refer the interested reader to the abundance of literature on the subject (see [7] and the references therein).

The algorithm proceeds by perturbing the curve c ₀(s), in a fashion to be described shortly, to obtain a finite collection of new curves, compute the associated cost (and a₁, a₂) for each new curve from (8), and then select the new estimate, c ₁(s), of c*(s) as a curve which gave the minimum cost over all the perturbations, provided that minimum is less than or equal to the current cost for c ₀(s). The curves are obtained by taking each of the K-2 distinct control points of c ₀(s) and moving them by a fixed amount h, in mm, in the vertical, horizontal, and diagonal directions. Thus, there are 8 possible moves for each of the K-2 distinct control points, resulting in at most 8(K-2) new curves for which the cost is computed. Then c ₁(s) is selected from among these curves and the process is iterated³ until the cost can no longer be reduced by making these types of perturbations to the control points of the current estimate c_k (s). The algorithm is outlined in Figure 2, where h is understood to have been chosen apriori.

 

Fig. 2. Anomaly Recovery Algorithm

Download Full Size | PDF

4 Numerical Experiments

All our experiments were conducted in Matlab using double precision, floating point arithmetic. We used Matlab’s Spline Toolbox to create and manipulate the B-splines. The region of interest is 3cm in width and 3cm in depth, and the region is discretized into a 31×31 pixel image. There were 10 sources across the top of the region and 10 receivers along the bottom, and data was collected at the modulation frequencies of 0 and 200 MHz. Thus, the matrix G had dimension 300×31².

In both examples, we assumed that the object was a homogeneous perturbation on a homogeneous background of unknown value. Hence, the matrices B₁ and B₂ in (6) were both taken to be m×1 matrices with entries of all ones. Therefore the two coefficients a₁ and a₂ were both estimated during the minimization of the cost function in (8). This is different from the work reported in [9] where the a₂ was taken to be zero both in forming the data and during the reconstruction phase. Recall that g(r) represents the space varying perturbation in the absorption coefficient, µ _a, relative to the background coefficient used to determine G. In our experiments, we used µ_a =.05/cm and µ′_s =10/cm to construct G.

To generate the data, c*(s) was taken to be a quadratic spline curve defined using 6 distinct control points and S* was the corresponding discrete indicator matrix. With a*₁,a*₂ as the true values of the anomaly and background, respectively, the noise-free data was formed as the matrix-vector product G _g with g in (6) and G described above. This generated a data vector of length 300, the first 100 points of which corresponded to 0 MHz, the second 100 to the in-phase at 200, and the 3rd to the quadrature at 200.

We added shot noise to the data as follows. For a particular source and receiver pair the noise component n_sr (ω) (at any given frequency ω, for either an in-phase or quadrature component) is given by

where v ₁ is a number chosen from a normal distribution with mean zero and variance one⁴. Here σ_sr (0) is proportional to the square root of the total signal strength at DC:

where we denote by ỹ_s (r) the noise-free scattered field corresponding to source s, receiver r, at DC and we use ${\stackrel{\mathit{~}}{y}}_s^{\mathit{\text{inc}}}$ (r) to denote the incident field at source s, receiver r at DC, and γ is the proportionality constant. We denote the 300-length, noise-contaminated data vector by y and the noise-free data by ỹ. The value of the signal to noise ratio that we give in the examples corresponds to

$\mathit{S}NR=10{\log }_{10}\frac{{\parallel \tilde{\mathbf{y}}\parallel }_2^2}{{\parallel \mathbf{n}\parallel }_2^2}.$.

Since the noise is not white, we must replace J(a₁, a₂) in (7) by

$J({\mathbf{a}}_1,{\mathbf{a}}_2)≔{\parallel {\Sigma }^{-1}(\mathbf{G}[{\mathbf{SB}}_1,\left(\mathbf{I}-\mathbf{S}\right){\mathbf{B}}_2]\left[\begin{array}{c}{\mathbf{a}}_1\\ {\mathbf{a}}_2\end{array}\right]-\mathbf{y})\parallel }_2,$,

where ∑^-1 is a diagonal matrix with 1/σ_sr (0) placed appropriately along the diagonal.

In our examples, we compare the solution computed using our algorithm with the solution obtained via the commonly used truncated singular value decomposition (TSVD) method [7]. Essentially, if G=${\sum }_{i=1}^n$ µ_i u _i ${\mathbf{v}}_i^T$ is the singular value decomposition of the matrix G, the TSVD regularized solution for truncation index j,j≤n is ${\mathbf{g}}^{ts\upsilon d}=\sum _{i=1}^j\frac{{\mathbf{u}}_i^T\mathbf{y}}{{\mu }_i}{\mathbf{v}}_i$. We refer to the best, or optimal, TSVD solution as the one corresponding to the index j which gives the smallest relative error in the 2-norm with respect to the true solution.

4.1 Example 1

The contour of the anomaly, b*(r), in this example is the solid green curve in Figure 3 a. We determine our discrete, noise-free data vector ỹ by using b*(r) to generate S, setting a₁=10, a₂=-1, forming g _true according to (6), and setting ỹ=Gg _true . Physically, these values correspond to a background perturbation in µ_a from the known value of .05/cm by -.01/cm and a perturbation in µ_a within the anomaly of .10/cm. The discretized true image is shown in Figure 3b. We then added a noise vector n to the data to generate the noisy data y=ỹ+n. The noise was determined as described above so that the SNR relative to the data was about 23dB (the proportionality constant γ was 1e-9).

Figure 3c shows the optimal TSVD solution using the noisy data. Notice that the TSVD reconstruction incorrectly shows the anomaly to be much larger than it should be, and that the computed value of the anomaly is less than 1/5 of what it should be. The relative 2-norm difference error between g^true and g^tsvd (that is, the relative root mean squared error (RMSE)) is ‖g ^true -g ^tsvd ‖2/‖g ^true ‖2=.84.

Since the TSVD solution was already available, we used this to initialize c ₀(s) for our algorithm. Thus, the blue curve in Figure 3a shows c ₀(s) to be an approximation to the outline of where the anomaly appears to be in the TSVD reconstruction, which we obtained by postprocessing the TSVD image manually. We set h, the amount to perturb the control points, to 1 mm in our algorithm and took λ=.005 in (9). The final estimate of the contour, b_k *(s), that is reached by the algorithm is the red curve in Figure 3b; that is, no other curve is locally more optimal than that curve. The total reconstruction obtained using our algorithm to determine both b_k *(s) and a₁,a₂ is displayed in Figure 3d. The relative 2-norm error between the true and our computed solution is .56, which is a significant improvement in the error over the TSVD solution. Notice that although the relative error value seems high, as evidenced by the Figure 3d the shape-based reconstruction very accurately reconstructs the value of g inside and outside the anomaly and fairly accurately estimates the extent of the anomaly (the reconstruction is off by only 10 pixels). Thus, the relative error measure does not accurately reflect the high quality of the reconstruction. Further, recall from the discussion at the beginning of this section that c*(s) was defined with 6 distinct control points while our reconstruction assumed only 5 distinct control points. We believe that if we use more control points for c ₀(s), we could have reconstructed the irregular contour more accurately. The trade-off for using more control points in the reconstruction is the time it takes to execute one outer iteration.

4.2 Example 2

For our second example, the contour of the anomaly is given by the blue curve in Figure 4a. The value of g(r) inside the anomaly, a₁ was 20 and the value of the perturbation on the background, a₂, was -2 (corresponding to perturbations in µ_a of .2/cm and -.02/cm, respectively). The true solution, illustrated in Figure 4b, and noisy data y were formed as described in example 1 with the noise computed for a 28 dB SNR (the proportionality constant γ is the same as in Example 1).

Figure 4c illustrates the TSVD solution. Notice that the value of the anomaly at best is 1/4 of what it should have been, and that the anomaly looks to be smeared from top to bottom, as in the first example. This latter phenomenon is to be expected due to the “averaging” nature of the TSVD method and the lack of data at the sides of the region of interest. The relative error between the true solution and TSVD solutions is .84.

Again, because the TSVD solution was available, we manually took an estimate of the outline of the reconstructed anomaly to serve as the initial guess at the contour, c ₀(s). The curve c ₀(s) is displayed as the blue curve in Figure 4a. As in the first example, we set the curve perturbation h to 1 mm. For this example, we took λ=.05.

The final estimate of the contour is shown in red in Figure 4a, while the reconstructed image is displayed in Figure 4d. Note the vast improvement of our reconstruction over that obtained by TSVD. We estimate a1 as 22.8 and a2 as -2. The relative error between the true solution and the estimate obtained by our algorithm was .48, almost 1/2 the TSVD error. Also, the shape of our reconstructed anomaly is only off by 8 pixels with respect to the shape of the true anomaly. Here again, we believe that if we use more control points for c ₀(s), we could have reconstructed the contour more accurately.

 

Fig. 3. From left to right, top to bottom: a) Contour curves for b*(r), c ₀(r), b_k *(r) b) True image of the perturbation in the absorption coefficient c) TSVD reconstruction of the absorption perturbation d) our reconstruction of the absorption using h=1 and λ=.005.

Download Full Size | PDF

 

Fig. 4. From left to right, top to bottom: a) Contour curves for b*(r), c ₀(r), b_k *(r) b) True image of the perturbation in the absorption coefficient c) TSVD reconstruction of the absorption perturbation d) our reconstruction of the absorption using h=1 and λ=.05.

Download Full Size | PDF

5 Conclusions and Future Work

In this work, we presented an effective technique for simultaneously solving the image formation and object characterization problems from DPDW data with a particular Gaussian noise model. The success of our method was based on the underlying model for the perturbation of the absorption coefficient. Our model formulates the problem in terms of only a small number of unknowns: namely, those that describe the B-spline contour of the anomaly and those that describe the values of the perturbations over both the anomaly and the background. Our examples illustrated that our reconstructions contain significantly more qualitative and quantitative information than the commonly used TSVD reconstructions. Further, although our algorithm was implemented in serial, it is in fact highly parallelizable in that the cost evaluations at any given iteration can be performed in parallel. Thus, it is computationally feasible to consider more complicated structures (i.e. those with more control points that would each need to be perturbed during an iteration). Finally, the model and algorithm presented here are readily adapted to account for non-linear scattering models that arise when the Born approximation is no longer valid. In the future we hope to report results for non-linear scattering models and to generalize this work for scattering problems and 3D reconstructions. The generalization of our approach to multiple anomalies is straightforward if initial estimates of boundaries of such anomalies can be found; more work needs to be done to generalize our method to find multiple anomalies with a single boundary as a starting guess.

Although we have used a fixed knot sequence and number of control points to describe our B-spline curve, we note that it is possible to add and delete knots (and correspondingly, control points) from our description of the contour as well as specify more complicated basis functions B ₁(r) and B ₂(r), thereby allowing reconstruction of more complex anomalies and backgrounds. The question of how to choose A appropriately is also an important issue. All these considerations were beyond the scope of this work and remain areas for future research.