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. Also accessible as abel.Transform.

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).

daun
the deconvolution method with Tikhonov regularization of Daun et al. and its extensions.

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.

Note: in PyAbel ≤0.8.4 the intensity distribution was off by a factor of π, please keep this in mind when comparing absolute intensities.
transform_options (dict) – Additional arguments passed to the individual transform functions. See the documentation for the individual transform method for options: basex, daun, direct, hansenlaw, linbasex, onion_bordas, onion_peeling, rbasex, three_point, two_point.
center_options (dict) – Additional arguments to be passed to the centering function, see abel.tools.center.center_image().
angular_integration_options (dict) – Additional arguments passed to the angular integration functions, see abel.tools.vmi.angular_integration_3D().
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:

symmetry_axis = 0 (vertical):

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

symmetry_axis = 1 (horizontal):

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

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).

daun *

Methods based on onion-peeling deconvolution using Tikhonov regularization described in

K. J. Daun, K. A. Thomson, F. Liu, G. J. Smallwood, “Deconvolution of axisymmetric flame properties using Tikhonov regularization”, Appl. Opt. 45, 4638–4646 (2006).

In addition to the original implicit step-functions basis (“onion peeling”) and the derivative regularization, linear and quadratic basis functions are implemented, as well as the $L_2$-norm Tikhonov regularization (like in basex) and non-negative least-squares solution.

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 (CH₂OH and CD₂OH)”, 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_3D().
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.
radial (numpy 1D array) – with method='linbasex': radial grid for Beta array
Beta (numpy 2D array) – with method='linbasex': 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]
projection (numpy 2D array) – with method='linbasex': 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.transform.set_basis_dir(basis_dir='', make=True)[source]

Changes the path to the directory for saving/loading cached basis sets that transform methods use by default.

Parameters:

basis_dir (str or None) – absolute or relative path. Passing '' (default) resets to the system-dependent default path, see default_basis_dir(). For the current working directory (as in PyAbel up to v0.8.4), use '.'. To disable basis-set caching on disk, use None.
make (bool) – create the directory if it does not exist (default: yes)

Return type:

None

abel.transform.get_basis_dir(make=False)[source]

Gets the path to the directory for saving/loading cached basis sets that transform methods use by default. If not changed by set_basis_dir(), it depends on the operating system, see default_basis_dir().

Parameters:: make (bool) – create the directory if it does not exist (default: no)
Returns:: path – absolute or relative path if disk caching is enabled, otherwise None
Return type:: str or None

abel.transform.default_basis_dir()[source]

Gets full path to the system-dependent default directory for saving/loading cached basis sets:

Linux (and other Unix-like):: ~/.cache/PyAbel (or $XDG_CACHE_HOME/PyAbel if set)
macOS:: /Users/<user>/Library/Caches/PyAbel
Windows:: <user profile>\AppData\Local\PyAbel\cache (or %LOCALAPPDATA%\PyAbel\cache if set). See important notes below.

Parameters:: None
Returns:: path – full path to the system-dependent default basis-sets directory
Return type:: str

Notes for MS Windows users

Python installed from Microsoft Store redirects subdirectory creation in AppData\Local to a “private per-user, per-app location” AppData\Local\Packages\Python...\LocalCache\Local (see Using Python on Windows / Known Issues). However, if AppData\Local\PyAbel\ already exists (for example, was manually created not from Python), apparently it should be usable.
Old Windows versions (2000, XP, Server 2003) by default don’t set the LOCALAPPDATA environment variable, so PyAbel will create and use the AppData\Local subtree in the user profile folder. This is probably fine, but not how it should be. To use the standard location, please do
```
set LOCALAPPDATA=%USERPROFILE%\Local Settings\Application Data
```
before starting Python. Or permanently set it in “Environment Variables” from Windows “System Properties”.

abel.transform.basis_dir_cleanup(basis_dir='', method=None)[source]

Deletes saved basis sets.

Parameters:

basis_dir (str or None) – path to the directory with saved basis sets. Use '' for the default directory, see get_basis_dir(). (For convenience, None can be passed to do nothing.)
method (str or list of str or None) – transform methods for which basis sets should be deleted. Can be a single string (see the method parameter in Transform) or a list of strings. Use 'all' to delete basis sets for all methods. None does nothing.

Return type:

None

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 or None) – path to the directory for saving / loading the basis sets. Use '' for the default directory. If None, the basis set will not be loaded from or 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 or None) – path to the directory for saving / loading the basis sets. Use '' for the default directory. If None, the basis sets will not be loaded from or 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.

Return type:

None

abel.basex.basis_dir_cleanup(basis_dir='')[source]

Utility function.

Deletes basis sets saved on disk.

Parameters:: basis_dir (str or None) – absolute or relative path to the directory with saved basis sets. Use '' for the default directory. None does nothing.
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.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 or None) – 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. Use '' for the default directory. If None, the operator array will not be loaded from or 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 or None) – 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. Use '' for the default directory. If None, the operator array will not be loaded from or 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 or None) – 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. Use '' for the default directory. If None, the operator array will not be loaded from or 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. Use '' for the default directory. For None, do not load or 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
Return type:: None

abel.dasch.basis_dir_cleanup(method, basis_dir='')[source]

Utility function.

Deletes deconvolution operator arrays saved on disk.

Parameters:

method (str) – Abel transform method 'onion_peeling', 'three_point', or 'two_point'
basis_dir (str or None) – absolute or relative path to the directory with saved deconvolution operator arrays. Use '' for the default directory. None does nothing.

Return type:

None

abel.daun module

abel.daun.daun_transform(data, reg=0.0, degree=0, dr=1.0, direction='inverse', basis_dir=None, verbose=True)[source]

Forward and inverse Abel transforms based on onion-peeling deconvolution using Tikhonov regularization described in

K. J. Daun, K. A. Thomson, F. Liu, G. J. Smallwood, “Deconvolution of axisymmetric flame properties using Tikhonov regularization”, Appl. Opt. 45, 4638–4646 (2006).

with additional basis-function types and regularization methods (see description).

This function operates on the “right side” of an image, that it, just one half of a cylindrically symmetric image, with the axial pixels located in the 0-th column.

Parameters:

data (m × n numpy array) – the image to be transformed. data[:, 0] should correspond to the central column of the image.
reg (float or tuple or str) – regularization for the inverse transform:

strength:
same as ('diff', strength)

('diff', strength):
Tikhonov regularization using the first-order difference operator (first-derivative approximation), as described in the original article

('L2', strength):
Tikhonov regularization using the $L_2$ norm, like in BASEX

('L2c', strength):
same as ('L2', strength), but with an intensity correction applied to compensate the drop near the symmetry axis

'nonneg':
non-negative least-squares solution.

Warning: this regularization method is very slow, typically taking up to a minute for a megapixel image.
degree (int) – degree of basis-function polynomials:

0:
rectangular functions (step-function approximation), corresponding to “onion peeling” from the original article

1:
triangular functions (piecewise linear approximation)

2:
piecewise quadratic functions (smooth approximation)

3:
piecewise cubic functions (cubic-spline approximation)
dr (float) – pixel size in the radial direction. This only affects the absolute scaling of the transformed image.
direction (str: 'forward' or 'inverse') – type of Abel transform to be performed
basis_dir (str, optional) – path to the directory for saving / loading the basis set (potentially useful only for degree = 3 and transform without regularization; time savings in other cases are small and might be negated by the disk-access overhead). Use '' for the default directory. If None (default), the basis set will not be loaded from or saved to disk.
verbose (bool) – determines whether progress report should be printed

Returns:

recon – the transformed (half) image

Return type:

m × n numpy array

abel.daun.get_bs_cached(n, degree=0, reg_type='diff', strength=0, direction='inverse', basis_dir=None, verbose=False)[source]

Internal function.

Gets the basis set and calculates the necessary transform matrix (notice that inverse direction with 'nonneg' regularization, as well as with strength = 0 for degree ≠ 3, gives the forward (triangular) matrix, to be used in solvers).

Parameters:

n (int) – half-width of the image in pixels, must include the axial pixel
degree (int) – polynomial degree for basis functions (0–3)
reg_type (None or str) – regularization type (None, 'diff', 'L2', 'L2c', 'nonneg')
strength (float) – Tikhonov regularization parameter (for reg_type = 'diff' and 'L2'/'L2c', ignored otherwise)
direction (str: 'forward' or 'inverse') – type of Abel transform to be performed
basis_dir (str or None) – path to the directory for saving / loading the basis set. Use '' for the default directory. If None, the basis sets will not be loaded from or saved to disk.
verbose (bool) – print some debug information

Returns:

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

Return type:

n × n numpy array

abel.daun.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 set
'inverse': only inverse transform

Return type:

None

abel.daun.basis_dir_cleanup(basis_dir='')[source]

Utility function.

Deletes basis sets saved on disk.

Parameters:: basis_dir (str or None) – absolute or relative path to the directory with saved basis sets. Use '' for the default directory. None does nothing.
Return type:: None

abel.direct module

abel.direct.direct_transform(fr, dr=None, r=None, direction='inverse', derivative=<function gradient>, int_func=<function trapezoid>, 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/np.trapezoid, 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.

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.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 a single image quadrant in the upper right orientation (Q0). Is not applicable to images with odd Legendre orders.

Parameters not described below are passed directly to linbasex_transform_full().

Parameters:

IM (numpy 2D array) – upper right quadrant of the image data, must be square in shape
return_Beta (bool) – in addition to the transformed image, return the radial, Beta and projections arrays
dr (any) – dummy variable for call compatibility with the other methods

Returns:

inv_IM (numpy 2D array) – upper right quadrant of the inverse Abel transformed image
radial (numpy 1D array) – (only if return_Beta = True) radii of each Newton sphere
Beta (numpy 2D array) – (only 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 radial is the speed distribution

Beta[1] vs radial is the anisotropy of each Newton sphere
projections (numpy 2D array) – (only if return_Beta = True) 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=<deprecated>, 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).

Lin-Basex 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.

This function operates on the whole image.

Parameters:

IM (numpy 2D array) – image data must have square shape of odd size
basis_dir (str or None) – path to the directory for saving / loading the basis sets. Use '' for the default directory. If None (default), the basis set will not be loaded from or saved to disk.
proj_angles (list of float) – 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 of int) – 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$ halfwidth equal to smoothing.
rcond (float) – (default 0.0005) scipy.linalg.lstsq() fit conditioning value. Use 0 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.
threshold (float) – threshold for normalization of higher-order Newton spheres (default 0.2): if Beta[0] < threshold, the associated Beta[j] for all j ⩾ 1 are set to zero
clip (int) – clip first vectors (smallest Newton spheres) to avoid singularities (default 0)
norm_range (tuple of int) – (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.
direction (str) – Abel transform direction. Only “inverse” is implemented.
verbose (bool) – print information about processing (normally used for debugging)

Returns:

inv_IM (numpy 2D array) – inverse Abel transformed image
radial (numpy 1D array) – radii of each Newton sphere
Beta (numpy 2D array) – 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 radial is the speed distribution

Beta[1] vs radial is the anisotropy of each Newton sphere
projections (numpy 2D array) – 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.

Warning

This function is deprecated and will be remove in the future. See issue #356.

For integrating the speed distribution and averaging the anisotropy, please use mean_beta().

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.mean_beta(radial, Beta, regions)[source]

Integrate normalized intensity (Beta[0]) and perform intensity-weighted averaging of anisotropy (Beta[1:]) over ranges of Newton spheres.

Parameters:

radial (numpy 1D array) – radii of Newton spheres
Beta (numpy 2D array) – speed and anisotropy distribution from linbasex_transform_full()
regions (list of tuple of int) – radial ranges [(min0, max0), (min1, max1), …]. Note that inclusion of regions where Beta[0] is below threshold set in linbasex_transform_full() will bias the mean anisotropies towards zero.

Returns:

Beta_mean – overall intensity (Beta_mean[0]) and mean anisotropy values (Beta_mean[1:]) in each region

Return type:

2D 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 or None) – path to the directory for saving / loading the basis. Use '' for the default directory. If None, the basis set will not be loaded from or saved to disk.
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
Return type:: None

abel.linbasex.basis_dir_cleanup(basis_dir='')[source]

Utility function.

Deletes basis sets saved on disk.

Parameters:: basis_dir (str or None) – relative or absolute path to the directory with saved basis sets. Use '' for the default directory. None does nothing.
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.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]

rBasex Abel transform for velocity-mapping images, operating in polar coordinates.

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). Use '' for the default directory. 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. Use '' for the default directory. If None, the basis set will not be loaded from or 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.

Return type:

None

abel.rbasex.basis_dir_cleanup(basis_dir='')[source]

Utility function.

Deletes basis sets saved on disk.

Parameters:: basis_dir (str or None) – absolute or relative path to the directory with saved basis sets. Use '' for the default directory. None does nothing.
Return type:: None

Image processing tools

abel.tools.analytical module

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

Bases: object

Base class for functions that have a known Abel transform (see GaussianAnalytical for a concrete example). Every such class should expose the following public attributes:

r

vector of positions along the $r$ axis

Type:: numpy array

func

the values of the original function (same shape as r for 1D functions, or same row size as r for 2D images)

Type:: numpy array

abel

the values of the Abel transform (same shape as func)

Type:: numpy array

mask_valid

mask (same shape as func) where the function is well smoothed/well behaved (no known artefacts in the inverse Abel reconstuction), typically excluding the origin, the domain boundaries, and function discontinuities, that can be used for unit testing.

Type:: numpy array

Parameters:

n (int) – number of points along the $r$ axis (saved to attribute n)
r_max (float) – maximum $r$ value (saved to attribute r_max)
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]

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

Bases: BaseAnalytical

Define a 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 for 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.

Parameters:

r (1D array) – array of positions along the r axis. Must start with 0.
A0 (float or 1D array) – height of the step. If 1D array, the height can be variable along the z axis
r0, r1 (float) – positions of the step along the r axis

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: 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
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₁ (r − r_0) + c₂ (r − r_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)
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]

class abel.tools.analytical.PiecewisePolynomial(n, r_max, ranges, symmetric=True)[source]

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

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

Bases: 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
sigma (float) – 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)
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]

class abel.tools.analytical.TransformPair(n, profile=5)[source]

Bases: 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=None, temperature=200)[source]

Bases: BaseAnalytical

Sample images, made up of Gaussian functions (or cubic splines, for 'O2').

Parameters:

n (integer) – image size n rows × n cols (must be odd for most purposes; even n values would result in half-pixel centering)
name (str) –

'Dribinski'
Sample test image used in the BASEX paper Rev. Sci. Instrum. 73, 2634 (2002), intensity function Eq. (16) (there are some missing negative exponents in the publication).

9 Gaussian peaks with alternating anisotropies (β = −1, 0, +2), plus 1 wide background Gaussian. Peak amplitudes are designed to produce comparable heights in the speed distribution, thus the peaks at small radii appear very intense in the image and its Abel transform.

'Gaussian'
Isotropic 2D Gaussian $\exp(-r^2 / \textbf{sigma}^2)$.

Its Abel transform is also a Gaussian with the same width: $\sqrt{\pi}\, \textbf{sigma} \exp(-r^2 / \textbf{sigma}^2)$.

'Gerber'
Artificial test image used in the lin-BASEX article Rev. Sci. Instrum. 84, 033101 (2013), Table I.

8 Gaussian peaks with various intensities and anisotropies up to 4th order (β₄).

'O2'
Synthetic image mimicking a velocity-map image of O⁺ from multiphoton photodissociation/ionization ${\rm O}_2 \xrightarrow{4h\nu} {\rm O} + {\rm O}^+ + e^-$ at $44\,444~\text{cm}^{-1}$ (225 nm)

Multiple peaks with various intensities and anisotropies; see, for example, J. Chem. Phys. 107, 2357 (1997).

