|
1.IntroductionDiffuse optical tomography (DOT) has been developed from simulations and a phantom experimental realm to in vivo experiments such as animal and human subject imaging.1, 2, 3, 4, 5, 6 DOT image quality depends on many practical issues, such as intrinsic tissue heterogeneity, imaging system noise, and improper deployment of sources and detectors on the imaging probe. Many research regarding enhancing the accuracy and reducing the computational cost of the forward model, developing advanced inverse methods and improving some experimental conditions (e.g., source detector calibration) have been reported. 5, 6, 7, 8, 9, 10, 11 Optimizations of source and detector (SD) arrangements have also been studied by some groups in recent years. 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22 The singular value analysis (SVA) technique was widely used for selecting optimized SD sets through evaluating the useful information contained in the weight matrix of a given imaging setup. 13, 14, 15, 16, 17, 18, 19, 20, 21 It can provide a generic estimation of the SD set's performance over the whole imaging domain. However, it cannot directly provide a quantitative analysis of the accuracy of reconstructed values under specified noise conditions, unless the inverse problem is explicitly solved. Besides, the selections of the thresholds for singular value selections depend on the experimental setup and measurement scheme. Such thresholds must be determined case by case.20, 21 The computational cost for SVA of Jacobian matrices in a Matlab® environment is high, especially when dealing with large numbers of matrices of large sizes. In our study, we introduce a rigorous and low computational cost method to select optimized SD sets by directly estimating the accuracy of reconstructed values, quantitatively. Thus, it can distinguish different SD sets’ performance on a quantitative base. In order to select optimized SD sets, we must investigate the performance of each SD set from an SD group in terms of spatial resolutions (location of the perturbation) and optical properties of the perturbations. The accuracy of the estimations of the locations and the optical properties of the perturbations depends on the noise level and also the sensitivity of the measurements to these parameters.23, 24 The Cramer–Rao lower bound25 (CRLB) was introduced by some groups to calculate the precision limit (accuracy) of the estimations of the perturbation depth (single parameter) based on the given forward model, noise model, and noise level.23, 24 In this paper, we adapted the CRLB method into the issue of optimization of SD arrangements in the DOT field. To verify the model we derived, we jointly estimated the precision limits of the depth and the perturbation value of a single target, which was embedded in different depths of a diffusive semi-infinite geometry and with different noise levels. Our estimation method is based on the diffusion equation in the Laplace domain, a Gaussian noise model, and the CRLB. The reliability of the CRLB values was examined through comparisons of the CRLB values and the corresponding sample variances of the same SD sets. Furthermore, we compared the effectiveness of CRLB method to the commonly used SVA method on the same SD sets. Reconstruction images from two simulated SD sets, one with the higher precision limits and the other with the lower precision limits, were presented in this paper to demonstrate the effectiveness of our method for selecting optimized SD arrangements in terms of image qualities. The proposed method will be applied to optimize the SD arrangements on a rotation probe, which will be used in our time domain (TD) DOT system.26, 27 Although all the simulations were conducted in the Laplace domain,28, 29 our method can be adapted to frequency domain and time domain as well. 2.Theory2.1.Forward formulaFor a given Laplace parameter p, the photon density at location r can be given by the first Born approximation,6 as follows: Eq. 1[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \langle {\rm \Phi (}{\bf r}{\rm, }p{\rm)}\rangle = \langle {\rm \Phi }_{{\rm bg}} {\rm (}{\bf r}{\rm, }p{\rm)}\rangle + {\rm \Phi }_{{\rm pert}} {\rm (}{\bf r}{\rm, }p{\rm)}{\rm .} \end{equation}\end{document}Eq. 2[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm \Phi }_{{\rm pert}} {\rm (}{\bf r}{\rm, }p{\rm)} = \int_V {G_{\rm 0} {\rm (}{\bf r}{\rm, }{\bf r}_{\rm d} {\rm, }p{\rm) \delta \mu }_{\rm a} {\rm (}{\bf r}{\rm)}} G_{\rm 0} {\rm (}{\bf r}_{\rm s} {\rm, }{\bf r}{\rm, }p{\rm)}d^3 {\bf r}. \end{equation}\end{document}Eq. 3[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\bm \Phi }_{{\rm pert}} = {\bm {\rm A \delta}}. \end{equation}\end{document}2.2.Cramér–Rao Lower BoundWhile using the same reconstruction method, different sets of measurements using the same SD set lead to different estimations of [TeX:] ${\bm \delta}$ mainly due to the system noises. We examined the accuracy of reconstructed [TeX:] ${\bm \delta}$ from two aspects: the accuracy of estimated reconstructed perturbation value Δμa and the accuracy of estimated perturbation center R. For the single target case in the imaging geometry, Δμa and R can be considered as random variables with standard deviations σ(Δμa) and σ(R). With other conditions held the same, the SD arrangements that result in lower values of σ(Δμa) and σ(R) after adequate repetitions of experiments have a higher possibility of obtaining better quality of reconstruction image and higher noise immunity. However, it is impractical to test all the SD sets from the entire SD group by conducting simulations or repetitive experiments on every set because of the time cost and the change of experimental conditions. Here we introduce the Cramer–Rao lower bound analysis, which combines the forward model solution, the noise model, and the specified noise level together,23, 24, 25 to calculate the lower bound for the deviations σ(Δμa) and σ(R). It is known that in the simplest form of the Cramér–Rao bound analysis, the estimator is assumed to be unbiased, and the estimator in diffuse optical imaging is generally biased, mainly due to the ill-posed inverse problem. However, the ill-posedness of the inverse problem can be minimized in certain situations, especially when a priori information is available. In our method, the target is assumed to be single and there are only two unknown parameters: the center location and the absorption coefficient perturbation. The estimators for them should be approximately unbiased when there are adequate source-detector pairs. Advanced instrumentation methods (e.g., time-resolved DOT) can also help alleviate the problem. Other issues, such as the accuracy of the forward model, are of less importance and their effects could be removed by calibrations. Thus, the CRLB analysis can be adapted to our problem. The following relationships exist: Eq. 4[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm \sigma (\Delta \mu }_{\rm a} {\rm)} \ge \sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }, \qquad {\rm \sigma (R)} \ge \sqrt {{\rm CRLB}_R }, \end{equation}\end{document}In a DOT system, the measurements of photon density obtained from the same source-detector pair obeys Gaussian statistics, with the mean value 〈Φk( [TeX:] ${\bm \theta}$ )〉 and the standard deviation σk( [TeX:] ${\bm \theta}$ ) equal to the noise strength. Subscript k denotes the k’th source-detector pair in one SD set. [TeX:] ${\bm \theta}$ is defined as [TeX:] ${\bm \theta}$ = [Δμa,R]T. Thus, the distribution of a measurement set [TeX:] ${\bm \Phi}$ = [Φ1,Φ2,Φ3, …,ΦM]T is denoted as [TeX:] ${\bm \Phi}$ ∼ N[〈 [TeX:] ${\bm \Phi}$ ( [TeX:] ${\bm \theta}$ )〉,C( [TeX:] ${\bm \theta}$ )], where 〈 [TeX:] ${\bm \Phi}$ ( [TeX:] ${\bm \theta}$ )〉 is the M × 1 mean value vector representing the independent measurements in one SD set and C( [TeX:] ${\bm \theta}$ ) is the diagonal M × M covariance matrix with the main diagonal element formulized as σk( [TeX:] ${\bm \theta}$ )2. The probability density function is written as25 Eq. 5[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \hspace*{-20pt}{\rm P(}{\bm \Phi }{\rm)} &=& \frac{1}{{(2{\rm \pi })^{M/2} {\rm det}^{1/2} {\bf C}({\bm \theta })}}\nonumber\\ &&\times\,{\rm exp}\left[ { - \frac{1}{2}[{\bm \Phi } - \langle {\bm \Phi }({\bm \theta })\rangle ]^T {\bf C}({\bm \theta })^{ - 1} [{\bm \Phi } - \langle {\bm \Phi }({\bm \theta })\rangle ]} \right], \end{eqnarray}\end{document}Eq. 6[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm \sigma }_k ({\bm \theta }) = a\langle {\rm \Phi }_k ({\bm \theta })\rangle, \end{equation}\end{document}Because we considered only one single target, Eq. 3 can be rewritten as Eq. 7[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\bm \Phi }_{{\rm pert}} = {\bm {\rm A}}_1 \cdot {\rm \Delta \mu }_{\rm a} . \end{equation}\end{document}To calculate the lower bounds [TeX:] $\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }$ and [TeX:] $\sqrt {{\rm CRLB}_R }$ for Gaussian observations, the Fisher information matrix (FIM) is introduced with typical elements,25 Eq. 8[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \displaystyle{\rm F}_{i,j}& \mathop = \limits^{{\rm def}}&\left\langle \left({\frac{\partial }{{\partial {\bm \theta }_{\rm i} }}\ln {\rm P(}{\bm \Phi }{\rm)}} \right)\left({\frac{\partial }{{\partial {\bm \theta }_{\rm j} }}\ln {\rm P(}{\bm \Phi }{\rm)}} \right)\Big| {\bm \theta }\right\rangle \nonumber\\ &=& \left[\displaystyle{\frac{{\partial \langle {\bm \Phi }({\bm \theta })\rangle }}{{\partial {\bm \theta }_{\rm i} }}} \right]^T {\bf C}({\bm \theta })^{ - 1} \left[\displaystyle{\frac{{\partial \langle {\bm \Phi }({\bm \theta })\rangle }}{{\partial {\bm \theta }_{\rm j} }}} \right]\nonumber\\ && \, + \displaystyle\frac{1}{2}{\rm tr}\left[ {{\bf C}({\bm \theta })^{ - 1} \frac{{\partial {\bf C}({\bm \theta })}}{{\partial {\bm \theta }_{\rm i} }}{\bf C}({\bm \theta })^{ - 1} \frac{{\partial {\bf C}({\bm \theta })}}{{\partial {\bm \theta }_{\rm j} }}} \right], \end{eqnarray}\end{document}
[TeX:]
\documentclass[12pt]{minimal}\begin{document}\begin{equation*}
{\bf F}\mathop = \limits \left[ \begin{array}{l}
F_{11}, \quad F_{12} \\[3pt]
F_{21}, \quad F_{22} \\
\end{array} \right]
\end{equation*}\end{document}
with
Eq. 9[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} {\rm F}_{{\rm 1,1}} &=& \sum_{k = 1}^{X_1, \ldots, X_n } {\left({2 + \frac{1}{{a^2 }}} \right)\frac{{\left({\sum\nolimits_{n \subset N_1 } {{A}_{k,n} } } \right)^2 }}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }}} {\rm, }\nonumber\\ \hspace*{-3pt}{F}_{{\rm 1,2}}\!=\!{F}_{{\rm 2,1}}\!&=&\! \sum_{k\! =\! 1}^{M} {\!\!\left(\! {2\!\! +\!\! \frac{1}{{a^2 }}}\! \right)\!\!\frac{\Delta\mu_{\rm a}{\left({\sum\nolimits_{n \subset N_1 } {{A}_{k,n} } } \right){\left(\partial_{\rm R}\sum\nolimits_{n \subset N_1 } {{A}_{k,n} }\right)}}}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }}} {\rm, }\nonumber\\ \hspace*{-12pt}F_{{\rm 2,2}} &=& \sum_{k = 1}^M {\left({2 + \frac{1}{{a^2 }}} \right)\frac{{{\rm \Delta \mu }_{\rm a}^{\rm 2} \left({\partial _{\rm R} \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 }}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }}} {\rm, } \end{eqnarray}\end{document}Eq. 10[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{l} \hspace*{-34pt}\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }\; \mathop = \limits^{{\rm def}}\; \sqrt {\displaystyle\frac{{F_{2,2} }}{{{\rm det}\left({\bf F} \right)}}} \vspace*{8pt}\\ \hspace*{16.5pt} = \sqrt {\left({2 + \displaystyle\frac{1}{{a^2 }}} \right)^{ - 1} \displaystyle\frac{{\sum\nolimits_{{\rm k} = {\rm 1}}^{\rm M} {\left[ {\left({\partial _{\rm R} \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 /\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 } \right]} }}{{\sum\nolimits_{\left\{ {k,j} \right\} \subset C} {\displaystyle\frac{{\left[ {\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right) - \left({\sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)} \right]^2 }}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 \langle {\rm \Phi }_j \left({\bm \theta } \right)\rangle ^2 }}} }}} . \\ \end{array} \end{equation}\end{document}Eq. 11[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{l} \hspace*{-25pt}\sqrt {{\rm CRLB}_R }\; \mathop = \limits^{{\rm def}}\; \sqrt {\displaystyle\frac{{F_{1,1} }}{{{\rm det}\left({\bf F} \right)}}}\hspace*{25pt} \vspace*{8pt}\\ \hspace*{18pt} = \sqrt {\left({2 + \displaystyle\frac{1}{{a^2 }}} \right)^{ - 1} \displaystyle\frac{{\sum\nolimits_{k = 1}^M {\left[ {\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 /\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 } \right]} }}{{{\rm \Delta \mu }_{\rm a}^{\rm 2} \sum\nolimits_{\left\{ {k,j} \right\} \subset C} {\displaystyle\frac{{\left[ {\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right) - \left({\sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)} \right]^2 }}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 \langle {\rm \Phi }_j ({\bm \theta })\rangle ^2 }}} }}} {\rm }{\rm .}\hspace*{-18pt} \\ \end{array} \end{equation}\end{document}[TeX:] $\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }$ represents the cross-coupling precision limit for Δμa when taking the influence of unknown R into account. [TeX:] $F_{11}^{ - 1/2}$ is the precision limit for Δμa when R is known. Same rules apply for [TeX:] $\sqrt {{\rm CRLB}_{R} }$ and [TeX:] $F_{22}^{ - 1/2}$ . In our case, Δμa and R affect the signal in a coupling way. Note that the true values of δμa and r are independent; however, the covariance of the random variables Δμa and R may not equal to zero. The lowest value of the covariance between Δμa and R can be obtained from the antidiagonal of F −1. 3.Simulation Results and DiscussionWe conducted simulation to compare the performance of different SD arrangements on a rotation probe shown in Fig. 1. The rotation center is at the center of the probe. There are 23 source positions and 12 detector positions on the probe plane for selection. The examined geometry was chosen to be semi-infinite, and the dimensions of the imaging domain were 3.5×3.5×3.6 cm. The imaging domain started from the probe surface. We discretized the imaging domain into voxels with the size of 0.5×0.5×0.4 cm. The total number of voxels was 441. The rotation probe and the imaging domain were coaxial. We set the angle of rotation at 40 deg. All possible combinations of three sources and four detectors (876,645 sets) from the SD group were examined. Some works14, 15 have suggested choosing identical quantities of sources and detectors to achieve better image qualities. Our proposed method can be used to compare selections of different number of sources and detectors, and there is no restriction either on the combinations of sources and detectors or on the number of sources and detectors. However, in our case, we focused on choosing better SD arrangements of given numbers of sources and detectors for an illustration. The Laplace parameters used in our model ranged from –0.4×109 S−1 to 0.6×109 S−1 in step size of 0.2×109 S−1. The absorption coefficient and reduced scattering coefficient of the homogenous background were set to be 0.02 and 6 cm−1, respectively, which were close to the optical properties for normal breast tissue.11 The refractive index of the background was 1.4, which was close to the tissue refractive index.27 We focused on investigating the precision limits of Δμa and depth L [simply replace ∂R with ∂L in Eqs. 10, 11] of the central targets, in order to select those SD sets that are more sensitive to the perturbations in the central area and have higher robustness in noisy environment. In our simulations, the central targets were defined as those with their centers along the axis passing through the center of the cylindrical rotation probe. A central target of 3×3×3 voxels with the absorption coefficient equal to 0.22 cm−1 and reduced scattering coefficient equal to 6 cm−1 was embedded into the homogeneous diffusive background at different depths. The corresponding Cramer–Rao lower bound was calculated for noise levels of 1, 2, and 3% of the signal. Figure 2 illustrates the parameters’ precision limits for all the examined SD sets. From Fig. 2, one can observe that with the same target depth in Figs. 2a and 2c, the ranking orders of SD sets, based on the corresponding precision limits of both parameters, are similar for different noise levels. The same rules apply for the cases of different target depths with the same noise level. Some small differences in the ranking of some local SD sets, with different noise levels or different target depths, are due to numerical errors from Matlab calculations of precision limits of those local SD sets. These precision limits from different local SD sets have very close values to each others and thus lead to small changes in the ranking of local SD sets, because we cannot avoid the numerical errors. From our observations, we deduce that those SD sets with lower precision limits in a certain environment (with specified noise level or with specified target depth in the central area), comparing to the rest of SD sets, have a more robust performance across the central area. This conclusion can be extended to the noncentral area: those SD sets with lower precision limits in a certain environment (with specified noise level or with specified target depth), comparing to the rest of SD sets, have a more robust performance across the z-axis (while keeping the XY position of the target constant). This is because of the composition of Eqs. 10, 11. Furthermore, for other imaging geometry, as long as the single volume target moves along the normal vector of the surface, the ranking of different SD sets based on the precision limits of perturbation depths and values remains similar under different noise conditions and with different target depths. In our case, we only focused on the central area because of the difficulty of imaging this area. To verify the correlation between the performance of the SD sets in terms of accuracy of reconstructed parameters and the precision limits calculated from the proposed method, a Levenberg–Marquardt algorithm30 with positivity constraint was used to reconstruct the optical properties of the imaging domain under simulated noise conditions. We randomly chose 500 sets with index numbers smaller than 600, 500 sets with index numbers between 465,000 and 475,000, and another 500 sets with index numbers larger than 870,000, according to the SD sets’ index numbers from Fig. 2, to conduct the investigations. Figure 3 shows the sample variances of the target perturbation values after 20 times’ repetitions of simulations with 1, 2, and 3% noise levels, respectively, and the CRLB values calculated from the proposed method for different SD sets. Figure 4 shows the sample variances of the target center and the corresponding CRLB values with the same simulation conditions as the data shown in Fig. 3. From these two figures, we observed that those SD sets with lower precision limits throughout all the tested target depths have a higher chance of achieving higher accuracy of reconstructed parameters. Thus, the signal sensitivity from these SD sets to the reconstructed parameters is higher. We also observed that there are some of sample variances smaller than the corresponding CRLB. This is mainly due to numerical errors from calculation. In addition, for those SD sets with close precision limits, the sample variances of the reconstructed parameters are not strictly following the same ranking of the precision limits. This is due to the fluctuations of the simulation data generated with simulated random noises and also the numerical errors from calculations. In a word, for those SD sets with similar precision limits, the performances of them in terms of the accuracy of the reconstructed parameters are comparable. The effectiveness of the CRLB method and the commonly used SVA method for selecting optimized SD sets were also compared. Figure 5 shows the SVA analysis of the total number of useful singular values above a threshold of 10−4 and the sum of them for the same group of SD sets as in Figs. 3 and 4. Note that the trend of useful singular value numbers of different SD sets is not affected by the value of threshold, although the actual numbers of useful singular value change when using different threshold values. 13, 14, 15, 16, 17, 18, 19, 20, 21 It is known that the SVA method is a widely used method for evaluating the performance of an imaging probe configuration. Its judgement is based on the number of useful singular values above a specified threshold. Figures 3, 4, 5 show that for two SD sets with significant different numbers of useful singular values, the SD set with larger numbers of useful singular values is more optimized than the other. Also from Figs. 3, 4, 5, we observed that for those SD sets with similar numbers of useful singular values, their performances are comparable. Similar observations are found in the CRLB case: for those SD sets with similar CRLB values of examined parameters, their performances are comparable. Nevertheless, the SVA method cannot provide quantitative analysis on particular parameters, or under different conditions. For example, one cannot estimate how low the precision limit of the estimated target position can be or whether the performance of the SD set under different noise conditions is stable or not. In contrast, the CRLB method provides quantitative analysis of the interested parameters. Besides, the computational time for the CRLB analysis is much lower than that for the SVA analysis. This was demonstrated through our simulations. By using a PC with Intel(R) Core(TM)2 CPU 6300 at 1.86 GHz, the computation time for calculating [TeX:] $\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }$ and [TeX:] $\sqrt {{\rm CRLB}_R }$ of one SD set was ∼0.005 s, while the computation time for SVA of the SD set's Jacobian matrix was ∼0.450 s, which is nine times longer than our method. To summarize, qualitatively, the accuracy of the SVA method and that of the CRLB method are comparable; quantitatively, only the CRLB method can provide meaningful estimations of particular parameters under different setting and conditions. To provide a direct image about the arrangements of the sources and detectors, we choose six SD sets for illustration. Figures 6a, 6b, 6c list three SD sets with the three lowest precision limits of very close values, and Figs. 6d, 6e, 6f list the other three SD sets with the three highest precision limits of very close values based on the data from Fig. 2. These six sets performed stable under other conditions as shown in Fig. 2. From Figs. 6a, 6b, 6c, one can see that the distances from the probe center to the sources or the detectors are close to each other within the same SD set; the sources and the detectors are distributed dispersedly from the XY center. From Figs. 6d, 6e, 6f, one can see that the distances from the probe center to the sources or the detectors differ from each other: some detectors are too close to the probe center. Besides, the sources are too close to each other as well as the detectors. Because we were observing the central area, we may infer that those SD sets with dispersed distributions of sources and detectors and proper distances (neither too far nor too close) from sources and detectors to the probe center perform better than others. Reconstructed images of different central targets with different noise levels at the YZ plane (X = 0) using SD sets 1 [Fig. 6a] and 6 [Fig. 6f] are presented in Fig. 7 to illustrate the differences in image qualities caused by the differences in SD arrangements. Figure 7 shows that the image quality (in terms of accuracy of reconstructed target depths, values and background artificial effects) is better using the SD set 1 than that using the SD set 6 for all the examined target depths and the examined noise levels. All images use the same color bar for the data range from 0.02 to 0.22 cm−1 for easy comparing. Detailed data of the reconstructed values is presented in Table 1. From our investigations, the SD sets with lower precision limits for all the conditions we applied can be chosen for the optimization purpose. Table 1Reconstruction values from Fig. 7.
4.ConclusionWe have introduced a rigorous and computationally efficient methodology for selecting optimized source and detector arrangements for DOT systems. Simulations were conducted on a rotation probe hosting three sources and four detectors to validate the effectiveness of the proposed method. The various performances of different source and detector sets were investigated based on precision limits of the target depth and target perturbation values. The SD sets corresponding to the lowest precision limits can be selected for designing an optimized DOT imaging probe, which leads to the best possible image quality. We also discussed the advantages of our method over the SVA method. Our method can be adapted to other imaging geometries, different numbers of sources and detectors and different combinations of sources and detectors. AppendicesAppendixThe derivations of Eqs. 8, 9, 10 presented are shown here. Because we only consider the single-target case, there is only one perturbation value, and Eq. 3 can be rewritten as [TeX:] ${\bm \Phi }_{{\rm pert}} = {\bm {\rm A}}_1 \cdot {\rm \Delta \mu }_{\rm a}$ . [TeX:] ${\bm {\rm A}}_1$ is a M × 1 vector with each element of it formulized as [TeX:] $\sum\nolimits_{n \subset N_1 } {A_{k,n} }$ , where n denotes the voxel index number, N 1 is the set of voxel index numbers of a single target of 3×3×3 voxels, k is the index of the k’th measurement from one SD set, and M is the total number of measurements from one SD set. By substituting Eqs. 1, 3, 6 into Eq. 7 and deriving from it, we get the first derivative of [TeX:] $\langle \Phi({\bm \theta})\rangle$ and [TeX:] $C({\bm \theta})$ with respect to ∂Δμa and ∂R: Eq. 12[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\partial \langle {\bm \Phi }({\bm \theta })\rangle }}{{\partial {\rm \Delta \mu }_{\rm a} }} = \frac{{\partial {\bm \Phi }_{{\rm pert}} ({\bm \theta })}}{{\partial {\rm \Delta \mu }_{\rm a} }} = {\bf {\rm A}}_1; \end{equation}\end{document}Eq. 13[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\partial \langle {\bm \Phi }({\bm \theta })\rangle }}{{\partial R}} = \frac{{\partial {\bm \Phi }_{{\rm pert}} ({\bm \theta })}}{{\partial R}} = {\rm \Delta \mu }_{\rm a} \cdot \partial _R {\bf {\rm A}}_1; \end{equation}\end{document}Eq. 14[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\partial {\bf C}({\bm \theta })}}{{\partial {\rm \Delta \mu }_{\rm a} }} = 2a^2 \cdot {\rm diag}\left[ {\langle {\bm \Phi }({\bm \theta })\rangle \cdot {\bf {\rm A}}_1 } \right]^T; \end{equation}\end{document}Eq. 15[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\partial {\bf C}({\bm \theta })}}{{\partial {R}}} = 2a^2 {\rm \Delta \mu }_{\rm a} \cdot {\rm diag}\left[ {\langle {\bm \Phi }({\bm \theta })\rangle \cdot \partial _R {\bm {\rm A}}_1 } \right]^T . \end{equation}\end{document}Eq. 16[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&\left[ {\displaystyle\frac{{\partial \langle {\bm \Phi }({\bm \theta })\rangle }}{{\partial {\rm \Delta \mu }_{\rm a} }}} \right]^T {\bf C}({\bm \theta })^{ - 1} \left[ {\displaystyle\frac{{\partial \langle {\bm \Phi }({\bm \theta })\rangle }}{{\partial R}}} \right] = {\rm \Delta \mu }_{\rm a} \cdot {\bm {\rm A}}_1 ^T {\bf C}({\bm \theta })^{ - 1} \partial _R {\bm {\rm A}}_1 \nonumber\\[6pt] &&\;\; = \sum_{k = 1}^M {\displaystyle\frac{{{\rm \Delta \mu }_{\rm a} \left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)}}{{a^2 \langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }};} \end{eqnarray}\end{document}Eq. 17[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \hspace*{-15pt}&&\frac{1}{2}{\rm tr}\left[ {{\bf C}\left({\bm \theta } \right)^{ - 1} \frac{{\partial {\bf C}\left({\bm \theta } \right)}}{{\partial {\rm \Delta \mu }_{\rm a} }}{\bf C}\left({\bm \theta } \right)^{ - 1} \frac{{\partial {\bf C}\left({\bm \theta } \right)}}{{\partial R}}} \right]\nonumber\\[6pt] \hspace*{-15pt}&&\quad = \sum_{k = 1}^M {\frac{{2{\rm \Delta \mu }_{\rm a} \left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)}}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }}} . \end{eqnarray}\end{document}Thus, Eq. 18[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&\displaystyle F_{{\rm 1,2}} = F_{{\rm 2,1}} = \sum_{k = 1}^M {\frac{{{\rm \Delta \mu }_{\rm a} \left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)}}{{a^2 \langle {\rm \Phi }_{\rm k} ({\bm \theta })\rangle ^2 }}} \nonumber\\[3pt] &&\quad\displaystyle + \sum_{k = 1}^M {\frac{{2{\rm \Delta \mu }_{\rm a} \left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)}}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }}} \nonumber\\[3pt] &&\quad\hspace*{-17pt}\displaystyle = \sum_{k = 1}^M {\left({2 + \frac{1}{{a^2 }}} \right)\frac{{{\rm \Delta \mu }_{\rm a} \left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)}}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }}}. \end{eqnarray}\end{document}Using Eqs. 7, 9, [TeX:] $\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }$ in Eq. 4 can be written as Eq. 19[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }\, \mathop = \limits^{{\rm def}}\, \sqrt {\frac{{F_{2,2} }}{{{\rm det}\left({\bf F} \right)}}} = \sqrt {\frac{{F_{2,2} }}{{F_{1,1} F_{2,2} - F_{1,2}^2 }}} . \end{equation}\end{document}
[TeX:]
\documentclass[12pt]{minimal}\begin{document}\begin{eqnarray*}
\displaystyle F_{1,1} F_{2,2} &=& \left({2 + \frac{1}{{a^2 }}} \right)^2 {\rm \Delta \mu }_{\rm a}^{\rm 2} \left[ {\sum_{k = 1}^M {\displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 }}{{\langle {\rm \Phi }_{\rm k} ({\bm \theta })\rangle ^2 }}} } \right]\left[ {\sum_{k = 1}^M {\displaystyle\frac{{\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 }}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }}} } \right] \\[6pt]
\displaystyle &=& \left({2 + \frac{1}{{a^2 }}} \right)^2 {\rm \Delta \mu }_{\rm a}^{\rm 2} \left\{ \begin{array}{l}
\displaystyle\sum_{k = 1}^M {\displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 \left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 }}{{\langle {\rm \Phi }_{\rm k} ({\bm \theta })\rangle ^4 }}} \\[10pt]
+ \displaystyle\sum_{\left\{ {k,j} \right\} \subset C} {\left[ {\displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 \left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)^2 }}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 \langle {\rm \Phi }_j ({\bm \theta })\rangle ^2 }} + \displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)^2 \left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 }}{{\langle {\rm \Phi }_j ({\bm \theta })\rangle ^2 \langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }}} \right]} \\
\end{array} \right\} \\
\end{eqnarray*}\end{document}
[TeX:]
\documentclass[12pt]{minimal}\begin{document}\begin{eqnarray*}
\displaystyle F_{1,2}^2 &=& \left({2 + \frac{1}{{a^2 }}} \right)^2 {\rm \Delta \mu }_{\rm a}^{\rm 2} \left[ {\sum_{k = 1}^M {\displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)}}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }}} } \right]^2 \\[6pt]
\displaystyle &=& \left({2 + \displaystyle\frac{1}{{a^2 }}} \right)^2 {\rm \Delta \mu }_{\rm a}^{\rm 2} \left\{ \begin{array}{l}
\displaystyle\sum_{k = 1}^M {\displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 \left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 }}{{\langle {\rm \Phi }_{\rm k} ({\bm \theta })\rangle ^4 }}} \\[10pt]
+ 2\sum_{\left\{ {k,j} \right\} \subset C} {\left[ {\displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)}}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 \langle {\rm \Phi }_j ({\bm \theta })\rangle ^2 }}} \right]} \\
\end{array} \right\}, \\
\end{eqnarray*}\end{document}
we obtain
[TeX:]
\documentclass[12pt]{minimal}\begin{document}\begin{eqnarray*}
\displaystyle F_{1,1} F_{2,2} - F_{1,2}^2 &=& \left({2 + \frac{1}{{a^2 }}} \right)^2 {\rm \Delta \mu }_{\rm a}^{\rm 2} \sum_{\left\{ {k,j} \right\} \subset C} {\left[ \begin{array}{l}
\displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 \left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)^2 }}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 \langle {\rm \Phi }_j ({\bm \theta })\rangle ^2 }} + \displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)^2 \left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)^2 }}{{\langle {\rm \Phi }_j ({\bm \theta })\rangle ^2 \langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 }} \\[18pt]
- 2\displaystyle\frac{{\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)}}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 \langle {\rm \Phi }_j ({\bm \theta })\rangle ^2 }} \\[3pt]
\end{array} \right]} \\[5pt]
\displaystyle &=& \left({2 + \displaystyle\frac{1}{{a^2 }}} \right)^2 {\rm \Delta \mu }_{\rm a}^{\rm 2} \sum_{\left\{ {k,j} \right\} \subset C} {\displaystyle\frac{{\left[ {\left({\sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right) - \left({\sum\nolimits_{n \subset N_1 } {A_{j,n} } } \right)\left({\partial _R \sum\nolimits_{n \subset N_1 } {A_{k,n} } } \right)} \right]^2 }}{{\langle {\rm \Phi }_k ({\bm \theta })\rangle ^2 \langle {\rm \Phi }_j ({\bm \theta })\rangle ^2 }}}, \\
\end{eqnarray*}\end{document}
where {k,j} is a subset of C and C includes all the subsets of combinations of two index numbers out of M index numbers. M is the total number of measurements from one SD set. Thus, we obtain Eqs. 10, 11 from the derivations above.ReferencesV. Ntziachristos, A. G. Yodh, M. Schnall, and
B. Chance,
“Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,”
Proc. Natl. Acad. Sci. U.S.A., 97 2767
–2772
(2000). https://doi.org/10.1073/pnas.040570597 Google Scholar
B. W. Pogue, S. P. Poplack, T. O. McBride, W. A. Wells, K. S. Osterman, U. L. Osterberg, and
K. D. Paulsen,
“Quantitative hemoglobin tomography with diffuse near-infrared spectroscopy: Pilot results in the breast,”
Radiology, 218
(1), 261
–266
(2001). https://doi.org/http://radiology.rsna.org/content/218/1/261.abstract Google Scholar
S. B. Colak, M. B. vander Mark, G. W. Hooft, J. H. Hoogenraad, E. S. Van Der Linden, and
F. A. Kuijpers,
“Clinical optical tomography and NIR spectroscopy for breast cancer detection,”
IEEE J. Sel. Top. Quantum Electron., 5 1143
–1158
(1999). https://doi.org/10.1109/2944.796341 Google Scholar
S. Fantini, M. Franceschini, E. Gratton, D. Hueber, W. Rosenfeld, D. Maulik, P. Stubblefield, and
M. Stankovic,
“Non-invasive optical mapping of the piglet brain in real time,”
Opt. Express, 4
(8), 308
–314
(1999). https://doi.org/10.1364/OE.4.000308 Google Scholar
D. A. Boas, T. Gaudette, and
S. R. Arridge,
“Simultaneous imaging and optode calibration with diffuse optical tomography,”
Opt. Express, 8
(5), 263
–270
(2001). https://doi.org/10.1364/OE.8.000263 Google Scholar
S. R. Arridge,
“Topical review: optical tomography in medical imaging,”
Inv. Prob., 15 R41
–R93
(1999). https://doi.org/10.1088/0266-5611/15/2/022 Google Scholar
A. Serdaroglua, B. Yazicia, and
K. Kwona,
“Optimum source design for detection of heterogeneities in diffuse optical imaging,”
Proc. SPIE, 6139 61391A
(2006). https://doi.org/10.1117/12.647209 Google Scholar
A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and
D. A. Boas,
“Tomographic optical breast imaging guided by three-dimensional mammography,”
App. Opt., 42
(25), 5181
–5190
(2003). https://doi.org/10.1364/AO.42.005181 Google Scholar
D. A. Boas, A. M. Dale, and
M. A. Franceschini,
“Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,”
NeuroImage, 23 S275
–S288
(2004). https://doi.org/10.1016/j.neuroimage.2004.07.011 Google Scholar
J. J. Stott, J. P. Culver, S. R. Arridge, and
D. A. Boas,
“Optode positional calibration in diffuse optical tomography,”
App. Opt., 42
(16), 3154
–3162
(2003). https://doi.org/10.1364/AO.42.003154 Google Scholar
M. Huang, T. Xie, N. Chen, and
Q. Zhu,
“Simultaneous reconstruction of absorption and scattering maps with ultrasound localization: Feasibility study using transmission geometry,”
Appl. Opt., 42
(19), 4102
–4114
(2003). https://doi.org/10.1364/AO.42.004102 Google Scholar
Q. Wang, Z. Liu, and
H. Jiang,
“Optimization and evaluation of a three-dimensional diffuse optical tomography system for brain imaging,”
J. X-Ray Sci. Technol., 15
(4), 223
–234
(2007). https://doi.org/http://iospress.metapress.com/content/w6g92h436l7000u7/ Google Scholar
H. Xu, H. Dehghani, B. W. Pogue, R. Springett, K. D. Paulsen, and
J. F. Dunn,
“Near-infrared imaging in the small animal brain: optimization of fiber positions,”
J. Biomed. Opt., 8
(1), 102
–110
(2003). https://doi.org/10.1117/1.1528597 Google Scholar
E. E. Graves, J. P. Culver, J. Ripoll, R. Weissleder, and
V. Ntziachristos,
“Singular-value analysis and optimization of experimental parameters in fluorescence molecular tomography,”
J. Opt. Soc. Am. A, 21
(2), 231
–241
(2004). https://doi.org/10.1364/JOSAA.21.000231 Google Scholar
P. K. Yalavarthy, H. Dehghani, B. W. Pogue, and
K. D. Paulsen,
“Critical computational aspects of near infrared circular tomographic imaging: Analysis of measurement number, mesh resolution and reconstruction basis,”
Opt. Express, 14
(13), 6113
–6127
(2006). https://doi.org/10.1364/OE.14.006113 Google Scholar
J. P. Culver, V. Ntziachristos, M. J. Holboke, and
A. G. Yodh,
“Optimization of optode arrangements for diffuse optical tomography: A singular-value analysis,”
Opt. Lett., 26
(10), 701
–703
(2001). https://doi.org/10.1364/OL.26.000701 Google Scholar
A. Joshi, W. Bangerth, and
E. M. Sevick-Muraca,
“Non-contact fluorescence optical tomography with scanning patterned illumination,”
Opt. Express, 14
(14), 6516
–6534
(2006). https://doi.org/10.1364/OE.14.006516 Google Scholar
E. E. Graves, J. Ripoll, R. Weissleder, and
V. Ntziachristos,
“A submillimeter resolution fluorescence molecular imaging system for small animal imaging,”
Med. Phys., 30
(5), 901
–911
(2003). https://doi.org/10.1118/1.1568977 Google Scholar
E. E. Graves, J. P. Culver, J. Ripoll, R. Weissleder, and
V. Ntziachristos,
“Singular-value analysis and optimization of experimental parameters in fluorescence molecular tomography,”
J. Opt. Soc. Am. A, 21
(2), 231
–241
(2004). https://doi.org/10.1364/JOSAA.21.000231 Google Scholar
T. Lasser and
V. Ntziachristos,
“Optimization of 360 projection fluorescence molecular tomography,”
Med. Image Anal., 11 389
–399
(2007). https://doi.org/10.1016/j.media.2007.04.003 Google Scholar
A. T. N. Kumar, S. B. Raymond, G. Boverman, D. A. Boas, and
B. J. Bacskai,
“Time resolved fluorescence tomography of turbid media based on lifetime contrast,”
Opt. Express, 14
(25), 12255
–12270
(2006). https://doi.org/10.1364/OE.14.012255 Google Scholar
B. W. Pogue, T. O. McBride, U. L. Osterberg, and
K. D. Paulsen,
“Comparison of imaging geometries for diffuse optical tomography of tissue,”
Opt. Express, 4
(8), 270
–286
(1999). https://doi.org/10.1364/OE.4.000270 Google Scholar
M. Boffety, M. Allain, A. Sentenac, M. Massonneau, and
R. Carminati,
“Analysis of the depth resolution limit of luminescence diffuse optical imaging,”
Opt. Lett., 33
(20), 2290
–2292
(2008). https://doi.org/10.1364/OL.33.002290 Google Scholar
A. Sentenac, C. Gu´erin, P. C. Chaumet, F. Drsek, H. Giovannini, N. Bertaux, and
M. Holschneider,
“Influence of multiple scattering on the resolution of an imaging system: A Cramer-Rao analysis,”
Opt. Express, 15
(3), 1340
–1347
(2007). https://doi.org/10.1364/OE.15.001340 Google Scholar
S. M. Kay,
“Cramer–Rao lower bound,”
Fundamentals of Statistical Signal Processing: Estimation Theory, 27
–81 Prentice Hall, Englewood Cliffs, NJ
(1993). Google Scholar
W. Mo and
N. Chen,
“Design of an advanced time-domain diffuse optical tomography system,”
IEEE J. Sel. Top. Quant., 16
(3), 581
–587
(2010). https://doi.org/http://ieeexplore.ieee.org/Xplore/login.jsp?url=http://ieeexplore.ieee.org/iel5/2944/4481213/05280376.pdf%3Farnumber%3D5280376&authDecision=-203 Google Scholar
W. Mo, T. S. S. Chan, L. Chen, and
N. Chen,
“Quantitative characterization of optical and physiological parameters in normal breasts using time-resolved spectroscopy: in vivo results of 19 Singapore women,”
J. Biomed. Opt., 14
(6), 064004
(2009). https://doi.org/10.1117/1.3257251 Google Scholar
H. Zhao, Y. Tanikawa, Y. Yamada, and
Y. Gao,
“A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomography,”
Opt. Express, 14
(16), 7109
–7124
(2006). https://doi.org/10.1364/OE.14.007109 Google Scholar
H. Zhao, F. Gao, Y. Tanikawa, K. Homma, and
Y. Yamada,
“Time-resolved optical tomographic imaging for the provision of both anatomical and functional information about biological tissue,”
Appl. Opt., 43 1905
–1916
(2005). https://doi.org/10.1364/AO.44.001905 Google Scholar
S. R. Arridge and
J. C. Schotland,
“Optical tomography: Forward and inverse problems,”
Inv. Prob., 25
(12), 123010
(2009). https://doi.org/10.1088/0266-5611/25/12/123010 Google Scholar
|