Repository: colour-science/colour

Issue:
Implement support for "Barten (1999)" contrast sensitivity function.
## References
- http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=837000
- http://car.france3.mars.free.fr/Formation%20INA%20HD/HDTV/HDTV%20%202007%20v35/SMPTE%20normes%20et%20confs/Contrastm.pdf
- http://www.umich.edu/~ners580/ners-bioe_481/lectures/pdfs/Barten_ch3s1-9.pdf
- https://www.smpte.org/sites/default/files/2014-05-06-EOTF-Miller-1-2-handout.pdf
- https://www.smpte.org/sites/default/files/23-1615-TS7-2-IProc02-Miller.pdf
- http://www.avsforum.com/forum/139-display-calibration/2020689-need-help-barten-contrast-sensitivity-function.html

Discussion:
coveralls: 
[![Coverage Status](https://coveralls.io/builds/21420529/badge)](https://coveralls.io/builds/21420529)

Coverage increased (+0.03%) to 97.477% when pulling **954ef16c64663a7695c0616ceed6a1ac529f5d26 on feature/csf_barten1999** into **3f05c910ba35923e9bdcd7fc63648b58653d0386 on develop**.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/21420529/badge)](https://coveralls.io/builds/21420529)

Coverage increased (+0.03%) to 97.477% when pulling **954ef16c64663a7695c0616ceed6a1ac529f5d26 on feature/csf_barten1999** into **3f05c910ba35923e9bdcd7fc63648b58653d0386 on develop**.

KelSolaar: Here is a test with a maximised spatial frequency *u* versus fixed *u = 15*:

![image](https://user-images.githubusercontent.com/99779/52169888-5ffc4980-27a4-11e9-9922-d1fc133deadc.png)

It is actually going way below 2% for low luminance values, fixed *u* was closer to *Figure 31*. I'm starting to wonder what settings were used to generate the figure.

Here is the code used to generate mine:

```python
# %%
from __future__ import division

import numpy as np
from scipy.optimize import fmin

from colour.plotting import colour_style, plot_multi_functions
from colour.contrast.barten1999 import (
    function_contrast_sensitivity_Barten1999, pupil_diameter_Barten1999,
    retinal_illuminance_Barten1999, sigma_Barten1999)

_ = colour_style()

# %%
settings_BT2246 = {
    'k': 3.0,
    'T': 0.1,
    'X_max': 12,
    'N_max': 15,
    'n': 0.03,
    'p': 1.2274 * 10**6,
    'phi_0': 3 * 10**-8,
    'u_0': 7,
}


def maximise_spatial_frequency(L):
    maximised_spatial_frequency = []
    for L_v in L:
        X_0 = 60
        d = pupil_diameter_Barten1999(L_v, X_0)
        sigma = sigma_Barten1999(0.5 / 60, 0.08 / 60, d)
        E = retinal_illuminance_Barten1999(L_v, d, True)
        CSF = lambda x: -function_contrast_sensitivity_Barten1999(
            u=x,
            sigma=sigma,
            X_0=X_0,
            E=E,
            **settings_BT2246)

        maximised_spatial_frequency.append(fmin(CSF, 0, disp=False)[0])

    return maximised_spatial_frequency


L = np.logspace(np.log10(0.01), np.log10(100), 100)

X_0 = 60
d = pupil_diameter_Barten1999(L, X_0)
sigma = sigma_Barten1999(0.5 / 60, 0.08 / 60, d)
E = retinal_illuminance_Barten1999(L, d, True)
u = maximise_spatial_frequency(L)

y1 = lambda x: 1 / function_contrast_sensitivity_Barten1999(
    u=u,
    sigma=sigma,
    X_0=X_0,
    E=E,
    **settings_BT2246) * 2 * (1/ 1.27)

y2 = lambda x: 1 / function_contrast_sensitivity_Barten1999(
    u=X_0 / 15,
    sigma=sigma,
    X_0=X_0,
    E=E,
    **settings_BT2246) * 2 * (1/ 1.27)


plot_multi_functions(
    {'Maximised U': y1, 'Fixed U': y2},
    samples=L,
    log_x=10,
    # log_y=10,
    bounding_box=[0.1, 100, 0.001, 0.025],
    **{'axes.grid.which': 'both'})

```
KelSolaar: It works! I was actually clipping the X-axis:

![image](https://user-images.githubusercontent.com/99779/52170276-6262a180-27ac-11e9-9b74-8bc5b5fcc9ea.png)

KelSolaar: @nick-shaw, @MichaelMauderer Any chances to look at that? Pretty sure you started @nick-shaw, should we merge?
KelSolaar: Fixed the _BIBLIOGRAPHY.bib_ conflict by rebasing on top of `develop`, it is not super elegant but at this stage and given this file acts almost like a binary file where changes can happen anywhere without much control, it is the path of least resistance.

Patch:
diff --git a/colour/__init__.py b/colour/__init__.py
index 2b84d6abb..63428af59 100644
--- a/colour/__init__.py
+++ b/colour/__init__.py
@@ -17,6 +17,7 @@ Sub-packages
 -   colorimetry: Core objects for colour computations.
 -   constants: *CIE* and *CODATA* constants.
 -   continuous: Base objects for continuous data representation.
+-   contrast: Objects for contrast sensitivity computation.
 -   corresponding: Corresponding colour chromaticities computations.
 -   difference: Colour difference computations.
 -   examples: Examples for the sub-packages.
@@ -126,6 +127,8 @@ from .corresponding import (BRENEMAN_EXPERIMENTS,
                             BRENEMAN_EXPERIMENTS_PRIMARIES_CHROMATICITIES,
                             CORRESPONDING_CHROMATICITIES_PREDICTION_MODELS,
                             corresponding_chromaticities_prediction)
+from .contrast import (CONTRAST_SENSITIVITY_METHODS,
+                       contrast_sensitivity_function)
 from .phenomena import (rayleigh_scattering, scattering_cross_section,
                         sd_rayleigh_scattering)
 from .notation import (MUNSELL_COLOURS, MUNSELL_VALUE_METHODS,
@@ -251,6 +254,7 @@ __all__ += [
     'CORRESPONDING_CHROMATICITIES_PREDICTION_MODELS',
     'corresponding_chromaticities_prediction'
 ]
+__all__ += ['CONTRAST_SENSITIVITY_METHODS', 'contrast_sensitivity_function']
 __all__ += [
     'rayleigh_scattering', 'scattering_cross_section', 'sd_rayleigh_scattering'
 ]
diff --git a/colour/contrast/__init__.py b/colour/contrast/__init__.py
new file mode 100644
index 000000000..54bc365b2
--- /dev/null
+++ b/colour/contrast/__init__.py
@@ -0,0 +1,141 @@
+# -*- coding: utf-8 -*-
+"""
+References
+----------
+-   :cite:`Barten1999` : Barten, P. G. (1999). Contrast Sensitivity of the
+    Human Eye and Its Effects on Image Quality. SPIE. doi:10.1117/3.353254
+-   :cite:`Barten2003` : Barten, P. G. J. (2003). Formula for the contrast
+    sensitivity of the human eye. In Y. Miyake & D. R. Rasmussen (Eds.),
+    Proceedings of SPIE (Vol. 5294, pp. 231-238). doi:10.1117/12.537476
+-   :cite:`Cowan2004` : Cowan, M., Kennel, G., Maier, T., Walker, B., & Cowan,
+    M. (2004). Constant Sensitivity Experiment to Determine the Bit Depth for
+    Digital Cinema, (September). Retrieved from
+    http://car.france3.mars.free.fr/Formation INA HD/HDTV/HDTV  2007 v35/\
+SMPTE normes et confs/Contrastm.pdf
+-   :cite:`InternationalTelecommunicationUnion2015` : International
+    Telecommunication Union. (2015). Report ITU-R BT.2246-4 - The present state
+    of ultra-high definition television BT Series Broadcasting service, 5,
+    1-92.
+"""
+
+from __future__ import absolute_import
+
+from colour.utilities import CaseInsensitiveMapping, filter_kwargs
+from .barten1999 import (optical_MTF_Barten1999, pupil_diameter_Barten1999,
+                         sigma_Barten1999, retinal_illuminance_Barten1999,
+                         maximum_angular_size_Barten1999,
+                         contrast_sensitivity_function_Barten1999)
+
+__all__ = [
+    'optical_MTF_Barten1999', 'pupil_diameter_Barten1999', 'sigma_Barten1999',
+    'retinal_illuminance_Barten1999', 'maximum_angular_size_Barten1999',
+    'contrast_sensitivity_function_Barten1999'
+]
+
+CONTRAST_SENSITIVITY_METHODS = CaseInsensitiveMapping({
+    'Barten 1999': contrast_sensitivity_function_Barten1999,
+})
+CONTRAST_SENSITIVITY_METHODS.__doc__ = """
+Supported contrast sensitivity methods.
+
+References
+----------
+:cite:`Barten1999`, :cite:`Barten2003`, :cite:`Cowan2004`,
+:cite:`InternationalTelecommunicationUnion2015`,
+
+CONTRAST_SENSITIVITY_METHODS : CaseInsensitiveMapping
+    **{'Barten 1999'}**
+"""
+
+
+def contrast_sensitivity_function(method='Barten 1999', **kwargs):
+    """
+    Returns the contrast sensitivity :math:`S` of the human eye according to
+    the contrast sensitivity function (CSF) described by given method.
+
+    Parameters
+    ----------
+    method : unicode, optional
+        **{'Barten 1999'}**,
+        Computation method.
+
+    Other Parameters
+    ----------------
+    E : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Retinal illuminance :math:`E` in Trolands.
+    N_max : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Maximum number of cycles :math:`N_{max}` over which the eye can
+        integrate the information.
+    T : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Integration time :math:`T` in seconds of the eye.
+    X_0 : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Angular size :math:`X_0` in degrees of the object in the x direction.
+    Y_0 : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Angular size :math:`Y_0` in degrees of the object in the y direction.
+    X_max : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Maximum angular size :math:`X_{max}` in degrees of the integration
+        area in the x direction.
+    Y_max : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Maximum angular size :math:`Y_{max}` in degrees of the integration
+        area in the y direction.
+    k : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Signal-to-noise (SNR) ratio :math:`k`.
+    n : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Quantum efficiency of the eye :math:`n`.
+    p : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Photon conversion factor :math:`p` in
+        :math:`photons\\div seconds\\div degrees^2\\div Trolands` that
+        depends on the light source.
+    phi_0 : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Spectral density :math:`\\phi_0` in :math:`seconds degrees^2` of the
+        neural noise.
+    sigma : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Standard deviation :math:`\\sigma` of the line-spread function
+        resulting from the convolution of the different elements of the
+        convolution process.
+    u : numeric
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Spatial frequency :math:`u`, the cycles per degree.
+    u_0 : numeric or array_like, optional
+        {:func:`colour.contrast.contrast_sensitivity_function_Barten1999`},
+        Spatial frequency :math:`u_0` in :math:`cycles\\div degrees` above
+        which the lateral inhibition ceases.
+
+    Returns
+    -------
+    ndarray
+        Contrast sensitivity :math:`S`.
+
+    References
+    ----------
+    :cite:`Barten1999`, :cite:`Barten2003`, :cite:`Cowan2004`,
+    :cite:`InternationalTelecommunicationUnion2015`,
+
+    Examples
+    --------
+    >>> contrast_sensitivity_function(u=4)  # doctest: +ELLIPSIS
+    360.8691122...
+    >>> contrast_sensitivity_function('Barten 1999', u=4)  # doctest: +ELLIPSIS
+    360.8691122...
+    """
+
+    function = CONTRAST_SENSITIVITY_METHODS[method]
+
+    S = function(**filter_kwargs(function, **kwargs))
+
+    return S
+
+
+__all__ += ['CONTRAST_SENSITIVITY_METHODS', 'contrast_sensitivity_function']
diff --git a/colour/contrast/barten1999.py b/colour/contrast/barten1999.py
new file mode 100644
index 000000000..39eae6831
--- /dev/null
+++ b/colour/contrast/barten1999.py
@@ -0,0 +1,462 @@
+# -*- coding: utf-8 -*-
+"""
+Barten (1999) Contrast Sensitivity Function
+===========================================
+
+Defines *Barten (1999)* contrast sensitivity function:
+
+-   :func:`colour.contrast.contrast_sensitivity_function_Barten1999`
+
+See Also
+--------
+`Barten (1999) Contrast Sensitivity Function Jupyter Notebook
+<http://nbviewer.jupyter.org/github/colour-science/colour-notebooks/\
+blob/master/notebooks/contrast/barten1999.ipynb>`_
+
+References
+----------
+-   :cite:`Barten1999` : Barten, P. G. (1999). Contrast Sensitivity of the
+    Human Eye and Its Effects on Image Quality. SPIE. doi:10.1117/3.353254
+-   :cite:`Barten2003` : Barten, P. G. J. (2003). Formula for the contrast
+    sensitivity of the human eye. In Y. Miyake & D. R. Rasmussen (Eds.),
+    Proceedings of SPIE (Vol. 5294, pp. 231-238). doi:10.1117/12.537476
+-   :cite:`Cowan2004` : Cowan, M., Kennel, G., Maier, T., Walker, B., & Cowan,
+    M. (2004). Constant Sensitivity Experiment to Determine the Bit Depth for
+    Digital Cinema, (September). Retrieved from
+    http://car.france3.mars.free.fr/Formation INA HD/HDTV/HDTV  2007 v35/\
+SMPTE normes et confs/Contrastm.pdf
+-   :cite:`InternationalTelecommunicationUnion2015` : International
+    Telecommunication Union. (2015). Report ITU-R BT.2246-4 - The present state
+    of ultra-high definition television BT Series Broadcasting service, 5,
+    1-92.
+"""
+
+from __future__ import division, unicode_literals
+
+import numpy as np
+
+from colour.utilities import as_float_array, as_float
+
+__author__ = 'Colour Developers'
+__copyright__ = 'Copyright (C) 2013-2019 - Colour Developers'
+__license__ = 'New BSD License - http://opensource.org/licenses/BSD-3-Clause'
+__maintainer__ = 'Colour Developers'
+__email__ = 'colour-science@googlegroups.com'
+__status__ = 'Production'
+
+__all__ = [
+    'optical_MTF_Barten1999', 'pupil_diameter_Barten1999', 'sigma_Barten1999',
+    'retinal_illuminance_Barten1999', 'maximum_angular_size_Barten1999',
+    'contrast_sensitivity_function_Barten1999'
+]
+
+
+def optical_MTF_Barten1999(u, sigma=0.01):
+    """
+    Returns the optical modulation transfer function (MTF) :math:`M_{opt}` of
+    the eye using *Barten (1999)* method.
+
+    Parameters
+    ----------
+    u : numeric or array_like
+        Spatial frequency :math:`u`, the cycles per degree.
+    sigma : numeric or array_like, optional
+        Standard deviation :math:`\\sigma` of the line-spread function
+        resulting from the convolution of the different elements of the
+        convolution process.
+
+    Returns
+    -------
+    numeric or array_like
+        Optical modulation transfer function (MTF) :math:`M_{opt}` of the eye.
+
+    References
+    ----------
+    :cite:`Barten1999`, :cite:`Barten2003`, :cite:`Cowan2004`,
+    :cite:`InternationalTelecommunicationUnion2015`,
+
+    Examples
+    --------
+    >>> optical_MTF_Barten1999(4, 0.01)  # doctest: +ELLIPSIS
+    0.9689107...
+    """
+
+    u = as_float_array(u)
+    sigma = as_float_array(sigma)
+
+    return as_float(np.exp(-2 * np.pi ** 2 * sigma ** 2 * u ** 2))
+
+
+def pupil_diameter_Barten1999(L, X_0=60, Y_0=None):
+    """
+    Returns the pupil diameter for given luminance and object or stimulus
+    angular size using *Barten (1999)* method.
+
+    Parameters
+    ----------
+    L : numeric or array_like
+        Average luminance :math:`L` in :math:`cd/m^2`.
+    X_0 : numeric or array_like, optional
+        Angular size of the object :math:`X_0` in degrees in the x direction.
+    Y_0 : numeric or array_like, optional
+        Angular size of the object :math:`X_0` in degrees in the y direction.
+
+    References
+    ----------
+    :cite:`Barten1999`, :cite:`Barten2003`, :cite:`Cowan2004`,
+    :cite:`InternationalTelecommunicationUnion2015`,
+
+    Examples
+    --------
+    >>> pupil_diameter_Barten1999(100, 60, 60)  # doctest: +ELLIPSIS
+    2.0777571...
+    """
+
+    L = as_float_array(L)
+    X_0 = as_float_array(X_0)
+    Y_0 = X_0 if Y_0 is None else as_float_array(Y_0)
+
+    return as_float(5 - 3 * np.tanh(0.4 * np.log(L * X_0 * Y_0 / 40 ** 2)))
+
+
+def sigma_Barten1999(sigma_0=0.5 / 60, C_ab=0.08 / 60, d=2.1):
+    """
+    Returns the standard deviation :math:`\\sigma` of the line-spread function
+    resulting from the convolution of the different elements of the convolution
+    process using *Barten (1999)* method.
+
+    The :math:`\\sigma` quantity depends on the pupil diameter :math:`d` of the
+    eye lens. For very small pupil diameters, :math:`\\sigma` increases
+    inversely proportionally with pupil size because of diffraction, and for
+    large pupil diameters, :math:`\\sigma` increases about linearly with pupil
+    size because of chromatic aberration and others aberrations.
+
+    Parameters
+    ----------
+    sigma_0 : numeric or array_like, optional
+        Constant :math:`\\sigma_{0}` in degrees.
+    C_ab : numeric or array_like, optional
+        Spherical aberration of the eye :math:`C_{ab}` in
+        :math:`degrees\\div mm`.
+    d : numeric or array_like, optional
+        Pupil diameter :math:`d` in millimeters.
+
+    Returns
+    -------
+    ndarray
+        Standard deviation :math:`\\sigma` of the line-spread function
+        resulting from the convolution of the different elements of the
+        convolution process.
+
+    Warnings
+    --------
+    This definition expects :math:`\\sigma_{0}` and :math:`C_{ab}` to be given
+    in degrees and :math:`degrees\\div mm` respectively. However, in the
+    literature, the values for :math:`\\sigma_{0}` and
+    :math:`C_{ab}` are usually given in :math:`arc min` and
+    :math:`arc min\\div mm` respectively, thus they need to be divided by 60.
+
+    References
+    ----------
+    :cite:`Barten1999`, :cite:`Barten2003`, :cite:`Cowan2004`,
+    :cite:`InternationalTelecommunicationUnion2015`,
+
+    Examples
+    --------
+    >>> sigma_Barten1999(0.5 / 60, 0.08 / 60, 2.1)  # doctest: +ELLIPSIS
+    0.0087911...
+    """
+
+    sigma_0 = as_float_array(sigma_0)
+    C_ab = as_float_array(C_ab)
+    d = as_float_array(d)
+
+    return as_float(np.sqrt(sigma_0 ** 2 + (C_ab * d) ** 2))
+
+
+def retinal_illuminance_Barten1999(
+        L, d=2.1, apply_stiles_crawford_effect_correction=True):
+    """
+    Returns the retinal illuminance :math:`E` in Trolands for given average
+    luminance :math:`L` and pupil diameter :math:`d` using *Barten (1999)*
+    method.
+
+    Parameters
+    ----------
+    L : numeric or array_like
+        Average luminance :math:`L` in :math:`cd/m^2`.
+    d : numeric or array_like, optional
+        Pupil diameter :math:`d` in millimeters.
+    apply_stiles_crawford_effect_correction : bool, optional
+        Whether to apply the correction for *Stiles-Crawford* effect.
+
+    Returns
+    -------
+    ndarray
+        Retinal illuminance :math:`E` in Trolands.
+
+    Notes
+    -----
+    -   This definition is for use with photopic viewing conditions and thus
+        corrects for the Stiles-Crawford effect by default, i.e. directional
+        sensitivity of the cone cells with lower response of cone cells
+        receiving light from the edge of the pupil.
+
+    References
+    ----------
+    :cite:`Barten1999`, :cite:`Barten2003`, :cite:`Cowan2004`,
+    :cite:`InternationalTelecommunicationUnion2015`,
+
+    Examples
+    --------
+    >>> retinal_illuminance_Barten1999(100, 2.1)  # doctest: +ELLIPSIS
+    330.4115803...
+    >>> retinal_illuminance_Barten1999(100, 2.1, False)  # doctest: +ELLIPSIS
+    346.3605900...
+    """
+
+    d = as_float_array(d)
+    L = as_float_array(L)
+
+    E = (np.pi * d ** 2) / 4 * L
+
+    if apply_stiles_crawford_effect_correction:
+        E *= (1 - (d / 9.7) ** 2 + (d / 12.4) ** 4)
+
+    return E
+
+
+def maximum_angular_size_Barten1999(u, X_0=60, X_max=12, N_max=15):
+    """
+    Returns the maximum angular size :math:`X` of the object considered using
+    *Barten (1999)* method.
+
+    Parameters
+    ----------
+    u : numeric
+        Spatial frequency :math:`u`, the cycles per degree.
+    X_0 : numeric or array_like, optional
+        Angular size :math:`X_0` in degrees of the object in the x direction.
+    X_max : numeric or array_like, optional
+        Maximum angular size :math:`X_{max}` in degrees of the integration
+        area in the x direction.
+    N_max : numeric or array_like, optional
+        Maximum number of cycles :math:`N_{max}` over which the eye can
+        integrate the information.
+
+    Returns
+    -------
+    numeric or ndarray
+        Maximum angular size :math:`X` of the object considered.
+
+    References
+    ----------
+    :cite:`Barten1999`, :cite:`Barten2003`, :cite:`Cowan2004`,
+    :cite:`InternationalTelecommunicationUnion2015`,
+
+    Examples
+    --------
+    >>> maximum_angular_size_Barten1999(4)  # doctest: +ELLIPSIS
+    3.5729480...
+    """
+
+    u = as_float_array(u)
+    X_0 = as_float_array(X_0)
+    X_max = as_float_array(X_max)
+    N_max = as_float_array(N_max)
+
+    return (1 / X_0 ** 2 + 1 / X_max ** 2 + u ** 2 / N_max ** 2) ** -0.5
+
+
+def contrast_sensitivity_function_Barten1999(u,
+                                             sigma=sigma_Barten1999(
+                                                 0.5 / 60, 0.08 / 60, 2.1),
+                                             k=3.0,
+                                             T=0.1,
+                                             X_0=60,
+                                             Y_0=None,
+                                             X_max=12,
+                                             Y_max=None,
+                                             N_max=15,
+                                             n=0.03,
+                                             p=1.2274 * 10 ** 6,
+                                             E=retinal_illuminance_Barten1999(
+                                                 20, 2.1),
+                                             phi_0=3 * 10 ** -8,
+                                             u_0=7):
+    """
+    Returns the contrast sensitivity :math:`S` of the human eye according to
+    the contrast sensitivity function (CSF) described by *Barten (1999)*.
+
+    Contrast sensitivity is defined as the inverse of the modulation threshold
+    of a sinusoidal luminance pattern. The modulation threshold of this pattern
+    is generally defined by 50% probability of detection. The contrast
+    sensitivity function or CSF gives the contrast sensitivity as a function of
+    spatial frequency. In the CSF, the spatial frequency is expressed in
+    angular units with respect to the eye. It reaches a maximum between 1 and
+    10 cycles per degree with a fall off at higher and lower spatial
+    frequencies.
+
+    Parameters
+    ----------
+    u : numeric
+        Spatial frequency :math:`u`, the cycles per degree.
+    sigma : numeric or array_like, optional
+        Standard deviation :math:`\\sigma` of the line-spread function
+        resulting from the convolution of the different elements of the
+        convolution process.
+    k : numeric or array_like, optional
+        Signal-to-noise (SNR) ratio :math:`k`.
+    T : numeric or array_like, optional
+        Integration time :math:`T` in seconds of the eye.
+    X_0 : numeric or array_like, optional
+        Angular size :math:`X_0` in degrees of the object in the x direction.
+    Y_0 : numeric or array_like, optional
+        Angular size :math:`Y_0` in degrees of the object in the y direction.
+    X_max : numeric or array_like, optional
+        Maximum angular size :math:`X_{max}` in degrees of the integration
+        area in the x direction.
+    Y_max : numeric or array_like, optional
+        Maximum angular size :math:`Y_{max}` in degrees of the integration
+        area in the y direction.
+    N_max : numeric or array_like, optional
+        Maximum number of cycles :math:`N_{max}` over which the eye can
+        integrate the information.
+    n : numeric or array_like, optional
+        Quantum efficiency of the eye :math:`n`.
+    p : numeric or array_like, optional
+        Photon conversion factor :math:`p` in
+        :math:`photons\\div seconds\\div degrees^2\\div Trolands` that
+        depends on the light source.
+    E : numeric or array_like, optional
+        Retinal illuminance :math:`E` in Trolands.
+    phi_0 : numeric or array_like, optional
+        Spectral density :math:`\\phi_0` in :math:`seconds degrees^2` of the
+        neural noise.
+    u_0 : numeric or array_like, optional
+        Spatial frequency :math:`u_0` in :math:`cycles\\div degrees` above
+        which the lateral inhibition ceases.
+
+    Returns
+    -------
+    ndarray
+        Contrast sensitivity :math:`S`.
+
+    Warnings
+    --------
+    This definition expects :math:`\\sigma_{0}` and :math:`C_{ab}` used in the
+    computation of :math:`\\sigma` to be given in degrees and
+    :math:`degrees\\div mm` respectively. However, in the literature, the
+    values for :math:`\\sigma_{0}` and :math:`C_{ab}` are usually given in
+    :math:`arc min` and :math:`arc min\\div mm` respectively, thus they need to
+    be divided by 60.
+
+    Notes
+    -----
+    -   The formula holds for bilateral viewing and for equal dimensions of
+        the object in x and y direction. For monocular vision, the contrast
+        sensitivity is a factor :math:`\\sqrt{2}` smaller.
+    -   *Barten (1999)* CSF default values for the :math:`k`,
+        :math:`\\sigma_{0}`, :math:`C_{ab}`, :math:`T`, :math:`X_{max}`,
+        :math:`N_{max}`, :math:`n`, :math:`\\phi_{0}` and :math:`u_0` constants
+        are valid for a standard observer with good vision and with an age
+        between 20 and 30 years.
+    -   The other constants have been filled using reference data from
+        *Figure 31* in :cite:`InternationalTelecommunicationUnion2015` but
+        must be adapted to the current use case.
+    -   The product of :math:`u`, the cycles per degree, and :math:`X_0`,
+        the number of degrees, gives the number of cycles :math:`P_c` in a
+        pattern. Therefore, :math:`X_0` can be made a variable dependent on
+        :math:`u` such as :math:`X_0 = P_c / u`.
+
+    References
+    ----------
+    :cite:`Barten1999`, :cite:`Barten2003`, :cite:`Cowan2004`,
+    :cite:`InternationalTelecommunicationUnion2015`,
+
+    Examples
+    --------
+    >>> contrast_sensitivity_function_Barten1999(4)  # doctest: +ELLIPSIS
+    360.8691122...
+
+    Reproducing *Figure 31* in \
+:cite:`InternationalTelecommunicationUnion2015` illustrating the minimum
+    detectable contrast according to *Barten (1999)* model with the assumed
+    conditions for UHDTV applications. The minimum detectable contrast
+    :math:`MDC` is then defined as follows::
+
+        :math:`MDC = 1 / CSF * 2 * (1 / 1.27)`
+
+    where :math:`2` is used for the conversion from modulation to contrast and
+    :math:`1 / 1.27` is used for the conversion from sinusoidal to rectangular
+    waves.
+
+    >>> from scipy.optimize import fmin
+    >>> settings_BT2246 = {
+    ...     'k': 3.0,
+    ...     'T': 0.1,
+    ...     'X_max': 12,
+    ...     'N_max': 15,
+    ...     'n': 0.03,
+    ...     'p': 1.2274 * 10 ** 6,
+    ...     'phi_0': 3 * 10 ** -8,
+    ...     'u_0': 7,
+    ... }
+    >>>
+    >>> def maximise_spatial_frequency(L):
+    ...     maximised_spatial_frequency = []
+    ...     for L_v in L:
+    ...         X_0 = 60
+    ...         d = pupil_diameter_Barten1999(L_v, X_0)
+    ...         sigma = sigma_Barten1999(0.5 / 60, 0.08 / 60, d)
+    ...         E = retinal_illuminance_Barten1999(L_v, d, True)
+    ...         maximised_spatial_frequency.append(
+    ...             fmin(lambda x: (
+    ...                     -contrast_sensitivity_function_Barten1999(
+    ...                         u=x,
+    ...                         sigma=sigma,
+    ...                         X_0=X_0,
+    ...                         E=E,
+    ...                         **settings_BT2246)
+    ...                 ), 0, disp=False)[0])
+    ...         return as_float(np.array(maximised_spatial_frequency))
+    >>>
+    >>> L = np.logspace(np.log10(0.01), np.log10(100), 10)
+    >>> X_0 = Y_0 = 60
+    >>> d = pupil_diameter_Barten1999(L, X_0, Y_0)
+    >>> sigma = sigma_Barten1999(0.5 / 60, 0.08 / 60, d)
+    >>> E = retinal_illuminance_Barten1999(L, d)
+    >>> u = maximise_spatial_frequency(L)
+    >>> (1 / contrast_sensitivity_function_Barten1999(
+    ...     u=u, sigma=sigma, E=E, X_0=X_0, Y_0=Y_0, **settings_BT2246)
+    ...  * 2 * (1/ 1.27))
+    ... # doctest: +ELLIPSIS
+    array([ 0.0207396...,  0.0134885...,  0.0096063...,  0.0077299...,  \
+0.0068983...,
+            0.0065057...,  0.0062712...,  0.0061198...,  0.0060365...,  \
+0.0059984...])
+    """
+
+    u = as_float_array(u)
+    k = as_float_array(k)
+    T = as_float_array(T)
+    X_0 = as_float_array(X_0)
+    Y_0 = X_0 if Y_0 is None else as_float_array(Y_0)
+    X_max = as_float_array(X_max)
+    Y_max = X_max if Y_max is None else as_float_array(Y_max)
+    N_max = as_float_array(N_max)
+    n = as_float_array(n)
+    p = as_float_array(p)
+    E = as_float_array(E)
+    phi_0 = as_float_array(phi_0)
+    u_0 = as_float_array(u_0)
+
+    M_opt = optical_MTF_Barten1999(u, sigma)
+
+    M_as = 1 / (maximum_angular_size_Barten1999(u, X_0, X_max, N_max) *
+                maximum_angular_size_Barten1999(u, Y_0, Y_max, N_max))
+
+    S = (M_opt / k) / np.sqrt(2 / T * M_as * (1 / (n * p * E) + phi_0 /
+                                              (1 - np.exp(-(u / u_0) ** 2))))
+
+    return as_float(S)
diff --git a/colour/examples/contrast/examples_contrast.py b/colour/examples/contrast/examples_contrast.py
new file mode 100644
index 000000000..02f0c8a95
--- /dev/null
+++ b/colour/examples/contrast/examples_contrast.py
@@ -0,0 +1,109 @@
+# -*- coding: utf-8 -*-
+"""
+Showcases contrast sensitivity computations.
+"""
+
+from pprint import pprint
+
+import numpy as np
+from scipy.optimize import fmin
+
+import colour
+from colour.utilities import as_float, message_box
+from colour.plotting import colour_style, plot_single_function
+
+message_box('Contrast Sensitivity Computations')
+
+colour_style()
+
+message_box(('Computing the contrast sensitivity for a spatial frequency "u" '
+             'of 4, an angular size "X_0" of 60 and a retinal illuminance "E" '
+             'of 65 using "Barten (1999)" method.'))
+pprint(colour.contrast_sensitivity_function(u=4, X_0=60, E=65))
+pprint(
+    colour.contrast.contrast_sensitivity_function_Barten1999(
+        u=4, X_0=60, E=65))
+
+print('\n')
+
+message_box(('Computing the minimum detectable contrast with the assumed '
+             'conditions for UHDTV applications as given in "ITU-R BT.2246-4"'
+             '"Figure 31" and using "Barten (1999)" method.'))
+
+settings_BT2246 = {
+    'k': 3.0,
+    'T': 0.1,
+    'X_max': 12,
+    'N_max': 15,
+    'n': 0.03,
+    'p': 1.2274 * 10 ** 6,
+    'phi_0': 3 * 10 ** -8,
+    'u_0': 7,
+}
+
+
+def maximise_spatial_frequency(L):
+    """
+    Maximises the spatial frequency :math:`u` for given luminance value.
+
+    Parameters
+    ----------
+    L : numeric or array_like
+        Luminance value at which to maximize the spatial frequency :math:`u`.
+
+    Returns
+    -------
+    numeric or ndarray
+        Maximised spatial frequency :math:`u`.
+    """
+
+    maximised_spatial_frequency = []
+    for L_v in L:
+        X_0 = 60
+        d = colour.contrast.pupil_diameter_Barten1999(L_v, X_0)
+        sigma = colour.contrast.sigma_Barten1999(0.5 / 60, 0.08 / 60, d)
+        E = colour.contrast.retinal_illuminance_Barten1999(L_v, d, True)
+        maximised_spatial_frequency.append(
+            fmin(lambda x: (
+                    -colour.contrast.contrast_sensitivity_function_Barten1999(
+                        u=x,
+                        sigma=sigma,
+                        X_0=X_0,
+                        E=E,
+                        **settings_BT2246)
+                ), 0, disp=False)[0])
+
+    return as_float(np.array(maximised_spatial_frequency))
+
+
+L = np.logspace(np.log10(0.01), np.log10(100), 100)
+X_0 = Y_0 = 60
+d = colour.contrast.barten1999.pupil_diameter_Barten1999(L, X_0, Y_0)
+sigma = colour.contrast.barten1999.sigma_Barten1999(0.5 / 60, 0.08 / 60, d)
+E = colour.contrast.barten1999.retinal_illuminance_Barten1999(L, d)
+u = maximise_spatial_frequency(L)
+
+pprint(1 / colour.contrast_sensitivity_function(
+    u=u, sigma=sigma, E=E, X_0=X_0, Y_0=Y_0, **settings_BT2246) * 2 *
+       (1 / 1.27))
+pprint(1 / colour.contrast.contrast_sensitivity_function_Barten1999(
+    u=u, sigma=sigma, E=E, X_0=X_0, Y_0=Y_0, **settings_BT2246) * 2 *
+       (1 / 1.27))
+
+plot_single_function(
+    lambda x: (
+        1 / colour.contrast.contrast_sensitivity_function_Barten1999(
+            u=u, sigma=sigma, E=E, X_0=X_0, Y_0=Y_0, **settings_BT2246)
+        * 2 * (1 / 1.27)),
+    samples=L,
+    log_x=10,
+    **{
+        'title':
+            'Examples of HVS Minimum Detectable Contrast Characteristics',
+        'x_label':
+            'Luminance ($cd/m^2$)',
+        'y_label':
+            'Minimum Detectable Contrast',
+        'axes.grid.which':
+            'both'
+    })
diff --git a/docs/colour.contrast.rst b/docs/colour.contrast.rst
new file mode 100644
index 000000000..617f56406
--- /dev/null
+++ b/docs/colour.contrast.rst
@@ -0,0 +1,44 @@
+Contrast Sensitivity
+====================
+
+.. contents:: :local:
+
+Contrast Sensitivity
+--------------------
+
+``colour``
+
+.. currentmodule:: colour
+
+.. autosummary::
+    :toctree: generated/
+
+    contrast_sensitivity_function
+    CONTRAST_SENSITIVITY_METHODS
+
+Barten (1999) Contrast Sensitivity Function
+-------------------------------------------
+
+``colour.contrast``
+
+.. currentmodule:: colour.contrast
+
+.. autosummary::
+    :toctree: generated/
+
+    contrast_sensitivity_function_Barten1999
+
+**Ancillary Objects**
+
+``colour.contrast``
+
+.. currentmodule:: colour.contrast
+
+.. autosummary::
+    :toctree: generated/
+
+    optical_MTF_Barten1999
+    pupil_diameter_Barten1999
+    sigma_Barten1999
+    retinal_illuminance_Barten1999
+    maximum_angular_size_Barten1999
diff --git a/docs/generated/colour.CONTRAST_SENSITIVITY_METHODS.rst b/docs/generated/colour.CONTRAST_SENSITIVITY_METHODS.rst
new file mode 100644
index 000000000..d3c7b27d0
--- /dev/null
+++ b/docs/generated/colour.CONTRAST_SENSITIVITY_METHODS.rst
@@ -0,0 +1,6 @@
+colour.CONTRAST\_SENSITIVITY\_METHODS
+=====================================
+
+.. currentmodule:: colour
+
+.. autodata:: CONTRAST_SENSITIVITY_METHODS
\ No newline at end of file
diff --git a/docs/generated/colour.contrast.contrast_sensitivity_function_Barten1999.rst b/docs/generated/colour.contrast.contrast_sensitivity_function_Barten1999.rst
new file mode 100644
index 000000000..8e2a7d19e
--- /dev/null
+++ b/docs/generated/colour.contrast.contrast_sensitivity_function_Barten1999.rst
@@ -0,0 +1,6 @@
+colour.contrast.function\_contrast\_sensitivity\_Barten1999
+===========================================================
+
+.. currentmodule:: colour.contrast
+
+.. autofunction:: contrast_sensitivity_function_Barten1999
\ No newline at end of file
diff --git a/docs/generated/colour.contrast.maximum_angular_size_Barten1999.rst b/docs/generated/colour.contrast.maximum_angular_size_Barten1999.rst
new file mode 100644
index 000000000..6d3a8a2a4
--- /dev/null
+++ b/docs/generated/colour.contrast.maximum_angular_size_Barten1999.rst
@@ -0,0 +1,6 @@
+colour.contrast.maximum\_angular\_size\_Barten1999
+==================================================
+
+.. currentmodule:: colour.contrast
+
+.. autofunction:: maximum_angular_size_Barten1999
\ No newline at end of file
diff --git a/docs/generated/colour.contrast.optical_MTF_Barten1999.rst b/docs/generated/colour.contrast.optical_MTF_Barten1999.rst
new file mode 100644
index 000000000..940ebf126
--- /dev/null
+++ b/docs/generated/colour.contrast.optical_MTF_Barten1999.rst
@@ -0,0 +1,6 @@
+colour.contrast.optical\_MTF\_Barten1999
+========================================
+
+.. currentmodule:: colour.contrast
+
+.. autofunction:: optical_MTF_Barten1999
\ No newline at end of file
diff --git a/docs/generated/colour.contrast.pupil_diameter_Barten1999.rst b/docs/generated/colour.contrast.pupil_diameter_Barten1999.rst
new file mode 100644
index 000000000..3a3ab4b18
--- /dev/null
+++ b/docs/generated/colour.contrast.pupil_diameter_Barten1999.rst
@@ -0,0 +1,6 @@
+colour.contrast.pupil\_diameter\_Barten1999
+===========================================
+
+.. currentmodule:: colour.contrast
+
+.. autofunction:: pupil_diameter_Barten1999
\ No newline at end of file
diff --git a/docs/generated/colour.contrast.retinal_illuminance_Barten1999.rst b/docs/generated/colour.contrast.retinal_illuminance_Barten1999.rst
new file mode 100644
index 000000000..801fa9705
--- /dev/null
+++ b/docs/generated/colour.contrast.retinal_illuminance_Barten1999.rst
@@ -0,0 +1,6 @@
+colour.contrast.retinal\_illuminance\_Barten1999
+================================================
+
+.. currentmodule:: colour.contrast
+
+.. autofunction:: retinal_illuminance_Barten1999
\ No newline at end of file
diff --git a/docs/generated/colour.contrast.sigma_Barten1999.rst b/docs/generated/colour.contrast.sigma_Barten1999.rst
new file mode 100644
index 000000000..cb8368670
--- /dev/null
+++ b/docs/generated/colour.contrast.sigma_Barten1999.rst
@@ -0,0 +1,6 @@
+colour.contrast.sigma\_Barten1999
+=================================
+
+.. currentmodule:: colour.contrast
+
+.. autofunction:: sigma_Barten1999
\ No newline at end of file
diff --git a/docs/generated/colour.contrast_sensitivity_function.rst b/docs/generated/colour.contrast_sensitivity_function.rst
new file mode 100644
index 000000000..0ada5b5d2
--- /dev/null
+++ b/docs/generated/colour.contrast_sensitivity_function.rst
@@ -0,0 +1,6 @@
+colour.function\_contrast\_sensitivity
+======================================
+
+.. currentmodule:: colour
+
+.. autofunction:: contrast_sensitivity_function
\ No newline at end of file


Test patch:
diff --git a/colour/contrast/tests/__init__.py b/colour/contrast/tests/__init__.py
new file mode 100644
index 000000000..40a96afc6
--- /dev/null
+++ b/colour/contrast/tests/__init__.py
@@ -0,0 +1,1 @@
+# -*- coding: utf-8 -*-
diff --git a/colour/contrast/tests/test_barten1999.py b/colour/contrast/tests/test_barten1999.py
new file mode 100644
index 000000000..bd2141249
--- /dev/null
+++ b/colour/contrast/tests/test_barten1999.py
@@ -0,0 +1,546 @@
+# -*- coding: utf-8 -*-
+"""
+Defines unit tests for :mod:`colour.contrast.barten1999` module.
+"""
+
+from __future__ import division, unicode_literals
+
+import numpy as np
+import unittest
+from itertools import permutations
+
+from colour.contrast import (optical_MTF_Barten1999, pupil_diameter_Barten1999,
+                             sigma_Barten1999, retinal_illuminance_Barten1999,
+                             maximum_angular_size_Barten1999,
+                             contrast_sensitivity_function_Barten1999)
+from colour.utilities import ignore_numpy_errors
+
+__author__ = 'Colour Developers'
+__copyright__ = 'Copyright (C) 2013-2019 - Colour Developers'
+__license__ = 'New BSD License - http://opensource.org/licenses/BSD-3-Clause'
+__maintainer__ = 'Colour Developers'
+__email__ = 'colour-science@googlegroups.com'
+__status__ = 'Production'
+
+__all__ = [
+    'TestOpticalMTFBarten1999', 'TestPupilDiameterBarten1999',
+    'TestSigmaBarten1999', 'TestRetinalIlluminanceBarten1999',
+    'TestMaximumAngularSizeBarten1999',
+    'TestContrastSensitivityFunctionBarten1999'
+]
+
+
+class TestOpticalMTFBarten1999(unittest.TestCase):
+    """
+    Defines :func:`colour.contrast.barten1999.optical_MTF_Barten1999`
+    definition unit tests methods.
+    """
+
+    def test_optical_MTF_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.optical_MTF_Barten1999`
+        definition.
+        """
+
+        np.testing.assert_almost_equal(
+            optical_MTF_Barten1999(4, 0.01), 0.968910791191297, decimal=7)
+
+        np.testing.assert_almost_equal(
+            optical_MTF_Barten1999(8, 0.01), 0.881323136669471, decimal=7)
+
+        np.testing.assert_almost_equal(
+            optical_MTF_Barten1999(4, 0.05), 0.454040738727245, decimal=7)
+
+    def test_n_dimensional_optical_MTF_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.optical_MTF_Barten1999`
+        definition n-dimensions support.
+        """
+
+        u = np.array([4, 8, 12])
+        sigma = np.array([0.01, 0.05, 0.1])
+        M_opt = np.array([9.6891079e-01, 4.2499056e-02, 4.5228645e-13])
+        np.testing.assert_almost_equal(
+            optical_MTF_Barten1999(u, sigma), M_opt, decimal=7)
+
+        u = np.tile(u, (6, 1))
+        sigma = np.tile(sigma, (6, 1))
+        M_opt = np.tile(M_opt, (6, 1))
+        np.testing.assert_almost_equal(
+            optical_MTF_Barten1999(u, sigma), M_opt, decimal=7)
+
+        u = np.reshape(u, (2, 3, 3))
+        sigma = np.reshape(sigma, (2, 3, 3))
+        M_opt = np.reshape(M_opt, (2, 3, 3))
+        np.testing.assert_almost_equal(
+            optical_MTF_Barten1999(u, sigma), M_opt, decimal=7)
+
+    @ignore_numpy_errors
+    def test_nan_optical_MTF_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.optical_MTF_Barten1999`
+        definition nan support.
+        """
+
+        cases = [-1.0, 0.0, 1.0, -np.inf, np.inf, np.nan]
+        cases = set(permutations(cases * 3, r=3))
+        for case in cases:
+            optical_MTF_Barten1999(np.array(case), np.array(case))
+
+
+class TestPupilDiameterBarten1999(unittest.TestCase):
+    """
+    Defines :func:`colour.contrast.barten1999.pupil_diameter_Barten1999`
+    definition unit tests methods.
+    """
+
+    def test_pupil_diameter_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.pupil_diameter_Barten1999`
+        definition.
+        """
+
+        np.testing.assert_almost_equal(
+            pupil_diameter_Barten1999(20, 60), 2.272517118855717, decimal=7)
+
+        np.testing.assert_almost_equal(
+            pupil_diameter_Barten1999(0.2, 600), 2.272517118855717, decimal=7)
+
+        np.testing.assert_almost_equal(
+            pupil_diameter_Barten1999(20, 60, 30),
+            2.459028745178825,
+            decimal=7)
+
+    def test_n_dimensional_pupil_diameter_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.pupil_diameter_Barten1999`
+        definition n-dimensions support.
+        """
+
+        L = np.array([0.2, 20, 100])
+        X_0 = np.array([60, 120, 240])
+        Y_0 = np.array([60, 30, 15])
+        d = np.array([5.92690125, 2.27251712, 2.07775713])
+        np.testing.assert_almost_equal(
+            pupil_diameter_Barten1999(L, X_0, Y_0), d, decimal=7)
+
+        L = np.tile(L, (6, 1))
+        X_0 = np.tile(X_0, (6, 1))
+        d = np.tile(d, (6, 1))
+        np.testing.assert_almost_equal(
+            pupil_diameter_Barten1999(L, X_0, Y_0), d, decimal=7)
+
+        L = np.reshape(L, (2, 3, 3))
+        X_0 = np.reshape(X_0, (2, 3, 3))
+        d = np.reshape(d, (2, 3, 3))
+        np.testing.assert_almost_equal(
+            pupil_diameter_Barten1999(L, X_0, Y_0), d, decimal=7)
+
+    @ignore_numpy_errors
+    def test_nan_pupil_diameter_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.pupil_diameter_Barten1999`
+        definition nan support.
+        """
+
+        cases = [-1.0, 0.0, 1.0, -np.inf, np.inf, np.nan]
+        cases = set(permutations(cases * 3, r=3))
+        for case in cases:
+            pupil_diameter_Barten1999(
+                np.array(case), np.array(case), np.array(case))
+
+
+class TestSigmaBarten1999(unittest.TestCase):
+    """
+    Defines :func:`colour.contrast.barten1999.sigma_Barten1999` definition unit
+    tests methods.
+    """
+
+    def test_sigma_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.sigma_Barten1999` definition.
+        """
+
+        np.testing.assert_almost_equal(
+            sigma_Barten1999(0.5 / 60, 0.08 / 60, 2.1),
+            0.008791157173231,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            sigma_Barten1999(0.75 / 60, 0.08 / 60, 2.1),
+            0.012809761902549,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            sigma_Barten1999(0.5 / 60, 0.16 / 60, 2.1),
+            0.010040141654601,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            sigma_Barten1999(0.5 / 60, 0.08 / 60, 2.5),
+            0.008975274678558,
+            decimal=7)
+
+    def test_n_dimensional_sigma_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.sigma_Barten1999` definition
+        n-dimensions support.
+        """
+
+        sigma_0 = np.array([0.25 / 60, 0.5 / 60, 0.75 / 60])
+        C_ab = np.array([0.04 / 60, 0.08 / 60, 0.16 / 60])
+        d = np.array([2.1, 2.5, 5.0])
+        sigma = np.array([0.00439558, 0.00897527, 0.01827643])
+        np.testing.assert_almost_equal(
+            sigma_Barten1999(sigma_0, C_ab, d), sigma, decimal=7)
+
+        sigma_0 = np.tile(sigma_0, (6, 1))
+        C_ab = np.tile(C_ab, (6, 1))
+        sigma = np.tile(sigma, (6, 1))
+        np.testing.assert_almost_equal(
+            sigma_Barten1999(sigma_0, C_ab, d), sigma, decimal=7)
+
+        sigma_0 = np.reshape(sigma_0, (2, 3, 3))
+        C_ab = np.reshape(C_ab, (2, 3, 3))
+        sigma = np.reshape(sigma, (2, 3, 3))
+        np.testing.assert_almost_equal(
+            sigma_Barten1999(sigma_0, C_ab, d), sigma, decimal=7)
+
+    @ignore_numpy_errors
+    def test_nan_sigma_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.sigma_Barten1999`
+        definition nan support.
+        """
+
+        cases = [-1.0, 0.0, 1.0, -np.inf, np.inf, np.nan]
+        cases = set(permutations(cases * 3, r=3))
+        for case in cases:
+            sigma_Barten1999(np.array(case), np.array(case), np.array(case))
+
+
+class TestRetinalIlluminanceBarten1999(unittest.TestCase):
+    """
+    Defines :func:`colour.contrast.barten1999.retinal_illuminance_Barten1999`
+    definition unit tests methods.
+    """
+
+    def test_retinal_illuminance_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.retinal_illuminance_Barten1999`
+        definition.
+        """
+
+        np.testing.assert_almost_equal(
+            retinal_illuminance_Barten1999(20, 2.1, True),
+            66.082316060529919,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            retinal_illuminance_Barten1999(20, 2.5, True),
+            91.815644777503664,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            retinal_illuminance_Barten1999(20, 2.1, False),
+            69.272118011654939,
+            decimal=7)
+
+    def test_n_dimensional_retinal_illuminance_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.retinal_illuminance_Barten1999`
+        definition n-dimensions support.
+        """
+
+        L = np.array([0.2, 20, 100])
+        d = np.array([2.1, 2.5, 5.0])
+        E = np.array([0.66082316, 91.81564478, 1493.69535909])
+        np.testing.assert_almost_equal(
+            retinal_illuminance_Barten1999(L, d), E, decimal=7)
+
+        L = np.tile(L, (6, 1))
+        d = np.tile(d, (6, 1))
+        E = np.tile(E, (6, 1))
+        np.testing.assert_almost_equal(
+            retinal_illuminance_Barten1999(L, d), E, decimal=7)
+
+        L = np.reshape(L, (2, 3, 3))
+        d = np.reshape(d, (2, 3, 3))
+        E = np.reshape(E, (2, 3, 3))
+        np.testing.assert_almost_equal(
+            retinal_illuminance_Barten1999(L, d), E, decimal=7)
+
+    @ignore_numpy_errors
+    def test_nan_retinal_illuminance_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.retinal_illuminance_Barten1999`
+        definition nan support.
+        """
+
+        cases = [-1.0, 0.0, 1.0, -np.inf, np.inf, np.nan]
+        cases = set(permutations(cases * 3, r=3))
+        for case in cases:
+            retinal_illuminance_Barten1999(np.array(case), np.array(case))
+
+
+class TestMaximumAngularSizeBarten1999(unittest.TestCase):
+    """
+    Defines :func:`colour.contrast.barten1999.maximum_angular_size_Barten1999`
+    definition unit tests methods.
+    """
+
+    def test_maximum_angular_size_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.\
+maximum_angular_size_Barten1999` definition.
+        """
+
+        np.testing.assert_almost_equal(
+            maximum_angular_size_Barten1999(4, 60, 12, 15),
+            3.572948005052482,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            maximum_angular_size_Barten1999(8, 60, 12, 15),
+            1.851640199545103,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            maximum_angular_size_Barten1999(4, 120, 12, 15),
+            3.577708763999663,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            maximum_angular_size_Barten1999(4, 60, 24, 15),
+            3.698001308168194,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            maximum_angular_size_Barten1999(4, 60, 12, 30),
+            6.324555320336758,
+            decimal=7)
+
+    def test_n_dimensional_maximum_angular_size_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.\
+maximum_angular_size_Barten1999` definition n-dimensions support.
+        """
+
+        u = np.array([4, 8, 12])
+        X_0 = np.array([60, 120, 240])
+        X_max = np.array([12, 14, 16])
+        N_max = np.array([15, 20, 25])
+        X = np.array([3.57294801, 2.46055152, 2.06581761])
+        np.testing.assert_almost_equal(
+            maximum_angular_size_Barten1999(u, X_0, X_max, N_max),
+            X,
+            decimal=7)
+
+        u = np.tile(u, (6, 1))
+        X_0 = np.tile(X_0, (6, 1))
+        X = np.tile(X, (6, 1))
+        np.testing.assert_almost_equal(
+            maximum_angular_size_Barten1999(u, X_0, X_max, N_max),
+            X,
+            decimal=7)
+
+        u = np.reshape(u, (2, 3, 3))
+        X_0 = np.reshape(X_0, (2, 3, 3))
+        X = np.reshape(X, (2, 3, 3))
+        np.testing.assert_almost_equal(
+            maximum_angular_size_Barten1999(u, X_0, X_max, N_max),
+            X,
+            decimal=7)
+
+    @ignore_numpy_errors
+    def test_nan_maximum_angular_size_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.\
+maximum_angular_size_Barten1999` definition nan support.
+        """
+
+        cases = [-1.0, 0.0, 1.0, -np.inf, np.inf, np.nan]
+        cases = set(permutations(cases * 3, r=3))
+        for case in cases:
+            maximum_angular_size_Barten1999(
+                np.array(case), np.array(case), np.array(case), np.array(case))
+
+
+class TestContrastSensitivityFunctionBarten1999(unittest.TestCase):
+    """
+    Defines :func:`colour.contrast.barten1999.\
+contrast_sensitivity_function_Barten1999` definition unit tests methods.
+    """
+
+    def test_contrast_sensitivity_function_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.\
+contrast_sensitivity_function_Barten1999` definition.
+        """
+
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=4,
+                sigma=0.01,
+                E=65,
+                X_0=60,
+                X_max=12,
+                Y_0=60,
+                Y_max=12,
+                p=1.2 * 10 ** 6),
+            352.761342126727020,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=8,
+                sigma=0.01,
+                E=65,
+                X_0=60,
+                X_max=12,
+                Y_0=60,
+                Y_max=12,
+                p=1.2 * 10 ** 6),
+            177.706338840717340,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=4,
+                sigma=0.02,
+                E=65,
+                X_0=60,
+                X_max=12,
+                Y_0=60,
+                Y_max=12,
+                p=1.2 * 10 ** 6),
+            320.872401634215750,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=4,
+                sigma=0.01,
+                E=130,
+                X_0=60,
+                X_max=12,
+                Y_0=60,
+                Y_max=12,
+                p=1.2 * 10 ** 6),
+            455.171315756946400,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=4,
+                sigma=0.01,
+                E=65,
+                X_0=120,
+                X_max=12,
+                Y_0=60,
+                Y_max=12,
+                p=1.2 * 10 ** 6),
+            352.996281545740660,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=4,
+                sigma=0.01,
+                E=65,
+                X_0=60,
+                X_max=24,
+                Y_0=60,
+                Y_max=12,
+                p=1.2 * 10 ** 6),
+            358.881580104493650,
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=4,
+                sigma=0.01,
+                E=65,
+                X_0=240,
+                X_max=12,
+                Y_0=60,
+                Y_max=12,
+                p=1.2 * 10 ** 6),
+            contrast_sensitivity_function_Barten1999(
+                u=4,
+                sigma=0.01,
+                E=65,
+                X_0=60,
+                X_max=12,
+                Y_0=240,
+                Y_max=12,
+                p=1.2 * 10 ** 6),
+            decimal=7)
+
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=4,
+                sigma=0.01,
+                E=65,
+                X_0=60,
+                X_max=12,
+                Y_0=60,
+                Y_max=12,
+                p=1.4 * 10 ** 6),
+            374.791328640476140,
+            decimal=7)
+
+    def test_n_dimensional_contrast_sensitivity_function_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.\
+contrast_sensitivity_function_Barten1999` definition n-dimensions support.
+        """
+
+        u = np.array([4, 8, 12])
+        sigma = np.array([0.01, 0.02, 0.04])
+        E = np.array([0.65, 90, 1500])
+        X_0 = np.array([60, 120, 240])
+        S = np.array([39.86810824, 142.24417024, 4.39800166])
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=u, sigma=sigma, E=E, X_0=X_0),
+            S,
+            decimal=7)
+
+        u = np.tile(u, (6, 1))
+        E = np.tile(E, (6, 1))
+        S = np.tile(S, (6, 1))
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=u, sigma=sigma, E=E, X_0=X_0),
+            S,
+            decimal=7)
+
+        u = np.reshape(u, (2, 3, 3))
+        E = np.reshape(E, (2, 3, 3))
+        S = np.reshape(S, (2, 3, 3))
+        np.testing.assert_almost_equal(
+            contrast_sensitivity_function_Barten1999(
+                u=u, sigma=sigma, E=E, X_0=X_0),
+            S,
+            decimal=7)
+
+    @ignore_numpy_errors
+    def test_nan_contrast_sensitivity_function_Barten1999(self):
+        """
+        Tests :func:`colour.contrast.barten1999.\
+contrast_sensitivity_function_Barten1999` definition nan support.
+        """
+
+        cases = [-1.0, 0.0, 1.0, -np.inf, np.inf, np.nan]
+        cases = set(permutations(cases * 3, r=3))
+        for case in cases:
+            contrast_sensitivity_function_Barten1999(
+                u=np.array(case),
+                sigma=np.array(case),
+                E=np.array(case),
+                X_0=np.array(case))
+
+
+if __name__ == '__main__':
+    unittest.main()