pytranskit.optrans.continuous package

base transform

class pytranskit.optrans.continuous.base.BaseMapper2D

Bases: BaseTransform

Base class for 2D optimal transport transform methods (e.g. CLOT, VOT2D).

Warning

This class should not be used directly. Use derived classes instead.

apply_forward_map(transport_map, sig1)

Appy forward transport map.

Parameters:

• transport_map (array, shape (2, height, width)) – Forward transport map.

• sig1 (array, shape (height, width)) – Signal to transform.

Returns:

sig0_recon – Reconstructed reference signal sig0.

Return type:

array, shape (height, width)

apply_inverse_map(transport_map, sig0)

Appy inverse transport map.

Parameters:

• transport_map (array, shape (2, height, width)) – Forward transport map. Inverse is computed in this function.

• sig0 (array, shape (height, width)) – Reference signal.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

class pytranskit.optrans.continuous.base.BaseTransform

Bases: object

Base class for optimal transport transform methods.

Warning

This class should not be used directly. Use derived classes instead.

apply_forward_map()

Placeholder for application of forward transport map. Subclasses should implement this method!

apply_inverse_map()

Placeholder for application of inverse transport map. Subclasses should implement this method!

forward()

Placeholder for forward transform. Subclasses should implement this method!

inverse()

Inverse transform.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

pytranskit.optrans.continuous.base.assert_equal_shape(a, b, names=None)

Throw a ValueError if a and b are not the same shape.

pytranskit.optrans.continuous.base.check_array(array, ndim=None, dtype='numeric', force_all_finite=True, force_strictly_positive=False)

Input validation on an array, list, or similar.

Parameters:

• array (object) – Input object to check/convert

• ndim (int or None (default=None)) – Number of dimensions that array should have. If None, the dimensions are not checked

• dtype (string, type, list of types or None (default='numeric')) – Data type of result. If None, the dtype of the input is preserved. If ‘numeric’, dtype is preserved unless array.dtype is object. If dtype is a list of types, conversion on the first type is only performed if the dtype of the input is not in the list.

• force_all_finite (boolean (default=True)) – Whether to raise an error on np.inf and np.nan in array

• force_strictly_positive (boolean (default=False)) – Whether to raise an error if any array elements are <= 0

Returns:

array_converted – The converted and validated array.

Return type:

object

pytranskit.optrans.continuous.base.griddata2d(img, f, order=1, fill_value=0.0)

Interpolate 2d scattered data

Parameters:

• img (2d array, shape (height, width)) – Image to interpolate.

• f (3d array, shape (2, height, width)) – Coordinates of at which img is defined. First dimension f[0] corresponds to y-coordinates, second dimension f[1] is x-coordinates.

• order (int (default=1)) – Order of the interpolation. Must be in the range 0-2.

• fill_value (float (default=0.)) – Value used for points outside the boundaries.

Returns:

Image interpolated on to regular gird defined by: x = 0:(width-1), y = 0:(height-1).

Return type:

out, 2d array, shape (height, width)

pytranskit.optrans.continuous.base.interp2d(img, f, order=1, fill_value=0.0)

Interpolation for 2D gridded data.

Parameters:

• img (2d array, shape (height, width)) – Image to interpolate. This function assumes a grid of sample points: x = 0:(width-1), y = 0:(height-1).

• f (3d array, shape (2, height, width)) – Coordinates of interpolated points. First dimension f[0] corresponds to y-coordinates, second dimension f[1] is x-coordinates.

• order (int (default=1)) – Order of the spline interpolation. Must be in the range 0-5.

• fill_value (float (default=0.)) – Value used for points outside the boundaries.

Returns:

Image interpolated at points defined by f.

Return type:

out, 2d array, shape (height, width)

cdt

class pytranskit.optrans.continuous.cdt.BaseTransform

Bases: object

Base class for optimal transport transform methods.

Warning

This class should not be used directly. Use derived classes instead.

apply_forward_map()

Placeholder for application of forward transport map. Subclasses should implement this method!

apply_inverse_map()

Placeholder for application of inverse transport map. Subclasses should implement this method!

forward()

Placeholder for forward transform. Subclasses should implement this method!

inverse()

Inverse transform.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

class pytranskit.optrans.continuous.cdt.CDT

Bases: BaseTransform

Cumulative Distribution Transform.

displacements_

Displacements u.

Type:

1d array

transport_map_

Transport map f.

Type:

1d array

References

[The cumulative distribution transform and linear pattern classification] (https://arxiv.org/abs/1507.05936)

apply_forward_map(transport_map, sig1)

Appy forward transport map.

Parameters:

• transport_map (1d array) – Forward transport map.

• sig1 (1d array) – Signal to transform.

Returns:

sig0_recon – Reconstructed reference signal sig0.

Return type:

1d array

apply_inverse_map(transport_map, sig0, x)

Apply inverse transport map.

Parameters:

• transport_map (1d array) – Forward transport map. Inverse is computed in this function.

• sig0 (1d array) – Reference signal.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

1d array

forward(x0, sig0, x1, sig1, rm_edge=False)

Forward transform.

Parameters:

• x0 (1d array) – Independent axis variable of reference signal (sig0).

• sig0 (1d array) – Reference signal.

• x1 (1d array) – Independent axis variable of the signal to transform (sig1).

• sig1 (1d array) – Signal to transform.

Returns:

• sig1_cdt (1d array) – CDT of input signal sig1 (new definition).

• sig1_hat (1d array) – old definition.

• xilde (1d array) – Independent axis variable in CDT space.

inverse(transport_map, sig0, x1)

Inverse transform.

Parameters:

• transport_map (1d array) – Forward transport map.

• sig0 (1d array) – Reference signal.

• x1 (1d array) – Independent axis variable of the signal to reconstruct.

Returns:

sig1_recon – Reconstructed signal.

Return type:

1d array

pytranskit.optrans.continuous.cdt.assert_equal_shape(a, b, names=None)

Throw a ValueError if a and b are not the same shape.

pytranskit.optrans.continuous.cdt.check_array(array, ndim=None, dtype='numeric', force_all_finite=True, force_strictly_positive=False)

Input validation on an array, list, or similar.

Parameters:

• array (object) – Input object to check/convert

• ndim (int or None (default=None)) – Number of dimensions that array should have. If None, the dimensions are not checked

• dtype (string, type, list of types or None (default='numeric')) – Data type of result. If None, the dtype of the input is preserved. If ‘numeric’, dtype is preserved unless array.dtype is object. If dtype is a list of types, conversion on the first type is only performed if the dtype of the input is not in the list.

• force_all_finite (boolean (default=True)) – Whether to raise an error on np.inf and np.nan in array

• force_strictly_positive (boolean (default=False)) – Whether to raise an error if any array elements are <= 0

Returns:

array_converted – The converted and validated array.

Return type:

object

pytranskit.optrans.continuous.cdt.interp(x, xp, fp, left=None, right=None, period=None)

One-dimensional linear interpolation for monotonically increasing sample points.

Returns the one-dimensional piecewise linear interpolant to a function with given discrete data points (xp, fp), evaluated at x.

Parameters:

• x (array_like) – The x-coordinates at which to evaluate the interpolated values.

• xp (1-D sequence of floats) – The x-coordinates of the data points, must be increasing if argument period is not specified. Otherwise, xp is internally sorted after normalizing the periodic boundaries with xp = xp % period.

• fp (1-D sequence of float or complex) – The y-coordinates of the data points, same length as xp.

• left (optional float or complex corresponding to fp) – Value to return for x < xp[0], default is fp[0].

• right (optional float or complex corresponding to fp) – Value to return for x > xp[-1], default is fp[-1].

• period (None or float, optional) –

  A period for the x-coordinates. This parameter allows the proper interpolation of angular x-coordinates. Parameters left and right are ignored if period is specified.

  New in version 1.10.0.

Returns:

y – The interpolated values, same shape as x.

Return type:

float or complex (corresponding to fp) or ndarray

Raises:

ValueError – If xp and fp have different length If xp or fp are not 1-D sequences If period == 0

See also

scipy.interpolate

Warning

The x-coordinate sequence is expected to be increasing, but this is not explicitly enforced. However, if the sequence xp is non-increasing, interpolation results are meaningless.

Note that, since NaN is unsortable, xp also cannot contain NaNs.

A simple check for xp being strictly increasing is:

np.all(np.diff(xp) > 0)

Examples

>>> xp = [1, 2, 3]
>>> fp = [3, 2, 0]
>>> np.interp(2.5, xp, fp)
1.0
>>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
array([3.  , 3.  , 2.5 , 0.56, 0.  ])
>>> UNDEF = -99.0
>>> np.interp(3.14, xp, fp, right=UNDEF)
-99.0

Plot an interpolant to the sine function:

>>> x = np.linspace(0, 2*np.pi, 10)
>>> y = np.sin(x)
>>> xvals = np.linspace(0, 2*np.pi, 50)
>>> yinterp = np.interp(xvals, x, y)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'o')
[<matplotlib.lines.Line2D object at 0x...>]
>>> plt.plot(xvals, yinterp, '-x')
[<matplotlib.lines.Line2D object at 0x...>]
>>> plt.show()

Interpolation with periodic x-coordinates:

>>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
>>> xp = [190, -190, 350, -350]
>>> fp = [5, 10, 3, 4]
>>> np.interp(x, xp, fp, period=360)
array([7.5 , 5.  , 8.75, 6.25, 3.  , 3.25, 3.5 , 3.75])

Complex interpolation:

>>> x = [1.5, 4.0]
>>> xp = [2,3,5]
>>> fp = [1.0j, 0, 2+3j]
>>> np.interp(x, xp, fp)
array([0.+1.j , 1.+1.5j])

pytranskit.optrans.continuous.cdt.signal_to_pdf(input, sigma=0.0, epsilon=1e-08, total=1.0)

Get the (smoothed) probability density function of a signal.

Performs the following operations: 1. Smooth sigma with a Gaussian filter 2. Normalize signal such that it sums to 1 3. Add epsilon to ensure signal is strictly positive 4. Re-normalize signal such that it sums to total

Parameters:

• input (ndarray) – Input array

• sigma (scalar) – Standard deviation for Gaussian kernel. The standard deviations of the Gaussian filter are given for each axis as a sequence, or as a single number, in which case it is equal for all axes.

• epsilon (scalar) – Offset to ensure that signal is strictly positive.

• total (scalar) – Value of the signal summation.

Returns:

pdf – Returned array of same shape as input

Return type:

ndarray

clot

class pytranskit.optrans.continuous.clot.CDT

Bases: BaseTransform

Cumulative Distribution Transform.

displacements_

Displacements u.

Type:

1d array

transport_map_

Transport map f.

Type:

1d array

References

