abel package

abel.transform module

class abel.transform.Transform(IM, direction='inverse', method='three_point', origin='none', symmetry_axis=None, use_quadrants=(True, True, True, True), symmetrize_method='average', angular_integration=False, transform_options={}, center_options={}, angular_integration_options={}, recast_as_float64=True, verbose=False, center=<deprecated>)[source]

Bases: object

Abel transform image class.

This class provides whole-image forward and inverse Abel transforms, together with preprocessing (centering, symmetrizing) and postprocessing (integration) functions.

Parameters:
  • IM (a N×M numpy array) – This is the image to be transformed

  • direction (str) – The type of Abel transform to be performed.

    forward

    A forward Abel transform takes a (2D) slice of a 3D image and returns the 2D projection.

    inverse (default)

    An inverse Abel transform takes a 2D projection and reconstructs a 2D slice of the 3D image.

  • method (str) – specifies which numerical approximation to the Abel transform should be employed (see below). The options are

    basex

    the Gaussian “basis set expansion” method of Dribinski et al. (2002).

    direct

    a naive implementation of the analytical formula by Roman Yurchuk.

    hansenlaw

    the recursive algorithm described by Hansen and Law (1985).

    linbasex

    the 1D projections of velocity-mapping images in terms of 1D spherical functions by Gerber et al. (2013).

    onion_bordas

    the algorithm of Bordas and co-workers (1996), re-implemented by Rallis, Wells and co-workers (2014).

    onion_peeling

    the onion peeling deconvolution as described by Dasch (1992).

    rbasex

    a method similar to pBasex by Garcia et al. (2004) for velocity-mapping images, but with analytical basis functions developed by Ryazanov (2012).

    three_point

    the three-point transform of Dasch (1992).

    two_point

    the two-point transform of Dasch (1992).

  • origin (tuple or str) – Before applying Abel transform, the image is centered around this point.

    If a tuple (float, float) is provided, this specifies the image origin in the (row, column) format. If a string is provided, an automatic centering algorithm is used:

    image_center

    The origin is assumed to be the center of the image.

    convolution

    The origin is found from autoconvolution of image projections along each axis.

    slice

    The origin is found by comparing slices in the horizontal and vertical directions.

    com

    The origin is calculated as the center of mass.

    gaussian

    The origin is found using a fit to a Gaussian function. This only makes sense if your data looks like a Gaussian.

    none (default)

    No centering is performed. An image with an odd number of columns must be provided.

  • symmetry_axis (None, int or tuple) – Symmetrize the image about the numpy axis 0 (vertical), 1 (horizontal), (0, 1) (both axes). Note that the Abel transform is always performed around the vertical axis. This parameter only affects how the image is modified before (and after) applying the Abel transform. For more information, see the “Quadrant combining” note below.

  • use_quadrants (tuple of 4 booleans) – select quadrants to be used in the analysis: (Q0, Q1, Q2, Q3). Quadrants are numbered counter-clockwide from upper right. See note below for description of quadrants. Default is (True, True, True, True), which uses all quadrants.

  • symmetrize_method (str) – Method used for symmetrizing the image.

    average

    Average the quadrants, in accordance with the symmetry_axis.

    fourier

    Axial symmetry implies that the Fourier components of the 2D projection should be real. Removing the imaginary components in reciprocal space leaves a symmetric projection.

    K. R. Overstreet, P. Zabawa, J. Tallant, A. Schwettmann, J. P. Shaffer, “Multiple scattering and the density distribution of a Cs MOT”, Optics Express 13, 9672–9682 (2005).

  • angular_integration (bool) – Integrate the image over angle to give the radial (speed) intensity distribution.

  • transform_options (tuple) – Additional arguments passed to the individual transform functions. See the documentation for the individual transform method for options.

  • center_options (tuple) – Additional arguments to be passed to the centering function, see abel.tools.center.center_image().

  • angular_integration_options (tuple (or dict)) – Additional arguments passed to the angular integration functions, see abel.tools.vmi.angular_integration().

  • recast_as_float64 (bool) – determines whether the input image should be recast to float64. Many images are imported in other formats (such as uint8 or uint16), and this does not always play well with the transorm algorithms. This should probably always be set to True (default).

  • verbose (bool) – determines whether non-critical output should be printed.

Note

Quadrant combining: The quadrants can be combined (averaged) using the use_quadrants keyword in order to provide better data quality.

The quadrants are numbered starting from Q0 in the upper right and proceeding counter-clockwise:

+--------+--------+
| Q1   * | *   Q0 |
|   *    |    *   |
|  *     |     *  |                                 AQ1 | AQ0
+--------o--------+ --([inverse] Abel transform)--> ----o----
|  *     |     *  |                                 AQ2 | AQ3
|   *    |    *   |
| Q2  *  | *   Q3 |          AQi == [inverse] Abel transform
+--------+--------+                 of quadrant Qi

Three cases are possible:

  1. symmetry_axis = 0 (vertical):

    Combine:  Q01 = Q0 + Q1, Q23 = Q2 + Q3
    inverse image   AQ01 | AQ01
                    -----o----- (left and right sides equivalent)
                    AQ23 | AQ23
    
  2. symmetry_axis = 1 (horizontal):

    Combine: Q12 = Q1 + Q2, Q03 = Q0 + Q3
    inverse image   AQ12 | AQ03
                    -----o----- (top and bottom equivalent)
                    AQ12 | AQ03
    
  3. symmetry_axis = (0, 1) (both):

    Combine: Q = Q0 + Q1 + Q2 + Q3
    inverse image   AQ | AQ
                    ---o---  (all quadrants equivalent)
                    AQ | AQ
    

Notes

As mentioned above, PyAbel offers several different approximations to the the exact Abel transform. All the methods should produce similar results, but depending on the level and type of noise found in the image, certain methods may perform better than others. Please see the Transform Methods section of the documentation for complete information.

The methods marked with a * indicate methods that generate basis sets. The first time they are run for a new image size, it takes seconds to minutes to generate the basis set. However, this basis set is saved to disk can be reloaded, meaning that future transforms are performed much more quickly.

basex *

The “basis set exapansion” algorithm describes the data in terms of gaussian-like functions, which themselves can be Abel-transformed analytically. With the default functions, centered at each pixel, this method also does not make any assumption about the shape of the data. This method is one of the de-facto standards in photoelectron/photoion imaging.

V. Dribinski, A. Ossadtchi, V. A. Mandelshtam, H. Reisler, “Reconstruction of Abel-transformable images: The Gaussian basis-set expansion Abel transform method”, Rev. Sci. Instrum. 73, 2634–2642 (2002).

direct
This method attempts a direct integration of the Abel-transform integral. It makes no assumptions about the data (apart from cylindrical symmetry), but it typically requires fine sampling to converge. Such methods are typically inefficient, but thanks to this Cython implementation (by Roman Yurchuk), this “direct” method is competitive with the other methods.
hansenlaw

This “recursive algorithm” produces reliable results and is quite fast (~0.1 s for a 1001×1001 image). It makes no assumptions about the data (apart from cylindrical symmetry). It tends to require that the data is finely sampled for good convergence.

E. W. Hansen, P.-L. Law, “Recursive methods for computing the Abel transform and its inverse”, J. Opt. Soc. Am. A 2, 510–520 (1985).

linbasex *

Velocity-mapping images are composed of projected Newton spheres with a common centre. The 2D images are usually evaluated by a decomposition into base vectors, each representing the 2D projection of a set of particles starting from a centre with a specific velocity distribution. Lin-BASEX evaluates 1D projections of VM images in terms of 1D projections of spherical functions, instead.

Th. Gerber, Yu. Liu, G. Knopp, P. Hemberger, A. Bodi, P. Radi, Ya. Sych, “Charged particle velocity map image reconstruction with one-dimensional projections of spherical functions”, Rev. Sci. Instrum. 84, 033101 (2013).

onion_bordas

The onion peeling method, also known as “back projection”, originates from C. Bordas, F. Paulig, “Photoelectron imaging spectrometry: Principle and inversion method”, Rev. Sci. Instrum. 67, 2257–2268 (1996).

The algorithm was subsequently coded in MatLab by C. E. Rallis, T. G. Burwitz, P. R. Andrews, M. Zohrabi, R. Averin, S. De, B. Bergues, B. Jochim, A. V. Voznyuk, N. Gregerson, B. Gaire, I. Znakovskaya, J. McKenna, K. D. Carnes, M. F. Kling, I. Ben-Itzhak, E. Wells, “Incorporating real time velocity map image reconstruction into closed-loop coherent control”, Rev. Sci. Instrum. 85, 113105 (2014), which was used as the basis of this Python port. See issue #56.

onion_peeling *
This is one of the most compact and fast algorithms, with the inverse Abel transform achieved in one Python code-line, PR #155. See also three_point is the onion peeling algorithm as described by Dasch (1992), reference below.
rbasex *

The pBasex method by G. A. Garcia, L. Nahon, I. Powis, “Two-dimensional charged particle image inversion using a polar basis function expansion”, Rev. Sci. Instrum. 75, 4989–2996 (2004) adapts the BASEX (“basis set expansion”) method to the specific case of velocity-mapping images by using a basis of 2D functions in polar coordinates, such that the reconstructed radial distributions are obtained directly from the expansion coefficients.

This method employs the same approach, but uses more convenient basis functions, which have analytical Abel transforms separable into radial and angular parts, developed in M. Ryazanov, “Development and implementation of methods for sliced velocity map imaging. Studies of overtone-induced dissociation and isomerization dynamics of hydroxymethyl radical (CH2OH and CD2OH)”, Ph.D. dissertation, University of Southern California, 2012 (ProQuest, USC).

three_point *

The “Three Point” Abel transform method exploits the observation that the value of the Abel inverted data at any radial position r is primarily determined from changes in the projection data in the neighborhood of r. This method is also very efficient once it has generated the basis sets.

C. J. Dasch, “One-dimensional tomography: a comparison of Abel, onion-peeling, and filtered backprojection methods”, Appl. Opt. 31, 1146–1152 (1992).

two_point *
Another Dasch method. Simple, and fast, but not as accurate as the other methods.

The following class attributes are available, depending on the calculation.

Returns:
  • transform (numpy 2D array) – the 2D forward/inverse Abel-transformed image.
  • angular_integration (tuple) – (radial-grid, radial-intensity) radial coordinates and the radial intensity (speed) distribution, evaluated using abel.tools.vmi.angular_integration().
  • residual (numpy 2D array) – residual image (not currently implemented).
  • IM (numpy 2D array) – the input image, re-centered (optional) with an odd-size width.
  • method (str) – transform method, as specified by the input option.
  • direction (str) – transform direction, as specified by the input option.
  • Beta (numpy 2D array) – with method=linbasex, transform_options=dict(return_Beta=True): Beta array coefficients of Newton-sphere spherical harmonics
    Beta[0] - the radial intensity variation

    Beta[1] - the anisotropy parameter variation

    …Beta[n] - higher-order terms up to legedre_orders=[0, ..., n]

  • radial (numpy 1D array) – with method=linbasex, transform_options=dict(return_Beta=True): radial grid for Beta array
  • projection (numpy 2D array) – with method=linbasex, transform_options=dict(return_Beta=True): radial projection profiles at angles proj_angles
  • distr (Distributions.Results object) – with method=rbasex: the object from which various radial distributions can be retrieved

abel.basex module

abel.basex.basex_transform(data, sigma=1.0, reg=0.0, correction=True, basis_dir='./', dr=1.0, verbose=True, direction='inverse')[source]

This function performs the BASEX (BAsis Set EXpansion) Abel transform. It works on a “right side” image. I.e., it works on just half of a cylindrically symmetric object, and data[0,0] should correspond to a central pixel. To perform a BASEX transform on a whole image, use

abel.Transform(image, method='basex', direction='inverse').transform

This BASEX implementation only works with images that have an odd-integer full width.

Parameters:
  • data (m × n numpy array) – the image to be transformed. data[:,0] should correspond to the central column of the image.

  • sigma (float) – width parameter for basis functions, see equation (14) in the article. Determines the number of basis functions (n/sigma rounded). Can be any positive number, but using sigma < 1 is not very meaningful and requires regularization.

  • reg (float) –

    regularization parameter, square of the Tikhonov factor.

    reg=0 means no regularization,

    reg=100 is a reasonable value for megapixel images.

    Forward transform requires regularization only if sigma < 1, and reg should be ≪ 1.

  • correction (boolean) – apply intensity correction in order to reduce method artifacts (intensity normalization and oscillations)

  • basis_dir (str) – path to the directory for saving / loading the basis sets. If None, the basis set will not be saved to disk.

  • dr (float) – size of one pixel in the radial direction. This only affects the absolute scaling of the transformed image.

  • verbose (boolean) – determines whether statements should be printed

  • direction (str: 'forward' or 'inverse') – type of Abel transform to be performed

Returns:

recon – the transformed (half) image

Return type:

m × n numpy array

