## Abstract

We evaluate the ultimate transverse spatial resolution that can be expected in Diffuse Optical Tomography, in the configuration of projection imaging. We show how such a performance can be approached using time-resolved measurements and reasonable assumptions, in the context of a linearized diffusion model.

© 2009 Optical Society of America

## 1. Introduction

Diffuse optical tomography (DOT) in the near infrared is the subject of numerous and fruitful researches for more than two decades, with a lot of promising applications in biomedical optics for oxygenation monitoring [1–9], neurosciences [10–15] or for mammography [16–21]. As the diffusion of light in tissue blurs very efficiently the spatial information, DOT implies the resolution of a complex ill-defined inverse problem. This usually leads to a poor spatial resolution. A way to improve this resolution is to couple DOT techniques with other imaging techniques like MRI, in order to inject an *a priori* information for the inversion procedure [22–24]. But this solution tends to nullify the low cost advantage of DOT techniques. In this paper we are going to investigate the ultimate resolution limits that should be reasonably expected from pure DOT techniques. We will focus on the transverse resolution of projection imaging, commonly used in mammography [18–21]: the breast is modelled by a slab of thickness d, and the source and the detector are face-to-face at both sides of the slab (see Fig. 1); the source-detector system translates over the whole plane, in order to get a projection image.

Within the frame of the diffusion approximation, the time-resolved (TR) Green function of this problem satisfies the diffusion Eq. (25):

where c is the speed of light in the medium, *D* the diffusion constant, *µ _{a}* the absorption coefficient, and

*is the position of a point-like source, located for instance at*

**r**_{s}*z*=

*l** to model an incident laser beam [25] (

*l**~

*3D*being the transport mean free path; let us recall that

*l**=

*µ*’

_{s}

^{−1}, where

*µ*’

_{s}is the reduced scattering coefficient). Boundary conditions can be accounted of in many way with quite the same results at large scale [26], taking for instance

*G*(

**)=0 on extrapolated planes**

*r**z*=−

*0.7l** and

*z*=

*d*+

*0.7l**. If

**=**

*r**lies on the exit surface at*

**r**_{d}*z*=

*d*,

*G*(

*) is directly linked to the transmittance and is proportional to the measured signal.*

**r**_{d}In most cases, inverse problems [27] are based on a linearization of the diffusion Eq. (1) with respect to the optical properties *D* and *µ _{a}* of the medium. For instance, when writing

*µ*(

_{a}**)=**

*r**µ*+

_{a0}*δµ*(

_{a}**) and if**

*r**G*is the Green function for the background

_{0}*µ*(

_{a}**)=**

*r**µ*, one can use the Rytov approximation:

_{a0}The aim of the inverse problem is to invert this linear set of equations and to reconstruct *δµ _{a}* from the contrast function

*C*(

*). The contrast function exhibits the information specifically linked to the inhomogeneity*

**r**_{d},t*δµ*. There are several tricks to experimentally get this function, using measurements at a reference position [18], or using the administration of a dye that induces high variations of

_{a}*µ*and low variations of

_{a}*D*[28]; another important example concerns brain functional imaging, where changes with respect to a reference state can also allow the derivation of a contrast function (see for instance references [12,13]); there is also the opportunity offered by diffusing wave spectroscopy to specifically vary the effective absorption coefficient µa with the control of the correlation time [29]. We will therefore consider in the following that this function is experimentally available, and that the only function we need to compute is the integral kernel in the right member of (2). To this end, we will assume that the background has well-known constant optical properties

*D*and

_{0}*µ*: the problem of the indeterminacy of these coefficients is left to the end of this paper. With such an assumption

_{a0}*G*can be easily obtained [25] and the integral kernel presents a translational symmetry with respect to transverse coordinates

_{0}*=(*

**ρ***x,y*), and (2) becomes:

With the technique of projection imaging, the source and the detector are face-to-face so that * ρ_{d}*=

*. The incident laser beam and the detector however have a spatial extension, that can be accounted for with convolution products over*

**ρ**s*and*

**ρ**_{d}*, leading to a contrast function C(ρ,t) that only depends on the central position of the source-detector system. These results can be extended to the continuous wave (CW) case through integration over time.*

**ρ**_{s}## 2. The ultimate spatial resolution

#### 2.1 Theoretical background

As said above, the inverse problem usually consists to invert an ill-posed linear set of equations *f*=*Lg*, and more precisely to optimize the cost function

where a Thikonov regularization term of strength *λ* is added in order to prevent large oscillations in the reconstructed image. The choice of *λ* is not obvious and is quite arbitrary. Several works are devoted to it, and we are going to discuss it in the present paper. With a linear system the minimum of this cost function is:

In this section we are going to deal with a very simple problem in which the inclusion *δµ _{a}* is supposed to be located at a well-defined position

*z*. In that case the contrast function (3) is a simple convolution product

*C*(

*)=*

**ρ***κ*⊗

_{z,t}*δ*(

_{µa}*), and*

**ρ***C*is a blurred version of

*δµ*. Presented in that way, the problem seems to be over-simplified, but it is an interesting way in order to probe the ultimate limits for the spatial resolution of inverse problems: the next section will be devoted to a more general case. A convolution product can be easily inverted in the Fourier space, and (5) directly leads to

_{a}*g*̃(

*k*)=

*W*̃

_{z,t}(

*k*)

*f*̃(

*k*), with:

W looks like the Wiener filter, extensively used in image processing [30]. In the Wiener filter however the parameter *λ* depends on the spatial frequency and also on an *a priori* knowledge of the power spectrum of the initial image. Here W gives the solution of an inverse problem, and *λ* is a Thikonov regularization parameter that needs to be fixed [31]. The deconvolution is characterized by the transfer function *h*̃_{z,t}=*W*̃_{z,t}
*κ*̃_{z,t}, which tends to *1* as *λ* tends to *0*. In this paper the spatial resolution will be defined as the Full Width at Half Maximum (FWHM) of this transfer function *h( ρ)*. Of course, a numerical simulation with

*λ*=

*0*will lead to unrealistically good results, but such results won’t live through a little bit of noise. In fact, the Thikonov parameter

*λ*allows

*W*to filter the noise at high spatial frequencies: if the noise is high, a high value is required for

*λ*. Let us quantify this. We consider an image of

*N*x

*N*pixels, with a distance

*Δx*between pixels and a uniform white noise

*b*. As the filter

_{ij}*W*is linear, we can focus on its action on this noise. If we introduce the discrete Fourier transform (DFT)

*b*̃

_{ij}of

*b*, we have for the variance of the noise after filtering:

_{ij}In the continuum limit, with Δ*k*=*2π*/(*N*Δ*x*), * k_{ij}*=(

*i*Δ

*k,j*Δ

*k*) and ∑→Δ

*k*

^{−2}∫

*dk*

^{2}we have

With diffusion processes, the kernel *κ* is nearly Gaussian [32], and it is remarkable that for a Gaussian kernel

the integral (8) can be explicitly calculated:

Note that *α* is directly linked to the FWHM of the kernel, *ie* to the spatial resolution *R* before deconvolution, with *α*=*R*
^{2}/(16ln 2). If we introduce the signal to noise ratio (SNR) of the data we can directly write from Eq. (7):

Let us note that this expression of *γ* does not depend on the definition of the ‘Signal’: here we only assume that the signal is not lowered by the deconvolution procedure, so that *SNR _{f}* is a lower bound for the SNR after deconvolution. From Eqs. (10) and (11) we have:

and *λ* depends on the SNR as expected. This formula enlighten the arbitrariness in the choice of *λ*, which appears through the choice of *SNR _{f}*. The user has to choose the SNR he tolerates. Let us also note that (12) depends on

*Δx*: the smaller

*Δx*is, the more numerous the implied measurements are. This is equivalent to a better SNR. For

*λ*=

*SNR*

^{−2}, the SNR after deconvolution will be

*SNR*≡

_{f}*SNR*

_{0}≈1,5