[The cumulative distribution transform and linear pattern classification] (https://arxiv.org/abs/1507.05936)

apply_forward_map(transport_map, sig1)

Appy forward transport map.

Parameters:

• transport_map (1d array) – Forward transport map.

• sig1 (1d array) – Signal to transform.

Returns:

sig0_recon – Reconstructed reference signal sig0.

Return type:

1d array

apply_inverse_map(transport_map, sig0, x)

Apply inverse transport map.

Parameters:

• transport_map (1d array) – Forward transport map. Inverse is computed in this function.

• sig0 (1d array) – Reference signal.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

1d array

forward(x0, sig0, x1, sig1, rm_edge=False)

Forward transform.

Parameters:

• x0 (1d array) – Independent axis variable of reference signal (sig0).

• sig0 (1d array) – Reference signal.

• x1 (1d array) – Independent axis variable of the signal to transform (sig1).

• sig1 (1d array) – Signal to transform.

Returns:

• sig1_cdt (1d array) – CDT of input signal sig1 (new definition).

• sig1_hat (1d array) – old definition.

• xilde (1d array) – Independent axis variable in CDT space.

inverse(transport_map, sig0, x1)

Inverse transform.

Parameters:

• transport_map (1d array) – Forward transport map.

• sig0 (1d array) – Reference signal.

• x1 (1d array) – Independent axis variable of the signal to reconstruct.

Returns:

sig1_recon – Reconstructed signal.

Return type:

1d array

class pytranskit.optrans.continuous.clot.CLOT(lr=0.01, momentum=0.0, decay=0.0, max_iter=300, tol=0.001, verbose=0)

Bases: object

Continuous Linear Optimal Transport Transform.

This uses Nesterov’s accelerated gradient descent to remove the curl in the initial mapping.

Parameters:

• lr (float (default=0.01)) – Learning rate.

• momentum (float (default=0.)) – Nesterov accelerated gradient descent momentum.

• decay (float (default=0.)) – Learning rate decay over each update.

• max_iter (int (default=300)) – Maximum number of iterations.

• tol (float (default=0.001)) – Stop iterating when change in cost function is below this threshold.

• verbose (int (default=1)) – Verbosity during optimization. 0=no output, 1=print cost, 2=print all metrics.

displacements_

Displacements u. First index denotes direction: displacements_[0] is y-displacements, and displacements_[1] is x-displacements.

Type:

array, shape (2, height, width)

transport_map_

Transport map f. First index denotes direction: transport_map_[0] is y-map, and transport_map_[1] is x-map.

Type:

array, shape (2, height, width)

displacements_initial_

Initial displacements computed using the method by Haker et al.

Type:

array, shape (2, height, width)

transport_map_initial_

Initial transport map computed using the method by Haker et al.

Type:

array, shape (2, height, width)

cost_

Value of cost function at each iteration.

Type:

list of float

curl_

Curl at each iteration.

Type:

list of float

References

[A continuous linear optimal transport approach for pattern analysis in image datasets] (https://www.sciencedirect.com/science/article/pii/S0031320315003507) [Optimal mass transport for registration and warping] (https://link.springer.com/article/10.1023/B:VISI.0000036836.66311.97)

apply_forward_map(transport_map, sig1)

Apply forward transport map.

Parameters:

• transport_map (array, shape (2, height, width)) – Forward transport map.

• sig1 (array, shape (height, width)) – Signal to transform.

Returns:

sig0_recon – Reconstructed reference signal sig0.

Return type:

array, shape (height, width)

apply_inverse_map(transport_map, sig0)

Apply inverse transport map.

Parameters:

• transport_map (array, shape (2, height, width)) – Forward transport map. Inverse is computed in this function.

• sig0 (array, shape (height, width)) – Reference signal.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

forward(sig0, sig1)

Forward transform.

Parameters:

• sig0 (array, shape (height, width)) – Reference image.

• sig1 (array, shape (height, width)) – Signal to transform.

Returns:

lot – LOT transform of input image sig1. First index denotes direction: lot[0] is y-LOT, and lot[1] is x-LOT.

Return type:

array, shape (2, height, width)

pytranskit.optrans.continuous.clot.assert_equal_shape(a, b, names=None)

Throw a ValueError if a and b are not the same shape.

pytranskit.optrans.continuous.clot.check_array(array, ndim=None, dtype='numeric', force_all_finite=True, force_strictly_positive=False)

Input validation on an array, list, or similar.

Parameters:

• array (object) – Input object to check/convert

• ndim (int or None (default=None)) – Number of dimensions that array should have. If None, the dimensions are not checked

• dtype (string, type, list of types or None (default='numeric')) – Data type of result. If None, the dtype of the input is preserved. If ‘numeric’, dtype is preserved unless array.dtype is object. If dtype is a list of types, conversion on the first type is only performed if the dtype of the input is not in the list.

• force_all_finite (boolean (default=True)) – Whether to raise an error on np.inf and np.nan in array

• force_strictly_positive (boolean (default=False)) – Whether to raise an error if any array elements are <= 0

Returns:

array_converted – The converted and validated array.

Return type:

object

pytranskit.optrans.continuous.clot.dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False)

Return the Discrete Cosine Transform of arbitrary type sequence x.

Parameters:

• x (array_like) – The input array.

• type ({1, 2, 3, 4}, optional) – Type of the DCT (see Notes). Default type is 2.

• n (int, optional) – Length of the transform. If n < x.shape[axis], x is truncated. If n > x.shape[axis], x is zero-padded. The default results in n = x.shape[axis].

• axis (int, optional) – Axis along which the dct is computed; the default is over the last axis (i.e., axis=-1).

• norm ({None, 'ortho'}, optional) – Normalization mode (see Notes). Default is None.

• overwrite_x (bool, optional) – If True, the contents of x can be destroyed; the default is False.

Returns:

y – The transformed input array.

Return type:

ndarray of real

See also

idct

Inverse DCT

Notes

For a single dimension array x, dct(x, norm='ortho') is equal to MATLAB dct(x).

There are, theoretically, 8 types of the DCT, only the first 4 types are implemented in scipy. ‘The’ DCT generally refers to DCT type 2, and ‘the’ Inverse DCT generally refers to DCT type 3.

Type I

There are several definitions of the DCT-I; we use the following (for norm=None)

\[y_k = x_0 + (-1)^k x_{N-1} + 2 \sum_{n=1}^{N-2} x_n \cos\left( \frac{\pi k n}{N-1} \right)\]

If norm='ortho', x[0] and x[N-1] are multiplied by a scaling factor of \(\sqrt{2}\), and y[k] is multiplied by a scaling factor f

\[\begin{split}f = \begin{cases} \frac{1}{2}\sqrt{\frac{1}{N-1}} & \text{if }k=0\text{ or }N-1, \\ \frac{1}{2}\sqrt{\frac{2}{N-1}} & \text{otherwise} \end{cases}\end{split}\]

New in version 1.2.0: Orthonormalization in DCT-I.

Note

The DCT-I is only supported for input size > 1.

Type II

There are several definitions of the DCT-II; we use the following (for norm=None)

\[y_k = 2 \sum_{n=0}^{N-1} x_n \cos\left(\frac{\pi k(2n+1)}{2N} \right)\]

If norm='ortho', y[k] is multiplied by a scaling factor f

\[\begin{split}f = \begin{cases} \sqrt{\frac{1}{4N}} & \text{if }k=0, \\ \sqrt{\frac{1}{2N}} & \text{otherwise} \end{cases}\end{split}\]

which makes the corresponding matrix of coefficients orthonormal (O @ O.T = np.eye(N)).

Type III

There are several definitions, we use the following (for norm=None)

\[y_k = x_0 + 2 \sum_{n=1}^{N-1} x_n \cos\left(\frac{\pi(2k+1)n}{2N}\right)\]

or, for norm='ortho'

\[y_k = \frac{x_0}{\sqrt{N}} + \sqrt{\frac{2}{N}} \sum_{n=1}^{N-1} x_n \cos\left(\frac{\pi(2k+1)n}{2N}\right)\]

The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up to a factor 2N. The orthonormalized DCT-III is exactly the inverse of the orthonormalized DCT-II.

Type IV

There are several definitions of the DCT-IV; we use the following (for norm=None)

\[y_k = 2 \sum_{n=0}^{N-1} x_n \cos\left(\frac{\pi(2k+1)(2n+1)}{4N} \right)\]

If norm='ortho', y[k] is multiplied by a scaling factor f

\[f = \frac{1}{\sqrt{2N}}\]

New in version 1.2.0: Support for DCT-IV.

References

[1]

‘A Fast Cosine Transform in One and Two Dimensions’, by J. Makhoul, IEEE Transactions on acoustics, speech and signal processing vol. 28(1), pp. 27-34, :doi:`10.1109/TASSP.1980.1163351` (1980).

[2]

Wikipedia, “Discrete cosine transform”, https://en.wikipedia.org/wiki/Discrete_cosine_transform

Examples

The Type 1 DCT is equivalent to the FFT (though faster) for real, even-symmetrical inputs. The output is also real and even-symmetrical. Half of the FFT input is used to generate half of the FFT output:

>>> from scipy.fftpack import fft, dct
>>> fft(np.array([4., 3., 5., 10., 5., 3.])).real
array([ 30.,  -8.,   6.,  -2.,   6.,  -8.])
>>> dct(np.array([4., 3., 5., 10.]), 1)
array([ 30.,  -8.,   6.,  -2.])

pytranskit.optrans.continuous.clot.griddata2d(img, f, order=1, fill_value=0.0)

Interpolate 2d scattered data

Parameters:

• img (2d array, shape (height, width)) – Image to interpolate.

• f (3d array, shape (2, height, width)) – Coordinates of at which img is defined. First dimension f[0] corresponds to y-coordinates, second dimension f[1] is x-coordinates.

• order (int (default=1)) – Order of the interpolation. Must be in the range 0-2.

• fill_value (float (default=0.)) – Value used for points outside the boundaries.

Returns:

Image interpolated on to regular gird defined by: x = 0:(width-1), y = 0:(height-1).

Return type:

out, 2d array, shape (height, width)

pytranskit.optrans.continuous.clot.idct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False)

Return the Inverse Discrete Cosine Transform of an arbitrary type sequence.

Parameters:

• x (array_like) – The input array.

• type ({1, 2, 3, 4}, optional) – Type of the DCT (see Notes). Default type is 2.

• n (int, optional) – Length of the transform. If n < x.shape[axis], x is truncated. If n > x.shape[axis], x is zero-padded. The default results in n = x.shape[axis].

• axis (int, optional) – Axis along which the idct is computed; the default is over the last axis (i.e., axis=-1).

• norm ({None, 'ortho'}, optional) – Normalization mode (see Notes). Default is None.

• overwrite_x (bool, optional) – If True, the contents of x can be destroyed; the default is False.

Returns:

idct – The transformed input array.

Return type:

ndarray of real

See also

dct

Forward DCT

Notes

For a single dimension array x, idct(x, norm='ortho') is equal to MATLAB idct(x).

‘The’ IDCT is the IDCT of type 2, which is the same as DCT of type 3.

IDCT of type 1 is the DCT of type 1, IDCT of type 2 is the DCT of type 3, and IDCT of type 3 is the DCT of type 2. IDCT of type 4 is the DCT of type 4. For the definition of these types, see dct.

Examples

The Type 1 DCT is equivalent to the DFT for real, even-symmetrical inputs. The output is also real and even-symmetrical. Half of the IFFT input is used to generate half of the IFFT output:

>>> from scipy.fftpack import ifft, idct
>>> ifft(np.array([ 30.,  -8.,   6.,  -2.,   6.,  -8.])).real
array([  4.,   3.,   5.,  10.,   5.,   3.])
>>> idct(np.array([ 30.,  -8.,   6.,  -2.]), 1) / 6
array([  4.,   3.,   5.,  10.])

pytranskit.optrans.continuous.clot.interp2d(img, f, order=1, fill_value=0.0)

Interpolation for 2D gridded data.

Parameters:

• img (2d array, shape (height, width)) – Image to interpolate. This function assumes a grid of sample points: x = 0:(width-1), y = 0:(height-1).

• f (3d array, shape (2, height, width)) – Coordinates of interpolated points. First dimension f[0] corresponds to y-coordinates, second dimension f[1] is x-coordinates.

• order (int (default=1)) – Order of the spline interpolation. Must be in the range 0-5.

• fill_value (float (default=0.)) – Value used for points outside the boundaries.

Returns:

Image interpolated at points defined by f.

Return type:

out, 2d array, shape (height, width)

pytranskit.optrans.continuous.clot.signal_to_pdf(input, sigma=0.0, epsilon=1e-08, total=1.0)

Get the (smoothed) probability density function of a signal.

Performs the following operations: 1. Smooth sigma with a Gaussian filter 2. Normalize signal such that it sums to 1 3. Add epsilon to ensure signal is strictly positive 4. Re-normalize signal such that it sums to total

Parameters:

• input (ndarray) – Input array

• sigma (scalar) – Standard deviation for Gaussian kernel. The standard deviations of the Gaussian filter are given for each axis as a sequence, or as a single number, in which case it is equal for all axes.

• epsilon (scalar) – Offset to ensure that signal is strictly positive.

• total (scalar) – Value of the signal summation.

Returns:

pdf – Returned array of same shape as input

Return type:

ndarray

radoncdt

class pytranskit.optrans.continuous.radoncdt.BaseTransform

Bases: object

Base class for optimal transport transform methods.

Warning

This class should not be used directly. Use derived classes instead.

apply_forward_map()

Placeholder for application of forward transport map. Subclasses should implement this method!

apply_inverse_map()

Placeholder for application of inverse transport map. Subclasses should implement this method!

forward()