abel.basex.basex_core_transform(rawdata, A)[source]

Internal function that does the actual BASEX transform. It requires that the transform matrix be passed.

Parameters:
  • rawdata (m × n numpy array) – right half (with the axis) of the input image.
  • A (n × n numpy array) – 2D array given by the transform-calculation function
Returns:

IM – the Abel-transformed image

Return type:

m × n numpy array

abel.basex.get_bs_cached(n, sigma=1.0, reg=0.0, correction=True, basis_dir='.', dr=1.0, verbose=False, direction='inverse')[source]

Internal function.

Gets BASEX basis sets, using the disk as a cache (i.e. load from disk if they exist, if not, calculate them and save a copy on disk) and calculates the transform matrix. To prevent saving the basis sets to disk, set basis_dir=None. Loaded/calculated matrices are also cached in memory.

Parameters:
  • n (int) – Abel transform will be performed on an n pixels wide area of the (half) image
  • sigma (float) – width parameter for basis functions
  • reg (float) – regularization parameter
  • correction (boolean) – apply intensity correction. Corrects wrong intensity normalization (seen for narrow basis sets), intensity oscillations (seen for broad basis sets), and intensity drop-off near r = 0 due to regularization.
  • basis_dir (str) – path to the directory for saving / loading the basis sets. If None, the basis sets will not be saved to disk.
  • dr (float) – pixel size. This only affects the absolute scaling of the output.
  • verbose (boolean) – determines whether statements should be printed
  • direction (str: 'forward' or 'inverse') – type of Abel transform to be performed
Returns:

A – matrix of the Abel transform (forward or inverse)

Return type:

n × n numpy array

abel.basex.cache_cleanup(select='all')[source]

Utility function.

Frees the memory caches created by get_bs_cached(). This is usually pointless, but might be required after working with very large images, if more RAM is needed for further tasks.

Parameters:

select (str) – selects which caches to clean:

all (default)

everything, including basis;

forward

forward transform;

inverse

inverse transform.

Returns:

Return type:

None

abel.basex.get_basex_correction(A, sigma, direction)[source]

Internal function.

The default BASEX basis and the way its projection is calculated leads to artifacts in the reconstructed distribution – incorrect overall intensity for sigma = 1, intensity oscillations for other sigma values, intensity fluctuations (and drop-off for reg > 0) near r = 0. This function generates the intensity correction profile from the BASEX result for a step function with a soft edge (to avoid ringing) aligned with the last basis function.

Parameters:
  • A (n × n numpy array) – matrix of the Abel transform
  • sigma (float) – basis width parameter
  • direction (str: 'forward' or 'inverse') – type of the Abel transform
Returns:

cor – intensity correction profile

Return type:

1 × n numpy array

abel.linbasex module

abel.linbasex.linbasex_transform(IM, basis_dir=None, proj_angles=[0, 1.5707963267948966], legendre_orders=[0, 2], radial_step=1, smoothing=0, rcond=0.0005, threshold=0.2, return_Beta=False, clip=0, norm_range=(0, -1), direction='inverse', verbose=False, dr=None)[source]

Wrapper function for linbasex to process supplied quadrant-image as a full-image.

PyAbel transform functions operate on the right side of an image. Here we follow the basex technique of duplicating the right side to the left re-forming the whole image.

Inverse Abel transform using 1d projections of images.

Th. Gerber, Yu. Liu, G. Knopp, P. Hemberger, A. Bodi, P. Radi, Ya. Sych, “Charged particle velocity map image reconstruction with one-dimensional projections of spherical functions”, Rev. Sci. Instrum. 84, 033101 (2013).

linbasex models the image using a sum of Legendre polynomials at each radial pixel, As such, it should only be applied to situations that can be adequately represented by Legendre polynomials, i.e., images that feature spherical-like structures. The reconstructed 3D object is obtained by adding all the contributions, from which slices are derived.

Parameters:
  • IM (numpy 2D array) – image data must be square shape of odd size

  • proj_angles (list) – projection angles, in radians (default \([0, \pi/2]\)) e.g. \([0, \pi/2]\) or \([0, 0.955, \pi/2]\) or \([0, \pi/4, \pi/2, 3\pi/4]\)

  • legendre_orders (list) – orders of Legendre polynomials to be used as the expansion

    • even polynomials [0, 2, …] gerade
    • odd polynomials [1, 3, …] ungerade
    • all orders [0, 1, 2, …].

    In a single photon experiment there are only anisotropies up to second order. The interaction of 4 photons (four wave mixing) yields anisotropies up to order 8.

  • radial_step (int) – number of pixels per Newton sphere (default 1)

  • smoothing (float) – convolve Beta array with a Gaussian function of 1/e 1/2 width smoothing.

  • rcond (float) – (default 0.0005) scipy.linalg.lstsq fit conditioning value. set rcond to zero to switch conditioning off. Note: In the presence of noise the equation system may be ill posed. Increasing rcond smoothes the result, lowering it beyond a minimum renders the solution unstable. Tweak rcond to get a “reasonable” solution with acceptable resolution.

  • clip (int) – clip first vectors (smallest Newton spheres) to avoid singularities (default 0)

  • norm_range (tuple) – (low, high) normalization of Newton spheres, maximum in range Beta[0, low:high]. Note: Beta[0, i] the total number of counts integrated over sphere i, becomes 1.

  • threshold (float) – threshold for normalization of higher order Newton spheres (default 0.2) Set all Beta[j], j>=1 to zero if the associated Beta[0] is smaller than threshold.

  • return_Beta (bool) – return the Beta array of Newton spheres, as the tuple: radial-grid, Beta for the case legendre_orders=[0, 2]

    Beta[0] vs radius -> speed distribution

    Beta[2] vs radius -> anisotropy of each Newton sphere

    see ‘Returns’.

  • direction (str) – “inverse” - only option for this method. Abel transform direction.

  • dr (None) – dummy variable for call compatibility with the other methods

  • verbose (bool) – print information about processing (normally used for debugging)

Returns:

  • inv_IM (numpy 2D array) – inverse Abel transformed image

  • radial, Beta, projections (tuple) – (if return_Beta=True)

    contributions of each spherical harmonic \(Y_{i0}\) to the 3D distribution contain all the information one can get from an experiment. For the case legendre_orders=[0, 2]:

    Beta[0] vs radius -> speed distribution

    Beta[1] vs radius -> anisotropy of each Newton sphere.

    projections : are the radial projection profiles at angles proj_angles

abel.linbasex.linbasex_transform_full(IM, basis_dir=None, proj_angles=[0, 1.5707963267948966], legendre_orders=[0, 2], radial_step=1, smoothing=0, rcond=0.0005, threshold=0.2, clip=0, return_Beta=False, norm_range=(0, -1), direction='inverse', verbose=False)[source]

Inverse Abel transform using 1d projections of images.

Th. Gerber, Yu. Liu, G. Knopp, P. Hemberger, A. Bodi, P. Radi, Ya. Sych, “Charged particle velocity map image reconstruction with one-dimensional projections of spherical functions”, Rev. Sci. Instrum. 84, 033101 (2013).

linbasex models the image using a sum of Legendre polynomials at each radial pixel, As such, it should only be applied to situations that can be adequately represented by Legendre polynomials, i.e., images that feature spherical-like structures. The reconstructed 3D object is obtained by adding all the contributions, from which slices are derived.

Parameters:
  • IM (numpy 2D array) – image data must be square shape of odd size

  • proj_angles (list) – projection angles, in radians (default \([0, \pi/2]\)) e.g. \([0, \pi/2]\) or \([0, 0.955, \pi/2]\) or \([0, \pi/4, \pi/2, 3\pi/4]\)

  • legendre_orders (list) – orders of Legendre polynomials to be used as the expansion

    • even polynomials [0, 2, …] gerade
    • odd polynomials [1, 3, …] ungerade
    • all orders [0, 1, 2, …].

    In a single photon experiment there are only anisotropies up to second order. The interaction of 4 photons (four wave mixing) yields anisotropies up to order 8.

  • radial_step (int) – number of pixels per Newton sphere (default 1)

  • smoothing (float) – convolve Beta array with a Gaussian function of 1/e 1/2 width smoothing.

  • rcond (float) – (default 0.0005) scipy.linalg.lstsq fit conditioning value. set rcond to zero to switch conditioning off. Note: In the presence of noise the equation system may be ill posed. Increasing rcond smoothes the result, lowering it beyond a minimum renders the solution unstable. Tweak rcond to get a “reasonable” solution with acceptable resolution.

  • clip (int) – clip first vectors (smallest Newton spheres) to avoid singularities (default 0)

  • norm_range (tuple) – (low, high) normalization of Newton spheres, maximum in range Beta[0, low:high]. Note: Beta[0, i] the total number of counts integrated over sphere i, becomes 1.

  • threshold (float) – threshold for normalization of higher order Newton spheres (default 0.2) Set all Beta[j], j>=1 to zero if the associated Beta[0] is smaller than threshold.

  • return_Beta (bool) – return the Beta array of Newton spheres, as the tuple: radial-grid, Beta for the case legendre_orders=[0, 2]

    Beta[0] vs radius -> speed distribution

    Beta[2] vs radius -> anisotropy of each Newton sphere

    see ‘Returns’.

  • direction (str) – “inverse” - only option for this method. Abel transform direction.

  • dr (None) – dummy variable for call compatibility with the other methods

  • verbose (bool) – print information about processing (normally used for debugging)

Returns:

  • inv_IM (numpy 2D array) – inverse Abel transformed image

  • radial, Beta, projections (tuple) – (if return_Beta=True)

    contributions of each spherical harmonic \(Y_{i0}\) to the 3D distribution contain all the information one can get from an experiment. For the case legendre_orders=[0, 2]:

    Beta[0] vs radius -> speed distribution

    Beta[1] vs radius -> anisotropy of each Newton sphere.

    projections : are the radial projection profiles at angles proj_angles

abel.linbasex.int_beta(Beta, radial_step=1, threshold=0.1, regions=None)[source]

Integrate beta over a range of Newton spheres.

Parameters:
  • Beta (numpy array) – Newton spheres
  • radial_step (int) – number of pixels per Newton sphere (default 1)
  • threshold (float) – threshold for normalisation of higher orders, 0.0 … 1.0.
  • regions (list of tuple radial ranges) – [(min0, max0), (min1, max1), …]
Returns:

Beta_in – integrated normalized Beta array [Newton sphere, region]

Return type:

numpy array

abel.linbasex.get_bs_cached(cols, basis_dir=None, legendre_orders=[0, 2], proj_angles=[0, 1.5707963267948966], radial_step=1, clip=0, verbose=False)[source]

load basis set from disk, generate and store if not available.

Checks whether file: linbasex_basis_{cols}_{legendre_orders}_{proj_angles}_{radial_step}_{clip}*.npy is present in basis_dir

Either, read basis array or generate basis, saving it to the file.

Parameters:
  • cols (int) – width of image
  • basis_dir (str) – path to the directory for saving / loading the basis
  • legendre_orders (list) – default [0, 2] = 0 order and 2nd order polynomials
  • proj_angles (list) – default [0, np.pi/2] in radians
  • radial_step (int) – pixel grid size, default 1
  • clip (int) – image edge clipping, default 0 pixels
  • verbose (boolean) – print information for debugging
Returns:

  • D (tuple (B, Bpol)) – of ndarrays B (pol, proj, cols, cols) Bpol (pol, proj)
  • file.npy (file) – saves basis to file name linbasex_basis_{cols}_{legendre_orders}_{proj_angles}_{radial_step}_{clip}.npy

abel.linbasex.cache_cleanup()[source]

Utility function.

Frees the memory caches created by get_bs_cached(). This is usually pointless, but might be required after working with very large images, if more RAM is needed for further tasks.

Parameters:None
Returns:
Return type:None

abel.rbasex module

abel.rbasex.rbasex_transform(IM, origin='center', rmax='MIN', order=2, odd=False, weights=None, direction='inverse', reg=None, out='same', basis_dir=None, verbose=False)[source]

This function takes the input image and outputs its forward or inverse Abel transform as an image and its radial distributions.

The origin, rmax, order, odd and weights parameters are passed to abel.tools.vmi.Distributions, so see its documentation for their detailed descriptions.

Parameters:
  • IM (m × n numpy array) – the image to be transformed

  • origin (tuple of int or str) – image origin, explicit in the (row, column) format, or as a location string (by default, the image center)

  • rmax (int or string) – largest radius to include in the transform (by default, the largest radius with at least one full quadrant of data)

  • order (int) – highest angular order present in the data, ≥ 0 (by default, 2). Working with very high orders (≳ 15) can result in excessive noise, especially at small radii and for narrow peaks.

  • odd (bool) – include odd angular orders (by default is False, but is enabled automatically if order is odd)

  • weights (m × n numpy array, optional) – weighting factors for each pixel. The array shape must match the image shape. Parts of the image can be excluded from analysis by assigning zero weights to their pixels. By default is None, which applies equal weight to all pixels.

  • direction (str: 'forward' or 'inverse') – type of Abel transform to be performed (by default, inverse)

  • reg (None or str or tuple (str, float), optional) – regularization to use for inverse Abel transform. None (default) means no regularization, a string selects a non-parameterized regularization method, and parameterized methods are selected by a tuple (method, strength). Available methods are:

    ('L2', strength):

    Tikhonov \(L_2\) regularization with strength as the square of the Tikhonov factor. This is the same as “Tikhonov regularization” used in BASEX, with almost identical effects on the radial distributions.

    ('diff', strength):

    Tikhonov regularization with the difference operator (approximation of the derivative) multiplied by the square root of strength as the Tikhonov matrix. This tends to produce less blurring, but more negative overshoots than 'L2'.

    ('SVD', strength):

    truncated SVD (singular value decomposition) with N = strength × rmax largest singular values removed for each angular order. This mimics the approach proposed (but in fact not used) in pBasex. Not recommended due to generally poor results.

    'pos':

    non-parameterized method, finds the best (in the least-squares sense) solution with non-negative \(\cos^n\theta \sin^m\theta\) terms (see cossin()). For order = 0, 1, and 2 (with odd = False) this is equivalent to \(I(r, \theta) \geqslant 0\); for higher orders this assumption is stronger than \(I \geqslant 0\) and corresponds to no interference between different multiphoton channels. Not implemented for odd orders > 1.

    Notice that this method is nonlinear, which also means that it is considerably slower than the linear methods and the transform operator cannot be cached.

    In all cases, strength = 0 provides no regularization. For the Tikhonov methods, strength ~ 100 is a reasonable value for megapixel images. For truncated SVD, strength must be < 1; strength ~ 0.1 is a reasonable value; strength ~ 0.5 can produce noticeable ringing artifacts. See the full description and examples there.

  • out (str or None) – shape of the output image:

    'same' (default):

    same shape and origin as the input

    'fold' (fastest):

    Q0 (upper right) quadrant (for odd=False) or right half (for odd=True) up to rmax, but limited to the largest input-image quadrant (or half)

    'unfold':

    like 'fold', but symmetrically “unfolded” to all 4 quadrants

    'full':

    all pixels with radii up to rmax

    'full-unique':

    the unique part of 'full': Q0 (upper right) quadrant for odd=False, right half for odd=True

    None:

    no image (recon will be None). Can be useful to avoid unnecessary calculations when only the transformed radial distributions (distr) are needed.

  • basis_dir (str, optional) – path to the directory for saving / loading the basis set (useful only for the inverse transform without regularization; time savings in other cases are small and might be negated by the disk-access overhead). If None (default), the basis set will not be loaded from or saved to disk.

  • verbose (bool) – print information about processing (for debugging), disabled by default

Returns:

  • recon (2D numpy array or None) – the transformed image. Is centered and might have different dimensions than the input image.
  • distr (Distributions.Results object) – the object from which various distributions for the transformed image can be retrieved, see abel.tools.vmi.Distributions.Results

abel.rbasex.get_bs_cached(Rmax, order=2, odd=False, direction='inverse', reg=None, valid=None, basis_dir=None, verbose=False)[source]

Internal function.

Gets the basis set (from cache or runs computations and caches them) and calculates the transform matrix. Loaded/calculated matrices are also cached in memory.

Parameters:
  • Rmax (int) – largest radius to be transformed
  • order (int) – highest angular order
  • odd (bool) – include odd angular orders
  • direction (str: 'forward' or 'inverse') – type of Abel transform to be performed
  • reg (None or str or tuple (str, float)) – regularization type and strength for inverse transform
  • valid (None or bool array) – flags to exclude invalid radii from transform
  • basis_dir (str, optional) – path to the directory for saving / loading the basis set. If None, the basis set will not be saved to disk.
  • verbose (bool) – print some debug information
Returns:

A – (Rmax + 1) × (Rmax + 1) matrices of the Abel transform (forward or inverse) for each angular order

Return type:

list of 2D numpy arrays

abel.rbasex.cache_cleanup(select='all')[source]

Utility function.

Frees the memory caches created by get_bs_cached(). This is usually pointless, but might be required after working with very large images, if more RAM is needed for further tasks.

Parameters:

select (str) – selects which caches to clean:

all (default)

everything, including basis;

forward

forward transform;

inverse

inverse transform.

Returns:

Return type:

None

abel.hansenlaw module

abel.hansenlaw.hansenlaw_transform(image, dr=1, direction='inverse', hold_order=0, **kwargs)[source]

Forward/Inverse Abel transformation using the algorithm from

E. W. Hansen, “Fast Hankel transform algorithm”, IEEE Trans. Acoust. Speech Signal Proc. 33, 666–671 (1985)

and

E. W. Hansen, P.-L. Law, “Recursive methods for computing the Abel transform and its inverse”, J. Opt. Soc. Am. A 2, 510–520 (1985).

This function performs the Hansen–Law transform on only one “right-side” image:

Abeltrans = abel.hansenlaw.hansenlaw_transform(image, direction='inverse')

Note

Image should be a right-side image, like this:

.         +--------      +--------+
.         |      *       | *      |
.         |   *          |    *   |  <---------- im
.         |  *           |     *  |
.         +--------      o--------+
.         |  *           |     *  |
.         |   *          |    *   |
.         |     *        | *      |
.         +--------      +--------+

In accordance with all PyAbel methods the image origin o is defined to be mid-pixel i.e. an odd number of columns, for the full image.

For the full image transform, use the abel.Transform.

Inverse Abel transform:

iAbel = abel.Transform(image, method='hansenlaw').transform

Forward Abel transform:

fAbel = abel.Transform(image, direction='forward', method='hansenlaw').transform
Parameters:
  • image (1D or 2D numpy array) – Right-side half-image (or quadrant). See figure below.
  • dr (float) – Sampling size, used for Jacobian scaling. Default: 1 (appliable for pixel images).
  • direction (string ‘forward’ or ‘inverse’) – forward or inverse Abel transform. Default: ‘inverse’.
  • hold_order (int 0 or 1) – The order of the hold approximation, used to evaluate the state equation integral. 0 assumes a constant intensity across a pixel (between grid points) for the driving function (the image gradient for the inverse transform, or the original image, for the forward transform). 1 assumes a linear intensity variation between grid points, which may yield a more accurate transform for some functions (see PR 211). Default: 0.
Returns:

aim – forward/inverse Abel transform half-image

Return type:

1D or 2D numpy array

abel.dasch module

abel.dasch.two_point_transform(IM, basis_dir='.', dr=1, direction='inverse', verbose=False)[source]

The two-point deconvolution method.

C. J. Dasch, “One-dimensional tomography: a comparison of Abel, onion-peeling, and filtered backprojection methods”, Appl. Opt. 31, 1146–1152 (1992).

Parameters:
  • IM (1D or 2D numpy array) – right-side half-image (or quadrant)
  • basis_dir (str) – path to the directory for saving / loading the “two-point” deconvolution operator array. Here, called basis_dir for consistency with the other true basis methods. If None, the operator array will not be saved to disk.
  • dr (float) – sampling size (=1 for pixel images), used for Jacobian scaling. The resulting inverse transform is simply scaled by 1/dr.
  • direction (str) – only the direction=”inverse” transform is currently implemented
  • verbose (bool) – trace printing
Returns:

inv_IM – the “two-point” inverse Abel transformed half-image

Return type:

1D or 2D numpy array

abel.dasch.three_point_transform(IM, basis_dir='.', dr=1, direction='inverse', verbose=False)[source]

The three-point deconvolution method.

C. J. Dasch, “One-dimensional tomography: a comparison of Abel, onion-peeling, and filtered backprojection methods”, Appl. Opt. 31, 1146–1152 (1992).

Parameters:
  • IM (1D or 2D numpy array) – right-side half-image (or quadrant)
  • basis_dir (str) – path to the directory for saving / loading the “three-point” deconvolution operator array. Here, called basis_dir for consistency with the other true basis methods. If None, the operator array will not be saved to disk.
  • dr (float) – sampling size (=1 for pixel images), used for Jacobian scaling. The resulting inverse transform is simply scaled by 1/dr.
  • direction (str) – only the direction=”inverse” transform is currently implemented
  • verbose (bool) – trace printing
Returns:

inv_IM – the “three-point” inverse Abel transformed half-image

Return type:

1D or 2D numpy array

abel.dasch.onion_peeling_transform(IM, basis_dir='.', dr=1, direction='inverse', verbose=False)[source]

The onion-peeling deconvolution method.

C. J. Dasch, “One-dimensional tomography: a comparison of Abel, onion-peeling, and filtered backprojection methods”, Appl. Opt. 31, 1146–1152 (1992).

Parameters:
  • IM (1D or 2D numpy array) – right-side half-image (or quadrant)
  • basis_dir (str) – path to the directory for saving / loading the “onion-peeling” deconvolution operator array. Here, called basis_dir for consistency with the other true basis methods. If None, the operator array will not be saved to disk.
  • dr (float) – sampling size (=1 for pixel images), used for Jacobian scaling. The resulting inverse transform is simply scaled by 1/dr.
  • direction (str) – only the direction=”inverse” transform is currently implemented
  • verbose (bool) – trace printing
Returns:

inv_IM – the “onion-peeling” inverse Abel transformed half-image

Return type:

1D or 2D numpy array

abel.dasch.dasch_transform(IM, D)[source]

Inverse Abel transform using the given deconvolution D-operator array.

Parameters:
  • IM (2D numpy array) – image data
  • D (2D numpy array) – deconvolution operator array, of shape (cols, cols)
Returns:

inv_IM – inverse Abel transform according to deconvolution operator D

Return type:

2D numpy array

abel.dasch.get_bs_cached(method, cols, basis_dir='.', verbose=False)[source]

Load Dasch method deconvolution operator array from cache, or disk. Generate and store if not available.

Checks whether method deconvolution array has been previously calculated, or whether the file {method}_basis_{cols}.npy is present in basis_dir.

Either, assign, read, or generate the deconvolution array (saving it to file).

Parameters:
  • method (str) – Abel transform method onion_peeling, three_point, or two_point
  • cols (int) – width of image
  • basis_dir (str or None) – path to the directory for saving or loading the deconvolution array. For None do not save the deconvolution operator array
  • verbose (boolean) – print information (mainly for debugging purposes)
Returns:

  • D (numpy 2D array of shape (cols, cols)) – deconvolution operator array for the associated method
  • file.npy (file) – saves D, the deconvolution array to file name: {method}_basis_{cols}.npy

abel.dasch.cache_cleanup()[source]

Utility function.

Frees the memory caches created by get_bs_cached(). This is usually pointless, but might be required after working with very large images, if more RAM is needed for further tasks.

Parameters:None
Returns:
Return type:None

abel.onion_bordas module

abel.onion_bordas.onion_bordas_transform(IM, dr=1, direction='inverse', shift_grid=True, **kwargs)[source]

Onion peeling (or back projection) inverse Abel transform.

This algorithm was adapted by Dan Hickstein from the original Matlab implementation, created by Chris Rallis and Eric Wells of Augustana University, and described in

C. E. Rallis, T. G. Burwitz, P. R. Andrews, M. Zohrabi, R. Averin, S. De, B. Bergues, B. Jochim, A. V. Voznyuk, N. Gregerson, B. Gaire, I. Znakovskaya, J. McKenna, K. D. Carnes, M. F. Kling, I. Ben-Itzhak, E. Wells, “Incorporating real time velocity map image reconstruction into closed-loop coherent control”, Rev. Sci. Instrum. 85, 113105 (2014).

The algorithm actually originates from

C. Bordas, F. Paulig, “Photoelectron imaging spectrometry: Principle and inversion method”, Rev. Sci. Instrum. 67, 2257–2268 (1996).

This function operates on the “right side” of an image. i.e. it works on just half of a cylindrically symmetric image. Unlike the other transforms, the image origin should be at the left edge, not mid-pixel. This corresponds to an even-width full image.

However, shift_grid=True (default) provides the typical behavior, where the image origin corresponds to the pixel center in the 0th column.

To perform a onion-peeling transorm on a whole image, use

abel.Transform(image, method='onion_bordas').transform
Parameters:
  • IM (1D or 2D numpy array) – right-side half-image (or quadrant)
  • dr (float) – sampling size (=1 for pixel images), used for Jacobian scaling. The resulting inverse transform is simply scaled by 1/dr.
  • direction (str) – only the inverse transform is currently implemented.
  • shift_grid (bool) – place the image origin on the grid (left edge) by shifting the image 1/2 pixel to the left.
Returns:

AIM – the inverse Abel transformed half-image

Return type:

1D or 2D numpy array

abel.direct module

abel.direct.direct_transform(fr, dr=None, r=None, direction='inverse', derivative=<function gradient>, int_func=<function trapz>, correction=True, backend='C', **kwargs)[source]

This algorithm performs a direct computation of the Abel transform integrals. When correction=False, the pixel at the lower bound of the integral (where y=r) is skipped, which causes a systematic error in the Abel transform. However, if correction=True is used, then an analytical transform transform is applied to this pixel, which makes the approximation that the function is linear across this pixel. With correction=True, the Direct method produces reasonable results.

The Direct method is implemented in both Python and a compiled C version using Cython, which is much faster. The implementation can be selected using the backend argument. If the C-backend is not available, you must re-install PyAbel with Numpy, Cython, and a C-compiler already installed.

By default, integration at all other pixels is performed using the Trapezoidal rule.

Parameters:
  • fr (1d or 2d numpy array) – input array to which direct/inverse Abel transform will be applied. For a 2d array, the first dimension is assumed to be the z axis and the second the r axis.
  • dr (float) – spatial mesh resolution (optional, default to 1.0)
  • r (1D ndarray) – the spatial mesh (optional). Unusually, direct_transform should, in principle, be able to handle non-uniform data. However, this has not been regorously tested.
  • direction (string) – Determines if a forward or inverse Abel transform will be applied. can be ‘forward’ or ‘inverse’.
  • derivative (callable) – a function that can return the derivative of the fr array with respect to r. (only used in the inverse Abel transform).
  • int_func (function) – This function is used to complete the integration. It should resemble np.trapz, in that it must be callable using axis=, x=, and dx= keyword arguments.
  • correction (boolean) – If False the pixel where the weighting function has a singular value (where r==y) is simply skipped, causing a systematic under-estimation of the Abel transform. If True, integration near the singular value is performed analytically, by assuming that the data is linear across that pixel. The accuracy of this approximation will depend on how the data is sampled.
  • backend (string) – There are currently two implementations of the Direct transform, one in pure Python and one in Cython. The backend paremeter selects which method is used. The Cython code is converted to C and compiled, so this is faster. Can be ‘C’ or ‘python’ (case insensitive). ‘C’ is the default, but ‘python’ will be used if the C-library is not available.
Returns:

out – with either the direct or the inverse abel transform.

Return type:

1d or 2d numpy array of the same shape as fr

abel.direct.is_uniform_sampling(r)[source]

Returns True if the array is uniformly spaced to within 1e-13. Otherwise False.

Image processing tools

abel.tools.analytical module

class abel.tools.analytical.BaseAnalytical(n, r_max, symmetric=True, **args)[source]

Bases: object

class abel.tools.analytical.StepAnalytical(n, r_max, r1, r2, A0=1.0, ratio_valid_step=1.0, symmetric=True)[source]

Bases: abel.tools.analytical.BaseAnalytical

Define a symmetric step function and calculate its analytical Abel transform. See examples/example_basex_step.py.

Parameters:
  • n (int) – number of points along the r axis
  • r_max (float) – range of the r interval
  • symmetric (boolean) – if True, the r interval is [-r_max, r_max] (and n should be odd), otherwise the r interval is [0, r_max]
  • r1, r2 (float) – bounds of the step function if r > 0 (symmetric function is constructed for r < 0)
  • A0 (float) – height of the step
  • ratio_valid_step (float) – in the benchmark take only the central ratio*100% of the step (exclude possible artefacts on the edges)
abel_step_analytical(r, A0, r0, r1)[source]

Forward Abel transform of a step function located between r0 and r1, with a height A0

A0 +                  +-------------+
   |                  |             |
   |                  |             |
 0 | -----------------+             +-------------
   +------------------+-------------+------------>
   0                  r0            r1           r axis
Parameters:
  • r1 (1D array) – vecor of positions along the r axis. Must start with 0.
  • r0, r1 (float) – positions of the step along the r axis
  • A0 (float or 1D array) – height of the step. If 1D array, the height can be variable along the z axis
Returns:

Return type:

1D array, if A0 is a float, a 2D array otherwise

sym_abel_step_1d(r, A0, r0, r1)[source]

Produces a symmetrical analytical transform of a 1D step

class abel.tools.analytical.Polynomial(n, r_max, r_1, r_2, c, r_0=0.0, s=1.0, reduced=False, symmetric=True)[source]

Bases: abel.tools.analytical.BaseAnalytical

Define a polynomial function and calculate its analytical Abel transform.

(See Polynomials for details and examples.)

Parameters:
  • n (int) – number of points along the r axis
  • r_max (float) – range of the r interval
  • symmetric (boolean) – if True, the r interval is [−r_max, +r_max] (and n should be odd), otherwise the r interval is [0, r_max]
  • r_1, r_2 (float) – r bounds of the polynomial function if r > 0; outside [r_1, r_2] the function is set to zero (symmetric function is constructed for r < 0)
  • c (numpy array) – polynomial coefficients in order of increasing degree: [c₀, c₁, c₂] means c₀ + c₁ r + c₂ r²
  • r_0 (float, optional) – origin shift: the polynomial is defined as c₀ + c₁ (rr_0) + c₂ (rr_0)² + …
  • s (float, optional) – r stretching factor (around r_0): the polynomial is defined as c₀ + c₁ (r/s) + c₂ (r/s)² + …
  • reduced (boolean, optional) – internally rescale the r range to [0, 1]; useful to avoid floating-point overflows for high degrees at large r (and might improve numerical accuracy)
class abel.tools.analytical.PiecewisePolynomial(n, r_max, ranges, symmetric=True)[source]

Bases: abel.tools.analytical.BaseAnalytical

Define a piecewise polynomial function (sum of Polynomials) and calculate its analytical Abel transform.

Parameters:
  • n (int) – number of points along the r axis

  • r_max (float) – range of the r interval

  • symmetric (boolean) – if True, the r interval is [−r_max, +r_max] (and n should be odd), otherwise the r interval is [0, r_max]

  • ranges (iterable of unpackable) –

    (list of tuples of) polynomial parameters for each piece:

    [(r_1_1st, r_2_1st, c_1st),
     (r_1_2nd, r_2_2nd, c_2nd),
     ...
     (r_1_nth, r_2_nth, c_nth)]
    

    according to Polynomial conventions. All ranges are independent (may overlap and have gaps, may define polynomials of any degrees) and may include optional Polynomial parameters

class abel.tools.analytical.GaussianAnalytical(n, r_max, sigma=1.0, A0=1.0, ratio_valid_sigma=2.0, symmetric=True)[source]

Bases: abel.tools.analytical.BaseAnalytical

Define a gaussian function and calculate its analytical Abel transform. See examples/example_basex_gaussian.py.

Parameters:
  • n (int) – number of points along the r axis
  • r_max (float) – range of the r interval
  • symmetric (boolean) – if True, the r interval is [-r_max, r_max] (and n should be odd), otherwise, the r interval is [0, r_max]
  • sigma (floats) – sigma parameter for the gaussian
  • A0 (float) – amplitude of the gaussian
  • ratio_valid_sigma (float) – in the benchmark take only the range 0 < r < ration_valid_sigma * sigma (exclude possible artefacts on the axis and the possibly clipped tail)
class abel.tools.analytical.TransformPair(n, profile=5)[source]

Bases: abel.tools.analytical.BaseAnalytical

Abel-transform pair analytical functions.

profiles 1–7: Table 1 of G. C.-Y. Chan, Gary M. Hieftje, “Estimation of confidence intervals for radial emissivity and optimization of data treatment techniques in Abel inversion”, Spectrochimica Acta B 61, 31–41 (2006).

See abel.tools.transform_pairs.

Returns:
  • r (numpy array) – vector of positions along the r axis: linspace(0, 1, n)
  • dr (float) – radial interval
  • func (numpy array) – values of the original function (same shape as r)
  • abel (numpy array) – values of the Abel transform (same shape as func)
  • label (str) – name of the curve
  • mask_valid (boolean array) – set all True. Used for unit tests
class abel.tools.analytical.SampleImage(n=361, name='dribinski', sigma=2, temperature=200)[source]

Bases: abel.tools.analytical.BaseAnalytical

Sample images, made up of Gaussian functions

Parameters:
  • n (integer) – image size n rows x n cols
  • name (str) – one of “dribinski” or “Ominus”
  • sigma (float) – Gaussian 1/e width (pixels)
  • temperature (float) – for ‘Ominus’ only anion levels have Boltzmann population weight (2J+1) exp(-177.1 h c 100/k/temperature)
image

image

Type:2D numpy array
name

sample image name

Type:str

abel.tools.center module

abel.tools.center.find_origin(IM, method='image_center', axes=(0, 1), verbose=False, **kwargs)[source]

Find the coordinates of image origin, using the specified method.

Parameters:
  • IM (2D np.array) – image data

  • method (str) – determines how the origin should be found. The options are:

    image_center

    the center of the image is used as the origin. The trivial result.

    com

    the origin is found as the center of mass.

    convolution

    the origin is found as the maximum of autoconvolution of the image projections along each axis.

    gaussian

    the origin is extracted by a fit to a Gaussian function. This is probably only appropriate if the data resembles a gaussian.

    slice

    the image is broken into slices, and these slices compared for symmetry.

  • axes (int or tuple of int) – find origin coordinates: 0 (vertical), or 1 (horizontal), or (0, 1) (both vertical and horizontal).

Returns:

out – coordinates of the origin of the image in the (row, column) format. For coordinates not in axes, the center of the image is returned.

Return type:

(float, float)

abel.tools.center.center_image(IM, method='com', odd_size=True, square=False, axes=(0, 1), crop='maintain_size', order=3, verbose=False, center=<deprecated>, **kwargs)[source]

Center image with the custom value or by several methods provided in find_origin() function.

Parameters:
  • IM (2D np.array) – The image data.

  • method (str or tuple of float) – either a tuple (float, float), the coordinate of the origin of the image in the (row, column) format, or a string to specify an automatic centering method:

    image_center

    the center of the image is used as the origin. The trivial result.

    com

    the origin is found as the center of mass.

    convolution

    the origin is found as the maximum of autoconvolution of the image projections along each axis.

    gaussian

    the origin is extracted from a fit to a Gaussian function. This is probably only appropriate if the data resembles a gaussian.

    slice

    the image is broken into slices, and these slices compared for symmetry.

  • odd_size (boolean) – if True, the returned image will contain an odd number of columns. Most of the transform methods require this, so it’s best to set this to True if the image will subsequently be Abel-transformed.

  • square (bool) – if True, the returned image will have a square shape.

  • crop (str) – determines how the image should be cropped. The options are:

    maintain_size

    return image of the same size. Some regions of the original image may be lost, and some regions may be filled with zeros.

    valid_region

    return the largest image that can be created without padding. All of the returned image will correspond to the original image. However, portions of the original image will be lost. If you can tolerate clipping the edges of the image, this is probably the method to choose.

    maintain_data

    the image will be padded with zeros such that none of the original image will be cropped.

    See set_center() for examples.

  • axes (int or tuple of int) – center image with respect to axis 0 (vertical), 1 (horizontal), or both axes (0, 1) (default). When specifying an explicit origin in method, unused coordinates can also be passed as None, for example, method=(row, None) or method=(None, col).

  • order (int) – interpolation order, see set_center() for details.

Returns:

out – centered image

Return type:

2D np.array

abel.tools.center.set_center(data, origin, crop='maintain_size', axes=(0, 1), order=3, verbose=False, center=<deprecated>)[source]

Move image origin to mid-point of image (rows // 2, cols // 2).

Parameters:
  • data (2D np.array) – the image data

  • origin (tuple of float) – (row, column) coordinates of the image origin. Coordinates set to None are ignored.

  • crop (str) – determines how the image should be cropped. The options are:

    maintain_size (default)

    return image of the same size. Some regions of the original image may be lost and some regions may be filled with zeros.

    valid_region

    return the largest image that can be created without padding. All of the returned image will correspond to the original image. However, portions of the original image will be lost. If you can tolerate clipping the edges of the image, this is probably the method to choose.

    maintain_data

    the image will be padded with zeros such that none of the original image will be cropped.

    Examples:

    _images/crop_options.svg
  • axes (int or tuple of int) – center image with respect to axis 0 (vertical), 1 (horizontal), or both axes (0, 1) (default).

  • order (int) – interpolation order (0–5, default is 3) for centering with fractional origin. Lower orders work faster; order = 0 (also implied for integer origin) means a whole-pixel shift without interpolation and is much faster.

  • verbose (bool) – print some information for debugging

Returns:

out – centered image

Return type:

2D np.array

abel.tools.center.find_origin_by_center_of_mass(IM, axes=(0, 1), verbose=False, round_output=False, **kwargs)[source]

Find image origin by calculating its center of mass.

Parameters:
  • IM (numpy 2D array) – image data
  • round_output (bool) – if True, the coordinates are rounded to integers; otherwise they are floats.
  • axes (int or tuple) – find origin coordinates: 0 (vertical), or 1 (horizontal), or (0, 1) (both vertical and horizontal).
Returns:

origin – (row, column)

Return type:

(float, float)

abel.tools.center.find_origin_by_convolution(IM, axes=(0, 1), projections=False, **kwargs)[source]

Find the image origin as the maximum of autoconvolution of its projections along each axis.

Code from the linbasex juptyer notebook.

Parameters:
  • IM (numpy 2D array) – image data
  • projections (bool) – if this parameter is True, the autoconvoluted projections along both axes will be returned after the origin.
  • axes (int or tuple) – find origin coordinates: 0 (vertical), or 1 (horizontal), or (0, 1) (both vertical and horizontal).
Returns:

origin – (row, column)

or (row, column), conv_0, conv_1

Return type:

(float, float)

abel.tools.center.find_origin_by_center_of_image(data, axes=(0, 1), verbose=False, **kwargs)[source]

Find image origin simply as its center, from its dimensions.

Parameters:
  • IM (numpy 2D array) – image data
  • axes (int or tuple) – has no effect
Returns:

origin – (row, column)

Return type:

(int, int)

abel.tools.center.find_origin_by_gaussian_fit(IM, axes=(0, 1), verbose=False, round_output=False, **kwargs)[source]

Find image origin by fitting the summation along rows and columns of the data to two 1D Gaussian functions.

Parameters:
  • IM (numpy 2D array) – image data
  • axes (int or tuple) – find origin coordinates: 0 (vertical), or 1 (horizontal), or (0, 1) (both vertical and horizontal).
  • round_output (bool) – if True, the coordinates are rounded to integers; otherwise they are floats.
Returns:

origin – (row, column)

Return type:

(float, float)

abel.tools.center.axis_slices(IM, radial_range=(0, -1), slice_width=10)[source]

Returns vertical and horizontal slice profiles, summed across slice_width.

Parameters:
  • IM (2D np.array) – image data
  • radial_range (tuple of float) – (rmin, rmax) range to limit data
  • slice_width (int) – width of the image slice, default 10 pixels
Returns:

top, bottom, left, right – image slices oriented in the same direction

Return type:

1D np.arrays shape (rmin:rmax, 1)

abel.tools.center.find_origin_by_slice(IM, axes=(0, 1), slice_width=10, radial_range=(0, -1), axis=<deprecated>, **kwargs)[source]

Find the image origin by comparing opposite sides.

Parameters:
  • IM (2D np.array) – the image data
  • slice_width (integer) – Sum together this number of rows (cols) to improve signal, default 10.
  • radial_range (tuple) – (rmin, rmax): radial range [rmin:rmax] for slice profile comparison.
  • axes (int or tuple) – find origin coordinates: 0 (vertical), or 1 (horizontal), or (0, 1) (both vertical and horizontal).
Returns:

origin – (row, column)

Return type:

(float, float)

abel.tools.center.find_center(IM, center='image_center', square=False, verbose=False, **kwargs)[source]

Deprecated function. Use find_origin() instead.

abel.tools.center.find_center_by_center_of_mass(IM, verbose=False, round_output=False, **kwargs)[source]

Deprecated function. Use find_origin_by_center_of_mass() instead.

abel.tools.center.find_center_by_convolution(IM, **kwargs)[source]

Deprecated function. Use find_origin_by_convolution() instead.

abel.tools.center.find_center_by_center_of_image(data, verbose=False, **kwargs)[source]

Deprecated function. Use find_origin_by_center_of_image() instead.

abel.tools.center.find_center_by_gaussian_fit(IM, verbose=False, round_output=False, **kwargs)[source]

Deprecated function. Use find_origin_by_gaussian_fit() instead.

abel.tools.center.find_image_center_by_slice(IM, slice_width=10, radial_range=(0, -1), axis=(0, 1), **kwargs)[source]

Deprecated function. Use find_origin_by_slice() instead.

abel.tools.circularize module

abel.tools.circularize.circularize_image(IM, method='lsq', origin=None, radial_range=None, dr=0.5, dt=0.5, smooth=<deprecated>, ref_angle=None, inverse=False, return_correction=False, tol=0, center=<deprecated>)[source]

Corrects image distortion on the basis that the structure should be circular.

This is a simplified radial scaling version of the algorithm described in J. R. Gascooke, S. T. Gibson, W. D. Lawrance, “A ‘circularisation’ method to repair deformations and determine the centre of velocity map images”, J. Chem. Phys. 147, 013924 (2017).

This function is especially useful for correcting the image obtained with a velocity-map-imaging spectrometer, in the case where there is distortion of the Newton sphere (ring) structure due to an imperfect electrostatic lens or stray electromagnetic fields. The correction allows the highest-resolution 1D photoelectron distribution to be extracted.

The algorithm splits the image into “slices” at many different angles (set by dt) and compares the radial intensity profile of adjacent slices. A scaling factor is found which aligns each slice profile with the previous slice. The image is then corrected using a spline function that smoothly connects the discrete scaling factors as a continuous function of angle.

This circularization algorithm should only be applied to a well-centered image, otherwise use the origin keyword (described below) to center it.

Parameters:
  • IM (numpy 2D array) – Image to be circularized.

  • method (str) – Method used to determine the radial correction factor to align slice profiles:

    argmax

    compare intensity-profile.argmax() of each radial slice. This method is quick and reliable, but it assumes that the radial intensity profile has an obvious maximum. The positioning is limited to the nearest pixel.

    lsq

    minimize the difference between a slice intensity-profile with its adjacent slice. This method is slower and may fail to converge, but it may be applied to images with any (circular) structure. It aligns the slices with sub-pixel precision.

  • origin (float tuple, str or None) – Pre-center image using abel.tools.center.center_image(). May be an explicit (row, column) tuple or a method name: 'com', 'convolution', 'gaussian;, 'image_center', 'slice'. None (default) assumes that the image is already centered.

  • radial_range (tuple or None) – Limit slice comparison to the radial range tuple (rmin, rmax), in pixels, from the image origin. Use to determine the distortion correction associated with particular peaks. It is recommended to select a region of your image where the signal-to-noise ratio is highest, with sharp persistent (in angle) features.

  • dr (float) – Radial grid size for the polar coordinate image, default = 0.5 pixel. This is passed to abel.tools.polar.reproject_image_into_polar().

    Small values may improve the distortion correction, which is often of sub-pixel dimensions, at the cost of reduced signal to noise for the slice intensity profile. As a general rule, dr should be significantly smaller than the radial “feature size” in the image.

  • dt (float) – Angular grid size. This sets the number of radial slices, given by \(2\pi/dt\). Default = 0.1, ~ 63 slices. More slices, using smaller dt, may provide a more detailed angular variation of the correction, at the cost of greater signal to noise in the correction function.

    Also passed to abel.tools.polar.reproject_image_into_polar().

  • smooth (float) – Deprecated, use tol instead. The relationship is smooth = Nangles × tol2, where Nangles is the number of slices (see dt).

  • ref_angle (None or float) – Reference angle for which radial coordinate is unchanged. Angle varies between \(-\pi\) and \(\pi\), with zero angle vertical.

    None uses numpy.mean() of the radial correction function, which attempts to maintain the same average radial scaling. This approximation is likely valid, unless you know for certain that a specific angle of your image corresponds to an undistorted image.

  • inverse (bool) – Apply an inverse Abel transform to the polar-coordinate image, to remove the background intensity. This may improve the signal-to-noise ratio, allowing the weaker intensity featured to be followed in angle.

    Note that this step is only for the purposes of allowing the algorithm to better follow peaks in the image. It does not affect the final image that is returned, except for (hopefully) slightly improving the precision of the distortion correction.

  • return_correction (bool) – Additional outputs, as describe below.

  • tol (float) – Root-mean-square (RMS) fitting tolerance for the spline function. At the default zero value, the spline interpolates between the discrete scaling factors. At larger values, a smoother spline is found such that its RMS deviation from the discrete scaling factors does not exceed this number. For example, tol=0.01 means 1% RMS tolerance for the radial scaling correction. At very large tolerances, the spline degenerates to a constant, the average of the discrete scaling factors.

    Typically, tol may remain zero (use interpolation), but noisy data may require some smoothing, since the found discrete scaling factors can have noticeable errors. To examine the relative scaling factors and how well they are represented by the spline function, use the option return_correction=True.

Returns:

  • IMcirc (numpy 2D array) – Circularized version of the input image, same size as input.
  • The following values are returned if return_correction=True
  • angles (numpy 1D array) – Mid-point angle (radians) of each image slice.
  • radial_correction (numpy 1D array) – Radial correction scale factor at each angular slice.
  • radial_correction_function (function(numpy 1D array)) – Function that may be used to evaluate the radial correction at any angle.

abel.tools.circularize.circularize(IM, radial_correction_function, ref_angle=None)[source]

Remap image from its distorted grid to the true cartesian grid.

Parameters:
  • IM (numpy 2D array) – Original image
  • radial_correction_function (function(numpy 1D array)) – A function returning the radial correction for a given angle. It should accept a numpy 1D array of angles.
abel.tools.circularize.correction(polarIMTrans, angles, radial, method)[source]

Determines a radial correction factors that align an angular slice radial intensity profile with its adjacent (previous) slice profile.

Parameters:
  • polarIMTrans (numpy 2D array) – Polar coordinate image, transposed \((\theta, r)\) so that each row is a single angle.

  • angles (numpy 1D array) – Angle coordinates for one row of polarIMTrans.

  • radial (numpy 1D array) – Radial coordinates for one column of polarIMTrans.

  • method (str) –

    argmax

    radial correction factor from position of maximum intensity.

    lsq

    least-squares determine a radial correction factor that will align a radial intensity profile with the previous, adjacent slice.

abel.tools.math module

abel.tools.math.gradient(f, x=None, dx=1, axis=-1)[source]

Return the gradient of 1 or 2-dimensional array. The gradient is computed using central differences in the interior and first differences at the boundaries.

Irregular sampling is supported (it isn’t supported by np.gradient)

Parameters:
  • f (1d or 2d numpy array) – Input array.
  • x (array_like, optional) – Points where the function f is evaluated. It must be of the same length as f.shape[axis]. If None, regular sampling is assumed (see dx)
  • dx (float, optional) – If x is None, spacing given by dx is assumed. Default is 1.
  • axis (int, optional) – The axis along which the difference is taken.
Returns:

out – Returns the gradient along the given axis.

Return type:

array_like

Notes

To-Do: implement smooth noise-robust differentiators for use on experimental data. http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/

abel.tools.math.gaussian(x, a, mu, sigma, c)[source]

Gaussian function

\(f(x)=a e^{-(x - \mu)^2 / (2 \sigma^2)} + c\)

Parameters:
  • x (1D np.array) – coordinate
  • a (float) – the height of the curve’s peak
  • mu (float) – the position of the center of the peak
  • sigma (float) – the standard deviation, sometimes called the Gaussian RMS width
  • c (float) – non-zero background
Returns:

out – the Gaussian profile

Return type:

1D np.array

abel.tools.math.guss_gaussian(x)[source]

Find a set of better starting parameters for Gaussian function fitting

Parameters:x (1D np.array) – 1D profile of your data
Returns:out – estimated value of (a, mu, sigma, c)
Return type:tuple of float
abel.tools.math.fit_gaussian(x)[source]

Fit a Gaussian function to x and return its parameters

Parameters:x (1D np.array) – 1D profile of your data
Returns:out – (a, mu, sigma, c)
Return type:tuple of float

abel.tools.polar module

abel.tools.polar.reproject_image_into_polar(data, origin=None, Jacobian=False, dr=1, dt=None)[source]

Reprojects a 2D numpy array (data) into a polar coordinate system, with the pole placed at origin and the angle measured clockwise from the upward direction. The resulting array has rows corresponding to the radial grid, and columns corresponding to the angular grid.

Parameters:
  • data (2D np.array) – the image array
  • origin (tuple or None) – (row, column) coordinates of the image origin. If None, the geometric center of the image is used.
  • Jacobian (bool) – Include r intensity scaling in the coordinate transform. This should be included to account for the changing pixel size that occurs during the transform.
  • dr (float) – radial coordinate spacing for the grid interpolation. Tests show that there is not much point in going below 0.5.
  • dt (float or None) – angular coordinate spacing (in radians). If None, the number of angular grid points will be set to the largest dimension (the height or the width) of the image.
Returns:

  • output (2D np.array) – the polar image (r, theta)
  • r_grid (2D np.array) – meshgrid of radial coordinates
  • theta_grid (2D np.array) – meshgrid of angular coordinates

Notes

Adapted from: https://stackoverflow.com/questions/3798333/image-information-along-a-polar-coordinate-system

abel.tools.polar.index_coords(data, origin=None)[source]

Creates x and y coordinates for the indices in a numpy array, relative to the origin, with the x axis going to the right, and the y axis going up.

Parameters:
  • data (numpy array) – 2D data. Only the array shape is used.
  • origin (tuple or None) – (row, column). Defaults to the geometric center of the image.
Returns:

x, y

Return type:

2D numpy arrays

abel.tools.polar.cart2polar(x, y)[source]

Transform Cartesian coordinates to polar.

Parameters:x, y (floats or arrays) – Cartesian coordinates
Returns:r, theta – Polar coordinates
Return type:floats or arrays
abel.tools.polar.polar2cart(r, theta)[source]

Transform polar coordinates to Cartesian.

Parameters:r, theta (floats or arrays) – Polar coordinates
Returns:x, y – Cartesian coordinates
Return type:floats or arrays

abel.tools.polynomial module

See Polynomials for details and examples.

class abel.tools.polynomial.Polynomial(r, r_min, r_max, c, r_0=0.0, s=1.0, reduced=False)[source]

Bases: object

Polynomial function and its Abel transform.

Supports multiplication and division by numbers.

Parameters:
  • r (numpy array) – r values at which the function is generated (and x values for its Abel transform); must be non-negative and in ascending order
  • r_min, r_max (float) – r domain: the function is defined as the polynomial on [r_min, r_max] and zero outside it; 0 ≤ r_min < r_maxmax r (r_max might exceed maximal r, but usually by < 1 pixel)
  • c (numpy array) – polynomial coefficients in order of increasing degree: [c₀, c₁, c₂] means c₀ + c₁ r + c₂ r²
  • r_0 (float, optional) – origin shift: the polynomial is defined as c₀ + c₁ (rr_0) + c₂ (rr_0)² + …
  • s (float, optional) – r stretching factor (around r_0): the polynomial is defined as c₀ + c₁ (r/s) + c₂ (r/s)² + …
  • reduced (boolean, optional) – internally rescale the r range to [0, 1]; useful to avoid floating-point overflows for high degrees at large r (and might improve numeric accuracy)
class abel.tools.polynomial.PiecewisePolynomial(r, ranges)[source]

Bases: abel.tools.polynomial.Polynomial

Piecewise polynomial function (sum of Polynomials) and its Abel transform.

Supports multiplication and division by numbers.

Parameters:
  • r (numpy array) – r values at which the function is generated (and x values for its Abel transform)

  • ranges (iterable of unpackable) –

    (list of tuples of) polynomial parameters for each piece:

    [(r_min_1st, r_max_1st, c_1st),
     (r_min_2nd, r_max_2nd, c_2nd),
     ...
     (r_min_nth, r_max_nth, c_nth)]
    

    according to Polynomial conventions. All ranges are independent (may overlap and have gaps, may define polynomials of any degrees) and may include optional Polynomial parameters

abel.tools.transform_pairs module

Analytical function Abel-transform pairs

profiles 1–7:
G. C.-Y. Chan, Gary M. Hieftje, “Estimation of confidence intervals for radial emissivity and optimization of data treatment techniques in Abel inversion”, Spectrochimica Acta B 61, 31–41 (2006), Table 1.

Note

the transform pair functions are more conveniently accessed through abel.tools.analytical.TransformPair:

func = abel.tools.analytical.TransformPair(n, profile=nprofile)

which sets the radial range r and provides attributes .func (source), .abel (projection), .r (radial range), .dr (step), .label (the profile name)

Parameters:r (floats or numpy 1D array of floats) – value or grid to evaluate the function pair: 0 < r < 1
returns:source, projection – source function profile (inverse Abel transform of projection), projection functon profile (forward Abel transform of source)
rtype:tuple of 1D numpy arrays of shape r
abel.tools.transform_pairs.a(n, r)[source]

coefficient

\[a_n = \sqrt{n^2 - r^2}\]
abel.tools.transform_pairs.profile1(r)[source]

profile1: C. J. Cremers, R. C. Birkebak, “Application of the Abel Integral Equation to Spectrographic Data”, Appl. Opt. 5, 1057–1064 (1966), Eq. (13).

\[ \begin{align}\begin{aligned}\epsilon(r) &= 0.75 + 12r^2 - 32r^3 & 0 \le r \le 0.25\\\epsilon(r) &= \frac{16}{27}(1 + 6r - 15r^2 + 8r^3) & 0.25 \lt r \le 1\\I(r) &= \frac{1}{108}(128a_1 +a_{0.25}) + \frac{2}{27}r^2 (283a_{0.25} - 112a_1) +\\& \,\,\,\, \frac{8}{9}r^2\left[4(1+r^2)\ln\frac{1+a_1}{r} - (4+31r^2)\ln\frac{0.25+a_{0.25}}{r}\right] & 0 \le r \le 0.25\\I(r) &= \frac{32}{27}\left[a_1 - 7a_1 r + 3r^2(1+r^2) \ln\frac{1+a_1}{r}\right] & 0.25 \lt r \le 1\end{aligned}\end{align} \]
_images/abel-1.svg
abel.tools.transform_pairs.profile2(r)[source]

profile2: C. J. Cremers, R. C. Birkebak, “Application of the Abel Integral Equation to Spectrographic Data”, Appl. Opt. 5, 1057–1064 (1966), Eq. (11).

\[ \begin{align}\begin{aligned}\epsilon(r) &= 1 - 3r^2 + 2r^3 & 0 \le r \le 1\\I(r) &= a_1\left(1-\frac{5}{2}r^2\right) + \frac{3}{2}r^4\ln\frac{1+a_1}{r} & 0 \le r \le 1\end{aligned}\end{align} \]
_images/abel-2.svg
abel.tools.transform_pairs.profile3(r)[source]

profile3: C. J. Cremers, R. C. Birkebak, “Application of the Abel Integral Equation to Spectrographic Data”, Appl. Opt. 5, 1057–1064 (1966), Eq. (12).

\[ \begin{align}\begin{aligned}\epsilon(r) &= 1-2r^2 & 0 \le r \le 0.5\\\epsilon(r) &= 2(1-r^2)^2 & 0.5 \lt r \le 1\\I(r) &= \frac{4a_1}{3}(1+2r^2)-\frac{2 a_{0.5}}{3}(1+8r^2) - 4r^2\ln\frac{1-a_1}{0.5+a_{0.5}} & 0 \le r \le 0.5\\I(r) &= \frac{4a_1}{3}(1+2r^2)-4r^2\ln\frac{1-a_1}{r} & 0.5 \lt r \le 1\end{aligned}\end{align} \]
_images/abel-3.svg
abel.tools.transform_pairs.profile4(r)[source]

profile4: R. Álvarez, A. Rodero, M. C. Quintero, “An Abel inversion method for radially resolved measurements in the axial injection torch”, Spectochim. Acta B 57, 1665–1680 (2002), Eq. (10).

Note

Published projection has misprints (“193.30083” instead of “196.30083” in both cases).

\[ \begin{align}\begin{aligned}\epsilon(r) &= 0.1 + 5.51r^2 - 5.25r^3 & 0 \le r \le 0.7\\\epsilon(r) &= -40.74 + 155.56r - 188.89r^2 + 74.07r^3 & 0.7 \lt r \le1\\I(r) &= 22.68862a_{0.7} - 14.811667a_1 + (217.557a_{0.7} - 196.30083a_1)r^2 +\\ & \,\,\, 155.56r^2\ln\frac{1 + a_1}{0.7 + a_{0.7}} + r^4\left(55.5525\ln\frac{1 + a_1}{r} - 59.49\ln\frac{0.7 + a_{0.7}}{r}\right) & 0 \le r \le 0.7\\I(r) &= -14.811667a_1 - 196.30083a_1 r^2 + r^2(155.56 + 55.5525r^2) \ln\frac{1 + a_1}{r} & 0.7 \lt r \le 1\end{aligned}\end{align} \]
_images/abel-4.svg
abel.tools.transform_pairs.profile5(r)[source]

profile5: M. J. Buie, J. T. P. Pender, J. P. Holloway, T. Vincent, P. L. G. Ventzek, M. L. Brake, J. Quant. Spectrosc. Radiat. Transf. 55, 231–243 (1996), Table 1, № 1.

\[ \begin{align}\begin{aligned}\epsilon(r) &= 1 & 0 \le r \le 1\\I(r) &= 2a_1 & 0 \le r \le 1\end{aligned}\end{align} \]
_images/abel-5.svg
abel.tools.transform_pairs.profile6(r)[source]

profile6: M. J. Buie, J. T. P. Pender, J. P. Holloway, T. Vincent, P. L. G. Ventzek, M. L. Brake, J. Quant. Spectrosc. Radiat. Transf. 55, 231–243 (1996), Table 1, № 7.

\[ \begin{align}\begin{aligned}\epsilon(r) &= (1-r^2)^{-\frac{3}{2}} \exp\left[1.1^2\left( 1 - \frac{1}{1-r^2}\right)\right] & 0 \le r \le 1\\I(r) &= \frac{\sqrt{\pi}}{1.1a_1} \exp\left[1.1^2\left( 1 - \frac{1}{1-r^2}\right)\right] & 0 \le r \le 1\end{aligned}\end{align} \]
_images/abel-6.svg
abel.tools.transform_pairs.profile7(r)[source]

profile7: M. J. Buie, J. T. P. Pender, J. P. Holloway, T. Vincent, P. L. G. Ventzek, M. L. Brake, J. Quant. Spectrosc. Radiat. Transf. 55, 231–243 (1996), Table 1, № 9 (divided by 2).

\[ \begin{align}\begin{aligned}\epsilon(r) &= \frac{1}{2}(1+10r^2-23r^4+12r^6) & 0 \le r \le 1\\I(r) &= \frac{8}{105}a_1(19 + 34r^2 - 125r^4 + 72r^6) & 0 \le r \le 1\end{aligned}\end{align} \]
_images/abel-7.svg

abel.tools.symmetry module

abel.tools.symmetry.get_image_quadrants(IM, reorient=True, symmetry_axis=None, use_quadrants=(True, True, True, True), symmetrize_method='average')[source]

Given an image (m,n) return its 4 quadrants Q0, Q1, Q2, Q3 as defined below.

Parameters:
  • IM (2D np.array) – Image data shape (rows, cols)

  • reorient (boolean) – Reorient quadrants to match the orientation of Q0 (top-right)

  • symmetry_axis (int or tuple) – can have values of None, 0, 1, or (0,1) and specifies no symmetry, vertical symmetry axis, horizontal symmetry axis, and both vertical and horizontal symmetry axes. Quadrants are added. See Note.

  • use_quadrants (boolean tuple) – Include quadrant (Q0, Q1, Q2, Q3) in the symmetry combination(s) and final image

  • symmetrize_method (str) – Method used for symmetrizing the image.

    average

    Simply average the quadrants.

    fourier

    Axial symmetry implies that the Fourier components of the 2-D projection should be real. Removing the imaginary components in reciprocal space leaves a symmetric projection.

    K. R. Overstreet, P. Zabawa, J. Tallant, A. Schwettmann, J. P. Shaffer, “Multiple scattering and the density distribution of a Cs MOT”, Optics Express 13, 9672–9682 (2005).

Returns:

Q0, Q1, Q2, Q3 – shape: (rows//2+rows%2, cols//2+cols%2) all oriented in the same direction as Q0 if reorient=True

Return type:

tuple of 2D np.arrays

Notes

The symmetry_axis keyword averages quadrants like this:

+--------+--------+
| Q1   * | *   Q0 |
|   *    |    *   |
|  *     |     *  |               cQ1 | cQ0
+--------o--------+ --(output) -> ----o----
|  *     |     *  |               cQ2 | cQ3
|   *    |    *   |
| Q2  *  | *   Q3 |          cQi == combined quadrants
+--------+--------+

symmetry_axis = None - individual quadrants
symmetry_axis = 0 (vertical) - average Q0+Q1, and Q2+Q3
symmetry_axis = 1 (horizontal) - average Q1+Q2, and Q0+Q3
symmetry_axis = (0, 1) (both) - combine and average all 4 quadrants

The end results look like this:

(0) symmetry_axis = None

    returned image   Q1 | Q0
                    ----o----
                     Q2 | Q3

(1) symmetry_axis = 0

    Combine:  Q01 = Q0 + Q1, Q23 = Q2 + Q3
    returned image   Q01 | Q01
                    -----o-----
                     Q23 | Q23

(2) symmetry_axis = 1

    Combine: Q12 = Q1 + Q2, Q03 = Q0 + Q3
    returned image   Q12 | Q03
                    -----o-----
                     Q12 | Q03

(3) symmetry_axis = (0, 1)

    Combine all quadrants: Q = Q0 + Q1 + Q2 + Q3
    returned image   Q | Q
                    ---o---  all quadrants equivalent
                     Q | Q
abel.tools.symmetry.put_image_quadrants(Q, original_image_shape, symmetry_axis=None)[source]

Reassemble image from 4 quadrants Q = (Q0, Q1, Q2, Q3) The reverse process to get_image_quadrants(reorient=True)

Note: the quadrants should all be oriented as Q0, the upper right quadrant

Parameters:
  • Q (tuple of np.array (Q0, Q1, Q2, Q3)) – Image quadrants all oriented as Q0 shape (rows//2+rows%2, cols//2+cols%2)

    +--------+--------+
    | Q1   * | *   Q0 |
    |   *    |    *   |
    |  *     |     *  |
    +--------o--------+
    |  *     |     *  |
    |   *    |    *   |
    | Q2  *  | *   Q3 |
    +--------+--------+
    
  • original_image_shape (tuple) – (rows, cols)

    reverses the padding added by get_image_quadrants() for odd-axis sizes

    odd row trims 1 row from Q1, Q0

    odd column trims 1 column from Q1, Q2

  • symmetry_axis (int or tuple) –

    impose image symmetry

    symmetry_axis = 0 (vertical)   - Q0 == Q1 and Q3 == Q2 symmetry_axis = 1 (horizontal) - Q2 == Q1 and Q3 == Q0

Returns:

IM

Reassembled image of shape (rows, cols):

symmetry_axis =

None             0              1           (0,1)

 Q1 | Q0        Q1 | Q1        Q1 | Q0       Q1 | Q1
----o----  or  ----o----  or  ----o----  or ----o----
 Q2 | Q3        Q2 | Q2        Q1 | Q0       Q1 | Q1

Return type:

np.array

abel.tools.vmi module

abel.tools.vmi.angular_integration(IM, origin=None, Jacobian=True, dr=1, dt=None)[source]

Angular integration of the image.

Returns the one-dimensional intensity profile as a function of the radial coordinate.

Note: the use of Jacobian=True applies the correct Jacobian for the integration of a 3D object in spherical coordinates.

Parameters:
  • IM (2D numpy.array) – the image data
  • origin (tuple or None) – image origin in the (row, column) format. If None, the geometric center of the image (rows // 2, cols // 2) is used.
  • Jacobian (bool) – Include \(r\sin\theta\) in the angular sum (integration). Also, Jacobian=True is passed to abel.tools.polar.reproject_image_into_polar(), which includes another value of r, thus providing the appropriate total Jacobian of \(r^2\sin\theta\).
  • dr (float) – radial grid spacing in pixels (default 1). dr=0.5 may reduce pixel granularity of the speed profile.
  • dt (float or None) – angular grid spacing in radians. If None, the number of theta values will be set to largest dimension (the height or the width) of the image, which should typically ensure good sampling.
Returns:

  • r (1D numpy.array) – radial coordinates
  • speeds (1D numpy.array) – integrated intensity array (vs radius).

abel.tools.vmi.average_radial_intensity(IM, **kwargs)[source]

Calculate the average radial intensity of the image, averaged over all angles. This differs form abel.tools.vmi.angular_integration() only in that it returns the average intensity, and not the integrated intensity of a 3D image. It is equivalent to calling abel.tools.vmi.angular_integration() with Jacobian=True and then dividing the result by 2π.

Parameters:
Returns:

  • r (1D numpy.array) – radial coordinates
  • intensity (1D numpy.array) – intensity profile as a function of the radial coordinate

abel.tools.vmi.radial_integration(IM, origin=None, radial_ranges=None)[source]

Intensity variation in the angular coordinate.

This function is the \(\theta\)-coordinate complement to abel.tools.vmi.angular_integration().

Evaluates intensity vs angle for defined radial ranges. Determines the anisotropy parameter for each radial range.

See examples/example_O2_PES_PAD.py.

Parameters:
  • IM (2D numpy.array) – the image data

  • origin (tuple or None) – image origin in the (row, column) format. If None, the geometric center of the image (rows // 2, cols // 2) is used.

  • radial_ranges (list of tuple ranges or int step) –

    tuple

    integration ranges [(r0, r1), (r2, r3), ...] evaluates the intensity vs angle for the radial ranges r0_r1, r2_r3, etc.

    int

    the whole radial range (0, step), (step, 2*step), ..

Returns:

  • Beta (array of tuples) – (beta0, error_beta_fit0), (beta1, error_beta_fit1), … corresponding to the radial ranges
  • Amplitude (array of tuples) – (amp0, error_amp_fit0), (amp1, error_amp_fit1), … corresponding to the radial ranges
  • Rmidpt (numpy float 1D array) – radial mid-point of each radial range
  • Intensity_vs_theta (2D numpy.array) – intensity vs angle distribution for each selected radial range
  • theta (1D numpy.array) – angle coordinates, referenced to vertical direction

abel.tools.vmi.anisotropy_parameter(theta, intensity, theta_ranges=None)[source]

Evaluate anisotropy parameter \(\beta\), for \(I\) vs \(\theta\) data:

\[I = \frac{\sigma_\text{total}}{4\pi} [ 1 + \beta P_2(\cos\theta) ],\]

where \(P_2(x)=\frac{3x^2-1}{2}\) is a 2nd-order Legendre polynomial.

J. Cooper, R. N. Zare, “Angular Distribution of Photoelectrons”, J. Chem. Phys. 48, 942–943 (1968).

Parameters:
  • theta (1D numpy array) – angle coordinates, referenced to the vertical direction.
  • intensity (1D numpy array) – intensity variation with angle
  • theta_ranges (list of tuples) – angular ranges over which to fit [(theta1, theta2), (theta3, theta4)]. Allows data to be excluded from fit, default include all data.
Returns:

  • beta (tuple of floats) – (anisotropy parameter, fit error)
  • amplitude (tuple of floats) – (amplitude of signal, fit error)

abel.tools.vmi.toPES(radial, intensity, energy_cal_factor, per_energy_scaling=True, photon_energy=None, Vrep=None, zoom=1)[source]

Convert speed radial coordinate into electron kinetic or electron binding energy. Return the photoelectron spectrum (PES).

This calculation uses a single scaling factor energy_cal_factor to convert the radial pixel coordinate into electron kinetic energy.

Additional experimental parameters: photon_energy will give the energy scale as electron binding energy, in the same energy units, while Vrep, the VMI lens repeller voltage (volts), provides for a voltage-independent scaling factor. i.e. energy_cal_factor should remain approximately constant.

The energy_cal_factor is readily determined by comparing the generated energy scale with published spectra. e.g. for O photodetachment, the strongest fine-structure transition occurs at the electron affinity \(EA = 11\,784.676(7)\) cm\(^{-1}\). Values for the ANU experiment are given below, see also examples/example_hansenlaw.py.

Parameters:
  • radial (numpy 1D array) – radial coordinates.
  • intensity (numpy 1D array) – intensity values, at the radial array.
  • energy_cal_factor (float) – energy calibration factor that will convert radius squared into energy. The units affect the units of the output. e.g. inputs in eV/pixel2, will give output energy units in eV. A value of \(1.148427\times 10^{-5}\) cm\(^{-1}/\)pixel2 applies for “examples/data/O-ANU1024.txt” (with Vrep = -98 volts).
  • per_energy_scaling (bool) – sets the intensity Jacobian. If True, the returned intensities correspond to an “intensity per eV” or “intensity per cm-1 “. If False, the returned intensities correspond to an “intensity per pixel”.
  • photon_energy (None or float) – measurement photon energy. The output energy scale is then set to electron-binding-energy in units of energy_cal_factor. The conversion from wavelength (nm) to photon_energy (cm−1) is \(10^{7}/\lambda\) (nm) e.g. 1.0e7/812.51 for “examples/data/O-ANU1024.txt”.
  • Vrep (None or float) – repeller voltage. Convenience parameter to allow the energy_cal_factor to remain constant, for different VMI lens repeller voltages. Defaults to None, in which case no extra scaling is applied. e.g. -98 (volts), for “examples/data/O-ANU1024.txt”.
  • zoom (float) – additional scaling factor if the input experimental image has been zoomed. Default 1.
Returns:

  • eKBE (numpy 1D array of floats) – energy scale for the photoelectron spectrum in units of energy_cal_factor. Note that the data is no longer on a uniform grid.
  • PES (numpy 1D array of floats) – the photoelectron spectrum, scaled according to the per_energy_scaling input parameter.

class abel.tools.vmi.Distributions(origin='center', rmax='MIN', order=2, odd=False, use_sin=True, weights=None, method='linear')[source]

Bases: object

Class for calculating various radial distributions.

Objects of this class hold the analysis parameters and cache some intermediate computations that do not depend on the image data. Multiple images can be analyzed (using the same parameters) by feeding them to the object:

distr = Distributions(parameters)
results1 = distr(image1)
results2 = distr(image2)

If analyses with different parameters are required, multiple objects can be used. For example, to analyze 4 quadrants independently:

distr0 = Distributions('ll', ...)
distr1 = Distributions('lr', ...)
distr2 = Distributions('ur', ...)
distr3 = Distributions('ul', ...)

for image in images:
    Q0, Q1, Q2, Q3 = ...
    res0 = distr0(Q0)
    res1 = distr1(Q1)
    res2 = distr2(Q2)
    res3 = distr3(Q3)

However, if all the quadrants have the same dimensions, it is more memory-efficient to flip them all to the same orientation and use a single object:

distr = Distributions('ll', ...)

for image in images:
    Q0, Q1, Q2, Q3 = ...
    res0 = distr(Q0)
    res1 = distr(Q1[:, ::-1])  # or np.fliplr
    res2 = distr(Q2[::-1, ::-1])  # or np.flip(Q2, (0, 1))
    res3 = distr(Q3[::-1, :])  # or np.flipud

More concise function to calculate distributions for single images (without caching) are also available, see harmonics(), Ibeta() below.

Parameters:
  • origin (tuple of int or str) – origin of the radial distributions (the pole of polar coordinates) within the image.

    (int, int):

    explicit row and column indices

    str:

    location string specifying the vertical and horizontal positions (in this order!) using the words from the following diagram:

                  left            center             right
    
       top/upper  [0, 0]---------[0, n//2]--------[0, n-1]
                  |                                      |
                  |                                      |
          center  [m//2, 0]    [m//2, n//2]    [m//2, n-1]
                  |                                      |
                  |                                      |
    bottom/lower  [m-1, 0]------[m-1, n//2]-----[m-1, n-1]
    

    The words can be abbreviated to their first letter each (such as 'top left''tl', the space is then not required).

    'center center'/'cc' can also be shortened to 'center'/'c'.

    Examples:

    'center' or 'cc' (default) for the full centered image

    'center left'/'cl' for the right image half, vertically centered

    'bottom left'/'bl' or 'lower left'/'ll' for the upper-right image quadrant

  • rmax (int or str) – largest radius to include in the distributions

    int:

    explicit value

    'hor':

    fitting inside horizontally

    'ver':

    fitting inside vertically

    'HOR':

    touching horizontally

    'VER':

    touching vertically

    'min':

    minimum of 'hor' and 'ver', the largest area with 4 full quadrants

    'max':

    maximum of 'hor' and 'ver', the largest area with 2 full quadrants

    'MIN' (default):

    minimum of 'HOR' and 'VER', the largest area with 1 full quadrant (thus the largest with the full 90° angular range)

    'MAX':

    maximum of 'HOR' and 'VER'

    'all':

    covering all pixels (might have huge errors at large r, since the angular dependences must be inferred from very small available angular ranges)

  • order (int) – highest order in the angular distributions, ≥ 0 (by default, 2). Requesting very high orders (≳ 15) can result in excessive noise, especially at small radii and for narrow peaks.

  • odd (bool) – include odd angular orders. By default is False, but is enabled automatically if order is odd. Notice that although odd orders can be extracted from the upper or lower image part alone, analyzing the whole image is more reliable.

  • use_sin (bool) – use \(|\sin \theta|\) weighting (enabled by default). This is the weight implied in spherical integration (for the total intensity, for example) and with respect to which the Legendre polynomials are orthogonal, so using it in the fitting procedure gives the most reasonable results even if the data deviates form the assumed angular behavior. It also reduces contributions from the centerline noise.

  • weights (m × n numpy array, optional) – in addition to the optional \(|\sin \theta|\) weighting (see use_sin above), use given weights for each pixel. The array shape must match the image shape.

    Parts of the image can be excluded from the fitting by assigning zero weights to their pixels.

    (Note: if use_sin=False, a reference to this array is cached instead of its content, so if you modify the array between creating the object and using it, the results will be surprising. However, if needed, you can pass a copy as weights=weights.copy().)

  • method (str) – numerical integration method used in the fitting procedure

    'nearest':

    each pixel of the image is assigned to the nearest radial bin. The fastest, but noisier (especially for high orders).

    'linear' (default):

    each pixel of the image is linearly distributed over the two adjacent radial bins. About twice slower than 'nearest', but smoother.

    'remap':

    the image is resampled to a uniform polar grid, then polar pixels are summed over all angles for each radius. The smoothest, but significantly slower and might have problems with rmax > 'MIN' and discontinuous weights.

class Results(r, cn, order, odd, valid=None)[source]

Bases: object

Class for holding the results of image analysis.

Distributions.image() returns an object of this class, from which various distributions can be retrieved using the methods described below, for example:

distr = Distributions(...)
res = distr(IM)
harmonics = res.harmonics()

All distributions are returned as 2D arrays with the rows (1st index) corresponding to particular terms of the expansion and the columns (2nd index) corresponding to the radii. Odd angular terms are included only when they are used (odd = True or order is odd), otherwise there are only 1 + order/2 rows. The terms can be easily separated like I, beta2, beta4 = res.Ibeta(). Python 3 users can also collect all \(\beta\) parameters as I, *beta = res.Ibeta() for any order. Alternatively, transposing the results as Ibeta = res.Ibeta().T allows accessing all terms \(\big(I(r), \beta_2(r), \beta_4(r), \dots\big)\) at particular radius r as Ibeta[r].

r

radii from 0 to rmax

Type:numpy array
order

highest order in the angular distributions

Type:int
odd

whether odd angular orders are present

Type:bool
orders

orders for all angular terms:

[0, 2, …, order] for odd = False,

[0, 1, 2, …, order] for odd = True

Type:list of int
sinpowers

sine powers \(m\) in the \(\cos^n\theta \cdot \sin^m\theta\) terms from cossin(); cosine powers \(n\) are given by orders (see above)

Type:list of int
valid

flags for each radius indicating whether it has valid data (radii that have zero weights for all pixels will have no valid data)

Type:bool array
cos()[source]

Radial distributions of \(\cos^n \theta\) terms (0 ≤ norder).

(You probably do not need them.)

Returns:cosn – radial dependences of the \(\cos^n \theta\) terms, ordered from the lowest to the highest power
Return type:(# terms) × (rmax + 1) numpy array
rcos()[source]

Same as cos(), but prepended with the radii row.

cossin()[source]

Radial distributions of \(\cos^n \theta \cdot \sin^m \theta\) terms (n + m = order, and n + m = order − 1 for odd orders, with m always even).

For order = 0:

\(\cos^0 \theta\) is the total intensity.

For order = 1:

\(\cos^0 \theta\) is the total intensity,

\(\cos^1 \theta\) is the antisymmetric component.

For order = 2

\(\sin^2 \theta\) corresponds to “perpendicular” (⟂) transitions,

\(\cos^2 \theta\) corresponds to “parallel” (∥) transitions.

For order = 4

\(\sin^4 \theta\) corresponds to ⟂,⟂,

\(\cos^2 \theta \cdot \sin^2 \theta\) corresponds to ∥,⟂ and ⟂,∥,

\(\cos^4 \theta\) corresponds to ∥,∥.

And so on.

Notice that higher orders can represent lower orders as well:

\(\sin^2 \theta + \cos^2 \theta= \cos^0 \theta \quad\) (⟂ + ∥ = 1),

\(\sin^4 \theta + \cos^2 \theta \cdot \sin^2 \theta = \sin^2 \theta \quad\) (⟂,⟂ + ∥,⟂ = ⟂,⟂ + ⟂,∥ = ⟂),

\(\cos^2 \theta \cdot \sin^2 \theta + \cos^4 \theta = \cos^2 \theta \quad\) (∥,⟂ + ∥,∥ = ⟂,∥ + ∥,∥ = ∥),

and so forth.

Returns:cosnsinm – radial dependences of the \(\cos^n \theta \cdot \sin^m \theta\) terms, ordered from lower to higher \(\cos \theta\) powers
Return type:(# terms) × (rmax + 1) numpy array
rcossin()[source]

Same as cossin(), but prepended with the radii row.

harmonics()[source]

Radial distributions of spherical harmonics (Legendre polynomials \(P_n(\cos \theta)\)).

Spherical harmonics are orthogonal with respect to integration over the full sphere:

\[\iint P_n P_m \,d\Omega = \int_0^{2\pi} \int_0^\pi P_n(\cos \theta) P_m(\cos \theta) \,\sin\theta d\theta \,d\varphi = 0\]

for nm; and \(P_0(\cos \theta)\) is the spherically averaged intensity.

Returns:Pn – radial dependences of the \(P_n(\cos \theta)\) terms
Return type:(# terms) × (rmax + 1) numpy array
rharmonics()[source]

Same as harmonics(), but prepended with the radii row.

Ibeta(window=1)[source]

Radial intensity and anisotropy distributions.

A cylindrically symmetric 3D intensity distribution can be expanded over spherical harmonics (Legendre polynomials \(P_n(\cos \theta)\)) as (including even and odd terms)

\[I(r, \theta, \varphi) \, d\Omega = \frac{1}{4\pi} I(r) \big[1 + \beta_1(r) P_1(\cos \theta) + \beta_2(r) P_2(\cos \theta) + \dots\big],\]

or, for distributions with top–bottom symmetry (only even terms),

\[I(r, \theta, \varphi) \, d\Omega = \frac{1}{4\pi} I(r) \big[1 + \beta_2(r) P_2(\cos \theta) + \beta_4(r) P_4(\cos \theta) + \dots\big],\]

where \(I(r)\) is the “radial intensity distribution” integrated over the full sphere:

\[I(r) = \int_0^{2\pi} \int_0^\pi I(r, \theta, \varphi) \,r^2 \sin\theta d\theta \,d\varphi,\]

and \(\beta_n(r)\) are the dimensionless “anisotropy parameters” describing relative contributions of each harmonic order (\(\beta_0(r) = 1\) by definition). In particular:

\(\beta_2 = 2\) for the \(\cos^2 \theta\) (∥) angular distribution,

\(\beta_2 = 0\) for the isotropic distribution,

\(\beta_2 = -1\) for the \(\sin^2 \theta\) (⟂) angular distribution.

The radial intensity distribution alone for data with arbitrary angular variations can be obtained by using weight='sin' and order=0.

Parameters:window (int) – window size in pixels for radial averaging of \(\beta\). Since anisotropy parameters are non-linear, the central moving average is applied to the harmonics (which are linear), and then \(\beta\) is calculated from them. In case of well separated peaks, setting window to the peak width will result in \(\beta\) values at peak centers equal to total peak anisotropies (beware of the background, however).
Returns:Ibeta – radial intensity distribution (0-th term) and radial dependences of anisotropy parameters (other terms)
Return type:(# terms) × (rmax + 1) numpy array
rIbeta(window=1)[source]

Same as Ibeta(), but prepended with the radii row.

image(IM)[source]

Analyze an image.

This method can be also conveniently accessed by “calling” the object itself:

distr = Distributions(...)
Ibeta = distr(IM).Ibeta()
Parameters:IM (m × n numpy array) – the image to analyze
Returns:results – the object with analysis results, from which various distributions can be retrieved, see Results
Return type:Distributions.Results object
abel.tools.vmi.harmonics(IM, origin='cc', rmax='MIN', order=2, **kwargs)[source]

Convenience function to calculate harmonic distributions for a single image. Equivalent to Distributions(...).image(IM).harmonics().

Notice that this function does not cache intermediate calculations, so using it to process multiple images is several times slower than through a Distributions object.

abel.tools.vmi.rharmonics(IM, origin='cc', rmax='MIN', order=2, **kwargs)[source]

Same as harmonics(), but prepended with the radii row.

abel.tools.vmi.Ibeta(IM, origin='cc', rmax='MIN', order=2, window=1, **kwargs)[source]

Convenience function to calculate radial intensity and anisotropy distributions for a single image. Equivalent to Distributions(...).image(IM).Ibeta(window).

Notice that this function does not cache intermediate calculations, so using it to process multiple images is several times slower than through a Distributions object.

abel.tools.vmi.rIbeta(IM, origin='cc', rmax='MIN', order=2, window=1, **kwargs)[source]

Same as Ibeta(), but prepended with the radii row.

abel.benchmark module

class abel.benchmark.Timent(skip=0, repeat=1, duration=0.0)[source]

Bases: object

Helper class for measuring execution times.

The constructor only initializes the timing-procedure parameters. Use the time() method to run it for particular functions.

Parameters:
  • skip (int) – number of “warm-up” iterations to perform before the measurements. Can be specified as a negative number, then abs(skip) “warm-up” iterations are performed, but if this took more than duration seconds, they are accounted towards the measured iterations.
  • repeat (int) – minimal number of measured iterations to perform. Must be positive.
  • duration (float) – minimal duration (in seconds) of the measurements.
time(func, *args, **kwargs)[source]

Repeatedly executes a function at least repeat times and for at least duration seconds (see above), then returns the average time per iteration. The actual number of measured iterations can be retrieved from Timent.count.

Parameters:
  • func (callable) – function to execute
  • *args, **kwargs (any, optional) – parameters to pass to func
Returns:

average function execution time

Return type:

float

Notes

The measurements overhead can be estimated by executing

Timent(...).time(lambda: None)

with a sufficiently large number of iterations (to avoid rounding errors due to the finite timer precision). In 2018, this overhead was on the order of 100 ns per iteration.

class abel.benchmark.AbelTiming(n=[301, 501], select='all', repeat=1, t_min=0.1, t_max=inf, verbose=True)[source]

Bases: object

Benchmark performance of different Abel implementations (basis generation, forward and inverse transforms, as applicable).

Parameters:
  • n (int or sequence of int) – array size(s) for the benchmark (assuming 2D square arrays (nn))

  • select (str or sequence of str) – methods to benchmark. Use 'all' (default) for all available or choose any combination of individual methods:

    select=['basex', 'direct_C', 'direct_Python', 'hansenlaw',
            'linbasex', 'onion_bordas, 'onion_peeling', 'two_point',
            'three_point']
    
  • repeat (int) – repeat each benchmark at least this number of times to get the average values

  • t_min (float) – repeat each benchmark for at least this number of seconds to get the average values

  • t_max (float) – do not benchmark methods at array sizes when this is expected to take longer than this number of seconds. Notice that benchmarks for the smallest size from n are always run and that the estimations can be off by a factor of 2 or so.

  • verbose (boolean) – determines whether benchmark progress should be reported (to stderr)

n

array sizes from the parameter n, sorted in ascending order

Type:list of int
bs, fabel, iabel

benchmark results — dictionaries for

bs
basis-set generation
fabel
forward Abel transform
iabel
inverse Abel transform

with methods as keys and lists of timings in milliseconds as entries. Timings correspond to array sizes in AbelTiming.n; for skipped benchmarks (see t_max) they are np.nan.

Type:dict of list of float

Notes

The results can be output in a nice format by simply print(AbelTiming(...)).

Keep in mind that most methods have \(O(n^2)\) memory and \(O(n^3)\) time complexity, so going from n = 501 to n = 5001 would require about 100 times more memory and take about 1000 times longer.

class abel.benchmark.DistributionsTiming(n=[301, 501], shape='half', rmax='MIN', order=2, weight=['none', 'sin', 'sin+array'], method='all', repeat=1, t_min=0.1)[source]

Bases: object

Benchmark performance of different VMI distributions implementations.

Parameters:
  • n (int or sequence of int) – array size(s) for the benchmark (assuming full images to be 2D square arrays (nn))

  • shape (str) – image shape:

    'Q':

    one quadrant ((n + 1)/2, (n + 1)/2)

    'half' (default):

    half image (n, (n + 1)/2), vertically centered

    'full':

    full image (n, n), centered

  • rmax (str or sequence of str) – 'MIN' (default) and/or 'all', see rmax in abel.tools.vmi.Distributions

  • order (int) – highest order in the angular distributions. Even number ≥ 0.

  • weight (str or sequence of str) – weighting to test. Use 'all' for all available or choose any combination of individual types:

    weight=['none', 'sin', 'array', 'sin+array']
    
  • method (str or sequence of str) – methods to benchmark. Use 'all' (default) for all available or choose any combination of individual methods:

    method=['nearest', 'linear', 'remap']
    
  • repeat (int) – repeat each benchmark at least this number of times to get the average values

  • t_min (float) – repeat each benchmark for at least this number of seconds to get the average values

n

array sizes from the parameter n

Type:list of int
results

benchmark results — multi-level dictionary, in which results[method][rmax][weight] is the list of timings in milliseconds corresponding to array sizes in DistributionsTiming.n. Each timing is a tuple (t1, t) with t1 corresponding to single-image (non-cached) performance, and t corresponding to batch (cached) performance.

Type:dict of dict of dict of list of tuple of float

Notes

The results can be output in a nice format by simply print(DistributionsTiming(...)).

abel.benchmark.is_symmetric(arr, i_sym=True, j_sym=True)[source]

Takes in an array of shape (n, m) and check if it is symmetric

Parameters:
  • arr (1D or 2D array)
  • i_sym (array) – symmetric with respect to the 1st axis
  • j_sym (array) – symmetric with respect to the 2nd axis
Returns:

  • a binary array with the symmetry condition for the corresponding quadrants.
  • The globa

Notes

If both i_sym = True and j_sym = True, the input array is checked for polar symmetry.

See issue #34 comment for the defintion of a center of the image.

abel.benchmark.absolute_ratio_benchmark(analytical, recon, kind='inverse')[source]

Check the absolute ratio between an analytical function and the result of a inverse Abel reconstruction.

Parameters:
  • analytical (one of the classes from analytical, initialized)
  • recon (1D ndarray) – a reconstruction (i.e. inverse abel) given by some PyAbel implementation