Skip to content

nabu.preproc

source package nabu.preproc

Classes

  • CCDFilter — Filtering applied on radios.

  • Log — Helper class to take -log(radios)

  • CTFPhaseRetrieval — This class implements the CTF formula of [1] in its regularised form which avoids the zeros of unregularised_filter_denominator (unreg_filter_denom is the so here named denominator/delta_beta of equation 8).

  • DistortionCorrection — A class for estimating and correcting image distortion.

  • DoubleFlatField — Init double flat field by summing a series of urls and considering the same subregion of them.

  • FlatFieldDataUrls — Initialize a flat-field normalization process with DataUrls.

  • PaganinPhaseRetrieval — Paganin Phase Retrieval for an infinitely distant point source. Formula (10) in [1].

  • VerticalShift — This class is used when a vertical translation (along the tomography rotation axis) occurred.

source class CCDFilter(radios_shape: tuple, kernel_size: int = 3, correction_type: str = 'median_clip', median_clip_thresh: float = 0.1, abs_diff=False, preserve_borders=False)

Filtering applied on radios.

Initialize a CCDCorrection instance.

Parameters

  • radios_shape : tuple — A tuple describing the shape of the radios stack, in the form (n_radios, n_z, n_x).

  • kernel_size : int — Size of the kernel for the median filter. Default is 3.

  • correction_type : str — Correction type for radios ("median_clip", "sigma_clip", ...)

  • median_clip_thresh : float, optional — Threshold for the median clipping method.

  • abs_diff : boolean — by default False: the correction is triggered when img - median > threshold. If equals True: correction is triggered for abs(img-media) > threshold

  • preserve borders : boolean — by default False: If equals True: the borders (width=1) are not modified.

Notes

A CCD correction is a process (usually filtering) taking place in the radios space. Available filters: - median_clip: if the value of the current pixel exceeds the median of adjacent pixels (a 3x3 neighborhood) more than a threshold, then this pixel value is set to the median value.

Methods

  • median_filter — Perform a median filtering on an image.

  • median_clip_mask — Compute a mask indicating whether a pixel is valid or not, according to the median-clip method.

  • median_clip_correction — Compute the median clip correction on one image.

  • dezinger_correction — Compute the median clip correction on a radios stack, and propagates the invalid pixels into vert and horiz directions.

source staticmethod CCDFilter.median_filter(img, kernel_size=3)

Perform a median filtering on an image.

source method CCDFilter.median_clip_mask(img, return_medians=False)

Compute a mask indicating whether a pixel is valid or not, according to the median-clip method.

Parameters

  • img : numpy.ndarray — Input image

  • return_medians : bool, optional — Whether to return the median values additionally to the mask.

source method CCDFilter.median_clip_correction(radio, output=None)

Compute the median clip correction on one image.

Parameters

  • radios : numpy.ndarray, optional — A radio image.

  • output : numpy.ndarray, optional — Output array

source method CCDFilter.dezinger_correction(radios, dark=None, nsigma=5, output=None)

Compute the median clip correction on a radios stack, and propagates the invalid pixels into vert and horiz directions.

Parameters

  • radios : numpy.ndarray — A radios stack.

  • dark : numpy.ndarray, optional — A dark image. Default is None. If not None, it is subtracted from the radios.

  • nsigma : float, optional — Number of standard deviations to use for the zinger detection. Default is 5.

  • output : numpy.ndarray, optional — Output array

Raises

  • ValueError

source class Log(radios_shape, clip_min=None, clip_max=None)

Helper class to take -log(radios)

Parameters

  • clip_min : float, optional — Before taking the logarithm, the values are clipped to this minimum.

  • clip_max : float, optional — Before taking the logarithm, the values are clipped to this maximum.

Methods

  • take_logarithm — Take the negative logarithm of a radios chunk. Processing is done in-place !

source method Log.take_logarithm(radios)

Take the negative logarithm of a radios chunk. Processing is done in-place !

Parameters

  • radios : array — Radios chunk.

source class CTFPhaseRetrieval(shape, geo_pars, delta_beta, padded_shape='auto', padding_mode='reflect', translation_vh=None, normalize_by_mean=False, lim1=1e-05, lim2=0.2, use_rfft=False, fftw_num_threads=None, fft_num_threads=None, logger=None)

This class implements the CTF formula of [1] in its regularised form which avoids the zeros of unregularised_filter_denominator (unreg_filter_denom is the so here named denominator/delta_beta of equation 8).

References

[1] B. Yu, L. Weber, A. Pacureanu, M. Langer, C. Olivier, P. Cloetens, and F. Peyrin, "Evaluation of phase retrieval approaches in magnified X-ray phase nano computerized tomography applied to bone tissue", Optics Express, Vol 26, No 9, 11110-11124 (2018)

Initialize a Contrast Transfer Function phase retrieval.

Parameters

  • geo_pars : GeoPars — the geometry description

  • delta_beta : float — the delta/beta ratio

  • padded_shape : str or tuple, optional — Padded image shape, in the form (num_rows, num_columns) i.e (vertical, horizontal). By default, it is twice the image shape.

  • padding_mode : str — Padding mode. It must be valid for the numpy.pad function

  • translation_vh : array, optional — Shift in the form (y, x). It is used to perform a translation of the image before applying the CTF filter.

  • normalize_by_mean : bool — Whether to divide the (padded) image with its mean before applying the CTF filter.

  • lim1 : float >0 — the regulariser strenght at low frequencies

  • lim2 : float >0 — the regulariser strenght at high frequencies

  • use_rfft : bool, optional — Whether to use real-to-complex (R2C) FFT instead of usual complex-to-complex (C2C).

  • fftw_num_threads : bool or None or int, optional — DEPRECATED - please use fft_num_threads instead.

  • fft_num_threads : bool or None or int, optional — Number of threads to use for FFT. If a number is provided: number of threads to use for FFT. You can pass a negative number to use N - fft_num_threads cores.

  • logger : optional — a logger object

Methods

source method CTFPhaseRetrieval.retrieve_phase(img, output=None)

Apply the CTF filter to retrieve the phase.

Parameters

  • img : np.ndarray — Projection image. It must have been already flat-fielded.

Returns

  • ph : numpy.ndarray — Phase image

source class DistortionCorrection(estimation_method='fft-correlation', estimation_kwargs=None, correction_method='interpn', correction_kwargs=None)

A class for estimating and correcting image distortion.

Initialize a DistortionCorrection object.

Parameters

  • estimation_method : str — Name of the method to use for estimating the distortion

  • estimation_kwargs : dict, optional — Named arguments to pass to the estimation method, in the form of a dictionary.

  • correction_method : str — Name of the method to use for correcting the distortion

  • correction_kwargs : dict, optional — Named arguments to pass to the correction method, in the form of a dictionary.

Methods

source method DistortionCorrection.estimate_distortion(image, reference_image)

source method DistortionCorrection.correct_distortion(image, coords)

source method DistortionCorrection.estimate_and_correct(image, reference_image)

source class DoubleFlatField(shape, result_url=None, sub_region=None, detector_corrector=None, input_is_mlog=True, output_is_mlog=False, average_is_on_log=False, sigma_filter=None, filter_mode='reflect', log_clip_min=None, log_clip_max=None)

Init double flat field by summing a series of urls and considering the same subregion of them.

Parameters

  • shape : tuple — Expected shape of radios chunk to process

  • result_url : url, optional — where the double-flatfield is stored after being computed, and possibly read (instead of re-computed) before processing the images.

  • sub_region : tuple, optional — If provided, this must be a tuple in the form (start_x, end_x, start_y, end_y). Each image will be cropped to this region. This is used to specify a chunk of files. Each of the parameters can be None, in this case the default start and end are taken in each dimension.

  • input_is_mlog : boolean, default True — the input is considred as minus logarithm of normalised radios

  • output_is_mlog : boolean, default True — the output is considred as minus logarithm of normalised radios

  • average_is_on_log : boolean, False — the minus logarithm of the data is averaged the clipping value that is applied prior to the logarithm

  • sigma_filter : optional — if given a high pass filter is applied by signal -gaussian_filter(signal,sigma,filter_mode)

  • filter_mode : optional, default 'reflect' — the padding scheme applied a the borders ( same as scipy.ndimage.filtrs.gaussian_filter)

Methods

source method DoubleFlatField.compute_double_flatfield(radios, recompute=False)

Read the radios and generate the "double flat field" by averaging and possibly other processing.

Parameters

  • radios : array — Input radios chunk.

  • recompute : bool, optional — Whether to recompute the double flatfield if already computed.

source method DoubleFlatField.get_double_flatfield(radios=None, compute=False)

Get the double flat field or a subregion of it.

Parameters

  • radios : array, optional — Input radios chunk

  • compute : bool, optional — Whether to compute the double flatfield anyway even if a dump file exists.

Raises

  • ValueError

source method DoubleFlatField.apply_double_flatfield(radios)

Apply the "double flatfield" filter on a chunk of radios. The processing is done in-place !

source class FlatFieldDataUrls(radios_shape: tuple, flats: dict, darks: dict, radios_indices=None, interpolation: str = 'linear', distortion_correction=None, nan_value=1.0, radios_srcurrent=None, flats_srcurrent=None, **chunk_reader_kwargs)

Bases : FlatField

Initialize a flat-field normalization process with DataUrls.

Parameters

  • radios_shape : tuple — A tuple describing the shape of the radios stack, in the form (n_radios, n_z, n_x).

  • flats : dict — Dictionary where the key is the flat index, and the value is a silx.io.DataUrl pointing to the flat.

  • darks : dict — Dictionary where the key is the dark index, and the value is a silx.io.DataUrl pointing to the dark.

  • radios_indices : array, optional — Array containing the radios indices. radios_indices[0] is the index of the first radio, and so on.

  • interpolation : str, optional — Interpolation method for flat-field. See below for more details.

  • distortion_correction : DistortionCorrection, optional — A DistortionCorrection object. If provided, it is used to correct flat distortions based on each radio.

  • nan_value : float, optional — Which float value is used to replace nan/inf after flat-field.

Other Parameters

The other named parameters are passed to ChunkReader(). Please read its documentation for more information.

Notes

Usually, when doing a scan, only one or a few darks/flats are acquired. However, the flat-field normalization has to be performed on each radio, although incoming beam can fluctuate between projections. The usual way to overcome this is to interpolate between flats. If interpolation="nearest", the first flat is used for the first radios subset, the second flat is used for the second radios subset, and so on. If interpolation="linear", the normalization is done as a linear function of the radio index.

source class PaganinPhaseRetrieval(shape, distance=0.5, energy=20, delta_beta=250.0, pixel_size=1e-06, padding='edge', use_rfft=True, use_R2C=None, fftw_num_threads=None, fft_num_threads=None)

Paganin Phase Retrieval for an infinitely distant point source. Formula (10) in [1].

Parameters

  • shape : int or tuple — Shape of each radio, in the format (num_rows, num_columns), i.e (size_vertical, size_horizontal). If an integer is provided, the shape is assumed to be square.

  • distance : float, optional — Propagation distance in meters.

  • energy : float, optional — Energy in keV.

  • delta_beta : float, optional — delta/beta ratio, where n = (1 - delta) + i*beta is the complex refractive index of the sample.

  • pixel_size : float or tuple, optional — Detector pixel size in meters. Default is 1e-6 (one micron) If a tuple is passed, the pixel size is set as (horizontal_size, vertical_size).

  • padding : str, optional — Padding method. Available are "zeros", "mean", "edge", "sym", "reflect". Default is "edge". Please refer to the "Padding" section below for more details.

  • use_rfft : bool, optional — Whether to use Real-to-Complex (R2C) transform instead of standard Complex-to-Complex transform, providing better performances

  • use_R2C : bool, optional — DEPRECATED, use use_rfft instead

  • fftw_num_threads : bool or None or int, optional — DEPRECATED - please use fft_num_threads

  • fft_num_threads : bool or None or int, optional — Number of threads for FFT. Default is to use all available threads. You can pass a negative number to use N - fft_num_threads cores.

Important

Mind the units! Distance and pixel size are in meters, and energy is in keV.

References

[1] D. Paganin Et Al, "Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object", Journal of Microscopy, Vol 206, Part 1, 2002

Notes

Padding methods

The phase retrieval is a convolution done in Fourier domain using FFT, so the Fourier transform size has to be at least twice the size of the original data. Mathematically, the data should be padded with zeros before being Fourier transformed. However, in practice, this can lead to artefacts at the edges (Gibbs effect) if the data does not go to zero at the edges. Apart from applying an apodization (Hamming, Blackman, etc), a common strategy to avoid these artefacts is to pad the data. In tomography reconstruction, this is usually done by replicating the last(s) value(s) of the edges ; but one can think of other methods:

  • "zeros": the data is simply padded with zeros.
  • "mean": the upper side of extended data is padded with the mean of the first row, the lower side with the mean of the last row, etc.
  • "edge": the data is padded by replicating the edges. This is the default mode.
  • "sym": the data is padded by mirroring the data with respect to its edges. See numpy.pad().
  • "reflect": the data is padded by reflecting the data with respect to its edges, including the edges. See numpy.pad().

Formulas

The radio is divided, in the Fourier domain, by the original "Paganin filter" [1].

.. math::

F = 1 + \frac{\delta}{\beta} \lambda D \pi |k|^2

where k is the wave vector.

Methods

source method PaganinPhaseRetrieval.compute_filter()

source method PaganinPhaseRetrieval.pad_with_values(data, top_val=0, bottom_val=0, left_val=0, right_val=0)

Pad the data into self.padded_data with values.

Parameters

  • data : numpy.ndarray — data (radio)

  • top_val : float or numpy.ndarray, optional — Value(s) to fill the top of the padded data with.

  • bottom_val : float or numpy.ndarray, optional — Value(s) to fill the bottom of the padded data with.

  • left_val : float or numpy.ndarray, optional — Value(s) to fill the left of the padded data with.

  • right_val : float or numpy.ndarray, optional — Value(s) to fill the right of the padded data with.

source method PaganinPhaseRetrieval.pad_data(data, padding_method=None)

Raises

  • ValueError

source method PaganinPhaseRetrieval.apply_filter(radio, padding_method=None, output=None)

source method PaganinPhaseRetrieval.lmicron_to_db(Lmicron)

Utility to convert the "Lmicron" parameter of PyHST to a value of delta/beta. Please see the doc of nabu.preproc.phase.lmicron_to_db()

source class VerticalShift(radios_shape, shifts)

This class is used when a vertical translation (along the tomography rotation axis) occurred.

These translations are meant "per projection" and can be due either to mechanical errors, or can be applied purposefully with known motor movements to smear rings artefacts.

The object is initialised with an array of shifts: one shift for each projection. A positive shifts means that the axis has moved in the positive Z direction. The interpolation is done taking for a pixel (y,x) the pixel found at (y+shft,x) in the recorded images.

The method apply_vertical_shifts performs the correctionson the radios.

Parameters

  • radios_shape : tuple — Shape of the radios chunk, in the form (n_radios, n_y, n_x)

  • shifts : sequence of floats — one shift for each projection

Notes

During the acquisition, there might be other translations, each of them orthogonal to the rotation axis. - A "horizontal" translation in the detector plane: this is handled directly in the Backprojection operation. - A translation along the beam direction: this one is of no concern for parallel-beam geometry

Methods

  • apply_vertical_shifts — Parameters

    radios: a sequence of np.array The input radios. If the optional parameter is not given, they are modified in-place iangles: a sequence of integers Must have the same lenght as radios. It contains the index at which the shift is found in self.shifts given by shifts argument in the initialisation of the object. output: a sequence of np.array, optional If given, it will be modified to contain the shifted radios. Must be of the same shape of radios.

source method VerticalShift.apply_vertical_shifts(radios, iangles=None, output=None)

Parameters

radios: a sequence of np.array The input radios. If the optional parameter is not given, they are modified in-place iangles: a sequence of integers Must have the same lenght as radios. It contains the index at which the shift is found in self.shifts given by shifts argument in the initialisation of the object. output: a sequence of np.array, optional If given, it will be modified to contain the shifted radios. Must be of the same shape of radios.