Optimization of source and detector configurations based on Cramer-Rao lower bound analysis

Ling Chen; Nanguang Chen

doi:10.1117/1.3549738

1 March 2011 Optimization of source and detector configurations based on Cramer-Rao lower bound analysis

Ling Chen, Nanguang Chen

Author Affiliations +

Journal of Biomedical Optics, Vol. 16, Issue 3, 035001 (March 2011). https://doi.org/10.1117/1.3549738

Abstract

Optimization of source and detector (SD) arrangements in a diffuse optical tomography system is helpful for improving measurements' sensitivity to localized changes in imaging domain and enhancing the capacity of noise resistance. We introduced a rigorous and computationally efficient methodology and adapt it into the diffuse optics field to realize the optimizations of SD arrangements. Our method is based on Cramer-Rao lower bound analysis, which combines the diffusion-forward model and a noise model together. This method can be used to investigate the performance of the SD arrangements through quantitative estimations of lower bounds of the standard variances of the reconstructed perturbation depths and values. More importantly, it provides direct estimations of parameters without solving the inverse problem. Simulations are conducted in the reflection geometry to validate the effectiveness of the method on selections of the optimized SD sets, with a fixed number of sources and detectors, from an SD group on a planar probe surface. The impacts of different noise levels and target perturbation depths are considered in the simulations. It is demonstrated that the SD sets selected by this method afford better reconstructed images. This methodology can be adapted to other probe surfaces and other imaging geometries.

1. Introduction

Diffuse 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 bound²⁵ (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. Theory

2.1.

Forward formula

For a given Laplace parameter p, the photon density at location r can be given by the first Born approximation,⁶ 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}

⟨ Φ (r, p) ⟩ = ⟨ Φ_{bg} (r, p) ⟩ + Φ_{pert} (r, p) .

where 〈Φ_bg(r,p)〉 is the background photon density of a homogenous tissue at the location r, and Φ_pert(r,p) is the scattering field caused by the perturbations in optical properties of the tissue. Assuming absorption variations only,

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}

Φ_{pert} (r, p) = \int_{V} G_{0} (r, r_{d}, p {) δ μ}_{a} (r) G_{0} (r_{s}, r, p) d^{3} r .

Rewritten Eq. 2 into a matrix form,

Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\bm \Phi }_{{\rm pert}} = {\bm {\rm A \delta}}. \end{equation}\end{document}

Φ_{pert} = A δ .

[TeX:] ${\bm \Phi }_{{\rm pert}} \!\!=\!\! \left[ {{\rm \Phi }_{{\rm pert1}} {\rm (}{\bf r}{\rm, }p{\rm),\Phi }_{{\rm pert2}} {\rm (}{\bf r}{\rm, }p{\rm),\Phi }_{{\rm pert3}} {\rm (}{\bf r}{\rm, }p{\rm),} \ldots {\rm, \Phi }_{{\rm pert}M} {\rm (}{\bf r}{\rm, }p{\rm)}} \right]^T$

Φ_{pert} = {[Φ_{pert 1} (r, p {), Φ}_{pert 2} (r, p {), Φ}_{pert 3} (r, p), ... {, Φ}_{pert M} (r, p)]}^{T}

is a M × 1 vector representing M measurements. Each element of [TeX:] ${\bm {\rm A}}$

A

is A _{i, j} = G ₀(r _j,r _di,p)G ₀(r _si,r _j,p)dv, where r _si is the source position of i’th measurement, r _di is the detector position of i’th measurement, r _j is the position of the j’th voxel, and dv is the voxel size. [TeX:] ${\bm \delta}$

δ

= [δμ_a1,δμ_a2,δμ_a3, …,δμ_aN]^T denotes the perturbations of the absorption coefficients at N voxel positions. Many works have been conducted to develop inverse methods for estimating [TeX:] ${\bm \delta}$

δ

within the tissue.^{6, 30}

2.2.

Cramér–Rao Lower Bound

While 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}

{σ (Δ μ}_{a}) \geq \sqrt{{CRLB}_{{Δ μ}_{a}}}, σ (R) \geq \sqrt{{CRLB}_{R}},

where [TeX:] $\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }$

\sqrt{{CRLB}_{{Δ μ}_{a}}}

and [TeX:] $\sqrt {{\rm CRLB}_R }$

\sqrt{{CRLB}_{R}}

are the lower bounds of the deviations σ(Δμ_a) and σ(R).

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}$ $Φ$ = [Φ₁,Φ₂,Φ₃, …,Φ_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}$ $θ$ )². The probability density function is written as²⁵

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}

\begin{matrix} P (Φ) & = & \frac{1}{{(2 π)}^{M / 2} \det^{1 / 2} C (θ)} \\ \times \exp [- \frac{1}{2} {[Φ - ⟨ Φ (θ) ⟩]}^{T} C {(θ)}^{- 1} [Φ - ⟨ Φ (θ) ⟩]], \end{matrix}

where M is the number of measurements in one measurement set. σ_k( [TeX:] ${\bm \theta}$

θ

) can be represented as

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}

σ_{k} (θ) = a ⟨ Φ_{k} (θ) ⟩,

where a denotes the noise level.

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}

Φ_{pert} = A_{1} \cdot {Δ μ}_{a} .

[TeX:] ${\bm {\rm A}}_1$

A_{1}

is a M × 1 vector with each element of it formulized as [TeX:] $\sum\nolimits_{n \subset N_1 } {A_{k,n} }$

\sum_{n \subset N_{1}} A_{k, n}

, where n denotes the voxel index number, N ₁ is the set of voxel index numbers of a single target with the size of a few 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.

To calculate the lower bounds [TeX:] $\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }$ $\sqrt{{CRLB}_{{Δ μ}_{a}}}$ and [TeX:] $\sqrt {{\rm CRLB}_R }$ $\sqrt{{CRLB}_{R}}$ for Gaussian observations, the Fisher information matrix (FIM) is introduced with typical elements,²⁵

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}

\begin{matrix} F_{i, j} & \overset{def}{=} & 〈(\frac{\partial}{\partial θ_{i}} \ln P (Φ)) (\frac{\partial}{\partial θ_{j}} \ln P (Φ)) | θ〉 \\ = & {[\frac{\partial ⟨ Φ (θ) ⟩}{\partial θ_{i}}]}^{T} C {(θ)}^{- 1} [\frac{\partial ⟨ Φ (θ) ⟩}{\partial θ_{j}}] \\ + \frac{1}{2} tr [C {(θ)}^{- 1} \frac{\partial C (θ)}{\partial θ_{i}} C {(θ)}^{- 1} \frac{\partial C (θ)}{\partial θ_{j}}], \end{matrix}

where 〈 · · · 〉 denotes the ensemble average. Substitute Eqs. 1, 6, 7 into Eq. 8, we obtain the FIM adapted to our problem,

[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}

F = [\begin{matrix} F_{11}, F_{12} \\ F_{21}, F_{22} \end{matrix}]

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}

\begin{matrix} F_{1, 1} & = & \sum_{k = 1}^{X_{1}, ..., X_{n}} (2 + \frac{1}{a^{2}}) \frac{{(\sum_{n \subset N_{1}} A_{k, n})}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{2}}, \\ F_{1, 2} = F_{2, 1} & = & \sum_{k = 1}^{M} (2 + \frac{1}{a^{2}}) \frac{Δ μ_{a} (\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}{{⟨ Φ_{k} (θ) ⟩}^{2}}, \\ F_{2, 2} & = & \sum_{k = 1}^{M} (2 + \frac{1}{a^{2}}) \frac{{Δ μ}_{a}^{2} {(\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{2}}, \end{matrix}

where ∂_R represents ∂/∂R, the subscript n in Eq. 9 denotes the n’th voxel. For multiparameter case, the Cramér–Rao lower bound states that the covariance matrix of [TeX:] ${\bm \theta}$

θ

satisfies [TeX:] ${\rm Cov}({\bm \theta })\break \ge {\bf F}^{ - 1}$

Cov (θ) \geq F^{- 1}

.²⁵ Thus, [TeX:] $\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }$

\sqrt{{CRLB}_{{Δ μ}_{a}}}

and [TeX:] $\sqrt {{\rm CRLB}_{R} }$

\sqrt{{CRLB}_{R}}

in Eq. 4 can be written as

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}

\begin{matrix} \sqrt{{CRLB}_{{Δ μ}_{a}}} \overset{def}{=} \sqrt{\frac{F_{2, 2}}{\det (F)}} \\ = \sqrt{{(2 + \frac{1}{a^{2}})}^{- 1} \frac{\sum_{k = 1}^{M} [{(\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}^{2} / {⟨ Φ_{k} (θ) ⟩}^{2}]}{\sum_{\{k, j\} \subset C} \frac{{[(\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{j, n}) - (\sum_{n \subset N_{1}} A_{j, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})]}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{2} {⟨ Φ_{j} (θ) ⟩}^{2}}}} . \end{matrix}

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}

\begin{matrix} \sqrt{{CRLB}_{R}} \overset{def}{=} \sqrt{\frac{F_{1, 1}}{\det (F)}} \\ = \sqrt{{(2 + \frac{1}{a^{2}})}^{- 1} \frac{\sum_{k = 1}^{M} [{(\sum_{n \subset N_{1}} A_{k, n})}^{2} / {⟨ Φ_{k} (θ) ⟩}^{2}]}{{Δ μ}_{a}^{2} \sum_{\{k, j\} \subset C} \frac{{[(\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{j, n}) - (\sum_{n \subset N_{1}} A_{j, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})]}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{2} {⟨ Φ_{j} (θ) ⟩}^{2}}}} . \end{matrix}

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. Recall that M is the total number of measurements from one SD set. Detailed derivations of Eqs. 10, 11 can be found in the Appendix.

[TeX:] $\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }$ $\sqrt{{CRLB}_{{Δ μ}_{a}}}$ represents the cross-coupling precision limit for Δμ_a when taking the influence of unknown R into account. [TeX:] $F_{11}^{ - 1/2}$ $F_{11}^{- 1 / 2}$ is the precision limit for Δμ_a when R is known. Same rules apply for [TeX:] $\sqrt {{\rm CRLB}_{R} }$ $\sqrt{{CRLB}_{R}}$ and [TeX:] $F_{22}^{ - 1/2}$ $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 ⁻¹.

3. Simulation Results and Discussion