'Ominus' or 'O-'
Simulate the photoelectron spectrum of O⁻ photodetachment ${}^3P_J \gets {}^2P_{3/2,1/2}$.

6 transitions, triplet neutral, and doublet anion.
sigma (float) – 1/e halfwidth of peaks in pixels, default values are: 2⋅r_max/180 for 'Dribinski', 2⋅r_max/500 for 'Ominus', r_max/3 for 'Gaussian', $\sqrt{2}$ (std. dev. = 1) for 'Gerber'.

For 'O2': HWHM of narrow peaks in pixels, default is 1.5 for any r_max.
temperature (float) – anion temperature in kelvins (default: 200) for 'Ominus': anion levels have Boltzmann population weight $(2J + 1) \exp[-hc \cdot 177.1~\text{cm}^{-1}/(k \cdot \textbf{temperature})]$

name

sample-image name

Type:: str

transform(tol=0.0048)[source]

Compute forward Abel transform of the image as an analytical Abel transform of its piecewise polynomial approximation (except 'Gaussian' and 'O2', which are computed exactly).

Parameters:: tol (float) – relative tolerance of the approximation (max. deviation divided by max. amplitude, default: 4.8e-3 ≲ 0.5%); the resulting Abel transform is somewhat more accurate
Returns:: abel – Abel-transformed image, also accessible as the abel attribute
Return type:: 2D numpy array

property image: Deprecated. Use func instead.

