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 \(\psi\) of every view

\[ \psi_{j+1} = \{\alpha (\mathbb{I} - \mathcal{O}) + \mathcal{F} [(1 + \alpha) \mathcal{O} - \alpha \mathbb{I}]\}(\psi_{j}) \]

where \(\mathcal{O}\) and \(\mathcal{F}\) are the overlap and Fourier update, respectively. \(\mathbb{I}\) is the idenity matrix and \(\alpha\) is a tuning parameter that blends pure alternating projections (\(\alpha=0\)) with the original DM algorithm (\(\alpha=1\)). By its very nature, \(\alpha=1\) is more exploratory but can be a bit too “aggressive” in some cases, whereas \(\alpha=0\) has a tendency to get stuck in local minima. In practice, \(\alpha=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 \(\mathcal{F}[\Psi]\), the model \(\Psi\) (transformed into the Fourier domain) is being updated such that its intensity is replaced by the measured data \(I\) while its phase remains unchanged

\[ \mathcal{F}[\Psi] = \Psi \cdot \frac{\sqrt{|I|} + \rho \cdot (|\Psi| - \sqrt{|I|})}{|\Psi|} \]

where \(\rho\) is a renormalisation factor defined as

\[\begin{split} \rho = \begin{cases} 1 & \text{if } b \geq \langle(|\Psi| - \sqrt{|I|})^2\rangle \\ \sqrt{b / \langle(|\Psi| - \sqrt{|I|})^2\rangle} & \text{if } b < \langle(|\Psi| - \sqrt{|I|})^2\rangle \\ 0 & \text{if } b = \text{None} \end{cases} \end{split}\]

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 \(\mathcal{O}[\psi]\), the exit wave model \(\psi\) (back-transformed into the real domain) is being used to update object \(O\) and probe \(P\)

\[ O = \frac{\sum_k P^* \cdot \psi_k}{\sum_k|P|^2} \]
\[ P = \frac{\sum_k O_k^* \cdot \psi_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

\[ \mathcal{O}[\psi]_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)
Contents