pynx.cdi: Coherent Diffraction Imaging

Description

This modules provides algorithm for 2D and 3D reconstruction of single objects using several algorithms:

  • hybrid input-output (HIO)
  • error reduction (ER)
  • relaxed averaged alternating reflections (RAAR),
  • maximum likelihood conjugate gradient minimization…
  • partial coherence

The calculations use an ‘operator’ approach, where operations on a given cdi dataset can be simply written by multiplying it by the corresponding operator (FT, projection,..) or by a series of operators.

API Reference

Note that the Python API uses an ‘operator’ approach, to enable writing complex operations in a more mathematical and natural way.

CDI base classes

This is the CDI base classes, which can be used with operators

class pynx.cdi.cdi.CDI(iobs, support=None, obj=None, mask=None, pixel_size_detector=None, wavelength=None, detector_distance=None)

Reconstruction class for two or three-dimensional coherent diffraction imaging data.

Constructor. All arrays are assumed to be centered in (0,0) to avoid fftshift

Parameters:
  • iobs

    2D/3D observed diffraction data (intensity). Assumed to be corrected and following Poisson statistics, will be converted to float32. Dimensions should be divisible by 4 and have a prime factor decomposition up to 7. The following special values are allowed:

    • values<=-1e30 are considered masked
    • -1e30 < values < 0 are observed but used as free pixels
  • support – initial support in real space (1 = inside support, 0 = outside)
  • obj – initial object. If None, will be initialized to a complex object, with uniform random 0<modulus<1 and 0<phase<2pi. The data is assumed to be centered in (0,0) to avoid fft shifts.
  • mask – mask for the diffraction data (0: valid pixel, >0: masked)
  • pixel_size_detector – detector pixel size (meters)
  • wavelength – experiment wavelength (meters)
  • detector_distance – detector distance (meters)
copy(copy_history_llk=False, copy_free_pixels=False)

Creates a copy (without any reference passing) of this object, unless copy_obj is False.

Parameters:
  • copy_history_llk – if True, the history, cycle number, llk, etc.. are copied to the new object.
  • copy_free_pixels – if True, copy the distribution of free pixels to the new CDI object. Otherwise the new CDI object does not have free pixels initialised.
Returns:

a copy of the object.

get_iobs(shift=False)

Get the observed intensity data array.

Parameters:shift – if True, the data array will be fft-shifted so that the center of the data is in the center of the array, rather than in the corner (the default).
Returns:the 2D or 3D CDI numpy data array
get_llk(noise=None, normalized=True)

Get the normalized log-likelihoods, which should converge to 1 for a statistically ideal fit.

Parameters:
  • noise – either ‘gaussian’, ‘poisson’ or ‘euclidian’, will return the corresponding log-likelihood.
  • normalized – if True, will return normalized values so that the llk from a statistically ideal model should converge to 1
Returns:

the log-likelihood, or if noise=None, a tuple of the three (poisson, gaussian, euclidian) log-likelihoods.

get_llkn()

Get the poisson normalized log-likelihood, which should converge to 1 for a statistically ideal fit to a Poisson-noise dataset.

Returns:the normalized log-likelihood, for Poisson noise.
get_obj(shift=False)

Get the object data array. This will automatically get the latest data, either from GPU or from the host memory, depending where the last changes were made.

Parameters:shift – if True, the data array will be fft-shifted so that the center of the data is in the center of the array, rather than in the corner (the default).
Returns:the 2D or 3D CDI numpy data array
get_support(shift=False)

Get the support array. This will automatically get the latest data, either from GPU or from the host memory, depending where the last changes were made.

Parameters:shift – if True, the data array will be fft-shifted so that the center of the data is in the center of the array, rather than in the corner (the default).
Returns:the 2D or 3D CDI numpy data array
get_x_y_z()

Get 1D arrays of x and y (z if 3d) coordinates, taking into account the pixel size. The arrays are centered at (0,0) - i.e. with the origin in the corner for FFT puroposes. x is an horizontal vector and y vertical, and (if 3d) z along third dimension.

Returns:a tuple (x, y) or (x, y, z) of 1D numpy arrays
in_object_space()
Returns:True if the current obj array is in object space, False otherwise.
init_free_pixels(ratio=0.05, island_radius=3, exclude_zone_center=0.05, verbose=False)

Random selection of ‘free’ pixels in the observed intensities, to be used as unbiased indicator. :param ratio: (approximate) relative number of pixels to be included in the free mask :param island_radius: free island radius, to avoid pixel correlation due to finit object size :param exclude_zone_center: the relative radius of the zone to be excluded near the center :param verbose: if True, be verbose :return: nothing. Free pixel values are modified as iobs_free = -iobs - 1

reset_history()

Reset history, and set current cycle to zero :return: nothing

save_data_cxi(filename, sample_name=None, experiment_id=None, instrument=None, process_parameters=None)

Save the diffraction data (observed intensity, mask) to an HDF% CXI file. :param filename: the file name to save the data to :param sample_name: optional, the sample name :param experiment_id: the string identifying the experiment, e.g.: ‘HC1234: Siemens star calibration tests’ :param instrument: the string identifying the instrument, e.g.: ‘ESRF id10’ :param process_parameters: a dictionary of parameters which will be saved as a NXcollection :return: Nothing. a CXI file is created

save_obj_cxi(filename, sample_name=None, experiment_id=None, instrument=None, note=None, crop=0, save_psf=False, process_notes=None, process_parameters=None, append=False, **kwargs)

Save the result of the optimisation (object, support) to an HDF5 CXI file. :param filename: the file name to save the data to :param sample_name: optional, the sample name :param experiment_id: the string identifying the experiment, e.g.: ‘HC1234: Siemens star calibration tests’ :param instrument: the string identifying the instrument, e.g.: ‘ESRF id10’ :param note: a string with a text note giving some additional information about the data, a publication… :param crop: integer, if >0, the object will be cropped arount the support, plus a margin of ‘crop’ pixels. :param save_psf: if True, also save the psf (if present) :param process_notes: a dictionary which will be saved in ‘/entry_N/process_1’. A dictionary entry can also be a ‘note’ as keyword and a dictionary as value - all key/values will then be saved as separate notes. Example: process={‘program’: ‘PyNX’, ‘note’:{‘llk’:1.056, ‘nb_photons’: 1e8}} :param process_parameters: a dictionary of parameters which will be saved as a NXcollection :param append: if True and the file already exists, the result will be saved in a new entry. :return: Nothing. a CXI file is created

set_free_pixel_mask(m, verbose=False, shift=False)

Set the free pixel mask. Assumes the free pixel mask is correct (excluding center, bad pixels) :param m: the boolean mask (1=masked, 0 otherwise) :param shift: if True, the data array will be fft-shifted so that the center of the stored data is

in the corner of the array. [default: the array is already shifted]
Returns:
set_iobs(iobs, shift=False)

Set the observed intensity data array.

Parameters:
  • iobs – the 2D or 3D CDI numpy data array (float32 numpy array)
  • shift – if True, the data array will be fft-shifted so that the center of the stored data is in the corner of the array. [default: the array is already shifted]
Returns:

nothing

set_mask(mask, shift=False)

Set the mask data array. Note that since the mask is stored by setting observed intensitied of masked pixels to negative values, it is not possible to un-mask pixels.

Parameters:
  • obj – the 2D or 3D CDI mask array
  • shift – if True, the data array will be fft-shifted so that the center of the stored data is in the corner of the array. [default: the array is already shifted]
Returns:

nothing

set_obj(obj, shift=False)

Set the object data array. Assumed to be in object (not Fourier) space

Parameters:
  • obj – the 2D or 3D CDI numpy data array (complex64 numpy array)
  • shift – if True, the data array will be fft-shifted so that the center of the stored data is in the corner of the array. [default: the array is already shifted]
Returns:

nothing

set_support(support, shift=False)

Set the support data array.

Parameters:
  • obj – the 2D or 3D CDI numpy data array (complex64 numpy array)
  • shift – if True, the data array will be fft-shifted so that the center of the stored data is in the corner of the array. [default: the array is already shifted]
Returns:

nothing

update_history(mode='llk', verbose=False, **kwargs)

Update the history record. :param mode: either ‘llk’ (will record new log-likelihood, nb_photons_calc, average value..), or ‘algorithm’ (will only update the algorithm) or ‘support’. :param verbose: if True, print some info about current process :param kwargs: other parameters to be recorded, e.g. support_threshold=xx, dt= :return: nothing

class pynx.cdi.cdi.OperatorCDI

Base class for an operator on CDI objects, not requiring a processing unit.

timestamp_increment(cdi)

Increment the timestamp counter corresponding to the processing language used (OpenCL or CUDA) Virtual, must be derived.

Parameters:w – the object (e.g. wavefront) the operator will be applied to.
Returns:
pynx.cdi.cdi.save_cdi_data_cxi(filename, iobs, wavelength=None, detector_distance=None, pixel_size_detector=None, mask=None, sample_name=None, experiment_id=None, instrument=None, note=None, iobs_is_fft_shifted=False, process_parameters=None)

Save the diffraction data (observed intensity, mask) to an HDF5 CXI file, NeXus-compatible. :param filename: the file name to save the data to :param iobs: the observed intensity :param wavelength: the wavelength of the experiment (in meters) :param detector_distance: the detector distance (in meters) :param pixel_size_detector: the pixel size of the detector (in meters) :param mask: the mask indicating valid (=0) and bad pixels (>0) :param sample_name: optional, the sample name :param experiment_id: the string identifying the experiment, e.g.: ‘HC1234: Siemens star calibration tests’ :param instrument: the string identifying the instrument, e.g.: ‘ESRF id10’ :param iobs_is_fft_shifted: if true, input iobs (and mask if any) have their origin in (0,0[,0]) and will be shifted back to centered-versions before being saved. :param process_parameters: a dictionary of parameters which will be saved as a NXcollection :return: Nothing. a CXI file is created

CDI Operators

This section lists the operators, which can be imported automatically using from pynx.cdi import *. When this import is done, either CUDA (preferably) or OpenCL operators will be imported. The documentation below corresponds to OpenCL operators, but this should be identical to CUDA operators.

CDI Runner class

The ‘runner’ class is used for command-line analysis scripts.

class pynx.cdi.runner.runner.CDIRunner(argv, params, cdi_runner_scan_class)

Class to process CDI data, with parameters from the command-line or from a text file.

Parameters:
  • argv – the command-line parameters
  • params – parameters for the optimization, with some default values.
  • ptycho_runner_scan_class – the class to use to run the analysis.
check_params()

Check if self.params includes a minimal set of valid parameters

Returns: Nothing. Will raise an exception if necessary

check_params_beamline()

Check if self.params includes a minimal set of valid parameters, specific to a beamline

Returns: Nothing. Will raise an exception if necessary

parse_arg(k, v=None)

Interprets one parameter. Will :param k: string with the name of the parameter :param v: value of the parameter, or None if a keyword parameter :return: True if parameter could be interpreted, False otherwise

parse_arg_beamline(k, v)

Parse argument in a beamline-specific way. This function only parses single arguments. If an argument is recognized and interpreted, the corresponding value is added to self.params

This method should be superseded in a beamline/instrument-specific child class.

Returns:True if the argument is interpreted, false otherwise
parse_argv()

Parses all the arguments given on a command line,

Returns: nothing

parse_parameters_file(filename)

Read parameters from a text file, written with one parameter per line as a python script :return: nothing. The parameters are accepted if understood, and stored in self.params

process_scans()

Run all the analysis on the supplied scan list

Returns:Nothing
exception pynx.cdi.runner.runner.CDIRunnerException
class pynx.cdi.runner.runner.CDIRunnerScan(params, scan)

Abstract class to handle CDI data. Must be derived to be used.

center_crop_data()

Auto-center diffraction pattern, and resize according to allowed FFT prime decomposition. Takes into account ROI if supplied. :return:

corr_flat_field()

Correct the image by a known flat field. Currently only applies correction for maxipix large pixels :return:

init_mask()

Load mask if the corresponding parameter has been set, or just initialize an array of 0. This must be called after iobs has been imported, and before cropping, so at the end of load_data()

Returns:Nothing. self.mask is created.
init_object()

Prepare the initial object. It will be random unless a file was supplied for input. :return: Nothing. self.object is created

init_support()

Prepare the initial support. Note that the mask should be initialized first. if self.params[‘support’] == ‘auto’, nothing is done, and the support will be created after the cdi object, to use GPU functions. :return: Nothing. self.support is created

load_data()

Loads data. This function only loads data (2D or 3D) from generic files, and should be derived for beamline-specific imports (from spec+individual images, etc..)

Returns: nothing. Loads data in self.iobs, and intializes mask and initial support.

prepare_cdi(free_pixel_mask=None)

Prepare CDI object from input parameters.

Parameters:free_pixel_mask – optionally supply a free pixel mask to keep it for successive runs
Returns:nothing. Creates self.cdi object.
prepare_data(init_mask=True, rebin_data=True, init_support=True, center_crop_data=True, corr_flat_field=True)

Prepare CDI data for processing.

Parameters:
  • init_mask – if True (the default), initialize the mask
  • rebin_data – if True (the default), rebin the data
  • init_support – if True (the default), initialize the support
  • center_crop_data – if True (the default), center & crop data
  • corr_flat_field – if True (the default), correct for the flat field
Returns:

prepare_processing_unit()

Prepare processing unit (CUDA, OpenCL, or CPU).

Returns: nothing. Creates self.processing_unit

rebin_data()

Will rebin the data (iobs and mask) and update pixel_size_detector according to self.params[‘rebin’]. this must be called at the end of load_data, after init_mask. If a ROI has been given, cropping will occur here. :return:

run(file_name=None)

Main work function. Will run selected algorithms according to parameters

Parameters:file_name – output file_name. If None, the result is not saved.
Returns:
run_algorithm(algo_string, file_name=None)

Run a single or suite of algorithms in a given run.

Parameters:
  • algo_string – a single or suite of algorithm steps to use, which can either correspond to a change of parameters, i.e. ‘beta=0.5’, ‘support_threshold=0.3’, ‘positivity=1’, or operators which should be applied to the cdi object, such as: ‘hio**20’: 20 cycles of HIO ‘er’: a single cycle or error reduction ‘detwinhio**20’: 20 cycles of hio after halving the support ‘er**20*hio**100’: 100 cycles of hio followed by 20 cycles of ER ‘(sup*er**10*raar**100)**10’: 10 times [100 cycles of RAAR + 10 ER + Support update]
  • file_name – output file_name. If None, the result is not saved.
save_result(file_name=None, verbose=True)

Save the results (object, support,…) at the end of a run :param file_name: the filename to save the data to. If its extension is either npz or cxi, this will

override params[‘output_format’]. Otherwise, the extension will be replaced. Note that the full filename will be appended with the llk value.
Returns:

Examples

# -*- coding: utf-8 -*-

# PyNX - Python tools for Nano-structures Crystallography
#   (c) 2017-present : ESRF-European Synchrotron Radiation Facility
#       authors:
#         Vincent Favre-Nicolin, favre@esrf.fr

import numpy as np
from scipy.fftpack import ifftshift, fftshift, fft2

# NB: It is possible to use the PYNX_PU environment variable to choose the GPU or language,
# e.g. using PYNX_PU=opencl or PYNX_PU=cuda or PYNX_PU=cpu or PYNX_PU=Titan, etc..

from pynx.utils.pattern import siemens_star
from pynx.cdi import *

# Test on a simulated pattern (2D)
n = 512

# Siemens-Star object
obj0 = siemens_star(n, nb_rays=18, r_max=60, nb_rings=3)

# Start from a slightly loose disc support
x, y = np.meshgrid(np.arange(-n // 2, n // 2, dtype=np.float32), np.arange(-n // 2, n // 2, dtype=np.float32))
r = np.sqrt(x ** 2 + y ** 2)
support = r < 65

iobs = abs(ifftshift(fft2(fftshift(obj0.astype(np.complex64))))) ** 2
iobs = np.random.poisson(iobs * 1e10 / iobs.sum())
mask = np.zeros_like(iobs, dtype=np.int16)
if True:
    # Mask some values in the central beam
    print("Removing %6.3f%% intensity" % (iobs[255:257, 255:257].sum() / iobs.sum() * 100))
    iobs[255:257, 255:257] = 0
    mask[255:257, 255:257] = 1

cdi = CDI(fftshift(iobs), obj=None, support=fftshift(support), mask=fftshift(mask), wavelength=1e-10,
          pixel_size_detector=55e-10)

cdi.init_free_pixels()

# Initial scaling of the object [ only useful if there are masked pixels !]
cdi = ScaleObj(method='F') * cdi

show = 40

# Do 100 cycles of RAAR
cdi = RAAR(calc_llk=20) ** 100 * cdi

if show > 0:
    cdi = ShowCDI(fig_num=1) * cdi

sup = SupportUpdate(threshold_relative=0.3, smooth_width=(2.0, 0.5, 800), force_shrink=False)

cdi = (sup * ER(calc_llk=20, show_cdi=show, fig_num=1) ** 5 * RAAR(calc_llk=20, show_cdi=show, fig_num=1) ** 40) ** 5 * cdi
cdi = DetwinRAAR(nb_cycle=10) * cdi
cdi = (sup * ER(calc_llk=20, show_cdi=show, fig_num=1) ** 5 * RAAR(calc_llk=20, show_cdi=show, fig_num=1) ** 40) ** 20 * cdi

# Finish with ML or ER
# cdi = ML(reg_fac=1e-2, calc_llk=5, show_cdi=show, fig_num=1) ** 100 * cdi
cdi = ER(calc_llk=5, show_cdi=show, fig_num=1) ** 100 * cdi

cdi = ShowCDI(fig_num=1) * cdi
# -*- coding: utf-8 -*-

# PyNX - Python tools for Nano-structures Crystallography
#   (c) 2017-present : ESRF-European Synchrotron Radiation Facility
#       authors:
#         Vincent Favre-Nicolin, favre@esrf.fr

import numpy as np
from scipy.ndimage import gaussian_filter
from scipy.fftpack import fftshift
import fabio
from matplotlib import pyplot as plt

from pynx.cdi import *

# CDI example from an experimental data set from id10 (courtesy of Yuriy Chushkin)

iobs = np.flipud(fabio.open("data/logo5mu_20sec.edf").data)
support = np.flipud(fabio.open("data/mask_logo5mu_20sec.edf").data)
n = len(iobs)
x, y = np.meshgrid(np.arange(0, n, dtype=np.float32), np.arange(0, n, dtype=np.float32))

# Mask specific to this dataset (from beamstop, after symmetrization of observed data)
mask = np.logical_or(iobs != 0, np.logical_or(abs(x - 258) > 30, abs(y - 258) > 30))
mask *= np.logical_or(iobs != 0, np.logical_or(abs(x - 246) > 30, abs(y) < 495))
mask *= np.logical_or(iobs != 0, np.logical_or(abs(x - 266) > 30, abs(y) > 20))
mask *= np.logical_or(iobs != 0, np.logical_or(abs(x - 10) > 30, abs(y - 270) > 5))
mask *= np.logical_or(iobs != 0, np.logical_or(abs(x - 498) > 30, abs(y - 240) > 5))
mask *= np.logical_or(iobs != 0, np.logical_or(abs(x) > 30, abs(y) > 20))
mask *= np.logical_or(iobs != 0, np.logical_or(abs(x-510) > 30, abs(y-510) > 20))
mask = (mask == 0)  # 0: OK, 1: masked

plt.figure(1, figsize=(8, 8))

# ========================= Try first from the known support (easy !) =============================
cdi = CDI(fftshift(iobs), obj=None, support=fftshift(support), mask=fftshift(mask), wavelength=1e-10,
          pixel_size_detector=55e-6)

# Set free mask for unbiased likelihood estimation
cdi.init_free_pixels()

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

# Do 4 * (50 cycles of HIO + 20 of ER), displaying object after each group of cycle
cdi = (ShowCDI(fig_num=1) * ER(calc_llk=10) ** 20 * HIO(calc_llk=20) ** 50) ** 4 * cdi

# cdi = (ShowCDI(fig_num=1) * ER(calc_llk=10) ** 20 * CF(calc_llk=20) ** 50) ** 4 * cdi

# cdi = (ShowCDI(fig_num=1) * ER(calc_llk=10) ** 20 * RAAR(calc_llk=20) ** 50) ** 4 * cdi

print("\n======================================================================\n")
print("This was too easy - start again from a loose support !\n")

# ========================= Try again from a loose support ========================================
support = np.flipud(fabio.open("data/mask_logo5mu_20sec.edf").data)
support = gaussian_filter(support.astype(np.float32), 4) > 0.2
cdi = CDI(fftshift(iobs), obj=None, support=fftshift(support), mask=fftshift(mask), wavelength=1e-10,
          pixel_size_detector=55e-6)

# Set free mask for unbiased likelihood estimation
cdi.init_free_pixels()

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

# Display every N cycles
show = 100

# Do 200 cycles of HIO, displaying object every N cycle and log-likelihood every 20 cycle
cdi = HIO(calc_llk=20, show_cdi=show, fig_num=2) ** 200 * cdi

# Support update operator
sup = SupportUpdate(threshold_relative=0.25, smooth_width=(2, 0.5, 800), force_shrink=False)

if True:
    # Do 40 cycles of HIO, then 5 of ER, update support, repeat
    cdi = (sup * ER(calc_llk=20, show_cdi=show, fig_num=2) ** 5
           * HIO(calc_llk=20, show_cdi=show, fig_num=2) ** 40) ** 20 * cdi
else:
    # Do 40 cycles of HIO, update support, repeat
    cdi = (sup * HIO(calc_llk=20, show_cdi=show, fig_num=2) ** 40) ** 20 * cdi

# Finish with ML or ER
# cdi = ML(reg_fac=1e-2, calc_llk=20, show_cdi=show, fig_num=1) ** 100 * cdi
cdi = ER(calc_llk=20, show_cdi=show, fig_num=2) ** 100 * cdi