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