CDI example loading data from a CXI file#
[1]:
%matplotlib widget
import numpy as np
import h5py
from numpy.fft import fftshift
from scipy.ndimage import center_of_mass
import matplotlib.pyplot as plt
# This imports all necessary operators. GPU will be auto-selected
from pynx.cdi import *
from pynx.utils.math import smaller_primes
Filename#
Change to your CXI file
[2]:
filename = "/Users/vincent/dev/pynx/pynx/cdi/examples/alignment_S2280.cxi"
Extract data#
[3]:
cxi = h5py.File(filename, 'r')
if '/entry_1/instrument_1/source_1/energy' in cxi:
nrj = cxi['/entry_1/instrument_1/source_1/energy'][()] / 1.60218e-16
wavelength = 12.384 / nrj * 1e-10
print(" CXI input: Energy = %8.2fkeV" % nrj)
else:
wavelength = None
if '/entry_1/instrument_1/detector_1/distance' in cxi:
detector_distance = cxi['/entry_1/instrument_1/detector_1/distance'][()]
print(" CXI input: detector distance = %8.2fm" % detector_distance)
else:
detector_distance = None
if '/entry_1/instrument_1/detector_1/x_pixel_size' in cxi:
pixel_size_detector = cxi['/entry_1/instrument_1/detector_1/x_pixel_size'][()]
print(" CXI input: detector pixel size = %8.2fum" % (pixel_size_detector * 1e6))
else:
pixel_size_detector = None
print(" CXI input: loading iobs")
if 'entry_1/instrument_1/detector_1/data' in cxi:
iobs = cxi['entry_1/instrument_1/detector_1/data'][()].astype(np.float32)
else:
iobs = cxi['entry_1/data_1/data'].value.astype(np.float32)
# This is the detector mask
if 'entry_1/instrument_1/detector_1/mask' in cxi:
mask = cxi['entry_1/instrument_1/detector_1/mask'][()].astype(np.int8)
nb = mask.sum()
print(" CXI input: loading mask, with %d pixels masked (%6.3f%%)" % (nb, nb * 100 / mask.size))
else:
mask = None
CXI input: Energy = 0.00keV
CXI input: loading iobs
CXI input: loading mask, with 1056784 pixels masked ( 1.544%)
Centre & crop data#
You can skip this if your data is already centred and if the arrays dimensions already have a prime factor decomposition with factors not larger than 7.
You can also change max_size to a larger value.
[4]:
max_size = 256
if iobs.ndim == 3:
nz0, ny0, nx0 = iobs.shape
# Find center of mass
z0, y0, x0 = center_of_mass(iobs)
print("Center of mass at:", z0, y0, x0)
iz0, iy0, ix0 = int(round(z0)), int(round(y0)), int(round(x0))
# Max symmetrical box around center of mass
nx = 2 * min(ix0, nx0 - ix0)
ny = 2 * min(iy0, ny0 - iy0)
nz = 2 * min(iz0, nz0 - iz0)
if max_size is not None:
nx = min(nx, max_size)
ny = min(ny, max_size)
nz = min(nz, max_size)
# Crop data to fulfill FFT size requirements
nz1, ny1, nx1 = smaller_primes((nz, ny, nx), maxprime=7, required_dividers=(2,))
print("Centering & reshaping data: (%d, %d, %d) -> (%d, %d, %d)" % (nz0, ny0, nx0, nz1, ny1, nx1))
iobs = iobs[iz0 - nz1 // 2:iz0 + nz1 // 2, iy0 - ny1 // 2:iy0 + ny1 // 2,
ix0 - nx1 // 2:ix0 + nx1 // 2]
if mask is not None:
mask = mask[iz0 - nz1 // 2:iz0 + nz1 // 2, iy0 - ny1 // 2:iy0 + ny1 // 2,
ix0 - nx1 // 2:ix0 + nx1 // 2]
else:
ny0, nx0 = iobs.shape
# Find center of mass
y0, x0 = center_of_mass(iobs)
print("Center of mass at:", y0, x0)
iy0, ix0 = int(round(y0)), int(round(x0))
# Max symmetrical box around center of mass
nx = 2 * min(ix0, nx0 - ix0)
ny = 2 * min(iy0, ny0 - iy0)
if max_size is not None:
nx = min(nx, max_size)
ny = min(ny, max_size)
nz = min(nz, max_size)
# Crop data to fulfill FFT size requirements
ny1, nx1 = smaller_primes((ny, nx), maxprime=7, required_dividers=(2,))
print("Centering & reshaping data: (%d, %d) -> (%d, %d)" % (ny0, nx0, ny1, nx1))
iobs = iobs[iy0 - ny1 // 2:iy0 + ny1 // 2, ix0 - nx1 // 2:ix0 + nx1 // 2]
if mask is not None:
mask = mask[iy0 - ny1 // 2:iy0 + ny1 // 2, ix0 - nx1 // 2:ix0 + nx1 // 2]
Center of mass at: 128.33750872519897 228.5596927477386 293.1825317435623
Centering & reshaping data: (257, 516, 516) -> (256, 256, 256)
Create & initialise the CDI object#
[5]:
cdi = CDI(fftshift(iobs), obj=None, support=None, mask=fftshift(mask), wavelength=wavelength,
pixel_size_detector=pixel_size_detector)
cdi = InitFreePixels() * cdi
Init the object support & a random object#
For Bragg CDI (without a beamstop) you can usually rely on auto-correlation to get a first estimate of the object support. Calling ScaleObj() is only necessary if a mask is used.
There are several options for the object random initialisation, use InitObjRandom? for details.
[6]:
cdi = AutoCorrelationSupport(threshold=0.1) * cdi
cdi = ShowCDI() * ScaleObj() * InitObjRandom(src="support",amin=0.8,amax=1, phirange=0.5) * cdi
Solve the object using RAAR & ER#
The important parameters below which can affect the reconstruction are the support update ones. You can try: * SupportUpdate(threshold_relative=0.12, method='max', smooth_width=(2,0.5,600)) * SupportUpdate(threshold_relative=0.3, method='rms', smooth_width=(2,0.5,600))
both should work in this case, but depending on datasets one may work better than the other
See the other support update options by typing SupportUpdate?
You can make more complex combinations of algorithms simply by multiplying the operators (ER, HIO, RAAR, SupportUpdate, CF,…) as you see fit.
Note that the recipe below should work on a majority of trials, but will sometimes fail. You can restart from a random object.
[7]:
sup = SupportUpdate(threshold_relative=0.3, method='rms', smooth_width=(2,0.5,600))
cdi = (sup * RAAR(beta=0.9, calc_llk=50, show_cdi=50)**20)**40 * cdi
cdi = (sup * ER(calc_llk=50, show_cdi=50)**20)**10 * cdi
RAAR # 0 LLK= 7.142[free= 1.273](p), nb photons=2.188722e+06, support:nb=210057 ( 1.252%) <obj>= 3.23 max= 1177.73, out=58.633% dt/cycle=1.2612s
RAAR # 50 LLK= 0.899[free= 0.354](p), nb photons=2.920514e+07, support:nb=198200 ( 1.181%) <obj>= 12.14 max= 1315.76, out=7.712% dt/cycle=0.0412s
RAAR #100 LLK= 0.989[free= 0.365](p), nb photons=3.000815e+07, support:nb=185884 ( 1.108%) <obj>= 12.71 max= 1291.80, out=8.122% dt/cycle=0.0411s
RAAR #150 LLK= 0.961[free= 0.363](p), nb photons=2.981901e+07, support:nb=184143 ( 1.098%) <obj>= 12.73 max= 1290.42, out=7.697% dt/cycle=0.0413s
RAAR #200 LLK= 1.012[free= 0.361](p), nb photons=3.022515e+07, support:nb=178240 ( 1.062%) <obj>= 13.02 max= 1271.01, out=8.549% dt/cycle=0.0427s
RAAR #250 LLK= 0.991[free= 0.353](p), nb photons=2.991549e+07, support:nb=175248 ( 1.045%) <obj>= 13.07 max= 1273.02, out=8.022% dt/cycle=0.0418s
RAAR #300 LLK= 1.070[free= 0.346](p), nb photons=3.047046e+07, support:nb=176571 ( 1.052%) <obj>= 13.14 max= 1254.09, out=9.101% dt/cycle=0.0421s
RAAR #350 LLK= 1.031[free= 0.338](p), nb photons=3.012230e+07, support:nb=163710 ( 0.976%) <obj>= 13.56 max= 1246.80, out=8.382% dt/cycle=0.0420s
RAAR #400 LLK= 1.141[free= 0.346](p), nb photons=3.102152e+07, support:nb=159442 ( 0.950%) <obj>= 13.95 max= 1233.33, out=9.646% dt/cycle=0.0410s
RAAR #450 LLK= 1.057[free= 0.326](p), nb photons=3.040137e+07, support:nb=154886 ( 0.923%) <obj>= 14.01 max= 1233.62, out=8.254% dt/cycle=0.0420s
RAAR #500 LLK= 1.152[free= 0.346](p), nb photons=3.095714e+07, support:nb=146458 ( 0.873%) <obj>= 14.54 max= 1213.34, out=10.470% dt/cycle=0.0410s
RAAR #550 LLK= 0.858[free= 0.314](p), nb photons=2.905919e+07, support:nb=121711 ( 0.725%) <obj>= 15.45 max= 1226.84, out=7.969% dt/cycle=0.0409s
RAAR #600 LLK= 0.834[free= 0.313](p), nb photons=2.987586e+07, support:nb=110876 ( 0.661%) <obj>= 16.42 max= 1228.06, out=6.306% dt/cycle=0.0412s
RAAR #650 LLK= 0.786[free= 0.293](p), nb photons=2.944555e+07, support:nb=111226 ( 0.663%) <obj>= 16.27 max= 1234.48, out=5.803% dt/cycle=0.0405s
RAAR #700 LLK= 0.828[free= 0.297](p), nb photons=3.001148e+07, support:nb=111183 ( 0.663%) <obj>= 16.43 max= 1233.24, out=5.931% dt/cycle=0.0416s
RAAR #750 LLK= 0.783[free= 0.293](p), nb photons=2.944706e+07, support:nb=111043 ( 0.662%) <obj>= 16.28 max= 1235.17, out=5.761% dt/cycle=0.0406s
ER #800 LLK= 0.827[free= 0.302](p), nb photons=2.999764e+07, support:nb=111272 ( 0.663%) <obj>= 16.42 max= 1233.02, out=5.961% dt/cycle=0.0409s
ER #850 LLK= 0.194[free= 0.203](p), nb photons=2.498174e+07, support:nb= 99204 ( 0.591%) <obj>= 15.87 max= 1238.38, out=3.292% dt/cycle=0.0354s
ER #900 LLK= 0.193[free= 0.201](p), nb photons=2.498246e+07, support:nb= 97997 ( 0.584%) <obj>= 15.97 max= 1235.91, out=3.334% dt/cycle=0.0359s
ER #950 LLK= 0.193[free= 0.201](p), nb photons=2.498050e+07, support:nb= 97633 ( 0.582%) <obj>= 16.00 max= 1234.94, out=3.333% dt/cycle=0.0364s
Save result in a CXI file#
Type cdi.save_obj_cxi? to get the full documentation about the function
[8]:
# cdi.save_obj_cxi("result-Pt.cxi")
[ ]: