1.

Introduction

Diffuse optical tomography (DOT) is a promising functional imaging modality that uses near infrared (NIR) light in the spectral range of 600 to 1000 nm.¹ Typically, in DOT, NIR light is injected at one location and optical projection measurements are made at other locations either by charge-coupled device (CCD) cameras or photomultiplier tubes (PMT). Generally, it is favored to have large data sets for meeting higher resolution demands and full sensitivity coverage of an imaged object.^2–7 On the other hand, it is discouraged to have dense sampling in light of scanning time and system compactness as will be described in the following. By scanning time, we mean the time associated with the source channel sequencing in a multichannel system where detection can be done in parallel, while for a single-channel system, it is associated with the total number of paired measurements. CCD cameras offer a wide-area detection scheme for acquiring optical projections, but they suffer from noise problems because of a relatively small dynamic range and a lower photo-detection sensitivity to small photon flux signals. PMT and fiber-based photodiode detection systems are less sensitive to the noise issues.^6,8 However, in such focused detection systems, where data are acquired by each detector channel individually, data acquisition time for a successful image reconstruction can be substantially increased. An increased scanning time can complicate the imaging problem due to patient movement and other physiological fluctuations. On the contrary, increasing the number of detection channels for parallel data acquisition certainly would result in a bulky imaging system and is not a cost-effective approach. Therefore, one needs a source–detector distribution optimization that can provide a sufficient sampling density while minimizing the number of channels.

The source–detector (SD) optimization problem has been a topic of wide interest in the optical imaging community, and various approaches have been investigated that do not necessarily solve an inverse problem. Among such approaches, singular value decomposition (SVD) analysis of the sensitivity matrix of the forward problem has been widely studied for designing an appropriate scanning scheme. By counting the number of singular values above a given threshold,^9–12 one can analyze different scanning schemes. Although SVD analysis is conceptually straightforward and can provide useful information, it lacks an insight on the information contained in the singular vectors,^6,13 and therefore, does not provide a way to distinguish the redundant SD pairs in a given imaging task. Similarly, Cramer-Rao lower bound (CRLB) analysis has been studied^13,14 along the same line as SVD analysis. The studies on assessment of data collection strategies in terms of data correlation has also been investigated^15–17 quite recently. Such kinds of optimization studies are certainly helpful for the reduction of data acquisition time. Although some of the optimization algorithms proposed in Refs. 15 and 16 are able to successfully reduce the number of measurements, they require a heuristic threshold parameter to obtain the number of independent measurements. More importantly, it how to adaptably reduce data measurements in a controlled manner at some particular reduction level has not been reported. Additionally, the effect of change in data correlation on neighboring measurements after deleting a particular least informative data point has not been considered, because the existing approaches are based on a single-step ranking process. Moreover, the previously studied approaches mostly focus on the data space properties in terms of the SD geometry optimization and only a few studies⁶ have considered the effect of sampling density on image space properties.

In this work, we quantitatively re-examine the SD optimization issue by introducing mathematical metrics to evaluate the correlation in the projection data, which we refer to as data space ranking in this work, and to evaluate the resolution characteristics or sampling density variation for a particular imaging geometry, which we call image space ranking (ISR).¹⁸ Both mathematical metrics are used to complement the information provided by the SVD analysis. Our approach deals with both data space ranking and ISR for optimizing the SD distribution via an iterative scheme that can accommodate various regularizations in its framework. Prior knowledge can also be incorporated in our scheme to further improve the SD optimization although it is not the scope of this work.

2.

Methods

2.3.

Metrics for DOT Geometry Optimization

The design optimization strategies for well-posed problems are generally based on the covariance matrix, where a suitable design is chosen based on some minimum norms (trace and determinant) of the covariance matrix of the estimate errors. However, such strategies are not suitable for ill-posed problems like DOT where estimators of an unknown parameter are most likely biased because of the regularization.²³ The bias introduced by a regularization does not depend on the noise level and, therefore, cannot be reduced by repeated or redundant measurements usually offered by a CCD-based detection system. In particular circumstances, this nonstochastic component of error, or bias, may be dominant over the stochastic component of measurement noise.²⁴ Therefore, when designing an acquisition scheme, the role of the regularization should also be considered in such a design optimization.

This work aims to suggest the mathematical metrics and an optimization methodology which can be used to objectively analyze the sampling scheme in question and second to use these metrics to select those measurements that produce the least correlation while containing sufficient information for the estimation of the unknown parameter. The methodology is inspired by the work,¹⁸ where the authors applied a ranking scheme to select structural health monitoring sensors. Here, we extend the idea to the DOT problem and propose metrics that utilized the sensitivity matrix which indicate how sensitive the measurement is at a given source–detector location to the change in the optical parameters. However, the system considered in Ref. 18 was an over-determined system, whereas we consider an under-determined and ill-conditioned system. In addition, we use the SVD, which is numerically more stable, instead of the eigenvalue decomposition for the derivation of rankings.

From linear algebra theory, the rank of a matrix is the number of linearly independent rows (or columns) it has, and so rank can be thought as a measure of redundancy. The key idea of the work is to weight each measurement location according to its contribution to the rank of the sensitivity matrix “$J$” to evaluate the effective linear independence of the measurements for a given SD geometry, and similar to weight each image voxel or nodes to its contribution to the rank of sensitivity matrix “$J$” to evaluate the effect of SD geometry in image space. SVD is one of the most widely used methods to efficiently determine the rank of a matrix. The singular vectors are orthonormal, therefore, the following relations hold:

Eq. (10)

$$V{V}^T=I,$$

Eq. (11)

$$U{U}^T=I,$$

Eq. (12)

$$V^T{J}^{T}JV=S^{T}S=S^2,$$

Eq. (13)

$$U^{T}J{J}^{T}U=S{S}^T=S^2.$$

2.3.1.

Data space ranking

Because the singular vectors are orthogonal, they will span an $n$-dimensional space. We form the matrix product “$G$,” where each row contains the squared contribution of the rows of “$J$” in terms of the coordinate system defined by “$V$”:

Eq. (14)

$$G=[JV]\circ [JV],$$

where the symbol $\circ$ is used to represent the Hadamard product. One can similarly rewrite Eq. (14) as follows:

Eq. (15)

$$G=[US]\circ [US].$$

Each column of “$G$” sums to the corresponding squared singular value of “$J$.” Thus, the $i^{\prime }$th column represents the contribution of the $i^{\prime }$th SD measurement pair to the associated squared singular value. If the right-hand side of Eq. (15) is postmultiplied by a matrix of inverse of the squared singular values, each direction in the $n$-dimensional space will have equal importance.¹⁸ However, because of the presence of smaller singular values, it is not straightforward to get the inversion.²⁵ One way to get the inversion is to add a small regularization factor ${\lambda }^2$ as follows:

Eq. (16)

$$F_{SS}=([US]\circ [US])\times {(S^2+{\lambda }^{2}I)}^{-1}.$$

By using the Hadamard product properties of a multiplication of diagonal matrix, Eq. (16) can be rewritten by

Eq. (17)

$$F_{SS}=[U\circ U]\times [S^2\times {(S^2+{\lambda }^{2}I)}^{-1}],$$

where the term in the right bracket is known as the filter factor in Tikhonov (Tikh) regularization ($F_{\mathrm{tikh}}$) sense²⁴

Eq. (18)

$$\begin{gathered} F_{\mathrm{tikh}}=S^2\times {(S^2+{\lambda }^{2}I)}^{-1} \\ {F}_{\mathrm{tikh},mn}=\{\begin{array}{ll}\frac{s_m^2}{s_m^2+{\lambda }_m^2},& m=n\\ 0,& m\ne n\end{array}. \end{gathered}$$

Other forms of the filter factors are also available in the literature.²⁴ For example, the filter factor for truncated singular value decomposition (TSVD) is given by Eq. (19)

Eq. (19)

$$F_{\mathrm{tsvd},ij}=\{\begin{array}{ll}1,& s_i>ϵ\\ 0,& s_i<ϵ\end{array}.$$

where “$ϵ$” corresponds to the truncate on singular values. The $i^{\prime }$th term in the $j^{\prime }$th column of the matrix $F_{SS}$ (referred to as the matrix of fractional squared singular values) represents the fractional contribution of the $i^{\prime }$th SD measurement to the $j^{\prime }$th squared singular value. Summation over each row in Eq. (17) would yield the contribution of the corresponding SD pair to the rank of the sensitivity matrix, and it is referred to as data space ranking (DSR) in this work. Note that the DSR can range between zero and one {0, 1}

Eq. (20)

$$\mathrm{DSR}=\sum _j{U}_{ij}\times \sqrt{F_{jj}}\circ U_{ij}\times \sqrt{F_{jj}},$$

where $F_{jj}$ is a selected Tikh or TSVD filter factor. For a full-rank well-conditioned matrix, where the rank of the matrix can be equal to the number of rows, we do not need to apply the filter factor for regularization and the DSR would produce a vector of ones. However, in the case of DOT where the sensitivity matrix is ill-conditioned and a regularization is necessary, the DSR would range between {0, 1} and provide a useful gadget to analyze the data properties.

2.3.2.

Image space ranking

Similarly, one can define an ISR by forming a matrix “$G$” from the squared contribution of the rows of “$J^T$” spanned by “$U$”

Eq. (21)

$$G=[J^{T}U]\circ [J^{T}U].$$

Rewriting Eq. (21), one can get

Eq. (22)

$$G=[V{S}^T]\circ [V{S}^T].$$

One can similarly introduce the matrix of fractional squared singular values

Eq. (23)

$$F_{SS}=[V\circ V]\times [S^2\times {(S^2+{\lambda }^{2}I)}^{-1}].$$

Summation over each row of Eq. (23) would yield the contribution of the image basis (voxels or nodes) to the rank of the sensitivity matrix, and it is referred to as ISR. The ISR ranges between zero and one, {0, 1}

Eq. (24)

$$\mathrm{ISR}=\sum _j{V}_{ij}\times \sqrt{F_{jj}}\circ V_{ij}\times \sqrt{F_{jj}}.$$

2.3.3.

Sparse selection of source–detector pairs

A schematic flowchart of the sparse selection algorithm of the SD pairs is shown in Fig. 1. The procedure begins with calculating a sensitivity matrix “$J$” for a set of densely sampled SD pair measurements. The maximum number of iterations is determined by the predefined reduction level of measurements. In each iteration step, the DSR is computed for all the remaining SD pair measurements and the SD pair corresponding to the lowest value of DSR is removed. The matrix “$J$” is then updated with the reduced number of SD pairs, and the above procedure is repeated until the desired number of iterations is reached. Ideally, a single SD pair is supposed to be removed during an iteration but more than one can be removed depending on the DSR values. One thing we would like to note is that the methodology can also be applied to the optimum SD design in a multimodality imaging situation,^26–28 where the region-of-interest (ROI) is known as a priori in the anatomical imaging domain. Under such situations, the unknown basis in the sensitivity matrix can be adapted to the ROI; either by simply selecting the basis (voxels or nodes) corresponding to the ROI only or by modifying the basis to a tissue-particular basis through a mapping function, e.g., for breast: fiber-glandular, adipose, and tumor regions.

Fig. 1

Iterative procedure to automatically select most informative sparse SD sensitivity maps.

3.

Results and Discussion

In order to show the effects of the regularizations on DSR, the DSR was plotted as a function of the measured channel in Fig. 2 for the rectangular geometry with 15 sources and 15 detectors. In the Tikh regularization method, the weight parameters were varied ($\lambda =0.5,1$). In the TSVD regularization method, the first 50 and 70 modes with higher singular values were selected. It is observed that the stronger the regularization is, the smaller the DSR becomes. However, the channel by channel DSR variation trend is not very dependent on the regularizations. This implies that a sparse selection of SD pair measurements can be effectively performed even with a suboptimal regularization in action.

Fig. 2

A plot showing the variation of DSR at different SD measurements for the selected regularization methods. (Tikh, TSVD).

Let us consider the mammographically compressed breast imaging scenario in 2-D to optimize the sampling scheme, noting that similar optimization can also be done for the circularly geometry. Field-of-view (FOV) that entails the total coverage in the imaging domain is shown in Fig. 3 for the rectangular geometry. We have varied the detector sampling distance by 1, 2, and 5 mm in the FOV limited by $X=-35\mathrm{mm}$ and $X=35\mathrm{mm}$ to show the effect of data correlation. We have considered the sampling distance of 1 mm to represent the CCD-based wide detection system and sampling distances of 2 and 5 mm to represent the photomultiplier-based focused type detection system. Figure 4 shows the channel by channel DSR and it is observed that the sampling interval of 5 mm provides the least correlated data measurements or the highest DSR. The right- and left-most detection channels present higher DSR values compared to others because the overlap with the neighboring channel is minimum. Similar optimization can also be performed for deciding the intervals of the source array with a single-detection channel.

Fig. 3

Illustration of scanning scheme where the single source is placed on the upper plane and the detector array on the bottom in a rectangular geometry.

Fig. 4

Plot of DSR in the FOV for a single source and various number of detection channels.

Now we consider all the possible combinations of sources and detectors for various scanning conditions summarized in Table 1 for the rectangular geometry. Figure 5 shows the results of DSR. Box plots are also presented in Fig. 6 to compare the DSR among various configurations. It is seen that the median value of DSR is the highest in the case of the $15\times 15$ transmission geometry and that the DSR values are generally skewed toward higher values. The first and third quartile values of the DSR distribution can also be useful. The DSR values of the reflectance case present the higher third quartile value and the maximum value, but the median value is lower compared to the $15\times 15$ transmission case.

Fig. 5

Plot of DSR in the FOV for 15 source channels and various number of detection channels.

Fig. 6

Box plot showing the variation in DSR for different scanning geometries of Fig. 5.

The ISR results are shown in Fig. 7. In Fig. 7(a), the log of ISR at each node is plotted. It can be conjectured that even though the DSR can be higher for the reflectance geometry, the depth resolution would be poorer. It is important to note that, for the case of a densely sampled scheme, ISR is similar to the sparse cases. This implies that the redundancy in data measurements does not essentially contribute to enhancing the image resolution in the case of an ill-conditioned modality like DOT. Figures 7(b) and 7(c) show the line profiles of ISR along the mid horizontal line and the diagonal, respectively, for a more quantitative visualization.

Fig. 7

(a) Log of ISR plotted at each node for transmission (Trans) and reflectance (Ref) type scanning geometries as described in Table 1, (b) ISR plot along the midline, and (c) ISR plot along the diagonal.

The results shown above confirmed that both DSR and ISR can be useful for analyzing sampling schemes in DOT. Now let us move on to the study of iterative sparse selection of the SD pairs using these metrics. As shown in Fig. 8, we considered various reduction rates of data sampling by 50% (120), 65% (84), and 80% (48) of the total number of measurements. The numbers in parenthesis represent the number of SD pairs. The same number of measurements was also randomly selected to compare with the proposed sparse selection methodology. Various noise levels (0%, 1%, 2%, and 5%) have been tested as explained in the method section. Figures 8(a) and 8(b) show the reconstructed absorption and scattering parameter images, respectively, where a single anomaly of radius 10 mm was added at the position of $x=20\mathrm{mm}$, $y=0\mathrm{mm}$. In all cases, the proposed sparse selection method outperformed the random sampling approach. In addition, the reconstructed contrast values of the absorption and scattering in the case of 50% reduction were comparable to those in the full measurement data set. The shape of the anomaly was also better preserved in the case of the proposed selection scheme. In Figs. 9(a) and 9(b), the line profiles across the middle of the reconstructed images are shown in the 1% noise case.

Fig. 8

Reconstructed images of a circular phantom with a single anamoly at various noise levels at various reduction of measurements selected by optimized (OPT) and random (RAN) methods. (a) The reconstructed absorption coefficient images and (b) the reconstructed scattering coefficient images.

Fig. 9

Line profiles across the middle of the reconstructed images for 1% noise case in Fig. 8 for (a) absorption coefficient images and (b) scattering coefficient images.

A similar study has been conducted for the case of two anomalies of radius 7.5 mm placed at the locations of ($x=20\mathrm{mm}$, $y=-10\mathrm{mm}$) and ($x=-20\mathrm{mm}$, $y=10\mathrm{mm}$), respectively, in the circular geometry. The reconstruction results are shown in Fig. 10. Similar observations and conclusions can be drawn in this example as well.

Fig. 10

Reconstructed images of a circular phantom with two anamlies. (a) Absorption coefficient images and (b) scattering coefficient images.

The rectangular geometry was also tested, where a single anomaly of radius 6.5 mm was placed at the position of ($x=0\mathrm{mm}$, $y=6\mathrm{mm}$). The same reduction rates by 50% (112), 65% (80), and 80% (45) have been investigated with respect to the total of 225 measurements. The reconstructed results are shown in Fig. 11.

Fig. 11

Reconstructed images of a rectanugular phantom with 1% noise level at various levels of measurement reduction. (a) Absorption coefficient images and (b) scattering coefficient images.

We further look into the relationship among the DSR, ISR, and the reconstructed results. In Fig. 12, box plots of the DSR distribution in the circular geometry corresponding to Fig. 10 are shown. The median values for the optimized selection scheme were consistently higher than that of the random method. It is presumed that, in the random measurements, the informative or least correlated measurements have also been partly removed from the data set. Figure 12 shows that, at 80% reduction level, we have almost all independent measurements as all the SD pairs cluster around the DSR value of one.

Fig. 12

Box plot of DSR at various number of measurements selected by the optimized (OPT) and the random (RAN) methods.

Deleting low-rank SD sensitivity maps improves the conditioning of the sensitivity matrix “$J$” since the condition number is defined by the ratio of the largest to the smallest singular values. As we keep iterating in the proposed sparse selection scheme, the DSR is increased by removing the redundant information and the sensitivity matrix approaches toward the highest possible rank. The condition number is linked to the rank of a matrix and the highest rank can be ensured by minimizing the condition number.³⁰ This implies that the proposed sparse selection scheme results in minimization of the condition number.

In Figs. 13 and 14, we show the ISR maps and line profiles for various sampling conditions. The optimally selected ones showed more uniform distribution of ISR in the image domain compared to the randomly selected ones.

Fig. 13

Log of ISR showing the resolution characteristics at various reduction of measurements.

Fig. 14

Line profiles of ISR at various number of measurement by the optimized (OPT) and the random (RAN) methods along the mid-central line in Fig. 13.

The bar plot of trace of the ISR is shown in Fig. 15. In an ideal condition, the trace of the ISR should be equal to the rank of the matrix. As the number of measurements decreases, the trace moves closer to the system rank. However, the sampling density will be poorer as the number of measurements decreases and would lead to a degraded image resolution.

Fig. 15

Trace of image space rank at various number of measurements.

For a quantitative assessment, the root-mean-square-error (RMSE) of the reconstructed absorption and scattering parameters of the anomaly region was calculated. The RMSE is given by Eq. (25)

Table 3

RMSE(×10−2) of the reconstructed absorption and scattering coefficient images at various numbers of measurements (NM).

One can also set a criterion to automatically abort the iterative process of sparse selection to the reasonable levels of DSR and ISR. For the cases where the ROI is known as a priori, the number of measurements can be significantly reduced. As earlier explained, the ROI can be either extracted from a multimodality imaging or determined numerically as support recovery in the joint sparsity formulation.³¹ To evaluate the performance of our methodology, we performed a simulation study for a mouse head model. The finite-element meshes were generated from the 2-D MRI image of a mouse head with the segmentation of the brain (lime) and remaining structures (navy) as shown in Fig. 16(a). We added two circular anomalies of radius of 2.5 mm in the brain region (currant) to mimic the tumor. Figure 16(b) shows the reconstructed absorption coefficient image when all the sources and detectors are used. Figure 16(c) shows the reconstructed image when the proposed sparse selection ($\mathrm{NM}=28$) has been applied for the whole volume. Figure 16(d) shows the reconstructed image when the same number of measurements was randomly selected. Figure 16(e) shows the reconstructed image by the proposed sparse selection scheme while limiting the unknown basis to the ROI.

Fig. 16

Reconstructed absorption coefficient images of a mouse head model with the 2-D segmented mouse head image shown in (a).

It is observed that the sparse selection scheme applied to the ROI only provides substantially enhanced image quality compared to that applied to the whole volume. It should be noted here that we did not utilize the prior information in the image reconstruction directly. We utilized a priori information only for designing the SD scheme.

Table 4 summarizes the RMSE calculation results of the reconstructed absorption coefficient images. The presented method is quite general in a sense that it can be applied to any configuration geometry and imaging domain.

Table 4

RMSE(×10−2) for mouse brain mesh at 50% reduced number of measurements by different selection methods.

4.

Conclusions

In this study, we proposed an iterative sparse selection scheme to optimally reduce the number of measurements in DOT. To accomplish the goal, we proposed two metrics, namely data space and ISRs. It was confirmed through the simulation studies that both measures provide useful criteria for deciding optimum sampling intervals among detection channels, scanning field of view, and measurement configuration. Our results demonstrated that a comparable contrast recovery is possible for the optimized sparse configuration of SD pairs compared to the dense configuration. We also showed that the proposed method would apply well to the ROI imaging.