Navigation

  • index
  • modules |
  • PyNX 2023.1.1 documentation »
  • CDI example: ESRF logo

CDI example: ESRF logo¶

[1]:
%matplotlib notebook
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 -O 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]:
################## 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

plt.figure()
cdi = (ER(show_cdi=20,calc_llk=20, fig_num=-1) ** 20 * HIO(show_cdi=50, calc_llk=20, fig_num=-1) ** 50) ** 4 * cdi
 HIO #  0 LLK= 4613.295[free=  0.000](p), nb photons=1.101198e+09, support:nb=  7577 ( 2.890%) <obj>=    381.23 max=   1707.42, out=15.593% dt/cycle=4.4058s
 HIO # 20 LLK= 2413.191[free=  0.000](p), nb photons=1.890848e+09, support:nb=  7577 ( 2.890%) <obj>=    499.55 max=   1647.40, out=3.097% dt/cycle=0.0154s
 HIO # 40 LLK= 4172.237[free=  0.000](p), nb photons=2.146526e+09, support:nb=  7577 ( 2.890%) <obj>=    532.25 max=   1668.80, out=5.347% dt/cycle=0.0003s
  ER # 60 LLK=  67.490[free=  0.000](p), nb photons=1.656973e+09, support:nb=  7577 ( 2.890%) <obj>=    467.64 max=   1706.49, out=0.484% dt/cycle=0.0037s
 HIO # 80 LLK= 1414.964[free=  0.000](p), nb photons=1.763803e+09, support:nb=  7577 ( 2.890%) <obj>=    482.48 max=   1671.53, out=2.781% dt/cycle=0.0195s
 HIO #100 LLK= 3402.531[free=  0.000](p), nb photons=2.023631e+09, support:nb=  7577 ( 2.890%) <obj>=    516.79 max=   1657.91, out=5.011% dt/cycle=0.0004s
  ER #120 LLK= 4767.501[free=  0.000](p), nb photons=2.269720e+09, support:nb=  7577 ( 2.890%) <obj>=    547.32 max=   1702.35, out=5.459% dt/cycle=0.0188s
 HIO #140 LLK=  60.848[free=  0.000](p), nb photons=1.668316e+09, support:nb=  7577 ( 2.890%) <obj>=    469.24 max=   1715.61, out=0.446% dt/cycle=0.0185s
 HIO #160 LLK= 2509.240[free=  0.000](p), nb photons=1.901482e+09, support:nb=  7577 ( 2.890%) <obj>=    500.95 max=   1660.06, out=3.722% dt/cycle=0.0210s
 HIO #180 LLK= 4046.212[free=  0.000](p), nb photons=2.133695e+09, support:nb=  7577 ( 2.890%) <obj>=    530.66 max=   1672.04, out=4.996% dt/cycle=0.0003s
  ER #200 LLK=  70.683[free=  0.000](p), nb photons=1.672874e+09, support:nb=  7577 ( 2.890%) <obj>=    469.88 max=   1726.76, out=0.505% dt/cycle=0.0003s
 HIO #220 LLK= 1432.781[free=  0.000](p), nb photons=1.770090e+09, support:nb=  7577 ( 2.890%) <obj>=    483.34 max=   1673.32, out=2.757% dt/cycle=0.0181s
 HIO #240 LLK= 3402.303[free=  0.000](p), nb photons=2.033601e+09, support:nb=  7577 ( 2.890%) <obj>=    518.07 max=   1670.75, out=4.964% dt/cycle=0.0003s
  ER #260 LLK= 4718.998[free=  0.000](p), nb photons=2.252937e+09, support:nb=  7577 ( 2.890%) <obj>=    545.29 max=   1707.74, out=5.682% dt/cycle=0.0189s

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.

[6]:
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
plt.figure()

cdi = (sup * HIO(calc_llk=0, positivity=True, fig_num=-1, 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=-1, show_cdi=25) ** 25)**20 * cdi

