OSA's Digital Library

Optics Express

Optics Express

  • Editor: J. H. Eberly
  • Vol. 7, Iss. 13 — Dec. 18, 2000
  • pp: 468–480
« Show journal navigation

Recovery of piecewise constant coefficients in optical diffusion tomography

V. Kolehmainen, M. Vauhkonen, J. P. Kaipio, and S. R. Arridge  »View Author Affiliations


Optics Express, Vol. 7, Issue 13, pp. 468-480 (2000)
http://dx.doi.org/10.1364/OE.7.000468


View Full Text Article

Acrobat PDF (933 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

In optical diffusion tomography the reconstruction of the absorbtion and scattering coefficients is conventionally carried out in a pixel basis. The resulting number of unknowns makes the associated inverse problem severely ill-posed. We have recently proposed a new approach in which the goal is to reconstruct boundaries of piece-wise constant tissue regions as well as the diffusion and absorption coefficients within these regions. This method assumes that there is a feasible initial guess on the domain boundaries. In this paper we propose an extension to this approach in which the initial estimate for the boundary and coefficient estimation is extracted from a conventional pixel based reconstruction using standard image processing operations. In the computation of the pixel based reconstruction the output least squares problem is augmented with an approximated total variation prior. The performance of the proposed approach is evaluated using simulated frequency domain data. It is shown that since the total variation type approach favors domains with constant coefficients it is well suited for the fixing of the starting point for the actual boundary and coefficient reconstruction method.

© Optical Society of America

1 Introduction

In an optical diffusion tomography experiment, light from a near infrared laser source is guided to the surface Ω of a highly scattering object Ω by optic fibers, and the transmitted light is measured on Ω by using detecting optic fibres and light sensitive detectors. The objective is to estimate the diffusion and absorbtion coefficients (κ, µ a) within Ω based on the set of transmission data. The reconstuction of (κ, µ a) is a nonlinear ill-posed inverse problem.

In this study we consider the reconstruction of κ and µ a in the case of an object domain Ω⊂ℝ2 which is known to consist of L+1 disjoint tissue regions Ak

Ω=k=0LAk,
(1)

see Fig. 1. The regions Ak are assumed to be bounded by smooth closed boundary curves {C , =1,…, L} and have constant diffusion and absorbtion coefficients {κk , µ a,k, k=0,1,…, L} within the regions. The outer boundary of the background region A 0 is Ω.

Fig. 1. An example of a piecewise constant object domain Ω=∪ k , Ak with constant coefficients {κk , µ a,k, k=0, 1, 2, 3}. The outer boundaries of regions A 1, A 2 and A 3 are denoted by C 1, C 2 and C 3, respectively. The background region A 0 has Ω as the outer boundary.

We have recently proposed [1

1. V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data,” Inverse Problems 15, 1375–1391 (1999). [CrossRef]

, 2

2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

] an approach for the reconstruction of the shapes and locations of the region boundaries {C } and the values {κk , µak } within the regions {Ak } based on the transmission data on ∂Ω. In the proposed method, the shapes and locations of the boundaries {C } are approximated with a set of shape coefficients γ. The FEM discretisation of the forward model, which is the diffusion approximation (DA) to the radiative transfer equation (RTE), was formulated as a mapping from the set of parameters γ, κ=(κ 0,…, κL )T and µ a=(µ a,0,…, µ a,L)T to the data on the boundary Ω. The inverse problem was then stated as an optimisation problem in which we sought to minimise the residual norm between the measured and computed data.

The performance of the method was evaluated using noisy simulated frequency domain data. In [2

2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

] we assumed that the number of regions was known a priori. The method was shown to work well when the initial estimate for the parameters was relatively good in the sense that the location and size of the different tissue regions {Ak } were approximately known a priori. If the location of the regions in the initial estimate were severely erroneous the method could fail to find some regions.

In this paper we show how the shape estimation method can be initialised in the absense of a priori information about the number and locations of the tissue regions. We propose a three stage approach to the recovery of piecewise constant coefficients. As the first stage, constant coefficients (κ˜ 0, µ˜ a,0) are fitted to the whole domain Ω using nonlinear least squares. This estimate is used as the initial estimate for the second stage, in which conventional pixelwise images are reconstructed.

In the second stage the inverse problem is stated as Tikhonov regularized problem. As the regularizing penalty functional we use an approximation for the total variation norm of the solution [3

3. D. Dobson and F. Santosa, “An image enhancement technique for electrical impedance tomography,” Inverse Problems 10, 317–334 (1994). [CrossRef]

, 4

4. J. P. Kaipio, V. Kolehmainen, E. Somersalo, and M. Vauhkonen, “Statistical inversion methods in electrical impedance tomography,” Inverse Problems (2000), in Press.

, 5

5. K. D. Paulsen and H. Jiang, “Enhanced frequency-domain optical image reconstruction in tissues through total variation minimization,” Appl. Opt. 35, 3447–3458 (1996). [CrossRef] [PubMed]

]. The regularizing functional is also called a prior due to its identification with the prior distribution in the statistical interpretation of inverse problems [4

4. J. P. Kaipio, V. Kolehmainen, E. Somersalo, and M. Vauhkonen, “Statistical inversion methods in electrical impedance tomography,” Inverse Problems (2000), in Press.

]. From the statistical point of view, the total variation norm is a feasible choice as a regularizing penalty functional when the images are assumed to consist of a few different levels of {κk , µ a,k} with relatively short boundary lines {C }. This is because the total variation penalty functional tends to favor such coefficient distributions which have domains with short boundaries and almost constant values within each region. Once approximate convergence of the second stage is achieved, the initial configuration of the boundaries {C } can be found from the pixelwise reconstructions using standard image processing techniques such as edge detection. An example is given in Section 4.

In the third stage, the final images are reconstructed using the shape estimation method that was proposed in [1

1. V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data,” Inverse Problems 15, 1375–1391 (1999). [CrossRef]

, 2

2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

]. To improve the convergence speed of the the earlier approach, we now also use line search. The performance of this three stage inversion scheme is evaluated by simulations in which the optical parameters (κ, µ a) are relevant for medical applications of optical tomography.

The rest of this paper is organised as follows. In Section 2 we discuss briefly the parametrisation γ of the region boundaries {C } and the diffusion approximation to the RTE. In Section 3 a short stage-by-stage description of the new three stage inversion approach is given and in Section 4 the performance of this approach is evaluated using noisy simulated data. In Section 5 we give conclusions and some suggestions for future work.

2 Forward Model

2.1 Representation of the boundaries {Cℓ}

In this study the domains {Ak } ∊ Ω are assumed to be simply connected and disjoint (i.e. CiCj =∅, ij), see Fig. 1. We also assume that the outer boundary of the region A 0, that is, Ω is known and has the property Ω∩C =∅ ∀ℓ. If the outer boundaries {C , =1, 2,…,L} of the regions {A } are sufficiently smooth, they can be approximated in the form

{C(S)}=(x(s)y(s))=n=1Nθ(γnxθnx(s)γnyθny(s)),=1,...,L
(2)

where θn are periodic and differentiable basis functions with period 1. In this study we express both coordinates of the boundary curves as Fourier series with respect to the curve parameter s, that is, we use basis functions of the form

θ1α(s)=1
θnα(s)=sin(2πn2(s+ϕs)),n=2,4,6,
θnα(s)=cos(2π(n12)(s+ϕs)),n=3,5,7,
(3)

where ϕs is the phase of the curve parameter and α denotes either x or y. Furthermore, let γ denote the vector of the shape coefficients, that is,

γ=(γ1x1,,γNθx1,γ1y1,,γNθy1,,,γ1xL,,γNθxL,γ1yL,,γNθyL)T.
(4)

For more details on the properties of the boundary representation (2) and (3), see [1

1. V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data,” Inverse Problems 15, 1375–1391 (1999). [CrossRef]

, 2

2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

].

2.2 Diffusion approximation to the RTE

Consider the following experiment: S optic fibers are placed on the source positions εjΩ on the boundary of the highly scattering (µ′ sµ a) object Ω and M optic fibers are placed in the detector positions ζ iΩ. Light from an intensity modulated near infrared laser source is guided to the object via one of the source fibers at εj , and the amount of transmitted light is measured on all the detector locations ζ i ,i=1,…,M. Then this process is repeated for all S source locations.

A well established model for the above experiment is the diffusion approximation (DA) of the radiative transfer equation [6

6. M. C. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, New York, 1967).

, 7

7. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). [CrossRef]

]. In the case of an intensity modulated light source, the diffusion approximation can be expressed as

