Compensation of optode sensitivity and position errors in diffuse optical tomography using the approximation error approach.

Diffuse optical tomography is highly sensitive to measurement and modeling errors. Errors in the source and detector coupling and positions can cause significant artifacts in the reconstructed images. Recently the approximation error theory has been proposed to handle modeling errors. In this article, we investigate the feasibility of the approximation error approach to compensate for modeling errors due to inaccurately known optode locations and coupling coefficients. The approach is evaluated with simulations. The results show that the approximation error method can be used to recover from artifacts in reconstructed images due to optode coupling and position errors.


Introduction
Diffuse optical tomography (DOT) is a non-invasive soft tissue imaging modality with applications for example in breast and brain imaging [1][2][3]. It involves the use of near-infrared light for probing and determining spatially distributed optical parameters. The image reconstruction problem in DOT is ill-posed. The ill-posedness means that even small errors in measurements or modeling can cause large errors in the reconstructions. In contact based DOT measurement setup, lasers are coupled to fiber optodes and these are used to direct light into the tissue. The transmitted light is collected using another set of optodes coupled to photo-multipliers. In practical experiments, the positions of the source and detector optodes are not always known accurately. There are also other uncertain model parameters such as coupling coefficients that include source strength, coupling losses of the optodes and detector efficiency or gain. It has previously been observed that small errors in source and detector positions and coupling coefficients can cause large artifacts in the reconstructed images [4,5].
Various approaches for reduction of errors caused by inaccurately known optode coupling have been developed. Difference imaging can reduce measurement artifacts to an extent in the reconstructed images [6]. Difference imaging requires a reference measurement with known optical properties and provides information about changes in the optical properties of the target between the reference and measurement states. However, this technique only works when the changes in the optical properties are not very large. A modified technique based on the same idea that differential changes are less susceptible to surface errors was developed in [7]. This technique uses a contrast agent that needs to be injected during imaging which may not always be appropriate. A similar technique [8] uses the differences between measurements at two separate wavelengths to reduce the coupling errors. A calibration method using a homogeneous diffusive phantom with source placed at the center and optodes placed equidistant from the center to calibrate relative optode sensitivities was used in [6]. Furthermore, a calibration method based on rotational symmetry of source and detector positions in a measurement setup was utilized in [5,9,10]. Difference imaging and calibration approaches assume that the source and detector coupling coefficients do not change after the calibration. This may not always be a realistic assumption since coupling coefficients are dependent on the statistical fluctuations of the source and detectors and may also change due to movement of the patient or the optode fibers. Therefore, approaches in which the errors caused by coupling coefficients are handled during the image reconstruction have been developed. A regularizing functional which suppresses variation of the optical parameters near the boundary was utilized in [11]. The weighing function method uses priori constraints at image boundary without actual evidence from data. Simultaneous image reconstruction and optode calibration was performed in [12] and was further extended to include phase errors in experiments [13]. An extension to time domain measurements can be found in [14]. This method solves the optode coupling coefficients along with the optical parameters. A Bayesian approach to solve the coupling coefficients along with the optical coefficients was used in [15].
While performing corrections for source and detector coupling related errors, the position errors are usually not considered. A calibration technique where the signal distortions due to optode position errors were approximated as coupling coefficients was presented in [4,16]. However, the effects due to position shifts are different from those of coupling losses and hence this approximation may only work for small distortions. To our knowledge modeling and estimation of optode position errors has not been performed.
Recently a theory known as the approximation error (AE) theory [17,18] has been developed to handle modeling errors. The approach is based on the Bayesian inversion paradigm, where all the unknowns including primary unknowns (absorption, scatter) as well as the unknown nuisance parameters such as the optode positions and coupling coefficients are modeled as random variables. The key idea in the approximation error approach is to represent the modeling errors caused by uncertainly known nuisance parameters as an additive modeling error noise in the observation model. The realization of the modeling error is obviously unknown and can not be computed without knowing the realizations of the unknowns and nuisance parameters. However, given the prior probability density models of all the uncertainly known parameters, one can carry out an approximate marginalization of the problem over the nuisance parameters and model reduction errors by using a Gaussian approximation for the joint statistics of the primary unknowns and the modelling error noise [19]. In DOT, the approximate marginalization by the approximation error approach has proven successful in recovery from modeling errors caused by coarse finite element discretization [17,20], geometric mis-modeling and domain truncation [21][22][23], using approximate physical models [24,25], and treatment of scatter as fixed nuisance parameter while estimating absorption image only [19].
In this work we propose the approximation error approach for recovery from modeling errors caused by inaccurately known source and detector coupling coefficients and locations. The remainder of the paper is organized as follows. In Section 2, we review the light transport model we use and describe the modeling of the coupling coefficients. We also review the usual image reconstruction using the conventional measurement error model and then describe how the approximation error approach can be used to compensate for the uncertainty in the optode coupling and locations. Results using simulated measurement data are presented in Section 3. Finally the conclusions are given in Section 4.

Diffusion approximation model
In a typical DOT measurement setup visible/near infra-red light is injected to an object from the object surface. Let Ω ⊂ R n , where n is the dimension of the medium (n = 2, 3), model this object domain. In a diffusive medium like soft tissue, the commonly used light transport model for DOT is the diffusion approximation (DA) for the radiative transport equation (RTE) [26]. In this paper the frequency domain version of the diffusion approximation model is used as the model for light propagation in tissues [27] where Φ i (r) := Φ i is the photon density for the i:th source, µ a (r) := µ a is the absorption coefficient, i is the imaginary unit, ω is the angular modulation frequency of the input signal, c is the speed of light in the medium and q 0,i (r) := q 0 is the source within object domain Ω operating at frequency ω. κ(r) := κ is the diffusion coefficient. The diffusion coefficient κ is given by κ(r) = 1/(n(µ a (r) + µ s (r))), where µ s (r) := µ s is the reduced scattering coefficient.
Let us for now assume that all sources and detectors are ideally coupled, i.e. there are no unknown amplitude loss or phase delays occurring in the source and detector optodes. Let us denote the location of the source optodes by m i ⊂ ∂ Ω, i = 1 . . . N s and the location of detector optodes by n j ⊂ ∂ Ω, j = 1 . . . N d . When the source is modeled as diffuse boundary source, q 0,i (r) = 0 and the source term can be written in the boundary condition wheren is the outward normal to the boundary at point r, q i is the strength of the boundary source at location m i , γ is dimension dependent constant (γ = 1/π when Ω ⊂ R 2 , γ = 1/2 when Ω ⊂ R 3 ) and α is a parameter governing the internal reflection at the boundary ∂ Ω. The measurable quantity, the exitance, Γ i, j at detector j under illumination from source i is defined by Let Γ ∈ C N s N d denote vector of complex valued measurement data corresponding to measurement between all source-detector pairs i, j with single indexation A typical frequency domain DOT measurement setup collects the amplitude and phase as measurement data where y ∈ R 2N s N d is data vector that contains the measured log amplitude and phase for all source-detector pairs. The source and detector locations m i and n j are surface patches of known length in 2D and area in 3D. We parameterize the locations by the center point of the source and detector optodes and use notation for the vector of source and detector location parameters. Using this notation, the observation model is where e ∈ R 2N s N d models the random noise in measurements, x = (µ a , µ s ) T ∈ R 2N n is discretized optical coefficients and the mapping A is typically based on the finite element method (FEM) solution of Eq. (1-3).

Source and detector coupling coefficients
Consider now a practical situation where there are coupling losses in the source and detector optodes. Following [13], we model the coupling losses in source q i in (2) by a complex valued multiplicative coupling coefficientŝ i ∈ C, leading to photon densitỹ whereŝ i = s i exp(iδ i ) is the source coupling coefficient with amplitude factor s i and phase factor δ i . Similarly the coupling losses in measurement optodes are modeled with multiplicative coupling coefficientŝ leading to exitance [13]: Let and define vector valued mapping g(ζ ) ∈ C N d N s such that Using these notations and taking the log transform of vector of data of the form (7) leads to the separation of coupling coefficients into additive components leading to observation model where ε 1 (ζ ) ∈ R 2N s N d is the discrepancy in the log transformed measurement compared to the ideal (no coupling losses) model (5) for given realization ζ . Notice that when ideal sources and detectors are assumed (no losses), we have s i = 1, δ i = 0 for all i, d j = 1, η j = 0 for all j and ε 1 ≡ 0, i.e., model (10) becomes equal to (5).