cdi = (sup * ER(calc_llk=25, positivity=True, fig_num=-1, show_cdi=25) ** 25)**8 * cdi
 HIO #300 LLK= 6769.820[free=  0.000](p), nb photons=2.218663e+09, support:nb= 14504 ( 5.533%) <obj>=    391.11 max=   1737.76, out=26.942% dt/cycle=0.0055s
 HIO #325 LLK= 4520.894[free=  0.000](p), nb photons=2.846647e+09, support:nb= 15996 ( 6.102%) <obj>=    421.85 max=   2161.88, out=6.585% dt/cycle=0.0167s
 HIO #350 LLK= 4544.086[free=  0.000](p), nb photons=2.998468e+09, support:nb= 15066 ( 5.747%) <obj>=    446.12 max=   2172.74, out=5.588% dt/cycle=0.0148s
 HIO #375 LLK= 5101.502[free=  0.000](p), nb photons=3.107423e+09, support:nb= 14910 ( 5.688%) <obj>=    456.52 max=   2176.52, out=5.782% dt/cycle=0.0148s
 HIO #400 LLK= 5341.136[free=  0.000](p), nb photons=3.229872e+09, support:nb= 14401 ( 5.494%) <obj>=    473.58 max=   2185.01, out=5.644% dt/cycle=0.0154s
 HIO #425 LLK= 5785.357[free=  0.000](p), nb photons=3.312973e+09, support:nb= 14601 ( 5.570%) <obj>=    476.34 max=   2183.78, out=5.748% dt/cycle=0.0169s
 HIO #450 LLK= 5928.857[free=  0.000](p), nb photons=3.395379e+09, support:nb= 14513 ( 5.536%) <obj>=    483.69 max=   2189.40, out=5.679% dt/cycle=0.0149s
 HIO #475 LLK= 6103.035[free=  0.000](p), nb photons=3.423146e+09, support:nb= 14680 ( 5.600%) <obj>=    482.89 max=   2192.24, out=5.653% dt/cycle=0.0150s
 HIO #500 LLK= 6176.540[free=  0.000](p), nb photons=3.450363e+09, support:nb= 14728 ( 5.618%) <obj>=    484.02 max=   2186.54, out=5.780% dt/cycle=0.0151s
 HIO #525 LLK= 6281.780[free=  0.000](p), nb photons=3.448123e+09, support:nb= 15171 ( 5.787%) <obj>=    476.74 max=   2176.77, out=5.894% dt/cycle=0.0172s
 HIO #550 LLK= 6303.238[free=  0.000](p), nb photons=3.478170e+09, support:nb= 15783 ( 6.021%) <obj>=    469.44 max=   2170.16, out=5.832% dt/cycle=0.0152s
 HIO #575 LLK= 6242.092[free=  0.000](p), nb photons=3.440241e+09, support:nb= 16578 ( 6.324%) <obj>=    455.54 max=   2162.80, out=5.670% dt/cycle=0.0154s
 HIO #600 LLK= 6161.415[free=  0.000](p), nb photons=3.391452e+09, support:nb= 17433 ( 6.650%) <obj>=    441.07 max=   2147.30, out=5.708% dt/cycle=0.0150s
 HIO #625 LLK= 6106.698[free=  0.000](p), nb photons=3.337689e+09, support:nb= 18235 ( 6.956%) <obj>=    427.83 max=   2139.47, out=6.032% dt/cycle=0.0172s
 HIO #650 LLK= 5894.200[free=  0.000](p), nb photons=3.325478e+09, support:nb= 18196 ( 6.941%) <obj>=    427.50 max=   2139.68, out=5.538% dt/cycle=0.0151s
 HIO #675 LLK= 5842.846[free=  0.000](p), nb photons=3.344531e+09, support:nb= 17755 ( 6.773%) <obj>=    434.02 max=   2138.82, out=5.417% dt/cycle=0.0150s
 HIO #700 LLK= 5921.910[free=  0.000](p), nb photons=3.380962e+09, support:nb= 17348 ( 6.618%) <obj>=    441.46 max=   2137.23, out=5.312% dt/cycle=0.0152s
 HIO #725 LLK= 5897.153[free=  0.000](p), nb photons=3.398094e+09, support:nb= 16595 ( 6.330%) <obj>=    452.51 max=   2132.60, out=5.443% dt/cycle=0.0152s
 HIO #750 LLK= 5978.066[free=  0.000](p), nb photons=3.450331e+09, support:nb= 16158 ( 6.164%) <obj>=    462.10 max=   2120.66, out=5.573% dt/cycle=0.0171s
 HIO #775 LLK= 5877.459[free=  0.000](p), nb photons=3.473927e+09, support:nb= 15150 ( 5.779%) <obj>=    478.85 max=   2115.43, out=5.477% dt/cycle=0.0152s
  ER #800 LLK= 5974.333[free=  0.000](p), nb photons=3.527506e+09, support:nb= 14474 ( 5.521%) <obj>=    493.67 max=   2109.35, out=5.400% dt/cycle=0.0153s
  ER #825 LLK= 131.406[free=  0.000](p), nb photons=2.851373e+09, support:nb=  9254 ( 3.530%) <obj>=    555.09 max=   2114.66, out=1.035% dt/cycle=0.0150s
  ER #850 LLK= 157.721[free=  0.000](p), nb photons=2.828425e+09, support:nb=  8669 ( 3.307%) <obj>=    571.20 max=   2091.28, out=1.070% dt/cycle=0.0167s
  ER #875 LLK= 173.950[free=  0.000](p), nb photons=2.785787e+09, support:nb=  8436 ( 3.218%) <obj>=    574.65 max=   2072.98, out=0.955% dt/cycle=0.0146s
  ER #900 LLK= 179.700[free=  0.000](p), nb photons=2.759463e+09, support:nb=  8263 ( 3.152%) <obj>=    577.89 max=   2060.09, out=0.954% dt/cycle=0.0146s
  ER #925 LLK= 184.696[free=  0.000](p), nb photons=2.729027e+09, support:nb=  8148 ( 3.108%) <obj>=    578.73 max=   2048.45, out=0.940% dt/cycle=0.0144s
  ER #950 LLK= 188.447[free=  0.000](p), nb photons=2.707151e+09, support:nb=  8056 ( 3.073%) <obj>=    579.69 max=   2039.58, out=0.950% dt/cycle=0.0170s
  ER #975 LLK= 191.420[free=  0.000](p), nb photons=2.686692e+09, support:nb=  7994 ( 3.049%) <obj>=    579.73 max=   2032.32, out=0.945% dt/cycle=0.0146s
[ ]:

Table of Contents

  • CDI example: ESRF logo
    • Get ESRF logo diffraction data
    • Optimisation from the known support
    • Optimisation from a loose support

This Page

  • Show Source

Quick search

Navigation

  • index
  • modules |
  • PyNX 2023.1.1 documentation »
  • CDI example: ESRF logo
© Copyright European Synchrotron Radiation Facility, Grenoble, France. Created using Sphinx 7.2.6.