*R*/Δ

*x*; for instance, for

*Δx*=

*1mm*and

*R~1cm*one should await a

*SNR*of about

_{f}*15*. If one wants

*SNR*=

_{f}*g*

*SNR*has to be fixed to:

_{0}, λwhere *g* is a SNR enhancement factor. Let us note that *SNR _{f}* is fixed by the user, and that we can have

*SNR*<

*SNR*: W is a filter and can reduce the noise. Let us now see how λ is linked to the resolution: the transfer function

_{f}*h*̃ is approximately equal to

*1*if

*κ*̃≫√

*λ*; for the Gaussian kernel (9), this defines a disc in the Fourier plane of radius ${k}_{\mathrm{max}}=\sqrt{-\mathrm{ln}\lambda /(2\alpha )}$. If

*R*is the spatial resolution before deconvolution, we therefore roughly get for the spatial resolution

*R*after deconvolution:

_{d}where the numerical constant *1.41* is the FWHM of the function *J*
_{1}(*πx*)/*x*. We have numerically checked this formula, as shown in Fig. 2: the ratio *R _{d}/R* was numerically evaluated with a convolution/deconvolution procedure with a gaussian kernel applied to a single point (blue curve), and the result is compared to formula (14) (red curve). This formula appears to be quite good for low values of

*λ*, and is still valid up to

*6*% at λ=10

^{−2}.

Let us note that for low values of the SNR (*ie* λ ≥ 0.1), we can have *R _{d}*>

*R*. This does not mean that the filter

*W*is useless in that case: the noise has to be filtered, and the filter

*W*does this job with a degradation of the resolution. In that case we are however outside the field of inversion procedures and resolution enhancement, and the fact that (14) overestimates R

_{d}is not a problem.

#### 2.2 Numerical evaluation of the resolution limit

One purpose of this paper is to compare the resolution limit of the TR versus the CW cases. It is indeed well-known [33] that the resolution of the kernel κ is improved for short values of the transit time *t*. The SNR however decreases: this limits the resolution gain expected from the deconvolution algorithm and introduces an optimal measurement time. In the CW case, the initial resolution in worse, but the SNR is better. The formulas (13–14) allow a quantitative evaluation for these questions. For this we only need an evaluation of the SNR. For a detector directly placed at the contact of the output surface, the measured background signal, expressed in photo-electrons, is roughly ½*η*
*G _{0}Δτ*, where

*Δτ*is the temporal resolution of the detector and:

where *E* is the total incident energy during the exposition, S is the surface of the detector, and *a* its quantum efficiency. In the CW case, the signal is simply ½*η G*. The fundamental limit for the SNR is the quantum noise, ie the square root of the signal in photo-electrons. Quantum noise corresponds to a fundamental limit because it is the only source of random noise that cannot be cancelled, and that therefore has to be considered for an investigation of the ultimate resolution. Furthermore, reaching this limit is not a great difficulty with nowadays technology. The contrast function needs two measurements, one with and another without the inclusion, so that the corresponding quantum noise is approximately

Let us recall that the background Green function G_{0} is not supposed to vary a lot with transverse coordinates, and this quantum noise therefore almost has a constant variance, what *a posteriori* justifies the white noise model used to derive Eq. (10). The contrast function can be roughly evaluated as *C*≈*κ _{z,t}δµ_{a}V*, where

*δµ*is the typical absorption contrast of the inclusion and

_{a}*V*its volume. This leads to:

with:

A lot of parameters linked to the experimental setup are integrated in the multiplicative constant *F* defined in Eq. (18). This constant quantifies the quality of the setup, and for instance increases with the incident light energy, the detector quantum efficiency, the detector surface or the object contrast. It can be seen as an adjustable factor that accounts for all these parameters. In the CW case, the *Δτ* term in Eq. (17) has to be dropped.

In order to increase the SNR, one should use a large detector, and a high source power too; but for safety reasons with biomedical applications, this last point implies a spreading of the incident energy with a large source. As said in the introduction, we account for such spatial extensions of the source and the detector using additional convolution products, that should increase the FWHM of the kernel *κ*, and therefore deteriorate the spatial resolution. Our numerical simulations however show that, when compared to a null diameter, a diameter of 5mm for the source and the detector induces very low changes to the kernel *κ* for an object located in the middle of the slab. We therefore introduce this diameter of 5mm through suitable convolutions. In the same way, the time integration window *Δτ* of the detector should be the highest as possible in order to maximize the SNR in Eq. (17). *Δτ* should also remain small enough compared to typical transit times in order to perform a correct measurement of time-resolved data. As we are going to deal with typical transit times of a few nanoseconds, the value *Δτ*=*100ps* appeared to us as a good compromise.

Equations (13),14,17,18) now allow us to estimate the spatial resolution *R _{d}* that can be achieved after deconvolution. We performed numerical simulations with the same parameters as in ref [34]. to mimic the breast: a thickness

*d*=

*6cm*for the slab, a refractive index

*n*=

*1.4*,

*µ*’

_{s0}=

*0.75mm*

^{−1}(with

*D*=

_{0}*1/[3(µ’*) and

_{s0}+µ_{a0})]*µ*=

_{a0}*0.0025mm*and

^{−1}, 0.005mm^{−1}, 0.0075mm^{−1}*0.01mm*. The Fig. 3 presents the result obtained in the TR case, with

^{−1}*µ*=

_{a0}*0.005mm*, for the resolution

^{−1}*R*as a function of the transit time. This function

_{d}*R*clearly exhibits a minimum at a transit time

_{d}(t)*t*(which is 782ps in this case) that is short enough to allow a small width of the kernel κ, but not too short for the SNR to allow a good deconvolution.

_{opt}With the TR measurements, we observed that *R _{d}*(

*t*) always presents such a minimum. The Fig. 4(a) presents this optimal value

*R*(in mm) as a function of

_{d,opt}*F*. It clearly decreases for high values of

*F, ie*for good experimental sensibilities, but this graph is in fact presented with a logarithmic scale for the

*F*axis: this optimal resolution in fact does not evolve so much for a given order of magnitude of

*F*. In the same way, it only slightly depends on the absorption coefficient.

At this stage, let us try to find a reasonable range for *F*: we can for instance assume an incident laser energy *E*~*10mJ*, a detector surface of *S*~*20mm ^{2}* with a quantum efficiency

*a~8*% (the detector can for instance be a time gated intensified CCD camera). This leads to

*η~6.10*. For a SNR enhancement factor

^{16}*g*=

*3*and

*δµaV~2.5mm*(

^{2}*δµ*[28] and

_{a}~0.005mm-1*V~500mm*), this gives

^{3}*F~10*. This is of course a rough estimation, and all those parameters can greatly change depending on the experiment: for instance an incident energy of only

^{8}*1mJ*will lower the factor

*F*from

*10*to 3.

^{8}*10*, and the corresponding optimal spatial resolution will be increased of about 15%. This model therefore account for a large class of experiments. In the following we will deal with the case

^{7}*F*=

*10*, but our conclusions can be extended to a large range of values for this parameter.

^{8}For an absorption coefficient of *µ _{a0}*=

*0.005mm*and

^{−1}*F*=

*10*, we therefore assert that the optimal resolution in the time-resolved case is of about

^{8}*9.3mm*at an optimal transit time

*t*=

_{opt}*780ps*. This corresponds to a FWHM of the kernel

*κ*of

*R*=

*14.5mm*, corrected to

*9.3mm*after the Wiener filter. The Fig. 4b presents the gain in the spatial resolution using TR instead of CW. This gain decreases with

*F*, since the SNR is lower in the TR case. Furthermore, it decreases when

*µ*increases: in fact, an increase of

_{a0}*µ*penalizes long transit times, so that the absorption acts like an early photon time gate. Figure 4b shows the comparative advantage of TR technologies, as the spatial resolution is of only ~11mm in the CW case for

_{a}*F*=

*10*.

^{8}As such values are often appearing in the literature, a comment is needed at this stage: a numerical simulation is presented on Fig. 5; we consider an image (Fig. 5a) with: a single point, 2 discs of 10mm diameter distant from 5mm, and 4 sets of bars; in each set, the bars of thickness *e* are separated by the same distance *e*, where e takes the values *8mm, 7mm, 6mm* and *5mm*, depending on the set. This image is then convoluted by the TR kernel *κ* at *t _{opt}*=

*780ps*(Fig. 5e), where the objects are supposed to be located at the center of the slab (

*z*=

*30mm*). A Gaussian noise is added according to (17). The Wiener filter W is then used for image restoration (Fig. 5f), and the improvement is quite clear. The definition we used for the spatial resolution is the FWHM of the transfer function h, that is the FWHM of the image of the single point in Fig. 5f. But the two discs can be separated in Fig. 5f, although they are separated from only 5mm: In fact, the definition we use for the spatial resolution refers to the separation distance between the centers of objects rather than between their outer surface [35]. For instance, the bars can be clearly seen in Fig. 5f, except for a separation distance of

*10mm*(

*e*=

*5mm*) for which they can be still distinguished: this corresponds to the resolution limit. So the main result here is that

*5mm*thick objects separated from

*5mm*,

*ie*from

*d/12*(where

*d*is the thickness of the slab), should be distinguishable. Note at this stage that

*5mm*does not correspond to the minimal size that can be seen: in fact there is no minimal size, and an object can be seen provided

*δµ*is not too small. The point-like object is visible in all the images of Fig. 5, but its contrast was 40 times higher than for the other objects.

_{a}VFigure 5c shows the convoluted image with added noise in the CW case, and Fig. 3d shows the restored image. The *e*=*5mm* bars in this last figure cannot be seen any more: there is still a periodic structure, but there are only 5 oscillations instead of 6; in fact, as we are at the resolution limit, the Fourier components corresponding to 6 oscillations are more filtered than spatial frequencies that are just a little bit lower, leading to this periodic structure with only 5 oscillations. One should say that this is still a good result, but the SNR in this case is of about 6000, while it was 300 in the TR case. As a consequence, the restoration implies a highly good knowledge of the kernel *κ*. For instance, we can use a Gaussian fit of this kernel with less than 2% error; the Fig. 5b corresponds to a restoration of Fig. 5c using this Gaussian fit instead of the true function in the Wiener filter *W* (Eq. (6), and only the *e*=*7mm* and *e*=*8mm* bars are visible in the resulting image. In fact, seeing bars with *e*=*8mm* in the CW case through a *6cm* thick slab is still a really good performance [34], which is seldom reached in the field. One outcome of the present paper is that such results can be improved in the TR case. Of course we were dealing in this section with objects located at a well-defined depth, in order to get an ultimate resolution. We have now to discuss how such a resolution can be obtained without such an assumption.

## 3. Restoration of a 3D projection image

In the previous section, we have seen how the Wiener filter allows us to reach a ~*9mm* spatial resolution starting from a *R*=*14.5mm* FWHM of the kernel κ in the TR case. The resolution gain due to image restoration is about *60*%, which is high, but is not so much. In fact, the deconvolution of a Gaussian of FWHM R by a Gaussian of FWHM equal to *80*% of *R* gives exactly a Gaussian with a FWHM equal to *60*% of *R*. This suggests that if we use, in the Wiener filter (6), a kernel that is *80*% narrower than the real kernel *κ*, this should not greatly affect the performances of the image restoration. In the TR case, *R*(*z,t*) roughly varies as [33]:

For a fixed transit time *t*, if *R* is the value of (19) at *z*=*d/2, 80*% of this value is reached at *z*=*d/5*. We therefore performed again the same simulation as in section *2.2*, except that we compute, in all cases, a kernel corresponding to position *z*=*d/5*=*12mm* (and to the optimal transit time *t*=*780ps* of section *2.2*) for the Wiener filter; this filter is then used to restore the projection image obtained for different ‘real’ positions z of the objects.

We plot in Fig. 6b the result obtained for *z*=*30mm*; this result can be compared to Fig. 5f, which is reproduced in 6a: the spatial resolution we obtain in 6b is not very different from what we obtained in the previous section; it increases of about 6%, and the quality of the two images is comparable (6b is only a bit noisier). But of course, the restoration can be used at any depth between *12* and *48mm*. Figure 6c and 6d thus corresponds to the cases *z*=*20mm* and *z*=*12mm*, leading to quite good results.

In fact, the kernel used for the restoration is not ought to be identical to the true kernel: the main requirement is that the transfer function *h* does not amplify high frequencies, *ie* that |*h*̃|≤1, and this is the case in the examples we considered above. If however the FWHM of the restoration function is higher than the FWHM of the real kernel, this is not the case and it can induce oscillations like in Fig. 5b. A problem can therefore appear for an object close to the surface, but simulations performed for *z*=*10mm* and *z*=*5mm* (Fig. 6e and 6f) show that this is not so much a problem. Furthermore, considering again Eq. (19), 80% of the value of *R* is reached when *D* is reduced by 35%: there is therefore a quite high tolerance concerning the value of *D*, and one can expect that small fluctuations of D won’t induce too much perturbations, provided they are not registered in the contrast function (Eq. (3): see discussion in section 1).

The results presented in this section are based on the fact that the resolution gain due to image restoration is not too high. In our simulations, this resolution gain is much higher in the CW case, thanks to a high SNR. This method we propose for the restoration of a 3D projection image will therefore be less efficient in the CW case, as shown on Fig. 7: we consider an object in the middle of the slab (*z*=*30mm*), and the restoration is performed with the correct kernel in Fig. 7a (which is the same as Fig. 5d) and with a *z*=*12mm* kernel in Fig. 7b. This last figure clearly shows a resolution loss with this method in the CW case.

In summary, the depth position of an object is not a compulsory data for image restoration in projection imaging. As a consequence, such a knowledge will be very difficult to obtain only from projection data: 3D reconstruction needs more information, like for instance measurements with different relative positions between the source and the detector [34]. A restoration of projection images implied in a tomographic reconstruction procedure can also be obtained using the method presented in this paper, so that those performances obtained for projection imaging can certainly be extended to 3D imaging. Anyway, projection imaging has many applications, especially for mammography, and it could be really interesting to experimentally reach such performance to improve the pertinence of the diagnosis, and perhaps to detect early stage tumours.

## 3. Conclusion

We have evaluated in a simple case the ultimate transverse spatial resolution that can be expected in Diffuse Optical Tomography, when we use projection imaging and time-resolved measurements. We enlighten the existence of an optimal transit time for these measurements, in order to obtain the best resolution. We show that this resolution can be almost reached without any knowledge of the depth position of the objects, and that a significant improvement can be expected compared to the performances commonly obtained.

## References and links

**1. **F. F. Jöbsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science **198(4323)**, 1264–1267 (1977).
[CrossRef]

**2. **J. R. Wilson, D. M. Mancini, K. McCully, N. Ferraro, V. Lanoce, and B. Chance, “Noninvasive detection of skeletal muscle underperfusion with near-infrared spectroscopy in patients with heart failure,” Circulation **80(6)**, 1668–1674 (1989).
[CrossRef]

**3. **T. Hamaoka, H. Iwane, T. Katsumura, T. Shimomitsu, N. Murase, S. Nishio, T. Osada, T. Sako, H. Higuchi, M. Miwa, and B. Chance, “The quantitative measures of muscle oxygenation by near infrared time-resolved spectroscopy,” Med. Sci. Sports Exerc. **28(Supplement)**, 62 (1996).

**4. **S. Takatani and J. Ling, “Optical oximetry sensors for whole blood and tissue,” IEEE Eng. Med. Biol. Mag. **13(3)**, 347–357 (1994).
[CrossRef]

**5. **R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini, “Compact tissue oximeter based on dual-wavelength multichannel time-resolved reflectance,” Appl. Opt. **38(16)**, 3670–3680 (1999).
[CrossRef]

**6. **E. M. C. Hillman, J. C. Hebden, M. Schweiger, H. Dehghani, F. E. W. Schmidt, D. T. Delpy, and S. R. Arridge, “Time resolved optical tomography of the human forearm,” Phys. Med. Biol. **46(4)**, 1117–1130 (2001).
[CrossRef]

**7. **A. Kienle and T. Glanzmann, “In vivo determination of the optical properties of muscle with time-resolved reflectance using a layered model,” Phys. Med. Biol. **44(11)**, 2689–2702 (1999).
[CrossRef]

**8. **R. A. De Blasi, S. Fantini, M. A. Franceschini, M. Ferrari, and E. Gratton, “Cerebral and muscle oxygen saturation measurement by frequency-domain near-infra-red spectrometer,” Med. Biol. Eng. Comput. **33(2)**, 228–230 (1995).
[CrossRef]

**9. **G. Yu, T. Durduran, and G. Lech, C. Zhou, B. Chance, E. R. Mohler 3rd, and A. G. Yodh, “Time-dependent blood flow and oxygenation in human skeletal muscles measured with noninvasive near-infrared diffuse optical spectroscopies,” J. Biomed. Opt. 10(2), 024027 (2005).

**10. **J. C. Hebden, A. Gibson, R. M. Yusof, N. Everdell, E. M. C. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek, and J. S. Wyatt, “Three-dimensional optical tomography of the premature infant brain,” Phys. Med. Biol. **47(23)**, 4155–4166 (2002).
[CrossRef]

**11. **A. P. Gibson, T. Austin, N. L. Everdell, M. Schweiger, S. R. Arridge, J. H. Meek, J. S. Wyatt, D. T. Delpy, and J. C. Hebden, “Three-dimensional whole-head optical tomography of passive motor evoked responses in the neonate,” Neuroimage **30(2)**, 521–528 (2006).
[CrossRef]

**12. **D. K. Joseph, T. J. Huppert, M. A. Franceschini, and D. A. Boas, “Diffuse optical tomography system to image brain activation with improved spatial resolution and validation with functional magnetic resonance imaging,” Appl. Opt. **45(31)**, 8142–8151 (2006).
[CrossRef]

**13. **M. A. Franceschini, V. Toronov, M. Filiaci, E. Gratton, and S. Fantini, “On-line optical imaging of the human brain with 160-ms temporal resolution,” Opt. Express **6(3)**, 49–57 (2000), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-6-3-49. [CrossRef]

**14. **D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head”, Op. Express **10**, 159–170 (2002), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-10-3-159

**15. **A. Liebert, H. Wabnitz, J. Steinbrink, H. Obrig, M. Möller, R. Macdonald, A. Villringer, and H. Rinneberg, “Time-resolved multidistance near-infrared spectroscopy of the adult head: intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons,” Appl. Opt. **43(15)**, 3037–3047 (2004).
[CrossRef]

**16. **A. Pifferi, J. Swartling, E. Chikoidze, A. Torricelli, P. Taroni, A. Bassi, S. Andersson-Engels, and R. Cubeddu, “Spectroscopic time-resolved diffuse reflectance and transmittance measurements of the female breast at different interfiber distances,” J. Biomed. Opt. **9(6)**, 1143–1151 (2004).
[CrossRef]

**17. **S. Srinivasan, B. W. Pogue, B. Brooksby, S. Jiang, H. Dehghani, C. Kogel, W. A. Wells, S. P. Poplack, and K. D. Paulsen, “Near-infrared characterization of breast tumors in vivo using spectrally-constrained reconstruction,” Technol. Cancer Res. Treat. **4(5)**, 513–526 (2005).

**18. **B. Wassermann, A. Kummrow, K. T. Moesta, D. Grosenick, J. Mucke, H. Wabnitz, M. Möller, R. Macdonald, P. M. Schlag, and H. Rinneberg, “In-vivo tissue optical properties derived by linear perturbation theory for edge-corrected time-domain mammograms”, Op. Express **13**, 8571–8583 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-21-8571 [CrossRef]

**19. **D. Grosenick, K. T. Moesta, M. Möller, J. Mucke, H. Wabnitz, B. Gebauer, C. Stroszczynski, B. Wassermann, P. M. Schlag, and H. Rinneberg, “Time-domain scanning optical mammography: I. Recording and assessment of mammograms of 154 patients,” Phys. Med. Biol. **50(11)**, 2429–2449 (2005).
[CrossRef]

**20. **D. Grosenick, H. Wabnitz, K. T. Moesta, J. Mucke, P. M. Schlag, and H. Rinneberg, “Time-domain scanning optical mammography: II. Optical properties and tissue parameters of 87 carcinomas,” Phys. Med. Biol. **50(11)**, 2451–2468 (2005).
[CrossRef]

**21. **P. Taroni, A. Torricelli, L. Spinelli, A. Pifferi, F Arpaia, G. Danesini, and R. Cubeddu, “Time-resolved optical mammography between 637 and 985 nm: clinical study on the detection and identification of breast lesions,” Phys. Med. Biol. **50(11)**, 2469–2488 (2005).
[CrossRef]

**22. **A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. A. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt. **42(25)**, 5181–5190 (2003).
[CrossRef]

**23. **P. K. Yalavarthy, B. W. Pogue, H. Dehghani, C. M. Carpenter, S. Jiang, and K. D. Paulsen, “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express **15(13)**, 8043–8058 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-13-8043. [CrossRef]

**24. **R. L. Barbour, H. L. Graber, J. W. Chang, S. L. S. Barbour, P. C. Koo, and R. Aronson, “MRI-guided optical tomography: Prospects and computation for a new imaging method,” IEEE Comput. Sci. Eng. **2(4)**, 63–77 (1995).
[CrossRef]

**25. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. **28(12)**, 2331–2336 (1989).
[CrossRef]

**26. **A. H. Hielscher, S. L. Jacques, L. Wang, and F. K. Tittel, “The influence of boundary conditions on the accuracy of diffusion theory in time-resolved reflectance spectroscopy of biological tissues,” Phys. Med. Biol. **40(11)**, 1957–1975 (1995).
[CrossRef]

**27. **S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol. **42(5)**, 841–853 (1997).
[CrossRef]

**28. **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(6)**, 2767–2772 (2000).
[CrossRef]

**29. **M. Cheikh, H. L. Nghiêm, D. Ettori, E. Tinet, S. Avrillier, and J. M. Tualle, “Time-resolved diffusing wave spectroscopy applied to dynamic heterogeneity imaging,” Opt. Lett. **31(15)**, 2311–2313 (2006).
[CrossRef]

**30. **A. Khireddine, K. Benmahammed, and W. Puech, “Digital image restoration by Wiener filter in 2D case,” Adv. Eng. Software **38(7)**, 513–516 (2007).
[CrossRef]

**31. **H. Niu, P. Guo, L. Ji, Q. Zhao, and T. Jiang, “Improving image quality of diffuse optical tomography with a projection-error-based adaptive regularization method,” Opt. Express **16(17)**, 12423–12434 (2008), http://www.opticsinfobase.org/abstract.cfm?URI=oe-16-17-12423. [CrossRef]

**32. **V. Chernomordik, R. Nossal, and A. H. Gandjbakhche, “Point spread functions of photons in time-resolved transillumination experiments using simple scaling arguments,” Med. Phys. **23(11)**, 1857–1861 (1996).
[CrossRef]

**33. **V. Chernomordik, A. Gandjbakhche, M. Lepore, R. Esposito, and I. Delfino, “Depth dependence of the analytical expression for the width of the point spread function (spatial resolution) in time-resolved transillumination,” J. Biomed. Opt. **6(4)**, 441–445 (2001).
[CrossRef]

**34. **S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express **16(7)**, 5048–5060 (2008), http://www.opticsinfobase.org/abstract.cfm?URI=oe-16-7-5048. [CrossRef]

**35. **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).
[CrossRef]