Diffuse optical tomography (DOT) is an emerging technology for noninvasive imaging of large human organs. The diffusion equation, the basis for most DOT reconstruction algorithms,1, 2, 3, 4 can be analytically modeled as a linear system under the Born or Rytov approximation.5, 6, 7 Image reconstruction is often a severely ill-posed problem due to the finite number of linearly independent measurements, noise, and often complicated structure of the heterogeneities in the tissue. Reconstruction algorithms play a very important role in DOT, as the image formation processes are nontrivial. For this reason, the reconstruction problem has received the attention of many researchers. One of the conventional ways is to formulate an optimization problem that calculates the forward problem iteratively and minimizes the residue of the measured photon density perturbations. Adaptive mesh refinement algorithms for finite element forward solver8, 9, 10, 11, 12, 13, 14, 15 and conjugate gradient methods16, 17, 18, 19, 20, 21, 22, 23 are popular in this approach, and many modified and improved algorithms have been proposed for such approaches. Another approach is to use the linear model obtained from the Born/Rytov approximation and employ mathematical methods like pseudo-inverse, regularization, optode calibration, and other linear transformation techniques to solve the reconstruction problem.24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34 The present work falls in the latter category.
The proposed method is a multistage reconstruction algorithm that reduces the ill-posedness of the inverse problem iteratively and thus improves the image quality. It is assumed that the heterogeneities have only absorptive contrast against the background tissue. The first stage uses a signal subspace method to identify the regions that definitely belong to the background medium. Thus, the region of interest becomes smaller, and the ill-posedness of the problem reduces. The second stage performs a truncated singular value decomposition (SVD)–based pseudo-inverse for the reduced inverse problem. The third stage uses an important property of the truncated SVD-based inversion to guess the regions that belong to the background tissue and thus shrinks the region of interest around the scatterers. The second and third stages are performed iteratively until the region of interest cannot be reduced further based on the property of the truncated SVD-based inversion.
The advantages of the proposed method are as follows:
1. As the results show, the method works well even for complicated geometries, which have multiple heterogeneities present in the domain or have boundaries that are not simply connected.
2. The method needs far fewer measurements than are typically needed by previously reported inverse methods.
3. The method is simple and easy to implement.
4. The method is fast and converges in a few iterations.
Although the proposed method has been presented for the absorptive contrast only in order to remain pertinent to the proposed method, it can be extended to the general case, where both absorptive and scattering contrasts are present. It can also be easily adapted to the Rytov model of diffusion equation.
We first present the mathematical formulation of DOT. Then, we present the proposed algorithm and results in subsequent sections. We present two examples that demonstrate the applicability of the algorithm for complicated structure of heterogeneities. The structures chosen are inspired by the examples presented in Refs. 15, 16, 17. We conclude the paper with discussion and conclusions.
Diffusion Modelis the modulation frequency of the incident photon wave; is the speed of light in the tissue; , , and are, respectively, the absorption coefficient, the diffusion coefficient, and the reduced scattering coefficient of the tissue at a point ; and is the representation of the isotropic point source. We will refer to as the amplitude of the source at . We mention that the photon density and the source term have the time harmonic term , which has been suppressed in the paper and is present implicitly in the physical sense. is assumed to be constant in the tissue sample, and the absorption coefficient is split into two parts: , where is the homogeneous part of , which is equal to the background absorption coefficient, while is its heterogeneous part, also referred to as the contrast at the location . Further, is the Green’s function corresponding to the problem, is the tissue domain of interest containing the heterogeneities, and . The wave number used for defining the Green’s function is .
We assume that sources and detectors are being used for measurement. Each detector is a linear sensor sensing the scattered photon density at the detector position . The tissue may be discretized into voxels, and the absorption coefficient can be assumed to be constant over a voxel, if the discretization is sufficiently fine.
Stage 1: Signal Subspace Method
We form a mapping:is a vector comprising the perturbation in photon density measured at , to , is a vector comprising the source terms to , and is a matrix whose ’th element has the value , where is the volume (area) of the ’th voxel in a problem of 3-D (2-D) geometry, and is the representative location of the ’th voxel.
The matrix can be understood as a linear transformation from the sources to the photon density perturbations measured at the detectors. Thus, the vector lies in the range of . Using the singular value decomposition (SVD)35 of the matrix , , we know that the range of is spanned by the vectors , . Here, vectors and are the left and right singular vectors, respectively, corresponding to the ’th singular value obtained by SVD. We form a matrix that contains the vectors , , which span the range of and therefore form the space in which lies. We specifically mention that during reconstruction, is formed by illuminating the domain using one source [with amplitude ] at a time (thus, has only one nonzero element, the nonzero element being 1) and measuring the perturbations . Then, the column in corresponding to the active source is given by . Subsequently, the matrix is formed by applying SVD on the matrix that contains the actually measured photon density perturbations.
From the definition of , the matrix can be written as:is a matrix of dimension having elements , is a diagonal matrix whose diagonal elements are , and is a matrix of dimension having elements . Equations 3, 4 can be combined as follows:
For simplicity, we denote the ’th column in corresponding to the ’th voxel as . Using Eq. 5, the matrix can be understood as a linear transformation that transforms the photon sources induced at the voxels’ location, , to the photon density perturbations measured at the detectors. It should be noted that photon sources are induced only at the heterogeneous voxels—say, . Thus, physically, the vector lies in a -dimensional vector subspace that is spanned by the vectors .
Combining the physical and mathematical perspectives of the subspace to which belongs,36 it is evident that every vector can be represented as a linear combination of the vectors in :
It is expected that for homogeneous voxels , finding such a combination of the vectors in should not be possible. However, there might be some homogeneous voxels, such that corresponding to them are numerically linear combinations of the vectors . Thus, although they do not contribute physically to , they numerically belong to the subspace spanned by the . Due to this, despite being homogeneous, they can be expressed as linear combinations of vectors in . We denote such voxels depicting an ambiguous behavior as , and the remaining voxels as .
Based on the residue in Eq. 6, we form an error metric as follows:is small for the voxels and , and large for the voxels . By choosing a proper threshold, say, , we can classify the voxels with as , which definitely belong to and have . The concept of stage 1 is summarized in Table 1 .
Summary of stage 1.
|Serialno.||Category||Currentinduced||Result for Eq. 7||Inference based on ε(rm)|
|1||Zero||Definitely homogeneous,classify as|
|2||Zero||May be heterogeneous|
|3||Nonzero||May be heterogeneous|
It should be noted that although this method is seemingly similar to the multiple signal classification (MUSIC) algorithm,36, 37 there are some very important differences. First, MUSIC is a method used for qualitative determination of the location of heterogeneities, while the present method does not attempt to locate the heterogeneities. It rather attempts to find out the voxels that cannot be confused as heterogeneity and thus reduce the region of interest. Second, the conventional formulation of MUSIC is applicable to point-like heterogeneities only, while the current method is applicable to extended heterogeneities as well. Third, MUSIC considers the noise space that is orthogonal to , while the current method considers the range itself.
The choice of a threshold value determines the severity of rejection of the voxels in this stage. However, the suitable threshold changes from problem to problem and setup to setup. Empirically, we found that if the SNR of the system is dB , a threshold value of gives good results. It might be argued that this is a very large value for the least-squares error method applied here. However, it is safer to take a reasonably large value so that the heterogeneous voxels may not get rejected while reducing the ill-posedness of the problem, especially in the presence of noise.
Stage 2: Truncated SVD-Based Pseudoinverse
Here, we construct another mapping:is an -dimensional vector containing the measurements corresponding to all the pairs of sources and detectors. It is thus a concatenation of all the vectors formed by using only one source at a time. contains the elements . The elements of are , , where corresponds to a measurement pair .
We rewrite the mapping 8 as follows by separating the real and imaginary parts of and in order to ensure that obtained by inversion is real-valued:
The relevance of inversion based on this formulation can be emphasized by considering that in the previous stage, we used the concept of measured vector belonging to the range subspace. Although we significantly reduced the number of unknowns by doing so, we were unable to distinguish between the ambiguous voxels and the heterogeneous voxels. Since this stage uses the actual photon density measurements, which correspond to particular vectors in the range subspace, inversion based on it should be able to give a better insight about the heterogeneity of the voxels not rejected.
At this point, it is important to discuss again the relevance of stage 1 on the method. It is evident that the reduction of definitely homogeneous voxels reduces the number of unknowns and thus reduces the ill-posedness of the problem. It also has an impact on the computational complexity of stage 2. Stage 1 uses and to determine and reject the definitely homogeneous voxels. The maximum dimensions of and are both much less than the total number of voxels in the domain and close to the number of detectors. Further, since the matrix does not change for any voxel, the pseudoinverse needs to be calculated only once in stage 1. In stage 2, the maximum dimension of is the number of voxels under consideration—say, . Further, stage 2 requires the calculation of SVD of the matrix , which has a computational complexity of . Evidently, the computational complexity of stage 1 is much less than stage 2. In addition, stage 1 rejects some voxels from further consideration and thus reduces the computational load on stage 2.
Stage 3: Iterative Rejection of Ambiguous Voxels
Since the dimension of is more than the dimension of , Eq. 9 is an underdetermined linear equation. Thus, there are infinite possible solutions to this equation. However, truncated pseudo-inverse of provides a solution among the infinite pool of solutions that has the following properties:
1. In the -dimensional space of (i.e., the domain of ), solution vectors lying in the subspace spanned by the first few right singular vectors (as many as the rank of ) are considered.
2. Among the solutions lying in the subspace discussed earlier, the solution with minimum length (or minimum Frobenius norm) is chosen as the solution of the equation.
We recall that . If all the heterogeneities have similar—say, positive—contrast , then the actual vector will have some negative components (for the heterogeneous voxels) and some zero components (for the homogeneous voxels). Since the truncated pseudo-inverse will find a solution with minimum length, the solution chosen, , is as close to the origin of the dimensional space of as possible. This implies that in comparison to the actual vector , is expected to be smooth vector with small negative and positive values close to zero. Further, the voxels belonging to the heterogeneities or near to them will have negative values for the retrieved , and the voxels away from the heterogeneities are very likely to have positive values for the retrieved . Consequently, the voxels whose retrieved values are positive are very unlikely to belong to . Thus, we add them to the set of rejected voxels , and perform the stages 2 and 3 iteratively until all the elements of are negative.
However, it should be noted that the value of the absorption coefficient of the background cannot be known a priori for biological samples, and the value used is the statistical average available from previous research. Thus, the actual absorption coefficient of the tissue may be slightly larger or smaller than the value used for reconstruction. Due to this, it is more reasonable that we do not reject the voxels based on the negative definiteness of , and instead, we allow some margin on the positive side. Accordingly, we may reject voxels for which , where and is preferably small. Typically, the value of may be chosen to be the relative tolerance level in the estimated statistical average of the absorption coefficient of the background tissue.
We consider two 2-D examples to illustrate the performance of the algorithm. In both examples, the background tissue is of size and has the parameters and . All the heterogeneities have contrast . We have used eight sources and eight detectors, placed in a regular circular fashion around the domain. The radius of this circular arrangement is . The measurement setup is shown in Fig. 1 . The modulation frequency of the photon density waves is .
The tissue has been discretized into square pixels of dimension for the forward problem and for the inverse problem. Different pixel sizes have been chosen in order to avoid the inverse crime. The forward problem has been solved using Eq. 8, and white Gaussian noise has been added to the calculated photon density perturbations. The solution of the forward problem is corrupted by adding 3% white Gaussian noise.
Since we assume the background tissue to extend infinitely, we use the homogeneous 2-D Green’s function , where is the Hankel function of the first kind and zero order. The values of and are chosen to be and 0.05, respectively. All the figures representing the heterogeneities and the reconstruction results plot the actual/reconstructed absorptive contrast for the complete investigation domain.
The first example consists of three circular heterogeneities present in the investigation domain [see Fig. 1]. These heterogeneities are located at , , and , respectively, and have radius each.
The plot of is plotted in Fig. 1, and the voxels rejected as definitely homogeneous are shown in black in Fig. 1. The final reconstruction result is shown in Fig. 1. The result is obtained after 13 iterations. For comparison, we provide the reconstruction result when only the truncated SVD-based pseudo-inverse is used for reconstruction in Fig. 1. It is evident that Fig. 1 presents a very blurred reconstruction and the three heterogeneities are not well-resolved.
The results presented in Figs. 1 and 1 assume that the absorption coefficient of the tissue is known to be exactly . However, as discussed before, often the exact absorption coefficient of the background tissue is not known in advance, and the statistical average needs to be used for calculation. We present two example cases to demonstrate the effectiveness of the algorithm in such cases. In the first case, the actual absorptive contrast is 5% more than the statistical average . The result for this case is obtained after 14 iterations and is presented in Fig. 1. In the second case, the actual absorptive contrast is 5% less than the statistical average . The result for this case is obtained after 18 iterations and is presented in Fig. 1.
The second example consists of an annular heterogeneity present in the investigation domain [see Fig. 2 ]. It is centered at and has inner radius of and outer radius of . The center has been chosen such that there is not symmetry between the source/detector arrangement and the heterogeneity. The plot of is shown in Fig. 2. The reconstruction result is obtained after 12 iterations and is shown in Fig. 2. The reconstruction result using only the truncated SVD-based pseudo-inverse is shown in Fig. 2. As seen in Fig. 2, the inner hollow of the annular tissue is not detected using the truncated SVD-based pseudo-inverse.
Next, we present results for the cases where the actual absorption coefficient of the background tissue is different from the statistical average used for computation. In the first case, the actual absorptive contrast is 5% more than the statistical average . The result for this case is obtained after 14 iterations and is presented in Fig. 2. In the second case, the actual absorptive contrast is 5% less than the statistical average . The result for this case is obtained after 10 iterations and is presented in Fig. 2.
Discussion and Conclusion
While the conventional approach would directly apply some inversion technique (like truncated SVD-based pseudo-inverse) to solve Eq. 8 or 9 for , we have reduced the number of unknowns in in stage 1 and then iteratively in stage 3. The role of stage 1 is to provide an initial guess and accelerate the algorithm. It does not improve the resolution by itself, although it reduces the ill-posedness of the subsequent formulation. The improvement in resolution and quality of imaging is specifically due to the stage 3. Thus, the quality of reconstruction will not be greatly affected if stage 1 is omitted. However, the algorithm will perform slowly and will typically need more iterations to converge. On the other hand, if stage 3 is omitted, the reconstruction using stage 1 and stage 2 is only marginally better than performing stage 2 alone.
The rejection of the voxels definitely belonging to the background has given us two advantages. First, we have reduced the computational intensity of the problem. Second, and more important, we have made the problem of inversion more definitive and guided. Effectively, we have reduced the ill-posedness of the overall inversion problem by reducing the number of unknowns iteratively (due to the rejection of voxels). We have shown the performance of our algorithm in the presence of multiple heterogeneities and heterogeneities with complicated structures that may contain boundaries that are not simply connected.
We also note that although the present analysis has been developed for absorptive contrast only and in the premise of Born approximation, it is easily extensible to the scattering contrast and/or Rytov approximation. Further, although the concept has been demonstrated for an infinitely extending background tissue, it is directly applicable to finite tissue boundaries by the use of the appropriate Green’s function governing the problem scenario.
This work was supported by the Singapore Ministry of Education under Grant No. R263000485112, Life Science Institute/NUS Grant No. R397000615712 and SERC/ Grant No. 0521010098.