·κΦj(ω)+μaΦj(ω)+iωcΦj(ω)=q0j,rΩ
Φj(ω)+2κϑΦj(ω)ν=gsjrΩ,
(5)

where Φ j is the photon density µ a is the absorption coefficient [mm-1], κ=(3(µ a+µ′ s))-1 [mm] is the diffusion coefficient, µ′ s is the reduced scattering coefficient [mm-1], c is the speed of light in the medium, q0 j is the distribution of internal light sources, g S j is the distribution of boundary sources, v is the boundary normal direction, ν is a coefficient due to the mismatch of the refractive indices in Ω and the surrounding medium and the subindex j denotes that the object is illuminated from the source location εj .

Two widely used source models within the DA framework are the collimated source (CS) model and the diffuse boundary source (DS) model [8

8. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element model for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22, 1779–1792 (1995). [CrossRef] [PubMed]

]. In the case of the CS model g Sj becomes zero, and in the case of the DS model qo j becomes zero.

The relation between the photon density and the measured exitance at the measurement location ζ iΩ can be expressed as

Γi,j(ω)=κ(ζi)ν·Φj(ω).
(6)

For further details on the diffusion approximation and source models, see e.g. [6

6. M. C. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, New York, 1967).

, 9

9. J. P. Kaltenbach and M. Kaschke, “Frequency- and Time-domain Modelling of Light Transport in Random Media,” in Medical Optical Tomography: Functional Imaging and Monitoring, G. Muller, B. Chance, R. Alfano, S. Arridge, J. Beuthan, E. Gratton, M. Kaschke, B. Masters, S. Svanberg, and P. van der Zee, eds., (SPIE, Bellingham, WA, 1993), pp. 65–86.

, 10

10. R. A. J. Groenhuis, H. A. Ferwerda, and J. J. T. Bosch, “Scattering and Absorption of Turbid Materials Determined from Reflection Measurements. Part 1: Theory,” Appl. Opt. 22, 2456–2462 (1983). [CrossRef] [PubMed]

, 7

7. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). [CrossRef]

, 8

8. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element model for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22, 1779–1792 (1995). [CrossRef] [PubMed]

, 11

11. R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

].

In this study, the discretisation of the forward model (5–6) is based on the finite element method. The details of the discretisation are given in the references [1

1. V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data,” Inverse Problems 15, 1375–1391 (1999). [CrossRef]

, 2

2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

, 12

12. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A Finite Element Approach for Modeling Photon Transport in Tissue,” Med. Phys. 20, 299–309 (1993). [CrossRef] [PubMed]

, 8

8. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element model for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22, 1779–1792 (1995). [CrossRef] [PubMed]

]. In the following discussion, the vector of image parameters will be denoted by f. The FEM based map P : f↦>z which maps the image parameter vector f to the data vector z on Ω, will be denoted by

z=𝓟(f).
(7)

In this paper the data vector is of the form

z=(re(Γ1(ω)),,re(Γs(ω)),im(Γ1(ω)),,im(Γs(ω))T,
(8)

where Γ j (ω)=(Γ1,j(ω),…, ΓM,j(ω)) is the vector of the complex measurements for the j:th source [13

13. H. Jiang, K. D. Paulsen, U. L. Osterberg, B. W. Pogue, and M. S. Patterson, “Optical Image Reconstruction Using Frequency-domain Data: Simulations and Experiments,” J. Opt. Soc. Am. A 13, 253–266 (1995). [CrossRef]

]. This splitting of the data vector means that the discrete mapping (7) ensures that the adjoint operator gives real updates to the solution.

3 Inverse Problem

3.1 Reconstruction of the best fitting constants µ˜a,0 and κ˜0

f̂=argminfzmeas𝓟(f)zmeas22,
(9)

where the parameter vector is f=(κ˜ 0, µ˜ a,0)T ∊ ℝ×ℝ. The inverse problem (9) is usually stable (well-posed) due to the small number of unknowns and thus the problem can be solved without additional regularization.

A well-established method for the solution of stable nonlinear optimization problems is the Gauss-Newton method, which can be written formally as

f(l+1)=f(l)+(J˜(l)TJ˜(l))1J˜lTz˜(l),
(10)

where J̃(l)=diag(zmeas1)J(l) is the (scaled) Jacobian matrix of P(f (l)) and (l)=zmeas1(z meas-P(f (l)).

The motivation behind this stage is to find a good initial estimate for the second step which is the computation of conventional pixelwise reconstructions with the total variation penalty functional.

3.2 Reconstruction in a local pixel basis with the total variation regularization

κ=m=1Npκmχm,μa=m=1Npμa,mχm,
(11)

where xm is the characteristic function of pixel Ω m . Within the approximation (11), the functions (κ, µ a) are identified with the vectors μa=(μa,1,,μa,Np)TNp, κ=(κ1,,κNp)TNp and the parameter vector for the inverse problem becomes

f=(κμa)Np×Np.
(12)

At this stage, the ill-posed inverse problem is regularised by using the generalised Tikhonov regularisation. The regularised least squares problem becomes

f̂=argminf{zmeas𝓟(f)zmeas22+Ψ(f)},
(13)

where ψ(f)>0 is regularizing penalty functional. In this study, ψ(f) is of the form

Ψ(f)=βκTV(κ)+βμaTV(μa),
(14)

where TV(·) refers to the total variation norm [3

3. D. Dobson and F. Santosa, “An image enhancement technique for electrical impedance tomography,” Inverse Problems 10, 317–334 (1994). [CrossRef]

, 4

4. J. P. Kaipio, V. Kolehmainen, E. Somersalo, and M. Vauhkonen, “Statistical inversion methods in electrical impedance tomography,” Inverse Problems (2000), in Press.

]. Here βκ,βμa are constant hyper-parameters whose choice is disscussed in Section 4. With the piecewise constant basis for the coefficients (11), the TV norm can be expressed for both κ,µ a as

TV(α)=j=1JdjΔjTα,
(15)

where α is either κ or µ a, dj is the length of the jth edge between the adjacent pixels Ωi1j and Ωi2j and ΔjNp is the vector

Δj=(0,,1,(i1j)0,,0,-1,(i2j)0,,0)T

The choice of the regularizing penalty functional (14) is based on the assumed properties of the true images (κ,µ a): It is well known that the total variation norm is a feasible prior for piecewise constant images which consist of a few constant levels {κk ,µ a,k} with relatively short well defined boundary lines {C } [3

3. D. Dobson and F. Santosa, “An image enhancement technique for electrical impedance tomography,” Inverse Problems 10, 317–334 (1994). [CrossRef]

, 4

4. J. P. Kaipio, V. Kolehmainen, E. Somersalo, and M. Vauhkonen, “Statistical inversion methods in electrical impedance tomography,” Inverse Problems (2000), in Press.

]. While the TV prior has found much attention in impedance tomography, it seems to have been used in optical tomography previously only in [5

5. K. D. Paulsen and H. Jiang, “Enhanced frequency-domain optical image reconstruction in tissues through total variation minimization,” Appl. Opt. 35, 3447–3458 (1996). [CrossRef] [PubMed]

], where positive correlation of the prior distributions of (κ, µ a) was implicitly assumed, and only the special case of identical contrast in absorption and scattering was examined.

The problem (13) is solved using the Gauss-Newton method. Here, we face the difficulty that the total variation penalty functional is non-differentiable. To overcome this problem, we approximate the L 1 norm in (15) with the smooth function

thτ(t)=1τlog(cosh(τt)),
(16)

W(f(l))=ςlm=12Np1fm(l),
(17)

where {ϛl} is a sequence of decreasing positive coefficients (i.e. ϛl→0 as l→∞), fm(l) are the elements of the vector f (l) and l denotes the iteration index [16

16. A. Fiacco and G. McCormick, Nonlinear Programming (SIAM, 1990). [CrossRef]

]. With these modifications, the Gauss-Newton algorithm becomes

f(l+1)=f(l)+(J˜(l)TJ˜(l)+HTVτ(l)+HW(l))1(J˜(l)Tz˜lgTVτ(l)gW(l)),
(18)

where gTVτ(l) and HTVτ(l) are gradient and Hessian of the approximated total variation functional (gW(l) and HW(l) respectively for the functional W(f (l)). The details for the computation of the map P(f (l)) and the Jacobian matrix J(l) in equations (10) and (18) can be found in [7

7. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). [CrossRef]

]. Details for the computation of gTVτ(l) and HTVτ(l) can be found in [17

17. S. Järvenpää, “A finite element model for the inverse conductivity problem,”, Phil. Lic. thesis, University of Helsinki, Finland (1996).

, 18

18. V. Kolehmainen, E. Somersalo, P. J. Vauhkonen, M. Vauhkonen, and J. P. Kaipio, “A Bayesian approach and total variation priors in 3D electrical impedance tomography,” In Proc 20th Ann Int Conf IEEE Eng Med Biol Soc, pp. 1028–1031 (Hong Kong, China, 1998).

]. It should be noted that the effect of the functional (17) vanishes as l→∞ and the algorithm (18) converges to the inequality constrained (i.e. f>0) minimum of the functional (13) [16

16. A. Fiacco and G. McCormick, Nonlinear Programming (SIAM, 1990). [CrossRef]

].

The initial estimate for the final stage, which is the actual shape estimation method, can be extracted from the local pixel base reconstructions by using a number of standard image processing tools. An example of feasible image processing steps is given in Section 4.

3.3 Shape estimation

At the final stage of our approach, we utilise a modified version of the shape estimation method which was proposed in [2

2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

]. The goal in the method is to find the representation of the boundary shapes γ, and the values κ=(κ 0,…,κL )T, µ a=(µ a,0,µ a,L)T which give the best match to the data in the output LS sense. Thus, the problem can be stated as

f̂=argminfzmeas𝓟(f)zmeas22,
(19)

where the parameters are

f=(γκμa)2LNθxL+1xL+1.
(20)

If the number of parameters is not unreasonably high, the inverse problem is not illposed and regularization is not needed. However, to secure numerical stability during the iteration we use a Levenberg-Marquardt type method, which can be formally written as

f(l+1)=f(l)+s(l)d(l)
(21)
d(l)=(J˜(l)TJ˜(l)+λI)1J˜(l)Tz˜(l),
(22)

Ξ(l)(s)zmeas𝓟(f(l)+sd(l))zmeas22
(23)

In [2

2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

], the optimisation problem (19) was solved in rescaled parameter space in order to avoid numerical instabilities which are due to the different magnitudes of the image parameters γ, κ, µ a. In this study, we use a different parameter scaling scheme: Instead of scaling the parameters by the initial estimate f (0), the parameters are rescaled by γ˜ =γ/max(|γ (0)|), κ˜=κ/max(|κ (0)|) and µ˜ a=µ a/max(|µa(0)), and the optimisation is performed in the scaled parameter space. This choice yields a numerically more stable problem and speeds up the convergence noticeably.

The details for the computation of the map P(f) and the Jacobian J in (22) can be found in [1

1. V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data,” Inverse Problems 15, 1375–1391 (1999). [CrossRef]

, 2

2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

]. We note here only that the phase ϕs of the curve parameter s in equation (3) has to be fixed when applying equations (21)(22). If the phase (ϕs is not fixed we obtain a one-dimensional null space N with respect to each curve C in the block Jγ of the Jacobian J.

4 Results

The synthetic measurement system that we use in this paper consists of 16 sources and 16 detectors placed in equiangular positions with the order ε 1, ζ1, ε2, ζ2,…, ε 16, ζ16 on the boundary Ω of a circular domain Ω which has the radius of 25mm. The optical coefficients within Ω were chosen in the range of interest in medical imaging [7

7. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). [CrossRef]

]. For the absorption coefficient µ a we used values in the range 0.025mm-1→0.05mm-1 and for the reduced scattering coefficient µ′ s we used values in the range 2mm-1→ 4mm-1. The diffusion coefficient is computed as κ=(3(µ a+µ′ s))-1. The tissue regions {Ak } were bounded by trigonometric boundaries.

The true target images are shown in the top row of Fig. 2. The regions A 1 and A 2 have significant contrast in both the diffusion and absorption coefficients with respect to the background A 0. The region A 3, which is located on the lower right side of Ω, has significant contrast in the absorption coefficient and only a weak contrast in the diffusion coefficient with respect to the surrounding region A 0.

In the generation of the synthetic data with the FEM model we discretised the object domain Ω to P=8600 disjoint triangular elements with D=4429 vertex nodes. The simulated data was computed using a single frequency ω=300MHz and a multiplicative random noise with standard deviation of 1% of the measured signal was added to the simulated data [14

14. S. R. Arridge, M. Hiraoka, and M. Schweiger, “Statistical Basis for the Determination of Optical Pathlength in Tissue,” Phys. Med. Biol. 40, 1539–1558 (1995). [CrossRef] [PubMed]

]. In the image reconstruction process we used a FEM model in which Ω was divided to P=5776 triangular elements with D=3017 vertex nodes.

In the recovery of the best matching constants we used the initial guess (κ˜ 0, µ˜ a,0)=(0.1104,0.0200). The convergence of the iteration was verified by monitoring the output residual norm and the behaviour of the image parameters. The iteration (10) reached full convergence in 12 steps, having minimum at (κ˜ 0, µ˜ a,0)=(0.1642,0.0301).

The second stage, which is the reconstruction into the local pixel basis, was started from the constant values (κ 0, µ a,0)=(0.1642,0.0301). The domain was discretised into 13 concentric layers of generic quadrilateral pixels, the total number of pixels being Np =496. Each of these pixels is a union of several triangular FE elements. This disretisation is essentially identical to the “Joshua tree” in [19

19. M. Cheney, D. Isaacson, J. Newell, S. Simske, and J. Goble, “NOSER: An algorithm for solving the inverse conductivity problem,” Int J Imaging Systems and Technology 2, 66–75 (1990). [CrossRef]

].

In general, there is not any automatic methods for the choice of two regularization parameters βκ and βμa in equations (13)(14), and thus these have to be chosen manually. In the proposed approach the parameters βκ and βμa are chosen large enough to produce images with flat details and adequately well-defined boundaries. This is adequate since the purpose of this stage is only to produce images which serve as a basis for the initialisation of the boundaries. The initial values of the coefficients are not as crucial as feasible approximate boundaries. We note that the feasible values of βκ and βμa depend on imaging parameters such as geometry of the domain Ω, data types, noise levels of the instrumentation and the number of pixels in the mesh. For example, for the given imaging parameters we have used fixed values βκ =50 and βμa=150.

The second stage includes the fixed hyperparameter τ in equation (16) and the sequence {ϛl} of the positivity constraint parameters. Consider first the choice of the fixed parameter τ. The function hτ (t) has the property that hτ (t)→|t| as τ increases. However, as τ increases the second derivatives of hτ (t), which are needed in the Hessian matrix ΗTVτ(l), tend to become very large, leading to numerical instabilities in the iteration (18). Thus, a tradeoff between the accuracy of the approximation (16) and the numerical stability of the iteration (18) has to be made. We have used the value τ=2.5. The choice of the decreasing sequence {ϛl} in iteration (18) is not crucial as long as the positivity constraint remains intact since the asymptotic estimate is not a affected by the positivity constraint [16

16. A. Fiacco and G. McCormick, Nonlinear Programming (SIAM, 1990). [CrossRef]

]. For the given geometry, we have used the sequence {ϛl+1=0.5 ϛ,l l≥1} with the initialisation ϛ0=10-4. The essential convergence of the iteration (18) was reached in 5 steps. The resulting images are shown in the second row of Fig. 2.

Fig. 2. A test case with three boundaries {C }. Region A 1 is located up, A 2 is located on the lower left side of Ω and A 3 on the lower right side. Left column shows images of κ([mm]) and right column shows images of µ a ([mm-1]). Rows from top to bottom: True distributions, pixelwise reconstructions with the approximated TV prior, initial estimates for the shape estimation and the bottom row shows the final estimates with the shape estimation method. The animation (349kB) shows the progress of the inversion approach step by step.

The initial estimate for the shape estimation method was extracted from the pixelwise reconstructions by standard image processing techniques. To make the image processing possible by standard tools, the pixelwise reconstructions were mapped from the “Joshua tree” discretisation to a regular Nr ×Nr grid in the domain G=[-25mm, 25mm]×[-25mm, 25mm] ⊂ℝ2. The number Nr was chosen such that the resolution of the regular grid was approximately equal to the resolution of the “Joshua tree” mesh. The left hand image in Fig. 3 shows the reconstructed µ a in the rectangular mesh. For the grid points which lay outside Ω, we attached values from a pixel on the outermost layer of the Joshua tree mesh. Next, the mapped images of (κ, µ a) were filtered with a band pass filter and then the edges were detected by looking for zero crossings from the filtered images [20

20. J. S. Lim, Two-dimensional signal and image processing (Prentice Hall, Englewood Cliffs, NJ, 1990).

, 21

21. Matlab: Image Processing toolbox user’s guide, 2 ed., The MathWorks Inc., 24 Prime Park Way, Natick, MA 01760-1500, USA.

]. The edge detection program returned binary images, having ones at the points where edges were detected, and zeros elsewhere. The right image in Fig. 3 shows the edges that were found from the µ a reconstruction in the left image of Fig. 3. Based on the reconstructed images and the detected edges, the distinct regions {Ak , k=1, 2, 3} in addition to the background A 0 were clearly detectable. These regions were marked and boundary curves were fitted to the detected boundary pixels by the LS method, see the animation. These boundary curves and the average values of the coefficients within these regions were then used as the initial guess for the final stage. The initial guess for the final stage is shown in third row of Fig. 2. The initial values of the coefficients {κk , µ a,k, k=0,…,3} can be found in Table 1.

Fig. 3. Left: The reconstructed µ a in a regular grid. Right: The edges that were detected from the µ a image.

Table 1. Values of {κk , µa,k} in the initial and in the final estimates of the shape estimation in Fig. 2. The initial estimate is shown in the 3rd row in Fig. 2, and the final estimate in the bottom row, respectively. The dimension of κk is [mm] and the dimension of µa,k is [mm-1].

table-icon
View This Table

does not depend on λ. The feasible range for λ depends on the norm of the Jacobian and thus also on the imaging parameters. For the given imaging parameters, we have used typically the value λ=1. Starting from the initial estimate in the third row of Fig. 2, the iteration (21)–(22) converged in 26 iterations. The final estimates for (κ, µ a) in the case of this example are shown in the bottom row of Fig. 2. The error statistics of the initial and reconstructed values of {κ,k , µ a,k} in the third stage are given in Table 1.

As can be seen from Fig. 2, the shapes of the boundaries of the regions A 1A 3 are found with relatively good accuracy. Also the diffusion and absorption coefficients {κk , µ a,k} are found with good accuracy resulting in final parameter estimation error norms in the range ~0–4%, see Table 1.

5 Conclusions

Based on the shape estimation method, which was proposed in [2

2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

, 1

1. V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data,” Inverse Problems 15, 1375–1391 (1999). [CrossRef]

], we proposed a new three stage scheme for the reconstruction of piecewise constant domains in optical diffusion tomography. In the proposed approach, the shape estimation method is initialised by a local pixel basis reconstruction that is carried out by generalized Tikhonov regularization with a total variation type prior. Standard image processing methods are then used to obtain initial approximations for the boundaries and the coefficients.

The performance of the new approach was evaluated using noisy simulated frequency domain data. It was previously shown that the actual boundary recovery approach is able to yield good estimates for both the boundaries and the coefficients once the initial guess for the iteration was good enough. It was now shown that the proposed approach for the initialization yield the initial configuration with adequate accuracy. Although the performance of the method was evaluated with frequency domain data, the method can be extended to time resolved data types in a straightforward manner. The performance of the method using different data types is a topic of further research.

Although the method assumes that the domains are piecewise constant, it may prove to be useful also in such cases in which the the diffusion and absorption may have small spatial variations within the tissue regions {Ak } and the largest contrast exist over the boundaries {C }. The latter approach has also been discussed in [23

23. M. E. Kilmer, E. L. Miller, D. A. Boas, D. H. Brooks, C. A. DiMarzio, and R. J. Gaudette, “Direct object localization and characterization from diffuse photon density data,” Proceedings of the SPIE Photonics West Meeting, Jan. 1999.

] wherein an alternative shape estimation method has been investigated. In that case, B-splines with a finite number of control points were used rather than trigonometric functions, and the space of control points was searched with a local greedy strategy and explicit perturbation of a linear forward model.

The method in the form that it was presented here is applicable to 2D situations. In this study we assumed a circular domain, although, in principle, it can also be adopted to other geometries, such as the 2D reflection geometry investigated in [23

23. M. E. Kilmer, E. L. Miller, D. A. Boas, D. H. Brooks, C. A. DiMarzio, and R. J. Gaudette, “Direct object localization and characterization from diffuse photon density data,” Proceedings of the SPIE Photonics West Meeting, Jan. 1999.

], by minor modifications to the FEM based forward operator P(f). In that case it may be necessary to include shape regularisation methods as required in [23

23. M. E. Kilmer, E. L. Miller, D. A. Boas, D. H. Brooks, C. A. DiMarzio, and R. J. Gaudette, “Direct object localization and characterization from diffuse photon density data,” Proceedings of the SPIE Photonics West Meeting, Jan. 1999.

] although such methods were not needed here.

However, it is well known that optical diffusion tomography should be considered a 3D problem [24

24. M. Schweiger and S. R. Arridge, “Comparison of 2D and 3D reconstruction algorithms in Optical Tomography,” Appl. Opt. 37, 7419–7428 (1998). [CrossRef]

]. The validity of a 3D FEM model for optical diffusion tomography has been shown in simulations in [24

24. M. Schweiger and S. R. Arridge, “Comparison of 2D and 3D reconstruction algorithms in Optical Tomography,” Appl. Opt. 37, 7419–7428 (1998). [CrossRef]

] and using measured data in [25

25. S. R. Arridge, J. C. Hebden, M. Schweiger, F. E. W. Schmidt, M. E. Fry, E. M. C. Hillman, H. Dehghani, and D. T. Delby, “A method for three-dimensional time-resolved optical tomography,” Int. J. Imaging Syst. Technol. 11, 2–11 (2000). [CrossRef]

].

Acknowledgements

This work was supported by the Saastamoinen Foundation, Academy of Finland and the Vilho, Yrjö and Kalle Väisälä fund. S. Arridge acknowledges the support of the Wellcome Trust.

References and links

1.

V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data,” Inverse Problems 15, 1375–1391 (1999). [CrossRef]

2.

V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, “Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography,” Phys Med Biol (2000), in Press.

3.

D. Dobson and F. Santosa, “An image enhancement technique for electrical impedance tomography,” Inverse Problems 10, 317–334 (1994). [CrossRef]

4.

J. P. Kaipio, V. Kolehmainen, E. Somersalo, and M. Vauhkonen, “Statistical inversion methods in electrical impedance tomography,” Inverse Problems (2000), in Press.

5.

K. D. Paulsen and H. Jiang, “Enhanced frequency-domain optical image reconstruction in tissues through total variation minimization,” Appl. Opt. 35, 3447–3458 (1996). [CrossRef] [PubMed]

6.

M. C. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, New York, 1967).

7.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). [CrossRef]

8.

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element model for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22, 1779–1792 (1995). [CrossRef] [PubMed]

9.

J. P. Kaltenbach and M. Kaschke, “Frequency- and Time-domain Modelling of Light Transport in Random Media,” in Medical Optical Tomography: Functional Imaging and Monitoring, G. Muller, B. Chance, R. Alfano, S. Arridge, J. Beuthan, E. Gratton, M. Kaschke, B. Masters, S. Svanberg, and P. van der Zee, eds., (SPIE, Bellingham, WA, 1993), pp. 65–86.

10.

R. A. J. Groenhuis, H. A. Ferwerda, and J. J. T. Bosch, “Scattering and Absorption of Turbid Materials Determined from Reflection Measurements. Part 1: Theory,” Appl. Opt. 22, 2456–2462 (1983). [CrossRef] [PubMed]

11.

R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

12.

S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A Finite Element Approach for Modeling Photon Transport in Tissue,” Med. Phys. 20, 299–309 (1993). [CrossRef] [PubMed]

13.

H. Jiang, K. D. Paulsen, U. L. Osterberg, B. W. Pogue, and M. S. Patterson, “Optical Image Reconstruction Using Frequency-domain Data: Simulations and Experiments,” J. Opt. Soc. Am. A 13, 253–266 (1995). [CrossRef]

14.

S. R. Arridge, M. Hiraoka, and M. Schweiger, “Statistical Basis for the Determination of Optical Pathlength in Tissue,” Phys. Med. Biol. 40, 1539–1558 (1995). [CrossRef] [PubMed]

15.

M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999). [CrossRef] [PubMed]

16.

A. Fiacco and G. McCormick, Nonlinear Programming (SIAM, 1990). [CrossRef]

17.

S. Järvenpää, “A finite element model for the inverse conductivity problem,”, Phil. Lic. thesis, University of Helsinki, Finland (1996).

18.

V. Kolehmainen, E. Somersalo, P. J. Vauhkonen, M. Vauhkonen, and J. P. Kaipio, “A Bayesian approach and total variation priors in 3D electrical impedance tomography,” In Proc 20th Ann Int Conf IEEE Eng Med Biol Soc, pp. 1028–1031 (Hong Kong, China, 1998).

19.

M. Cheney, D. Isaacson, J. Newell, S. Simske, and J. Goble, “NOSER: An algorithm for solving the inverse conductivity problem,” Int J Imaging Systems and Technology 2, 66–75 (1990). [CrossRef]

20.

J. S. Lim, Two-dimensional signal and image processing (Prentice Hall, Englewood Cliffs, NJ, 1990).

