Tomography with Synchrotron Radiation Alvaro R. De Pierro, - - PowerPoint PPT Presentation

Tomography with Synchrotron Radiation

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil IPUC, Florian´

• polis, April 2014

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Work by

◮ Eduardo Xavier Miqueles, Brazilian Synchrotron Light

Laboratory (LNLS)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Work by

◮ Eduardo Xavier Miqueles, Brazilian Synchrotron Light

Laboratory (LNLS)

◮ Elias Helou, ICMC-USP, Brazil

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Work by

◮ Eduardo Xavier Miqueles, Brazilian Synchrotron Light

Laboratory (LNLS)

◮ Elias Helou, ICMC-USP, Brazil ◮ ARDP, ICMC-USP, Brazil

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Work by

◮ Eduardo Xavier Miqueles, Brazilian Synchrotron Light

Laboratory (LNLS)

◮ Elias Helou, ICMC-USP, Brazil ◮ ARDP, ICMC-USP, Brazil

Thanks to the Nuclear Instrumentation Laboratory (Federal University of Rio de Janeiro), Brazil, and the Brazilian Synchrotron Light Laboratory (LNLS), that provided the real data. Articles and Software: Physics in Medicine & Biology (2010), IEEE Transactions on Medical Imaging (2011), Studies in Applied Mathematics (2011), Computer Physics Communications (2011), Journal of Inverse and Ill-posed Problems (2012), Journal of Physics, RAFT.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Do You Want to Go Rafting? http://www.raftist.com/

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Synchrotron at the LNLS: a Light Source

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is a Light Source (http://www.lightsources.org/)

Accelerators producing intense beams of X-rays, ultra-violet and infrared light, making possible both basic and applied research in fields ranging from physics to biology and technology, which are not possible with more conventional equipment. Light refers to visible light, but also includes light with wavelengths that we cannot see, such as: radio waves, microwaves, infrared, ultraviolet, X-rays, and gamma rays.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is a Light Source

As a general rule,

◮ short-wavelength (hard) X-rays are most useful for probing

atomic structure

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is a Light Source

As a general rule,

◮ short-wavelength (hard) X-rays are most useful for probing

atomic structure

◮ long-wavelength (soft) X-rays and ultraviolet light are good

choices for studying chemical reactions

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is a Light Source

As a general rule,

◮ short-wavelength (hard) X-rays are most useful for probing

atomic structure

◮ long-wavelength (soft) X-rays and ultraviolet light are good

choices for studying chemical reactions

◮ Infrared is ideally suited to studying atomic vibrations in

molecules and solids

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is a Light Source

As a general rule,

◮ short-wavelength (hard) X-rays are most useful for probing

atomic structure

◮ long-wavelength (soft) X-rays and ultraviolet light are good

choices for studying chemical reactions

◮ Infrared is ideally suited to studying atomic vibrations in

molecules and solids

◮ long wavelength end (terahertz waves), it is also useful for

certain types of electronic structure experiments

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is a Light Source

As a general rule,

◮ short-wavelength (hard) X-rays are most useful for probing

atomic structure

◮ long-wavelength (soft) X-rays and ultraviolet light are good

choices for studying chemical reactions

◮ Infrared is ideally suited to studying atomic vibrations in

molecules and solids

◮ long wavelength end (terahertz waves), it is also useful for

certain types of electronic structure experiments X-rays is the range of the electromagnetic spectrum known as synchrotron light.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Ring

Electron gun (1-2). Electrons enter a circular-shaped booster ring (3), accelerated to relativistic speeds. Finally ,they enter another ring, a storage ring (4), where they circulate for hours. Special bending magnets help them keep to their circular path. Synchrotron light is produced when the electrons change direction around the ring when manipulated by bending magnets (5), entering beamlines (gates) (6).

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Gate (Beamline) at the LNLS

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Beamline (closer look)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What Can Be Done in a Light Source?

◮ biosciences (protein crystallography, cell biology)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What Can Be Done in a Light Source?

◮ biosciences (protein crystallography, cell biology) ◮ medical research (microbiology, disease mechanisms,

high-resolution imaging and cancer radiation therapy)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What Can Be Done in a Light Source?

◮ biosciences (protein crystallography, cell biology) ◮ medical research (microbiology, disease mechanisms,

high-resolution imaging and cancer radiation therapy)

◮ chemical and environmental sciences (toxicology, atmospheric

research, clean combustion and cleaner industrial production technologies)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What Can Be Done in a Light Source?

◮ biosciences (protein crystallography, cell biology) ◮ medical research (microbiology, disease mechanisms,

high-resolution imaging and cancer radiation therapy)

◮ chemical and environmental sciences (toxicology, atmospheric

research, clean combustion and cleaner industrial production technologies)

◮ agriculture (plant genomics, soil studies, animal and plant

imaging)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What Can Be Done in a Light Source?

◮ biosciences (protein crystallography, cell biology) ◮ medical research (microbiology, disease mechanisms,

high-resolution imaging and cancer radiation therapy)

◮ chemical and environmental sciences (toxicology, atmospheric

research, clean combustion and cleaner industrial production technologies)

◮ agriculture (plant genomics, soil studies, animal and plant

imaging)

◮ minerals exploration (rapid analysis of drill core samples,

comprehensive characterization of ores for ease of mineral processing)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What Can Be Done in a Light Source?

◮ biosciences (protein crystallography, cell biology) ◮ medical research (microbiology, disease mechanisms,

high-resolution imaging and cancer radiation therapy)

◮ chemical and environmental sciences (toxicology, atmospheric

research, clean combustion and cleaner industrial production technologies)

◮ agriculture (plant genomics, soil studies, animal and plant

imaging)

◮ minerals exploration (rapid analysis of drill core samples,

comprehensive characterization of ores for ease of mineral processing)

◮ Etc .........................

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Tomographic Problems

◮ X-rays Fluorescence Computed Tomography

(XFCT)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Tomographic Problems

◮ X-rays Fluorescence Computed Tomography

(XFCT)

◮ Polychromaticity (WP)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Tomographic Problems

◮ X-rays Fluorescence Computed Tomography

(XFCT)

◮ Polychromaticity (WP) ◮ High Resolution CT for Nanomaterials (WP+)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

XFCT

X-rays Fluorescence Computed Tomography

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Problem and the Data

◮ We want to reconstruct the concentration

distribution of a heavy metal (Copper, Zinc, Iron,..), or other element like Iodine, inside a body.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Problem and the Data

◮ We want to reconstruct the concentration

distribution of a heavy metal (Copper, Zinc, Iron,..), or other element like Iodine, inside a body.

◮ This concentration distribution could indicate

malignancy in a tissue, for example. Another application is determination of 3D rock structure in mineralogy.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Problem and the Data

◮ We want to reconstruct the concentration

distribution of a heavy metal (Copper, Zinc, Iron,..), or other element like Iodine, inside a body.

◮ This concentration distribution could indicate

malignancy in a tissue, for example. Another application is determination of 3D rock structure in mineralogy.

◮ Irradiation by high intensity monochromatic

synchrotron X rays at a specific energy of the element stimulates fluorescence emission (data).

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

X-Rays Fluorescence Computed Tomography (XFCT)

Aims at reconstructing fluorescence emitted by the body when bombarded by high intensity X-rays at a given energy (monochromatic).

θ

ξ ξ t

x

FLUORESCENCE DETECTOR OBJECT SOURCE TRANSMISSION DETECTOR

Ω

τ γ γ

1 2

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Generalized Attenuated Radon Transform

And the model is

d(t, θ) = RWf (t, θ) =

• x·ξ=t f (x)W (x, θ)dx

where f (x) is the emission (fluorescence) density at x, µ is the fluorescence attenuation, λ is the attenuation of the X-rays,

W (x, θ) = ωλ(x, θ)ωµ(x, θ), Dh(x, θ) =

• R h(x + qξ⊥)dq

ωµ(x, θ) =

• Γ e−Dµ(x,θ+γ)dγ , probability of reaching

the detector from x ωλ(x, θ) = e−Dλ(x,θ+π), survival probability of reaching x

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Inversion of Rω Alternatives:

◮ Just forget about W and make corrections ”a

posteriori”.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Inversion of Rω Alternatives:

◮ Just forget about W and make corrections ”a

posteriori”.

◮ Discretize and solve an optimization model. Too

computationally intensive . Not our option (so far).

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Inversion of Rω Alternatives:

◮ Just forget about W and make corrections ”a

posteriori”.

◮ Discretize and solve an optimization model. Too

computationally intensive . Not our option (so far).

◮ Approximate by a scaled Radon Inverse and

Iterate.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The Inversion of Rω Alternatives:

◮ Just forget about W and make corrections ”a

posteriori”.

◮ Discretize and solve an optimization model. Too

computationally intensive . Not our option (so far).

◮ Approximate by a scaled Radon Inverse and

Iterate.

◮ Try to find an analytic inverse. A nice

mathematical tour.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Iterative Inversion

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

First Option: Iterated Inversion

We have a reasonable (fast, accurate if there is not too much noise) inverse for R, so, let us try a fixed point iteration !!!!

f (k+1) = f (k) + e(k) = (I − 1

aR−1RW)f (k) + R−1 a d ,

e(k) =

R−1(d−RW f (k)) a

And what is a?

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Iterated Inversion

Clearly, convergence depends on how close 1

aR−1RW is to the

identity, equivalently, how well R−1 approximates R−1

W and this

depends on the attenuation. If it is too large, it will not work. Chang (IEEE TNS 78) suggested for SPECT a reasonable value for a as the average attenuation given by

a(x) = 1

2π

2π

0 W (x, θ)dθ.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

A Contraction Constant

The Contraction constant for Kf = 1 a

• R−1(R − RW )f − (1 − a)f
• = f − 1

aR−1RW f (1) is given by c(κa) = sup

u∈R2

sup

θ∈[0,2π]

• 1 −

1 2a∞ [W (u, θ) + W (u, θ + π)]

• (2)

Miqueles-DP, “On the Iterative Inversion of Generalized Attenuated Radon Transforms”, Journal of Inverse and Ill-posed Problems, 1569-3945, December 2012.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

And

Sequences c(κ) and c(κa) for different (increasing) values of attenuation µ, meaning that, we have a reasonable computable value measuring convergence rate and ill-conditioning.

2 4 6 8 10 0.6 0.7 0.8 0.9 1

j (a)

c(κ) c(κa) 2 4 6 8 10 0.2 0.4 0.6 0.8 1

j (b)

c(κ) c(κa) Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Classical Methods in a Continuous Setting: EM

The Expectation Maximization (EM) can also be applied to the linear part of our problem, assuming a known attenuation function. For the EM we have the following iteration (continuous version):

f (k+1)(x) = f (k)(x)BWd(k)(x) BWe(x) ,

where d(k)(t, θ) = d(t, θ)/RW f (k)(t, θ), BW d(x) = 2π W (x, θ)d(x · ξ, θ)dθ is the attenuated backprojection, and e = 1 in V.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Unknown Attenuation: Iteration Once Again

New Problem: Given d ∈ V, find {f , µ} ∈ U such that Y(f , µ) = RW (µ)f − d = 0 ∈ V. (3) Iterate: f (k+1) = L

• d, f (k), µ(k)

, µ(k+1) = N

• d, f (k+1), µ(k)

. L stands for an approximate inversion of RW given µ(k) and N for the application of (say) Newton‘s method to the equation for f (k+1) given. Miqueles-DP, “Iterative Reconstruction in X-ray Fluorescence Tomography Based on Radon Inversion”, IEEE Transactions on Medical Imaging, 30, 2, 2011.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Simulated Data

Figure below shows a 32x32 representation of functions {f , µ, λ} and the simulated attenuated Radon transform with 80 projections views and 60 rays per view.

Simulated data for XFCT. From left to right: density function f , fluorescence attenuation µ, transmission attenuation λ and attenuated Radon transform.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Real Data

A microscopic sample with a distribution of Copper and Zinc inside. For the Copper sample, each projection view had 23 rays, while 20 rays for the Zinc sample. The total number of views was 60 for both samples. Figure below shows the functions {Rµ, RW f }.

Real data. From left to right: transmission data for Cu sample, XFCT data for Cu sample, transmission data for Zn sample and XFCT data for Zn sample.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Simulated Data

AKT × EM with µ = λ and for iterations {1, 2, 3, 4, 20} (left to right). For each block, the EM reconstruction is shown in the first row and the AKT reconstruction in the second. Simulated case (32x32). No noise.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Real Data

AKT × EM with µ = λ and for iterations {1, 2, 3, 4, 20} (left to right). For each block, the EM reconstruction is shown in the first row and the AKT reconstruction in the second. Cu sample (60x60)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Real Data

AKT × EM with µ = λ and for iterations {1, 2, 3, 4, 20} (left to right). For each block, the EM reconstruction is shown in the first row and the AKT reconstruction in the second. Zn sample (60x60).

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Simulated and Real Data

0.05 0.10 0.2 0.4 0.2 0.4

Iterates {µ(0), µ(1)}, for AV, using AKT for f (1). Initial guess µ(0) was obtained using FBP of transmission data.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Analytic Approach

An Analytic Inverse

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse

Following Fokas, the spectral analysis of the differential equation: η · ∇u(x) + a(x, η)u(x) = f (x) η = η(κ) ∈ C2, κ ∈ C. allows us to to write the solution in terms of the GART. a = 0 leads to the Radon Transform, a(x) = µ(x) to the Attenuated Radon Transform and a(x, η) will be determined for the XFCT. In

• ur case

η(κ) = 1 2i 1 κ + κ

• , 1

2 1 κ − κ

• .

(4) η = o(κ) and each component of η is analytic in κ with a pole in zero.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse

Changing variables x → (z, ¯ z), with z ∈ C, defined by z = v · x and v = v(κ) ∈ C2. z ¯ z

• = Gx,

G = vT ¯ vT

• ,

J = −1 1

• (5)

The previous equation can be rewriten as (η · v)∂zu + (η · ¯ v)∂¯

zu + a(x)u(x) = f (x).

(6)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse

We choose vector {η, v} ∈ C 2 such that η · v = 0 and η · ¯ v = j(λ) = det G = −v · J¯

• v. Put v = −Jη and denote

η(κ) = (c(κ), b(κ))T, then j(κ) = c(κ)b(κ) − b(κ)c(κ) = 2iJ(κ), J(κ) = Imag

• c(κ)b(κ)
• .

(7) Choosing η and v so that η · v = 0 and η · ¯ v = j(κ) we get j(κ)∂¯

zu + a(x, η)u(x) = f (x).

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse

Multiplying by an Euler factor eq(x) we decouple the equation

• btaining two d-bar equations

∂¯

z

• u(x)eq(x)

= f (x) j(κ) eq(x), ∂¯

zq(x) = a(x)

j(κ) Define the singularity set S = {κ ∈ C: j(κ) = 0}

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse

And we can use for each equation the following (Fokas-Iserles, J.R.Soc.Interface, 3, 45-54, 2006.),

Lemma

For all κ ∈ S, the solution of the ∂¯

zˆ

u(x) = g(x)/j(κ) is given by ˆ u(x; κ) . = ∂−1

¯ z

g(x) j(κ)

• = α(κ)

2πi

• R2

g(y)dy v(κ) · (y − x). α(κ) = sign J(η). for g(x) = f (x)

j(κ) eq(x) or g(x) = a(x) j(κ) .

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

A Riemann-Hilbert problem

What is a scalar inhomogeneous Riemann-Hilbert (RH) problem? Given a closed contour S and H¨

• lder continuous functions f and g
• n S, find a sectionally analytic function Φ with finite degree at

infinity (Φ(z) ∼ cmzm + O(zm−1) as z → ∞, cm = 0, z / ∈ S) such that Φ+(t) = g(t)Φ−(t) + f (t) In our case g(t) = 1.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse: a Riemann-Hilbert problem

S determines a curve, dividing the complex plane into two regions R+ and R−, where d-bar equations determine u± for κ ∈ R±. The solution for all κ is determined by the jump J (x) = u+(x) − u−(x)

• n the curve S. Since u is an analytic function of κ ∈ S, there exist

z0 ∈ C and δ > 0 such that S is homotopic to a circle centered at z0 and radius δ. Assuming without loss of generality that δ = 1, the solution of our Riemann-Hilbert problem is given by (Ablowitz) u(x; κ) = 1 2πi

• |z−z0|=1

J (x) z − κdz = 1 2π 2π J (x)eiθ −1 κ + O 1 κ2

• dθ

= 1 κh(x) + O 1 κ2

• ,

h(x) = −1 2π 2π J (x)eiθdθ

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse

Therefore, from the original equation, and with the boundary condition u(x) = O(κ−1) as κ → ∞, we have f (x) = 1 κη(κ) · ∇h(x) + a(x)O 1 κ

• + O

1 κ2

• ,

κ → ∞ It only remains to compute the jump function J = J (x) in order to evaluate h in the above equation. And after many, many, ......., many, too many, ........ calculations .......... http://onlinelibrary.wiley.com/doi/10.1111/j.1467- 9590.2011.00527.x/abstract

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse

f (x) = 1 2π 2π iO(η, ξθ)

• eDa(x,θ)m{Ra, RW f }(x · ξθ, θ)
• dθ

(8) . = IηRW f (x) (9) giving rise to the inverse operator Iη. Where m{r, d} = e− r

2

• hc(r)H
• hc(r)e

r 2 d

• + hs(r)H
• hs(r)e

r 2 d

• (10)

where hc(r) = cos( 1

2H r) and hs(r) = sin( 1 2H r), and

O(η, ξθ) = [D(η)ξθ] · ∇ − i[D(η)ξ⊥

θ ] · ∇

(11) and matrix D(η) = diag (η1(κ)/κ, iη2(κ)/κ).

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse

Better

f (x) = 1 4π 2π ∂t

• eDa(x,θ)m{Ra, RWf }
• (x·ξθ, θ)dθ

(12)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Extending to XFCT. Single angle

First step is considering the case of a fixed angle and inverting the corresponding operator Rγ for a fixed γ ∈ Γ. The inspiration is the nonrealistic case where γ = π, the exponentials are “parallel“ and the solution of the problem is trivial, just considering a(x, θ) = λ(x, θ + π) + µ(x, θ + π) and Fokas approach applies in a straightforward manner. This suggests the necessity of considering a rotation of angle γ for the next step towards a generalization.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Extending to XFCT. Single angle. Rotation

For every given point x, the line ℓ(x, v) = {x + sv : s ≥ 0} can be mapped to the line ℓ(x, −v) through the rotation operator φx(y) = 2x − y (13) and also can be mapped to the line ℓ(x, Gγv) (for a fixed angle γ, being G T

γ = (ξγ, ξ⊥ γ ) a 2x2 rotation matrix) through the following

rotation ψγ,x(y) = Gγ(y − x) + x. (14)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Extending to XFCT. Single angle. Rotation

Rotations φx and ψγ,x The attenuation that we have to consider will be derived from aγ,x(y) = λ(φx(y)) + bγ,x(y) (15) with bγ,x a function defined by bγ,x(y) = µ(ψγ,x(y)). (16) And after several lemmas and lots of calculations, we get the inverse R−1

γ

for a fixed angle γ.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Series Inversion from the fixed angle approximation

Theorem

If the fluorescence attenuation map µ satisfies the inequality max

(x,θ) |1 − α(x, θ, β∗)| < 1, with α defined before, the inverse

• perator R−1

xfct is given by the Neumann series

R−1

xfct = 1

m

∞

• k=0
• I − 1

mR−1

β∗ Rxfct

k R−1

β∗

(17) with I the identity operator, β∗ = 1

2(γ1 + γ2), m = γ2 − γ1 and

R−1

β∗ from before.

In practical experiments, the angle section Γ, is symmetrically chosen to verify Γ ⊆ [0, π]. So, the optimal angle β∗ is π

2 and the condition above

is always satisfied since there is a minimum in the amount of scattered photons at π

2 and therefore the total fluorescence attenuation (the

divergent beam transform of µ) is stationary at this angle.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Series Inversion from the fixed angle approximation

Now, using the same change of variables as before, we define the function a, for fixed but arbitrary values of x, by a(x, η) = λ(φx(x)) + bη,x(x), (18) where x = φx(x) and b = bη,x(x) such that Dbη,x(x, η) = − ln 1 m

• Γ(x)

e−Dµ(x,Gγη)dγ, (19) Since Dµ is a positive function and

• Γ dγ e−Dµ < m, the above

logarithm is well defined, i.e., Db > 0.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse

If (t, ρ) is the change of variables in x = tξ + ρξ⊥, m = γ2 − γ1 ωxfct(x, θ) = e−Dλ(x,θ+π)

• Γ

e−Dµ(x,ξθ+γ)dγ. and p(t, θ) = Rλ(t, θ) + Rb(t, θ) and Rb defined by. Rb(x·ξθ, θ) = − ln 1 m2

• Γ

e−Dµ(x,θ+γ+π)dγ

Γ

e−Dµ(x,θ+γ)dγ

• .

and m{r, d} = e− r

2

• hc(r)H
• hc(r)e

r 2 d

• + hs(r)H
• hs(r)e

r 2 d

• (20)

with hc(r) = cos( 1

2H r) and hs(r) = sin( 1 2H r).

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

An Analytic Inverse: Finally

f (x) = 1 4π 2π ∂t

• mω−1

xfct(x, θ)m

• p, 1

mRxfctf

• (x · ξ, θ)
• dθ

(21) = 1 4π 2π ∂t

• ω−1

xfct(x, θ)m{p, Rxfctf }(x · ξ, θ)

• dθ

(22) = R−1

xfctRxfctf (x)

(23) Miqueles-DP, “On the analytic inversion of the Generalized Attenuated Radon Transform for X-ray Fluorescence Computed Tomography, Studies in Applied Mathematics, 127,4, 2011.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Simulated Data

Simulated data: {f1, λ1, µ1, Rxfctf1}, (256 × 256), sinograms

• btained with M = 360 views and N = 400 rays per view.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Simulated Data

Simulated data: {f2, λ2, µ2, Rxfctf2}, (80 × 80), sinograms

• btained with M = 360 views and N = 400 rays per view.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Simulated Data

xfct inversion. From left to right: true density map, R−1

xfctd and

(mRβ∗)−1d for f1

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Simulated Data

xfct inversion. From left to right: true density map, R−1

xfctd and

(mRβ∗)−1d for f1

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Some Experiments: Real Data

Sequence of partial sums of the approximating series for real data using µ = λ.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

So far, formulas are true if the incoming rays are monochromatic. If there is a monochromatic filter this is half true (± 5Kev). With white light (the whole spectrum), this is not true at all. So,.............................................................................................

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity p(η) =

• E
• Ω(η)

S(ε)f (x)Wµε,λ(x, η)dxdε

XFCT polychromatic sinogram p, given the emission f and linear attenuation coefficients µ and λ. Weight function W given by

Wµε,λ(x, θ) = e−Dµε(x,θ+π)

Γ e−Dλ(x,θ+γ)dγ,

and

Dh(x, θ) =

• R h(x + qξ⊥)dq

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

Equivalently

p(η) =

• Ω(η) f (x)ω(x, η)dx

(a) ω(x, η) =

• Γ e−Dλ(x,θ+γ)

E S(ε)e−Dµε(x,θ+π)dε

(b)

(a) Generalized Attenuated Radon Transform (GART) (b) Weight function for polychromatic XFCT.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity r(η) =

• E S(ε)e−Rµε(η)dε

Model for function µε = µ(·, ε) with S(ε) the energy spectrum. The equation above should be solved for µε, given the sinogram q(η) = − log r(η). Observe that

• E S(ε) = 1.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

Uz = Uz(ε) linear mass attenuation for Calcium in the range [1Kev,30Kev]. The jump close to 4 indicates a K level energy for that element.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

Contours for indicator functions {vj} at region Rj. Assume that µε(x) = J

j=1 Uzj(ε)vj(x). Region R4 is the main body of the

• bject1 with air within R8. Hence, vj(x) = 1Rj(x) and there is

Zinc, Nickel, Manganese, Iron, Potassium, Chlorine and Copper within regions {R5, R6, R7, R3, R1, R2}, respectively.

1We use 1A(x) as the indicator function over the set A. Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

(a)-(c): Functions µ(ε60, x), µ(ε180, x) and µ(ε380, x) respectively; where ε60 = 4Kev, ε180 = 10Kev and ε380 = 20Kev (d) Sinogram q(η) = − log r(η) (e) Sinogram r(η)

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

Map function R−1[− log r] with r given Constrast within each region is averaged by the FBP method.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

P[µε1, . . . , µεN](η) ≡

N

• j=1

S(εj)e−Rµεj (η)∆εj = r(η) Discrete model. Given r, we want to recover function µεj, ∀ j. We consider, for a given j, the equation P[µεj] = r

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

error= exp [−Rµε] − r/r . Uzj. S(ε) There exists τ ∈ E (MVT) such that P(µε) = exp [−Rµτ(η)]. τ (simulations) in the neighborhood of ε∗ = argmax{S(ε), ε ∈ E}. Two local minima, the first one very close to the maximum point

• f S (dashed line), i.e., near ε = 4.5Kev. Assumption:τ = ε∗

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

∂µε ∂ε < 0 since U′ z(ε) < 0. This is also true in general, i.e., the

higher the energy, the lower the attenuation. From this fact, and also from numerical experiments, we observed that

∂µε(x) ∂ε

≈ −Cµε(x) for some constant C. Let us assume C = 1. It follows from the mean value theorem that for ℓ ∈ (0, 1) µε(x) − µτ(x) = ∂µρ(x) ∂ε (ε − τ), ρ = (1 − ℓ)ε + ℓτ µε(x) − µτ(x) ≈ −µε(x)(ε − τ))

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

We approximate P as follows: e−Rµτ = e−RµεeR(µε−µτ) implying e−Rµτ ≈ e−Rµεe−(ε−τ)Rµε and therefore eRµτ ≈ e[1+(ε−τ)]Rµε. Hence, P[µε] ≈ ¶[µε] = eg(ε)Rµε with g(ε) = 1 + (ε − τ). Our new system of equations is ¶[µε] ≈ r,

• r equivalently, H[µε] ≡ ¶[µε] − r ≈ 0. We use Newton’s method

to solve it. The Fr´ echet derivative of H at u ∈ U, denoted by H′[u], is the linear transformation H′[u] = −g(ε)¶[u]R. This yields H′[u] · h = −g(ε)¶[u]Rh for any h ∈ U.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

(1) µ(k+1)

ε

(x) = µ(k)

ε (x) + z(k)(x),

ε ∈ E (2) H′ µ(k)

ε

• · zk = −H
• µ(k)

ε

• (3) z(k) =

1 g(ε)R−1

H[µ(k)

ε ]/¶[µ(k) ε ]

• (⋆) µ(k+1)

ε

= µ(k)

ε

+

1 g(ε)R−1

1 − eg(ε)Rµ(k)

ε r

• Sequence of steps to develop the reconstruction algorithm (⋆), for

a fixed ε. The final equation (⋆) Algorithm 1, is equivalent to a fixed point method for the Fredholm equation u(x) = ¯ u(x) + Du(x) associated to ¶[u] = r, where ¯ u = J−1[r], Du = J−1(r − ¶[u]) and J = H′(u). The convergence is guaranteed whenever J−1 is a good approximation for the inverse of P.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

Simulations, 256 × 256 phantom µε, and a polychromatic sinogram q = − log(r) obtained with 300 views and 300 rays at each view. The results are shown below. The stopping criterion was µ(k+1)

ε

− µ(k)

ε ≤ σ

with σ = 10−3.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Polychromaticity

Comparison between original phantom µε (left) and reconstructed µ(k)

ε

(right). (a) ε = ε100 = 5.95Kev and k = 30. (b) ε = ε80 = 4.95Kev and k = 15 . (c) ε = ε200 = 10.95Kev and k = 3

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Work in Progress

High Resolution CT for Nanomaterials

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

A Sinogram

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

A Reconstruction. What is it?

2000x2000 pixels, 1000 views, 2048 rays each, pixel size 1.64

• micrometers. White light, a ”narrow“ band around 20Kev.

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

A Toothpick!!!!!!!!!!!!!!

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is left: too many things that can be done

◮ Developing fast forward and backprojections for iterative

methods

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is left: too many things that can be done

◮ Developing fast forward and backprojections for iterative

methods

◮ Phase contrast tomography

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is left: too many things that can be done

◮ Developing fast forward and backprojections for iterative

methods

◮ Phase contrast tomography ◮ And for energies over (say) 30Kev?

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is left: too many things that can be done

◮ Developing fast forward and backprojections for iterative

methods

◮ Phase contrast tomography ◮ And for energies over (say) 30Kev? ◮ Etc,

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is left: too many things that can be done

◮ Developing fast forward and backprojections for iterative

methods

◮ Phase contrast tomography ◮ And for energies over (say) 30Kev? ◮ Etc, ◮ Etc,

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

What is left: too many things that can be done

◮ Developing fast forward and backprojections for iterative

methods

◮ Phase contrast tomography ◮ And for energies over (say) 30Kev? ◮ Etc, ◮ Etc, ◮ Etc......

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Thanks

Thanks

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

The First Team

Figure:

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation

Questions?

My best students

Alvaro R. De Pierro, University of S˜ ao Paulo, Department of Applied Mathematics and Statistics, Brazil Tomography with Synchrotron Radiation