Model based reconstruction for magnetic particle imaging in 2D and 3D

(1)

The Basic Model The MPI Core Operator Reconstruction in 2D and 3D Reconstruction Algorithm

Model based reconstruction for magnetic particle imaging in 2D and 3D

Tom M ¨arz, Andreas Weinmann

Canazei, September 2017

1 / 28

(2)

Overview

Review: The Basic Model

The MPI Core Operator

Reconstruction Formulae in 2D and 3D

Reconstruction Algorithm in 2D and 3D

(3)

Basic Principles

(courtesy of Knopp et al. 2010)

Data measured in MPI:voltageu(t) induced in the recording coils.

By Farraday’s law of induction,

u(t)=−d dtΦ(t), where the magnetic fluxΦ(t)is Φ(t) =µ0

Z

R³

R(x)(H(x,t) +M(x,t))dx, µ0magnet. permeability.

• The fluxΦ(t)is caused by the applied fieldH(x,t)and the magnetization responseM(x,t).

• R(x)∈R³^×³is the sensitivity pattern of the three recording coil pairs.

3 / 28

(4)

Basic Principles

Data measured in MPI:voltageu(t)given by

u(t)=−d

dtΦ(t), Φ(t) =µ0

Z

R³

R(x)(H(x,t) +M(x,t))dx,

with magnetic fluxΦ(t),applied fieldH(x,t),magnetization response M(x,t).

By the Langevin theory of paramagnetism for superparamagnetic nanoparticles,

M(x,t) =ρ(x)mL |H(x,t)|

Hsat

! H(x,t)

|H(x,t)|, L(x) =coth(x)−1 x, withL. . . Langevin function,m. . . magnetic moment of a single particle, Hsat. . . saturation parameter.

Signal to reconstruct in MPI:concentration of the particlesρ(x).

(5)

Simplification

Data:voltageu(t)given by u(t) =−µ0

d dt

Z

R³

R(x)(H(x,t) +M(x,t))dx.

Since the applied fieldHdoes not depend onρ,consider the data s(t) =−µ₀d

dt Z

R³

R(x)M(x,t)dx. In the interesting field of view,R is almost constant; hence

s(t) =−µ₀Rd dt

Z

R³

M(x,t)dx.

Reconstructthe particle densityρ(x)given via M(x,t) =ρ(x)mL |H(x,t)|

H_sat

! H(x,t)

|H(x,t)|.

Hence,

s(t) =−µ0mRd dt

Z

R³

ρ(x)L |H(x,t)| H_sat

! H(x,t)

|H(x,t)|dx.

5 / 28

(6)

Simplification

Data:voltageu(t)given by u(t)=−µ0

d dt

Z

R³

R(x)(H(x,t) +M(x,t))dx.

Since the applied fieldHdoes not depend onρ,consider thedata s(t)=−µ₀d

dt Z

R³

R(x)M(x,t)dx.

In the interesting field of view,R is almost constant; hence s(t) =−µ₀Rd

dt Z

R³

M(x,t)dx.

H_sat

! H(x,t)

|H(x,t)|.

Hence,

Z

R³

! H(x,t)

|H(x,t)|dx.

(7)

Simplification

d dt

Z

R³

R(x)(H(x,t) +M(x,t))dx.

dt Z

R³

R(x)M(x,t)dx.

In the interesting field of view,R is almost constant; hence s(t) =−µ0Rd

dt Z

R³

M(x,t)dx.

H_sat

! H(x,t)

|H(x,t)|.

Hence,

Z

R³

! H(x,t)

|H(x,t)|dx.

5 / 28

(8)

Simplification

d dt

Z

R³

R(x)(H(x,t) +M(x,t))dx.

dt Z

R³

R(x)M(x,t)dx.

In the interesting field of view,R is almost constant; hence s(t) =−µ0Rd

dt Z

R³

M(x,t)dx.

Reconstructthe particle densityρ(x)given via M(x,t) =ρ(x)mL |H(x,t)|

Hsat

! H(x,t)

|H(x,t)|.

Hence,

Z

ρ(x)L |H(x,t)|

H_sat

! H(x,t)

|H(x,t)|dx.

(9)

Simplification

Simplified problem: Reconstructρ(x)from voltage datas(t)related via s(t)

|{z}

data

= −µ0mR

| {z } constants

d dt

Z

R³

ρ(x)

|{z}

signal

L |H(x,t)|

Hsat

! H(x,t)

|H(x,t)|

| {z } dx.

The applied fieldHcosists of a static fieldH^Sand a dynamic fieldH^D: H(x,t) =H^S(x) +H^D(x,t) =Gx+P I(t)

In the interesting region,

H^S(x) =Gx=gdiag(−1,−1,2)x, H^D(x,t) =P I(t), g. . .nominal gradient of the static field,I(t). . .current in the coils, P. . .almost constant sensitivity profile of the drive field coils. Thefield free pointr(t)is given byH(r(t),t) =0. Then,

Gr(t) =−P I(t) ⇒ H(x,t) =Gx−P I(t) =−G(r(t)−x). Hence(Goodwill, Connolly),

s(t) =µ₀mRd dt

Z

R³

ρ(x)L |G(r(t)−x)| Hsat

! G(r(t)−x)

|G(r(t)−x)| dx.

6 / 28

(10)

Simplification

|{z}

data

= −µ0mR

| {z } constants

d dt

Z

R³

ρ(x)

|{z}

signal

L |H(x,t)|

Hsat

! H(x,t)

|H(x,t)|

| {z } dx.

The applied fieldHcosists of a static fieldH^Sand a dynamic fieldH^D: H(x,t) =H^S(x) +H^D(x,t)

=Gx+P I(t) In the interesting region,

H^S(x) =Gx=gdiag(−1,−1,2)x, H^D(x,t) =P I(t), g. . .nominal gradient of the static field,I(t). . .current in the coils, P. . .almost constant sensitivity profile of the drive field coils. Thefield free pointr(t)is given byH(r(t),t) =0. Then,

Z

R³

! G(r(t)−x)

|G(r(t)−x)| dx.

(11)

Simplification

|{z}

data

= −µ0mR

| {z } constants

d dt

Z

R³

ρ(x)

|{z}

signal

L |H(x,t)|

Hsat

! H(x,t)

|H(x,t)|

| {z } dx.

The applied fieldHcosists of a static fieldH^Sand a dynamic fieldH^D: H(x,t) =H^S(x) +H^D(x,t)=Gx+P I(t)

H^S(x) =Gx=gdiag(−1,−1,2)x, H^D(x,t) =P I(t), g. . .nominal gradient of the static field,I(t). . .current in the coils, P. . .almost constant sensitivity profile of the drive field coils.

Thefield free pointr(t)is given byH(r(t),t) =0. Then,

Z

R³

! G(r(t)−x)

|G(r(t)−x)| dx.

6 / 28

(12)

Simplification

|{z}

data

= −µ0mR

| {z } constants

d dt

Z

R³

ρ(x)

|{z}

signal

L |H(x,t)|

Hsat

! H(x,t)

|H(x,t)|

| {z } dx.

Thefield free pointr(t)is given byH(r(t),t) =0.

Then,

Z

R³

! G(r(t)−x)

|G(r(t)−x)| dx.

(13)

Simplification

|{z}

data

= −µ0mR

| {z } constants

d dt

Z

R³

ρ(x)

|{z}

signal

L |H(x,t)|

Hsat

! H(x,t)

|H(x,t)|

| {z } dx.

The applied fieldHcosists of a static fieldH^Sand a dynamic fieldH^D: H(x,t)=H^S(x) +H^D(x,t) =Gx+P I(t)

Thefield free pointr(t)is given byH(r(t),t) =0. Then, Gr(t) =−P I(t)

⇒ H(x,t) =Gx−P I(t) =−G(r(t)−x). Hence(Goodwill, Connolly),

Z

R³

! G(r(t)−x)

|G(r(t)−x)| dx.

6 / 28

(14)

Simplification

|{z}

data

= −µ0mR

| {z } constants

d dt

Z

R³

ρ(x)

|{z}

signal

L |H(x,t)|

Hsat

! H(x,t)

|H(x,t)|

| {z } dx.

The applied fieldHcosists of a static fieldH^Sand a dynamic fieldH^D: H(x,t)=H^S(x) +H^D(x,t) =Gx+P I(t)

Gr(t) =−P I(t) ⇒ H(x,t) =Gx−P I(t) =−G(r(t)−x).

Hence(Goodwill, Connolly), s(t) =µ₀mRd

dt Z

R³

! G(r(t)−x)

|G(r(t)−x)| dx.

(15)

Simplification

|{z}

data

= −µ0mR

| {z } constants

d dt

Z

R³

ρ(x)

|{z}

signal

L |H(x,t)|

Hsat

! H(x,t)

|H(x,t)|

| {z } dx.

Gr(t) =−P I(t) ⇒ H(x,t) =Gx−P I(t) =−G(r(t)−x).

Hence(Goodwill, Connolly), s(t) =µ₀mRd

dt Z

R³

ρ(x)L |G(r(t)−x)|

Hsat

! G(r(t)−x)

|G(r(t)−x)| dx.

6 / 28

(16)

II. The MPI Core Operator

(or, getting rid of particular trajectories)

(17)

Decomposing the model

Problem: Reconstructρ(x)from voltage datas(t)related via s(t)

|{z}

data

= −µ0mR

| {z } constants

d dt

Z

R³

ρ(x)

|{z}

signal

L |G(r(t)−x)|

Hsat

! G(r(t)−x)

|G(r(t)−x)|

| {z }

kernel

dx.

From a mathematical viewpoint, by transformation,

¯s(t) = d dt

Z

R³

¯ρ(bx)L |br(t)−bx| h

!

br(t)−bx

|br(t)−bx| dbx,

s(t) = d dt

Z

R³

ρ(x)L |r(t)−x| h

! r(t)−x

|r(t)−x| dx, Or,

s(t) =∇_rΦ(r) ˙r(t), where Φ(r) =Z

Rⁿ

ρ(x)L |r−x| h

! r−x

|r−x| dx.

8 / 28

(18)

Decomposing the model

Problem: Reconstructρ(x)from voltage datas(t)related via s(t)

|{z}

data

= −µ0mR

| {z } constants

d dt

Z

R³

ρ(x)

|{z}

signal

L |G(r(t)−x)|

Hsat

! G(r(t)−x)

|G(r(t)−x)|

| {z }

kernel

dx.

From a mathematical viewpoint, by transformation,

¯s(t) = d dt

Z

R³

¯ρ(bx)L |br(t)−bx| h

!

br(t)−bx

|br(t)−bx| dbx,

s(t) = d dt

Z

R³

ρ(x)L |r(t)−x|

h

! r(t)−x

|r(t)−x| dx, Or,

s(t) =∇_rΦ(r) ˙r(t), where Φ(r) = Z

Rⁿ

ρ(x)L |r−x|

h

! r−x

|r−x| dx.

(19)

Decomposing the model

Problem: Reconstructρ(x)from voltage datas(t)related via s(t) =∇_rΦ(r) ˙r(t), where Φ(r) =Z

Rⁿ

ρ(x)L |r−x|

h

! r−x

|r−x| dx.

=⇒ The signalsonly depends on the locationrand the velocityr˙of the field free point, and not on the particular trajectory.

• Application: Plugging together different trajectories, overlapping fields of view,. . .

• Mathematically: view MPI as an operator ρ→A_h[ρ](r,v)

whereAh[ρ]is a function on phase space, linear in the velocityv, Ah[ρ](r)v =∇_rΦ(r)v =

Z

Rⁿ

ρ(x)∇_r L |r−x| h

! r−x

|r−x|

! dx·v.

The MPI core operatorAhis independent of a particular trajectory.

9 / 28

(20)

Decomposing the model

Problem: Reconstructρ(x)from voltage datas(t)related via s(t) =∇_rΦ(r) ˙r(t), where Φ(r) =Z

Rⁿ

ρ(x)L |r−x|

h

! r−x

|r−x| dx.

Z

Rⁿ

ρ(x)∇_r L |r−x| h

! r−x

|r−x|

! dx·v.

(21)

Decomposing the model

Problem: Reconstructρ(x)from voltage datas(t)related via s(t) =∇_rΦ(r)r(t),˙ where Φ(r) =Z

Rⁿ

ρ(x)L |r−x|

h

! r−x

|r−x| dx.

Z

Rⁿ

ρ(x)∇_r L |r−x|

h

! r−x

|r−x|

! dx·v.

9 / 28

(22)

Decomposing the model

Problem: Reconstructρ(x)from voltage datas(t)related via s(t) =∇_rΦ(r)r(t),˙ where Φ(r) =Z

Rⁿ

ρ(x)L |r−x|

h

! r−x

|r−x| dx.

Z

Rⁿ

h

! r−x

|r−x|

! dx·v.

(23)

III. Reconstruction in 2D and 3D

10 / 28

(24)

Idealization limit h → 0

The MPI core operatorAhis given by Ah[ρ](r)v =

Z

Rⁿ

h

! r−x

|r−x|

! dx·v.

What happens ifh→0?

• Physical meaning: e.g., temperature decreases, or, particle size increases.

• In 1D: kernel tends to Dirac pulse.

• Idealized operator without blurring part.

(25)

Idealization limit h → 0

The MPI core operatorAhis given by Ah[ρ](r)=

Z

Rⁿ

h

! r−x

|r−x|

! dx.

Theorem (

M ¨arz, W., IPI, 2016

)

Letαh[ρ](r) =traceAh[ρ](r)and let α[ρ](r) =limh→0αh[ρ](r). Then, α[ρ](r) =Z

Rⁿ

ρ(x)κ(r−x)dx, ρin BV0(Ω). In dimension n>1,

κ(r−x) :=divr

r−x

|r−x|

!

= n−1

|r−x|.

• In 1D, we have a Dirac-kernelκ(r−x) =2δ(r−x).A peaking property was conjectured for nD which is not true by the theorem.

• Striking point: the traceα[ρ]already contains all information onρ. This has not been realized before.

12 / 28

(26)

Idealization limit h → 0

Z

Rⁿ

h

! r−x

|r−x|

! dx.

Theorem (

M ¨arz, W., IPI, 2016

)

Rⁿ

κ(r−x) :=divr

r−x

|r−x|

!

= n−1

|r−x|.

• Striking point: the traceα[ρ]already contains all information onρ. This has not been realized before.

(27)

Idealization limit h → 0

Z

Rⁿ

h

! r−x

|r−x|

! dx.

Theorem (

M ¨arz, W., IPI, 2016

)

Rⁿ

κ(r−x) :=divr

r−x

|r−x|

!

= n−1

|r−x|.

• Striking point:the traceα[ρ]already contains all information onρ.

This has not been realized before.

12 / 28

(28)

Relation to the Laplace equation

We consider the MPI core operatorAhand its trace α_h[ρ](r) =traceA_h[ρ](r), together with the idealization limit

α[ρ](r) = lim

h→0αh[ρ](r), α[ρ](r) =Z

Rⁿ

ρ(x)κ(r−x)dx.

Corollary (

M ¨arz, W., IPI, 2016

)

We have the following dimension-dependent relations of the kernelκto the fundamental solutionΦ(r,x)of the Laplace equ.−∆rΦ =δ(r−x),

in 3D, κ(r−x) =8πΦ(r,x) with Φ(r,x) = 1 4π|r−x|; in 2D, κ(r−x) =−2π∇_rΦ(r,x)· r−x

|r−x|; with Φ(r,x) =−1

2πlog(|r−x|).

d²

(29)

Reconstruction Formula

Corollary (Reconstruction Formula for the Idealized Case,

M ¨arz, W., IPI, 2016

)

Consider the idealized MPI core operator A₀[ρ](r) = lim

h→0A_h[ρ](r) and

data F(r) =A₀[ρ](r)at each point r given for the idealized scenario.

Then,

ρ=κ⁻¹◦traceA0[ρ],

whereκis the dimension-dependent convolution kernel from above.

That means, we take the pointwise trace and then deconvolve w.r.t. the dimension-dependentκ. In particular, in 3D,

ρ= 1

8π∆◦traceA0[ρ],

14 / 28

(30)

Reconstruction Formula

Corollary (Reconstruction Formula for the Idealized Case,

M ¨arz, W., IPI, 2016

)

h→0A_h[ρ](r) and

Then,

That means, we take the pointwise trace and then deconvolve w.r.t. the dimension-dependentκ.

In particular, in 3D, ρ= 1

(31)

Reconstruction Formula

Corollary (Reconstruction Formula for the Idealized Case,

M ¨arz, W., IPI, 2016

)

h→0A_h[ρ](r) and

Then,

That means, we take the pointwise trace and then deconvolve w.r.t. the dimension-dependentκ. In particular, in 3D,

ρ= 1

14 / 28

(32)

Ill-Posedness

Corollary (Ill-Posedness

M ¨arz, W., IPI, 2016

)

Even the idealized MPI problem is ill-posed in 2D and 3D. Depending on the dimension the degree of ill-posedness, i.e., the order of gained Sobolev smoothness of the forward operator, is one in 2D, and two in 3D.

Remark.This is not the case in 1D whereκis a Dirac distribution.

(33)

Ill-Posedness

Corollary (Ill-Posedness

M ¨arz, W., IPI, 2016

)

Even the idealized MPI problem is ill-posed in 2D and 3D. Depending on the dimension the degree of ill-posedness, i.e., the order of gained Sobolev smoothness of the forward operator, is one in 2D, and two in 3D.

Remark.This is not the case in 1D whereκis a Dirac distribution.

15 / 28

(34)

Non-idealized situation

We consider the non-idealized situationh>0 now.

Theorem (

M ¨arz, W., IPI, 2016

)

The traceαh[ρ]of the MPI core operator Ah[ρ]is given by α_h[ρ](r) =Z

Rⁿ

ρ(x)κ_h(r−x)dx, ρ∈BV0(Ω),where the convolution kernelκhis given by

κ_h(y) = 1 hf |y|

h

!

, with f(z) =L⁰(z) +L(z)n−1

z , (1)

Lthe Langevin function;

f is analytic with only purely imaginary singularities; near zero, f has the power series expansion

f(z) =

∞

X

k=0

2^2k⁺²B2k+2

(2k+2)! (2k+n)z^2k, B_lthe l-th Bernoulli number, with a convergence radius ofπ. Thusκhis analytic.

(35)

Non-idealized situation

We consider the non-idealized situationh>0 now.

Theorem (

M ¨arz, W., IPI, 2016

)

The traceαh[ρ]of the MPI core operator Ah[ρ]is given by α_h[ρ](r) =Z

Rⁿ

ρ(x)κ_h(r−x)dx, ρ∈BV0(Ω),where the convolution kernelκhis given by

κ_h(y) = 1 hf |y|

h

!

, with f(z) =L⁰(z) +L(z)n−1

z , (1)

Lthe Langevin function; f is analytic with only purely imaginary singularities; near zero, f has the power series expansion

f(z) =

∞

X

k=0

2^2k⁺²B2k+2

(2k+2)! (2k+n)z^2k, B_lthe l-th Bernoulli number, with a convergence radius ofπ. Thusκhis analytic.

16 / 28

(36)

Reconstruction Formula

Corollary (Reconstruction Formula for the the Non-Idealized Case,

M ¨arz, W., IPI, 2016

)

Consider the MPI core operator Ah[ρ](r)and suppose that data F(r) =Ah[ρ](r)at each point r is given. Then,

κ_h∗ρ=traceAh[ρ],

whereκhis the analytic convolution kernelκh(y) = ¹_hf_|_y_|

h

with f(z) =L⁰(z) +L(z)ⁿ⁻_z¹from the slide before.

That means, we take the pointwise trace and then deconvolve w.r.t.κh.

(37)

Reconstruction Formula

Corollary (Reconstruction Formula for the the Non-Idealized Case,

M ¨arz, W., IPI, 2016

)

Consider the MPI core operator Ah[ρ](r)and suppose that data F(r) =Ah[ρ](r)at each point r is given. Then,

κ_h∗ρ=traceAh[ρ],

whereκhis the analytic convolution kernelκh(y) = ¹_hf_|_y_|

h

with f(z) =L⁰(z) +L(z)ⁿ⁻_z¹from the slide before.

That means, we take the pointwise trace and then deconvolve w.r.t.κh.

17 / 28

(38)

Ill-Posedness

Corollary (Ill-Posedness

M ¨arz, W., IPI, 2016

)

The non-idealized MPI problem is severely ill-posed in the following sense: there are no two spaces H^s,H^tin the (Hilbert-)Sobolev scale, such that the traceα_hof the MPI operator induces an isomorphism α_h:H^s→H^t between these two spaces.

(39)

There is a particularly nice relation between the tracesα[ρ]andαh[ρ]of the idealized and the non-idealized MPI operators in 3D.

Theorem (

M ¨arz, W., IPI, 2016

)

In 3D, we have that

αh[ρ] =−∆κh

8π ∗α[ρ].

This tells us that in 3D the non-idealizedαh[ρ]is a massively smoothed version of the idealizedα[ρ]. Recall that

κ_h(y) = 1 hf |y|

h

!

, with f(z) =

∞

X

k=0

2^2k⁺²B_2k+2

(2k+2)! (2k+n)z^2k.

19 / 28

(40)

IV. Reconstruction Algorithm

(41)

Discretization and Given Data

We use the derived reconstruction formulae to design a reconstruction algorithm for MPI in 2D/3D.

Measured data: Samples sk =s(tk) of time datas(t)associated with a scan trajectoryr(t).

• Recall:r(t)is related with the electrical currentI(t)via r(t) =−G⁻¹P I(t),

G=g diag(−1,−1,2),g. . .nominal gradient of the static field,P. . . sensitivity profile of the drive field coils.

• The measured data are a discrete sampling of the MPI core operator applied toρ

s_k =s(t_k) =A_h[ρ](r_k)v_k,

at locationrk =r(tk),with the trajectory having tangentvk =v(tk) at timetk,for finitely many measurements indexed byk.

• Note: the presented approach is independent of the particular trajectory type employed.

21 / 28

(42)

Discretization and Given Data

Measured data: Samples sk =s(tk) of time datas(t)associated with ascan trajectoryr(t).

• Recall:r(t)is related withthe electrical currentI(t)via r(t)=−G⁻¹PI(t),

s_k =s(t_k) =A_h[ρ](r_k)v_k,

at locationrk =r(tk),with the trajectory having tangentvk =v(tk) at timetk,for finitely many measurements indexed byk.

(43)

Discretization and Given Data

Measured data: Samples sk =s(tk) of time datas(t)associated with a scan trajectoryr(t).

• Recall:r(t)is related with the electrical currentI(t)via r(t) =−G⁻¹P I(t),

s_k =s(t_k) =A_h[ρ](r_k)v_k,

at locationr_k =r(t_k),with the trajectory having tangentv_k =v(t_k) at timetk,for finitely many measurements indexed byk.

21 / 28

(44)

Discretization and Given Data

Measured data: Samples sk =s(tk) of time datas(t)associated with a scan trajectoryr(t).

• Recall:r(t)is related with the electrical currentI(t)via r(t) =−G⁻¹P I(t),

s_k =s(t_k) =A_h[ρ](r_k)v_k,

at locationr_k =r(t_k),with the trajectory having tangentv_k =v(t_k) at timetk,for finitely many measurements indexed byk.

• Note: the presented approach is independent of the particular

(45)

Major Algorithmic Steps

Measured dataare a discrete sampling of the MPI core operator applied toρ

sk =s(tk) =Ah[ρ](rk)vk,

at locationrk =r(tk),with the trajectory having tangentvk =v(tk).

Reconstruction formula:

κ_h◦ρ=traceA_h[ρ],

Our scheme may be subdivided intotwo major steps.

Step 1:Deriving trace data on a spacial grid from the raw input. We obtain a grid functionurepresenting trace data, i.e.,

u≈traceA_h[ρ] in each pixel (grid cell).(Details follow.)

Step 2:Reconstruction of the signal from the derived trace data by deconvolution (ill-posed). Regularized solution of the problem

Find ρ in κ_h∗ρ=u, givenu.(Details follow.)

22 / 28

(46)

Major Algorithmic Steps

κ_h◦ρ=traceA_h[ρ], Our scheme may be subdivided intotwo major steps.

u≈traceA_h[ρ]

in each pixel (grid cell).(Details follow.)

(47)

Major Algorithmic Steps

κ_h◦ρ=traceA_h[ρ], Our scheme may be subdivided intotwo major steps.

u≈traceA_h[ρ]

in each pixel (grid cell).(Details follow.)

22 / 28

(48)

Example.

Noisy time signal – cut-out Intermediate: trace after spatial fitting.

red: component 1, green: component 2

(49)

Trace Data on a Spacial Grid from the Raw Input

Measured data are samples of the MPI core operator applied toρ, s_k =s(t_k) =A_h[ρ](r_k)v_k,

at locationrk =r(tk),tangentvk =v(tk).Reconstruction formula:

κh◦ρ=traceAh[ρ].

Step 1: Deriving trace data on a spacial grid from the raw input.

Find a grid functionu

u≈traceAh[ρ]defined on a grid ofN1×N2cells.

• On each cell,uis constant; each celliis represented by its center pointxi.

• A time samplet_k belongs to celli,ifr(t_k)is in this cell; For each cell i,we collect the signal datas(t_k_i)from samplest_k_i belonging to the celliand gather them in a matrixS_i. Accordingly, we collect the velocity vectorsr(t˙ _k_i) =v(t_k_i)and gather them in a matrixV_i.

• We obtain the matrix fitting problem w.r.t.Ai,

Ai Vi =Ah[ρ](xi)Vi =Si. (2) which we solve by least squares fitting. Then,ui :=traceAi.

24 / 28

(50)

Signal Reconstruction from the Trace Data by Deconvolution

Reconstruction formula:

κ_h◦ρ=traceA_h[ρ].

Step 1 yields a grid functionu

u≈traceAh[ρ]defined on a grid ofN1×N2cells.

Find ρ in κh∗ρ=u, givenuby classical Tychonov regularization

ρ=arg min

bρ µkDbρk²₂+kKhbρ−uk²₂. (3) We solve the corresponding discrete Euler-Lagrange equation

−µLρ+Kh(Khρ−u) =0. (4) Here,−L =D^TDis the five point stencil discretization of the Laplacian

(51)

Reconstruction Algorithm - Summary

input : Time dependent samplessk =s(tk)along trajectoryrk =r(tk)with tangentvk = ˙r(tk)at timestk;regularization parameterµ.

output: Reconstructed particle densityρ. fork ←1toKdo

// Associate time samples with pixel grid.

ifr_k in cell ithen

V(i)←[V(i), v_k]; //Append tangent direction.

S(i)←[S(i), sk]; // Append data value.

end end

fori←1toIdo

// For each cell fit trace data using (2). [Qi,Ri]←QR(V_i^T);

Ai←SiQiR_i^T; ui←traceAi; end

// Regularized deconvolution of the trace data using (3) by // solving (4) with conjugate gradients (CG).

ρ=CG(−µL+K_h², Khu);

26 / 28

(52)

Summary

• We have reviewed the MPI model.

• We have extracted the MPI core operator.

• We have analyzed the idealized situation and the non-idealized situation.

• We have obtained reconstruction formulae for both cases based on matrix traces of the MPI core operator.

• We have seen that even the idealized MPI problem is ill-posed in 2D and 3D, which contrasts the 1D situation.

• We have seen that the MPI problem is severly ill-posed.

• We have derived a reconstruction algorithm based on the reconstruction formulae.

(53)

Some References

B. Gleich and J ¨urgen Weizenecker.

Tomographic imaging using the nonlinear response of magnetic particles.

Nature, 435:1214–1217, 2005.

T. Knopp and Thorsten Buzug.

Magnetic Particle Imaging: An Introduction to Imaging Principles and Scanner Instrumentation.

Springer, 2012.

P. Goodwill and S. Conolly.

The X-space formulation of the magnetic particle imaging process: 1-D signal, resolution, bandwidth, SNR, SAR, and magnetostimulation.

IEEE Transactions on Medical Imaging, 29:1851–1859, 2010.

P. Goodwill and S. Conolly.

Multidimensional X-space magnetic particle imaging.

IEEE Transactions on Medical Imaging, 30:1581–1590, 2011.

T. M ¨arz and A. Weinmann.

Model-based reconstruction for magnetic particle imaging in 2D and 3D.

Inverse Problems and Imaging, 10:1087 – 1110, 2016.