We 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 works^{14, 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×10⁹ S⁻¹ to 0.6×10⁹ S⁻¹ in step size of 0.2×10⁹ S⁻¹. The absorption coefficient and reduced scattering coefficient of the homogenous background were set to be 0.02 and 6 cm⁻¹, respectively, which were close to the optical properties for normal breast tissue.¹¹ The refractive index of the background was 1.4, which was close to the tissue refractive index.²⁷

Fig. 1

XY plane of rotation probe: the circles on the left indicate the 23 source positions for selection and the asterisks on the right represent the 12 detector positions for selection. (Color online only.)

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⁻¹ and reduced scattering coefficient equal to 6 cm⁻¹ 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.

Fig. 2

Normalized precision limits of perturbation depth and perturbation value for all examined SD sets. Each point of data in (a–d) represents the normalized precision limit for one SD set under specified conditions. The normalization formulas and the ranges of the precision limits are indicated in (a–d). (a) The normalized precision limits of perturbation depth with target depth = 1 cm and noise level equal to 1, 2, and 3%, respectively. (b) Normalized precision limits of perturbation depth with noise level = 1% and target depth equal to 1.0, 1.8, and 2.6 cm, respectively. (c) Normalized precision limits of perturbation value with target depth = 1 cm and noise level equal to 1, 2, and 3%, respectively. (d) Normalized precision limits of perturbation value with noise level = 1% and target depth equal to 1.0, 1.8, and 2.6 cm, respectively.

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 algorithm³⁰ 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.

Fig. 3

Sample variances of the reconstructed target perturbation values versus the corresponding CRLB values for the selected 1500 SD sets. Each black asterisk represents the sample variance of the reconstructed target perturbation value generated from each SD set after 20 times’ repetitions of simulations. Each red asterisk represents the CRLB value of the same SD set. The data unit for all the data presented in the images is centimeters to the –1. The images are represented in columns based on noise level from 1 to 3%, and in rows based on target depth equal to 1.0, 1.8, and 2.6 cm, respectively. The target absorption perturbation value is 0.2 cm⁻¹. (Color online only.)

Fig. 4

Sample variances of the reconstructed target depths versus the corresponding CRLB values for the selected 1500 SD sets. Each black asterisk represents the sample variance of the reconstructed target depth generated from each SD set after 20 times’ repetitions of simulations. Each red asterisk represents the CRLB value of the same SD set. The data unit for all the data presented in the images is centimeters. The figures are represented in columns based on noise level from 1 to 3%, and in rows based on target depth equal to 1.0, 1.8, and 2.6 cm, respectively. The target absorption perturbation value is 0.2 cm⁻¹. (Color online only.)

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⁻⁴ 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} } }$ $\sqrt{{CRLB}_{{Δ μ}_{a}}}$ and [TeX:] $\sqrt {{\rm CRLB}_R }$ $\sqrt{{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.

Fig. 5

(a) The total number of the useful singular values above 10⁻⁴ for the selected 1500 SD sets. Each asterisk represents the total number of the useful singular values for each SD set. (b) The sum of the useful singular values above 10⁻⁴ for the selected 1500 SD sets. Each asterisk represents the sum of the useful singular values for each SD set.

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⁻¹ 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.

Fig. 6

Six SD sets with different precision limits of the perturbation depths and perturbation values. The solid circles from the 23 source positions indicate the chosen source positions while the solid circles from the 12 detector positions indicate the chosen detector positions for three sources and four detectors’ case. (a–c) represent three of the best SD sets with the lowest precision limits. (d–f) represent three of the worst SD sets with the highest precision limits. (Color online only.)

Fig. 7

Reconstruction images of targets in the YZ plane (X = 0): Target depth is (a,b) 1.0 cm, (c,d) 1.8 cm, and (e,f) 2.6 cm. (a,c,e) using SD set 1. (b,d,f) using SD set 6. Images in the first column show the original depths and values of the target. Images in the second column are reconstructed from a signal without noise. Images in the third column are reconstructed from a signal with 1% noise. Images in the fourth column are reconstructed from a signal with 2% noise. Images in the fifth column are reconstructed from a signal with 3% noise.

Table 1

Reconstruction values from Fig. 7.

		Noise level
		0%		1%		2%		3%
Target depth (cm)	SD set	μbg/cm	μtg/cm	μbg/cm	μtg/cm	μbg/cm	μtg/cm	μbg/cm	μtg/cm
1.0	SD set #1	0.0289±0.0246	0.1957±0.0352	0.0289±0.0244	0.1974±0.0482	0.0305±0.0250	0.1882±0.0514	0.0319±0.0317	0.1887±0.0511
	SD set #6	0.0326±0.0277	0.1924±0.0420	0.0324±0.0272	0.1956±0.0394	0.0338±0.0299	0.1863±0.0336	0.0349±0.0238	0.1482±0.0340
1.8	SD set #1	0.0343±0.0305	0.1710±0.0251	0.0342±0.0318	0.1949±0.0307	0.0402±0.0534	0.2138±0.0420	0.0436±0.0582	0.2167±0.0629
	SD set #6	0.0357±0.0296	0.1574±0.0265	0.0380±0.0315	0.1585±0.0309	0.0435±0.0355	0.1200±0.0421	0.0489±0.0457	0.1257±0.0291
2.6	SD set #1	0.0325±0.0301	0.1636±0.0214	0.0409±0.0408	0.1733±0.0244	0.0499±0.0458	0.1610±0.0139	0.0430±0.0391	0.1680±0.0196
	SD set #6	0.0340±0.0320	0.1517±0.0177	0.0396±0.0353	0.1458±0.0175	0.0427±0.0378	0.1440±0.0130	0.0467±0.0418	0.1520±0.0133
True values		0.22	0.02	0.22	0.02	0.22	0.02	0.22	0.22

4. Conclusion

We 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.

Appendices

Appendix

The 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}$ $Φ_{pert} = A_{1} \cdot {Δ μ}_{a}$ . [TeX:] ${\bm {\rm A}}_1$ $A_{1}$ is a M × 1 vector with each element of it formulized as [TeX:] $\sum\nolimits_{n \subset N_1 } {A_{k,n} }$ $\sum_{n \subset N_{1}} A_{k, n}$ , where n denotes the voxel index number, N ₁ 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})$ $C (θ)$ 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}

\frac{\partial ⟨ Φ (θ) ⟩}{\partial {Δ μ}_{a}} = \frac{\partial Φ_{pert} (θ)}{\partial {Δ μ}_{a}} = A_{1};

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}

\frac{\partial ⟨ Φ (θ) ⟩}{\partial R} = \frac{\partial Φ_{pert} (θ)}{\partial R} = {Δ μ}_{a} \cdot \partial_{R} A_{1};

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}

\frac{\partial C (θ)}{\partial {Δ μ}_{a}} = 2 a^{2} \cdot diag {[⟨ Φ (θ) ⟩ \cdot A_{1}]}^{T};

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}

\frac{\partial C (θ)}{\partial R} = 2 a^{2} {Δ μ}_{a} \cdot diag {[⟨ Φ (θ) ⟩ \cdot \partial_{R} A_{1}]}^{T} .

Using Eqs. 12 13 14 15 here, we get

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}

\begin{matrix} {[\frac{\partial ⟨ Φ (θ) ⟩}{\partial {Δ μ}_{a}}]}^{T} C {(θ)}^{- 1} [\frac{\partial ⟨ Φ (θ) ⟩}{\partial R}] = {Δ μ}_{a} \cdot A_{1}^{T} C {(θ)}^{- 1} \partial_{R} A_{1} \\ = \sum_{k = 1}^{M} \frac{{Δ μ}_{a} (\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}{a^{2} {⟨ Φ_{k} (θ) ⟩}^{2}}; \end{matrix}

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}

\begin{matrix} \frac{1}{2} tr [C {(θ)}^{- 1} \frac{\partial C (θ)}{\partial {Δ μ}_{a}} C {(θ)}^{- 1} \frac{\partial C (θ)}{\partial R}] \\ = \sum_{k = 1}^{M} \frac{2 {Δ μ}_{a} (\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}{{⟨ Φ_{k} (θ) ⟩}^{2}} . \end{matrix}

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}

\begin{matrix} F_{1, 2} = F_{2, 1} = \sum_{k = 1}^{M} \frac{{Δ μ}_{a} (\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}{a^{2} {⟨ Φ_{k} (θ) ⟩}^{2}} \\ + \sum_{k = 1}^{M} \frac{2 {Δ μ}_{a} (\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}{{⟨ Φ_{k} (θ) ⟩}^{2}} \\ = \sum_{k = 1}^{M} (2 + \frac{1}{a^{2}}) \frac{{Δ μ}_{a} (\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}{{⟨ Φ_{k} (θ) ⟩}^{2}} . \end{matrix}

The other elements in the Fisher matrix are also derived using Eqs. 12 13 14 15. Thus, we obtain Eq. 9.

Using Eqs. 7, 9, [TeX:] $\sqrt {{\rm CRLB}_{{\rm \Delta \mu }_{\rm a} } }$ $\sqrt{{CRLB}_{{Δ μ}_{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}

\sqrt{{CRLB}_{{Δ μ}_{a}}} \overset{def}{=} \sqrt{\frac{F_{2, 2}}{\det (F)}} = \sqrt{\frac{F_{2, 2}}{F_{1, 1} F_{2, 2} - F_{1, 2}^{2}}} .

Because*****

[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}

\begin{matrix} F_{1, 1} F_{2, 2} & = & {(2 + \frac{1}{a^{2}})}^{2} {Δ μ}_{a}^{2} [\sum_{k = 1}^{M} \frac{{(\sum_{n \subset N_{1}} A_{k, n})}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{2}}] [\sum_{k = 1}^{M} \frac{{(\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{2}}] \\ = & {(2 + \frac{1}{a^{2}})}^{2} {Δ μ}_{a}^{2} \{\begin{matrix} \sum_{k = 1}^{M} \frac{{(\sum_{n \subset N_{1}} A_{k, n})}^{2} {(\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{4}} \\ + \sum_{\{k, j\} \subset C} [\frac{{(\sum_{n \subset N_{1}} A_{k, n})}^{2} {(\partial_{R} \sum_{n \subset N_{1}} A_{j, n})}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{2} {⟨ Φ_{j} (θ) ⟩}^{2}} + \frac{{(\sum_{n \subset N_{1}} A_{j, n})}^{2} {(\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}^{2}}{{⟨ Φ_{j} (θ) ⟩}^{2} {⟨ Φ_{k} (θ) ⟩}^{2}}] \end{matrix}\} \end{matrix}

[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}

\begin{matrix} F_{1, 2}^{2} & = & {(2 + \frac{1}{a^{2}})}^{2} {Δ μ}_{a}^{2} {[\sum_{k = 1}^{M} \frac{(\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}{{⟨ Φ_{k} (θ) ⟩}^{2}}]}^{2} \\ = & {(2 + \frac{1}{a^{2}})}^{2} {Δ μ}_{a}^{2} \{\begin{matrix} \sum_{k = 1}^{M} \frac{{(\sum_{n \subset N_{1}} A_{k, n})}^{2} {(\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{4}} \\ + 2 \sum_{\{k, j\} \subset C} [\frac{(\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n}) (\sum_{n \subset N_{1}} A_{j, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{j, n})}{{⟨ Φ_{k} (θ) ⟩}^{2} {⟨ Φ_{j} (θ) ⟩}^{2}}] \end{matrix}\}, \end{matrix}

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}

\begin{matrix} F_{1, 1} F_{2, 2} - F_{1, 2}^{2} & = & {(2 + \frac{1}{a^{2}})}^{2} {Δ μ}_{a}^{2} \sum_{\{k, j\} \subset C} [\begin{matrix} \frac{{(\sum_{n \subset N_{1}} A_{k, n})}^{2} {(\partial_{R} \sum_{n \subset N_{1}} A_{j, n})}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{2} {⟨ Φ_{j} (θ) ⟩}^{2}} + \frac{{(\sum_{n \subset N_{1}} A_{j, n})}^{2} {(\partial_{R} \sum_{n \subset N_{1}} A_{k, n})}^{2}}{{⟨ Φ_{j} (θ) ⟩}^{2} {⟨ Φ_{k} (θ) ⟩}^{2}} \\ - 2 \frac{(\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n}) (\sum_{n \subset N_{1}} A_{j, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{j, n})}{{⟨ Φ_{k} (θ) ⟩}^{2} {⟨ Φ_{j} (θ) ⟩}^{2}} \end{matrix}] \\ = & {(2 + \frac{1}{a^{2}})}^{2} {Δ μ}_{a}^{2} \sum_{\{k, j\} \subset C} \frac{{[(\sum_{n \subset N_{1}} A_{k, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{j, n}) - (\sum_{n \subset N_{1}} A_{j, n}) (\partial_{R} \sum_{n \subset N_{1}} A_{k, n})]}^{2}}{{⟨ Φ_{k} (θ) ⟩}^{2} {⟨ Φ_{j} (θ) ⟩}^{2}}, \end{matrix}

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.

References

1.

V. 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

2.

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

3.

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

4.

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

5.

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

6.

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

7.

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

8.

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

9.

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

10.

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

11.

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

12.

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

13.

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

14.

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

15.

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

16.

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

17.

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

18.

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

19.

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

20.

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

21.

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

22.

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

23.

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

24.

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

25.

S. M. Kay, “Cramer–Rao lower bound,” Fundamentals of Statistical Signal Processing: Estimation Theory, 27 –81 Prentice Hall, Englewood Cliffs, NJ (1993). Google Scholar

26.

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

27.

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

28.

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

29.

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

30.

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

Citation Download Citation

Ling Chen and Nanguang Chen "Optimization of source and detector configurations based on Cramer-Rao lower bound analysis," Journal of Biomedical Optics 16(3), 035001 (1 March 2011). https://doi.org/10.1117/1.3549738

Published: 1 March 2011

Access the abstract

JOURNAL ARTICLE
11 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 11 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Sensors

Absorption

Diffuse optical imaging

Image quality

Imaging systems

Computer simulations

Target detection

1.

Introduction

2.

Theory