property abel: Abel transform of the image, computed (with default accuracy) only if necessary; see transform() for details.

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:
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
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.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(IM, 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(IM, 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 = N_angles × tol², where N_angles 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.
ref_angle (None or float) – Reference angle at which the radial correction function is renormalized to unity. If None, the angular average is used for renormalization.

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.

Returns:

radcorr – radial correction factors for angles

Return type:

numpy 1D array

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.guess_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.math.guss_gaussian(x)[source]: Deprecated function. Use guess_gaussian() instead.

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.BasePolynomial[source]

Bases: object

Abstract base class for polynomials. Implements multiplication and division by numbers. (Addition and subtraction of polynomials are not implemented because they are meaningful only for polynomials generated on the same grid. Use Piecewise... classes for sums of polynomials.)

func

values of the original function

Type:: numpy array

abel

values of the Abel transform

Type:: numpy array

copy()[source]: Return an independent copy.

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

Bases: BasePolynomial

Polynomial function and its Abel transform.

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_max ≲ max r (r_max might exceed maximal r, but usually by < 1 pixel; negative r_min or r_max are allowed for convenience but are interpreted as 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₁ (r − r_0) + c₂ (r − r_0)² + …
s (float, optional) – r stretching factor (around r_0): the polynomial is defined as c₀ + c₁ ((r − r_0)/s) + c₂ ((r − r_0)/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: BasePolynomial

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

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

p

Polynomial objects corresponding to each piece

Type:: list of Polynomial

copy()[source]: Make an independent copy.

class abel.tools.polynomial.SPolynomial(r, cos, r_min, r_max, c, r_0=0.0, s=1.0)[source]

Bases: BasePolynomial

Bivariate polynomial function $\sum_{mn} c_{mn} r^m \cos^n\theta$ in spherical coordinates and its Abel transform.

Parameters:

r, cos (numpy array) – r and cos θ values at which the function is generated; r must be non-negative. Arrays for generating a 2D image can be conveniently prepared by the rcos() function. On the other hand, the radial dependence alone (for a single cosine power) can be obtained by supplying a 1D r array and a cos array filled with ones.
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_max ≲ max r (r_max might exceed maximal r, but usually by < 1 pixel; negative r_min or r_max are allowed for convenience but are interpreted as 0)
c (2D numpy array) – polynomial coefficients for r and cos θ powers: c[m, n] is the coefficient for the $r^m \cos^n\theta$ term. This array can be conveniently constructed using Angular tools.
r_0 (float, optional) – r domain shift: the polynomial is defined in powers of (r − r_0) instead of r
s (float, optional) – r stretching factor (around r_0): the polynomial is defined in powers of (r − r_0)/s instead of r

class abel.tools.polynomial.PiecewiseSPolynomial(r, cos, ranges)[source]

Bases: BasePolynomial

Piecewise bivariate polynomial function (sum of SPolynomials) in spherical coordinates and its Abel transform.

Parameters:

r, cos (numpy array) – r and cos θ values at which the function is generated
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 SPolynomial conventions. All ranges are independent (may overlap and have gaps, may define polynomials of any degrees) and may include optional SPolynomial parameters (r_0, s).

abel.tools.polynomial.rcos(rows=None, cols=None, shape=None, origin=None)[source]

Create arrays with polar coordinates $r$ and $\cos\theta$: either from a pair of Cartesian arrays (rows, cols) with row and column values for each point or for a uniform grid with given dimensions and origin (shape, origin).

Parameters:

rows, cols (numpy array) – arrays with respectively row and column values for each point. Must have identical shapes (the output arrays will have the same shape), but might contain any values. For example, can be 2D arrays with integer pixel coordinates, or with floating-point numbers for sampling at subpixel points or on a distorted grid, or 1D arrays for sampling along some curve.
shape (tuple of int) – (rows, cols) – create output arrays of given shape, with values corresponding to a uniform pixel grid.
origin (tuple of float, optional) – position of the origin ($r = 0$) in the output array. By default, the center of the array is used (center of the middle pixel for odd-sized dimensions; even-sized dimensions will have a corresponding half-pixel shift).

Returns:

r (numpy array) – radii $r = \sqrt{\text{row}^2 + \text{col}^2}$ for each point
cos (numpy array) – cosines of the polar angle $\cos\theta = -\text{row}/r$ for each point (by convention, $\cos\theta = 0$ at $r = 0$)

class abel.tools.polynomial.Angular(c)[source]

Bases: object

Class helping to define angular dependences for SPolynomial and PiecewiseSPolynomial.

Supports arithmetic operations (addition, subtraction, multiplication of objects; multiplication and division by numbers) and outer product with radial coefficients (any list-like object). For example:

[3, 0, -1] * (Angular.cos(4) + Angular.sin(4) / 2)

represents $(3 - r^2)\big(\cos^4\theta + (\sin^4\theta) / 2\big)$, producing

[[ 1.5  0.  -3.  0.   4.5]
 [ 0.   0.  -0.  0.   0. ]
 [-0.5  0.   1.  0.  -1.5]]

which can be supplied as the coefficient matrix to SPolynomial. Likewise, a list of ranges for PiecewiseSPolynomial can be prepared as an outer product with a list of (r_min, r_max, coeffs) tuples (with optional other SPolynomial parameters), where 1D coeffs contain radial coefficients for a polynomial segment.

Parameters:: c (float or iterable of float) – list of coefficients: Angular([c₀, c₁, c₂, ...]) means $c_0 \cos^0\theta + c_1 \cos^1\theta + c_2 \cos^2\theta + \dots$; Angular(a) represents the isotropic distribution a⋅cos⁰ θ

c

coefficients for $\cos^n\theta$ powers, passed at instantiation directly (see above) or converted from other representations by the methods below.

Type:: numpy array

classmethod cos(n)[source]: Cosine power: Angular.cos(n) means $\cos^n\theta$.

classmethod sin(n)[source]: Sine power: Angular.sin(n) means $\sin^n\theta$ (n must be even).

classmethod cossin(m, n)[source]: Product of cosine and sine powers: Angular.cossin(m, n) means $\cos^m\theta \cdot \sin^n\theta$ (n must be even).

classmethod legendre(c)[source]

Weighted sum of Legendre polynomials in cos θ: Angular.legendre([c₀, c₁, c₂, ...]) means $c_0 P_0(\cos\theta) + c_1 P_1(\cos\theta) + c_2 P_2(\cos\theta) + \dots$

This method is intended to be called like

Angular.legendre([1, β₁, β₂, ...])

where $\beta_i$ are so-called anisotropy parameters. However, if you really need a single polynomial $P_n(\cos\theta)$, this can be easily achieved by

Angular.legendre([0] * n + [1])

class abel.tools.polynomial.ApproxGaussian(tol=0.0048)[source]

Bases: object

Piecewise quadratic approximation (non-negative and continuous but not exactly smooth) of the unit-amplitude, unit-SD Gaussian function $\exp(-r^2/2)$, equal to it at endpoints and midpoint of each piece. The forward Abel transform of this approximation will typically have a better relative accuracy than the approximation itself.

Parameters:

tol (float) – absolute approximation tolerance (maximal deviation). Some reference values yielding the best accuracy for certain number of segments:

tol	Better than	Segments
3.7e-2	5%	3
1.4e-2	2%	5
4.8e-3	0.5%	7 (default)
0.86e-3	0.1%	13
0.99e-4	0.01%	27
0.95e-5	0.001%	59

ranges

list of parameters (r_min, r_max, [c₀, c₁, c₂], r_0, s) that can be passed directly to PiecewisePolynomial or, after “multiplication” by Angular, to PiecewiseSPolynomial

Type:: lists of tuple

norm

the integral $\int_{-\infty}^{+\infty} f(r)\,dr$ for normalization (equals $\sqrt{2\pi}$ for the ideal Gaussian function, but slightly differs for the approximation)

Type:: float

scaled(A=1.0, r_0=0.0, sigma=1.0)[source]

Parameters for piecewise polynomials corresponding to the shifted and scaled Gaussian function $A \exp\big([(r - r_0)/\sigma]^2 / 2\big)$.

(Useful numbers: a Gaussian normalized to unit integral, that is the standard normal distribution, has $A = 1/\sqrt{2\pi}$; however, see norm above. A Gaussian with FWHM = 1 has $\sigma = 1/\sqrt{8\ln2}$.)

Parameters:

A (float) – amplitude
r_0 (float) – peak position
sigma (float) – standard deviation

Returns:

ranges – parameters for the piecewise polynomial approximating the shifted and scaled Gaussian

Return type:

list of tuple

abel.tools.polynomial.bspline(spl)[source]

Convert SciPy B-spline to PiecewisePolynomial parameters.

Parameters:: spl (tuple or BSpline or UnivariateSpline) – scipy.interpolate B-spline representation, such as splrep() results, BSpline object (result of make_interp_spline(), for example) or UnivariateSpline object
Returns:: ranges – list of parameters (r_min, r_max, coeffs, r_0) that can be passed directly to PiecewisePolynomial or, after “multiplication” by Angular, to PiecewiseSPolynomial
Return type:: list of tuple

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 < r \le 1\\I(r) &= \frac{1}{108}(128a_1 + a_{0.25}) + \frac{2}{27}r^2 (283a_{0.25} - 112a_1) + {}\\& \quad \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 < r \le 1\end{aligned}\end{align} \]

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} \]

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 < 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 < r \le 1\end{aligned}\end{align} \]

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

