pirt.interp - Interpolation functions

The interp module implements several functions for interpolation, implemented in Numba.

Interpolation by warping and projecting

pirt.warp(data, samples, order='linear', spline_type=0.0)

Interpolate (sample) data at the positions specified by samples (pixel coordinates).

Parameters:
data : array (float32 or float64)

Data to interpolate, can be 1D, 2D or 3D.

samples : tuple with numpy arrays

Each array specifies the sample position for one dimension (in x-y-z order). Can also be a stacked array as in skimage’s warp() (in z-y-x order).

order : integer or string

Order of interpolation. Can be 0:’nearest’, 1:’linear’, 2: ‘quadratic’, 3:’cubic’.

spline_type : float or string

Only for cubic interpolation. Specifies the type of spline. Can be ‘Basis’, ‘Hermite’, ‘Cardinal’, ‘Catmull-rom’, ‘Lagrange’, ‘Lanczos’, ‘quadratic’, or a float, specifying the tension parameter for the Cardinal spline. See the docs of get_cubic_spline_coefs() for more information.

Returns:
result : array

The result is of the same type as the data array, and of the same shape of the samples arrays, which can be of any shape. This flexibility makes this function suitable as a base function for higher level “sampling functions”.

Notes

The input data can have up to three dimensions. It can be of any dtype, but float32 or float64 is recommended in general.

An order of interpolation of 2 would naturally correspond to quadratic interpolation. However, due to its uneven coefficients it reques the same support (and speed) as a cubic interpolant. This implementation adds the two quadratic polynomials. Note that you can probably better use order=3 with a Catmull-Rom spline, which corresponds to the linear interpolation of the two quadratic polynomials.

It can be shown (see Thevenaz et al. 2000 “Interpolation Revisited”) that interpolation using a Cardinal spline is equivalent to interpolating-B-spline interpolation.

pirt.awarp(data, samples, order='linear', spline_type=0.0)

Interpolation in anisotropic array. Like warp(), but the samples are expressed in world coordimates.

pirt.project(data, samples)

Interpolate data to the positions specified by samples (pixel coordinates).

In contrast to warp(), the project() function applies forward deformations, moving the pixels/voxels to the given locations, rather than getting the pixel values from the given locations. Although this may feel closer to how one would like to think about deformations, this function is slower and has no options to determine the interpolation, because there is no interpolation, but splatting.

Parameters:
data : array (float32 or float64)

Data to interpolate, can be 1D, 2D or 3D.

samples : tuple with numpy arrays

Each array specifies the sample position for one dimension (in x-y-z order). In contrast to warp(), each array must have the same shape as data. Can also be a stacked array as in skimage’s warp() (in z-y-x order).

Returns:
result : array

The result is of the same type and shape as the data array.

pirt.aproject(data, samples)

Interpolation in anisotropic array. Like project(), but the samples are expressed in world coordimates.

Deforming

pirt.deform_backward(data, deltas, order=1, spline_type=0.0)

Interpolate data according to the deformations specified in deltas. Deltas should be a tuple of numpy arrays similar to ‘samples’ in the warp() function. They represent the relative sample positions expressed in world coordinates.

pirt.deform_forward(data, deltas)

Like deform_backward(), but applied to project (forward deformation).

High level functions

pirt.resize(data, new_shape, order=3, spline_type=0.0, prefilter=False, extra=False)

Resize the data to the specified new shape.

Parameters:
data : numpy array

The data to rezize.

new_shape : tuple

The new shape of the data (z-y-x order).

order : {0,1,3} or {‘nearest’, ‘linear’, ‘cubic’}

The interpolation order to use.

spline_type : float or string

Only for cubic interpolation. Specifies the type of spline. Can be ‘Basis’, ‘Hermite’, ‘Cardinal’, ‘Catmull-rom’, ‘Lagrange’, ‘Lanczos’, ‘quadratic’, or a float, specifying the tension parameter for the Cardinal spline. See the docs of get_cubic_spline_coefs() for more information.

prefilter : bool

Whether to apply (discrete Gaussian diffusion) anti-aliasing (when downampling). Default False.

extra : bool

Whether to extrapolate the data a bit. In this case, each datapoint is seen as spanning a space equal to the distance between the data points. This is the method used when you resize an image using e.g. paint.net or photoshop. If False, the first and last datapoint are exactly on top of the original first and last datapoint (like scipy.ndimage.zoom). Default False.

pirt.imresize(data, factor, order=3)

Convenience function to resize the image data (1D, 2D or 3D).

This function uses pirt.resize() with ‘prefilter’ and ‘extra’ set to True. This makes it more suitble for generic image resizing. Use pirt.resize() for more fine-grained control.

Parameters:
data : numpy array

The data to rezize.

new_shape : tuple

The new shape of the data (z-y-x order).

order : {0,1,3} or {‘nearest’, ‘linear’, ‘cubic’}

The interpolation order to use.

pirt.zoom(data, factor, order=3, spline_type=0.0, prefilter=False, extra=False)

Resize the data with the specified factor. The default behavior is the same as scipy.ndimage.zoom(), but three times faster.

Parameters:
data : numpy array

