3.1.

Stage 1: Signal Subspace Method

We form a mapping:

The matrix ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ can be understood as a linear transformation from the sources to the photon density perturbations measured at the detectors. Thus, the vector $\overline{\mathbf{V}}$ lies in the range of ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ . Using the singular value decomposition (SVD)³⁵ of the matrix ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ , ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ ${\overline{\mathbf{v}}}_j={\sigma }_j{\overline{\mathbf{u}}}_j$ , we know that the range of ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ is spanned by the vectors ${\overline{\mathbf{u}}}_j$ , ${\sigma }_j\ne 0$ . Here, vectors ${\overline{\mathbf{u}}}_j$ and ${\overline{\mathbf{v}}}_j$ are the left and right singular vectors, respectively, corresponding to the $j$ ’th singular value ${\sigma }_j$ obtained by SVD. We form a matrix $\stackrel{\mathrm{̿}}{\mathbf{U}}$ that contains the vectors ${\overline{\mathbf{u}}}_j$ , ${\sigma }_j\ne 0$ , which span the range of ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ and therefore form the space in which $\overline{\mathbf{V}}$ lies. We specifically mention that during reconstruction, ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ is formed by illuminating the domain using one source [with amplitude $Q\left({\mathbf{r}}_s\right)=\kappa$ ] at a time (thus, $\overline{\mathbf{Q}}$ has only one nonzero element, the nonzero element being 1) and measuring the perturbations $\overline{\mathbf{V}}$ . Then, the column in ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ corresponding to the active source is given by $\overline{\mathbf{V}}$ . Subsequently, the matrix $\stackrel{\mathrm{̿}}{\mathbf{U}}$ is formed by applying SVD on the matrix ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ that contains the actually measured photon density perturbations.

From the definition of ${\mathbf{\Phi }}_{\mathit{scat}}(d,s)$ , the matrix ${\stackrel{\mathrm{̿}}{\mathbf{\Phi }}}_{\mathbf{scat}}$ can be written as:

For simplicity, we denote the $m$ ’th column in ${\stackrel{\mathrm{̿}}{\mathbf{G}}}_{\mathbf{scat}}$ corresponding to the $m$ ’th voxel as $\overline{\mathbf{G}}\left({\mathbf{r}}_m\right)$ . Using Eq. 5, the matrix ${\stackrel{\mathrm{̿}}{\mathbf{G}}}_{\mathbf{scat}}$ can be understood as a linear transformation that transforms the photon sources induced at the voxels’ location, $\stackrel{\mathrm{̿}}{\mathbf{O}}{\stackrel{\mathrm{̿}}{\mathbf{G}}}_{\mathbf{inc}}\overline{\mathbf{Q}}$ , to the photon density perturbations measured at the detectors. It should be noted that photon sources are induced only at the heterogeneous voxels—say, ${\mathbf{r}}_m^{\left(\mathrm{het}\right)}$ . Thus, physically, the vector $\overline{\mathbf{V}}$ lies in a $D$ -dimensional vector subspace that is spanned by the vectors $\overline{\mathbf{G}}\left({\mathbf{r}}_m^{\left(\mathrm{het}\right)}\right)$ .

Combining the physical and mathematical perspectives of the subspace to which $\overline{\mathbf{V}}$ belongs,³⁶ it is evident that every vector $\overline{\mathbf{G}}\left({\mathbf{r}}_m^{\left(\mathrm{het}\right)}\right)$ can be represented as a linear combination of the vectors in $\stackrel{\mathrm{̿}}{\mathbf{U}}$ :

