2.1.

Formulation of the Forward Problem

In this paper, we assume frequency-resolved measurements; the light sources are amplitude modulated by a sinusoid of a specified frequency $\nu$ . Although we focus on this setup, the proposed algorithm is general enough to accommodate time-resolved and continuous-wave measurements (with an appropriate change in forward model). The widely used forward model for the photon fluence in a highly scattering medium is the Helmholtz frequency-domain diffusion equation^{4, 15, 31}

Functional activations give rise to small perturbations in ${\mu }_a\left(\mathbf{x}\right)$ due to changes in local oxy- and deoxyhemoglobin concentrations. Denoting the perturbations in the absorption coefficient and the corresponding perturbation in the fluence as $\mathrm{\Delta }{\mu }_a\left(\mathbf{x}\right)$ and $\mathrm{\Delta }\varphi \left(\mathbf{x}\right)$ , respectively, and by assuming $\parallel \mathrm{\Delta }\varphi \parallel \ll \parallel \varphi \parallel$ , we obtain

Eq. 3

2.2.

Noise Process

It is experimentally observed that the noise at the $j$ th detector due to the $i$ th source [defined by $\eta (i,j)$ ] follows a Gaussian process with zero mean and standard deviation proportional to the amplitude of the corresponding baseline signal ${\varphi }_{i,j}$ . (We assume source-detector distances of about 25 to $30\mathrm{mm}$ , including the biological noise contribution.) Collecting the measurements corresponding to the different source-detector pairs into a single vector, and denoting $f\left(\mathbf{x}\right)=\mathrm{\Delta }{\mu }_a\left(\mathbf{x}\right)$ , we obtain

2.3.

Simulation of the Forward Problem

We use Monte Carlo simulations of the optical transport problem, discussed in Refs. 16, 32, to generate the frequency-domain light bundles. The modulation frequency $\nu$ is chosen as $200\mathrm{MHz}$ . To generate realistic simulations, we use the segmentation of a magnetization-prepared rapid acquisition gradient echo (MP-RAGE) MRI brain scan along with scattering and absorbtion coefficients from Ref. 32 to obtain the light bundles. We insert a small bloblike perturbation (see Fig. 3) on the cortical region (10% of the nonstimulated absorption coefficient) and obtain the corresponding perturbed fluence measurements. We refer to the bloblike perturbation in the absorption coefficient as $f_{orig}$ and the corresponding optical perturbation as ${\mathbf{y}}_{orig}$ .

Fig. 3

Reconstruction in the absence of noise. The range of the original signal is scaled to be in the range 0 to 100; the same scale factor is applied to the reconstructions. The voxel dimensions in the above figure are $4{\mathrm{mm}}^3$ . The reconstructions are constrained to the brain (gray matter and white matter) and are displayed over the mask of the brain (gray matter and white matter). The nonbrain regions are marked as “nb” in the color bar. The top row consists of three slices of the original and level-set reconstructions each. Note that the level-set reconstructions are in close agreement with the original object. The corresponding slices of the reconstructions with SIRT and CG are displayed on the bottom row. Note that these reconstructions are distributed over a larger region and are oscillatory. The white arrows correspond to regions with negative amplitude. This effect is predominantly visible in the CG reconstructions. We used 64 iterations of CG and 400 iterations of SIRT. For the new algorithm, we used $\mu =2e-3$ , $\zeta =0.2$ , $\lambda =0.6$ .

Based on previous experimental results, we set the standard deviation of the noise process to be 2% of the baseline signal. Because the perturbation in the optical signal due to functional activation is often much smaller in amplitude as compared to the baseline, this corresponds to a very noisy scenario [the signal to noise ratio (SNR) of the perturbation signal is approximately $-15\mathrm{dB}$ ]. At such a high noise level, it is difficult to obtain reasonable reconstructions. Hence, it is common practice to average $N$ successive measurements to improve the SNR. With $N$ averages, the standard deviation of the noise process will decrease by a factor of $\sqrt{N}$ . In our experimental setup,^{16, 33} we time-domain multiplex the signals from individual channels and hence the time resolution is $20\mathrm{ms}$ per channel: we typically use around 50 to 100 averages; the resulting temporal resolution is still sufficient for hemodynamic studies.

