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
[ ]: