Projectional Engines

Projectional Engines#

Learning more about projectional engines/algorithms such as the popular difference map algorithm.

In PtyPy we have implemented a generalised version of the difference map algorithm with the following update for the exit wave $ψ$ of every view

ψ_{j + 1} = {α (I - O) + F [(1 + α) O - α I]} (ψ_{j})

where $O$ and $F$ are the overlap and Fourier update, respectively. $I$ is the idenity matrix and $α$ is a tuning parameter that blends pure alternating projections ( $α = 0$ ) with the original DM algorithm ( $α = 1$ ). By its very nature, $α = 1$ is more exploratory but can be a bit too “aggressive” in some cases, whereas $α = 0$ has a tendency to get stuck in local minima. In practice, $α = 0.95$ seems to be a good and robust compromise. In our moonflower example, we can adjust the alpha parameter accordingly

p.engines.engine00 = u.Param()
p.engines.engine00.name = "DM"
p.engines.engine00.numiter = 80
p.engines.engine00.alpha = 0.95

The Fourier update#

During the Fourier update $F [Ψ]$ , the model $Ψ$ (transformed into the Fourier domain) is being updated such that its intensity is replaced by the measured data $I$ while its phase remains unchanged

F [Ψ] = Ψ \cdot \frac{\sqrt{| I |} + ρ \cdot (| Ψ | - \sqrt{| I |})}{| Ψ |}

where $ρ$ is a renormalisation factor defined as

\begin{array}{r} ρ = {\begin{cases} 1 & if b \geq ⟨ (| Ψ | - \sqrt{| I |})^{2} ⟩ \\ \sqrt{b / ⟨ (| Ψ | - \sqrt{| I |})^{2} ⟩} & if b < ⟨ (| Ψ | - \sqrt{| I |})^{2} ⟩ \\ 0 & if b = None \end{cases} \end{array}

with $b$ being the so called power bound, which can be modified using the fourier_power_bound engine parameter. For Poisson-sampled data, the theoretical value for this power bound parameter is $0.25$

p.engines.engine00.fourier_power_bound = 0.25

but it often makes sense to lower or higher this power bound parameter based on the experimental conditions. If the power bound is too low, it might lead to over-fitting. If the power bound is too high, the Fourier update would no longer be enforced.

The overlap update#

During the Overlap update $O [ψ]$ , the exit wave model $ψ$ (back-transformed into the real domain) is being used to update object $O$ and probe $P$

O = \frac{\sum_{k} P^{*} \cdot ψ_{k}}{\sum_{k} | P |^{2}}

P = \frac{\sum_{k} O_{k}^{*} \cdot ψ_{k}}{\sum_{k} | O_{k} |^{2}}

where $\sum_{k}$ is the summation over all views and $^{*}$ denotes the complex conjugate. Based on these updates, the new model is being updated such that

O [ψ]_{k} = P \cdot O_{k}

for all views $k$ . Since object and probe update depend on each other, both are repeated several times to reach convergence. In PtyPy, termination of this inner loop is triggered by reaching overlap_max_iterations or if the root mean square difference between the current and previous estimate of the probe is smaller than overlap_converge_factor. Additionally, the behaviour of the overlap update can be controlled by delaying the first update of the probe using probe_update_first or by swapping the order of the object and probe update using update_object_first which is True by default.

p.engines.engine00.overlap_converge_factor = 0.05
p.engines.engine00.overlap_max_iterations = 10
p.engines.engine00.update_object_first = True
p.engines.engine00.probe_update_start = 2

Challenge

Modify different DM engine parameters and observe their impact on the outcome of the reconstruction.

import ptypy
import ptypy.utils as u

p = u.Param()
p.verbose_level = "interactive"
p.io = u.Param()
p.io.rfile = None
p.io.autosave = u.Param(active=False)
p.io.interaction = u.Param(active=False)

# Live-plotting
p.io.autoplot = u.Param()
p.io.autoplot.active=True
p.io.autoplot.threaded = False
p.io.autoplot.layout = "jupyter"
p.io.autoplot.interval = 10

p.scans = u.Param()
p.scans.MF = u.Param()
p.scans.MF.name = "Full"
p.scans.MF.data= u.Param()
p.scans.MF.data.name = "MoonFlowerScan"
p.scans.MF.data.shape = 128
p.scans.MF.data.num_frames = 200
p.scans.MF.data.save = None
p.scans.MF.data.density = 0.2
p.scans.MF.data.photons = 1e8
p.scans.MF.data.psf = 0.

# Define reconstruction engines
p.engines = u.Param()

# Difference Map (DM) engine
p.engines.engine00 = u.Param()
p.engines.engine00.name = "DM"
p.engines.engine00.numiter = 80

# General update
p.engines.engine00.alpha = 0.95

# Fourier update
p.engines.engine00.fourier_power_bound = 0.25

# Overlap update
p.engines.engine00.overlap_converge_factor = 0.05
p.engines.engine00.overlap_max_iterations = 10
p.engines.engine00.update_object_first = True
p.engines.engine00.probe_update_start = 2

P = ptypy.core.Ptycho(p,level=5)

Projectional Engines

Contents

Projectional Engines#

The Fourier update#

The overlap update#