We assumed the source-detector configuration in Ref. 16 for the simulations (see Fig. 2 ). We denote the noise corrupted version of ${\mathbf{y}}_{\mathrm{orig}}$ as ${\mathbf{y}}_{\mathrm{measured}}$ . We then use the reconstruction algorithms to obtain the estimate of the signal (denoted by $\widehat{f}$ ) from ${\mathbf{y}}_{\mathrm{measured}}$ .

Fig. 2

Schematic of the optical probe: 16 sources arranged in a circular pattern and 4 detectors in the center used to obtain 64 source-detector pairs.

2.4.

Reconstruction

By collecting the voxel values of $f\left(\mathbf{x}\right)$ into a vector $\mathbf{f}$ , the forward model (9) can be expressed as a matrix equation $\mathbf{y}=\mathsf{A}\mathbf{f}+\mathit{\eta }$ . Each voxel in the reconstructed image would correspond to a column of $\mathsf{A}$ . Similarly, a measurement corresponds to a row in $\mathsf{A}$ . In the context of sparse sampling, it is shown that if the submatrix of $\mathsf{A}$ , obtained by selecting any of its $M$ columns is invertible, an activation pattern with $M$ active voxels can be uniquely reconstructed with probability 1.³⁴ Unfortunately, the value of $M$ , for which the condition number of all $M$ column submatrices of $\mathsf{A}$ is reasonably low (well-conditioned matrix is required to ensure robust reconstructions in the presence of noise), is fairly small for the NIRSI forward model $(M\le 16)$ . Hence, we have to enforce additional constraints that are relevant to the functional imaging context (e.g., the activation pattern being smooth in the region of support, the boundary of the activations being smooth) to further reduce the degrees of freedom. We would now be searching for an activation pattern in the class of reconstructions that satisfy these constraints; we hypothesize that these constraints ensure a reasonably well-conditioned inversion problem.

The standard $l_1$ minimization approach to sparse problems^{23, 24} does not offer a convenient way to account for the additional properties of brain activations. Instead, we propose an alternate algorithm using variational principles. Assuming the activations to be support limited to subregions of the brain specified by $\Omega \subseteq D$ (see Fig. 1), we rewrite the forward model as

Note that $f$ and its support $\Omega$ in 9 are unknowns. We pose the estimation of these quantities as a numerical optimization problem; we define the cost function as

Eq. 10

The second and third terms in Eq. 10 are similar to the standard Tikhonov regularization terms used in inverse problems to make the inverse problem well-posed.¹⁴ The important difference with the standard setting is that these integrals are evaluated over $\Omega$ as compared to ${\mathbb{R}}^3$ (the entire volume) in the standard scheme. This restriction ensures that we would not be smoothing across the boundaries between the activated and the nonactivated regions, thus avoiding the blurring of the edges of the activation regions. The second term penalizes the amplitude of the reconstructions, thus eliminating reconstructions with high amplitudes. The third term imposes a penalty on the roughness of $f$ on its support. If $\lambda \to \infty$ , $\mathbf{\nabla }f$ would only be supported on $\mathrm{d}\Omega$ , the boundary of $\Omega$ ; we would be reconstructing a piecewise-constant function $f$ . The last term in Eq. 10 constrains the volume of the activations; it ensures that the estimation of the support $\Omega$ is well-posed. The parameters $\lambda$ , $\mu$ , and $\zeta$ control the strength of the additional penalties relative to the data consistency term. Note that in the absence of the second and third terms, the cost function is simply an $l_0$ minimization problem; one tries to find an $f$ with the smallest possible support.

This approach of simultaneously estimating $\Omega$ and the function values is an extension of Ref. 25 to three dimensions. This approach is also conceptually similar to Refs. 26 and 35, 36, 37, 38. A similar shape-based approach was used in the context of diffuse optical imaging in Ref. 39. They modeled the shape of the activation as an ellipsoid while the function values inside and outside the activation are approximated by harmonics. Our approach is more general because it can account for perturbations of arbitrary shape and topology. Moreover, we model the function values as arbitrary smooth functions.

The main difference between the new method and the one in Ref. 25 is the inclusion of additional Tikhonov regularization terms. In addition, the last term in Eq. 10 is also different from Ref. 25; we minimize the volume of $\Omega$ rather than the surface area of the boundary $\mathrm{d}\Omega$ . The main reason to change this term is to improve the numerical stability of the algorithm. The functional activations typically occupy a very small region of space ($\approx 10$ to 15 voxels). Hence, the approximation of the curvature term using simple finite difference operators is not good enough for this application. Furthermore, the new term minimizes the volume of the support that is desirable to ensure sparse reconstructions.

