PyNX

  • Changelog
  • Installation
  • Tutorials
  • Scripts Reference
  • API Reference

Table of Contents

  • Changelog
  • Installation
  • Tutorials
  • Scripts Reference
  • API Reference
On this page
  • Get ESRF logo diffraction data
  • Optimisation from the known support
  • Optimisation from a loose support

CDI example: ESRF logo¶

[1]:
%matplotlib ipympl
import os
import numpy as np
from numpy.fft import fftshift
from scipy.ndimage import gaussian_filter
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm

# This imports all necessary operators. GPU will be auto-selected
from pynx.cdi import *

Get ESRF logo diffraction data¶

This dataset was recorded on at id10@ESRF, courtesy of Yuriy Chushkin

[2]:
if not os.path.exists('logo5mu_20sec.npz'):
    os.system('curl -OL https://zenodo.org/record/3451855/files/logo5mu_20sec.npz')

print(list(np.load('logo5mu_20sec.npz').keys()))
iobs = np.load('logo5mu_20sec.npz')['iobs']
mask = np.load('logo5mu_20sec.npz')['mask']
support = np.load('logo5mu_20sec.npz')['support']

plt.figure(1, figsize=(9,4))
plt.subplot(121)
plt.imshow(iobs, norm=LogNorm())
plt.colorbar()
plt.subplot(122)
plt.imshow(mask)
plt.tight_layout()
['mask', 'support', 'iobs']

Optimisation from the known support¶

First create the CDI object.

Set the initial object array using random values from the existing mask.

then optimise it using 4 series with 50 cycles of HIO and 20 of ER algorithms.

It converges very easily as the support is known.

[3]:
# Create a figure here so it is displayed before computin begins in the next cell
fig_num = plt.figure().number
[4]:
################## Try first from the known support (too easy !) ################################################
cdi = CDI(fftshift(iobs), obj=None, support=fftshift(support), mask=fftshift(mask), wavelength=1e-10,
          pixel_size_detector=55e-6)

# Initial object
cdi = InitObjRandom(src="support", amin=0, amax=1, phirange=0) * cdi

# Initial scaling, required by mask
cdi = ScaleObj(method='F') * cdi

cdi = (ER(show_cdi=20,calc_llk=20, fig_num=fig_num) ** 20 * HIO(show_cdi=50, calc_llk=20, fig_num=fig_num) ** 50) ** 4 * cdi
 HIO #  0 LLK= 4578.038[free=  0.000](p), nb photons=1.105158e+09, support:nb=  7577 ( 2.890%) <obj>=    381.91 max=   1707.01, out=15.829% dt/cycle=0.5506s
 HIO # 20 LLK= 2442.971[free=  0.000](p), nb photons=1.893612e+09, support:nb=  7577 ( 2.890%) <obj>=    499.92 max=   1646.07, out=3.084% dt/cycle=0.0006s
 HIO # 40 LLK= 4128.753[free=  0.000](p), nb photons=2.142394e+09, support:nb=  7577 ( 2.890%) <obj>=    531.74 max=   1673.40, out=5.265% dt/cycle=0.0007s
  ER # 60 LLK=  66.098[free=  0.000](p), nb photons=1.657782e+09, support:nb=  7577 ( 2.890%) <obj>=    467.75 max=   1709.32, out=0.477% dt/cycle=0.0083s
 HIO # 80 LLK= 1409.898[free=  0.000](p), nb photons=1.761791e+09, support:nb=  7577 ( 2.890%) <obj>=    482.20 max=   1669.11, out=2.770% dt/cycle=0.0078s
 HIO #100 LLK= 3419.759[free=  0.000](p), nb photons=2.022614e+09, support:nb=  7577 ( 2.890%) <obj>=    516.66 max=   1664.03, out=4.991% dt/cycle=0.0005s
  ER #120 LLK= 4762.980[free=  0.000](p), nb photons=2.267921e+09, support:nb=  7577 ( 2.890%) <obj>=    547.10 max=   1707.47, out=5.724% dt/cycle=0.0004s
 HIO #140 LLK=  62.978[free=  0.000](p), nb photons=1.671479e+09, support:nb=  7577 ( 2.890%) <obj>=    469.68 max=   1721.63, out=0.460% dt/cycle=0.0074s
 HIO #160 LLK= 2521.567[free=  0.000](p), nb photons=1.903893e+09, support:nb=  7577 ( 2.890%) <obj>=    501.27 max=   1661.22, out=3.604% dt/cycle=0.0004s
 HIO #180 LLK= 4045.329[free=  0.000](p), nb photons=2.139451e+09, support:nb=  7577 ( 2.890%) <obj>=    531.38 max=   1676.58, out=4.786% dt/cycle=0.0004s
  ER #200 LLK=  73.400[free=  0.000](p), nb photons=1.675225e+09, support:nb=  7577 ( 2.890%) <obj>=    470.21 max=   1729.36, out=0.521% dt/cycle=0.0004s
 HIO #220 LLK= 1441.742[free=  0.000](p), nb photons=1.774170e+09, support:nb=  7577 ( 2.890%) <obj>=    483.89 max=   1674.28, out=2.685% dt/cycle=0.0077s
 HIO #240 LLK= 3434.888[free=  0.000](p), nb photons=2.036689e+09, support:nb=  7577 ( 2.890%) <obj>=    518.46 max=   1669.36, out=4.936% dt/cycle=0.0007s
  ER #260 LLK= 4745.007[free=  0.000](p), nb photons=2.253117e+09, support:nb=  7577 ( 2.890%) <obj>=    545.31 max=   1708.72, out=5.922% dt/cycle=0.0007s

Optimisation from a loose support¶

This will use the SupportUpdate operator in addition to ER and HIO.

We use positivity to facilitate convergence towards the solution. This example is difficult to get a good convergence because the object is divided in many small parts, and it is easy to get a wrong support

We can either start from: * the real support expanded by a gaussian convolution. * from a large circle.

The second choice is more difficult as the initial support is symmetrical, and there is an ambiguity between equivalent twin solutions. For this we perform a few cycles with a halved-support to break the symmetry

The critical parameter here is the threshold value. It is either possible to give it relative the the rms value of the density inside the support, the average or the maximum value (in this case, the threshold has to be set lower than for rms or average methods). In this example rms with threshold=0.28 works well (about half executions give a good result when starting from a smoothed initial support (much fewer for a large circle)).

If the resulting object does not look good, just re-execute the cell.

[5]:
if True:
    tmp = np.arange(-len(iobs)/2, len(iobs)/2)
    y, x = np.meshgrid(tmp, tmp, indexing='ij')
    r = np.sqrt(x**2+y**2)
    sup = r<70
else:
    sup = gaussian_filter(support.copy().astype(np.float32), 6) > 0.1

if False:
    plt.figure(figsize=(5,3))
    plt.imshow(sup, origin='lower')
    plt.title('Initial support')

cdi = CDI(fftshift(iobs), obj=None, support=fftshift(sup), mask=fftshift(mask), wavelength=1e-10,
          pixel_size_detector=55e-6)
# cdi.init_free_pixels()  # We will not use free log-likelihood in this example

# Initial object
cdi = InitObjRandom(src="support", amin=0, amax=1, phirange=0) * cdi

# Initial scaling, required by mask
cdi = ScaleObj(method='F') * cdi

# Support update operator
sup = SupportUpdate(threshold_relative=0.28, method='rms', smooth_width=(2,0.5,600),
                    force_shrink=False, post_expand=None)
#sup = SupportUpdate(threshold_relative=0.07, method='max', smooth_width=(2,0.5,600),
#                    force_shrink=False, post_expand=None)

# Start with HIO, detwin with a halved support, then HIO and ER,
# with a support update every 25 cycles
cdi = (sup * HIO(calc_llk=0, positivity=True, fig_num=fig_num, show_cdi=25) ** 25)**4 * cdi

cdi = DetwinHIO(positivity=True, detwin_axis=0)**20 * cdi

cdi = (sup * HIO(calc_llk=25, positivity=True, fig_num=fig_num, show_cdi=25) ** 25)**20 * cdi

cdi = (sup * ER(calc_llk=25, positivity=True, fig_num=fig_num, show_cdi=25) ** 25)**8 * cdi
 HIO #300 LLK= 6922.383[free=  0.000](p), nb photons=2.125585e+09, support:nb= 15094 ( 5.758%) <obj>=    375.26 max=   1713.41, out=28.770% dt/cycle=0.0026s
 HIO #325 LLK= 4113.302[free=  0.000](p), nb photons=2.764542e+09, support:nb= 16079 ( 6.134%) <obj>=    414.65 max=   2190.33, out=5.876% dt/cycle=0.0091s
 HIO #350 LLK= 4362.989[free=  0.000](p), nb photons=2.904820e+09, support:nb= 15213 ( 5.803%) <obj>=    436.97 max=   2192.30, out=5.103% dt/cycle=0.0066s
 HIO #375 LLK= 4764.158[free=  0.000](p), nb photons=3.060104e+09, support:nb= 14473 ( 5.521%) <obj>=    459.82 max=   2208.43, out=5.194% dt/cycle=0.0068s
 HIO #400 LLK= 5152.735[free=  0.000](p), nb photons=3.188035e+09, support:nb= 13875 ( 5.293%) <obj>=    479.34 max=   2206.85, out=5.367% dt/cycle=0.0068s
 HIO #425 LLK= 5523.166[free=  0.000](p), nb photons=3.248677e+09, support:nb= 13570 ( 5.177%) <obj>=    489.29 max=   2179.74, out=5.525% dt/cycle=0.0068s
 HIO #450 LLK= 5608.852[free=  0.000](p), nb photons=3.306351e+09, support:nb= 12955 ( 4.942%) <obj>=    505.19 max=   2160.48, out=5.445% dt/cycle=0.0085s
 HIO #475 LLK= 5864.203[free=  0.000](p), nb photons=3.370520e+09, support:nb= 12763 ( 4.869%) <obj>=    513.89 max=   2144.07, out=5.818% dt/cycle=0.0067s
 HIO #500 LLK= 6011.030[free=  0.000](p), nb photons=3.424409e+09, support:nb= 12490 ( 4.765%) <obj>=    523.61 max=   2127.83, out=5.742% dt/cycle=0.0070s
 HIO #525 LLK= 6354.836[free=  0.000](p), nb photons=3.422021e+09, support:nb= 12916 ( 4.927%) <obj>=    514.73 max=   2095.74, out=6.003% dt/cycle=0.0070s
 HIO #550 LLK= 6186.282[free=  0.000](p), nb photons=3.445023e+09, support:nb= 12637 ( 4.821%) <obj>=    522.12 max=   2090.94, out=5.724% dt/cycle=0.0081s
 HIO #575 LLK= 6410.061[free=  0.000](p), nb photons=3.491656e+09, support:nb= 13227 ( 5.046%) <obj>=    513.79 max=   2092.56, out=5.775% dt/cycle=0.0085s
 HIO #600 LLK= 6315.852[free=  0.000](p), nb photons=3.513132e+09, support:nb= 13609 ( 5.191%) <obj>=    508.08 max=   2091.54, out=5.680% dt/cycle=0.0068s
 HIO #625 LLK= 6267.222[free=  0.000](p), nb photons=3.492241e+09, support:nb= 13848 ( 5.283%) <obj>=    502.18 max=   2085.39, out=5.585% dt/cycle=0.0069s
 HIO #650 LLK= 6111.112[free=  0.000](p), nb photons=3.460939e+09, support:nb= 13461 ( 5.135%) <obj>=    507.06 max=   2078.57, out=5.382% dt/cycle=0.0068s
 HIO #675 LLK= 6397.715[free=  0.000](p), nb photons=3.443083e+09, support:nb= 14196 ( 5.415%) <obj>=    492.48 max=   2061.92, out=5.703% dt/cycle=0.0070s
 HIO #700 LLK= 6275.630[free=  0.000](p), nb photons=3.373595e+09, support:nb= 14846 ( 5.663%) <obj>=    476.70 max=   2051.03, out=5.627% dt/cycle=0.0085s
 HIO #725 LLK= 6010.485[free=  0.000](p), nb photons=3.391134e+09, support:nb= 14388 ( 5.489%) <obj>=    485.48 max=   2073.13, out=5.179% dt/cycle=0.0068s
 HIO #750 LLK= 5979.675[free=  0.000](p), nb photons=3.383946e+09, support:nb= 13907 ( 5.305%) <obj>=    493.28 max=   2061.53, out=5.011% dt/cycle=0.0068s
 HIO #775 LLK= 5810.552[free=  0.000](p), nb photons=3.294052e+09, support:nb= 13165 ( 5.022%) <obj>=    500.21 max=   2009.77, out=4.972% dt/cycle=0.0068s
  ER #800 LLK= 5043.545[free=  0.000](p), nb photons=2.781448e+09, support:nb= 11484 ( 4.381%) <obj>=    492.14 max=   1808.12, out=4.434% dt/cycle=0.0084s
  ER #825 LLK=  68.021[free=  0.000](p), nb photons=1.991054e+09, support:nb=  6404 ( 2.443%) <obj>=    557.59 max=   1753.39, out=0.807% dt/cycle=0.0066s
  ER #850 LLK=  74.080[free=  0.000](p), nb photons=1.863145e+09, support:nb=  5354 ( 2.042%) <obj>=    589.91 max=   1680.19, out=1.133% dt/cycle=0.0084s
  ER #875 LLK=  68.358[free=  0.000](p), nb photons=1.698009e+09, support:nb=  4498 ( 1.716%) <obj>=    614.41 max=   1601.49, out=0.881% dt/cycle=0.0064s
  ER #900 LLK=  67.290[free=  0.000](p), nb photons=1.643536e+09, support:nb=  4328 ( 1.651%) <obj>=    616.23 max=   1583.36, out=0.599% dt/cycle=0.0064s
  ER #925 LLK=  69.028[free=  0.000](p), nb photons=1.634482e+09, support:nb=  4300 ( 1.640%) <obj>=    616.53 max=   1580.82, out=0.557% dt/cycle=0.0064s
  ER #950 LLK=  69.461[free=  0.000](p), nb photons=1.632347e+09, support:nb=  4291 ( 1.637%) <obj>=    616.78 max=   1580.13, out=0.549% dt/cycle=0.0063s
  ER #975 LLK=  69.641[free=  0.000](p), nb photons=1.631749e+09, support:nb=  4290 ( 1.637%) <obj>=    616.73 max=   1579.98, out=0.545% dt/cycle=0.0079s
[ ]:

© Copyright European Synchrotron Radiation Facility, Grenoble, France.

Created using Sphinx 5.3.0.