Repository: colour-science/colour

Issue:
Implement support for spectral data "Piecewise Cubic Hermite Interpolation".

Discussion:

Patch:
diff --git a/colour/algebra/__init__.py b/colour/algebra/__init__.py
index f46d14502..77bfc6f5e 100644
--- a/colour/algebra/__init__.py
+++ b/colour/algebra/__init__.py
@@ -6,17 +6,19 @@ from __future__ import absolute_import
 from .coordinates import *  # noqa
 from . import coordinates
 from .extrapolation import Extrapolator1d
-from .interpolation import (LinearInterpolator1d,
-                            SplineInterpolator,
-                            SpragueInterpolator)
+from .interpolation import (LinearInterpolator,
+                            SpragueInterpolator,
+                            CubicSplineInterpolator,
+                            PchipInterpolator)
 from .matrix import is_identity
 from .random import random_triplet_generator
 
 __all__ = []
 __all__ += coordinates.__all__
 __all__ += ['Extrapolator1d']
-__all__ += ['LinearInterpolator1d',
-            'SplineInterpolator',
-            'SpragueInterpolator']
+__all__ += ['LinearInterpolator',
+            'SpragueInterpolator',
+            'CubicSplineInterpolator',
+            'PchipInterpolator']
 __all__ += ['is_identity']
 __all__ += ['random_triplet_generator']
diff --git a/colour/algebra/extrapolation.py b/colour/algebra/extrapolation.py
index 148074866..3e60ad80a 100644
--- a/colour/algebra/extrapolation.py
+++ b/colour/algebra/extrapolation.py
@@ -74,10 +74,10 @@ class Extrapolator1d(object):
     --------
     Extrapolating a single numeric variable:
 
-    >>> from colour.algebra import LinearInterpolator1d
+    >>> from colour.algebra import LinearInterpolator
     >>> x = np.array([3, 4, 5])
     >>> y = np.array([1, 2, 3])
-    >>> interpolator = LinearInterpolator1d(x, y)
+    >>> interpolator = LinearInterpolator(x, y)
     >>> extrapolator = Extrapolator1d(interpolator)
     >>> extrapolator(1)
     -1.0
@@ -91,7 +91,7 @@ class Extrapolator1d(object):
 
     >>> x = np.array([3, 4, 5])
     >>> y = np.array([1, 2, 3])
-    >>> interpolator = LinearInterpolator1d(x, y)
+    >>> interpolator = LinearInterpolator(x, y)
     >>> extrapolator = Extrapolator1d(interpolator, method='Constant')
     >>> extrapolator(np.array([0.1, 0.2, 8, 9]))
     array([ 1.,  1.,  3.,  3.])
@@ -100,7 +100,7 @@ class Extrapolator1d(object):
 
     >>> x = np.array([3, 4, 5])
     >>> y = np.array([1, 2, 3])
-    >>> interpolator = LinearInterpolator1d(x, y)
+    >>> interpolator = LinearInterpolator(x, y)
     >>> extrapolator = Extrapolator1d(interpolator, method='Constant', left=0)
     >>> extrapolator(np.array([0.1, 0.2, 8, 9]))
     array([ 0.,  0.,  3.,  3.])
diff --git a/colour/algebra/interpolation.py b/colour/algebra/interpolation.py
index 83aa98bb4..4830120da 100644
--- a/colour/algebra/interpolation.py
+++ b/colour/algebra/interpolation.py
@@ -7,10 +7,12 @@ Interpolation
 
 Defines classes for interpolating variables.
 
--   :class:`LinearInterpolator1d`: 1-D function linear interpolation.
--   :class:`SplineInterpolator`: 1-D function cubic spline interpolation.
+-   :class:`LinearInterpolator`: 1-D function linear interpolation.
 -   :class:`SpragueInterpolator`: 1-D function fifth-order polynomial
     interpolation.
+-   :class:`CubicSplineInterpolator`: 1-D function cubic spline interpolation.
+-   :class:`PchipInterpolator`: 1-D function piecewise cube Hermite
+    interpolation.
 """
 
 from __future__ import division, unicode_literals
@@ -32,12 +34,13 @@ __maintainer__ = 'Colour Developers'
 __email__ = 'colour-science@googlegroups.com'
 __status__ = 'Production'
 
-__all__ = ['LinearInterpolator1d',
-           'SplineInterpolator',
-           'SpragueInterpolator']
+__all__ = ['LinearInterpolator',
+           'SpragueInterpolator',
+           'CubicSplineInterpolator',
+           'PchipInterpolator']
 
 
-class LinearInterpolator1d(object):
+class LinearInterpolator(object):
     """
     Linearly interpolates a 1-D function.
 
@@ -74,7 +77,7 @@ class LinearInterpolator1d(object):
     ...               27.8007,
     ...               86.0500])
     >>> x = np.arange(len(y))
-    >>> f = LinearInterpolator1d(x, y)
+    >>> f = LinearInterpolator(x, y)
     >>> # Doctests ellipsis for Python 2.x compatibility.
     >>> f(0.5)  # doctest: +ELLIPSIS
     7.64...
@@ -225,31 +228,6 @@ class LinearInterpolator1d(object):
             raise ValueError('"{0}" is above interpolation range.'.format(x))
 
 