2.5.

Algorithm

In Sec. 2.4, we have seen that the reconstruction of the activation pattern involves the minimization of Eq. 10

3.1.

Validation Using Simulated Data

The Monte Carlo setup to simulate the image formation is discussed in Sec. 2.3. The baseline measurements are computed using Monte Carlo simulations on the head model, derived from segmentations of a MPRAGE head scan. Similarly, the perturbed measurements are then obtained by running the scheme on the head model with the specified perturbation in the absorption coefficients. The perturbations are shown in Figs. 3a, 3b, 3c and Figs. 4a, 4b, 4c , respectively. A noise process of standard deviation proportional to the magnitude of the baseline signal is then added to the perturbed signal. The difference between the perturbed measurements and the baseline signal would correspond to the perturbed optical signal. We then perform the reconstructions using the three algorithms and the results are compared to the original perturbations.

Fig. 4

Reconstruction in the presence of noise. The range of the original signal is scaled to be in the range 0 to 100; the same scale factor is applied to the reconstructions as well. The voxel dimensions in the above figure are $4{\mathrm{mm}}^3$ . We considered a noise process of standard deviation that is 2% of the baseline signal. We considered 50 averages of this signal and reconstructed the signal. Note that the SIRT and TCG give very smooth reconstructions. Also note that their amplitudes (from the color bar) are not in good agreement with the original. On the other hand, the level-set algorithm gave reasonable reconstructions indicating its improved resolving capability. Here we used $\mu =2e-3$ , $\zeta =0.3$ , $\lambda =0.8$ . Note that the optimal parameter set depends on the noise level and hence is different from the noise-free case considered in Fig. 3.

The optimal parameters (number of iterations in TCG and SIRT and the variables $\lambda$ , $\mu$ , and $\zeta$ in the new algorithm) are determined by comparing the reconstructions to the original, for a typical perturbation. This perturbation is different from the ones considered in the subsequent experiments, which are reported in this paper. This optimal parameter set, so derived, was then used for all experiments at this noise level. Since the parameters depend on the noise level, the optimal parameter set has to be determined for each noise level. Because this set is a function of the noise levels, the experiment had to be repeated for different noise levels (or number of averages).

The algorithm is implemented in MATLAB, with a C implementation of the support update. The decay of the criterion in Eq. 10 as a function of the number of CG iterations in the cases of single-blob and two-blob perturbations (Figs. 3 and 4) are shown in Fig. 5 . It is seen that the criterion decreases rapidly after every support update. The level-set algorithm on the average took $\sim 1\min$ to converge in comparison to $\sim 10\mathrm{s}$ for the TCG on a Pentium 4, $3.2\mathrm{GHz}$ computer.

Fig. 5

Decay of the criterion (10) as a function of the number of CG iterations. The dotted lines mark the locations of the level-set update, where the support of the activations are changed. Each support update will add or delete new voxels to the current estimate of the support. The newly added voxels are initialized by zeros to start the CG iteration. This results in a slight increase in the value of the criterion [mainly due to a high value of the third term in Eq. 10] after some of the support updates (at the steps where new voxels are added). Note that in these steps, the criterion decays to a value smaller than the one before the support update in a few iterations. Both the CG and the level-set iterations are terminated when the relative decrease in the criterion falls below a specified threshold.

3.1.1.

Qualitative comparisons

In this section, we perform qualitative comparisons of the three algorithms: SIRT, TCG, and level set. First, we compare them in the absence of noise. A few slices of the original and the reconstructed activations are overlaid on the mask of the brain (gray-matter and white-matter regions) in Fig. 3. Note that the SIRT and CG reconstructions are distributed over a much larger region of the mask. There are also some regions on these estimates where the reconstructions are oscillatory (the white arrows in the figure show regions of negative amplitudes). The noisy-looking reconstructions of the CG could be due to the oscillatory nature of the singular vectors of the system matrix $\mathsf{A}$ . These oscillatory reconstructions were also reported in Ref. 16. On the other hand, the new approach gives reconstructions that are similar to the original signal.

To understand the flexibility and limitations of the algorithm, we study the performance of the algorithm for two different activations in Fig. 6 . In the first case, we consider two blobs with different amplitudes. It is seen that the level-set algorithm provides a reasonably good recovery of the amplitudes and the shape of the reconstructions. We then considered a single activation that is deeper in the brain. In this case, the reconstructions are biased to the cortical surface as expected. This is due to the exponential drop in sensitivity of the optical technique with depth. With standard reconstruction schemes, this bias is predominant and often leads to artifacts as seen in Fig. 9.

• 3. Support centroid error (SCE): We choose the distance between the centroids of the original activation and that of the detected activation as the error measure. This approach uses the classification from the earlier metric.

Fig. 6

Reconstruction for various test activation patterns. Two adjacent slices of each data set are displayed. We consider an activation pattern with two blobs that have different amplitudes on the left. It is seen that the level-set reconstruction scheme give a reasonably good reconstruction of the activations. On the right, we test the reconstructions for a deeper activation. It is seen that the level-set reconstructions are biased to the surface of the cortex. Also note that the amplitude of the reconstructions are much lower than the original value. This is because light bundles decay exponentially with depth. Here, we used the same noise level and same parameters as in Fig. 4.

In Fig. 4, we study the case where there are two activated regions in the presence of noise. Note that the SIRT and TCG algorithms give a uniformly smooth reconstruction; these algorithms are not capable of resolving the two activations. On the other hand, the new algorithm results in two distinct activations, although the shapes are distorted. This clearly indicates that the new approach can resolve adjacent structures better.

3.1.2.

Quantitative comparisons

In this section, we compare the different algorithms by comparing the estimates with the original signal $f_{orig}$ . Here, we considered the case with only one activation. We use different metrics for comparison:

• 1. Mean squared error (MSE): The normalized mean squared error in estimation is given by

  Eq. 23

  $$\mathrm{MSE}=\frac{{\int }_{{\mathbb{R}}^3}\mid f_{\mathrm{orig}}\left(\mathbf{x}\right)-f_{\mathrm{est}}\left(\mathbf{x}\right){\mid }^2\mathrm{d}\mathbf{x}}{{\int }_{{\mathbb{R}}^3}\mid f_{\mathrm{orig}}\left(\mathbf{x}\right){\mid }^2\mathrm{d}\mathbf{x}}.$$
  While this error term is widely used and is simple, its value can be misleading in our setting. The level-set algorithm produces piecewise smooth reconstructions as compared to smooth reconstructions with SIRT and CG. Because the optical light bundles are smooth and few in number, it is not possible to precisely localize the functional activations. The MSE would indicate larger errors for the level-set method, even when the piecewise smooth reconstructions are only slightly off. On the other hand, this metric gives small values while comparing smooth functions, even when the locations are wrong.

• 2. Support error (SE): We use the standard MATLAB (Mathworks, Natick, New Jersey) implementation of the $k$ -means algorithm to classify the reconstructions into two regions (activated and background). Examples of the segmentations of CG and SIRT reconstructions are shown in Fig. 7 . We then compare the support of the estimated activations to the original support using the metric

  Eq. 24

  $$\mathrm{SE}=\frac{{\int }_{{\mathbb{R}}^3}\mid {\chi }_{{\Omega }_{\mathrm{orig}}}\left(\mathbf{x}\right)-{\chi }_{{\Omega }_{\mathrm{est}}}\left(\mathbf{x}\right){\mid }^2\mathrm{d}\mathbf{x}}{{\int }_{{\mathbb{R}}^3}\mid {\chi }_{{\Omega }_{\mathrm{orig}}}\left(\mathbf{x}\right){\mid }^2\mathrm{d}\mathbf{x}},$$
  where the characteristic function of a region $\Omega$ is defined as

  Eq. 25

  $${\chi }_{\Omega }\left(\mathbf{x}\right)=\{\begin{array}{cc}\hfill 1,\hfill & \hfill \mathrm{if}\mathbf{x}∊\Omega ,\hfill \\ \hfill 0,\hfill & \hfill \mathrm{otherwise}.\hfill \end{array}$$
  Note that the above metric counts the misclassified voxels in the reconstruction and hence ignores any amplitude changes. This metric can be argued to make better sense for the comparison of reconstructions of functional activations.

Fig. 7

Examples of segmentations of a slice of the reconstructions using $k$ -means algorithm. The error between the original support and the segmentation is used as a criterion for comparison. See Sec. 3.1.2 for details.

Fig. 9

Reconstruction of perturbation in oxyhemoglobin $\left(\mathrm{\Delta }{\mathrm{HbO}}_2\right)$ and deoxyhemoglobin $\left(\mathrm{\Delta }\mathrm{Hb}\right)$ concentrations (in nanomoles/liter) in comparison with the fMRI maps. Because many averages were used to suppress the noise, we used the same parameter set as in Fig. 3. The reconstructions are constrained to the brain (gray matter and white matter). Three slices overlaid on the mask of the brain are displayed. As expected, the NIRSI reconstructions are biased to the cortical surface due to the exponential decay of sensitivity with depth. However, in comparision to CG and SIRT, the level-set reconstructions more closely resemble the fMRI activation maps.

Using these error measures, we quantitatively compare the performance of the algorithms (level-set, CG, and SIRT) as a function of the number of averages for the single activation case (Fig. 3). Note that the SNR is proportional to the square root of the number of averages. The results are displayed in Fig. 8 . Note that the new algorithm fares better under both SE and SCE, in clear agreement to the visual comparisons. However, the MSE indicates a higher error in comparison to the others. As we have explained before, this is due to the inherent weakness of MSE to compare piecewise constant functions. Note that the errors of any of the algorithms does not decay to zero even when there is no noise, indicating the insufficiency of the data (ill-posed nature of the problem).

Fig. 8

Quantitative comparison of the performance of the algorithms to noise. We plot the reconstruction error, evaluated using different metrics, as a function of the SNR of the measurements. Note that the level-set algorithm fares significantly better than the other under the SE and SCE criteria, in clear agreement with the visual comparisons. The MSE criterion indicates worse results. As discussed below, this is mainly due to the limitation of the MSE in comparing two piecewise constant functions.

3.2.

Experimental Results

We now study the performance of the algorithm to experimental data. The optical probe was designed with 16 pairs of 400-mm-diameter core plastic-clad multimode silica source fibers and 4 detector fiber bundles. The topology of the probe is shown in Fig. 1. A MRI visible marker is attached adjacent to each of the 16 optical source fibers and 4 detector fiber bundles so that accurate source and detector positions could be estimated from MR images. Experimental fMRI and NIRSI data were obtained from a normal subject, who was scanned using a Siemens Allegra 3T system. The study was approved by the University of Illinois at Urbana-Champaign Institutional Review Board. Written informed consent was obtained from the volunteer before the study began.

The optical signals were recorded using a near-infrared spectrometer (Imagent, ISS, Champaign, Illinois). Data acquisition was synchronized with the fMRI measurements using the transistor-transistor logic (TTL)-trigger signal from the MR scanner, which also triggered the beginning of the visual stimulation paradigm. The optical sources were laser diodes (690 and $830\mathrm{nm}$ ) that are amplitude modulated at $150\mathrm{MHz}$ and time-multiplexed. Light reaching the detectors was amplified by photomultiplier tubes and converted into ac, dc, and phase signals for each of the source-detector combinations, or channels, at each wavelength.

A pair of nonmagnetic goggles (Resonance Technology, Inc., North Ridge, California) with liquid crystal display screens were placed in front of the subjects’ eyes inside the birdcage head coil. The visual stimulation paradigm consisted of five blocks each with 28.8-s fixation followed by $19.2\mathrm{s}$ of a black-and-white checkerboard pattern flashing on for $50\mathrm{ms}$ and off for $1.95\mathrm{s}$ .

Automatic estimation of source and detector positions from the MR images, incorporating deformation of the optical probe, is performed using custom-developed image processing algorithms written in MATLAB. The optical data is averaged over 80 activation-relaxation trials to produce the estimate of the signal difference between the active and baseline conditions.¹⁶ The level-set reconstructions of the perturbation in absorption coefficients are performed at the two wavelengths. The perturbation in the oxy- and the deoxyhemoglobin concentrations are estimated using the Lambert-Beer law.

The estimated perturbation maps, overlaid on the brain surface, are shown in Fig. 9 . One would expect a systematic increase and decrease in the oxy- and deoxyhemoglobin concentrations, respectively, at the locations corresponding to the fMRI activations. It is seen from the reconstructions that this is roughly the trend in the NIRSI reconstructions. However, unlike the fMRI image, both CG and SIRT reconstructions indicate activations in a shallow region of the brain without separation between left and right sides of the visual cortex. On the contrary, the level-set images show two separate activated regions on the left and right sides of the brain, which are deeper and similar to the fMRI images. All optical images show stronger activation in the left hemisphere because in this location the activation region (as revealed by fMRI) is slightly closer to the surface of the head.