Placeholder for forward transform. Subclasses should implement this method!

inverse()

Inverse transform.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

class pytranskit.optrans.continuous.radoncdt.CDT

Bases: BaseTransform

Cumulative Distribution Transform.

displacements_

Displacements u.

Type:

1d array

transport_map_

Transport map f.

Type:

1d array

References

[The cumulative distribution transform and linear pattern classification] (https://arxiv.org/abs/1507.05936)

apply_forward_map(transport_map, sig1)

Appy forward transport map.

Parameters:

• transport_map (1d array) – Forward transport map.

• sig1 (1d array) – Signal to transform.

Returns:

sig0_recon – Reconstructed reference signal sig0.

Return type:

1d array

apply_inverse_map(transport_map, sig0, x)

Apply inverse transport map.

Parameters:

• transport_map (1d array) – Forward transport map. Inverse is computed in this function.

• sig0 (1d array) – Reference signal.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

1d array

forward(x0, sig0, x1, sig1, rm_edge=False)

Forward transform.

Parameters:

• x0 (1d array) – Independent axis variable of reference signal (sig0).

• sig0 (1d array) – Reference signal.

• x1 (1d array) – Independent axis variable of the signal to transform (sig1).

• sig1 (1d array) – Signal to transform.

Returns:

• sig1_cdt (1d array) – CDT of input signal sig1 (new definition).

• sig1_hat (1d array) – old definition.

• xilde (1d array) – Independent axis variable in CDT space.

inverse(transport_map, sig0, x1)

Inverse transform.

Parameters:

• transport_map (1d array) – Forward transport map.

• sig0 (1d array) – Reference signal.

• x1 (1d array) – Independent axis variable of the signal to reconstruct.

Returns:

sig1_recon – Reconstructed signal.

Return type:

1d array

class pytranskit.optrans.continuous.radoncdt.RadonCDT(theta=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179]))

Bases: BaseTransform

Radon Cumulative Distribution Transform.

Parameters:

theta (1d array (default=np.arange(180))) – Radon transform projection angles.

displacements_

Displacements u.

Type:

array, shape (t, len(theta))

transport_map_

Transport map f.

Type:

array, shape (t, len(theta))

References

[The Radon cumulative distribution transform and its application to image classification] (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4871726/)

apply_forward_map(transport_map, sig1)

Appy forward transport map.

Parameters:

• transport_map (array, shape (t, len(theta))) – Forward transport map.

• sig1 (2d array, shape (height, width)) – Signal to transform.

Returns:

sig0_recon – Reconstructed reference signal sig0.

Return type:

array, shape (height, width)

apply_inverse_map(transport_map, sig0, x_range)

Appy inverse transport map.

Parameters:

• transport_map (2d array, shape (t, len(theta))) – Forward transport map. Inverse is computed in this function.

• sig0 (array, shape (height, width)) – Reference signal.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

forward(x0_range, sig0, x1_range, sig1, rm_edge=False)

Forward transform.

Parameters:

• sig0 (array, shape (height, width)) – Reference image.

• sig1 (array, shape (height, width)) – Signal to transform.

Returns:

rcdt – Radon-CDT of input image sig1.

Return type:

array, shape (t, len(theta))

inverse(transport_map, sig0, x1_range)

Inverse transform.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

pytranskit.optrans.continuous.radoncdt.assert_equal_shape(a, b, names=None)

Throw a ValueError if a and b are not the same shape.

pytranskit.optrans.continuous.radoncdt.check_array(array, ndim=None, dtype='numeric', force_all_finite=True, force_strictly_positive=False)

Input validation on an array, list, or similar.

Parameters:

• array (object) – Input object to check/convert

• ndim (int or None (default=None)) – Number of dimensions that array should have. If None, the dimensions are not checked

• dtype (string, type, list of types or None (default='numeric')) – Data type of result. If None, the dtype of the input is preserved. If ‘numeric’, dtype is preserved unless array.dtype is object. If dtype is a list of types, conversion on the first type is only performed if the dtype of the input is not in the list.

• force_all_finite (boolean (default=True)) – Whether to raise an error on np.inf and np.nan in array

• force_strictly_positive (boolean (default=False)) – Whether to raise an error if any array elements are <= 0

Returns:

array_converted – The converted and validated array.

Return type:

object

pytranskit.optrans.continuous.radoncdt.iradon(radon_image, theta=None, output_size=None, filter_name='ramp', interpolation='linear', circle=True, preserve_range=True)

Inverse radon transform.

Reconstruct an image from the radon transform, using the filtered back projection algorithm.

Parameters:

• radon_image (array) – Image containing radon transform (sinogram). Each column of the image corresponds to a projection along a different angle. The tomography rotation axis should lie at the pixel index radon_image.shape[0] // 2 along the 0th dimension of radon_image.

• theta (array_like, optional) – Reconstruction angles (in degrees). Default: m angles evenly spaced between 0 and 180 (if the shape of radon_image is (N, M)).

• output_size (int, optional) – Number of rows and columns in the reconstruction.

• filter_name (str, optional) – Filter used in frequency domain filtering. Ramp filter used by default. Filters available: ramp, shepp-logan, cosine, hamming, hann. Assign None to use no filter.

• interpolation (str, optional) – Interpolation method used in reconstruction. Methods available: ‘linear’, ‘nearest’, and ‘cubic’ (‘cubic’ is slow).

• circle (boolean, optional) – Assume the reconstructed image is zero outside the inscribed circle. Also changes the default output_size to match the behaviour of radon called with circle=True.

• preserve_range (bool, optional) – Whether to keep the original range of values. Otherwise, the input image is converted according to the conventions of img_as_float. Also see https://scikit-image.org/docs/dev/user_guide/data_types.html

Returns:

