Translator Disclaimer
24 October 2016 Sampling scheme optimization for diffuse optical tomography based on data and image space rankings
Author Affiliations +
Abstract
We present a methodology for the optimization of sampling schemes in diffuse optical tomography (DOT). The proposed method exploits singular value decomposition (SVD) of the sensitivity matrix, or weight matrix, in DOT. Two mathematical metrics are introduced to assess and determine the optimum source–detector measurement configuration in terms of data correlation and image space resolution. The key idea of the work is to weight each data measurement, or rows in the sensitivity matrix, and similarly to weight each unknown image basis, or columns in the sensitivity matrix, according to their contribution to the rank of the sensitivity matrix, respectively. The proposed metrics offer a perspective on the data sampling and provide an efficient way of optimizing the sampling schemes in DOT. We evaluated various acquisition geometries often used in DOT by use of the proposed metrics. By iteratively selecting an optimal sparse set of data measurements, we showed that one can design a DOT scanning protocol that provides essentially the same image quality at a much reduced sampling.

1.

Introduction

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

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

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

2.

Methods

2.1.

Review of Diffuse Optical Tomography Forward and Inverse Problem

For a given distribution of optical parameters, absorption coefficient μa(r), and reduced scattering coefficient μs(r)=(1g)μs(r), the light propagation in biological tissues is usually modeled by

Eq. (1)

D(r)ϕ(r,ω)+[μa(r)+(iωcm(r))ϕ(r,ω)]ϕ(r,ω)=q0(r,ω),
where D(r)=13[μa(r)+μs(r)] is the diffusion coefficient, cm(r) is the speed of the light in the medium, q0(r,ω) is the input source term, and the ϕ(r,ω) is the generated photon fluence rate at position “r” with a modulation frequency “ω.”2,3,19 The forward model data ϕi,jc are generated by sampling in the data space with the ith source and the jth detector positions and is represented by

Eq. (2)

ϕi,jc=ϕi,jc(μa,D).
In an inverse problem, for a given experimental data ϕi,jm, we wish to find optical parameters μ={μa,D} (for both absorption and diffusion coefficients) by minimizing the objective function as follows:

Eq. (3)

O(ϕi,jm,μ)=argminμ(ϕi,jmϕi,jc)2.
The objective function O(ϕi,jm,μ) measures the compatibility between the estimated optical parameters μ={μa,D} and the experimentally measured data ϕi,jm. Let μ={μa,D} be an initial estimate of the optical parameters of the perfect model prior to inversion, and ϕi,jc(μ) be the corresponding data. A linearization of this optimization problem can be formulated by the use of a perturbation approach. Consider perturbing the prior initial optical parameters μ by a small amount Δμ.20 Equation (3) can be rewritten as Eq. (4) with δ=ϕi,jm(μ)ϕi,jc(μ) representing the compatibility in the initial estimate

Eq. (4)

O(ϕi,jm,μ)=argminμ[δ+ϕi,jc(μ)ϕi,jc(μ+Δμ)]2.
The linearization step is carried out by expanding ϕi,jc(μ+Δμ) through a Taylor series and neglecting higher order terms

Eq. (5)

ϕi,jc(μ+Δμ)=ϕi,jc(μ)+JΔμ+,
where Jmn=ϕc(μ)mμn|μ m={1,2,, number of SD pair measurements}, n={1,2,, number of unknowns} is the sensitivity or weight matrix and has four parts
Jmn=[J1=lnIDJ2=lnIμaJ3=ϕDJ4=ϕμa],
Each row in the matrix “J,” also known as sensitivity maps, relates the change in the boundary data either by Born approximation [amplitude (I) and phase] or by Rytov approximation [ln(I) and phase]3 to the small perturbation in the optical parameters (μa, D) of the reconstructed image basis (voxels or nodes). In the formulation of data and ISRs, which will be described later, only the J2 part of the Jacobean matrix is considered without loss of generality. Substituting Eq. (5) into Eq. (4), the linearized objective function can be written by

Eq. (6)

O(ϕi,jm,μ)argminμδJΔμ2.
The minimizer of Eq. (6) is given by

Eq. (7)

Δμ=(JTJ)1JTδ,

Eq. (8)

Δμ=(JTJ+αI)1JTδ.
Usually, in DOT, the square matrix JTJ in (7) turns out to be ill-conditioned due to the presence of extremely small singular values; one way to stabilize the inversion in such cases is to add the diagonal term as shown in (8). We exploited an adjoint formulation to calculate the sensitivity maps, which is known to be computationally robust compared to direct methods.21 The initial guess on the bulk optical properties for calculating the sensitivity matrix is obtained by a calibration procedure22 in the imaging domain whose analytical solutions are readily available.

2.2.

Singular Value Decomposition

A technique which is of particular interest for investigating ill-conditioned and/or rank-deficient problems is the SVD. In SVD, a matrix is factorized by a set of three matrices

Eq. (9)

J=USVT,
where U and V are the orthonormal matrices and S is a diagonal matrix of singular values of J. The columns in U and V are known as modes that span two different spaces: data space and image or model space, respectively. The magnitude of the singular values provides significance of the corresponding data and image space modes.

2.3.

Metrics for DOT Geometry Optimization

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

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

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

Eq. (10)

VVT=I,

Eq. (11)

UUT=I,

Eq. (12)

VTJTJV=STS=S2,

Eq. (13)

UTJJTU=SST=S2.

2.3.1.

Data space ranking

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

Eq. (14)

G=[JV][JV],
where the symbol is used to represent the Hadamard product. One can similarly rewrite Eq. (14) as follows:

Eq. (15)

G=[US][US].
Each column of “G” sums to the corresponding squared singular value of “J.” Thus, the ith column represents the contribution of the ith SD measurement pair to the associated squared singular value. If the right-hand side of Eq. (15) is postmultiplied by a matrix of inverse of the squared singular values, each direction in the n-dimensional space will have equal importance.18 However, because of the presence of smaller singular values, it is not straightforward to get the inversion.25 One way to get the inversion is to add a small regularization factor λ2 as follows:

Eq. (16)

FSS=([US][US])×(S2+λ2I)1.
By using the Hadamard product properties of a multiplication of diagonal matrix, Eq. (16) can be rewritten by

Eq. (17)

FSS=[UU]×[S2×(S2+λ2I)1],
where the term in the right bracket is known as the filter factor in Tikhonov (Tikh) regularization (Ftikh) sense24

Eq. (18)

Ftikh=S2×(S2+λ2I)1Ftikh,mn={sm2sm2+λm2,m=n0,mn.

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

Eq. (19)

Ftsvd,ij={1,si>ϵ0,si<ϵ.
where “ϵ” corresponds to the truncate on singular values. The ith term in the jth column of the matrix FSS (referred to as the matrix of fractional squared singular values) represents the fractional contribution of the ith SD measurement to the jth squared singular value. Summation over each row in Eq. (17) would yield the contribution of the corresponding SD pair to the rank of the sensitivity matrix, and it is referred to as data space ranking (DSR) in this work. Note that the DSR can range between zero and one {0, 1}

Eq. (20)

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

2.3.2.

Image space ranking

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

Eq. (21)

G=[JTU][JTU].
Rewriting Eq. (21), one can get

Eq. (22)

G=[VST][VST].

One can similarly introduce the matrix of fractional squared singular values

Eq. (23)

FSS=[VV]×[S2×(S2+λ2I)1].
Summation over each row of Eq. (23) would yield the contribution of the image basis (voxels or nodes) to the rank of the sensitivity matrix, and it is referred to as ISR. The ISR ranges between zero and one, {0, 1}

Eq. (24)

ISR=jVij×FjjVij×Fjj.

2.3.3.

Sparse selection of source–detector pairs

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

Fig. 1

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

JBO_21_10_106004_f001.png

2.4.

Numerical Study

We used two-dimensional (2-D) imaging geometries for simplicity in this work, but it can be extended to a three-dimensional (3-D) problem straightforwardly. In the case of breast DOT imaging, two most commonly considered imaging geometries include: a mammographically compressed breast with the sources being placed at one plane and the detectors on the other,29 or a circular distribution of the sources and the detectors around the breast volume.8 We employed both the rectangular imaging geometry and the circular imaging geometry in this work as summarized in Table 1. For a circular geometry, the sources and the detectors were placed around the object support in an equiangular fashion.

Table 1.

Description of geometries for which SD optimization analysis was performed.

TypesNo. of sourcesNo. of detectorsNo. of projections
Rectangular geometryReflectance1515210
Sparse transmission1515225
First dense transmission1535525
Second dense transmission15701050
Circular geometryTomographic1616240
Mouse head geometryTomographic8856

In all the simulation cases, we fixed the optical properties by: absorption coefficient μa=0.01  mm1, scattering coefficient μs=1  mm1, and refractive index=1.3 for the homogenous background, and absorption coefficient μa=0.02  mm1, scattering coefficient μs=2  mm1, and refractive index=1.3 for the tumor anomalies. One percent (1%) of Gaussian noise in the amplitude and in the phase angle was added for all the cases unless otherwise specified. The simulations were conducted assuming a Gaussian source intensity profile having the full-width-half-maximum of 2 mm.

We adopted a nonlinear iterative algorithm under the Rytov approximation and a dual-mesh scheme to avoid inverse crime for image reconstruction in all the simulation cases. Meshes were constructed as summarized in Table 2.

Table 2

Description of FEM meshes constructed for numerical study.

FEM triangular meshCircular geometryRectangular geometryMouse head geometry
One anomalyTwo anomaly
Coarse meshNodes17853321668
Elements341864001265
Fine meshNodes4882706545602187
Elements951713,83388364206

3.

Results and Discussion

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

Fig. 2

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

JBO_21_10_106004_f002.png

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

Fig. 3

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

JBO_21_10_106004_f003.png

Fig. 4

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

JBO_21_10_106004_f004.png

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

Fig. 5

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

JBO_21_10_106004_f005.png

Fig. 6

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

JBO_21_10_106004_f006.png

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

Fig. 7

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

JBO_21_10_106004_f007.png

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

Fig. 8

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

JBO_21_10_106004_f008.png

Fig. 9

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

JBO_21_10_106004_f009.png

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

Fig. 10

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

JBO_21_10_106004_f010.png

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

Fig. 11

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

JBO_21_10_106004_f011.png

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

Fig. 12

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

JBO_21_10_106004_f012.png

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

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

Fig. 13

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

JBO_21_10_106004_f013.png

Fig. 14

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

JBO_21_10_106004_f014.png

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

Fig. 15

Trace of image space rank at various number of measurements.

JBO_21_10_106004_f015.png

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

Eq. (25)

RMSE=1Ni=1N(μiaμir)2,
where N is the number of nodes or voxels within the anomaly region, and μia and μir are the actual and reconstructed optical coefficients (absorption or scattering coefficients), respectively. Table 3 summarizes the results of the previous simulation studies.

Table 3

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

Sampling schemeNMOptimizedRandom
CaseNoise (%)μaμs′μaμs′
Circular (1)02400.438540.9904
01200.446641.66830.483945.2165
0840.480548.29610.576354.5054
0480.585456.54650.771768.0097
12400.454537.7276
11200.471741.79400.520045.4893
1840.470946.95540.574549.5964
1480.589258.42210.736061.3935
22400.518038.4695
21200.645545.89010.801545.1149
2840.657751.98960.714160.1785
2480.683366.54850.770067.1893
52400.572057.0065
51200. 590558.82070.641559.0057
5840.660257.87080.620369.1217
5480.746268.93020.730177.0911
Circular (2)12400.477540.0851
11200.478845.25530.548446.5737
1840.522346.27480.619056.4978
1480.608257.66130.737966.0918
Rectangular12250.555046.2080
11120.526446.45140.624348.8919
1800.651853.66060.714054.2063
1450.675164.91960.741966.4600

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

Fig. 16

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

JBO_21_10_106004_f016.png

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

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

Table 4

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

Full (56)Optimized (28)Random (28)ROI optimized (28)
Mouse brain0.56320.66970.70390.6123

4.

Conclusions

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

Acknowledgments

This work was supported in part by the NRF Grant No. 2013M2A2A9043476 and the R&D Convergence Program of NST (National Research Council of Science and Technology) of Korea (Grant No. CAP-13-3-KERI).

References

1. 

Y. Yamada and S. Okawa, “Diffuse optical tomography: present status and its future,” Opt. Rev., 21 (3), 185 –205 (2014). http://dx.doi.org/10.1007/s10043-014-0028-7 Google Scholar

2. 

H. Dehghani et al., “Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng., 25 711 –732 (2009). http://dx.doi.org/10.1002/cnm.v25:6 Google Scholar

3. 

A. P. Gibson, J. C. Hebden and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol., 50 (4), R1 –R43 (2005). http://dx.doi.org/10.1088/0031-9155/50/4/R01 Google Scholar

4. 

Q. Fang et al., “Combined optical and X-ray tomosynthesis breast imaging,” Radiology, 258 (1), 89 –97 (2011). http://dx.doi.org/10.1148/radiol.10082176 Google Scholar

5. 

H. Dehghani et al., “Multiwavelength three-dimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results,” Appl. Opt., 42 (1), 135 –145 (2003). http://dx.doi.org/10.1364/AO.42.000135 Google Scholar

6. 

F. Leblond, K. M. Tichauer and B. W. Pogue, “Singular value decomposition metrics show limitations of detector design in diffuse fluorescence tomography,” Biomed. Opt. Express, 1 (5), 1514 –1531 (2010). http://dx.doi.org/10.1364/BOE.1.001514 Google Scholar

7. 

J. P. Culver et al., “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys., 30 (2), 235 –247 (2003). http://dx.doi.org/10.1118/1.1534109 Google Scholar

8. 

Y. Zhao et al., “Portable, parallel 9-wavelength near-infrared spectral tomography (NIRST) system for efficient characterization of breast cancer within the clinical oncology infusion suite,” Biomed. Opt. Express, 7 (6), 2186 –2201 (2016). http://dx.doi.org/10.1364/BOE.7.002186 Google Scholar

9. 

J. P. Culver et al., “Optimization of optode arrangements for diffuse optical tomography: a singular-value analysis,” Opt. Lett., 26 (10), 701 –703 (2001). http://dx.doi.org/10.1364/OL.26.000701 Google Scholar

10. 

H. Xu et al., “Near-infrared imaging in the small animal brain: optimization of fiber positions,” J. Biomed. Opt., 8 (1), 102 –110 (2003). http://dx.doi.org/10.1117/1.1528597 Google Scholar

11. 

P. K. Yalavarthy et al., “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). http://dx.doi.org/10.1364/OE.14.006113 Google Scholar

12. 

E. E. Graves et al., “Singular-value analysis and optimization of experimental parameters in fluorescence molecular tomography,” J. Opt. Soc. Am. A., 21 (2), 231 –241 (2004). http://dx.doi.org/10.1364/JOSAA.21.000231 Google Scholar

13. 

V. Pera, D. H. Brooks and M. Niedre, “On the use of the Cramér-Rao lower bound for diffuse optical imaging system design,” J. Biomed. Opt., 19 (2), 025002 (2014). http://dx.doi.org/10.1117/1.JBO.19.2.025002 Google Scholar

14. 

L. Chen and N. Chen, “Optimization of source and detector configurations based on Cramer-Rao lower bound analysis,” J. Biomed. Opt., 16 (3), 035001 (2011). http://dx.doi.org/10.1117/1.3549738 Google Scholar

15. 

D. Karkala and P. K. Yalavarthy, “Data-resolution based optimization of the data-collection strategy for near infrared diffuse optical tomography,” Med. Phys., 39 (8), 4715 –4725 (2012). http://dx.doi.org/10.1118/1.4736820 Google Scholar

16. 

C. B. Shaw and P. K. Yalavarthy, “Incoherence-based optimal selection of independent measurements in diffuse optical tomography measurements in diffuse optical tomography,” J. Biomed. Opt., 19 (3), 036017 (2014). http://dx.doi.org/10.1117/1.JBO.19.3.036017 Google Scholar

17. 

S. Sabir et al., “Adaptive selection of minimally correlated data for optimization of source-detector configuration in diffuse optical tomography,” Proc. SPIE, 9701 97010V (2016). http://dx.doi.org/10.1117/12.2211580 Google Scholar

18. 

D. C. Kammer, “Sensor placement for on-orbit modal identification and correlation of large space structure,” J. Guid. Control, Dyn., 14 (2), 251 –259 (1991). http://dx.doi.org/10.2514/3.20635 Google Scholar

19. 

S. Arridge and J. Schotland, “Optical tomography: forward and inverse problems,” Inverse Problems, 25 (12), 123010 (2009). http://dx.doi.org/10.1088/0266-5611/25/12/123010 Google Scholar

20. 

A. Curtis, “Optimal design of focused experiments and surveys,” Geophys. J. Int., 139 (1), 205 –215 (1999). http://dx.doi.org/10.1046/j.1365-246X.1999.00947.x Google Scholar

21. 

S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part 2: finite-element-method calculations,” Appl. Opt., 34 (34), 8026 –8037 (1995). http://dx.doi.org/10.1364/AO.34.008026 Google Scholar

22. 

B. W. Pogue et al., “Calibration of near-infrared frequency-domain tissue spectroscopy for absolute absorption coefficient quantitation in neonatal head-simulating phantoms,” J. Biomed. Opt., 5 (2), 185 –193 (2000). http://dx.doi.org/10.1117/1.429985 Google Scholar

23. 

L. Horesh, E. Haber and L. Tenorio, “Optimal experimental design for the large-scale nonlinear Ill-posed problem of impedance imaging,” Large-Scale Inverse Problems and Quantification of Uncertainty, 273 –290 John Wiley & Sons, Ltd., Chichester, United Kingdom (2010). Google Scholar

24. 

R. C. Aster, B. Borchers and C. H. Thurber, Parameter Estimation and Inverse Problems, 2nd ed.Elsevier Academic Press, New York (2013). Google Scholar

25. 

A. Jin, B. Yazici and V. Ntziachristos, “Light illumination and detection patterns for fluorescence diffuse optical tomography based on compressive sensing,” IEEE Trans. Image Process., 23 (6), 2609 –2624 (2014). http://dx.doi.org/10.1109/TIP.2014.2300756 Google Scholar

26. 

P. K. Yalavarthy et al., “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express, 15 (13), 8043 –8058 (2007). http://dx.doi.org/10.1364/OE.15.008043 Google Scholar

27. 

Y. Zhao et al., “Optimization of image reconstruction for magnetic resonance imaging-guided near-infrared diffuse optical spectroscopy in breast,” J. Biomed. Opt., 20 (5), 056009 (2015). http://dx.doi.org/10.1117/1.JBO.20.5.056009 Google Scholar

28. 

B. A. Brooksby et al., “Near-infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstsructing heterogeneities,” IEEE J. Sel. Top. Quantum Electron., 9 (2), 199 –209 (2003). http://dx.doi.org/10.1109/JSTQE.2003.813304 Google Scholar

29. 

H. Y. Ban et al., “Heterodyne frequency-domain multispectral diffuse optical tomography of breast cancer in the parallel-plane transmission geometry,” Med. Phys., 43 (7), 4383 –4395 (2016). http://dx.doi.org/10.1118/1.4953830 Google Scholar

30. 

S. Borguet and O. Léonard, “The fisher information matrix as a relevant tool for sensor selection in engine health monitoring,” Int. J. Rotating Mach., 2008 1 –10 (2008). http://dx.doi.org/10.1155/2008/784749 Google Scholar

31. 

O. Lee et al., “Compressive diffuse optical tomography: noniterative exact reconstruction using joint sparsity,” IEEE Trans. Med. Imaging, 30 (5), 1129 –1142 (2011). http://dx.doi.org/10.1109/TMI.2011.2125983 Google Scholar

Biography

Sohail Sabir received his MS degree in nuclear and quantum engineering from KAIST, Republic of Korea, in 2014. He is currently a PhD candidate in the Department of Nuclear and Quantum Engineering at KAIST. His research interests include radiation safety in radiological imaging, optical imaging, and compressive sensing approaches in medical imaging.

Changhwan Kim received his BS and MS degrees in the Department of Nuclear and Quantum Engineering from KAIST, Republic of Korea, in 2013 and 2015, respectively. Currently, his PhD studies focus on the field of imaging physics including diffuse optical tomography and x-ray computed tomography.

Sanghoon Cho received his BS degree in applied physics from Dankook University, Republic of Korea, in 2015. He is a master’s degree candidate in the Nuclear and Quantum Engineering Department, KAIST, Republic of Korea. His current research interests are devoted to the study of signal processing in diffuse optical tomography for diagnostic imaging applications, nuclear medicine for biological and functional imaging, and scatter correction in cone-beam computed tomography imaging.

Duchang Heo received his PhD in electronics engineering from Korea University, Republic of Korea, in 2003. He joined Samsung Electronics Ltd., in 2004, where he was involved in research on broadband light-wave source for a network system and on a picoprojector incorporated in cellular phones. He joined KERI, Republic of Korea, in 2008, and he is conducting research on femtosecond fiber lasers and diffuse optical tomography.

Kee Hyun Kim received his BS degree in electrical and computer engineering from Seoul National University in 2008 and his MS degree in biomedical engineering from Seoul National University in 2010. In 2015, he became a research engineer in KERI, Republic of Korea. His current research interests are devoted to the study of instrumentation and signal processing in diffuse optical tomography for diagnostic imaging applications.

Jong Chul Ye received his BSc and MSc degrees from Seoul National University, Republic of Korea, and his PhD from Purdue University, West Lafayette. He joined KAIST, Republic of Korea, in 2004, where he is currently KAIST endowed chair professor and the professor of the Department of Bio/Brain Engineering. His current research interests include compressed sensing and statistical signal processing for various imaging modalities.

Seungryong Cho received his BS and MS degrees in physics from Seoul National University in 1995 and 1997, respectively, and his PhD in medical physics from the University of Chicago, USA, in 2009. He is an associate professor in the Department of Nuclear and Quantum Engineering, KAIST, Republic of Korea. His research interests include radiation-based imaging for medicine and industry, multimodality imaging, radiation therapy physics, image reconstruction algorithms, and image processing techniques.

© 2016 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2016/$25.00 © 2016 SPIE
Sohail Sabir, Changhwan Kim, Sanghoon Cho, Duchang Heo, Kee Hyun Kim, Jong Chul Ye, and Seungryong Cho "Sampling scheme optimization for diffuse optical tomography based on data and image space rankings," Journal of Biomedical Optics 21(10), 106004 (24 October 2016). https://doi.org/10.1117/1.JBO.21.10.106004
Published: 24 October 2016
JOURNAL ARTICLE
12 PAGES


SHARE
Advertisement
Advertisement
Back to Top