-if is_scipy_installed():
-    from scipy.interpolate import interp1d
-
-    class SplineInterpolator(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):
-            # TODO: Implements proper wrapper to ensure return values
-            # consistency and avoid having to cast to numeric in
-            # :meth:`SpectralPowerDistribution.interpolate` method.
-            super(SplineInterpolator, self).__init__(
-                kind='cubic', *args, **kwargs)
-else:
-    warning(('"scipy.interpolate.interp1d" interpolator is unavailable, using '
-             '"LinearInterpolator1d" instead!'))
-
-    SplineInterpolator = LinearInterpolator1d
-
-
 class SpragueInterpolator(object):
     """
     Constructs a fifth-order polynomial that passes through :math:`y` dependent
@@ -273,7 +251,7 @@ class SpragueInterpolator(object):
 
     See Also
     --------
-    LinearInterpolator1d
+    LinearInterpolator
 
     Notes
     -----
@@ -523,3 +501,29 @@ class SpragueInterpolator(object):
 
         if above_interpolation_range.any():
             raise ValueError('"{0}" is above interpolation range.'.format(x))
+
+
+if is_scipy_installed():
+    from scipy.interpolate import PchipInterpolator, interp1d
+
+    class CubicSplineInterpolator(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):
+            # TODO: Implements proper wrapper to ensure return values
+            # consistency and avoid having to cast to numeric in
+            # :meth:`SpectralPowerDistribution.interpolate` method.
+            super(CubicSplineInterpolator, self).__init__(
+                kind='cubic', *args, **kwargs)
+else:
+    warning(('"scipy.interpolate.PchipInterpolator" and '
+             '"scipy.interpolate.interp1d" interpolators are not available, '
+             'using "LinearInterpolator" instead!'))
+
+    PchipInterpolator = CubicSplineInterpolator = LinearInterpolator
diff --git a/colour/colorimetry/spectrum.py b/colour/colorimetry/spectrum.py
index 37412d73f..bff22dd7d 100644
--- a/colour/colorimetry/spectrum.py
+++ b/colour/colorimetry/spectrum.py
@@ -26,9 +26,10 @@ import numpy as np
 
 from colour.algebra import (
     Extrapolator1d,
-    LinearInterpolator1d,
-    SplineInterpolator,
-    SpragueInterpolator)
+    LinearInterpolator,
+    SpragueInterpolator,
+    CubicSplineInterpolator,
+    PchipInterpolator)
 from colour.utilities import (
     ArbitraryPrecisionMapping,
     is_iterable,
@@ -1587,11 +1588,8 @@ class SpectralPowerDistribution(object):
         """
 
         extrapolator = Extrapolator1d(
-            LinearInterpolator1d(self.wavelengths,
-                                 self.values),
-            method=method,
-            left=left,
-            right=right)
+            LinearInterpolator(self.wavelengths, self.values),
+            method=method, left=left, right=right)
 
         spd_shape = self.shape
         for i in np.arange(spd_shape.start,
@@ -1618,7 +1616,7 @@ class SpectralPowerDistribution(object):
         shape : SpectralShape, optional
             Spectral shape used for interpolation.
         method : unicode, optional
-            {None, 'Sprague', 'Cubic Spline', 'Linear'},
+            {None, 'Cubic Spline', 'Linear', 'Pchip', 'Sprague'},
             Enforce given interpolation method.
 
         Returns
@@ -1651,9 +1649,11 @@ class SpectralPowerDistribution(object):
         -------
         -   If *scipy* is not unavailable the *Cubic Spline* method will
             fallback to legacy *Linear* interpolation.
+        -   *Cubic Spline* interpolator requires at least 3 wavelengths
+            :math:`\lambda_n` for interpolation.
         -   *Linear* interpolator requires at least 2 wavelengths
             :math:`\lambda_n` for interpolation.
-        -   *Cubic Spline* interpolator requires at least 3 wavelengths
+        -   *Pchip* interpolator requires at least 2 wavelengths
             :math:`\lambda_n` for interpolation.
         -   Sprague (1880) interpolator requires at least 6 wavelengths
             :math:`\lambda_n` for interpolation.
@@ -1683,13 +1683,6 @@ class SpectralPowerDistribution(object):
 
         Non uniform data is using *Cubic Spline* interpolation by default:
 
-        >>> data = {
-        ...     510: 49.67,
-        ...     520: 69.59,
-        ...     530: 81.73,
-        ...     540: 88.19,
-        ...     550: 86.26,
-        ...     560: 77.18}
         >>> spd = SpectralPowerDistribution('Spd', data)
         >>> spd[511] = 31.41
         >>> spd.interpolate(SpectralShape(steps=1))  # doctest: +ELLIPSIS
@@ -1699,18 +1692,23 @@ class SpectralPowerDistribution(object):
 
         Enforcing *Linear* interpolation:
 
-        >>> data = {
-        ...     510: 49.67,
-        ...     520: 69.59,
-        ...     530: 81.73,
-        ...     540: 88.19,
-        ...     550: 86.26,
-        ...     560: 77.18}
         >>> spd = SpectralPowerDistribution('Spd', data)
-        >>> spd.interpolate(SpectralShape(steps=1), method='Linear')  # noqa  # doctest: +ELLIPSIS
+        >>> spd.interpolate(  # doctest: +ELLIPSIS
+        ...     SpectralShape(steps=1),
+        ...     method='Linear')
         <...SpectralPowerDistribution object at 0x...>
         >>> spd[515]  # doctest: +ELLIPSIS
         array(59.63...)
+
+        Enforcing *Pchip* interpolation:
+
+        >>> spd = SpectralPowerDistribution('Spd', data)
+        >>> spd.interpolate(  # doctest: +ELLIPSIS
+        ...     SpectralShape(steps=1),
+        ...     method='Pchip')
+        <...SpectralPowerDistribution object at 0x...>
+        >>> spd[515]  # doctest: +ELLIPSIS
+        array(58.8173260...)
         """
 
         spd_shape = self.shape
@@ -1732,25 +1730,28 @@ class SpectralPowerDistribution(object):
 
         if method is None:
             if is_uniform:
-                interpolator = SpragueInterpolator(wavelengths, values)
+                interpolator = SpragueInterpolator
             else:
-                interpolator = SplineInterpolator(wavelengths, values)
+                interpolator = CubicSplineInterpolator
+        elif method == 'cubic spline':
+            interpolator = CubicSplineInterpolator
+        elif method == 'linear':
+            interpolator = LinearInterpolator
+        elif method == 'pchip':
+            interpolator = PchipInterpolator
         elif method == 'sprague':
             if is_uniform:
-                interpolator = SpragueInterpolator(wavelengths, values)
+                interpolator = SpragueInterpolator
             else:
                 raise RuntimeError(
                     ('"Sprague" interpolator can only be used for '
                      'interpolating functions having a uniformly spaced '
                      'independent variable!'))
-        elif method == 'cubic spline':
-            interpolator = SplineInterpolator(wavelengths, values)
-        elif method == 'linear':
-            interpolator = LinearInterpolator1d(wavelengths, values)
         else:
             raise ValueError(
                 'Undefined "{0}" interpolator!'.format(method))
 
+        interpolator = interpolator(wavelengths, values)
         self.__data = SpectralMapping(
             [(wavelength, float(interpolator(wavelength)))
              for wavelength in shape])
@@ -3351,7 +3352,7 @@ class TriSpectralPowerDistribution(object):
         shape : SpectralShape, optional
             Spectral shape used for interpolation.
         method : unicode, optional
-            {None, 'Sprague', 'Cubic Spline', 'Linear'},
+            {None, 'Cubic Spline', 'Linear', 'Pchip', 'Sprague'},
             Enforce given interpolation method.
 
         Returns
@@ -3403,71 +3404,43 @@ class TriSpectralPowerDistribution(object):
         >>> tri_spd = TriSpectralPowerDistribution('Tri Spd', data, mapping)
         >>> tri_spd.interpolate(SpectralShape(steps=1))  # doctest: +ELLIPSIS
         <...TriSpectralPowerDistribution object at 0x...>
-        >>> tri_spd[515]
-        array([ 60.30332087,  93.27163315,  13.86051361])
+        >>> tri_spd[515]  # doctest: +ELLIPSIS
+        array([ 60.3033208...,  93.2716331...,  13.8605136...])
 
         Non uniform data is using *Cubic Spline* interpolation by default:
 
-        >>> x_bar = {
-        ...     510: 49.67,
-        ...     520: 69.59,
-        ...     530: 81.73,
-        ...     540: 88.19,
-        ...     550: 89.76,
-        ...     560: 90.28}
-        >>> y_bar = {
-        ...     510: 90.56,
-        ...     520: 87.34,
-        ...     530: 45.76,
-        ...     540: 23.45,
-        ...     550: 15.34,
-        ...     560: 10.11}
-        >>> z_bar = {
-        ...     510: 12.43,
-        ...     520: 23.15,
-        ...     530: 67.98,
-        ...     540: 90.28,
-        ...     550: 91.61,
-        ...     560: 98.24}
         >>> data = {'x_bar': x_bar, 'y_bar': y_bar, 'z_bar': z_bar}
         >>> mapping = {'x': 'x_bar', 'y': 'y_bar', 'z': 'z_bar'}
         >>> tri_spd = TriSpectralPowerDistribution('Tri Spd', data, mapping)
         >>> tri_spd[511] = np.array([31.41, 95.27, 15.06])
         >>> tri_spd.interpolate(SpectralShape(steps=1))  # doctest: +ELLIPSIS
         <...TriSpectralPowerDistribution object at 0x...>
-        >>> tri_spd[515]
-        array([  21.47104053,  100.64300155,   18.8165196 ])
+        >>> tri_spd[515]  # doctest: +ELLIPSIS
+        array([  21.4710405...,  100.6430015...,   18.8165196...])
 
         Enforcing *Linear* interpolation:
 
-        >>> x_bar = {
-        ...     510: 49.67,
-        ...     520: 69.59,
-        ...     530: 81.73,
-        ...     540: 88.19,
-        ...     550: 89.76,
-        ...     560: 90.28}
-        >>> y_bar = {
-        ...     510: 90.56,
-        ...     520: 87.34,
-        ...     530: 45.76,
-        ...     540: 23.45,
-        ...     550: 15.34,
-        ...     560: 10.11}
-        >>> z_bar = {
-        ...     510: 12.43,
-        ...     520: 23.15,
-        ...     530: 67.98,
-        ...     540: 90.28,
-        ...     550: 91.61,
-        ...     560: 98.24}
         >>> data = {'x_bar': x_bar, 'y_bar': y_bar, 'z_bar': z_bar}
         >>> mapping = {'x': 'x_bar', 'y': 'y_bar', 'z': 'z_bar'}
         >>> tri_spd = TriSpectralPowerDistribution('Tri Spd', data, mapping)
-        >>> tri_spd.interpolate(SpectralShape(steps=1), method='Linear')  # noqa  # doctest: +ELLIPSIS
+        >>> tri_spd.interpolate(  # doctest: +ELLIPSIS
+        ...     SpectralShape(steps=1),
+        ...     method='Linear')
+        <...TriSpectralPowerDistribution object at 0x...>
+        >>> tri_spd[515]  # doctest: +ELLIPSIS
+        array([ 59.63...,  88.95...,  17.79...])
+
+        Enforcing *Pchip* interpolation:
+
+        >>> data = {'x_bar': x_bar, 'y_bar': y_bar, 'z_bar': z_bar}
+        >>> mapping = {'x': 'x_bar', 'y': 'y_bar', 'z': 'z_bar'}
+        >>> tri_spd = TriSpectralPowerDistribution('Tri Spd', data, mapping)
+        >>> tri_spd.interpolate(  # doctest: +ELLIPSIS
+        ...     SpectralShape(steps=1),
+        ...     method='Pchip')
         <...TriSpectralPowerDistribution object at 0x...>
-        >>> tri_spd[515]
-        array([ 59.63,  88.95,  17.79])
+        >>> tri_spd[515]  # doctest: +ELLIPSIS
+        array([ 58.8173260...,  89.4355596...,  16.4545683...])
         """
 
         for i in self.__mapping.keys():
diff --git a/colour/colorimetry/tristimulus.py b/colour/colorimetry/tristimulus.py
index e3d3eaa32..a901fc0cc 100644
--- a/colour/colorimetry/tristimulus.py
+++ b/colour/colorimetry/tristimulus.py
@@ -19,8 +19,13 @@ from __future__ import division, unicode_literals
 
 import numpy as np
 
-from colour.algebra import SplineInterpolator, SpragueInterpolator
+from colour.algebra import (
+    CubicSplineInterpolator,
+    LinearInterpolator,
+    PchipInterpolator,
+    SpragueInterpolator)
 from colour.colorimetry import STANDARD_OBSERVERS_CMFS, ones_spd
+from colour.utilities import is_string
 
 __author__ = 'Colour Developers'
 __copyright__ = 'Copyright (C) 2013 - 2015 - Colour Developers'
@@ -111,7 +116,8 @@ def spectral_to_XYZ(spd,
 
 def wavelength_to_XYZ(wavelength,
                       cmfs=STANDARD_OBSERVERS_CMFS.get(
-                          'CIE 1931 2 Degree Standard Observer')):
+                          'CIE 1931 2 Degree Standard Observer'),
+                      method=None):
     """
     Converts given wavelength :math:`\lambda` to *CIE XYZ* tristimulus values
     using given colour matching functions.
@@ -128,6 +134,9 @@ def wavelength_to_XYZ(wavelength,
         Wavelength :math:`\lambda` in nm.
     cmfs : XYZ_ColourMatchingFunctions, optional
         Standard observer colour matching functions.
+    method : unicode, optional
+        {None, 'Cubic Spline', 'Linear', 'Pchip', 'Sprague'},
+        Enforce given interpolation method.
 
     Returns
     -------
@@ -136,22 +145,60 @@ def wavelength_to_XYZ(wavelength,
 
     Raises
     ------
+    RuntimeError
+        If Sprague (1880) interpolation method is forced with a
+        non-uniformly spaced independent variable.
     ValueError
-        If wavelength :math:`\lambda` is not contained in the colour matching
-        functions domain.
+        If the interpolation method is not defined or if wavelength
+        :math:`\lambda` is not contained in the colour matching functions
+        domain.
 
     Notes
     -----
     -   Output *CIE XYZ* tristimulus values are in domain [0, 1].
     -   If *scipy* is not unavailable the *Cubic Spline* method will fallback
         to legacy *Linear* interpolation.
+    -   Sprague (1880) interpolator cannot be used for interpolating
+        functions having a non-uniformly spaced independent variable.
+
+    Warning
+    -------
+    -   If *scipy* is not unavailable the *Cubic Spline* method will fallback
+        to legacy *Linear* interpolation.
+    -   *Cubic Spline* interpolator requires at least 3 wavelengths
+        :math:`\lambda_n` for interpolation.
+    -   *Linear* interpolator requires at least 2 wavelengths :math:`\lambda_n`
+        for interpolation.
+    -   *Pchip* interpolator requires at least 2 wavelengths :math:`\lambda_n`
+        for interpolation.
+    -   Sprague (1880) interpolator requires at least 6 wavelengths
+        :math:`\lambda_n` for interpolation.
 
     Examples
     --------
+    Uniform data is using Sprague (1880) interpolation by default:
+
     >>> from colour import CMFS
     >>> cmfs = CMFS.get('CIE 1931 2 Degree Standard Observer')
-    >>> wavelength_to_XYZ(480)  # doctest: +ELLIPSIS
+    >>> wavelength_to_XYZ(480, cmfs)  # doctest: +ELLIPSIS
     array([ 0.09564  ,  0.13902  ,  0.812950...])
+    >>> wavelength_to_XYZ(480.5, cmfs)  # doctest: +ELLIPSIS
+    array([ 0.0914287...,  0.1418350...,  0.7915726...])
+
+    Enforcing *Cubic Spline* interpolation:
+
+    >>> wavelength_to_XYZ(480.5, cmfs, 'Cubic Spline')  # doctest: +ELLIPSIS
+    array([ 0.0914288...,  0.1418351...,  0.7915729...])
+
+    Enforcing *Linear* interpolation:
+
+    >>> wavelength_to_XYZ(480.5, cmfs, 'Linear')  # doctest: +ELLIPSIS
+    array([ 0.0914697...,  0.1418482...,  0.7917337...])
+
+    Enforcing *Pchip* interpolation:
+
+    >>> wavelength_to_XYZ(480.5, cmfs, 'Pchip')  # doctest: +ELLIPSIS
+    array([ 0.0914280...,  0.1418341...,  0.7915711...])
     """
 
     cmfs_shape = cmfs.shape