• reconstructed (ndarray) – Reconstructed image. The rotation axis will be located in the pixel with indices (reconstructed.shape[0] // 2, reconstructed.shape[1] // 2).

• .. versionchanged:: 0.19 – In iradon, filter argument is deprecated in favor of filter_name.

References

[1]

AC Kak, M Slaney, “Principles of Computerized Tomographic Imaging”, IEEE Press 1988.

[2]

B.R. Ramesh, N. Srinivasa, K. Rajgopal, “An Algorithm for Computing the Discrete Radon Transform With Some Applications”, Proceedings of the Fourth IEEE Region 10 International Conference, TENCON ‘89, 1989

Notes

It applies the Fourier slice theorem to reconstruct an image by multiplying the frequency domain of the filter with the FFT of the projection data. This algorithm is called filtered back projection.

pytranskit.optrans.continuous.radoncdt.match_shape2d(a, b)

Crop array B such that it matches the shape of A.

Parameters:

• a (2d array) – Array of desired size.

• b (2d array) – Array to crop. Shape must be larger than (or equal to) the shape of array a.

Returns:

b_crop – Cropped version of b, with the same shape as a.

Return type:

2d array

pytranskit.optrans.continuous.radoncdt.radon(image, theta=None, circle=True, *, preserve_range=False)

Calculates the radon transform of an image given specified projection angles.

Parameters:

• image (array_like) – Input image. The rotation axis will be located in the pixel with indices (image.shape[0] // 2, image.shape[1] // 2).

• theta (array_like, optional) – Projection angles (in degrees). If None, the value is set to np.arange(180).

• circle (boolean, optional) – Assume image is zero outside the inscribed circle, making the width of each projection (the first dimension of the sinogram) equal to min(image.shape).

• preserve_range (bool, optional) – Whether to keep the original range of values. Otherwise, the input image is converted according to the conventions of img_as_float. Also see https://scikit-image.org/docs/dev/user_guide/data_types.html

Returns:

radon_image – Radon transform (sinogram). The tomography rotation axis will lie at the pixel index radon_image.shape[0] // 2 along the 0th dimension of radon_image.

Return type:

ndarray

References

[1]

AC Kak, M Slaney, “Principles of Computerized Tomographic Imaging”, IEEE Press 1988.

[2]

B.R. Ramesh, N. Srinivasa, K. Rajgopal, “An Algorithm for Computing the Discrete Radon Transform With Some Applications”, Proceedings of the Fourth IEEE Region 10 International Conference, TENCON ‘89, 1989

Notes

Based on code of Justin K. Romberg (https://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html)

pytranskit.optrans.continuous.radoncdt.signal_to_pdf(input, sigma=0.0, epsilon=1e-08, total=1.0)

Get the (smoothed) probability density function of a signal.

Performs the following operations: 1. Smooth sigma with a Gaussian filter 2. Normalize signal such that it sums to 1 3. Add epsilon to ensure signal is strictly positive 4. Re-normalize signal such that it sums to total

Parameters:

• input (ndarray) – Input array

• sigma (scalar) – Standard deviation for Gaussian kernel. The standard deviations of the Gaussian filter are given for each axis as a sequence, or as a single number, in which case it is equal for all axes.

• epsilon (scalar) – Offset to ensure that signal is strictly positive.

• total (scalar) – Value of the signal summation.

Returns:

pdf – Returned array of same shape as input

Return type:

ndarray

radoncdt3D

class pytranskit.optrans.continuous.radoncdt3D.BaseTransform

Bases: object

Base class for optimal transport transform methods.

Warning

This class should not be used directly. Use derived classes instead.

apply_forward_map()

Placeholder for application of forward transport map. Subclasses should implement this method!

apply_inverse_map()

Placeholder for application of inverse transport map. Subclasses should implement this method!

forward()

Placeholder for forward transform. Subclasses should implement this method!

inverse()

Inverse transform.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

class pytranskit.optrans.continuous.radoncdt3D.CDT

Bases: BaseTransform

Cumulative Distribution Transform.

displacements_

Displacements u.

Type:

1d array

transport_map_

Transport map f.

Type:

1d array

References

[The cumulative distribution transform and linear pattern classification] (https://arxiv.org/abs/1507.05936)

apply_forward_map(transport_map, sig1)

Appy forward transport map.

Parameters:

• transport_map (1d array) – Forward transport map.

• sig1 (1d array) – Signal to transform.

Returns:

sig0_recon – Reconstructed reference signal sig0.

Return type:

1d array

apply_inverse_map(transport_map, sig0, x)

Apply inverse transport map.

Parameters:

• transport_map (1d array) – Forward transport map. Inverse is computed in this function.

• sig0 (1d array) – Reference signal.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

1d array

forward(x0, sig0, x1, sig1, rm_edge=False)

Forward transform.

Parameters:

• x0 (1d array) – Independent axis variable of reference signal (sig0).

• sig0 (1d array) – Reference signal.

• x1 (1d array) – Independent axis variable of the signal to transform (sig1).

• sig1 (1d array) – Signal to transform.

Returns:

• sig1_cdt (1d array) – CDT of input signal sig1 (new definition).

• sig1_hat (1d array) – old definition.

• xilde (1d array) – Independent axis variable in CDT space.

inverse(transport_map, sig0, x1)

Inverse transform.

Parameters:

• transport_map (1d array) – Forward transport map.

• sig0 (1d array) – Reference signal.

• x1 (1d array) – Independent axis variable of the signal to reconstruct.

Returns:

sig1_recon – Reconstructed signal.

Return type:

1d array

class pytranskit.optrans.continuous.radoncdt3D.RadonCDT3D(Npoints=1024, use_gpu=False)

Bases: BaseTransform

3D Radon Cumulative Distribution Transform.

Parameters:

• Npoints (scaler, number of radon projections) –

• use_gpu (boolean, use GPU if True) –

apply_inverse_map(transport_map, sig0, x_range)

Appy inverse transport map.

Parameters:

• transport_map (2d array, shape (t, N)) – Forward transport map. Inverse is computed in this function.

• sig0 (array, shape (height, width, depth)) – Reference signal.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

dRamp(x, h, d=3)

forward(x0_range, sig0, x1_range, sig1, rm_edge=False)

Forward transform.

Parameters:

• sig0 (array, shape (height, width, depth)) – Reference image.

• sig1 (array, shape (height, width, depth)) – Signal to transform.

Returns:

rcdt – 3D Radon-CDT of input image sig1.

Return type:

array, shape (t, N)

get_rotation_matrix(phi, theta, forward=True)

Calculates the Rotation matrix for, Input:

phi: Rotation along the x-axis (1,0,0) in degrees theta: Rotation along the y-axis (0,1,0) in degrees

output:

M: The rotation matrix

inverse(transport_map, sig0, x1)

Inverse transform.

Returns:

sig1_recon – Reconstructed signal sig1.

Return type:

array, shape (height, width)

iradon_3D(img3Dhat, phis, thetas)

radon_3D(img3D, N)

rotate3D(img3D, phi, theta, forward=True)

Rotates an image around the x and y axis with, Inputs:

img3D: Input image numpy 3D phi: Angle rotation around x theta: Angle rotation around y

Returns:

Rotated image (numpy 3D)

Return type:

img3D_rot

sample_sphere(num_pts, return_type='spherical')

This function “uniformly” samples a sphere on num_pts Inputs:

num_pts= number of points to sample return_type = return points in ‘spherical’ or ‘cartesian’ coordinates

pytranskit.optrans.continuous.radoncdt3D.assert_equal_shape(a, b, names=None)

Throw a ValueError if a and b are not the same shape.

pytranskit.optrans.continuous.radoncdt3D.check_array(array, ndim=None, dtype='numeric', force_all_finite=True, force_strictly_positive=False)

Input validation on an array, list, or similar.

Parameters:

• array (object) – Input object to check/convert

• ndim (int or None (default=None)) – Number of dimensions that array should have. If None, the dimensions are not checked

• dtype (string, type, list of types or None (default='numeric')) – Data type of result. If None, the dtype of the input is preserved. If ‘numeric’, dtype is preserved unless array.dtype is object. If dtype is a list of types, conversion on the first type is only performed if the dtype of the input is not in the list.

• force_all_finite (boolean (default=True)) – Whether to raise an error on np.inf and np.nan in array

• force_strictly_positive (boolean (default=False)) – Whether to raise an error if any array elements are <= 0

Returns:

array_converted – The converted and validated array.

Return type:

object

pytranskit.optrans.continuous.radoncdt3D.fft(x, n=None, axis=-1, overwrite_x=False)

Return discrete Fourier transform of real or complex sequence.

The returned complex array contains y(0), y(1),..., y(n-1), where

y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum().

Parameters:

• x (array_like) – Array to Fourier transform.

• n (int, optional) – Length of the Fourier transform. If n < x.shape[axis], x is truncated. If n > x.shape[axis], x is zero-padded. The default results in n = x.shape[axis].

• axis (int, optional) – Axis along which the fft’s are computed; the default is over the last axis (i.e., axis=-1).

• overwrite_x (bool, optional) – If True, the contents of x can be destroyed; the default is False.

Returns:

z –

with the elements:

[y(0),y(1),..,y(n/2),y(1-n/2),...,y(-1)]        if n is even
[y(0),y(1),..,y((n-1)/2),y(-(n-1)/2),...,y(-1)]  if n is odd

where:

y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n), j = 0..n-1

Return type:

complex ndarray

See also

ifft

Inverse FFT

rfft

FFT of a real sequence

Notes

The packing of the result is “standard”: If A = fft(a, n), then A[0] contains the zero-frequency term, A[1:n/2] contains the positive-frequency terms, and A[n/2:] contains the negative-frequency terms, in order of decreasingly negative frequency. So ,for an 8-point transform, the frequencies of the result are [0, 1, 2, 3, -4, -3, -2, -1]. To rearrange the fft output so that the zero-frequency component is centered, like [-4, -3, -2, -1, 0, 1, 2, 3], use fftshift.

Both single and double precision routines are implemented. Half precision inputs will be converted to single precision. Non-floating-point inputs will be converted to double precision. Long-double precision inputs are not supported.

This function is most efficient when n is a power of two, and least efficient when n is prime.

Note that if x is real-valued, then A[j] == A[n-j].conjugate(). If x is real-valued and n is even, then A[n/2] is real.

If the data type of x is real, a “real FFT” algorithm is automatically used, which roughly halves the computation time. To increase efficiency a little further, use rfft, which does the same calculation, but only outputs half of the symmetrical spectrum. If the data is both real and symmetrical, the dct can again double the efficiency by generating half of the spectrum from half of the signal.

Examples

>>> from scipy.fftpack import fft, ifft
>>> x = np.arange(5)
>>> np.allclose(fft(ifft(x)), x, atol=1e-15)  # within numerical accuracy.
True

pytranskit.optrans.continuous.radoncdt3D.fftshift(x, axes=None)

Shift the zero-frequency component to the center of the spectrum.

This function swaps half-spaces for all axes listed (defaults to all). Note that y[0] is the Nyquist component only if len(x) is even.

Parameters:

• x (array_like) – Input array.

• axes (int or shape tuple, optional) – Axes over which to shift. Default is None, which shifts all axes.

Returns:

y – The shifted array.

Return type:

ndarray

See also

ifftshift

The inverse of fftshift.

Examples

>>> freqs = np.fft.fftfreq(10, 0.1)
>>> freqs
array([ 0.,  1.,  2., ..., -3., -2., -1.])
>>> np.fft.fftshift(freqs)
array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.])

Shift the zero-frequency component only along the second axis:

>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0.,  1.,  2.],
       [ 3.,  4., -4.],
       [-3., -2., -1.]])
