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 wholeimage 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 velocitymapping images in terms of 1D spherical functions by Gerber et al. (2013).
onion_bordas
the algorithm of Bordas and coworkers (1996), reimplemented by Rallis, Wells and coworkers (2014).
onion_peeling
the onion peeling deconvolution as described by Dasch (1992).
rbasex
a method similar to pBasex by Garcia et al. (2004) for velocitymapping images, but with analytical basis functions developed by Ryazanov (2012).
three_point
the threepoint transform of Dasch (1992).
two_point
the twopoint 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 counterclockwide 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 asuint8
oruint16
), and this does not always play well with the transorm algorithms. This should probably always be set toTrue
(default).verbose (bool) – determines whether noncritical 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 counterclockwise:
+++  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 gaussianlike functions, which themselves can be Abeltransformed 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 defacto standards in photoelectron/photoion imaging.
V. Dribinski, A. Ossadtchi, V. A. Mandelshtam, H. Reisler, “Reconstruction of Abeltransformable images: The Gaussian basisset expansion Abel transform method”, Rev. Sci. Instrum. 73, 2634–2642 (2002).
daun
*Methods based on onionpeeling 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 stepfunctions 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 nonnegative leastsquares solution.direct
This method attempts a direct integration of the Abeltransform 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
*Velocitymapping 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. LinBASEX 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 onedimensional 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. BenItzhak, E. Wells, “Incorporating real time velocity map image reconstruction into closedloop 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 codeline, 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, “Twodimensional 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 velocitymapping 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 overtoneinduced dissociation and isomerization dynamics of hydroxymethyl radical (CH_{2}OH and CD_{2}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, “Onedimensional tomography: a comparison of Abel, onionpeeling, 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 Abeltransformed image.
angular_integration (tuple) – (radialgrid, radialintensity) 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, recentered (optional) with an oddsize 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 arrayBeta (numpy 2D array) – with
method='linbasex'
: coefficients of Newtonsphere spherical harmonicsBeta[0] — the radial intensity variation
Beta[1] — the anisotropy parameter variation
… Beta[n] — higherorder terms up to legedre_orders = [0, …, n]
projection (numpy 2D array) – with
method='linbasex'
: radial projection profiles at angles proj_anglesdistr (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 systemdependent default path, seedefault_basis_dir()
. For the current working directory (as in PyAbel up to v0.8.4), use'.'
. To disable basisset caching on disk, useNone
.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, seedefault_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 systemdependent default directory for saving/loading cached basis sets:
 Linux (and other Unixlike):
~/.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 systemdependent default basissets directory
 Return type:
str
Notes for MS Windows users
Python installed from Microsoft Store redirects subdirectory creation in
AppData\Local
to a “private peruser, perapp location”AppData\Local\Packages\Python...\LocalCache\Local
(see Using Python on Windows / Known Issues). However, ifAppData\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 theAppData\Local
subtree in the user profile folder. This is probably fine, but not how it should be. To use the standard location, please doset 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, seeget_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 inTransform
) 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, useabel.Transform(image, method='basex', direction='inverse').transform
This BASEX implementation only works with images that have an oddinteger 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. IfNone
, 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 transformcalculation function
 Returns:
IM – the Abeltransformed 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 dropoff 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. IfNone
, 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 dropoff 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 twopoint deconvolution method.
C. J. Dasch, “Onedimensional tomography: a comparison of Abel, onionpeeling, and filtered backprojection methods”, Appl. Opt. 31, 1146–1152 (1992).
 Parameters:
IM (1D or 2D numpy array) – rightside halfimage (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. IfNone
, 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 halfimage
 Return type:
1D or 2D numpy array
 abel.dasch.three_point_transform(IM, basis_dir='', dr=1, direction='inverse', verbose=False)[source]¶
The threepoint deconvolution method.
C. J. Dasch, “Onedimensional tomography: a comparison of Abel, onionpeeling, and filtered backprojection methods”, Appl. Opt. 31, 1146–1152 (1992).
 Parameters:
IM (1D or 2D numpy array) – rightside halfimage (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. IfNone
, 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 halfimage
 Return type:
1D or 2D numpy array
 abel.dasch.onion_peeling_transform(IM, basis_dir='', dr=1, direction='inverse', verbose=False)[source]¶
The onionpeeling deconvolution method.
C. J. Dasch, “Onedimensional tomography: a comparison of Abel, onionpeeling, and filtered backprojection methods”, Appl. Opt. 31, 1146–1152 (1992).
 Parameters:
IM (1D or 2D numpy array) – rightside halfimage (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. IfNone
, 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 halfimage
 Return type:
1D or 2D numpy array
 abel.dasch.dasch_transform(IM, D)[source]¶
Inverse Abel transform using the given deconvolution Doperator 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
, ortwo_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. ForNone
, do not load or save the deconvolution operator arrayverbose (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 onionpeeling 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 basisfunction 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 0th 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 firstorder difference operator (firstderivative 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'
:nonnegative leastsquares solution.
Warning: this regularization method is very slow, typically taking up to a minute for a megapixel image.
degree (int) – degree of basisfunction polynomials:
 0:
rectangular functions (stepfunction 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 (cubicspline 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 performedbasis_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 diskaccess overhead). Use
''
for the default directory. IfNone
(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) – halfwidth 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 performedbasis_dir (str or None) – path to the directory for saving / loading the basis set. Use
''
for the default directory. IfNone
, 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.direct module¶
 abel.direct.direct_transform(fr, dr=None, r=None, direction='inverse', derivative=<function gradient>, int_func=<function trapz>, correction=True, backend='C', **kwargs)[source]¶
This algorithm performs a direct computation of the Abel transform integrals. When correction=False, the pixel at the lower bound of the integral (where y=r) is skipped, which causes a systematic error in the Abel transform. However, if correction=True is used, then an analytical transform transform is applied to this pixel, which makes the approximation that the function is linear across this pixel. With correction=True, the Direct method produces reasonable results.
The Direct method is implemented in both Python and a compiled C version using Cython, which is much faster. The implementation can be selected using the backend argument. If the Cbackend is not available, you must reinstall PyAbel with Numpy, Cython, and a Ccompiler 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 nonuniform data. However, this has not been regorously tested.
direction (string) – Determines if a forward or inverse Abel transform will be applied. can be ‘forward’ or ‘inverse’.
derivative (callable) – a function that can return the derivative of the fr array with respect to r. (only used in the inverse Abel transform).
int_func (function) – This function is used to complete the integration. It should resemble np.trapz, in that it must be callable using axis=, x=, and dx= keyword arguments.
correction (boolean) – If False the pixel where the weighting function has a singular value (where r==y) is simply skipped, causing a systematic underestimation 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 Clibrary 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.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 “rightside” image:
Abeltrans = abel.hansenlaw.hansenlaw_transform(image, direction='inverse')
Note
Image should be a rightside image, like this:
+ ++  *  *   *  *  < im  *  *  + o+  *  *   *  *   *  *  + ++
In accordance with all PyAbel methods the image origin
o
is defined to be midpixel 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) – Rightside halfimage (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
orinverse
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 halfimage
 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 sphereBeta (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 onedimensional projections of spherical functions”, Rev. Sci. Instrum. 84, 033101 (2013).
LinBasex models the image using a sum of Legendre polynomials at each radial pixel. As such, it should only be applied to situations that can be adequately represented by Legendre polynomials, i.e., images that feature sphericallike 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. IfNone
(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 singlephoton experiment there are only anisotropies up to second order. The interaction of 4 photons (fourwave 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 illposed. 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 higherorder 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 intensityweighted 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_dirEither, 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. IfNone
, 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.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. BenItzhak, E. Wells, “Incorporating real time velocity map image reconstruction into closedloop 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 midpixel. This corresponds to an evenwidth 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 onionpeeling transorm on a whole image, use
abel.Transform(image, method='onion_bordas').transform
 Parameters:
IM (1D or 2D numpy array) – rightside halfimage (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 halfimage
 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 velocitymapping 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 nonparameterized 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'
:nonparameterized method, finds the best (in the leastsquares sense) solution with nonnegative \(\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 (forodd=True
) up to rmax, but limited to the largest inputimage quadrant (or half)'unfold'
:like
'fold'
, but symmetrically “unfolded” to all 4 quadrants'full'
:all pixels with radii up to rmax
'fullunique'
:the unique part of
'full'
: Q0 (upper right) quadrant forodd=False
, right half forodd=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 diskaccess overhead). Use
''
for the default directory. IfNone
(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 performedreg (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. IfNone
, 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
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 asr
for 2D images) 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
 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 floatingpoint 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
Polynomial
s) 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 optionalPolynomial
parameterssymmetric (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
Abeltransform 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 halfpixel 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 linBASEX article Rev. Sci. Instrum. 84, 033101 (2013), Table I.
8 Gaussian peaks with various intensities and anisotropies up to 4th order (β_{4}).
'O2'
Synthetic image mimicking a velocitymap 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¶
sampleimage 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.8e3 ≲ 0.5%); the resulting Abel transform is somewhat more accurate
 Returns:
abel – Abeltransformed 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), or1
(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 toTrue
if the image will subsequently be Abeltransformed.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 asNone
, for example,method=(row, None)
ormethod=(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 midpoint 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 wholepixel 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), or1
(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), or1
(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), or1
(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), or1
(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 velocitymapimaging 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 highestresolution 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 wellcentered 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 intensityprofile.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 intensityprofile 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 subpixel precision.
origin (float tuple, str or None) – Precenter 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 signaltonoise 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 subpixel 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^{2}, 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
usesnumpy.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 polarcoordinate image, to remove the background intensity. This may improve the signaltonoise 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) – Rootmeansquare (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) – Midpoint 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
leastsquares 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 2dimensional 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
ToDo: implement smooth noiserobust differentiators for use on experimental data. http://www.holoborodko.com/pavel/numericalmethods/numericalderivative/smoothlownoisedifferentiators/
 abel.tools.math.gaussian(x, a, mu, sigma, c)[source]¶

\(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) – nonzero 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/imageinformationalongapolarcoordinatesystem
 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.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
 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 nonnegative 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 floatingpoint 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
Polynomial
s) 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 optionalPolynomial
parameters
 p¶
Polynomial
objects corresponding to each piece Type:
list of Polynomial
 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 nonnegative. 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 usingAngular
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
SPolynomial
s) 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 optionalSPolynomial
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 floatingpoint 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 oddsized dimensions; evensized dimensions will have a corresponding halfpixel 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
andPiecewiseSPolynomial
.Supports arithmetic operations (addition, subtraction, multiplication of objects; multiplication and division by numbers) and outer product with radial coefficients (any listlike 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 forPiecewiseSPolynomial
can be prepared as an outer product with a list of(r_min, r_max, coeffs)
tuples (with optional otherSPolynomial
parameters), where 1Dcoeffs
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 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 socalled 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 (nonnegative and continuous but not exactly smooth) of the unitamplitude, unitSD 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.7e2
5%
3
1.4e2
2%
5
4.8e3
0.5%
7 (default)
0.86e3
0.1%
13
0.99e4
0.01%
27
0.95e5
0.001%
59
 ranges¶
list of parameters
(r_min, r_max, [c₀, c₁, c₂], r_0, s)
that can be passed directly toPiecewisePolynomial
or, after “multiplication” byAngular
, toPiecewiseSPolynomial
 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 Bspline to
PiecewisePolynomial
parameters. Parameters:
spl (tuple or BSpline or UnivariateSpline) –
scipy.interpolate
Bspline representation, such assplrep()
results,BSpline
object (result ofmake_interp_spline()
, for example) orUnivariateSpline
object Returns:
ranges – list of parameters
(r_min, r_max, coeffs, r_0)
that can be passed directly toPiecewisePolynomial
or, after “multiplication” byAngular
, toPiecewiseSPolynomial
 Return type:
list of tuple
abel.tools.transform_pairs module¶
Analytical function Abeltransform 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.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 (topright)
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 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).
 Returns:
Q0, Q1, Q2, Q3 – shape: (
rows // 2 + rows % 2, cols // 2 + cols % 2
) all oriented in the same direction as Q0 ifreorient=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()
withreorient=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 sizesodd 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 onedimensional radial intensity profile by angular integration or averaging of the image, treated either as a twodimensional distribution or as a central slice of a cylindrically symmetric threedimensional 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 twodimensional 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 threedimensional cylindrically symmetric object represented by the image as its central slice. When applied to the inverse Abel transform of a velocitymapping 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 twodimensional 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 threedimensional 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 onedimensional 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 forJacobian=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()
orangular_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 toabel.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 callingabel.tools.vmi.angular_integration()
withJacobian=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()
oraverage_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 rangesr0_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 midpoint 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^21}{2}\) is a 2ndorder 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 voltageindependent 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_{2}^{−} photodetachment, the origin band occurs at the electron affinity \(EA = 3613\) cm^{−1}. Values for the ANU experiment are given below, see also examples/example_hansenlaw.py.
 Parameters:
radial (numpy 1D array) – radial coordinates.
intensity (numpy 1D array) – intensity values, at the radial array.
energy_cal_factor (float) – energy calibration factor that will convert radius squared into energy. The units affect the units of the output. e.g. inputs in eV/pixel^{2}, will give output energy units in eV. A value of \(1.204\times 10^{5}\) cm^{−1}/pixel^{2} applies for “examples/data/O2ANU1024.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^{−1}”. IfFalse
, 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 electronbindingenergy in units of energy_cal_factor. The conversion from wavelength (nm) to photon_energy (cm^{−1}) is \(10^{7}/\lambda\) (nm) e.g.
1.0e7/454.5
for “examples/data/O2ANU1024.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/O2ANU1024.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 memoryefficient 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, n1]     center [m//2, 0] [m//2, n//2] [m//2, n1]     bottom/lower [m1, 0][m1, n//2][m1, n1]
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 upperright 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 asweights=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 likeI, beta2, beta4 = res.Ibeta()
. Python 3 users can also collect all \(\beta\) parameters asI, *beta = res.Ibeta()
for any order. Alternatively, transposing the results asIbeta = res.Ibeta().T
allows accessing all terms \(\big(I(r), \beta_2(r), \beta_4(r), \dots\big)\) at particular radius r asIbeta[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 byorders
(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
 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
 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'
andorder=0
. Parameters:
window (int) – window size in pixels for radial averaging of \(\beta\). Since anisotropy parameters are nonlinear, 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 (0th term) and radial dependences of anisotropy parameters (other terms)
 Return type:
(# terms) × (rmax + 1) numpy array
 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.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 timingprocedure parameters. Use the
time()
method to run it for particular functions. Parameters:
skip (int) – number of “warmup” iterations to perform before the measurements. Can be specified as a negative number, then
abs(skip)
“warmup” 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
basisset 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 arenp.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 inabel.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 — multilevel dictionary, in which
results[method][rmax][weight]
is the list of timings in milliseconds corresponding to array sizes inDistributionsTiming.n
. Each timing is a tuple (t_{1}, t_{∞}) with t_{1} corresponding to singleimage (noncached) 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