@@ -163,9 +210,34 @@ def wavelength_to_XYZ(wavelength,
 
     if wavelength not in cmfs:
         wavelengths, values, = cmfs.wavelengths, cmfs.values
-        interpolator = (SpragueInterpolator
-                        if cmfs.is_uniform() else
-                        SplineInterpolator)
+
+        if is_string(method):
+            method = method.lower()
+
+        is_uniform = cmfs.is_uniform()
+
+        if method is None:
+            if is_uniform:
+                interpolator = SpragueInterpolator
+            else:
+                interpolator = CubicSplineInterpolator
+        elif method == 'cubic spline':
+            interpolator = CubicSplineInterpolator
+        elif method == 'linear':
+            interpolator = LinearInterpolator
+        elif method == 'pchip':
+            interpolator = PchipInterpolator
+        elif method == 'sprague':
+            if is_uniform:
+                interpolator = SpragueInterpolator
+            else:
+                raise RuntimeError(
+                    ('"Sprague" interpolator can only be used for '
+                     'interpolating functions having a uniformly spaced '
+                     'independent variable!'))
+        else:
+            raise ValueError(
+                'Undefined "{0}" interpolator!'.format(method))
 
         interpolators = [interpolator(wavelengths, values[..., i])
                          for i in range(values.shape[-1])]
diff --git a/colour/examples/algebra/examples_interpolation.py b/colour/examples/algebra/examples_interpolation.py
index 9f94f50e4..d2edba760 100644
--- a/colour/examples/algebra/examples_interpolation.py
+++ b/colour/examples/algebra/examples_interpolation.py
@@ -16,7 +16,7 @@ from colour.utilities.verbose import message_box
 message_box('Interpolation Computations')
 
 message_box(('Comparing Sprague (1880) and "Cubic Spline" recommended '
-             'interpolation methods.'))
+             'interpolation methods to "Pchip" method.'))
 
 uniform_spd_data = {
     340: 0.0000,
@@ -78,11 +78,16 @@ base_spd = colour.SpectralPowerDistribution(
 uniform_interpolated_spd = colour.SpectralPowerDistribution(
     'Uniform - Sprague Interpolation',
     uniform_spd_data)
+uniform_pchip_interpolated_spd = colour.SpectralPowerDistribution(
+    'Uniform - Pchip Interpolation',
+    uniform_spd_data)
 non_uniform_interpolated_spd = colour.SpectralPowerDistribution(
     'Non Uniform - Cubic Spline Interpolation',
     non_uniform_spd_data)
 
 uniform_interpolated_spd.interpolate(colour.SpectralShape(steps=1))
+uniform_pchip_interpolated_spd.interpolate(colour.SpectralShape(steps=1),
+                                           method='Pchip')
 non_uniform_interpolated_spd.interpolate(colour.SpectralShape(steps=1))
 
 shape = base_spd.shape
@@ -97,6 +102,10 @@ pylab.plot(uniform_interpolated_spd.wavelengths,
            uniform_interpolated_spd.values,
            label=uniform_interpolated_spd.name,
            linewidth=2)
+pylab.plot(uniform_pchip_interpolated_spd.wavelengths,
+           uniform_pchip_interpolated_spd.values,
+           label=uniform_pchip_interpolated_spd.name,
+           linewidth=2)
 pylab.plot(non_uniform_interpolated_spd.wavelengths,
            non_uniform_interpolated_spd.values,
            label=non_uniform_interpolated_spd.name,
diff --git a/colour/notation/munsell.py b/colour/notation/munsell.py
index 9a39cc6f7..93106af11 100644
--- a/colour/notation/munsell.py
+++ b/colour/notation/munsell.py
@@ -53,7 +53,7 @@ except ImportError:
 
 from colour.algebra import (
     Extrapolator1d,
-    LinearInterpolator1d,
+    LinearInterpolator,
     cartesian_to_cylindrical)
 from colour.colorimetry import ILLUMINANTS, luminance_ASTMD153508
 from colour.constants import (
@@ -197,7 +197,7 @@ def _munsell_value_ASTMD153508_interpolator():
     munsell_values = np.arange(0, 10, 0.001)
     if _MUNSELL_VALUE_ASTM_D1535_08_INTERPOLATOR_CACHE is None:
         _MUNSELL_VALUE_ASTM_D1535_08_INTERPOLATOR_CACHE = Extrapolator1d(
-            LinearInterpolator1d(
+            LinearInterpolator(
                 luminance_ASTMD153508(munsell_values),
                 munsell_values))
 
@@ -624,8 +624,8 @@ def munsell_specification_to_xyY(specification):
     else:
         Y_minus = luminance_ASTMD153508(value_minus)
         Y_plus = luminance_ASTMD153508(value_plus)
-        x = LinearInterpolator1d((Y_minus, Y_plus), (x_minus, x_plus))(Y)
-        y = LinearInterpolator1d((Y_minus, Y_plus), (y_minus, y_plus))(Y)
+        x = LinearInterpolator((Y_minus, Y_plus), (x_minus, x_plus))(Y)
+        y = LinearInterpolator((Y_minus, Y_plus), (y_minus, y_plus))(Y)
 
     return np.array([x, y, Y / 100])
 
@@ -826,7 +826,7 @@ def xyY_to_munsell_specification(xyY):
             theta_differences_indexes]
 
         hue_angle_difference_new = Extrapolator1d(
-            LinearInterpolator1d(
+            LinearInterpolator(
                 theta_differences,
                 hue_angles_differences))(0) % 360
         hue_angle_new = (hue_angle_current + hue_angle_difference_new) % 360
@@ -891,7 +891,7 @@ def xyY_to_munsell_specification(xyY):
 
         rho_bounds = rho_bounds[rhos_bounds_indexes]
         chroma_bounds = chroma_bounds[rhos_bounds_indexes]
-        chroma_new = LinearInterpolator1d(rho_bounds, chroma_bounds)(rho_input)
+        chroma_new = LinearInterpolator(rho_bounds, chroma_bounds)(rho_input)
 
         specification_current = [hue_current, value, chroma_new, code_current]
         x_current, y_current, Y_current = np.ravel(
@@ -1305,7 +1305,7 @@ def hue_to_hue_angle(hue, code):
     """
 
     single_hue = ((17 - code) % 10 + (hue / 10) - 0.5) % 10
-    return LinearInterpolator1d(
+    return LinearInterpolator(
         (0, 2, 3, 4, 5, 6, 8, 9, 10),
         (0, 45, 70, 135, 160, 225, 255, 315, 360))(single_hue)
 
@@ -1339,7 +1339,7 @@ def hue_angle_to_hue(hue_angle):
     (3.2160000..., 4)
     """
 
-    single_hue = LinearInterpolator1d(
+    single_hue = LinearInterpolator(
         (0, 45, 70, 135, 160, 225, 255, 315, 360),
         (0, 2, 3, 4, 5, 6, 8, 9, 10))(hue_angle)
 
@@ -1826,15 +1826,15 @@ def xy_from_renotation_ovoid(specification):
             specification).lower()
 
         if interpolation_method == 'linear':
-            x = LinearInterpolator1d((lower_hue_angle, upper_hue_angle),
-                                     (x_minus, x_plus))(hue_angle)
-            y = LinearInterpolator1d((lower_hue_angle, upper_hue_angle),
-                                     (y_minus, y_plus))(hue_angle)
+            x = LinearInterpolator((lower_hue_angle, upper_hue_angle),
+                                   (x_minus, x_plus))(hue_angle)
+            y = LinearInterpolator((lower_hue_angle, upper_hue_angle),
+                                   (y_minus, y_plus))(hue_angle)
         elif interpolation_method == 'radial':
-            theta = LinearInterpolator1d((lower_hue_angle, upper_hue_angle),
-                                         (theta_minus, theta_plus))(hue_angle)
-            rho = LinearInterpolator1d((lower_hue_angle, upper_hue_angle),
-                                       (rho_minus, rho_plus))(hue_angle)
+            theta = LinearInterpolator((lower_hue_angle, upper_hue_angle),
+                                       (theta_minus, theta_plus))(hue_angle)
+            rho = LinearInterpolator((lower_hue_angle, upper_hue_angle),
+                                     (rho_minus, rho_plus))(hue_angle)
 
             x = rho * np.cos(np.radians(theta)) + x_grey
             y = rho * np.sin(np.radians(theta)) + y_grey
@@ -1903,7 +1903,7 @@ def LCHab_to_munsell_specification(LCHab):
     else:
         code = 8
 
-    hue = LinearInterpolator1d((0, 36), (0, 10))(Hab % 36)
+    hue = LinearInterpolator((0, 36), (0, 10))(Hab % 36)
     if hue == 0:
         hue = 10
 
@@ -1986,8 +1986,8 @@ def maximum_chroma_from_renotation(hue, value, code):
         L9 = luminance_ASTMD153508(9)
         L10 = luminance_ASTMD153508(10)
 
-        max_chroma = min(LinearInterpolator1d((L9, L10), (ma_limit_mcw, 0))(L),
-                         LinearInterpolator1d((L9, L10), (ma_limit_mccw, 0))(
+        max_chroma = min(LinearInterpolator((L9, L10), (ma_limit_mcw, 0))(L),
+                         LinearInterpolator((L9, L10), (ma_limit_mccw, 0))(
                              L))
     return max_chroma
 
@@ -2065,9 +2065,9 @@ def munsell_specification_to_xy(specification):
             x = x_minus
             y = y_minus
         else:
-            x = LinearInterpolator1d((chroma_minus, chroma_plus),
-                                     (x_minus, x_plus))(chroma)
-            y = LinearInterpolator1d((chroma_minus, chroma_plus),
-                                     (y_minus, y_plus))(chroma)
+            x = LinearInterpolator((chroma_minus, chroma_plus),
+                                   (x_minus, x_plus))(chroma)
+            y = LinearInterpolator((chroma_minus, chroma_plus),
+                                   (y_minus, y_plus))(chroma)
 
         return x, y
diff --git a/colour/plotting/__init__.py b/colour/plotting/__init__.py
index 04dbb252b..0686b6164 100644
--- a/colour/plotting/__init__.py
+++ b/colour/plotting/__init__.py
@@ -12,9 +12,9 @@ from .common import (
     DEFAULT_FIGURE_HEIGHT,
     DEFAULT_FIGURE_SIZE,
     DEFAULT_FONT_SIZE,
-    DEFAULT_PARAMETERS,
     DEFAULT_COLOUR_CYCLE,
     DEFAULT_HATCH_PATTERNS,
+    DEFAULT_PARAMETERS,
     DEFAULT_PLOTTING_ILLUMINANT,
     DEFAULT_PLOTTING_OECF,
     ColourParameter,
@@ -91,9 +91,9 @@ __all__ += [
     'DEFAULT_FIGURE_HEIGHT',
     'DEFAULT_FIGURE_SIZE',
     'DEFAULT_FONT_SIZE',
-    'DEFAULT_PARAMETERS',
     'DEFAULT_COLOUR_CYCLE',
     'DEFAULT_HATCH_PATTERNS',
+    'DEFAULT_PARAMETERS',
     'DEFAULT_PLOTTING_ILLUMINANT',
     'DEFAULT_PLOTTING_OECF',
     'ColourParameter',
diff --git a/colour/plotting/colorimetry.py b/colour/plotting/colorimetry.py
index 049493214..730c22501 100644
--- a/colour/plotting/colorimetry.py
+++ b/colour/plotting/colorimetry.py
@@ -256,8 +256,6 @@ def multi_spd_plot(spds,
         y_limit_min.append(min(values))
         y_limit_max.append(max(values))
 
-        matplotlib.pyplot.rc("axes", color_cycle=["r", "g", "b", "y"])
-
         if use_spds_colours:
             XYZ = spectral_to_XYZ(spd, cmfs, illuminant) / 100
             if normalise_spds_colours:
diff --git a/colour/plotting/common.py b/colour/plotting/common.py
index 1b0eae504..4dd9bfd4e 100644
--- a/colour/plotting/common.py
+++ b/colour/plotting/common.py
@@ -53,9 +53,9 @@ __all__ = ['PLOTTING_RESOURCES_DIRECTORY',
            'DEFAULT_FIGURE_HEIGHT',
            'DEFAULT_FIGURE_SIZE',
            'DEFAULT_FONT_SIZE',
-           'DEFAULT_PARAMETERS',
            'DEFAULT_COLOUR_CYCLE',
            'DEFAULT_HATCH_PATTERNS',
+           'DEFAULT_PARAMETERS',
            'DEFAULT_PLOTTING_ILLUMINANT',
            'DEFAULT_PLOTTING_OECF',
            'ColourParameter',
@@ -123,22 +123,6 @@ DEFAULT_FONT_SIZE : numeric
 if 'Qt4Agg' in matplotlib.get_backend():
     DEFAULT_FONT_SIZE = 10
 
-DEFAULT_PARAMETERS = {
-    'figure.figsize': DEFAULT_FIGURE_SIZE,
-    'font.size': DEFAULT_FONT_SIZE,
-    'axes.titlesize': DEFAULT_FONT_SIZE * 1.25,
-    'axes.labelsize': DEFAULT_FONT_SIZE * 1.25,
-    'legend.fontsize': DEFAULT_FONT_SIZE * 0.9,
-    'xtick.labelsize': DEFAULT_FONT_SIZE,
-    'ytick.labelsize': DEFAULT_FONT_SIZE}
-"""
-Default plotting parameters.
-
-DEFAULT_PARAMETERS : dict
-"""
-
-pylab.rcParams.update(DEFAULT_PARAMETERS)
-
 DEFAULT_COLOUR_CYCLE = ('r', 'g', 'b', 'c', 'm', 'y', 'k')
 """
 Default colour cycle for plots.
@@ -155,6 +139,23 @@ DEFAULT_HATCH_PATTERNS : tuple
 {'\\\\', 'o', 'x', '.', '*', '//'}
 """
 
+DEFAULT_PARAMETERS = {
+    'figure.figsize': DEFAULT_FIGURE_SIZE,
+    'font.size': DEFAULT_FONT_SIZE,
+    'axes.titlesize': DEFAULT_FONT_SIZE * 1.25,
+    'axes.labelsize': DEFAULT_FONT_SIZE * 1.25,
+    'legend.fontsize': DEFAULT_FONT_SIZE * 0.9,
+    'xtick.labelsize': DEFAULT_FONT_SIZE,
+    'ytick.labelsize': DEFAULT_FONT_SIZE,
+    'axes.color_cycle': DEFAULT_COLOUR_CYCLE}
+"""
+Default plotting parameters.
+
+DEFAULT_PARAMETERS : dict
+"""
+
+pylab.rcParams.update(DEFAULT_PARAMETERS)
+
 DEFAULT_PLOTTING_ILLUMINANT = ILLUMINANTS.get(
     'CIE 1931 2 Degree Standard Observer').get('D65')
 """


Test patch:
diff --git a/colour/algebra/tests/tests_extrapolation.py b/colour/algebra/tests/tests_extrapolation.py
index 9565bb613..00f7a4262 100644
--- a/colour/algebra/tests/tests_extrapolation.py
+++ b/colour/algebra/tests/tests_extrapolation.py
@@ -16,7 +16,7 @@ else:
     import unittest
 from itertools import permutations
 
-from colour.algebra import Extrapolator1d, LinearInterpolator1d
+from colour.algebra import Extrapolator1d, LinearInterpolator
 from colour.utilities import ignore_numpy_errors
 
 __author__ = 'Colour Developers'
@@ -62,14 +62,14 @@ class TestExtrapolator1d(unittest.TestCase):
         """
 
         extrapolator = Extrapolator1d(
-            LinearInterpolator1d(
+            LinearInterpolator(
                 np.array([5, 6, 7]),
                 np.array([5, 6, 7])))
         np.testing.assert_almost_equal(extrapolator((4, 8)), (4., 8.))
         self.assertEqual(extrapolator(4), 4.)
 
         extrapolator = Extrapolator1d(
-            LinearInterpolator1d(
+            LinearInterpolator(
                 np.array([3, 4, 5]),
                 np.array([1, 2, 3])),
             method='Constant')
@@ -78,7 +78,7 @@ class TestExtrapolator1d(unittest.TestCase):
         self.assertEqual(extrapolator(0.1), 1.)
 
         extrapolator = Extrapolator1d(
-            LinearInterpolator1d(
+            LinearInterpolator(
                 np.array([3, 4, 5]),
                 np.array([1, 2, 3])),
             method='Constant',
@@ -88,7 +88,7 @@ class TestExtrapolator1d(unittest.TestCase):
         self.assertEqual(extrapolator(0.1), 0.)
 
         extrapolator = Extrapolator1d(
-            LinearInterpolator1d(
+            LinearInterpolator(
                 np.array([3, 4, 5]),
                 np.array([1, 2, 3])),
             method='Constant',
@@ -108,7 +108,7 @@ class TestExtrapolator1d(unittest.TestCase):
         cases = set(permutations(cases * 3, r=3))
         for case in cases:
             extrapolator = Extrapolator1d(
-                LinearInterpolator1d(np.array(case), np.array(case)))
+                LinearInterpolator(np.array(case), np.array(case)))
             extrapolator(case[0])
 
 
diff --git a/colour/algebra/tests/tests_interpolation.py b/colour/algebra/tests/tests_interpolation.py
index ba5113322..0b85c20f8 100644
--- a/colour/algebra/tests/tests_interpolation.py
+++ b/colour/algebra/tests/tests_interpolation.py
@@ -17,7 +17,7 @@ else:
 from itertools import permutations
 
 from colour.algebra import (
-    LinearInterpolator1d,
+    LinearInterpolator,
     SpragueInterpolator)
 from colour.utilities import ignore_numpy_errors
 
@@ -362,7 +362,7 @@ SPRAGUE_INTERPOLATED_POINTS_DATA_A_10_SAMPLES = (
 class TestLinearInterpolator1d(unittest.TestCase):
     """
     Defines
-    :func:`colour.algebra.interpolation.LinearInterpolator1d` class units
+    :func:`colour.algebra.interpolation.LinearInterpolator` class units
     tests methods.
     """
 
@@ -375,7 +375,7 @@ class TestLinearInterpolator1d(unittest.TestCase):
                                'y')
 
         for attribute in required_attributes:
-            self.assertIn(attribute, dir(LinearInterpolator1d))
+            self.assertIn(attribute, dir(LinearInterpolator))
 
     def test_required_methods(self):
         """
@@ -385,18 +385,18 @@ class TestLinearInterpolator1d(unittest.TestCase):
         required_methods = ()
 
         for method in required_methods:
-            self.assertIn(method, dir(LinearInterpolator1d))
+            self.assertIn(method, dir(LinearInterpolator))
 
     def test___call__(self):
         """
         Tests
-        :func:`colour.algebra.interpolation.LinearInterpolator1d.__call__`
+        :func:`colour.algebra.interpolation.LinearInterpolator.__call__`
         method.
         """
 
         steps = 0.1
         x = np.arange(len(POINTS_DATA_A))
-        linear_interpolator = LinearInterpolator1d(x, POINTS_DATA_A)
+        linear_interpolator = LinearInterpolator(x, POINTS_DATA_A)
 
         for i, value in enumerate(
                 np.arange(0, len(POINTS_DATA_A) - 1 + steps, steps)):
@@ -414,7 +414,7 @@ class TestLinearInterpolator1d(unittest.TestCase):
     def test_nan__call__(self):
         """
         Tests
-        :func:`colour.algebra.interpolation.LinearInterpolator1d.__call__`
+        :func:`colour.algebra.interpolation.LinearInterpolator.__call__`
         method nan support.
         """
 
@@ -422,7 +422,7 @@ class TestLinearInterpolator1d(unittest.TestCase):
         cases = set(permutations(cases * 3, r=3))
         for case in cases:
             try:
-                linear_interpolator = LinearInterpolator1d(
+                linear_interpolator = LinearInterpolator(
                     np.array(case), np.array(case))
                 linear_interpolator(case[0])
             except ValueError:
diff --git a/colour/colorimetry/tests/tests_spectrum.py b/colour/colorimetry/tests/tests_spectrum.py
index f32ef44ce..5677a5226 100644
--- a/colour/colorimetry/tests/tests_spectrum.py
+++ b/colour/colorimetry/tests/tests_spectrum.py
@@ -2434,6 +2434,20 @@ class TestSpectralPowerDistribution(unittest.TestCase):
             rtol=0.0000001,
             atol=0.0000001)
 
+        np.testing.assert_almost_equal(
+            self.__spd.clone().interpolate(
+                SpectralShape(steps=1),
+                method='Linear')[410],
+            np.array(0.0643),
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            self.__spd.clone().interpolate(
+                SpectralShape(steps=1),
+                method='Pchip')[410],
+            np.array(0.06439937984496125),
+            decimal=7)
+
     def test_align(self):
         """
         Tests
@@ -2971,6 +2985,20 @@ class TestTriSpectralPowerDistribution(unittest.TestCase):
                 rtol=0.0000001,
                 atol=0.0000001)
 
+        np.testing.assert_almost_equal(
+            self.__tri_spd.clone().interpolate(
+                SpectralShape(steps=1),
+                method='Linear')[411],
+            np.array([0.050334, 0.001404, 0.24018]),
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            self.__tri_spd.clone().interpolate(
+                SpectralShape(steps=1),
+                method='Pchip')[411],
+            np.array([0.04895501, 0.00136229, 0.23349933]),
+            decimal=7)
+
     def test_align(self):
         """
         Tests
diff --git a/colour/colorimetry/tests/tests_tristimulus.py b/colour/colorimetry/tests/tests_tristimulus.py
index 33fecf517..97a72d384 100644
--- a/colour/colorimetry/tests/tests_tristimulus.py
+++ b/colour/colorimetry/tests/tests_tristimulus.py
@@ -208,6 +208,30 @@ class TestWavelength_to_XYZ(unittest.TestCase):
             np.array([0.44575583, 0.18184213, 0.]),
             decimal=7)
 
+        np.testing.assert_almost_equal(
+            wavelength_to_XYZ(
+                480.5,
+                CMFS.get('CIE 2012 2 Degree Standard Observer'),
+                'Cubic Spline'),
+            np.array([0.07773422, 0.18148028, 0.7337162]),
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            wavelength_to_XYZ(
+                480.5,
+                CMFS.get('CIE 2012 2 Degree Standard Observer'),
+                'Linear'),
+            np.array([0.07779856, 0.18149335, 0.7340129]),
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            wavelength_to_XYZ(
+                480.5,
+                CMFS.get('CIE 2012 2 Degree Standard Observer'),
+                'Pchip'),
+            np.array([0.07773515, 0.18148048, 0.73372294]),
+            decimal=7)
+
     def test_n_dimensional_wavelength_to_XYZ(self):
         """
         Tests