Phase Space Analysis in Medical Imaging Chris Stucchio Courant - - PowerPoint PPT Presentation

▶
phase space analysis in medical imaging
SMART_READER_LITE
LIVE PREVIEW

Phase Space Analysis in Medical Imaging Chris Stucchio Courant - - PowerPoint PPT Presentation

Aug 21, 2023 278 likes •1.46k views

Phase Space Analysis in Medical Imaging Chris Stucchio Courant Inst. and Trading Games, Inc. Collaboration with L. Greengard. Magnetic Resonance Imaging Excellent soft tissue contrast. No radiation. 2003 Nobel Prize (Lauterberger,

slide-1
SLIDE 1

Phase Space Analysis in Medical Imaging

Chris Stucchio Courant Inst. and Trading Games, Inc. Collaboration with L. Greengard.

slide-2
SLIDE 2

Magnetic Resonance Imaging

Excellent soft tissue contrast. No radiation. 2003 Nobel Prize (Lauterberger, Mansfield). Damadian maybe deserves credit too?

slide-3
SLIDE 3

MY LATERAL SPINE

slide-4
SLIDE 4

Objectives and Challenges

Show radiologist accurate pictures Quantify anatomical features

GOALS

Noise Artifacts Ambiguity

CHALLENGES

slide-5
SLIDE 5

Accurate Pictures

slide-6
SLIDE 6

Segment Anatomical Features

Separate into distinct regions

slide-7
SLIDE 7

Identification

Label the segmented regions

SKULL STARBOARD OVAL BRAIN PART PORT OVAL BRAIN PART GREAT BIG BRAIN TUMOUR

slide-8
SLIDE 8

Diagnosis

Draw conclusion from image data

SKULL STARBOARD OVAL BRAIN PAR PORT OVAL BRAIN PART GREAT BIG BRAIN TUMOUR

PATIENT IS SICK

slide-9
SLIDE 9

How an MRI works

slide-10
SLIDE 10

How an MRI works

Big Magnet: 1-2 Tesla Nucleus of atoms has spin Level Splitting: magnetic field breaks spin symmetry

slide-11
SLIDE 11

How an MRI Works

Excited state decays to ground, emits radiation. Measuring the radiation gives information on

  • bject.

SPIN UP (EXCITED STATE) SPIN DOWN (GROUND STATE)

slide-12
SLIDE 12

How an MRI works

Bloch Equation (macroscopic model): M(t) is magnetization, B(t) the magnetic field. M0(x) = Cρ(x) M(x, 0) = M0(x)

∂t M(x, t) = γ M × B(x, t) − P1,2 M T2 − P3( M(x, t) − M0(x)) T1

slide-13
SLIDE 13

HOW AN MRI WORKS

COORDINATE SYSTEM

Z Y X

slide-14
SLIDE 14

How an MRI works

Hit system with weak RF pulse (excitation): Rotates spins from z-direction into x-y plane

  • B(x, t) = [0, f(t)w(xz), 0]
  • M(x, t) ×

B(x, t) = [0, 0, M0(x)] × [0, f(t)w(x3), 0] = [−M0(x)f(t)w(x3), 0, 0]

slide-15
SLIDE 15

How an MRI works

Switch off excitation pulse, use probe field: X-Y components decoupled from Z component Substitution: M(t) =

Mx(t) + i My(t)

  • B(t) = [B0 +

G(t) · [x1, x2, 0]T ] z

slide-16
SLIDE 16

How an MRI works

Bloch Equation: M(x, t0) = ρ(x)w(xz)h(t0) ∂tM(x, t) =

  • −iγ(B0 + G(t) · x) − 1

T2

  • M(x, t)
slide-17
SLIDE 17

How an MRI works

Use RF receiver coils measure emission in the sample.

S(t) ∼

  • M(x, t)dx + noise
slide-18
SLIDE 18

How an MRI works

Solution: Simplify:

M(x, t) = ρ(x)w(x3)e−iγB0te−iγ(

R t

t0 G(t′)dt′)·xe−t/T2

  • k(t)

= γ t

t0

G(t′)dt′ M(x, t) → eiγB0tM(x, t)

slide-19
SLIDE 19

How an MRI works

Solution: Signal:

M(x, t) = ρ(x)w(x3)e−ik(t)·xe−t/T2

S(t) ∼ e−t/T2

  • ρ(x)w(x3)e−ik(t)·xdx + noise
slide-20
SLIDE 20

How an MRI works

Signal: An MRI measures the Continuous Fourier Transform of the density.

S(t) ∼ e−t/T2 ˆ ρ(k(t)) + noise

slide-21
SLIDE 21

Image Reconstruction

ACCURATE PICTURES

slide-22
SLIDE 22

Fourier Inversion

Hugely ill posed problem. Then:

∃f(x) = 0, [ρ + f](k1,...,N) = ˆ ρ(k1,...,N)

Given ˆ ρ(k1), . . . , ˆ ρ(kN), find ρ(x)

slide-23
SLIDE 23

Fourier Inversion

Fourier’s Theorem. Assume Cartesian sampling. Best approximation to density in norm

ρ(x) ≈

  • n

ˆ ρ(2π n)e−i2π

nx

L2([0, 1]2)

slide-24
SLIDE 24

Fourier Inversion

Fourier Transform not convergent pointwise Regularization discards information

ρ(x) ≈

  • n

ˆ ρ(2π n)w( n)e−i2π

nx

slide-25
SLIDE 25

FOURIER INVERSION

slide-26
SLIDE 26

OTHER ARTIFACTS

SMALL CURVATURE POSES PROBLEMS

slide-27
SLIDE 27

Current Solution

Reconstruct image using regularized discrete Fourier transform: Clean up regularized image in x-domain. Segment/identify based on cleaned up image.

ρ(x) ≈

  • n

ˆ ρ(2π n)w( n)e−i2π

nx

slide-28
SLIDE 28

Segmentation

OUTLINING THE IMPORTANT FEATURES

slide-29
SLIDE 29

Segmentation

Segmentation by anatomy/composition - outline the cancerous part Segmentation by perception - draw the same

  • utlines as a human

Image-space segmentation - separate based on image boundaries

GOALS

slide-30
SLIDE 30

Image boundaries

Image boundaries are places where image composition changes sharply. In medical images, this happens at discontinuities

  • f image.

Not true in other modalities.

slide-31
SLIDE 31

Discontinuities

Want to find discontinuities of an image. Image domain methods fail due to artifacts. Want to find discontinuities from raw MRI data, i.e. from samples of Fourier transform of image.

slide-32
SLIDE 32

Discontinuities

Simple model: a 1-d function with a discontinuity: If we localize on high frequencies, we can extract edges.

  • eikxf(x)dx = eikx0 f(x+

0 ) − f(x− 0 )

ik + O(k−2)

slide-33
SLIDE 33

1D Edge Detection

Laplace Filters, Gradient Filters, Concentration Kernels, etc.

STATE OF THE ART: CONCENTRATION KERNELS, C.F. TADMOR/GELB/ETC

slide-34
SLIDE 34

2D Edge Detectors

Tensor Products Radial Variables DFT−1[h( k)ˆ ρ( k)] R2 = R ⊗ R

slide-35
SLIDE 35

RESULT OF HIGH FREQUENCY FILTERS

EDGE DETECTOR RESOLUTION IS HALF THAT OF IMAGE

slide-36
SLIDE 36

RESULT OF HIGH FREQUENCY FILTERS

EDGE DETECTOR RESOLUTION IS HALF THAT OF IMAGE

slide-37
SLIDE 37

RESULT OF HIGH FREQUENCY FILTERS

EDGE DETECTOR RESOLUTION IS HALF THAT OF IMAGE

slide-38
SLIDE 38

RESULT OF HIGH FREQUENCY FILTERS

EDGE DETECTOR RESOLUTION IS HALF THAT OF IMAGE

slide-39
SLIDE 39

RESULT OF HIGH FREQUENCY FILTERS

EDGE DETECTOR RESOLUTION IS HALF THAT OF IMAGE

slide-40
SLIDE 40

RESULT OF HIGH FREQUENCY FILTERS

EDGE DETECTOR RESOLUTION IS HALF THAT OF IMAGE

slide-41
SLIDE 41

Problems

Noisy Does not separate regions Not obvious how to “fill in the holes”

slide-42
SLIDE 42

Problems

Noisy Does not separate regions Not obvious how to “fill in the holes”

slide-43
SLIDE 43

Problems

Noisy Does not separate regions Not obvious how to “fill in the holes”

slide-44
SLIDE 44

Boundary Reconstruction

Combinatorial methods Active Contours/Snakes/Level Sets Bayesian Methods

slide-45
SLIDE 45

Combinatorial Methods

Delaunay-methods: find the Crust of a point-set. Start with Delaunay graph. If a disk touches both ends of an edge in the Delaunay graph also touches a third vertex, then delete the edge.

(AMENTA, BERN, DEY, KUMAR, EPPSTEIN)

slide-46
SLIDE 46

Combinatorial Methods

WIN

slide-47
SLIDE 47

Combinatorial Methods

  • FIG. 6.

FAIL

slide-48
SLIDE 48

Fundamental requirements: Sensitive to noise:

Combinatorial Methods

sample spacing ≤ O(curve separation)

slide-49
SLIDE 49

Active Contours/Snakes

Start with small circle Expand circle, stopping at edges. Try to maintain curve smoothness.

slide-50
SLIDE 50

Level Set Segmentation

Don’t study contour directly - study level sets of auxiliary function instead. E(x) is result of edge detectors. ∂tφ(x, t) = |∇φ(x, t)| 1 + αE(x)f(φ(x, t) − 1)/2 + 2∆φ(x, t) + regularization

slide-51
SLIDE 51

LEVEL SET SEGMENTATION

slide-52
SLIDE 52

LEVEL SET SEGMENTATION

slide-53
SLIDE 53

LEVEL SET SEGMENTATION

slide-54
SLIDE 54

LEVEL SET SEGMENTATION

slide-55
SLIDE 55

LEVEL SET SEGMENTATION

slide-56
SLIDE 56

LEVEL SET SEGMENTATION

slide-57
SLIDE 57

LEVEL SET SEGMENTATION

slide-58
SLIDE 58

LEVEL SET SEGMENTATION

slide-59
SLIDE 59

LEVEL SET SEGMENTATION

slide-60
SLIDE 60

LEVEL SET SEGMENTATION

slide-61
SLIDE 61

2 dimensions is not 1 dimension “done twice”

slide-62
SLIDE 62

Parameterize images

ρ(x) =  

M−1

  • j=0

ρj1γj(x)   + ρtex(x)

slide-63
SLIDE 63

Parametric Model

Segmentation problem: Reconstruction problem: Find: M, γj(t) Find: M, γj(t), ρj, ρtex(x)

slide-64
SLIDE 64

Edges are the singular support of the function: Singular support is set of points

Singular Support

∀λ > 0, sup |

k|≥kr

  • ei

k·xρ(x)χ((x − x0)λ)dx

  • = O(kr

−3/2)

  • x0 = γj(t)
slide-65
SLIDE 65

Singular support extends to wavefront in higher dimensions Wavefront is set of surfels

Wavefront Set

∀λ > 0, sup

r≥kr

  • eirk0·xρ(x)χ((x − x0)λ)dx
  • = O(kr

−3/2)

( x0, k0) = (γj(t), ±Nj(t))

slide-66
SLIDE 66

2D is not 1D squared

singular support ⊂ RN

wavefront ⊂ RN × (SN−1/{±1})

Pxwavefront = singular support

slide-67
SLIDE 67

2D is not 1D squared

R1 × R1 = R2 × (S1/{±1})

slide-68
SLIDE 68

Wavefront Detection

slide-69
SLIDE 69

Wavefront Detectors

What does the Fourier transform of an edge look like?

slide-70
SLIDE 70

Wavefront Detectors

Calculate with Green’s Theorem

  • 1γj(

k) =

Ωj

ei

k·xdx1dx2 = Ωj

∂x1F2( k, x) − ∂x2F1( k, x)dx1dx2 =

  • S1 F(

k, γj(t)) · dγj(t) dt dt = 1 i| k|2

  • S1 ei

k·γj(t)

k⊥ · γ

j(t)dt

(HAT TIP: EUGENE SORETS)

slide-71
SLIDE 71

Wavefront Detectors

Phase stationary when k · γ′(t) = 0

M−1

  • j=0

ρj 1γj( k) =

M−1

  • j=0

ρj   ei

k·γ(tj( k))

  • k
  • 3/2

√π

  • κj(tj(

k)) + ei

k·γ(tj(−k))

  • k
  • 3/2

√π

  • κj(tj(−

k))    + O(kr

5/2)

(2.2)

tj( k) satisfies k · γ′(tj( k)) = 0

slide-72
SLIDE 72

Wavefront Detectors

Ray encodes location of edges with normals pointing in direction Localizing on this region yields surfels in the wavefront pointing in direction

  • k = kr

kθ

  • kθ
  • kθ
slide-73
SLIDE 73

DIRECTIONAL FILTERS

slide-74
SLIDE 74

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-75
SLIDE 75

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-76
SLIDE 76

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-77
SLIDE 77

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-78
SLIDE 78

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-79
SLIDE 79

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-80
SLIDE 80

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-81
SLIDE 81

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-82
SLIDE 82

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-83
SLIDE 83

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-84
SLIDE 84

WAVEFRONT FILTERS

ARROWS ARE TANGENTIAL TO THE EDGE

slide-85
SLIDE 85

SKIN OF LOWER BACK

slide-86
SLIDE 86

ZERO CURVATURE EXAMPLE

SPURIOUS EDGES ARE NOT DETECTED

slide-87
SLIDE 87

ZERO CURVATURE EXAMPLE

SPURIOUS EDGES ARE NOT DETECTED

slide-88
SLIDE 88

NOISE SENSITIVITY

slide-89
SLIDE 89

NOISE SENSITIVITY

slide-90
SLIDE 90

Analysis

slide-91
SLIDE 91

Assumptions

MRI measures Fourier transform of density Image piecewise constant plus smooth part The image boundaries are smooth Curvature bounded above and below The boundaries are separated from each other Minimum edge contrast

NOT SATISFIED IN PRACTICE

slide-92
SLIDE 92

How it works

eik·γj(tj(

k))

|k|3/2 √π

  • κj(tj(

k)) + O(|k|−5/2) Start with asymptotic expansion

slide-93
SLIDE 93

How it works

Drop higher order terms and apply directional filter eik·γj(tj(

k))

|k|3/2 √π

  • κj(tj(

k)) V(kθ)|k|1/2W(|k|)

slide-94
SLIDE 94

How it works

Then inverse Fourier Transform α

−α

∞ eik·γj(tj(kθ))e−ik·x k3/2

r

√π

  • κj(tj(kθ))

V(kθk3/2

r

W(krdkrdkθ

slide-95
SLIDE 95

How it works

Change variables tj(α)

tj(−α)

∞ eik·(γj(t)−x) k3/2

r

√π

  • κj(t)

V(kθ(t)k3/2

r

W(krdkrdt

slide-96
SLIDE 96

How it works

And evaluate inner integral tj(α)

tj(−α)

eik·γj(t) πκj(t)V(kθ(t)k3/2

r

ˇ W(Nj(t) · [γj(t) − x])dkrdt

slide-97
SLIDE 97

Proof of Correctness

SCHEMATIC PLOT OF THE INTEGRAND

slide-98
SLIDE 98

Proof of Correctness

Fast decay in normal direction Polynomial decay in tangential direction Parabolic scaling:

k domain: width = O(

  • length)

x domain: width = O(length2)

slide-99
SLIDE 99

Theorem

A directional filter will extract at least one surfel near the point where the tangent of an edge equals the direction of the filter. It will not extract surfels far from the edge. The theorem only applies to unrealistic parameter

  • choices. Algorithm still works on phantoms,

however.

slide-100
SLIDE 100

Segmentation with Surfels

slide-101
SLIDE 101

Combinatorial Reconstruction

Goal: combinatorial reconstruction of curves from scattered surfels How can tangential information help?

slide-102
SLIDE 102

Combinatorial Reconstruction

Points can only be connected in tangential direction.

slide-103
SLIDE 103

Reconstruction Algorithm:

Connect all points close to each other, but not within forbidden region. Prune graph, connecting only nearest tangential neighbors within the graph. Result is polygon with same topology as original curve. Then smooth polygon.

slide-104
SLIDE 104

Reconstruction Algorithm

Proven to work. Proof is an exercise in elementary calculus. Can filter uncorrelated noise via geometric constraints.

CURVE RECONSTRUCTION FROM POINTS AND TANGENTS.

  • L. GREENGARD AND C. STUCCHIO

ARXIV.ORG/ABS/0903.1817

sample spacing = O(

  • curve separation)
slide-105
SLIDE 105

FILTERING UNCORRELATED NOISE

slide-106
SLIDE 106

SEGMENTED PHANTOM

OVERSIMPLIFIED GEOMETRY

slide-107
SLIDE 107

Surfel Segmentation

Can prove segmentation algorithm correct by plugging output of wavefront theorem into input of curve reconstruction theorem. Combinatorial curve reconstruction only works for simplified geometry.

slide-108
SLIDE 108

Open problems

Build a level-set based segmentation algorithm that uses surfel data. Clean up the surfel data (Bayesian tricks)

slide-109
SLIDE 109

Reconstruction

slide-110
SLIDE 110

Reconstruction

Assume segmentation problem is approximately solved. Obvious idea: compute Fourier transform of discontinuities, subtract off, leaving only smooth part of function. Then manually draw discontinuities back.

slide-111
SLIDE 111

Fail

slide-112
SLIDE 112

Fourier Extension

Best approximation to low frequency data: High frequency data missing, but we can approximate:

ˆ ρmeas(k)

M−1

  • j=1

ρj 1γj(k) =

M−1

  • j=1

ρj 1 i|k|2

  • S1 eik·γj(t)k⊥ · γ′

j(t)dt

slide-113
SLIDE 113

Fourier Extension

Smooth Transition between them to avoid artifacts: ˆ ρreconstructed(k) = LPF(k)ˆ ρmeas(k) + HPF(k)

M−1

  • j=1

ρj 1γj(k)

slide-114
SLIDE 114

Fourier Extension

slide-115
SLIDE 115

Fourier Extension

slide-116
SLIDE 116

Conclusion

The wavefront of an image has more information than it’s singular support Surfels can be extracted directly from raw data Effectively segments and reconstructs phantoms Still needed: good geometric algorithms for surfel reconstruction