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.
    Retrieved from https://www.filmlight.ltd.uk/pdf/whitepapers/\
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'
__copyright__ = 'Copyright (C) 2013-2019 - Colour Developers'
__license__ = 'New BSD License - https://opensource.org/licenses/BSD-3-Clause'
__maintainer__ = 'Colour Developers'
__email__ = 'colour-science@googlegroups.com'
__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. padding_args : dict, optional 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 padding_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, padding_args=None, dtype=DEFAULT_FLOAT_DTYPE): self._x_p = None self._y_p = None self._x = None self._y = None self._window = None self._padding_args = {'pad_width': (window, window), 'mode': 'reflect'} self._dtype = dtype self.x = x self.y = y self.window = window self.padding_args = padding_args 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_p = np.pad( 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: self._y_p = np.pad(self._y, **self._padding_args) @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 def padding_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._padding_args @padding_args.setter def padding_args(self, value): """ Setter for the **self.padding_args** property. """ if value is not None: assert isinstance(value, Mapping), ( '"{0}" attribute: "{1}" type is not a "Mapping" instance!' ).format('padding_args', value) self._padding_args = value # 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. padding_args : dict, optional 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', 'Colour_Correct.cube') >>> LUT = colour.read_LUT(path) >>> 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', 'Colour_Correct.cube') >>> LUT = colour.read_LUT(path) >>> 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', 'Colour_Correct.cube') >>> LUT = colour.read_LUT(path) >>> 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', 'Colour_Correct.cube') >>> LUT = colour.read_LUT(path) >>> 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)