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.1 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 studied13,14 along the same line as SVD analysis. The studies on assessment of data collection strategies in terms of data correlation has also been investigated15–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 studies6 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).18 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.
Review of Diffuse Optical Tomography Forward and Inverse Problem
For a given distribution of optical parameters, absorption coefficient , and reduced scattering coefficient , the light propagation in biological tissues is usually modeled by2,3,19 The forward model data are generated by sampling in the data space with the th source and the th detector positions and is represented by 20 Equation (3) can be rewritten as Eq. (4) with representing the compatibility in the initial estimate 3 to the small perturbation in the optical parameters (, ) of the reconstructed image basis (voxels or nodes). In the formulation of data and ISRs, which will be described later, only the part of the Jacobean matrix is considered without loss of generality. Substituting Eq. (5) into Eq. (4), the linearized objective function can be written by 21 The initial guess on the bulk optical properties for calculating the sensitivity matrix is obtained by a calibration procedure22 in the imaging domain whose analytical solutions are readily available.
Singular Value Decomposition
A technique which is of particular interest for investigating ill-conditioned and/or rank-deficient problems is the SVD. In SVD, a matrix is factorized by a set of three matrices
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.23 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.24 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,18 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 “” 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 “” 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:
Data space ranking
Because the singular vectors are orthogonal, they will span an -dimensional space. We form the matrix product “,” where each row contains the squared contribution of the rows of “” in terms of the coordinate system defined by “”:18 However, because of the presence of smaller singular values, it is not straightforward to get the inversion.25 One way to get the inversion is to add a small regularization factor as follows: 24
Other forms of the filter factors are also available in the literature.24 For example, the filter factor for truncated singular value decomposition (TSVD) is given by Eq. (19)
Image space ranking
Similarly, one can define an ISR by forming a matrix “” from the squared contribution of the rows of “” spanned by “”
One can similarly introduce the matrix of fractional squared singular values
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 “” 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 “” 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.
We used two-dimensional (2-D) imaging geometries for simplicity in this work, but it can be extended to a three-dimensional (3-D) problem straightforwardly. In the case of breast DOT imaging, two most commonly considered imaging geometries include: a mammographically compressed breast with the sources being placed at one plane and the detectors on the other,29 or a circular distribution of the sources and the detectors around the breast volume.8 We employed both the rectangular imaging geometry and the circular imaging geometry in this work as summarized in Table 1. For a circular geometry, the sources and the detectors were placed around the object support in an equiangular fashion.
Description of geometries for which SD optimization analysis was performed.
|Types||No. of sources||No. of detectors||No. of projections|
|First dense transmission||15||35||525|
|Second dense transmission||15||70||1050|
|Mouse head geometry||Tomographic||8||8||56|
In all the simulation cases, we fixed the optical properties by: absorption coefficient , scattering coefficient , and for the homogenous background, and absorption coefficient , scattering coefficient , and for the tumor anomalies. One percent (1%) of Gaussian noise in the amplitude and in the phase angle was added for all the cases unless otherwise specified. The simulations were conducted assuming a Gaussian source intensity profile having the full-width-half-maximum of 2 mm.
We adopted a nonlinear iterative algorithm under the Rytov approximation and a dual-mesh scheme to avoid inverse crime for image reconstruction in all the simulation cases. Meshes were constructed as summarized in Table 2.
Description of FEM meshes constructed for numerical study.
|FEM triangular mesh||Circular geometry||Rectangular geometry||Mouse head geometry|
|One anomaly||Two anomaly|
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 (). 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.
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 and 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.
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 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 transmission case.
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.
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 , . 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.
A similar study has been conducted for the case of two anomalies of radius 7.5 mm placed at the locations of (, ) and (, ), 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.
The rectangular geometry was also tested, where a single anomaly of radius 6.5 mm was placed at the position of (, ). 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.
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.
Deleting low-rank SD sensitivity maps improves the conditioning of the sensitivity matrix “” 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.30 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.
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.
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 summarizes the results of the previous simulation studies.
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.31 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 () 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.
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.
RMSE(×10−2) for mouse brain mesh at 50% reduced number of measurements by different selection methods.
|Full (56)||Optimized (28)||Random (28)||ROI optimized (28)|
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.
This work was supported in part by the NRF Grant No. 2013M2A2A9043476 and the R&D Convergence Program of NST (National Research Council of Science and Technology) of Korea (Grant No. CAP-13-3-KERI).
Sohail Sabir received his MS degree in nuclear and quantum engineering from KAIST, Republic of Korea, in 2014. He is currently a PhD candidate in the Department of Nuclear and Quantum Engineering at KAIST. His research interests include radiation safety in radiological imaging, optical imaging, and compressive sensing approaches in medical imaging.
Changhwan Kim received his BS and MS degrees in the Department of Nuclear and Quantum Engineering from KAIST, Republic of Korea, in 2013 and 2015, respectively. Currently, his PhD studies focus on the field of imaging physics including diffuse optical tomography and x-ray computed tomography.
Sanghoon Cho received his BS degree in applied physics from Dankook University, Republic of Korea, in 2015. He is a master’s degree candidate in the Nuclear and Quantum Engineering Department, KAIST, Republic of Korea. His current research interests are devoted to the study of signal processing in diffuse optical tomography for diagnostic imaging applications, nuclear medicine for biological and functional imaging, and scatter correction in cone-beam computed tomography imaging.
Duchang Heo received his PhD in electronics engineering from Korea University, Republic of Korea, in 2003. He joined Samsung Electronics Ltd., in 2004, where he was involved in research on broadband light-wave source for a network system and on a picoprojector incorporated in cellular phones. He joined KERI, Republic of Korea, in 2008, and he is conducting research on femtosecond fiber lasers and diffuse optical tomography.
Kee Hyun Kim received his BS degree in electrical and computer engineering from Seoul National University in 2008 and his MS degree in biomedical engineering from Seoul National University in 2010. In 2015, he became a research engineer in KERI, Republic of Korea. His current research interests are devoted to the study of instrumentation and signal processing in diffuse optical tomography for diagnostic imaging applications.
Jong Chul Ye received his BSc and MSc degrees from Seoul National University, Republic of Korea, and his PhD from Purdue University, West Lafayette. He joined KAIST, Republic of Korea, in 2004, where he is currently KAIST endowed chair professor and the professor of the Department of Bio/Brain Engineering. His current research interests include compressed sensing and statistical signal processing for various imaging modalities.
Seungryong Cho received his BS and MS degrees in physics from Seoul National University in 1995 and 1997, respectively, and his PhD in medical physics from the University of Chicago, USA, in 2009. He is an associate professor in the Department of Nuclear and Quantum Engineering, KAIST, Republic of Korea. His research interests include radiation-based imaging for medicine and industry, multimodality imaging, radiation therapy physics, image reconstruction algorithms, and image processing techniques.