>>> np.fft.fftshift(freqs, axes=(1,))
array([[ 2.,  0.,  1.],
       [-4.,  3.,  4.],
       [-1., -3., -2.]])

pytranskit.optrans.continuous.radoncdt3D.ifft(x, n=None, axis=-1, overwrite_x=False)

Return discrete inverse Fourier transform of real or complex sequence.

The returned complex array contains y(0), y(1),..., y(n-1), where

y(j) = (x * exp(2*pi*sqrt(-1)*j*np.arange(n)/n)).mean().

Parameters:

• x (array_like) – Transformed data to invert.

• n (int, optional) – Length of the inverse Fourier transform. If n < x.shape[axis], x is truncated. If n > x.shape[axis], x is zero-padded. The default results in n = x.shape[axis].

• axis (int, optional) – Axis along which the ifft’s are computed; the default is over the last axis (i.e., axis=-1).

• overwrite_x (bool, optional) – If True, the contents of x can be destroyed; the default is False.

Returns:

ifft – The inverse discrete Fourier transform.

Return type:

ndarray of floats

See also

fft

Forward FFT

Notes

Both single and double precision routines are implemented. Half precision inputs will be converted to single precision. Non-floating-point inputs will be converted to double precision. Long-double precision inputs are not supported.

This function is most efficient when n is a power of two, and least efficient when n is prime.

