1. Introduction

Diffuse optical imaging is an emerging medical imaging modality which uses near-infrared (NIR) and visible red light to probe tissue physiology. Information about optical properties within the tissue being imaged is drawn from measurements of light that has traveled through the tissue.

In diffuse optical tomography, the 3D spatial distributions of optical properties (typically the scattering coefficient µ _s and the absorption coefficient µ _a), or changes in the optical properties are reconstructed. In brain activation imaging, reconstructed changes in µ _a at selected wavelengths are used to calculate changes in the concentrations of oxygenated (HbO2) and deoxygenated hemoglobin (HbR) [1, 2, 3]. The technique should not be confused with optical topography, in which optical measurements are used to produce a 2D image of HbO₂ and HbR changes in the activated brain areas.

In order to reconstruct changes in the optical properties of tissue, a model which predicts light propagation between the source and detector fibers is needed. This requires a computational model and a compatible anatomical model. Light propagation in tissue is thought to be accurately modeled by the Radiative Transfer Equation (RTE) [4]. The Monte Carlo method is a flexible, albeit computationally heavy, way of implementing the RTE [5]. The diffusion approximation (DA) to the RTE is considered fairly accurate in most tissues, and can be efficiently implemented using the finite element method [6]. Drawbacks of using the DA include the inability to correctly model light propagation in tissues with low scattering such as the cerebrospinal fluid (CSF), and inaccuracy close to the light sources, where the light propagation is not completely diffusive.

In brain activation imaging, linear methods are often used to reconstruct changes in optical absorption [7, 8]. In reconstructing changes in µ _a and µ _s using linear methods, it is assumed that the change in the optical properties is small and the background optical properties can be considered to be nearly time-invariant. An initial estimate of the distribution of optical properties is also needed when non-linear methods are used. A variety of background models have been used in the literature, with varying degrees of sophistication. They include homogeneous models with either simple geometric (for example semi-infinite slab [9]) or realistic surface shape [3], layered models [10] and realistic models based on MR images [7].

In our recent study [11], we found that accurate anatomical a priori information from MR imaging can significantly improve the spatial accuracy and robustness of reconstructing absorption changes due to brain activation in the neonatal human brain. Use of MR imaging based prior anatomical information has also been studied by other groups [7, 12, 13, 14]. However, MR image data of the subject is often not available, and creating a new MR imaging based model for each subject is time consuming. The latter is especially true in the case of infants, as automatic segmentation tools often perform poorly on infant MRI data. For this reason, we suggest that an atlas model which provides generic anatomic a priori information without requiring MR imaging of the subject could be a valuable aid in optical imaging of the infant brain. Several probabilistic atlases of the adult brain have been introduced (see e.g. [15, 16]), and they provide an effective way to use prior information of the location and variability of brain tissues and structures in image processing applications such as in brain segmentation [17]. These models, however, generally disregard the superficial tissues which are important for image reconstruction in optical imaging.

Registration of optical imaging data to a common stereotactic space such as the Montreal Neurological Institute (MNI) space has been proposed previously [18, 19, 20]. The methodology can be used to project optical imaging data measured at the surface of the head to a position in a brain atlas, which facilitates comparisons between subjects, and allows comparing the results of optical imaging to other modalities.

The use of generic anatomical information derived from a single individual other than the subject has previously been suggested for image reconstruction in optical brain activation imaging by others in the case of adult subjects [21], and by us in the case of infant subjects [23]. In the case of topographic (2D) reconstruction in adult subjects, use of a non-probabilistic model based on a hard segmentation of an averaged MR imaging based head atlas has been proposed [22]. However, no previously reported reconstruction methods utilize a probabilistic atlas in the calculation of the forward problem of diffuse optical imaging.

In this paper, we present to our knowledge the first implementation of a probabilistic head atlas intended for reconstructing activations in the infant brain from optical imaging data. We propose to use the atlas for solving the forward problem which is required for image reconstruction. Our atlas is based on a library of anatomical MR images of infants. It can be warped to correspond to the surface shape of the individual being imaged, and provides a map of the probability of each tissue type at each location inside the head. The atlas also retains the structural information from each MR image in the library, which can be exploited in calculating the forward problem (see Section 2.4). Construction of the atlas is presented, and its use in reconstructing brain activations is tested using Monte Carlo simulation. Monte Carlo was chosen as the computational model due to its ability to model arbitrary optical properties and its straightforward applicability to arbitrarily shaped anatomical structures.

2. Methods

2.1. Material

We used a set of seven anatomical MR images of heads of healthy full term neonates. The gestational ages of the infants at the time of imaging were between 39 weeks, 6 days and 41 weeks, 4 days (median 40 weeks). The data were originally collected for other studies at the Helsinki University Central Hospital (HUCH) and the Massachusetts General Hospital (MGH). In the HUCH data, the in-plane resolution was 0.94 mm and the slice thickness was 1.0 mm. For the MGH data, the in-plane resolution was 0.86 mm and the slice thickness 1.5 mm.

2.2. Construction of the probabilistic atlas

The MR images were segmented into scalp, skull, cerebrospinal fluid (CSF), and gray (GM) and white brain matter (WM). The segmentation was done by one of the authors (JH) using an in-house programmed interactive segmentation software, which features tools such as thresholding, region growth, dilation-erosion and painting by hand. Special care was taken to accurately segment the regions of CSF between the skull and the brain. This work-intensive approach was taken because results of automatic segmentation tools were not satisfactory.

In order to form the atlas, the MR images were coregistered together. This was done using only the outer surfaces and anatomical landmarks on the surface of the head. The internal structure of the head was not used in the coregistration. This approach will accurately match the head surface but not the internal structures. Since the purpose of the atlas is to provide generic knowledge of the internal structure of the head when the only surface shape is known, this result is desirable. Furthermore, when using the atlas, it has to be deformed to correspond to the head of the subject, for which only surface shape measurements are available. The anatomical landmarks considered which we found most useful in matching the orientations of the head models were the nasion and the right and left periauricular points.

The coregistration was done in two stages. First an affine transformation was used to set the orientation of the head models and the scale the same, and to match the anatomical landmark points. After this, elastic deformation was used to accurately match the outer shape of the head surface. Care was taken to preserve the thickness of the scalp, skull and the CSF layer, and the shape of the brain surface. The affine and elastic transformations of the MR images were performed using software developed at Helsinki University of Technology and the VTT Technical Research Centre of Finland [24]. Fig. 1 A-C shows the deformation of one head model to correspond to the surface shape of another. Fig. 1 also shows the anatomical variation between normal individuals, and illustrates how the shape of large pools of CSF varies between normal full term infants, even after the surface of the scalp and anatomical landmarks have been matched. This suggests that the probabilistic anatomical information provided by a probabilistic atlas may be more helpful than simply using the deformed head model of one subject to aid reconstructions from optical measurements conducted in another.

 

Fig. 1. Illustration of deformation between head models. In A and B, axial slices of segmented head models from two individuals. Due to difference in the orientations of the two models, the level of the slice in the axial direction is approximately the same in both cases in the anterior part of the head, and different in the posterior part of the head. In C, an axial slice of head model shown in B, coregistered and deformed to the surface shape of the model in A, at the same level as the slice in A. In D, an example of a layered model with same surface shape as the model shown in A.

Download Full Size | PDF

When creating the atlas, one MR image was selected as the target image to which the other images were coregistered. One of the images was left out of the atlas and used as the ’subject’ in simulated measurements of brain activations. The atlas thus consists of six MR images. Both the selection of the target image for creation of the atlas and the selection of the subject image were made randomly, in the beginning of the process.

2.3. Models

The performance of the atlas in reconstructing brain activations was compared to results obtained using simpler tissue models. These included a homogeneous model and layered models. The model which was used as the simulated ’subject’ is called the reference model. For completeness, reconstructions were also performed using the reference model.

All the anatomical models were deformed to correspond to the surface shape of the reference model. As when creating the atlas, only the surface shape and anatomical landmarks were used in deforming the atlas. We assumed accurate knowledge of the position of the anatomical landmarks and the surface shape of the reference model. All anatomical models were treated in the same way and had exactly the same surface shape. This allowed accurate comparison of the localization of the reconstructed perturbations obtained using the different anatomical models. Experimental errors which would in reality be present in the knowledge of the surface shape and positions of the optodes were not included in the simulation.

In Fig. 2, a transaxial slice of the reference and the probabilities of each of the five tissue classes in the atlas at the same transaxial level are shown. It can be seen that the surface-based coregistration results in a rather good coregistration of the brain. The CSF-filled Sylvian fissure between the frontal and temporal lobes coincides quite well in the different images of the atlas, as can be seen from the probability map for CSF. The reference model is also well, but not perfectly coregistered with the atlas.

The layered models consisted of five layers representing the scalp, the skull, the CSF, the GM and the WM. Three different layered models were created. In the first and the second models, the thickness of the tissue layers was based on atlas information. The tissue model was divided into thin shells, and the most common tissue type for each shell was calculated from the images in the atlas. Only the region underlying the optode array (see Section 2.6) was considered when calculating the most common tissue type. The CSF was not the most common tissue type in any of the shells. In the first layered model (named LAY), the shell with the most CSF voxels (between the skull and the GM) was set to CSF. In the second layered model (named LAYnoC) the CSF was not included. The third layered model (named LAYAh), we used ad hoc thicknesses for each tissue class. Fig. 1 D shows an axial slice of the LAY model. The thicknesses of the tissue layers in each layered model are given in Table 1. From the table it can be seen that the ad hoc and atlas based layered models differ considerably.

 

Fig. 2. A transaxial slice through the reference model, and probabilities for different tissue types of the atlas at the same level. To ease comparison between the reference model and the atlas, contours of the surface of the head and the CSF space are superimposed in red. Note that the surface shape for all models and all images in the atlas is the same.

Download Full Size | PDF

Table 1. Thickness of tissue layers^* in the layer models (mm). All models have central region of WM, thus no thickness is given for WM. LAY=atlas based layered model with CSF, LAYnoC=atlas based without CSF, LAYAh=ad hoc layered model.

View Table | View all tables in this article

It can be seen from the transaxial slices in Fig. 1 that there is large variation in the thickness of the scalp and the skull between different locations on the head. Under the optode grid used in this study, the thickness of the skull can vary from 1 mm to more than 4 mm. For CSF, the variation is even larger. This can cause problems when a model with tissue layers with fixed thickness is used. The probability maps of tissue types shown in Fig. 2 show that there is also variation between individuals, but this variation is smaller. This average thickness of the scalp and skull under the optode array calculated from the atlas were 2.7 mm and 2.7 mm, and for the reference model they were 2.4 mm and 2.2 mm, respectively. The figures given have been calculated from (voxel based) MR images, and are approximative.

The optical properties of different tissue classes in the reference model used for creating the simulated data were obtained from literature [25, 26]. We chose optical properties for light with wavelength 800 nm. The optical properties needed by the Monte Carlo model are the absorption and scattering coefficients µ _a and µ _s, the anisotropy factor g and the refractive index n. The reduced scattering coefficient is given by µ′_s=(1-g)µ _s.

When performing the reconstruction, optical parameters with the anatomical models were set in the following manner. The optical parameters used in the simulation are the reference set of parameters, called REF. The REF parameters were used with the anatomical reference model, the atlas, and all layered models. We also created a deviated set of optical parameters (called DEV). Error was introduced into the DEV optical parameters by sampling them from a normal distribution of the logarithm of µ _a and µ _s, centered at the correct (REF) parameters. Standard deviation for the logarithm of µ _a and µ _s were set at log(1.8). The DEV parameters were used for reconstruction using the atlas. We present reconstruction results for a representative DEV parameter set. Finally, we studied the case of a homogeneous model, with homogeneous parameters (named HOM). The REF, DEV and HOM parameter sets are tabulated in Table 2.

