# Source code for colour.algebra.interpolation

# -*- coding: utf-8 -*-
"""
Interpolation
=============

Defines classes and definitions for interpolating variables.

-   :class:colour.KernelInterpolator: 1-D function generic interpolation with
arbitrary kernel.
-   :class:colour.NearestNeighbourInterpolator: 1-D function
nearest-neighbour interpolation.
-   :class:colour.LinearInterpolator: 1-D function linear interpolation.
-   :class:colour.SpragueInterpolator: 1-D function fifth-order polynomial
interpolation using *Sprague (1880)* method.
-   :class:colour.CubicSplineInterpolator: 1-D function cubic spline
interpolation.
-   :class:colour.PchipInterpolator: 1-D function piecewise cube Hermite
interpolation.
-   :class:colour.NullInterpolator: 1-D function null interpolation.
-   :func:colour.lagrange_coefficients: Computation of
*Lagrange Coefficients*.
-   :func:colour.algebra.table_interpolation_trilinear: Trilinear
interpolation with table.
-   :func:colour.algebra.table_interpolation_tetrahedral: Tetrahedral
interpolation with table.
-   :attr:colour.TABLE_INTERPOLATION_METHODS: Supported table interpolation
methods.
-   :func:colour.table_interpolation: Interpolation with table using given
method.

References
----------
-   :cite:Bourkeb : Bourke, P. (n.d.). Trilinear Interpolation. Retrieved
from http://paulbourke.net/miscellaneous/interpolation/
-   :cite:Burger2009b : Burger, W., & Burge, M. J. (2009). Principles of
Digital Image Processing. London: Springer London.
doi:10.1007/978-1-84800-195-4
-   :cite:CIETC1-382005f : CIE TC 1-38. (2005). 9.2.4 Method of
interpolation for uniformly spaced independent variable. In CIE 167:2005
Recommended Practice for Tabulating Spectral Data for Use in Colour
Computations (pp. 1-27). ISBN:978-3-901-90641-1
-   :cite:CIETC1-382005h : CIE TC 1-38. (2005). Table V. Values of the
c-coefficients of Equ.s 6 and 7. In CIE 167:2005 Recommended Practice for
Tabulating Spectral Data for Use in Colour Computations (p. 19).
ISBN:978-3-901-90641-1
-   :cite:Fairman1985b : Fairman, H. S. (1985). The calculation of weight
factors for tristimulus integration. Color Research & Application, 10(4),
199-203. doi:10.1002/col.5080100407
-   :cite:Kirk2006 : Kirk, R. (2006). Truelight Software Library 2.0.
FL-TL-TN-0057-SoftwareLib.pdf
-   :cite:Westland2012h : Westland, S., Ripamonti, C., & Cheung, V. (2012).
Interpolation Methods. In Computational Colour Science Using MATLAB
(2nd ed., pp. 29-37). ISBN:978-0-470-66569-5
-   :cite:Wikipedia2005b : Wikipedia. (2005). Lanczos resampling. Retrieved
October 14, 2017, from https://en.wikipedia.org/wiki/Lanczos_resampling
-   :cite:Wikipedia2003a : Wikipedia. (2003). Lagrange polynomial -
Definition. Retrieved January 20, 2016, from https://en.wikipedia.org/\
wiki/Lagrange_polynomial#Definition
"""

from __future__ import division, unicode_literals

import itertools
import numpy as np
import scipy.interpolate
from collections import OrderedDict, Mapping
from six.moves import reduce

from colour.constants import DEFAULT_FLOAT_DTYPE, DEFAULT_INT_DTYPE
from colour.utilities import (CaseInsensitiveMapping, as_float_array, as_float,
closest_indexes, interval, is_integer,
is_numeric, runtime_warning, tsplit)

__author__ = 'Colour Developers'
__maintainer__ = 'Colour Developers'
__status__ = 'Production'

__all__ = [
'kernel_nearest_neighbour', 'kernel_linear', 'kernel_sinc',
'kernel_lanczos', 'kernel_cardinal_spline', 'KernelInterpolator',
'NearestNeighbourInterpolator', 'LinearInterpolator',
'SpragueInterpolator', 'CubicSplineInterpolator', 'PchipInterpolator',
'NullInterpolator', 'lagrange_coefficients',
'vertices_and_relative_coordinates', 'table_interpolation_trilinear',
'table_interpolation_tetrahedral', 'TABLE_INTERPOLATION_METHODS',
'table_interpolation'
]

[docs]def kernel_nearest_neighbour(x):
"""
Returns the *nearest-neighbour* kernel evaluated at given samples.

Parameters
----------
x : array_like
Samples at which to evaluate the *nearest-neighbour* kernel.

Returns
-------
ndarray
The *nearest-neighbour* kernel evaluated at given samples.

References
----------
:cite:Burger2009b

Examples
--------
>>> kernel_nearest_neighbour(np.linspace(0, 1, 10))
array([1, 1, 1, 1, 1, 0, 0, 0, 0, 0])
"""

return np.where(np.abs(x) < 0.5, 1, 0)

[docs]def kernel_linear(x):
"""
Returns the *linear* kernel evaluated at given samples.

Parameters
----------
x : array_like
Samples at which to evaluate the *linear* kernel.

Returns
-------
ndarray
The *linear* kernel evaluated at given samples.

References
----------
:cite:Burger2009b

Examples
--------
>>> kernel_linear(np.linspace(0, 1, 10))  # doctest: +ELLIPSIS
array([ 1.        ,  0.8888888...,  0.7777777...,  \
0.6666666...,  0.5555555...,
0.4444444...,  0.3333333...,  0.2222222...,  \
0.1111111...,  0.        ])
"""

return np.where(np.abs(x) < 1, 1 - np.abs(x), 0)

[docs]def kernel_sinc(x, a=3):
"""
Returns the *sinc* kernel evaluated at given samples.

Parameters
----------
x : array_like
Samples at which to evaluate the *sinc* kernel.
a : int, optional
Size of the *sinc* kernel.

Returns
-------
ndarray
The *sinc* kernel evaluated at given samples.

References
----------
:cite:Burger2009b

Examples
--------
>>> kernel_sinc(np.linspace(0, 1, 10))  # doctest: +ELLIPSIS
array([  1.0000000...e+00,   9.7981553...e-01,   9.2072542...e-01,
8.2699334...e-01,   7.0531659...e-01,   5.6425327...e-01,
4.1349667...e-01,   2.6306440...e-01,   1.2247694...e-01,
3.8981718...e-17])
"""

assert a >= 1, '"a" must be equal or superior to 1!'

return np.where(np.abs(x) < a, np.sinc(x), 0)

[docs]def kernel_lanczos(x, a=3):
"""
Returns the *lanczos* kernel evaluated at given samples.

Parameters
----------
x : array_like
Samples at which to evaluate the *lanczos* kernel.
a : int, optional
Size of the *lanczos* kernel.

Returns
-------
ndarray
The *lanczos* kernel evaluated at given samples.

References
----------
:cite:Wikipedia2005b

Examples
--------
>>> kernel_lanczos(np.linspace(0, 1, 10))  # doctest: +ELLIPSIS
array([  1.0000000...e+00,   9.7760615...e-01,   9.1243770...e-01,
8.1030092...e-01,   6.8012706...e-01,   5.3295773...e-01,
3.8071690...e-01,   2.3492839...e-01,   1.0554054...e-01,
3.2237621...e-17])
"""

assert a >= 1, '"a" must be equal or superior to 1!'

return np.where(np.abs(x) < a, np.sinc(x) * np.sinc(x / a), 0)

[docs]def kernel_cardinal_spline(x, a=0.5, b=0.0):
"""
Returns the *cardinal spline* kernel evaluated at given samples.

Notable *cardinal spline* :math:a and :math:b parameterizations:

-   *Catmull-Rom*: :math:(a=0.5, b=0)
-   *Cubic B-Spline*: :math:(a=0, b=1)
-   *Mitchell-Netravalli*: :math:(a=\\cfrac{1}{3}, b=\\cfrac{1}{3})

Parameters
----------
x : array_like
Samples at which to evaluate the *cardinal spline* kernel.
a : int, optional
:math:a control parameter.
b : int, optional
:math:b control parameter.

Returns
-------
ndarray
The *cardinal spline* kernel evaluated at given samples.

References
----------
:cite:Burger2009b

Examples
--------
>>> kernel_cardinal_spline(np.linspace(0, 1, 10))  # doctest: +ELLIPSIS
array([ 1.        ,  0.9711934...,  0.8930041...,  \
0.7777777...,  0.6378600...,
0.4855967...,  0.3333333...,  0.1934156...,  \
0.0781893...,  0.        ])
"""

x_abs = np.abs(x)
y = np.where(
x_abs < 1,
(-6 * a - 9 * b + 12) * x_abs ** 3 + (6 * a + 12 * b - 18) * x_abs ** 2
- 2 * b + 6,
(-6 * a - b) * x_abs ** 3 + (30 * a + 6 * b) * x_abs ** 2 +
(-48 * a - 12 * b) * x_abs + 24 * a + 8 * b,
)
y[x_abs >= 2] = 0

return 1 / 6 * y

[docs]class KernelInterpolator(object):
"""
Kernel based interpolation of a 1-D function.

The reconstruction of a continuous signal can be described as a linear
convolution operation. Interpolation can be expressed as a convolution of
the given discrete function :math:g(x) with some continuous interpolation
kernel :math:k(w):

:math:\\hat{g}(w_0) = [k * g](w_0) = \
\\sum_{x=-\\infty}^{\\infty}k(w_0 - x)\\cdot g(x)

Parameters
----------
x : array_like
Independent :math:x variable values corresponding with :math:y
variable.
y : array_like
Dependent and already known :math:y variable values to
interpolate.
window : int, optional
Width of the window in samples on each side.
kernel : callable, optional
Kernel to use for interpolation.
kernel_args : dict, optional
Arguments to use when calling the kernel.
Arguments to use when padding :math:y variable values with the
:func:np.pad definition.
dtype : type
Data type used for internal conversions.

Attributes
----------
x
y
window
kernel
kernel_args

Methods
-------
__call__

References
----------
:cite:Burger2009b, :cite:Wikipedia2005b

Examples
--------
Interpolating a single numeric variable:

>>> y = np.array([5.9200, 9.3700, 10.8135, 4.5100,
...               69.5900, 27.8007, 86.0500])
>>> x = np.arange(len(y))
>>> f = KernelInterpolator(x, y)
>>> f(0.5)  # doctest: +ELLIPSIS
6.9411400...

Interpolating an *array_like* variable:

>>> f([0.25, 0.75])  # doctest: +ELLIPSIS
array([ 6.1806208...,  8.0823848...])

Using a different *lanczos* kernel:

>>> f = KernelInterpolator(x, y, kernel=kernel_sinc)
>>> f([0.25, 0.75])  # doctest: +ELLIPSIS
array([ 6.5147317...,  8.3965466...])

Using a different window size:

>>> f = KernelInterpolator(
...     x,
...     y,
...     window=16,
...     kernel=kernel_lanczos,
...     kernel_args={'a': 16})
>>> f([0.25, 0.75])  # doctest: +ELLIPSIS
array([ 5.3961792...,  5.6521093...])
"""

[docs]    def __init__(self,
x,
y,
window=3,
kernel=kernel_lanczos,
kernel_args=None,
dtype=DEFAULT_FLOAT_DTYPE):
self._x_p = None
self._y_p = None

self._x = None
self._y = None
self._window = None
self._dtype = dtype

self.x = x
self.y = y
self.window = window

self._kernel = None
self.kernel = kernel
self._kernel_args = {}
self.kernel_args = kernel_args

self._validate_dimensions()

@property
def x(self):
"""
Getter and setter property for the independent :math:x variable.

Parameters
----------
value : array_like
Value to set the independent :math:x variable with.

Returns
-------
array_like
Independent :math:x variable.
"""

return self._x

@x.setter
def x(self, value):
"""
Setter for the **self.x** property.
"""

if value is not None:
value = np.atleast_1d(value).astype(self._dtype)

assert value.ndim == 1, (
'"x" independent variable must have exactly one dimension!')

value_interval = interval(value)

if value_interval.size != 1:
runtime_warning(('"x" independent variable is not uniform, '
'unpredictable results may occur!'))

self._x = value

if self._window is not None:
self._x, (self._window, self._window),
'linear_ramp',
end_values=(
np.min(self._x) - self._window * value_interval[0],
np.max(self._x) + self._window * value_interval[0]))

@property
def y(self):
"""
Getter and setter property for the dependent and already known
:math:y variable.

Parameters
----------
value : array_like
Value to set the dependent and already known :math:y variable
with.

Returns
-------
array_like
Dependent and already known :math:y variable.
"""

return self._y

@y.setter
def y(self, value):
"""
Setter for the **self.y** property.
"""

if value is not None:
value = np.atleast_1d(value).astype(self._dtype)

assert value.ndim == 1, (
'"y" dependent variable must have exactly one dimension!')

self._y = value

if self._window is not None:

@property
def window(self):
"""
Getter and setter property for the window.

Parameters
----------
value : int
Value to set the window with.

Returns
-------
int
Window.
"""

return self._window

@window.setter
def window(self, value):
"""
Setter for the **self.window** property.
"""

if value is not None:
assert is_integer(value), '"window" must be an integer!'

assert value >= 1, '"window" must be equal or superior to 1!'

self._window = value

# Triggering "self._x_p" update.
if self._x is not None:
self.x = self._x

# Triggering "self._y_p" update.
if self._y is not None:
self.y = self._y

@property
def kernel(self):
"""
Getter and setter property for the kernel callable.

Parameters
----------
value : callable
Value to set the kernel callable.

Returns
-------
callable
Kernel callable.
"""

return self._kernel

@kernel.setter
def kernel(self, value):
"""
Setter for the **self.kernel** property.
"""

if value is not None:
assert hasattr(
value,
'__call__'), ('"{0}" attribute: "{1}" is not callable!'.format(
'kernel', value))

self._kernel = value

@property
def kernel_args(self):
"""
Getter and setter property for the kernel call time arguments.

Parameters
----------
value : dict
Value to call the interpolation kernel with.

Returns
-------
dict
Kernel call time arguments.
"""

return self._kernel_args

@kernel_args.setter
def kernel_args(self, value):
"""
Setter for the **self.kernel_args** property.
"""

if value is not None:
assert isinstance(value, (dict, OrderedDict)), (
'"{0}" attribute: "{1}" type is not "dict" or "OrderedDict"!'
).format('kernel_args', value)

self._kernel_args = value

@property
"""
Getter and setter property for the kernel call time arguments.

Parameters
----------
value : dict
Value to call the interpolation kernel with.

Returns
-------
dict
Kernel call time arguments.
"""

"""
"""