The data to rezize.

factor : scalar or tuple

The resize factor, optionally for each dimension (z-y-z order).

order : {0,1,3} or {‘nearest’, ‘linear’, ‘cubic’}

The interpolation order to use.

spline_type : float or string

Only for cubic interpolation. Specifies the type of spline. Can be ‘Basis’, ‘Hermite’, ‘Cardinal’, ‘Catmull-rom’, ‘Lagrange’, ‘Lanczos’, ‘quadratic’, or a float, specifying the tension parameter for the Cardinal spline. See the docs of get_cubic_spline_coefs() for more information.

prefilter : bool

Whether to apply (discrete Gaussian diffusion) anti-aliasing (when downampling). Default False.

extra : bool

Whether to extrapolate the data a bit. In this case, each datapoint is seen as spanning a space equal to the distance between the data points. This is the method used when you resize an image using e.g. paint.net or photoshop. If False, the first and last datapoint are exactly on top of the original first and last datapoint (like numpy.zoom). Default False.

pirt.imzoom(data, factor, order=3)

Convenience function to resize the image data (1D, 2D or 3D) with the specified factor.

This function uses pirt.interp.resize() with ‘prefilter’ and ‘extra’ set to True. This makes it more suitble for generic image resizing. Use pirt.resize() for more fine-grained control.

Parameters:
data : numpy array

The data to rezize.

factor : scalar or tuple

The resize factor, optionally for each dimension (z-y-x order).

order : {0,1,3} or {‘nearest’, ‘linear’, ‘cubic’}

The interpolation order to use.

Helper functions

pirt.get_cubic_spline_coefs(t, spline_type='Catmull-Rom')

Calculates the coefficients for a cubic spline and returns them as a tuple. t is the ratio between “left” point and “right” point on the lattice. If performance matters, consider using get_lut() instead.

spline_type can be (case insensitive):

<number between -1 and 1>: Gives a Cardinal spline with the specified number as its tension parameter. A Cardinal spline is a type of Hermite spline, where the tangents are calculated using points p0 and p3; the coefficients can be directly applied to p0 p1 p2 p3.

‘Catmull-Rom’ or ‘Cardinal0’: is a cardinal spline a tension of 0. An interesting note: if we would create two quadractic splines, that would fit a polynomial “f(t) = a*t*t + b*t + c” to the first three and last three knots respectively, and if we would then combine the two using linear interpolation, we would obtain a catmull-rom spline. I don’t know whether this is how the spline was designed, or if this is a “side-effect”.

‘B-pline’: Basis spline. Not an interpolating spline, but an approximating spline. Here too, the coeffs can be applied to p0-p3.

‘Hermite’: Gives the Hermite coefficients to be applied to p1 m1 p2 m2 (in this order), where p1 p2 are the closest knots and m1 m2 the tangents.

‘Lagrange’: The Lagrange spline or Lagrange polynomial is an interpolating spline. It is the same as Newton polynomials, but implemented in a different manner (wiki). The Newton implementation should be more efficient, but this implementation is much simpler, and is very similar to the B-spline implementation (only the coefficients are different!). Also, when for example interpolating an image, coefficients are reused often and can be precalculated to enhance speed.

‘Lanczos’: Lanczos interpolation (windowed sync funcion). Note that this is not really a spline, and that sum of the coefficients is not exactly 1. Often used in audio processing. The Lanczos spline is very similar to the Cardinal spline with a tension of -0.25.

‘quadratic’: Quadratic interpolation with a support of 4, essentially the addition of the two quadratic polynoms.

‘linear’: Linear interpolation. Effective support is 2. Added for completeness and testing.

‘nearest’: Nearest neighbour interpolation. Added for completeness and testing.

pirt.meshgrid(*args)

Meshgrid implementation for 1D, 2D, 3D and beyond.

meshgrid(nx, ny) will create meshgrids with the specified shapes.

meshgrid(ndarray) will create meshgrids corresponding to the shape of the given array (which must have 2 or 3 dimension).

meshgrid([2,3,4], [1,2,3], [4,1,2]) uses the supplied values to create the grids. These lists can also be numpy arrays.

Returns a tuple with the grids in x-y-z order, with arrays of type float32.

pirt.make_samples_absolute(samples)

Note: this function is intended for sampes that represent a deformation; the number of dimensions of each array should match the number of arrays.

Given a tuple of arrays that represent a relative deformation expressed in world coordinates (x,y,z order), returns a tuple of sample arrays that represents the absolute sample locations in pixel coordinates. It is assumed that the sampling of the data is the same as for the sample arrays. The origin property is not used.

This process can also be done with relative ease by adding a meshgrid and then using awarp() or aproject(). But by combining it in one step, this process becomes faster and saves memory. Note that the deform_*() functions use this function.

Sampling a 2D slice in a 3D volume

class pirt.interp.SliceInVolume(self, pos, normal=None, previous=None)

Defines a slice in a volume.

The two span vectors are in v and u respectively. In other words, vec1 is up, vec2 is right.

convert_local_to_global(self, p2d, p3d)

Convert local 2D points to global 3D points. UNTESTED