Table 2. Optical properties of tissue types [25, 26]. Parameter sets REF, DEV and HOM are given. The absorption (µ_a) and reduced scattering coefficients (µ′_s), vary between sets of optical parameters, while the anisotropy factor g and the refractive index n remain the same for all sets.

View Table | View all tables in this article

2.4. Perturbation Monte Carlo

We used the Monte Carlo (MC) method to generate the simulated measurements. The perturbation Monte Carlo (pMC) [28, 29] method was used for reconstructing the perturbations from the simulated measurement data. The simulations were carried out in time domain (TD). The data types considered were the signal intensity and mean time of flight of photons.

In the MC method, light propagation is modelled by tracing paths of individual photon packets as they travel through a tissue model. In pMC, the photon paths are also used to obtain three-dimensional sensitivity maps which provide a linear estimate of the change in the optical signal in response to a change in the optical parameters at a specific location. Rules of photon propagation are described in previous literature [5] and shall not be discussed in detail.

The atlas we describe in Section 2.2 consists of individual MR images, and probability of a tissue type at each location is defined by these constituent images. When used with Monte Carlo, it is possible to run each photon packet using one randomly selected constituent image of the atlas. This is the approach we selected in this study. The sensitivity maps thus obtained contain averaged information from the constituent images.

In the present study, we only consider perturbations in absorption. The motivation for this choice was that in literature of NIRS studies of brain activations, it is commonly assumed that changes occur mainly in absorption, and that the concentrations of oxyhemoglobin and deoxyhemoglobin can be estimated from changes in the absorption coefficient for different wavelengths [7, 30, 31]. The sensitivity to changes in optical absorption can be derived from the expressions for signal intensity W and photon mean flight time 〈t〉 as obtained from MC simulation. These are given by [11, 30]

where w_p and t_p are the intensity contribution and flight time of the photon packet p, and µ _a(r) and l _r,p are the absorption coefficient and the path length of the photon packet p in the region r of the tissue, respectively. The sum in the exponential function goes over all regions r of the tissue model. In our implementation, the tissue model is divided into voxels.

The sensitivity of the signal to a change in µ _a in a region ρ is obtained by differentiation. We obtain [11, 30, 32]

where 〈l_ρ〉 is the average weighted photon pathlength in the region ρ and 〈l_ρt〉 is the average pathlength of the photon packets in the region ρ, weighted with the total flight time and the final weight of the photon packets at the detector.

All quantities in Eqs. 3 and 4 can be accumulated by MC simulation to obtain the derivatives required for the pMC sensitivity maps. With these maps, absorption changes can reconstructed without repeatedly calculating the computationally costly MC forward solution. In using stored maps of sensitivity, we make the approximation that the sensitivity profile is not significantly altered by the changes induced to the optical properties. Our previous results indicate that this linear approximation is sufficiently accurate for image reconstruction in the case of localized µ _a changes such as those studied in this paper [11]. Work by Liebert et al. suggests that linear methods may also be appropriate for estimating absorption changes over larger volumes in the brain [8].

The reconstruction algorithms consists of minimizing the objective function

where the vector Δµ⃗_a represents the change in µ _a in different voxels of the imaged volume.

A(Δµ⃗_a)=σ^-2 _s,d∑_s,d(f _s,d[Δµ⃗_a]-M _s,d) is the sum over all source-detector pairs (s,d) of the squared difference between model-predicted (f _s,d[Δµ⃗_a) and measured (M _s,d) difference data, weighted by the reciprocal of the variance σ _s,d in the measured (simulated) signal. Although all measurements are used, the weighting reduces the significance of measurements with long source-detector distances. Initial error in both signal intensity and mean time of flight were scaled to one, giving equal weight to both data types.

The functions R ₁(Δµ⃗_a) and R ₂(Δµ⃗_a) are regularization functions. R ₁(Δµ⃗_a)‖Δµ⃗_a‖² ₂ is simply the norm of the solution vector, or Tikhonov regularization. The second regularization function is given by R ₂=‖LΔµ⃗_a‖² ₂, where L is defined by

$L_{ij}=\{\begin{array}{cc}\multicolumn{1}{c}{nn,}& \multicolumn{1}{c}{\mathrm{if}i=j}\\ \multicolumn{1}{c}{-1,}& \multicolumn{1}{c}{\mathrm{if}j\mathrm{is}a\mathrm{neighbour}\mathrm{of}i,}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{\mathrm{otherwise}}\end{array}$

nn is the number of neighbours of voxel i (26-neighbourhood was used). It is used to reduce the high frequency noise present in the unregularized reconstructions. γ ₁ and γ ₂ are regularization parameters. Since optimal regularization depends on the contrast-to-noise ratio of the measured signal, they were chosen empirically for the four different scenarios (with varying number and contrast of absorptive perturbations) presented in Section 2.6. For each scenario, the same regularization parameters were used with all anatomical models.

We performed the reconstructions with an iterative algorithm using a line search. The search direction was the gradient of Φ with respect to µ⃗_a, divided by the diagonal of the Hessian matrix (thus elements of the search direction are given by $\genfrac{}{}{0.1ex}{}{\partial A\left(\Delta \overrightarrow{{\mu }_{\text{a}}}\right)}{\partial \overrightarrow{{\mu }_{\text{a}}}\left(r\right)}⁄\genfrac{}{}{0.1ex}{}{{\partial }^{2}A\left(\Delta \overrightarrow{{\mu }_{\text{a}}}\right)}{{\partial \overrightarrow{{\mu }_{\text{a}}}\left(r\right)}^2}$. The diagonal elements of the Hessian matrix are obtained by differentiating Eqs. 1 and 2 a second time, and necessary terms can be accumulated as in the case of the first derivatives (Eqs. 3 and 4).

2.5. Methods for assessing quality of reconstructions

The quality of reconstructions obtained using the different models for the optical background were assessed using qualitative assessment of images drawn from the reconstructions, and quantitative parameters calculated from the reconstructions. Both in the qualitative and the quantitative appraisal, the reconstructed perturbations are compared to the simulated perturbations, called the target perturbations.

Images were drawn of axial slices through the center of the target perturbations. In addition, tangential maps [11] representing a 2D distribution of µ _a changes on the cortex are shown. The tangential maps were generated by calculating for each voxel of the tissue surface the mean value of the 3D reconstruction data in the direction normal to the surface of the head between selected depth limits. The superficial and deep limits of this calculation were set to approximately correspond to the superficial and deep ends of the simulated perturbations. Due to the fact that the reconstructed perturbations are located more superficially than the target perturbation in all cases (as shown in the axial slices), the superficial limit was set 1 mm and the deep end 2 mm more superficially than the superficial and deep ends of the target perturbation, respectively. For reconstructions from multiple perturbations, contrast profiles are presented. The purpose of these profiles is to highlight how the combination of three target perturbations is reconstructed using the different background models. The contrast profiles were calculated as the mean contrast in the tangential maps defined above in a band which goes through the target perturbations, and has the width of the target perturbations.

The quantitative parameters calculated from the reconstruction results were the localization error of the perturbations, the peak contrast (PC), the integrated contrast (IC), and the contrast-to-noise ratio (CNR). The localization error was calculated as the difference between the position of the center-of-mass of the reconstructed perturbation and the target perturbation. The continuous region around the peak of the reconstructed perturbation with reconstructed contrast above 60 % of the peak value was considered.

The peak contrast was defined as the highest contrast within the region of the target perturbation, thus excluding erroneous peaks close to tissue surface. The integrated contrast was defined as the volume integral of reconstructed contrast within the region of the target perturbation. Both the peak and integrated contrasts are given as percentage of the values calculated from the target perturbations.

The CNR was calculated as the mean contrast within the region of the target perturbation, divided by the standard deviation of the background region. The background was defined as the volume below the optodes to a depth of 5 cm, excluding the regions of the target perturbations enlarged by increasing their original radius by a factor of 1.5 to allow for the spatial blurring by the diffuse imaging method [11].

2.6. Simulations

We simulated imaging of localized activations in the neonatal brain. The brain activations were located in the inferior part of the frontal lobe and the superior part of the temporal lobe in a region which contains such functionally important areas as Broca’s area and the auditory cortex. The brain activations were represented by absorption changes in the cerebral cortex. The perturbations were restricted to brain tissue, and therefore regions of CSF altered their otherwise spherical form.

The cases of a single absorptive perturbation and three absorptive perturbations were studied. The single perturbation was placed near to a region of CSF, at depth of approximately 10 mm from tissue surface (center of the perturbation). Of the three absorptive perturbations, the most anterior one (perturbation A) was set apart and deeper from tissue surface than the two more posterior perturbations (B and C), which were close to each other and sligthly more superficial. The depth (distance from tissue surface to center of perturbation) was 12.5 mm for perturbation A, and 9 mm for perturbations B and C. Perturbation A had a radius of 4.5 mm, and perturbations B and C had radii 4.2 mm.

In both cases presented above, perturbations with contrasts 20% and 11% above the background absorption coefficient in the GM were considered. The higher contrast case is slightly higher than would be expected in a real measurement, while the smaller contrast corresponds to what might be expected during brain activation, based on animal studies [33]. Two different contrasts with same noise level were used to illustrate the effect of contrast-to-noise ratio on the reconstructed image quality. The noise (standard deviation) in the simulated signal intensity and photon flight time were approximately 0.2 % and 0.2 ps at 12.5 mm source-detector separation and 0.7 % and 1.4 ps at 25 mm. Only photon noise was considered. The noise-level is similar to what may be obtained in a practical measurement with a near-perfect instrument [11].

A dense optode grid of 17 source optodes and 18 detector optodes was used (see Fig. 3). The shortest inter-optode distance in the optode grid is 8.6 mm. The choice of optode grid was motivated by our previous study, which suggested that a dense optode array should be used to achieve good spatial resolution [11].

The performance of the atlas and other models (see Section 2.3) were compared in reconstructing the simulated brain activations.

3. Results

3.1. Case of one absorptive perturbation

The results from the high contrast case are shown in Fig. 4, and for the low contrast case in Fig. 5.

 

Fig. 3. Arrangement of optodes. The 17 crosses stand for light sources, the 18 circles for detectors.

Download Full Size | PDF

 

Fig. 4. Results from reconstructing a single absorptive perturbation, high contrast case. In A, axial slices through center of the target perturbation. In B, tangential maps. On top row, the target perturbation, and the reconstruction obtained with the reference model. On second row, reconstructions obtained with the atlas using reference (left) and deviated (right) optical parameters. On third row, reconstructions obtained using the atlas based layered models with a CSF layer (left) and without a CSF layer (right). On fourth row, reconstructions obtained with the ’ad hoc’ layered model and the homogeneous model.

Download Full Size | PDF

 

Fig. 5. Results from reconstructing a single absorptive perturbation, low contrast case. In A, axial slices through center of the target perturbation. In B, tangential maps. Cases as in Fig. 4.

Download Full Size | PDF

The localization error, recovered contrast (PC and IC) and the CNR of the reconstructions in the case of a single perturbation are shown in Table 3.

For the high contrast case, the reconstruction was succesful using all the anatomical background models. The best localization accuracy was obtained using the reference model and the atlas. The results obtained with the atlas and the reference (REF) and deviated (DEV) optical parameters were similar. The layered models with a CSF layer (both atlas based and ’ad hoc’) produced good reconstructions, with a slight deterioration in image quality compared to the atlas and reference models. The reconstruction using the homogeneous model was the least succesful.

In the low contrast case, in which the change in signal due to the simulated perturbation was close to the noise level, the reference and atlas models performed well, whereas the results obtained using the layer models were less succesfull. The ’ad hoc’ layered model (LAYAh) performed better than the atlas based layer models. The reconstruction using the homogeneous model failed.

3.2. Case of three absorptive perturbations

The results for the case with high simulated contrast are shown in Fig. 6, and for the case with low simulated contrast case in Fig. 7. One should note that due to the slightly deeper location of perturbation A compared to perturbations B and C necessitated the use of different limits for the depth at which the tangential maps were calculated for the anterior (left) and posterior (right) parts. The width of the depth range was kept constant, and the depth limits were gently increased from the anterior part of perturbation B to the posterior part of perturbation A. In Fig. 8, contrast profiles (see Section 2.5) for the high and low contrast cases are shown superimposed on the profile calculated from an image of the target perturbations.

Due to the fusing of adjacent perturbations in the reconstructions, localization errors could not be calculated for the case of three perturbations. The contrasts PC and IC and CNR are

 

Fig. 6. Results from reconstructing three absorptive perturbations, high contrast case. In A, axial slices through center of the target perturbation. In B, tangential maps. On top row, the target perturbations, and the reconstruction obtained with the reference model. On second row, reconstructions obtained with the atlas using reference (left) and deviated (right) optical parameters. On third row, reconstructions obtained using the atlas based layered models with a CSF layer (left) and without a CSF layer (right). On fourth row, reconstructions obtained with the ’ad hoc’ layered model and the homogeneous model.

Download Full Size | PDF

 

Fig. 7. Results from reconstructing three absorptive perturbations, low contrast case. In A, axial slices through center of the target perturbation. In B, tangential maps. Cases as in Fig. 6.

Download Full Size | PDF

Table 3. Localization error, contrast and contrast-to-noise ratio for reconstructions from the single perturbation, with the high and low target contrasts. PC=Peak Contrast, IC=Integrated Contrast, CNR=Contrast-to-Noise Ratio. Contrasts given as percentrage of target contrast. On the row defining anatomical model: LAY=Atlas based layered model, LAYnoC=Atlas based layered model, no CSF, LAYAh=Ad hoc layered model. On the row defining optical parameters: REF=reference parameters, DEV=deviated parameters, HOM=homogeneous parameters.

View Table | View all tables in this article

given in Table 4.

Table 4. Localization error, contrast and contrast-to-noise ratio for reconstructions from three separate perturbations, with the high and low target contrasts. PC=Peak Contrast, IC=Integrated Contrast, CNR=Contrast-to-Noise Ratio. Contrasts given as percentrage of target contrast. On the row defining anatomical model: LAY=Atlas based layered model, LAYnoC=Atlas based layered model, no CSF, LAYAh=Ad hoc layered model. On the row defining optical parameters: REF=reference parameters, DEV=deviated parameters, HOM=homogeneous parameters. A, B and C refer to the target perturbations from posterior to anterior (left to right in Figs. 4 through 7.

View Table | View all tables in this article

The perturbations B and C dominated the reconstructions from the case with three perturbations. The perturbation A was reconstructed with reduced contrast, and in some cases appeared as a ’tail’ in the reconstruction, without appearing as a peak of its own. As can be seen in the Fig. 6, the result obtained with the ’ad hoc’ layered model shows a very strong ’tail’ that continues past the anterior perturbation. While similar tail also appears in other reconstructions, its contrast beyond the anterior perturbation is low.

As shown in Figs. 6 and 7, the contrast recovery for perturbations B and C was good using all except the homogeneous model, and they were fused together in all cases. As seen in the profile graphs in Fig. 8, the atlas models provide a more balanced symmetrical reconstruction of these two perturbations while the layered and homogeneous models reconstruct the central perturbation B as larger in contrast than perturbation C, especially in the high contrast case. One can also see, that the deeper perturbation A appears as a clear increase or plateau in contrast in the reconstructions obtained with the reference and atlas models. In reconstructions obtained with the layered and homogeneous models, the perturbation A is either reconstructed as a ’tail’ in contrast sloping downwards and left (layered models with CSF layer, both atlas based and ’ad hoc’), has very low contrast (layered model with no CSF layer) or both (homogeneous model). In reconstructions using the homogeneous model, perturbation A can be seen as separate, but it is reconstructed too posteriorly (to the right in Figs. 6 through 8). Thus the atlas with the reference (REF) optical parameters provided superior reconstruction quality in the cases with three perturbations, and the performance of the the atlas with the deviated (DEV) optical properties was also good.

 

Fig. 8. Contrast profiles calculated from the tangential maps in the case of three absorptive perturbations. On the y-axis, mean contrast in the band containing the target perturbations (10^-3 mm^-1). On the x-axis, location in mm. Cases as in Fig. 6.

Download Full Size | PDF

For optimal reconstruction results, the high an low contrast cases required different levels of regularization. In order to allow good contrast recovery for the high contrast perturbations, the Tikhonov regularization had to be reduced, which increased the background noise in the reconstructions. For this reason, the CNR values for the low contrast case are similar to and in some cases superior to the CNR values for the high contrast case.

4. Discussion

We have explored the feasibility of using atlas based anatomical a priori information to improve reconstruction of optical activation images in infants.

Our results from the case with a single absorptive perturbation show that localization accuracy can be improved by the use of the atlas (see Table 3). Our results also show that more accurate and sophisticated optical background models make the reconstruction more resistant to the effects of noise, and give more accurate and consistent results.

The layered model based on the atlas worked well in many cases, although often the ’ad hoc’ atlas model produced higher contrast. Reconstructions obtained using the homogeneous model had poor contrast and CNR, and poor localization accuracy in the depth direction, but in most cases quality of the image was reasonably good in the tangential direction.

In this study, errors in defining the surface shape of the head and position of the optodes were not taken into account. These error sources may weaken the behaviour of the atlas. However, we believe that the smooth probabilistic nature of the anatomical information of the atlas makes it resistant to these problems to some degree.

A challenge in using atlas based anatomical models is finding the correct optical parameters to be used with the model. While we found incorrect parameter values in the atlas model to reduce the quality of reconstructions, the results were better than those obtained using a simplified models. The background optical properties can at least in theory be determined using nonlinear image reconstruction techniques. However, currently nonlinear Monte Carlo may be computationally too heavy for this.

While the assumption that the human head is made of a number of tissue types with clearly defined optical parameters is itself an approximation, we believe that a model consisting of such tissue types can be made to correspond optically to the physical situation with sufficient accuracy. The probabilistic nature of the atlas model will probably help in making such a match.

5. Conclusions

We conclude that the use of probabilistic atlas information to aid in the reconstruction of brain activation images from diffuse optical imaging data is feasible. Our results show that when individual anatomical information is not available, generic anatomic information can be more helpful than a simple layered model of the head, even with careful selection of the thickness of the tissue layers. Our atlas used in this study consisted of only six subjects, which was found sufficient to demonstrate the concept. Adding MR images of more healthy subjects will make the atlas more representative of the subject population.

If using an anatomical atlas for reconstruction of brain activations is not possible, our results suggest that using a layered model may be advantageous. However, care should be taken when selecting the thickness of the tissue layers.

In this study, all the neonates used were normal subjects. In clinical situations, the brain may be abnormal and this may also be reflected in the anatomy. Therefore, the use of the probabilistic atlas to support reconstruction in abnormal cases is not advised and at the very least requires further study. Use of the atlas in case of subjects whose anatomy is normal but represents an extreme of the normal variation should also be studied further.

Future challenges for the use and development of atlas models for reconstructing brain activation images from optical imaging data include developing accurate and convenient methods for defining the surface shape of a subject and positions of anatomical landmarks and optodes, methods for determining the optical parameters of tissue classes using optical measurements, and creating atlases based on a larger set of MR images for different age groups.