This profile has a small discontinuity at $r = 0.7$.

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 < r \le1\\I(r) &= 22.68862a_{0.7} - 14.811667a_1 + {}\\ &\quad (217.557a_{0.7} - 196.30083a_1)r^2 + 155.56r^2\ln\frac{1 + a_1}{0.7 + a_{0.7}} + {}\\ &\quad 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 + {}\\ &\quad r^2(155.56 + 55.5525r^2) \ln\frac{1 + a_1}{r} & 0.7 < r \le 1\end{aligned}\end{align} \]

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.

Note

This profile is discontinuous (and its projection is not smooth) at $r = 1$, which can cause different problems in different methods, in particular, depending on their assumptions where the singularity is located within the last pixel.

\[ \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} \]

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} \]

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} \]

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() with 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---    ---o----    ---o---    ---o---
    Q2 | Q3    Q2 | Q2     Q1 | Q0    Q1 | Q1

Return type:

np.array

abel.tools.vmi module

abel.tools.vmi.radial_intensity(kind, IM, origin=None, dr=1, dt=None)[source]

Calculate the one-dimensional radial intensity profile by angular integration or averaging of the image, treated either as a two-dimensional distribution or as a central slice of a cylindrically symmetric three-dimensional distribution.

Parameters:

kind (str) – operation to perform:

'int2D':
integration in 2D over polar angles

'int3D':
integration in 3D over solid angles

'avg2D':
averaging in 2D over polar angles

'avg3D':
averaging in 3D over solid angles
IM (2D numpy.array) – the image data
origin (tuple of float or None) – image origin in the (row, column) format. If None, the geometric center of the image (rows // 2, cols // 2) is used.
dr (float) – radial grid spacing in pixels (default 1). dr=0.5 may reduce pixel granularity of the radial 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
intensity (1D numpy.array) – intensity profile as a function of the radial coordinate

abel.tools.vmi.angular_integration_2D(IM, origin=None, dr=1, dt=None)[source]

Angular integration of the image as a two-dimensional object.

Equivalent to radial_intensity('int2D', IM, origin, dr, dt).

abel.tools.vmi.angular_integration_3D(IM, origin=None, dr=1, dt=None)[source]

Angular integration of the three-dimensional cylindrically symmetric object represented by the image as its central slice. When applied to the inverse Abel transform of a velocity-mapping image, this yields the speed distribution.

Equivalent to radial_intensity('int3D', IM, origin, dr, dt).

abel.tools.vmi.average_radial_intensity_2D(IM, origin=None, dr=1, dt=None)[source]

Calculate the average radial intensity of the image as a two-dimensional object.

Equivalent to radial_intensity('avg2D', IM, origin, dr, dt).

abel.tools.vmi.average_radial_intensity_3D(IM, origin=None, dr=1, dt=None)[source]

Calculate the average radial intensity of the three-dimensional cylindrically symmetric object represented by the image as its central slice.

Equivalent to radial_intensity('avg3D', IM, origin, dr, dt).

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.

Warning

This function behaves incorrectly: misses a factor of π for 3D integration, with Jacobian=True, and for Jacobian=False returns the average (over polar angles) multiplied by 2π instead of integrating. It is currently deprecated and is provided only for backward compatibility, but will be removed in the future.

Please use radial_intensity(), angular_integration_2D() or angular_integration_3D().

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π.

Warning

This function is currently deprecated and is provided only for backward compatibility, but will be removed in the future.

Please use radial_intensity(), average_radial_intensity_2D() or average_radial_intensity_3D().

Parameters:

IM (2D numpy.array) – the image data
kwargs – additional keyword arguments to be passed to abel.tools.vmi.angular_integration()

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.average_radial_intensity_3D().

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 or int) –

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

int
radial step. Evaluates the intensity vs angle for the whole radial range (0, step), (step, 2*step), ..

Returns:

Beta (list of tuples) – (beta0, error_beta_fit0), (beta1, error_beta_fit1), … corresponding to the radial ranges
Amplitude (list of tuples) – (amp0, error_amp_fit0), (amp1, error_amp_fit1), … corresponding to the radial ranges
Rmidpt (list of float) – radial mid-point of each radial range
Intensity_vs_theta (list of 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 or None) – angular ranges over which to fit [(theta1, theta2), (theta3, theta4)]. Allows data to be excluded from fit; default (None) is to 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 origin band occurs at the electron affinity $EA = 3613$ cm⁻¹. 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/pixel², will give output energy units in eV. A value of $1.204\times 10^{-5}$ cm⁻¹/pixel² applies for “examples/data/O2-ANU1024.txt” (with Vrep = −2200 volts).
per_energy_scaling (bool) – sets the intensity Jacobian. If True, the returned intensities correspond to an “intensity per eV” or “intensity per cm⁻¹”. 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⁻¹) is $10^{7}/\lambda$ (nm) e.g. 1.0e7/454.5 for “examples/data/O2-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. -2200 (volts), for “examples/data/O2-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 ≤ n ≤ order).

(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 n ≠ m; 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 (n, n))

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 (n, n))
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 (t₁, t_∞) with t₁ 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

Return type:

a binary array with the symmetry condition for the corresponding quadrants

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