If the data type of x is real, a “real IFFT” algorithm is automatically used, which roughly halves the computation time.

Examples

>>> from scipy.fftpack import fft, ifft
>>> import numpy as np
>>> x = np.arange(5)
>>> np.allclose(ifft(fft(x)), x, atol=1e-15)  # within numerical accuracy.
True

pytranskit.optrans.continuous.radoncdt3D.signal_to_pdf(input, sigma=0.0, epsilon=1e-08, total=1.0)

Get the (smoothed) probability density function of a signal.

Performs the following operations: 1. Smooth sigma with a Gaussian filter 2. Normalize signal such that it sums to 1 3. Add epsilon to ensure signal is strictly positive 4. Re-normalize signal such that it sums to total

Parameters:

• input (ndarray) – Input array

• sigma (scalar) – Standard deviation for Gaussian kernel. The standard deviations of the Gaussian filter are given for each axis as a sequence, or as a single number, in which case it is equal for all axes.

• epsilon (scalar) – Offset to ensure that signal is strictly positive.

• total (scalar) – Value of the signal summation.

Returns:

pdf – Returned array of same shape as input

Return type:

ndarray

scdt

SCDT.py

Class for computing the Signed Cumulative Distribution Transform (SCDT).

Authors: Sumati Thareja D.M. Rocio Ivan Medri Akram Aldroubi Gustavo Rohde

Based on paper: ArXiv paper link

class pytranskit.optrans.continuous.scdt.SCDT(reference, x0=None)

Bases: object

Signed Cumulative Distribution Transform (SCDT)

Parameters:

• reference (1D array representing the reference density) –

• x0 (domain of reference) –

• reference_CDF (reference CDF) –

• xtilde (Domain of the reference's inverse CDF) –

• reference_CDT_inverse (inverse CDF of reference) –

calc_scdt(sig1, t1, s0, t0)

gen_inverse(f, dom_f, dom_gf)

gen_inverse calculates the generalized inverse of the function f input:

f: A one dimensional function represented by a vector dom_f: The domain of the function f dom_gf: The domain of the generalizd inverse of f

output:

The generalized inverse of f

istransform(Iposhat, Ineghat, Masspos, Massneg)

istrasform calculates the inverse of the stransform. input:

The 4 components of the transport transform for signed signals

output:

The original signal

itransform(Ihat)

itransform calculates the inverse of the CDT. It receives a transport displacement and the reference, and finds the one dimensional distribution I from it. input:

Transport displacement map The reference used for calculating the CDT

output:

I: The original distribution

stransform(I, x=None)

stransform calculates the transport transform (CDT) of a signal I for signals that may change sign input:

The original density I x -> domain of the density I

output:

The 4 components of the transform of signed signals: the CDT of the positive and the negative part of I, and the total masses of the positive and the negative part of I

transform(I)

transform calculates the transport map that morphs the one-dimensional distribution I into the reference. input:

I: A one dimensional distributions of size self.dim

output:

The CDT transformation of I

pytranskit.optrans.continuous.scdt.interp(x, xp, fp, left=None, right=None, period=None)

One-dimensional linear interpolation for monotonically increasing sample points.

Returns the one-dimensional piecewise linear interpolant to a function with given discrete data points (xp, fp), evaluated at x.

Parameters:

• x (array_like) – The x-coordinates at which to evaluate the interpolated values.

• xp (1-D sequence of floats) – The x-coordinates of the data points, must be increasing if argument period is not specified. Otherwise, xp is internally sorted after normalizing the periodic boundaries with xp = xp % period.

• fp (1-D sequence of float or complex) – The y-coordinates of the data points, same length as xp.

• left (optional float or complex corresponding to fp) – Value to return for x < xp[0], default is fp[0].

• right (optional float or complex corresponding to fp) – Value to return for x > xp[-1], default is fp[-1].

• period (None or float, optional) –

  A period for the x-coordinates. This parameter allows the proper interpolation of angular x-coordinates. Parameters left and right are ignored if period is specified.

  New in version 1.10.0.

Returns:

y – The interpolated values, same shape as x.

Return type:

float or complex (corresponding to fp) or ndarray

Raises:

ValueError – If xp and fp have different length If xp or fp are not 1-D sequences If period == 0

See also

scipy.interpolate

Warning

The x-coordinate sequence is expected to be increasing, but this is not explicitly enforced. However, if the sequence xp is non-increasing, interpolation results are meaningless.

Note that, since NaN is unsortable, xp also cannot contain NaNs.

A simple check for xp being strictly increasing is:

np.all(np.diff(xp) > 0)

Examples

>>> xp = [1, 2, 3]
>>> fp = [3, 2, 0]
>>> np.interp(2.5, xp, fp)
1.0
>>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
array([3.  , 3.  , 2.5 , 0.56, 0.  ])
>>> UNDEF = -99.0
>>> np.interp(3.14, xp, fp, right=UNDEF)
-99.0

Plot an interpolant to the sine function:

>>> x = np.linspace(0, 2*np.pi, 10)
>>> y = np.sin(x)
>>> xvals = np.linspace(0, 2*np.pi, 50)
>>> yinterp = np.interp(xvals, x, y)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'o')
[<matplotlib.lines.Line2D object at 0x...>]
>>> plt.plot(xvals, yinterp, '-x')
[<matplotlib.lines.Line2D object at 0x...>]
>>> plt.show()

Interpolation with periodic x-coordinates:

>>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
>>> xp = [190, -190, 350, -350]
>>> fp = [5, 10, 3, 4]
>>> np.interp(x, xp, fp, period=360)
array([7.5 , 5.  , 8.75, 6.25, 3.  , 3.25, 3.5 , 3.75])

Complex interpolation:

>>> x = [1.5, 4.0]
>>> xp = [2,3,5]
>>> fp = [1.0j, 0, 2+3j]
>>> np.interp(x, xp, fp)
array([0.+1.j , 1.+1.5j])