1 of 24

1

Bayesian phase retrieval for image reconstruction using fast Fourier transforms in Stan��

�Brian Ward, Bob Carpenter, and David Barmherzig�June 23, 2023

2 of 24

2

Source

(X-Ray)

Grid of detectors

Specimen

Reference

3 of 24

3

Ideal

Given a reference 𝑅, data 𝑌 = | 𝓕( 𝑋 + 𝑅 ) |2

where 𝓕 is an oversampled Fourier transform

Recover the source image 𝑋

​

In practice

Measurement error (low photon counts)

Missing data (beamstop)

4 of 24

4

Specimen

Reference

𝑋 + 𝑅

5 of 24

5

| 𝓕( 𝑋 + 𝑅 ) |2

​

6 of 24

6

Measurement model

​

Observed photon flux 𝑌̃ is distributed

​

𝑌̃ ~ Poisson( 𝑁𝑝𝑌 / 𝑌̅ )

​

where

𝑌̅ is the average value of 𝑌

𝑁𝑝 is the average photon flux per pixel

​

​

In practice, 𝑁𝑝 must be small (<10)

7 of 24

7

𝑁𝑝 = 1

8 of 24

8

9 of 24

9

Missing data

​

A beamstop prevents damage to sensors, occludes lowest frequencies

​

​

​

​

​

​

​

​

​

We use a beamstop of size 25x25

10 of 24

10

Final Data

11 of 24

11

Aside: How much data is missing

​

Take 𝓕( 𝑋 + 𝑅 ), occlude the highest frequencies, and invert*

12 of 24

12

13 of 24

13

14 of 24

14

MLE solution

15 of 24

15

Sampling

16 of 24

16

17 of 24

17

18 of 24

​

  • Barmherzig, D. A., & Sun, J. (2022). Towards practical holographic coherent diffraction imaging via maximum likelihood estimation. Opt. Express, 30(5), 6886–6906. doi:10.1364/OE.445015

​

18

19 of 24

19

Bonus Slides

20 of 24

20

Numerical comparison

Metric

MLE

Posterior Mean

Draw #100

RMSE

0.212

0.217

0.276

Structural Similarity (SSIM)

0.371

0.408

0.223

RMSE: Lower is better

SSIM: Higher is better

21 of 24

21

Type

complex_matrix

22 of 24

22

Speed

Model can be more vectorized than written on previous slides

​

The slowest part of this model is the FFT

​

  • Eigen’s ~10 times slower than MatLab (FFTW)

​

  • ~20% faster just by updating to Eigen 3.4

23 of 24

23

Varying 𝑁𝑝

24 of 24

24

Varying prior standard deviation