Posterior model
In the Bayesian approach to inverse problems, the principle is that all unknowns and measured quantities are considered as random variables and the uncertainty of their values are encoded into their probability distribution models [17,18,20,23]. The information on the unknown parameters (and the related uncertainty) given the measurements is embedded in the posterior model, given by the Bayes theorem as Here π(y|x, ζ , ξ , e) is the likelihood density modeling the probability of the measurement y when the realizations of x, ζ , ξ and e are given. π(x, ζ , ξ , e) is the prior density model and it models our information on the unknown parameters before the actual measurements.
To handle the uncertainty related to the unknown but uninteresting measurement errors e, the posterior (11) is usually marginalized with respect e as π(x, ζ , ξ |y) = π(x, ζ , ξ , e|y)de.
For details how the marginalization integration is obtained analytically in the case of the additive error model y = A(x) + e with Gaussian statistics for e, see section 2.3 in [19]. The posterior model π(x, ζ , ξ |y) is a probability density on a high-dimensional space. In computation of point estimate(s) from the posterior model, the most common choice is the maximum a posteriori (MAP) estimate. In principle, one could attempt to compute MAP estimate for all the unknown model parameters However, this would lead to computationally extensive and complicated problem, and to our knowledge the simultaneous estimation of (x, ξ , ζ ) has not been performed (for estimation of x and coupling coefficients ζ , see [12,13,15]). Alternatively, one could treat the uncertainty in the values of nuisance parameters (ξ , ζ ) in a similar manner than treating the uncertainty in e, that is, by marginalizing the posterior model as and then compute estimate for the primary unknowns from the posterior π(x|y). However, the solution of the integration (14) has closed form solution only in case of purely linear and Gaussian models. In case of non-linear problems such as DOT, the solution of (14) would require Markov chain Monte Carlo integration which would in most cases be infeasible for practical applications.
The key idea in the approximation error approach is to find approximationπ(x|y) for the posterior model (14) such that the marginalization over the uncertainty in the values of (ξ , ζ ) is carried out approximately but in a computationally feasible way.
Before introducing the approximation error approach for treating the uncertainty in the optode coupling and location parameters (ξ , ζ ), we first review the standard DOT reconstruction approach where ξ = ξ 0 and ζ = ζ 0 are treated as known and fixed conditioning variables. In most of DOT reconstruction schemes, the optode parameters ζ and ξ are treated as known deterministic parameters and independent of the optical coefficients. Within the Bayesian framework, any known parameter (whether measured or otherwise known) is interpreted as a conditioning variable. In the case at hand, if we take ζ and ξ to be known, we would have x and e as the only unknowns, and would obtain the conditional density Given the observation model (10) with ξ = ξ 0 and ζ = ζ 0 the observation model becomes y = A(x, ξ 0 ) + ε 1 (ζ 0 ) + e. Using Gaussian prior models for the unknown optical parameters x and the random measurement noise e wherex ∈ R 2N n andē ∈ R 2N s N d are the means, and Γ x ∈ R 2N n ×2N n and Γ e ∈ R 2N s N d ×2N s N d are the covariance matrices, and marginalizing over the unknown measurement errors e as π(x|y, ζ = ζ 0 , ξ = ξ 0 ) = π(x, e|y, ζ = ζ 0 , ξ = ξ 0 )de the posterior model becomes [19] π(x|y, The MAP estimate corresponding to the posterior (17) is obtained as where the Cholesky factors are L T x L x = Γ −1 x and L T e L e = Γ −1 e . If the coupling coefficients ζ 0 are not estimated by the pre-calibration techniques [5,6,9,10], they are typically fixed to the ideal "no losses" values, i.e., ε 1 (ζ 0 ) ≡ 0. In addition, if the random measurement errors have zero mean (ē = 0), the MAP estimate (18) has the usual form that is similar to the Tikhonov regularized least squares estimation, except that here the terms L e , L x andx would be derived from the related probabilistic models. We refer to the solution of (19) as the MAP estimate with the conventional measurement error model (CEM) approach.

Approximation error model
Consider now a practical situation where the detector and source locations and coupling coefficient may not be accurately known. Let ξ 0 and ζ 0 be the fixed realizations of the optode locations and coupling parameters that are to be used in estimation of the optical properties x. Evidently, if these fixed values are erroneous, these errors will lead to artefacts in the estimate of x.
To recover from these errors, we write the accurate measurement model (10) as where ε 1 and ε 2 are approximation errors that describe the discrepancy between the accurate model and the target model in which the optode parameters have the fixed values ζ 0 and ξ 0 . The measurement model (20) is called the approximation error model. In addition to marginalizing the problem over the uninteresting and unknown measurement noise e, the objective in the approximation error approach is to use measurement model (20) and treat the uncertainty related to the values of (ξ , ζ ) by carrying out an approximate marginalization of the posterior over the noise processes (ε 1 , ε 2 ). Following [19], we make Gaussian approximations for the density of ε 1 (ζ ) and the joint density of (x, ε 2 (x, ξ )) and denote the total additive noise by Further, if we make a technical approximation that x and ε 2 are mutually uncorrelated, we get an approximation that is referred to as the enhanced error model [17,18] π(x|y) ∝ exp{− 1 2 where the meann and covariance Γ n arē n =ε 1 +ε 2 +ē, Γ n = Γ ε 1 + Γ ε 2 + Γ e (22) andε 1 ,ε 2 and Γ ε 1 , Γ ε 2 are the means and covariances of the Gaussian approximations The MAP estimate with the enhanced error model is obtained as Where L T n L n = Γ −1 n . In the following sections, we refer to the solution of (24) as the MAP estimate with the approximation error model (AEM) approach.
Notice that the estimate (24) leads to a similar form of minimization problem than the MAP estimate with the conventional error model (19) (only the noise mean and Cholesky factorization of the noise precision matrix are different due to the different overall error model). Thus, the estimate for x can be efficiently computed by using the existing optimization codes for conventional measurement error model or Tikhonov regularized least squares.

Estimation of approximation error statistics
In the computation of MAP estimate (24), one has to have estimates for the meansε 1 andε 2 and the covariances Γ ε 1 and Γ ε 2 in equation (22). Closed form solutions for these are only available in purely linear and Gaussian case. In case of non-linear models, the means and covariances of ε 1 and ε 2 have to be estimated numerically by a simple Monte Carlo integration procedure. The following gives the outline of this estimation procedure.
2.5.1. Estimation ofε 1 and Γ ε 1 For the estimation of the mean and covariance for the optode coupling approximation error ε 1 (ζ ), we specify prior models π(s) and π(δ ) for the vectors of amplitude and phase coupling coefficients of the sources and prior models π(d) and π(η) for the amplitude and phase coupling coefficients of the detectors, respectively. The prior models are used for drawing sets of 1 := ε 1 (ζ (ℓ) ) by equations (8)(9)) and estimate the mean and covariance of ε 1 as from prior models π(x) and π(ξ ) = π(m)π(n). The samples are used to generate samples of ε 2 as ε (ℓ) The mean and covariances are estimated as Notice that, whereas the estimation of mean and covariance of ε 1 is computationally fast since it requires only evaluations of the mapping g(ζ ) (equations (8-9)), the estimation of the mean and covariance of ε 2 is computationally somewhat intensive as 2M forward solutions need to be evaluated. However, these computations need to be done only once for a fixed measurement setup and this estimation can be done off-line.

Simulation of the measurement data
In the numerical studies, the domain Ω ⊂ R 2 was a disk with radius r = 25 mm. The measurement setup consisted of N s = 16 sources and N d = 16 detectors. The source and detector optodes were modeled as 1 mm wide surface patches located at equi-spaced angular intervals on the boundary ∂ Ω. With this setup, the vector of DOT measurements (4) was y ∈ R 512 . A target with background optical properties µ a = 0.01mm −1 , µ s = 1mm −1 containing an absorption inclusion with µ a = 0.02mm −1 and scatter inclusion with µ s = 2mm −1 was constructed. The simulated measurement data was generated using FE approximation of the DA in a mesh with 33806 nodes and 67098 triangular elements. Random measurement noise e, that was drawn from a zero-mean Gaussian distribution where the standard deviations σ e,k were specified as 1% of the absolute value of simulated noise free measurement data, was added to the simulated measurement data.
In the computation of the MAP estimates (19) and (24) a FE mesh with 26075 nodes and 51636 elements was used. The MAP estimation problems were solved by a Gauss-Newton algorithm with an explicit line search algorithm [28]. The meanē = 0 and covariance Γ e were assumed known.

Prior models and computation of approximation error statistics
A proper (integrable) Gaussian smoothness prior was used as the prior model for the primary unknowns x = (µ a , µ s ) T . The same prior model was used both in the construction of the enhanced error model for ε 2 and all the MAP estimates based on the conventional measurement error model (19) and the approximation error model (24).
The absorption and scatter images µ a and µ s were modeled as mutually independent Gaussian random fields with a joint prior model In the construction of the mean vectorsμ a ,μ s and covariances Γ µ a and Γ µ s , the random field, say f (i.e., either µ a or µ s ), is considered in the form where f in is a spatially inhomogeneous parameter with zero mean, and f bg is a spatially constant (background) parameter with non-zero mean. For the latter, we can write f bg = qI, where I ∈ R N n is a vector of ones and q is a scalar random variable with distribution q ∼ N ( f * , σ 2 bg, f ). In the construction of Γ in, f , the approximate correlation length can be adjusted to match the size of the expected inhomogeneities and the marginal variances of f k :s are tuned based on the expected contrast of the inclusions. We model the distributions f in and f bg as mutually independent, that is, the background is mutually independent with the inhomogeneities. Thus, we havef [17,19,20] for further details, and see [29] for an alternative construction of a proper smoothness prior. The parameters in the prior model π(x) were selected as follows. The mean for background absorption and scatter were set as µ a, * = 0.01mm −1 and µ s, * = 1mm −1 and the standard deviations σ bg,µ a and σ bg,µ s of the background values were chosen such that 2 s.t.d. limits equaled 25% of the mean values µ a, * and µ s, * . In the construction of Γ in, f the correlation length for both µ a and µ s was set as 16 mm. The marginal standard deviations were set to equal values in each pixel and σ in,µ a and σ in,µ s were chosen such that 2 s.t.d. limits equaled 50% of the mean values µ a, * and µ s, * . Thus, the overall marginal standard deviations ( i.e., square root of diagonal elements of Γ µ a and Γ µ s ) were such that 2σ µ a = 0.0056mm −1 and 2σ µ s = 0.56mm −1 . This gives overall 2 s.t.d. intervals µ a ∈ [0.0044, 0.0156]mm −1 and µ s ∈ [0.44, 1.56]mm −1 , i.e., the values of absorption and scatter are expected to lie within theses intervals with prior probability of 95%.
In the prior model for the optode coupling parameters ζ = (s, δ , d, η) T , all the optode parameters were considered as mutually independent We modeled the amplitude parameters by uniform prior distributions π(s i ) = U(s min , 1), π(d j ) = U(d min , 1) between a selected minimum value and the ideal value 1, and for the phase parameters we used uniform prior models between ideal value of zero and a selected maximum phase shift π(δ i ) = U(0, δ max ), π(η j ) = U(0, η max ).
In the construction of the prior model π(ξ ) for the optode location vector ξ = (m, n) T , we modeled the locations of the optodes mutually independent and used parameterization ξ k = ξ 0,k + δ θ k where ξ 0,k is the angular location of the center point of optode k in the fixed parameterization ξ 0 that is to be used in the inverse problems and δ θ k ∼ U(−δ θ max , δ θ max ) models an error in the angular location.
To study the robustness of the approximation error model with respect specification of the prior model, the mean and covariance of the approximation error ε 1 (ζ ) were estimated using four different values of (s min , d min , δ max , η max ) in the prior models. Similarly, the mean and covariance of ε 2 (x, ξ ) was estimated using four different settings for δ θ max . The different values that were used are tabulated in Table 3.2. In the estimation of the statistics of ε 1 , equations (25)(26), N = 10000 random samples of the optode coupling parameters were used. In the construction of the statistics of ε 2 , equations (29-30), M = 2000 random samples of x and the optode location vectors ξ were used. The correlation structure of the covariance matrices Γ ε 1 and Γ ε 2 corresponding to prior parameters in the first row of Table 3.2 are displayed in Fig. 1.  between [a, b]. The perturbation δ θ in the location is given in degrees. The width of one source or detector fiber corresponds to 2.3 • and the separation of adjacent optodes in the equiangular case ζ 0 is 11.25 • . (1 • is equivalent to 0.44mm along the domain boundary ∂ Ω)

Coupling coefficients
Optode angular locations

Case 1: Separated and combined approximation errors
To study the effect of unknown optode coupling coefficients and locations, we considered the cases of i) pure coupling errors, ii) pure location errors and iii) combination of the coupling and location errors. The results are shown in Fig. 2. The first column in Fig. 2 shows the true absorption and scatter images (µ a top, µ s bottom). The second column in Fig. 2 shows the MAP estimate of µ a and µ s with conventional measurement error model (CEM), equation (19), when there are no optode coupling or location errors present. This estimate gives the reference estimate with the conventional model in the ideal case that there are no optode errors present.
Next, we generated a data where optode coupling losses were present but the optode locations were known exactly. The realization of ζ that was used for simulating the measurement data with coupling losses was drawn from the prior model π(ζ ) using the parameters in the first row of Table 3.2. The images on the first and second row of column three in Fig. 2 show the conventional MAP estimate (19) when the optode coupling ζ 0 is modeled (incorrectly) as ideal (measurement model y = A(x, ξ ) + e). The images on the first and second row of column 4 in Fig. 2 show the MAP estimate (24) with the approximation error model y = A(x, ξ ) + ε 1 + e where the approximation error ε 1 due to unknown optode coupling coefficients has been accounted for.
The third and fourth row in Fig 2 show results from a case where the optode coupling is exactly known (ideal coupling) but the optode locations are inaccurately known. For the simulation of the measurement data, a realization ξ of optode locations was drawn from the prior π(ξ ) with the parameters given in the first row of Table 3.2, and measurement data was simulated as y = A(x, ξ ) + e. The third column shows the conventional MAP estimate (19) using the incorrect fixed realization ξ 0 that corresponds to the equispaced locations (i.e., measurement model y = A(x, ξ 0 ) + e). The fourth column shows the MAP estimate (24) using the approximation error model y = A(x, ξ 0 ) + ε 2 (ξ ) + e where the approximation error ε 2 due to poorly known locations is taken into account.
Finally, the fifth and sixth row in Fig. 2 show a case where both, optode coupling and locations, are poorly known. The simulated measurement data was generated by drawing realizations of ξ and ζ from the their prior models using the parameters in the first row of Table 3.2 and then the data was computed as y = A(x, ξ ) + ε 1 (ζ ) + e. The third column shows the conventional MAP estimate (19) using the inexact fixed realizations (ξ 0 , ζ 0 ), where ξ 0 correspond to equispaced locations and ζ 0 to ideal coupling (i.e., measurement model y = A(x, ξ 0 ) + e). The fourth column shows the MAP estimate (24) using the approximation error model y = A(x, ξ 0 ) + ε 1 + ε 2 + e where both approximation errors are taken into account.
Evidently, as can be seen from the third column in Fig. 2, the conventional MAP estimates contain errors and distortions when optode coupling, optode locations or both are inaccurately known. Notice that the errors are especially severe in the absorption images, which is usually the parameter of main interest in medical applications. The MAP estimates with the approximation error model are basically free of these artefacts in all three situations and very similar to the reference MAP estimates in the second column that are conventional estimates in the ideal case that optode coupling and locations are exactly known. Thus, the reconstruction errors due to modeling errors in optode coupling and locations were efficiently removed using the approximation error approach.
To study the performance of the approximation error model in the ideal case when there are no actual modelling errors present in the data, we computed the reconstructions with the approximation error model from the data that was used for computing the reference reconstruction in the second column of Fig. 2. The reconstructions are shown in Fig. 3. The first column shows the estimate using the pure optode coupling error model y = A(x, ξ ) + ε 1 + e, the second column shows the reconstruction using the optode location error model y = A(x, ξ 0 ) + ε 2 (ξ ) + e and the third column show the estimate with the combined coupling and location error model y = A(x, ξ 0 ) + ε 1 + ε 2 + e. As can be seen, the approximation error model estimates are similar to the reference estimate in the second column of Fig. 2, which corresponds to correct noise model in this case. The difference in the estimates against the reference estimate is slightly larger in the cases of optode coupling error model and combined model than in the pure optode location error model. This discrepancy arises from the selection of the prior models for the optode coupling parameters; with the prior models used in the present case, the mean of optode coupling error ε 1 is non-zero and consequently the noise realization n = e + ε 1 with ε 1 = 0 has relatively low probability density with respect the actual noise model. The results indicate that the approximation error model performs also robustly in the ideal case when there are no optode coupling or location errors present in the measurement data. The CPU times for the estimates in Fig. 2 are listed the figure caption. The reference estimate in the second column corresponding to the ideal case of no modelling errors present has the shortest t CPU of all the estimates. However, the computation times of the AEM estimates in the fourth column are only moderately longer than in the ideal case. More importantly, comparing the computation times of the CEM estimates and AEM estimates in the third and fourth columns, the AEM estimate has shorter computation time than the CEM estimate in all three cases. This arises from faster convergence of the minimization when the noise model that is employed in the MAP functional is more realistic.