It is expected that for homogeneous voxels $(\delta {\mu }_a=0)$ , finding such a combination $\overline{a}$ of the vectors in $\stackrel{\mathrm{̿}}{\mathbf{U}}$ should not be possible. However, there might be some homogeneous voxels, such that $\overline{\mathbf{G}}\left({\mathbf{r}}_m\right)$ corresponding to them are numerically linear combinations of the vectors $\overline{\mathbf{G}}\left({\mathbf{r}}_m^{\left(\mathrm{het}\right)}\right)$ . Thus, although they do not contribute physically to $\overline{\mathbf{V}}$ , they numerically belong to the subspace spanned by the $\overline{\mathbf{G}}\left({\mathbf{r}}_m^{\left(\mathrm{het}\right)}\right)$ . Due to this, despite being homogeneous, they can be expressed as linear combinations of vectors in $\stackrel{\mathrm{̿}}{\mathbf{U}}$ . We denote such voxels depicting an ambiguous behavior as ${\mathbf{r}}_m^{\left(\mathrm{amb}\right)}$ , and the remaining voxels as ${\mathbf{r}}_m^{\left(\mathrm{hom}\right)}$ .

Based on the residue in Eq. 6, we form an error metric as follows:

Table 1

Summary of stage 1.

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 $\stackrel{\mathrm{̿}}{\mathbf{U}}$ , while the current method considers the range $\stackrel{\mathrm{̿}}{\mathbf{U}}$ itself.

The choice of a threshold value ${\epsilon }_t$ 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 $N$ dB $[N=20{\log }_{10}\left(\mathrm{SNR}\right)]$ , a threshold value of $-20\log \left({\epsilon }_t\right)=N∕6$ 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.

3.2.

Stage 2: Truncated SVD-Based Pseudoinverse

Here, we construct another mapping:

We rewrite the mapping 8 as follows by separating the real and imaginary parts of $\overline{\mathbf{V}}$ and $\stackrel{\mathrm{̿}}{\mathbf{W}}$ in order to ensure that $\overline{\mathbf{O}}$ 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 $\overline{\mathbf{V}}$ 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 $\overline{\mathbf{G}}\left({\mathbf{r}}_m\right)$ and $\stackrel{\mathrm{̿}}{\mathbf{U}}$ to determine and reject the definitely homogeneous voxels. The maximum dimensions of $\overline{\mathbf{G}}\left({\mathbf{r}}_m\right)$ and $\stackrel{\mathrm{̿}}{\mathbf{U}}$ are both much less than the total number of voxels in the domain and close to the number of detectors. Further, since the matrix $\stackrel{\mathrm{̿}}{\mathbf{U}}$ does not change for any voxel, the pseudoinverse needs to be calculated only once in stage 1. In stage 2, the maximum dimension of ${\stackrel{\mathrm{̿}}{\mathbf{W}}}^B$ is the number of voxels under consideration—say, $M^{\prime }$ . Further, stage 2 requires the calculation of SVD of the matrix ${\stackrel{\mathrm{̿}}{\mathbf{W}}}^B$ , which has a computational complexity of $O\left(2D{M}^{\prime 2}\right)$ . 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.

3.3.

Stage 3: Iterative Rejection of Ambiguous Voxels

Since the dimension of $\overline{\mathbf{O}}$ is more than the dimension of ${\overline{\mathbf{V}}}_{SD}^B$ , Eq. 9 is an underdetermined linear equation. Thus, there are infinite possible solutions to this equation. However, truncated pseudo-inverse of ${\stackrel{\mathrm{̿}}{\mathbf{W}}}^B$ provides a solution among the infinite pool of solutions that has the following properties:

We recall that $O\left({\mathbf{r}}_m\right)=-\delta {\mu }_a\left({\mathbf{r}}_m\right)∕\kappa$ . If all the heterogeneities have similar—say, positive—contrast $\left[\delta {\mu }_a\left({\mathbf{r}}_m^{\left(\mathrm{het}\right)}\right)\right]$ , then the actual vector $\overline{\mathbf{O}}$ 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, ${\overline{\mathbf{O}}}^{\mathit{inv}}$ , is as close to the origin of the $M$ dimensional space of $\overline{\mathbf{O}}$ as possible. This implies that in comparison to the actual vector $\overline{\mathbf{O}}$ , ${\overline{\mathbf{O}}}^{\mathit{inv}}$ 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 $O^{\mathit{inv}}\left({\mathbf{r}}_m\right)$ , and the voxels away from the heterogeneities are very likely to have positive values for the retrieved $O^{\mathit{inv}}\left({\mathbf{r}}_m\right)$ . Consequently, the voxels whose retrieved values $O^{\mathit{inv}}\left({\mathbf{r}}_m\right)$ are positive are very unlikely to belong to $\left\{{\mathbf{r}}_m^{\left(\mathrm{het}\right)}\right\}$ . Thus, we add them to the set of rejected voxels $\left\{{\mathbf{r}}_m^{\left(\mathrm{rej}\right)}\right\}$ , and perform the stages 2 and 3 iteratively until all the elements of ${\overline{\mathbf{O}}}^{\mathit{inv}}$ 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 $O^{\mathit{inv}}\left({\mathbf{r}}_m\right)$ , and instead, we allow some margin on the positive side. Accordingly, we may reject voxels for which $O^{\mathit{inv}}\left({\mathbf{r}}_m\right)>c\max \left(\left|O^{\mathit{inv}}\left({\mathbf{r}}_m\right)\right|\right)$ , where $c∊[0,1)$ and is preferably small. Typically, the value of $c$ may be chosen to be the relative tolerance level in the estimated statistical average of the absorption coefficient ${\mu }_a$ of the background tissue.

4.1.

Example 1

The first example consists of three circular heterogeneities present in the investigation domain [see Fig. 1]. These heterogeneities are located at $(-0.5,-0.6)\mathrm{cm}$ , $(-0.6,0.7)\mathrm{cm}$ , and $(0.7,0.3)\mathrm{cm}$ , respectively, and have radius $0.4\mathrm{cm}$ each.

The plot of $20\log \left[\epsilon \left({\mathbf{r}}_m\right)\right]$ 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 ${\mu }_{a0}=0.02{\mathrm{cm}}^{-1}$ . 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 ${\mu }_{a0}=0.02{\mathrm{cm}}^{-1}$ . 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 ${\mu }_{a0}=0.02{\mathrm{cm}}^{-1}$ . The result for this case is obtained after 18 iterations and is presented in Fig. 1.

4.2.

Example 2

The second example consists of an annular heterogeneity present in the investigation domain [see Fig. 2 ]. It is centered at $(0.1,0.2)\mathrm{cm}$ and has inner radius of $0.7\mathrm{cm}$ and outer radius of $1\mathrm{cm}$ . The center has been chosen such that there is not symmetry between the source/detector arrangement and the heterogeneity. The plot of $20\log \left[\epsilon \left({\mathbf{r}}_m\right)\right]$ 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.

Fig. 2

An annular heterogeneity [background tissue properties (statistical estimates): ${\mu }_s^{\prime }=10{\mathrm{cm}}^{-1}$ and ${\mu }_{a0}=0.02{\mathrm{cm}}^{-1}$ ; 3% additive white noise]. (a) The tissue region of interest and the distribution of heterogeneities. The absorptive contrast for the heterogeneities is $\delta {\mu }_a=0.2{\mathrm{cm}}^{-1}$ . (b) The plot of $20\log \epsilon \left({\mathbf{r}}_m\right)$ used to reject the definitely homogeneous voxels in stage 1. The threshold used for stage 1 is $20\log {\epsilon }_t=-5$ . (c) Multistage reconstruction when the background ${\mu }_{a0}$ is the same as the statistical estimate. (d) Truncated SVD-based reconstruction when the background ${\mu }_{a0}$ is the same as the statistical estimate. (e) Multistage reconstruction with 5% underestimated background ${\mu }_{a0}$ . (f) Multistage reconstruction with 5% overestimated background ${\mu }_{a0}$ .

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 ${\mu }_{a0}=0.02{\mathrm{cm}}^{-1}$ . 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 ${\mu }_{a0}=0.02{\mathrm{cm}}^{-1}$ . The result for this case is obtained after 10 iterations and is presented in Fig. 2.