21.

Matlab: Image Processing toolbox user’s guide, 2 ed., The MathWorks Inc., 24 Prime Park Way, Natick, MA 01760-1500, USA.

22.

D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” J. Soc. Indust. Appl. Math. 11, 431–441 (1963). [CrossRef]

23.

M. E. Kilmer, E. L. Miller, D. A. Boas, D. H. Brooks, C. A. DiMarzio, and R. J. Gaudette, “Direct object localization and characterization from diffuse photon density data,” Proceedings of the SPIE Photonics West Meeting, Jan. 1999.

24.

M. Schweiger and S. R. Arridge, “Comparison of 2D and 3D reconstruction algorithms in Optical Tomography,” Appl. Opt. 37, 7419–7428 (1998). [CrossRef]

25.

S. R. Arridge, J. C. Hebden, M. Schweiger, F. E. W. Schmidt, M. E. Fry, E. M. C. Hillman, H. Dehghani, and D. T. Delby, “A method for three-dimensional time-resolved optical tomography,” Int. J. Imaging Syst. Technol. 11, 2–11 (2000). [CrossRef]

26.

C. Brechbühler, G. Gerig, and O. Kübler, “Parametrization of closed surfaces for 3-D shape description,” Computer Vision and Image Understanding 61, 154–170 (1995). [CrossRef]