Case 2: Magnitude of errors in (ζ , ξ ) and sensitivity with respect the prior model
In the first test case, we used the same prior models π(ζ ) and π(ξ ) for estimation of the approximation error statistics and drawing the realizations of ζ and ξ that were used in simulating the measurement data with optode coupling or/and location errors. Basically, this case corresponds to a situation in which we know the actual prior probability distribution of the nuisance parameters. To investigate the impact of incorrect prior models of (ζ , ξ ) and how large errors the approximation error approach can tolerate in these parameters, we performed a test case where we simulated measurement data and constructed the approximation error statistics for ε 1 and ε 2 using the four different uniform prior distributions that are listed in Table 3.2. The results for the case of optode coupling errors are shown in Fig. 4 and for the case of optode location errors in Fig. 5. Both of these figures show a 4 × 4 table of MAP estimates (24) with the approximation error model such that in each image pair the left image shows µ a and the right image shows µ s . The support of the uniform prior models that were used for estimation of the approximation error statistics grows wider column wise from left to right and the support of the uniform prior that were used in the simulation of the measurement data increases from top to bottom. In the estimates indicated with small arrows the approximation error was trained with the same prior distribution that was used in the simulation of the measurement data, i.e., the image pairs on the diagonal of the 4 × 4 table correspond to the case that the actual prior distribution of ζ or ξ is known. We can see that using a prior that has a too restricted (too narrow) support compared to the actual distribution of the optode parameters leads to less efficient recovery from the modeling errors. However, training the approximation error statistics with a wider prior model than the actual distribution of the uncertainly known parameter does not seem to lead to deterioration in the recovery from the modeling error. Thus, we can say that it is safe to overestimate the uncertainty (up to a limit). In each of the image pairs, µ a is on the left and µ s on the right. The data at each of the four rows was generated using the different prior distributions given in Table 1. The AE statistics at each of the four columns was trained using the prior distributions in Table 1. The arrows denote pairs (µ a , µ s ) where the approximation error statistics was trained using the same prior that was used in simulation of the data with optode coupling errors.
When the magnitude of the errors were incremented from the largest values used in Figures 4  and 5, the recovery from the errors started to deteriorate gradually. However, the approximation error model was still able to give partial recovery even in the rather extreme case of having coupling errors randomly drawn from s i , d j ∼ U(0, 1.5) for the amplitude (zero corresponds to complete coupling failure!), and phase errors drawn from δ i , η j ∼ U(0, 3π), and location errors drawn from δ θ ∼ U(−20 • , +20 • ).

Conclusion
We have shown the feasibility of the approximation error approach for recovery from reconstruction artifacts caused by unknown source and detector optode coupling and positions. The In each of the image pairs, µ a is on the left and µ s on the right. The data at each of the four rows was generated using the prior distributions given in Table 1. The AE statistics in each of the four columns was trained using the prior distributions given in Table 1. The arrows denote pairs (µ a , µ s ) where the approximation error statistics was trained using the same prior distibution that was used in simulation of the data with misplaced optodes.
approach was tested with 2D simulations with various coupling and position errors. The results show that the approximation error approach can recover from pure and combined optode coupling and location errors. We also studied the performance using different magnitude of the errors and sensitivity with respect the specification of the prior model in the estimation of the approximation error statistics. The approach is robust for reasonably large range of errors as long as the prior model in the construction of the error model does not correspond to a too optimistic assumptions on the actual uncertainty related to the marginalized variables, in this paper, the optode locations and the coupling coefficients. Based on these results, we suggest that the approximation error approach could provide an efficient modeling protocol in practical cases where the optode coupling and locations are poorly known.