if value is not None:
assert isinstance(value, Mapping), (
'"{0}" attribute: "{1}" type is not a "Mapping" instance!'

# Triggering "self._y_p" update.
if self._y is not None:
self.y = self._y

[docs]    def __call__(self, x):
"""
Evaluates the interpolator at given point(s).

Parameters
----------
x : numeric or array_like
Point(s) to evaluate the interpolant at.

Returns
-------
float or ndarray
Interpolated value(s).
"""

x = np.atleast_1d(x).astype(self._dtype)

xi = as_float(self._evaluate(x))

return xi

def _evaluate(self, x):
"""
Performs the interpolator evaluation at given points.

Parameters
----------
x : ndarray
Points to evaluate the interpolant at.

Returns
-------
ndarray
Interpolated points values.
"""

self._validate_dimensions()
self._validate_interpolation_range(x)

x_interval = interval(self._x)[0]
x_f = np.floor(x / x_interval)

windows = (x_f[:, np.newaxis] + np.arange(-self._window + 1,
self._window + 1))
clip_l = min(self._x_p) / x_interval
clip_h = max(self._x_p) / x_interval
windows = np.clip(windows, clip_l, clip_h) - clip_l
windows = np.around(windows).astype(DEFAULT_INT_DTYPE)

return np.sum(
self._y_p[windows] * self._kernel(
x[:, np.newaxis] / x_interval - windows -
min(self._x_p) / x_interval, **self._kernel_args),
axis=-1)

def _validate_dimensions(self):
"""
Validates variables dimensions to be the same.
"""

if len(self._x) != len(self._y):
raise ValueError(
('"x" independent and "y" dependent variables have different '
'dimensions: "{0}", "{1}"').format(
len(self._x), len(self._y)))

def _validate_interpolation_range(self, x):
"""
Validates given point to be in interpolation range.
"""

below_interpolation_range = x < self._x[0]
above_interpolation_range = x > self._x[-1]

if below_interpolation_range.any():
raise ValueError('"{0}" is below interpolation range.'.format(x))

if above_interpolation_range.any():
raise ValueError('"{0}" is above interpolation range.'.format(x))

[docs]class NearestNeighbourInterpolator(KernelInterpolator):
"""
A nearest-neighbour interpolator.

Other Parameters
----------------
x : array_like
Independent :math:x variable values corresponding with :math:y
variable.
y : array_like
Dependent and already known :math:y variable values to
interpolate.
window : int, optional
Width of the window in samples on each side.
Arguments to use when padding :math:y variable values with the
:func:np.pad definition.
dtype : type
Data type used for internal conversions.
"""

[docs]    def __init__(self, *args, **kwargs):
kwargs['kernel'] = kernel_nearest_neighbour
if 'kernel_args' in kwargs:
del kwargs['kernel_args']

super(NearestNeighbourInterpolator, self).__init__(*args, **kwargs)

[docs]class LinearInterpolator(object):
"""
Linearly interpolates a 1-D function.

Parameters
----------
x : array_like
Independent :math:x variable values corresponding with :math:y
variable.
y : array_like
Dependent and already known :math:y variable values to
interpolate.
dtype : type
Data type used for internal conversions.

Attributes
----------
x
y

Methods
-------
__call__

Notes
-----
-   This class is a wrapper around *numpy.interp* definition.

Examples
--------
Interpolating a single numeric variable:

>>> y = np.array([5.9200, 9.3700, 10.8135, 4.5100,
...               69.5900, 27.8007, 86.0500])
>>> x = np.arange(len(y))
>>> f = LinearInterpolator(x, y)
>>> # Doctests ellipsis for Python 2.x compatibility.
>>> f(0.5)  # doctest: +ELLIPSIS
7.64...

Interpolating an *array_like* variable:

>>> f([0.25, 0.75])
array([ 6.7825,  8.5075])
"""

[docs]    def __init__(self, x, y, dtype=DEFAULT_FLOAT_DTYPE):
self._x = None
self._y = None
self._dtype = dtype

self.x = x
self.y = y

self._validate_dimensions()

@property
def x(self):
"""
Getter and setter property for the independent :math:x variable.

Parameters
----------
value : array_like
Value to set the independent :math:x variable with.

Returns
-------
array_like
Independent :math:x variable.
"""

return self._x

@x.setter
def x(self, value):
"""
Setter for the **self.x** property.
"""

if value is not None:
value = np.atleast_1d(value).astype(self._dtype)

assert value.ndim == 1, (
'"x" independent variable must have exactly one dimension!')

self._x = value

@property
def y(self):
"""
Getter and setter property for the dependent and already known
:math:y variable.

Parameters
----------
value : array_like
Value to set the dependent and already known :math:y variable
with.

Returns
-------
array_like
Dependent and already known :math:y variable.
"""

return self._y

@y.setter
def y(self, value):
"""
Setter for the **self.y** property.
"""

if value is not None:
value = np.atleast_1d(value).astype(self._dtype)

assert value.ndim == 1, (
'"y" dependent variable must have exactly one dimension!')

self._y = value

[docs]    def __call__(self, x):
"""
Evaluates the interpolating polynomial at given point(s).

Parameters
----------
x : numeric or array_like
Point(s) to evaluate the interpolant at.

Returns
-------
float or ndarray
Interpolated value(s).
"""

x = np.atleast_1d(x).astype(self._dtype)

xi = as_float(self._evaluate(x))

return xi

def _evaluate(self, x):
"""
Performs the interpolating polynomial evaluation at given points.

Parameters
----------
x : ndarray
Points to evaluate the interpolant at.

Returns
-------
ndarray
Interpolated points values.
"""

self._validate_dimensions()
self._validate_interpolation_range(x)

return np.interp(x, self._x, self._y)

def _validate_dimensions(self):
"""
Validates variables dimensions to be the same.
"""

if len(self._x) != len(self._y):
raise ValueError(
('"x" independent and "y" dependent variables have different '
'dimensions: "{0}", "{1}"').format(
len(self._x), len(self._y)))

def _validate_interpolation_range(self, x):
"""
Validates given point to be in interpolation range.
"""

below_interpolation_range = x < self._x[0]
above_interpolation_range = x > self._x[-1]

if below_interpolation_range.any():
raise ValueError('"{0}" is below interpolation range.'.format(x))

if above_interpolation_range.any():
raise ValueError('"{0}" is above interpolation range.'.format(x))

[docs]class SpragueInterpolator(object):
"""
Constructs a fifth-order polynomial that passes through :math:y dependent
variable.

*Sprague (1880)* method is recommended by the *CIE* for interpolating
functions having a uniformly spaced independent variable.

Parameters
----------
x : array_like
Independent :math:x variable values corresponding with :math:y
variable.
y : array_like
Dependent and already known :math:y variable values to
interpolate.
dtype : type
Data type used for internal conversions.

Attributes
----------
x
y

Methods
-------
__call__

Notes
-----
-   The minimum number :math:k of data points required along the
interpolation axis is :math:k=6.

References
----------
:cite:CIETC1-382005f, :cite:Westland2012h

Examples
--------
Interpolating a single numeric variable:

>>> y = np.array([5.9200, 9.3700, 10.8135, 4.5100,
...               69.5900, 27.8007, 86.0500])
>>> x = np.arange(len(y))
>>> f = SpragueInterpolator(x, y)
>>> f(0.5)  # doctest: +ELLIPSIS
7.2185025...

Interpolating an *array_like* variable:

>>> f([0.25, 0.75])  # doctest: +ELLIPSIS
array([ 6.7295161...,  7.8140625...])
"""

SPRAGUE_C_COEFFICIENTS = np.array([
[884, -1960, 3033, -2648, 1080, -180],
[508, -540, 488, -367, 144, -24],
[-24, 144, -367, 488, -540, 508],
[-180, 1080, -2648, 3033, -1960, 884],
])
"""
Defines the coefficients used to generate extra points for boundaries
interpolation.

SPRAGUE_C_COEFFICIENTS : array_like, (4, 6)

References
----------
:cite:CIETC1-382005h
"""

[docs]    def __init__(self, x, y, dtype=DEFAULT_FLOAT_DTYPE):
self._xp = None
self._yp = None

self._x = None
self._y = None
self._dtype = dtype

self.x = x
self.y = y

self._validate_dimensions()

@property
def x(self):
"""
Getter and setter property for the independent :math:x variable.

Parameters
----------
value : array_like
Value to set the independent :math:x variable with.

Returns
-------
array_like
Independent :math:x variable.
"""

return self._x

@x.setter
def x(self, value):
"""
Setter for the **self.x** property.
"""

if value is not None:
value = np.atleast_1d(value).astype(self._dtype)

assert value.ndim == 1, (
'"x" independent variable must have exactly one dimension!')

value_interval = interval(value)[0]

xp1 = value[0] - value_interval * 2
xp2 = value[0] - value_interval
xp3 = value[-1] + value_interval
xp4 = value[-1] + value_interval * 2

self._xp = np.concatenate(((xp1, xp2), value, (xp3, xp4)))

self._x = value

@property
def y(self):
"""
Getter and setter property for the dependent and already known
:math:y variable.

Parameters
----------
value : array_like
Value to set the dependent and already known :math:y variable
with.

Returns
-------
array_like
Dependent and already known :math:y variable.
"""

return self._y

@y.setter
def y(self, value):
"""
Setter for the **self.y** property.
"""

if value is not None:
value = np.atleast_1d(value).astype(self._dtype)

assert value.ndim == 1, (
'"y" dependent variable must have exactly one dimension!')

assert len(value) >= 6, (
'"y" dependent variable values count must be normalised to'
'domain [6:]!')

yp1 = np.ravel((np.dot(self.SPRAGUE_C_COEFFICIENTS[0],
np.array(value[0:6]).reshape(
(6, 1)))) / 209)[0]
yp2 = np.ravel((np.dot(self.SPRAGUE_C_COEFFICIENTS[1],
np.array(value[0:6]).reshape(
(6, 1)))) / 209)[0]
yp3 = np.ravel((np.dot(self.SPRAGUE_C_COEFFICIENTS[2],
np.array(value[-6:]).reshape(
(6, 1)))) / 209)[0]
yp4 = np.ravel((np.dot(self.SPRAGUE_C_COEFFICIENTS[3],
np.array(value[-6:]).reshape(
(6, 1)))) / 209)[0]

self._yp = np.concatenate(((yp1, yp2), value, (yp3, yp4)))

self._y = value

[docs]    def __call__(self, x):
"""
Evaluates the interpolating polynomial at given point(s).

Parameters
----------
x : numeric or array_like
Point(s) to evaluate the interpolant at.

Returns
-------
numeric or ndarray
Interpolated value(s).
"""

return self._evaluate(x)

def _evaluate(self, x):
"""
Performs the interpolating polynomial evaluation at given point.

Parameters
----------
x : numeric
Point to evaluate the interpolant at.

Returns
-------
float
Interpolated point values.
"""

x = as_float_array(x)

self._validate_dimensions()
self._validate_interpolation_range(x)

i = np.searchsorted(self._xp, x) - 1
X = (x - self._xp[i]) / (self._xp[i + 1] - self._xp[i])

r = self._yp

a0p = r[i]
a1p = ((2 * r[i - 2] - 16 * r[i - 1] + 16 * r[i + 1] -
2 * r[i + 2]) / 24)  # yapf: disable
a2p = ((-r[i - 2] + 16 * r[i - 1] - 30 * r[i] + 16 * r[i + 1] -
r[i + 2]) / 24)  # yapf: disable
a3p = ((-9 * r[i - 2] + 39 * r[i - 1] - 70 * r[i] + 66 * r[i + 1] -
33 * r[i + 2] + 7 * r[i + 3]) / 24)
a4p = ((13 * r[i - 2] - 64 * r[i - 1] + 126 * r[i] - 124 * r[i + 1] +
61 * r[i + 2] - 12 * r[i + 3]) / 24)
a5p = ((-5 * r[i - 2] + 25 * r[i - 1] - 50 * r[i] + 50 * r[i + 1] -
25 * r[i + 2] + 5 * r[i + 3]) / 24)

y = (a0p + a1p * X + a2p * X ** 2 + a3p * X ** 3 + a4p * X ** 4 +
a5p * X ** 5)

return y

def _validate_dimensions(self):
"""
Validates variables dimensions to be the same.
"""

if len(self._x) != len(self._y):
raise ValueError(
('"x" independent and "y" dependent variables have different '
'dimensions: "{0}", "{1}"').format(
len(self._x), len(self._y)))

def _validate_interpolation_range(self, x):
"""
Validates given point to be in interpolation range.
"""

below_interpolation_range = x < self._x[0]
above_interpolation_range = x > self._x[-1]

if below_interpolation_range.any():
raise ValueError('"{0}" is below interpolation range.'.format(x))

if above_interpolation_range.any():
raise ValueError('"{0}" is above interpolation range.'.format(x))

class CubicSplineInterpolator(scipy.interpolate.interp1d):
"""
Interpolates a 1-D function using cubic spline interpolation.

Notes
-----
-   This class is a wrapper around *scipy.interpolate.interp1d* class.
"""

def __init__(self, *args, **kwargs):
super(CubicSplineInterpolator, self).__init__(
kind='cubic', *args, **kwargs)

[docs]class PchipInterpolator(scipy.interpolate.PchipInterpolator):
"""
Interpolates a 1-D function using Piecewise Cubic Hermite Interpolating
Polynomial interpolation.

Attributes
----------
y

Notes
-----
-   This class is a wrapper around *scipy.interpolate.PchipInterpolator*
class.
"""

[docs]    def __init__(self, x, y, *args, **kwargs):
super(PchipInterpolator, self).__init__(x, y, *args, **kwargs)

self._y = y

@property
def y(self):
"""
Getter and setter property for the dependent and already known
:math:y variable.

Parameters
----------
value : array_like
Value to set the dependent and already known :math:y variable
with.

Returns
-------
array_like
Dependent and already known :math:y variable.
"""

return self._y

[docs]class NullInterpolator(object):
"""
Performs 1-D function null interpolation, i.e. a call within given
tolerances will return existing :math:y variable values and default
if outside tolerances.

Parameters
----------
x : ndarray
Independent :math:x variable values corresponding with :math:y
variable.
y : ndarray
Dependent and already known :math:y variable values to
interpolate.
absolute_tolerance : numeric, optional
Absolute tolerance.
relative_tolerance : numeric, optional
Relative tolerance.
default : numeric, optional
Default value for interpolation outside tolerances.
dtype : type
Data type used for internal conversions.

Attributes
----------
x
y
relative_tolerance
absolute_tolerance
default

Methods
-------
__call__

Examples
--------
>>> y = np.array([5.9200, 9.3700, 10.8135, 4.5100,
...               69.5900, 27.8007, 86.0500])
>>> x = np.arange(len(y))
>>> f = NullInterpolator(x, y)
>>> f(0.5)
nan
>>> f(1.0)  # doctest: +ELLIPSIS
9.3699999...
>>> f = NullInterpolator(x, y, absolute_tolerance=0.01)
>>> f(1.01)  # doctest: +ELLIPSIS
9.3699999...
"""

[docs]    def __init__(self,
x,
y,
absolute_tolerance=10e-7,
relative_tolerance=10e-7,
default=np.nan,
dtype=DEFAULT_FLOAT_DTYPE):
self._x = None
self._y = None
self._absolute_tolerance = None
self._relative_tolerance = None
self._default = None
self._dtype = dtype

self.x = x
self.y = y
self.absolute_tolerance = absolute_tolerance
self.relative_tolerance = relative_tolerance
self.default = default

self._validate_dimensions()

@property
def x(self):
"""
Getter and setter property for the independent :math:x variable.

Parameters
----------
value : array_like
Value to set the independent :math:x variable with.

Returns
-------
array_like
Independent :math:x variable.
"""

return self._x

@x.setter
def x(self, value):
"""
Setter for the **self.x** property.
"""

if value is not None:
value = np.atleast_1d(value).astype(self._dtype)

assert value.ndim == 1, (
'"x" independent variable must have exactly one dimension!')

self._x = value

@property
def y(self):
"""
Getter and setter property for the dependent and already known
:math:y variable.

Parameters
----------
value : array_like
Value to set the dependent and already known :math:y variable
with.

Returns
-------
array_like
Dependent and already known :math:y variable.
"""

return self._y

@y.setter
def y(self, value):
"""
Setter for the **self.y** property.
"""

if value is not None:
value = np.atleast_1d(value).astype(self._dtype)

assert value.ndim == 1, (
'"y" dependent variable must have exactly one dimension!')

self._y = value

@property
def relative_tolerance(self):
"""
Getter and setter property for the relative tolerance.

Parameters
----------
value : numeric
Value to set the relative tolerance with.

Returns
-------
numeric
Relative tolerance.
"""

return self._relative_tolerance

@relative_tolerance.setter
def relative_tolerance(self, value):
"""
Setter for the **self.relative_tolerance** property.
"""

if value is not None:
assert is_numeric(value), (
'"relative_tolerance" variable must be a "numeric"!')

self._relative_tolerance = value

@property
def absolute_tolerance(self):
"""
Getter and setter property for the absolute tolerance.

Parameters
----------
value : numeric
Value to set the absolute tolerance with.

Returns
-------
numeric
Absolute tolerance.
"""

return self._absolute_tolerance

@absolute_tolerance.setter
def absolute_tolerance(self, value):
"""
Setter for the **self.absolute_tolerance** property.
"""

if value is not None:
assert is_numeric(value), (
'"absolute_tolerance" variable must be a "numeric"!')

self._absolute_tolerance = value

@property
def default(self):
"""
Getter and setter property for the default value for call outside
tolerances.

Parameters
----------
value : numeric
Value to set the default value with.

Returns
-------
numeric
Default value.
"""

return self._default

@default.setter
def default(self, value):
"""
Setter for the **self.default** property.
"""

if value is not None:
assert is_numeric(value), (
'"default" variable must be a "numeric"!')

self._default = value

[docs]    def __call__(self, x):
"""
Evaluates the interpolator at given point(s).

Parameters
----------
x : numeric or array_like
Point(s) to evaluate the interpolant at.

Returns
-------
float or ndarray
Interpolated value(s).
"""

x = np.atleast_1d(x).astype(self._dtype)

xi = as_float(self._evaluate(x))

return xi

def _evaluate(self, x):
"""
Performs the interpolator evaluation at given points.

Parameters
----------
x : ndarray
Points to evaluate the interpolant at.

Returns
-------
ndarray
Interpolated points values.
"""

self._validate_dimensions()
self._validate_interpolation_range(x)

indexes = closest_indexes(self._x, x)
values = self._y[indexes]
values[~np.isclose(
self._x[indexes],
x,
rtol=self._absolute_tolerance,
atol=self._relative_tolerance)] = self._default

return values

def _validate_dimensions(self):
"""
Validates variables dimensions to be the same.
"""

if len(self._x) != len(self._y):
raise ValueError(
('"x" independent and "y" dependent variables have different '
'dimensions: "{0}", "{1}"').format(
len(self._x), len(self._y)))

def _validate_interpolation_range(self, x):
"""
Validates given point to be in interpolation range.
"""

below_interpolation_range = x < self._x[0]
above_interpolation_range = x > self._x[-1]

if below_interpolation_range.any():
raise ValueError('"{0}" is below interpolation range.'.format(x))

if above_interpolation_range.any():
raise ValueError('"{0}" is above interpolation range.'.format(x))

[docs]def lagrange_coefficients(r, n=4):
"""
Computes the *Lagrange Coefficients* at given point :math:r for degree
:math:n.

Parameters
----------
r : numeric
Point to get the *Lagrange Coefficients* at.
n : int, optional
Degree of the *Lagrange Coefficients* being calculated.

Returns
-------
ndarray

References
----------
:cite:Fairman1985b, :cite:Wikipedia2003a

Examples
--------
>>> lagrange_coefficients(0.1)
array([ 0.8265,  0.2755, -0.1305,  0.0285])
"""

r_i = np.arange(n)
L_n = []
for j in range(len(r_i)):
basis = [(r - r_i[i]) / (r_i[j] - r_i[i]) for i in range(len(r_i))
if i != j]
L_n.append(reduce(lambda x, y: x * y, basis))  # noqa

return np.array(L_n)

def vertices_and_relative_coordinates(V_xyz, table):
"""
Computes the vertices coordinates and indexes relative :math:V_{xyzr}
coordinates from given :math:V_{xyzr} values and interpolation table.

Parameters
----------
V_xyz : array_like
:math:V_{xyz} values to transform to indexes relative
:math:V_{xyzr} values.
table : array_like
4-Dimensional (NxNxNx3) interpolation table.

Returns
-------
tuple
Vertices coordinates and indexes relative :math:V_{xyzr} coordinates.

Examples
--------
>>> import os
>>> import colour
>>> path = os.path.join(
...     os.path.dirname(__file__),'..', 'io', 'luts', 'tests', 'resources',
...     'iridas_cube', 'ColourCorrect.cube')
>>> table = LUT.table
>>> prng = np.random.RandomState(4)
>>> V_xyz = colour.algebra.random_triplet_generator(3, random_state=prng)
>>> print(V_xyz)  # doctest: +ELLIPSIS
[[ 0.9670298...  0.7148159...  0.9762744...]
[ 0.5472322...  0.6977288...  0.0062302...]
[ 0.9726843...  0.2160895...  0.2529823...]]
>>> vertices, V_xyzr = vertices_and_relative_coordinates(V_xyz, table)
>>> print(vertices)
[[[ 0.833311  0.833311  0.833311]
[ 0.349416  0.657749  0.041083]
[ 0.797894 -0.035412 -0.035412]]
<BLANKLINE>
[[ 0.833311  0.833311  1.249963]
[ 0.340435  0.743769  0.340435]
[ 0.752767 -0.028479  0.362144]]
<BLANKLINE>
[[ 0.707102  1.110435  0.707102]
[ 0.344991  1.050213 -0.007621]
[ 0.633333  0.316667  0.      ]]
<BLANKLINE>
[[ 0.519714  0.744729  0.744729]
[ 0.314204  1.120871  0.314204]
[ 0.732278  0.315626  0.315626]]
<BLANKLINE>
[[ 1.06561   0.648957  0.648957]
[ 0.589195  0.589195  0.139164]
[ 1.196841 -0.053117 -0.053117]]
<BLANKLINE>
[[ 1.        0.666667  1.      ]
[ 0.594601  0.594601  0.369586]
[ 1.162588 -0.050372  0.353948]]
<BLANKLINE>
[[ 0.894606  0.894606  0.66959 ]
[ 0.663432  0.930188  0.12992 ]
[ 1.038439  0.310899 -0.05287 ]]
<BLANKLINE>
[[ 1.249966  1.249966  1.249966]
[ 0.682749  0.991082  0.374416]
[ 1.131225  0.29792   0.29792 ]]]
>>> print(V_xyzr)  # doctest: +ELLIPSIS
[[ 0.9010895...  0.1444479...  0.9288233...]
[ 0.6416967...  0.0931864...  0.0186907...]
[ 0.9180530...  0.6482684...  0.7589470...]]
"""

V_xyz = np.clip(V_xyz, 0, 1)
table = as_float_array(table)

V_xyz = np.reshape(V_xyz, (-1, 3))

# Indexes computations where i_m is the maximum index value on a given
# table axis, i_f and i_c respectively the floor and ceiling
# indexes encompassing a given V_xyz value.
i_m = np.array(table.shape[0:-1]) - 1
i_f = np.floor(V_xyz * i_m).astype(DEFAULT_INT_DTYPE)
i_c = np.clip(i_f + 1, 0, i_m)

# Relative to indexes V_xyz values.
V_xyzr = i_m * V_xyz - i_f

i_f_c = i_f, i_c

# Vertices computations by indexing table with the i_f and i_c
# indexes. 8 encompassing vertices are computed for a given V_xyz value
# forming a cube around it:
vertices = np.array([
table[i_f_c[i[0]][..., 0], i_f_c[i[1]][..., 1], i_f_c[i[2]][..., 2]]
for i in itertools.product(*zip([0, 0, 0], [1, 1, 1]))
])

return vertices, V_xyzr

[docs]def table_interpolation_trilinear(V_xyz, table):
"""
Performs trilinear interpolation of given :math:V_{xyz} values using
given interpolation table.

Parameters
----------
V_xyz : array_like
:math:V_{xyz} values to interpolate.
table : array_like
4-Dimensional (NxNxNx3) interpolation table.

Returns
-------
ndarray
Interpolated :math:V_{xyz} values.

References
----------
:cite:Bourkeb

Examples
--------
>>> import os
>>> import colour
>>> path = os.path.join(
...     os.path.dirname(__file__),'..', 'io', 'luts', 'tests', 'resources',
...     'iridas_cube', 'ColourCorrect.cube')
>>> table = LUT.table
>>> prng = np.random.RandomState(4)
>>> V_xyz = colour.algebra.random_triplet_generator(3, random_state=prng)
>>> print(V_xyz)  # doctest: +ELLIPSIS
[[ 0.9670298...  0.7148159...  0.9762744...]
[ 0.5472322...  0.6977288...  0.0062302...]
[ 0.9726843...  0.2160895...  0.2529823...]]
>>> table_interpolation_trilinear(V_xyz, table)  # doctest: +ELLIPSIS
array([[ 1.0120664...,  0.7539146...,  1.0228540...],
[ 0.5075794...,  0.6479459...,  0.1066404...],
[ 1.0976519...,  0.1785998...,  0.2299897...]])
"""

V_xyz = as_float_array(V_xyz)

vertices, V_xyzr = vertices_and_relative_coordinates(V_xyz, table)

vertices = np.moveaxis(vertices, 0, 1)
x, y, z = [f[:, np.newaxis] for f in tsplit(V_xyzr)]

weights = np.moveaxis(
np.transpose(
[(1 - x) * (1 - y) * (1 - z), (1 - x) * (1 - y) * z,
(1 - x) * y * (1 - z), (1 - x) * y * z, x * (1 - y) * (1 - z),
x * (1 - y) * z, x * y * (1 - z), x * y * z]), 0, -1)

xyz_o = np.reshape(np.sum(vertices * weights, 1), V_xyz.shape)

return xyz_o

[docs]def table_interpolation_tetrahedral(V_xyz, table):
"""
Performs tetrahedral interpolation of given :math:V_{xyz} values using
given interpolation table.

Parameters
----------
V_xyz : array_like
:math:V_{xyz} values to interpolate.
table : array_like
4-Dimensional (NxNxNx3) interpolation table.

Returns
-------
ndarray
Interpolated :math:V_{xyz} values.

References
----------
:cite:Kirk2006

Examples
--------
>>> import os
>>> import colour
>>> path = os.path.join(
...     os.path.dirname(__file__),'..', 'io', 'luts', 'tests', 'resources',
...     'iridas_cube', 'ColourCorrect.cube')
>>> table = LUT.table
>>> prng = np.random.RandomState(4)
>>> V_xyz = colour.algebra.random_triplet_generator(3, random_state=prng)
>>> print(V_xyz)  # doctest: +ELLIPSIS
[[ 0.9670298...  0.7148159...  0.9762744...]
[ 0.5472322...  0.6977288...  0.0062302...]
[ 0.9726843...  0.2160895...  0.2529823...]]
>>> table_interpolation_tetrahedral(V_xyz, table)  # doctest: +ELLIPSIS
array([[ 1.0196197...,  0.7674062...,  1.0311751...],
[ 0.5105603...,  0.6466722...,  0.1077296...],
[ 1.1178206...,  0.1762039...,  0.2209534...]])
"""

V_xyz = as_float_array(V_xyz)

vertices, V_xyzr = vertices_and_relative_coordinates(V_xyz, table)

vertices = np.moveaxis(vertices, 0, -1)
V000, V001, V010, V011, V100, V101, V110, V111 = tsplit(vertices)
x, y, z = [r[:, np.newaxis] for r in tsplit(V_xyzr)]

xyz_o = np.select([
np.logical_and(x > y, y > z),
np.logical_and(x > y, x > z),
np.logical_and(x > y, np.logical_and(y <= z, x <= z)),
np.logical_and(x <= y, z > y),
np.logical_and(x <= y, z > x),
np.logical_and(x <= y, np.logical_and(z <= y, z <= x)),
], [
(1 - x) * V000 + (x - y) * V100 + (y - z) * V110 + z * V111,
(1 - x) * V000 + (x - z) * V100 + (z - y) * V101 + y * V111,
(1 - z) * V000 + (z - x) * V001 + (x - y) * V101 + y * V111,
(1 - z) * V000 + (z - y) * V001 + (y - x) * V011 + x * V111,
(1 - y) * V000 + (y - z) * V010 + (z - x) * V011 + x * V111,
(1 - y) * V000 + (y - x) * V010 + (x - z) * V110 + z * V111,
])

xyz_o = np.reshape(xyz_o, V_xyz.shape)

return xyz_o

TABLE_INTERPOLATION_METHODS = CaseInsensitiveMapping({
'Trilinear': table_interpolation_trilinear,
'Tetrahedral': table_interpolation_tetrahedral,
})
TABLE_INTERPOLATION_METHODS.__doc__ = """
Supported table interpolation methods.

References
----------
:cite:Bourkeb, :cite:Kirk2006

TABLE_INTERPOLATION_METHODS : CaseInsensitiveMapping
**{'Trilinear', 'Tetrahedral'}**
"""

[docs]def table_interpolation(V_xyz, table, method='Trilinear'):
"""
Performs interpolation of given :math:V_{xyz} values using given
interpolation table.

Parameters
----------
V_xyz : array_like
:math:V_{xyz} values to interpolate.
table : array_like
4-Dimensional (NxNxNx3) interpolation table.
method : unicode, optional
**{'Trilinear', 'Tetrahedral'}**,
Interpolation method.

Returns
-------
ndarray
Interpolated :math:V_{xyz} values.

References
----------
:cite:Bourkeb, :cite:Kirk2006

Examples
--------
>>> import os
>>> import colour
>>> path = os.path.join(
...     os.path.dirname(__file__),'..', 'io', 'luts', 'tests', 'resources',
...     'iridas_cube', 'ColourCorrect.cube')
>>> table = LUT.table
>>> prng = np.random.RandomState(4)
>>> V_xyz = colour.algebra.random_triplet_generator(3, random_state=prng)
>>> print(V_xyz)  # doctest: +ELLIPSIS
[[ 0.9670298...  0.7148159...  0.9762744...]
[ 0.5472322...  0.6977288...  0.0062302...]
[ 0.9726843...  0.2160895...  0.2529823...]]
>>> table_interpolation(V_xyz, table)  # doctest: +ELLIPSIS
array([[ 1.0120664...,  0.7539146...,  1.0228540...],
[ 0.5075794...,  0.6479459...,  0.1066404...],
[ 1.0976519...,  0.1785998...,  0.2299897...]])
>>> table_interpolation(V_xyz, table, method='Tetrahedral')
... # doctest: +ELLIPSIS
array([[ 1.0196197...,  0.7674062...,  1.0311751...],
[ 0.5105603...,  0.6466722...,  0.1077296...],
[ 1.1178206...,  0.1762039...,  0.2209534...]])
"""

return TABLE_INTERPOLATION_METHODS.get(method)(V_xyz, table)