27.

A. Kelemen, G. Szekely, and G. Gerig, “Three-dimensional model-based segmentation,” IEEE Trans Med Imaging 18, 828–839 (1995). [CrossRef]

OCIS Codes
(100.3010) Image processing : Image reconstruction techniques
(100.3190) Image processing : Inverse problems

ToC Category:
Focus Issue: Diffuse optical tomography

History
Original Manuscript: October 16, 2000
Published: December 18, 2000

Citation
V. Kolehmainen, M. Vauhkonen, Jari Kaipio, and Simon Arridge, "Recovery of piecewise constant coefficients in optical diffusion tomography," Opt. Express 7, 468-480 (2000)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-7-13-468


Sort:  Journal  |  Reset  

References

  1. V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio, "Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data," Inverse Problems 15, 1375-1391 (1999). [CrossRef]
  2. V. Kolehmainen, S. R. Arridge, M. Vauhkonen, and J. P. Kaipio, "Simultaneous reconstruction of internal tissue region boundaries and coefficients in optical diffusion tomography," Phys Med Biol (2000), in Press.
  3. D. Dobson and F. Santosa, "An image enhancement technique for electrical impedance tomography," Inverse Problems 10, 317-334 (1994). [CrossRef]
  4. J. P. Kaipio, V. Kolehmainen, E. Somersalo, and M. Vauhkonen, "Statistical inversion methods in electrical impedance tomography," Inverse Problems (2000), in Press.
  5. K. D. Paulsen and H. Jiang, "Enhanced frequency-domain optical image reconstruction in tissues through total variation minimization," Appl. Opt. 35, 3447-3458 (1996). [CrossRef] [PubMed]
  6. M. C. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, New York, 1967).
  7. S. R. Arridge, "Optical tomography in medical imaging," Inverse Problems 15, R41-R93 (1999). [CrossRef]
  8. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, "The finite element model for the propagation of light in scattering media: Boundary and source conditions," Med. Phys. 22, 1779-1792 (1995). [CrossRef] [PubMed]
  9. J. P. Kaltenbach and M. Kaschke, "Frequency- and Time-domain Modelling of Light Transport in Random Media," in Medical Optical Tomography: Functional Imaging and Monitoring, G. Muller, B. Chance, R. Alfano, S. Arridge, J. Beuthan, E. Gratton, M. Kaschke, B. Masters, S. Svanberg, and P. van der Zee, eds., (SPIE, Bellingham, WA, 1993), pp. 65-86.
  10. R. A. J. Groenhuis, H. A. Ferwerda, and J. J. T. Bosch, "Scattering and Absorption of Turbid Materials Determined from Reflection Measurements. Part 1: Theory," Appl. Opt. 22, 2456-2462 (1983). [CrossRef] [PubMed]
  11. R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, "Boundary conditions for the diffusion equation in radiative transfer," J. Opt. Soc. Am. A 11, 2727-2741 (1994). [CrossRef]
  12. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, "A Finite Element Approach for Modeling Photon Transport in Tissue," Med. Phys. 20, 299-309 (1993). [CrossRef] [PubMed]
  13. H. Jiang, K. D. Paulsen, U. L. Osterberg, B. W. Pogue, and M. S. Patterson, "Optical Image Reconstruction Using Frequency-domain Data: Simulations and Experiments," J. Opt. Soc. Am. A 13, 253-266 (1995). [CrossRef]
  14. S. R. Arridge, M. Hiraoka, and M. Schweiger, "Statistical Basis for the Determination of Optical Pathlength in Tissue," Phys. Med. Biol. 40, 1539-1558 (1995). [CrossRef] [PubMed]
  15. M. Schweiger and S. R. Arridge, "Application of temporal filters to time resolved data in optical tomography," Phys. Med. Biol. 44, 1699-1717 (1999). [CrossRef] [PubMed]
  16. A. Fiacco and G. McCormick, Nonlinear Programming (SIAM, 1990). [CrossRef]
  17. S. Järvenpää, "A finite element model for the inverse conductivity problem,", Phil. Lic. thesis, University of Helsinki, Finland (1996).
  18. V. Kolehmainen, E. Somersalo, P. J. Vauhkonen, M. Vauhkonen, and J. P. Kaipio, "A Bayesian approach and total variation priors in 3D electrical impedance tomography," In Proc 20th Ann Int Conf IEEE Eng Med Biol Soc, pp. 1028-1031 (Hong Kong, China, 1998).
  19. M. Cheney, D. Isaacson, J. Newell, S. Simske, and J. Goble, "NOSER: An algorithm for solving the inverse conductivity problem," Int J Imaging Systems and Technology 2, 66-75 (1990). [CrossRef]
  20. J. S. Lim, Two-dimensional signal and image processing (Prentice Hall, Englewood Cliffs, NJ, 1990).
  21. Matlab: Image Processing toolbox user's guide, 2 ed., The MathWorks Inc., 24 Prime Park Way, Natick, MA 01760-1500, USA.
  22. D. W. Marquardt, "An algorithm for least-squares estimation of nonlinear parameters," J. Soc. Indust. Appl. Math. 11, 431-441 (1963). [CrossRef]
  23. M. E. Kilmer, E. L. Miller, D. A. Boas, D. H. Brooks, C. A. DiMarzio and R. J. Gaudette, "Direct object localization and characterization from diffuse photon density data," Proceedings of the SPIE Photonics West Meeting, Jan. 1999.
  24. M. Schweiger and S. R. Arridge, "Comparison of 2D and 3D reconstruction algorithms in Optical Tomography," Appl. Opt. 37, 7419-7428 (1998). [CrossRef]
  25. S. R. Arridge, J. C. Hebden, M. Schweiger, F. E. W. Schmidt, M. E. Fry, E. M. C. Hillman, H. Dehghani, and D. T. Delby, "A method for three-dimensional time-resolved optical tomography," Int. J. Imaging Syst. Technol. 11, 2-11 (2000). [CrossRef]
  26. C. Brechbühler, G. Gerig, and O. Kübler, "Parametrization of closed surfaces for 3-D shape description," Computer Vision and Image Understanding 61, 154-170 (1995). [CrossRef]
  27. A. Kelemen, G. Szekely, and G. Gerig, "Three-dimensional model-based segmentation," IEEE Trans Med Imaging 18, 828-839 (1995). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1. Fig. 2. Fig. 3.
 

Multimedia

Multimedia FilesRecommended Software
» Media 1: GIF (349 KB)      QuickTime

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited