"""
Jakob and Hanika (2019) - Reflectance Recovery
==============================================
Define the objects for reflectance recovery, i.e., spectral upsampling, using
*Jakob and Hanika (2019)* method:
- :func:`colour.recovery.sd_Jakob2019`
- :func:`colour.recovery.find_coefficients_Jakob2019`
- :func:`colour.recovery.XYZ_to_sd_Jakob2019`
- :class:`colour.recovery.LUT3D_Jakob2019`
References
----------
- :cite:`Jakob2019` : Jakob, W., & Hanika, J. (2019). A Low-Dimensional
Function Space for Efficient Spectral Upsampling. Computer Graphics Forum,
38(2), 147-155. doi:10.1111/cgf.13626
"""
from __future__ import annotations
import struct
from pathlib import Path
import numpy as np
from scipy.interpolate import RegularGridInterpolator
from scipy.optimize import minimize
from colour.algebra import smoothstep_function, spow
from colour.colorimetry import (
MultiSpectralDistributions,
SpectralDistribution,
SpectralShape,
handle_spectral_arguments,
intermediate_lightness_function_CIE1976,
sd_to_XYZ_integration,
)
from colour.constants import DTYPE_INT_DEFAULT
from colour.difference import JND_CIE1976
from colour.hints import (
ArrayLike,
Callable,
NDArrayFloat,
Tuple,
)
from colour.models import RGB_Colourspace, RGB_to_XYZ, XYZ_to_Lab, XYZ_to_xy
from colour.utilities import (
as_float_array,
as_float_scalar,
domain_range_scale,
full,
index_along_last_axis,
is_tqdm_installed,
message_box,
optional,
to_domain_1,
tsplit,
zeros,
)
if is_tqdm_installed():
from tqdm import tqdm
else: # pragma: no cover
from unittest import mock
tqdm = mock.MagicMock()
__author__ = "Colour Developers"
__copyright__ = "Copyright 2013 Colour Developers"
__license__ = "BSD-3-Clause - https://opensource.org/licenses/BSD-3-Clause"
__maintainer__ = "Colour Developers"
__email__ = "colour-developers@colour-science.org"
__status__ = "Production"
__all__ = [
"SPECTRAL_SHAPE_JAKOB2019",
"StopMinimizationEarlyError",
"sd_Jakob2019",
"error_function",
"dimensionalise_coefficients",
"lightness_scale",
"find_coefficients_Jakob2019",
"XYZ_to_sd_Jakob2019",
"LUT3D_Jakob2019",
]
SPECTRAL_SHAPE_JAKOB2019: SpectralShape = SpectralShape(360, 780, 5)
"""Spectral shape for *Jakob and Hanika (2019)* method."""
class StopMinimizationEarlyError(Exception):
"""
The exception used to stop :func:`scipy.optimize.minimize` once the
value of the minimized function is small enough. *SciPy* doesn't currently
offer a better way of doing it.
Attributes
----------
- :attr:`~colour.recovery.jakob2019.StopMinimizationEarlyError.coefficients`
- :attr:`~colour.recovery.jakob2019.StopMinimizationEarlyError.error`
"""
def __init__(self, coefficients: ArrayLike, error: float) -> None:
self._coefficients = as_float_array(coefficients)
self._error = as_float_scalar(error)
@property
def coefficients(self) -> NDArrayFloat:
"""
Getter property for the *Jakob and Hanika (2019)* exception
coefficients.
Returns
-------
:class:`numpy.ndarray`
*Jakob and Hanika (2019)* exception coefficients.
"""
return self._coefficients
@property
def error(self) -> float:
"""
Getter property for the *Jakob and Hanika (2019)* exception error
value.
Returns
-------
:class:`float`
*Jakob and Hanika (2019)* exception coefficients.
"""
return self._error
[docs]
def sd_Jakob2019(
coefficients: ArrayLike, shape: SpectralShape = SPECTRAL_SHAPE_JAKOB2019
) -> SpectralDistribution:
"""
Return a spectral distribution following the spectral model given by
*Jakob and Hanika (2019)*.
Parameters
----------
coefficients
Dimensionless coefficients for *Jakob and Hanika (2019)* reflectance
spectral model.
shape
Shape used by the spectral distribution.
Returns
-------
:class:`colour.SpectralDistribution`
*Jakob and Hanika (2019)* spectral distribution.
References
----------
:cite:`Jakob2019`
Examples
--------
>>> from colour.utilities import numpy_print_options
>>> with numpy_print_options(suppress=True):
... sd_Jakob2019([-9e-05, 8.5e-02, -20], SpectralShape(400, 700, 20))
... # doctest: +ELLIPSIS
SpectralDistribution([[ 400. , 0.3143046...],
[ 420. , 0.4133320...],
[ 440. , 0.4880034...],
[ 460. , 0.5279562...],
[ 480. , 0.5319346...],
[ 500. , 0.5 ...],
[ 520. , 0.4326202...],
[ 540. , 0.3373544...],
[ 560. , 0.2353056...],
[ 580. , 0.1507665...],
[ 600. , 0.0931332...],
[ 620. , 0.0577434...],
[ 640. , 0.0367011...],
[ 660. , 0.0240879...],
[ 680. , 0.0163316...],
[ 700. , 0.0114118...]],
SpragueInterpolator,
{},
Extrapolator,
{'method': 'Constant', 'left': None, 'right': None})
"""
c_0, c_1, c_2 = as_float_array(coefficients)
wl = shape.wavelengths
U = c_0 * wl**2 + c_1 * wl + c_2
R = 1 / 2 + U / (2 * np.sqrt(1 + U**2))
name = f"{coefficients!r} (COEFF) - Jakob (2019)"
return SpectralDistribution(R, wl, name=name)
def error_function(
coefficients: ArrayLike,
target: ArrayLike,
cmfs: MultiSpectralDistributions,
illuminant: SpectralDistribution,
max_error: float | None = None,
additional_data: bool = False,
) -> (
Tuple[float, NDArrayFloat]
| Tuple[float, NDArrayFloat, NDArrayFloat, NDArrayFloat, NDArrayFloat]
):
"""
Compute :math:`\\Delta E_{76}` between the target colour and the colour
defined by given spectral model, along with its gradient.
Parameters
----------
coefficients
Dimensionless coefficients for *Jakob and Hanika (2019)* reflectance
spectral model.
target
*CIE L\\*a\\*b\\** colourspace array of the target colour.
cmfs
Standard observer colour matching functions.
illuminant
Illuminant spectral distribution.
max_error
Raise ``StopMinimizationEarlyError`` if the error is smaller than this.
The default is *None* and the function doesn't raise anything.
additional_data
If *True*, some intermediate calculations are returned, for use in
correctness tests: R, XYZ and Lab.
Returns
-------
:class:`tuple` or :class:`tuple`
Tuple of computed :math:`\\Delta E_{76}` error and gradient of error,
i.e., the first derivatives of error with respect to the input
coefficients or tuple of computed :math:`\\Delta E_{76}` error,
gradient of error, computed spectral reflectance, *CIE XYZ* tristimulus
values corresponding to ``R`` and *CIE L\\*a\\*b\\** colourspace array
corresponding to ``XYZ``.
Raises
------
StopMinimizationEarlyError
Raised when the error is below ``max_error``.
"""
target = as_float_array(target)
c_0, c_1, c_2 = as_float_array(coefficients)
wv = np.linspace(0, 1, len(cmfs.shape))
U = c_0 * wv**2 + c_1 * wv + c_2
t1 = np.sqrt(1 + U**2)
R = 1 / 2 + U / (2 * t1)
t2 = 1 / (2 * t1) - U**2 / (2 * t1**3)
dR = np.array([wv**2 * t2, wv * t2, t2])
XYZ = sd_to_XYZ_integration(R, cmfs, illuminant, shape=cmfs.shape) / 100
dXYZ = np.transpose(
sd_to_XYZ_integration(dR, cmfs, illuminant, shape=cmfs.shape) / 100
)
XYZ_n = sd_to_XYZ_integration(illuminant, cmfs)
XYZ_n /= XYZ_n[1]
XYZ_XYZ_n = XYZ / XYZ_n
XYZ_f = intermediate_lightness_function_CIE1976(XYZ, XYZ_n)
dXYZ_f = np.where(
XYZ_XYZ_n[..., None] > (24 / 116) ** 3,
1 / (3 * spow(XYZ_n[..., None], 1 / 3) * spow(XYZ[..., None], 2 / 3)) * dXYZ,
(841 / 108) * dXYZ / XYZ_n[..., None],
)
def intermediate_XYZ_to_Lab(
XYZ_i: NDArrayFloat, offset: float | None = 16
) -> NDArrayFloat:
"""
Return the final intermediate value for the *CIE Lab* to *CIE XYZ*
conversion.
"""
return np.array(
[
116 * XYZ_i[1] - offset,
500 * (XYZ_i[0] - XYZ_i[1]),
200 * (XYZ_i[1] - XYZ_i[2]),
]
)
Lab_i = intermediate_XYZ_to_Lab(XYZ_f)
dLab_i = intermediate_XYZ_to_Lab(dXYZ_f, 0)
error = np.sqrt(np.sum((Lab_i - target) ** 2))
if max_error is not None and error <= max_error:
raise StopMinimizationEarlyError(coefficients, error)
derror = np.sum(dLab_i * (Lab_i[..., None] - target[..., None]), axis=0) / error
if additional_data:
return error, derror, R, XYZ, Lab_i
else:
return error, derror
def dimensionalise_coefficients(
coefficients: ArrayLike, shape: SpectralShape
) -> NDArrayFloat:
"""
Rescale the dimensionless coefficients to given spectral shape.
A dimensionless form of the reflectance spectral model is used in the
optimisation process. Instead of the usual spectral shape, specified in
nanometers, it is normalised to the [0, 1] range. A side effect is that
computed coefficients work only with the normalised range and need to be
rescaled to regain units and be compatible with standard shapes.
Parameters
----------
coefficients
Dimensionless coefficients.
shape
Spectral distribution shape used in calculations.
Returns
-------
:class:`numpy.ndarray`
Dimensionful coefficients, with units of
:math:`\\frac{1}{\\mathrm{nm}^2}`, :math:`\\frac{1}{\\mathrm{nm}}`
and 1, respectively.
"""
cp_0, cp_1, cp_2 = tsplit(coefficients)
span = shape.end - shape.start
c_0 = cp_0 / span**2
c_1 = cp_1 / span - 2 * cp_0 * shape.start / span**2
c_2 = cp_0 * shape.start**2 / span**2 - cp_1 * shape.start / span + cp_2
return np.array([c_0, c_1, c_2])
def lightness_scale(steps: int) -> NDArrayFloat:
"""
Create a non-linear lightness scale, as described in *Jakob and Hanika
(2019)*. The spacing between very dark and very bright (and saturated)
colours is made smaller, because in those regions coefficients tend to
change rapidly and a finer resolution is needed.
Parameters
----------
steps
Samples/steps count along the non-linear lightness scale.
Returns
-------
:class:`numpy.ndarray`
Non-linear lightness scale.
Examples
--------
>>> lightness_scale(5) # doctest: +ELLIPSIS
array([ 0. , 0.0656127..., 0.5 , 0.9343872..., \
1. ])
"""
linear = np.linspace(0, 1, steps)
return smoothstep_function(smoothstep_function(linear))
[docs]
def find_coefficients_Jakob2019(
XYZ: ArrayLike,
cmfs: MultiSpectralDistributions | None = None,
illuminant: SpectralDistribution | None = None,
coefficients_0: ArrayLike = zeros(3),
max_error: float = JND_CIE1976 / 100,
dimensionalise: bool = True,
) -> Tuple[NDArrayFloat, float]:
"""
Compute the coefficients for *Jakob and Hanika (2019)* reflectance
spectral model.
Parameters
----------
XYZ
*CIE XYZ* tristimulus values to find the coefficients for.
cmfs
Standard observer colour matching functions, default to the
*CIE 1931 2 Degree Standard Observer*.
illuminant
Illuminant spectral distribution, default to
*CIE Standard Illuminant D65*.
coefficients_0
Starting coefficients for the solver.
max_error
Maximal acceptable error. Set higher to save computational time.
If *None*, the solver will keep going until it is very close to the
minimum. The default is ``ACCEPTABLE_DELTA_E``.
dimensionalise
If *True*, returned coefficients are dimensionful and will not work
correctly if fed back as ``coefficients_0``. The default is *True*.
Returns
-------
:class:`tuple`
Tuple of computed coefficients that best fit the given colour and
:math:`\\Delta E_{76}` between the target colour and the colour
corresponding to the computed coefficients.
References
----------
:cite:`Jakob2019`
Examples
--------
>>> XYZ = np.array([0.20654008, 0.12197225, 0.05136952])
>>> find_coefficients_Jakob2019(XYZ) # doctest: +ELLIPSIS
(array([ 1.3723791...e-04, -1.3514399...e-01, 3.0838973...e+01]), \
0.0141941...)
"""
XYZ = as_float_array(XYZ)
coefficients_0 = as_float_array(coefficients_0)
cmfs, illuminant = handle_spectral_arguments(
cmfs, illuminant, shape_default=SPECTRAL_SHAPE_JAKOB2019
)
def optimize(
target_o: ArrayLike, coefficients_0_o: ArrayLike
) -> Tuple[NDArrayFloat, float]:
"""Minimise the error function using *L-BFGS-B* method."""
try:
result = minimize(
error_function,
coefficients_0_o,
(target_o, cmfs, illuminant, max_error),
method="L-BFGS-B",
jac=True,
)
return result.x, result.fun
except StopMinimizationEarlyError as error:
return error.coefficients, error.error
xy_n = XYZ_to_xy(sd_to_XYZ_integration(illuminant, cmfs))
XYZ_g = full(3, 0.5)
coefficients_g = zeros(3)
divisions = 3
while divisions < 10:
XYZ_r = XYZ_g
coefficient_r = coefficients_g
keep_divisions = False
coefficients_0 = coefficient_r
for i in range(1, divisions):
XYZ_i = (XYZ - XYZ_r) * i / (divisions - 1) + XYZ_r
Lab_i = XYZ_to_Lab(XYZ_i)
coefficients_0, error = optimize(Lab_i, coefficients_0)
if error > max_error:
break
else:
XYZ_g = XYZ_i
coefficients_g = coefficients_0
keep_divisions = True
else:
break
if not keep_divisions:
divisions += 2
target = XYZ_to_Lab(XYZ, xy_n)
coefficients, error = optimize(target, coefficients_0)
if dimensionalise:
coefficients = dimensionalise_coefficients(coefficients, cmfs.shape)
return coefficients, error
[docs]
def XYZ_to_sd_Jakob2019(
XYZ: ArrayLike,
cmfs: MultiSpectralDistributions | None = None,
illuminant: SpectralDistribution | None = None,
optimisation_kwargs: dict | None = None,
additional_data: bool = False,
) -> Tuple[SpectralDistribution, float] | SpectralDistribution:
"""
Recover the spectral distribution of given *CIE XYZ* tristimulus values
using *Jakob and Hanika (2019)* method.
Parameters
----------
XYZ
*CIE XYZ* tristimulus values to recover the spectral distribution from.
cmfs
Standard observer colour matching functions, default to the
*CIE 1931 2 Degree Standard Observer*.
illuminant
Illuminant spectral distribution, default to
*CIE Standard Illuminant D65*.
optimisation_kwargs
Parameters for :func:`colour.recovery.find_coefficients_Jakob2019`
definition.
additional_data
If *True*, ``error`` will be returned alongside the recovered spectral
distribution.
Returns
-------
:class:`tuple` or :class:`colour.SpectralDistribution`
Tuple of recovered spectral distribution and :math:`\\Delta E_{76}`
between the target colour and the colour corresponding to the computed
coefficients or recovered spectral distribution.
References
----------
:cite:`Jakob2019`
Examples
--------
>>> from colour import (
... CCS_ILLUMINANTS,
... MSDS_CMFS,
... SDS_ILLUMINANTS,
... XYZ_to_sRGB,
... )
>>> from colour.colorimetry import sd_to_XYZ_integration
>>> from colour.utilities import numpy_print_options
>>> XYZ = np.array([0.20654008, 0.12197225, 0.05136952])
>>> cmfs = (
... MSDS_CMFS["CIE 1931 2 Degree Standard Observer"]
... .copy()
... .align(SpectralShape(360, 780, 10))
... )
>>> illuminant = SDS_ILLUMINANTS["D65"].copy().align(cmfs.shape)
>>> sd = XYZ_to_sd_Jakob2019(XYZ, cmfs, illuminant)
>>> with numpy_print_options(suppress=True):
... sd # doctest: +ELLIPSIS
SpectralDistribution([[ 360. , 0.4893773...],
[ 370. , 0.3258214...],
[ 380. , 0.2147792...],
[ 390. , 0.1482413...],
[ 400. , 0.1086169...],
[ 410. , 0.0841255...],
[ 420. , 0.0683114...],
[ 430. , 0.0577144...],
[ 440. , 0.0504267...],
[ 450. , 0.0453552...],
[ 460. , 0.0418520...],
[ 470. , 0.0395259...],
[ 480. , 0.0381430...],
[ 490. , 0.0375741...],
[ 500. , 0.0377685...],
[ 510. , 0.0387432...],
[ 520. , 0.0405871...],
[ 530. , 0.0434783...],
[ 540. , 0.0477225...],
[ 550. , 0.0538256...],
[ 560. , 0.0626314...],
[ 570. , 0.0755869...],
[ 580. , 0.0952675...],
[ 590. , 0.1264265...],
[ 600. , 0.1779272...],
[ 610. , 0.2649393...],
[ 620. , 0.4039779...],
[ 630. , 0.5832105...],
[ 640. , 0.7445440...],
[ 650. , 0.8499970...],
[ 660. , 0.9094792...],
[ 670. , 0.9425378...],
[ 680. , 0.9616376...],
[ 690. , 0.9732481...],
[ 700. , 0.9806562...],
[ 710. , 0.9855873...],
[ 720. , 0.9889903...],
[ 730. , 0.9914117...],
[ 740. , 0.9931801...],
[ 750. , 0.9945009...],
[ 760. , 0.9955066...],
[ 770. , 0.9962855...],
[ 780. , 0.9968976...]],
SpragueInterpolator,
{},
Extrapolator,
{'method': 'Constant', 'left': None, 'right': None})
>>> sd_to_XYZ_integration(sd, cmfs, illuminant) / 100 # doctest: +ELLIPSIS
array([ 0.2066217..., 0.1220128..., 0.0513958...])
"""
XYZ = to_domain_1(XYZ)
cmfs, illuminant = handle_spectral_arguments(
cmfs, illuminant, shape_default=SPECTRAL_SHAPE_JAKOB2019
)
optimisation_kwargs = optional(optimisation_kwargs, {})
with domain_range_scale("ignore"):
coefficients, error = find_coefficients_Jakob2019(
XYZ, cmfs, illuminant, **optimisation_kwargs
)
sd = sd_Jakob2019(coefficients, cmfs.shape)
sd.name = f"{XYZ} (XYZ) - Jakob (2019)"
if additional_data:
return sd, error
else:
return sd
[docs]
class LUT3D_Jakob2019:
"""
Define a class for working with pre-computed lookup tables for the
*Jakob and Hanika (2019)* spectral upsampling method. It allows significant
time savings by performing the expensive numerical optimisation ahead of
time and storing the results in a file.
The file format is compatible with the code and *\\*.coeff* files in the
supplemental material published alongside the article. They are directly
available from
`Colour - Datasets <https://github.com/colour-science/colour-datasets>`__
under the record *4050598*.
Attributes
----------
- :attr:`~colour.recovery.LUT3D_Jakob2019.size`
- :attr:`~colour.recovery.LUT3D_Jakob2019.lightness_scale`
- :attr:`~colour.recovery.LUT3D_Jakob2019.coefficients`
- :attr:`~colour.recovery.LUT3D_Jakob2019.interpolator`
Methods
-------
- :meth:`~colour.recovery.LUT3D_Jakob2019.__init__`
- :meth:`~colour.recovery.LUT3D_Jakob2019.generate`
- :meth:`~colour.recovery.LUT3D_Jakob2019.RGB_to_coefficients`
- :meth:`~colour.recovery.LUT3D_Jakob2019.RGB_to_sd`
- :meth:`~colour.recovery.LUT3D_Jakob2019.read`
- :meth:`~colour.recovery.LUT3D_Jakob2019.write`
References
----------
:cite:`Jakob2019`
Examples
--------
>>> import os
>>> import colour
>>> from colour import CCS_ILLUMINANTS, SDS_ILLUMINANTS, MSDS_CMFS
>>> from colour.colorimetry import sd_to_XYZ_integration
>>> from colour.models import RGB_COLOURSPACE_sRGB
>>> from colour.utilities import numpy_print_options
>>> cmfs = (
... MSDS_CMFS["CIE 1931 2 Degree Standard Observer"]
... .copy()
... .align(SpectralShape(360, 780, 10))
... )
>>> illuminant = SDS_ILLUMINANTS["D65"].copy().align(cmfs.shape)
>>> LUT = LUT3D_Jakob2019()
>>> LUT.generate(RGB_COLOURSPACE_sRGB, cmfs, illuminant, 3, lambda x: x)
>>> path = os.path.join(
... colour.__path__[0],
... "recovery",
... "tests",
... "resources",
... "sRGB_Jakob2019.coeff",
... )
>>> LUT.write(path) # doctest: +SKIP
>>> LUT.read(path) # doctest: +SKIP
>>> RGB = np.array([0.70573936, 0.19248266, 0.22354169])
>>> with numpy_print_options(suppress=True):
... LUT.RGB_to_sd(RGB, cmfs.shape) # doctest: +ELLIPSIS
SpectralDistribution([[ 360. , 0.7666803...],
[ 370. , 0.6251547...],
[ 380. , 0.4584310...],
[ 390. , 0.3161633...],
[ 400. , 0.2196155...],
[ 410. , 0.1596575...],
[ 420. , 0.1225525...],
[ 430. , 0.0989784...],
[ 440. , 0.0835782...],
[ 450. , 0.0733535...],
[ 460. , 0.0666049...],
[ 470. , 0.0623569...],
[ 480. , 0.06006 ...],
[ 490. , 0.0594383...],
[ 500. , 0.0604201...],
[ 510. , 0.0631195...],
[ 520. , 0.0678648...],
[ 530. , 0.0752834...],
[ 540. , 0.0864790...],
[ 550. , 0.1033773...],
[ 560. , 0.1293883...],
[ 570. , 0.1706018...],
[ 580. , 0.2374178...],
[ 590. , 0.3439472...],
[ 600. , 0.4950548...],
[ 610. , 0.6604253...],
[ 620. , 0.7914669...],
[ 630. , 0.8738724...],
[ 640. , 0.9213216...],
[ 650. , 0.9486880...],
[ 660. , 0.9650550...],
[ 670. , 0.9752838...],
[ 680. , 0.9819499...],
[ 690. , 0.9864585...],
[ 700. , 0.9896073...],
[ 710. , 0.9918680...],
[ 720. , 0.9935302...],
[ 730. , 0.9947778...],
[ 740. , 0.9957312...],
[ 750. , 0.9964714...],
[ 760. , 0.9970543...],
[ 770. , 0.9975190...],
[ 780. , 0.9978936...]],
SpragueInterpolator,
{},
Extrapolator,
{'method': 'Constant', 'left': None, 'right': None})
"""
[docs]
def __init__(self) -> None:
self._interpolator: RegularGridInterpolator = RegularGridInterpolator(
np.array([]), np.array([])
)
self._size: int = 0
self._lightness_scale: NDArrayFloat = np.array([])
self._coefficients: NDArrayFloat = np.array([])
@property
def size(self) -> int:
"""
Getter property for the *Jakob and Hanika (2019)* interpolator
size, i.e., the sample count on one side of the 3D table.
Returns
-------
:class:`int`
*Jakob and Hanika (2019)* interpolator size.
"""
return self._size
@property
def lightness_scale(self) -> NDArrayFloat:
"""
Getter property for the *Jakob and Hanika (2019)* interpolator
lightness scale.
Returns
-------
:class:`numpy.ndarray`
*Jakob and Hanika (2019)* interpolator lightness scale.
"""
return self._lightness_scale
@property
def coefficients(self) -> NDArrayFloat:
"""
Getter property for the *Jakob and Hanika (2019)* interpolator
coefficients.
Returns
-------
:class:`numpy.ndarray`
*Jakob and Hanika (2019)* interpolator coefficients.
"""
return self._coefficients
@property
def interpolator(self) -> RegularGridInterpolator:
"""
Getter property for the *Jakob and Hanika (2019)* interpolator.
Returns
-------
:class:`scipy.interpolate.RegularGridInterpolator`
*Jakob and Hanika (2019)* interpolator.
"""
return self._interpolator
def _create_interpolator(self):
"""
Create a :class:`scipy.interpolate.RegularGridInterpolator` class
instance for read or generated coefficients.
"""
samples = np.linspace(0, 1, self._size)
axes = ([0, 1, 2], self._lightness_scale, samples, samples)
self._interpolator = RegularGridInterpolator(
axes, self._coefficients, bounds_error=False
)
def generate(
self,
colourspace: RGB_Colourspace,
cmfs: MultiSpectralDistributions | None = None,
illuminant: SpectralDistribution | None = None,
size: int = 64,
print_callable: Callable = print,
):
"""
Generate the lookup table data for given *RGB* colourspace, colour
matching functions, illuminant and given size.
Parameters
----------
colourspace
The *RGB* colourspace to create a lookup table for.
cmfs
Standard observer colour matching functions, default to the
*CIE 1931 2 Degree Standard Observer*.
illuminant
Illuminant spectral distribution, default to
*CIE Standard Illuminant D65*.
size
The resolution of the lookup table. Higher values will decrease
errors but at the cost of a much longer run time. The published
*\\*.coeff* files have a resolution of 64.
print_callable
Callable used to print progress and diagnostic information.
Examples
--------
>>> from colour import MSDS_CMFS, SDS_ILLUMINANTS
>>> from colour.models import RGB_COLOURSPACE_sRGB
>>> from colour.utilities import numpy_print_options
>>> cmfs = (
... MSDS_CMFS["CIE 1931 2 Degree Standard Observer"]
... .copy()
... .align(SpectralShape(360, 780, 10))
... )
>>> illuminant = SDS_ILLUMINANTS["D65"].copy().align(cmfs.shape)
>>> LUT = LUT3D_Jakob2019()
>>> print(LUT.interpolator) # doctest: +ELLIPSIS
<scipy...RegularGridInterpolator object at 0x...>
>>> LUT.generate(RGB_COLOURSPACE_sRGB, cmfs, illuminant, 3)
======================================================================\
=========
* \
*
* "Jakob et al. (2018)" LUT Optimisation \
*
* \
*
======================================================================\
=========
<BLANKLINE>
Optimising 27 coefficients...
<BLANKLINE>
>>> print(LUT.interpolator)
... # doctest: +ELLIPSIS
<scipy.interpolate...RegularGridInterpolator object at 0x...>
"""
cmfs, illuminant = handle_spectral_arguments(
cmfs, illuminant, shape_default=SPECTRAL_SHAPE_JAKOB2019
)
shape = cmfs.shape
xy_n = XYZ_to_xy(sd_to_XYZ_integration(illuminant, cmfs))
# It could be interesting to have different resolutions for lightness
# and chromaticity, but the current file format doesn't allow it.
lightness_steps = size
chroma_steps = size
self._lightness_scale = lightness_scale(lightness_steps)
self._coefficients = np.empty(
[3, chroma_steps, chroma_steps, lightness_steps, 3]
)
cube_indexes = np.ndindex(3, chroma_steps, chroma_steps)
total_coefficients = chroma_steps**2 * 3
# First, create a list of all the fully bright colours with the order
# matching cube_indexes.
samples = np.linspace(0, 1, chroma_steps)
ij = np.reshape(
np.transpose(np.meshgrid([1], samples, samples, indexing="ij")),
(-1, 3),
)
chromas = np.concatenate(
[
ij,
np.roll(ij, 1, axis=1),
np.roll(ij, 2, axis=1),
]
)
message_box(
'"Jakob et al. (2018)" LUT Optimisation',
print_callable=print_callable,
)
print_callable(f"\nOptimising {total_coefficients} coefficients...\n")
def optimize(
ijkL: ArrayLike, coefficients_0: ArrayLike, chroma: NDArrayFloat
) -> NDArrayFloat:
"""
Solve for a specific lightness and stores the result in the
appropriate cell.
"""
i, j, k, L = tsplit(ijkL, dtype=DTYPE_INT_DEFAULT)
RGB = self._lightness_scale[L] * chroma
XYZ = RGB_to_XYZ(RGB, colourspace, xy_n)
coefficients, _error = find_coefficients_Jakob2019(
XYZ, cmfs, illuminant, coefficients_0, dimensionalise=False
)
self._coefficients[i, L, j, k, :] = dimensionalise_coefficients(
coefficients, shape
)
return coefficients
with tqdm(total=total_coefficients) as progress:
for ijk, chroma in zip(cube_indexes, chromas):
progress.update()
# Starts from somewhere in the middle, similarly to how
# feedback works in "colour.recovery.\
# find_coefficients_Jakob2019" definition.
L_middle = lightness_steps // 3
coefficients_middle = optimize(
np.hstack([ijk, L_middle]), zeros(3), chroma
)
# Down the lightness scale.
coefficients_0 = coefficients_middle
for L in reversed(range(L_middle)):
coefficients_0 = optimize(
np.hstack([ijk, L]), coefficients_0, chroma
)
# Up the lightness scale.
coefficients_0 = coefficients_middle
for L in range(L_middle + 1, lightness_steps):
coefficients_0 = optimize(
np.hstack([ijk, L]), coefficients_0, chroma
)
self._size = size
self._create_interpolator()
def RGB_to_coefficients(self, RGB: ArrayLike) -> NDArrayFloat:
"""
Look up a given *RGB* colourspace array and return corresponding
coefficients. Interpolation is used for colours not on the table grid.
Parameters
----------
RGB
*RGB* colourspace array.
Returns
-------
:class:`numpy.ndarray`
Corresponding coefficients that can be passed to
:func:`colour.recovery.jakob2019.sd_Jakob2019` to obtain a spectral
distribution.
Raises
------
RuntimeError
If the pre-computed lookup table has not been generated or read.
Examples
--------
>>> from colour import MSDS_CMFS, SDS_ILLUMINANTS
>>> from colour.models import RGB_COLOURSPACE_sRGB
>>> cmfs = (
... MSDS_CMFS["CIE 1931 2 Degree Standard Observer"]
... .copy()
... .align(SpectralShape(360, 780, 10))
... )
>>> illuminant = SDS_ILLUMINANTS["D65"].copy().align(cmfs.shape)
>>> LUT = LUT3D_Jakob2019()
>>> LUT.generate(RGB_COLOURSPACE_sRGB, cmfs, illuminant, 3, lambda x: x)
>>> RGB = np.array([0.70573936, 0.19248266, 0.22354169])
>>> LUT.RGB_to_coefficients(RGB) # doctest: +ELLIPSIS
array([ 1.5013448...e-04, -1.4679754...e-01, 3.4020219...e+01])
"""
if len(self._interpolator.grid) != 0:
RGB = as_float_array(RGB)
value_max = np.max(RGB, axis=-1)
chroma = RGB / (value_max[..., None] + 1e-10)
i_m = np.argmax(RGB, axis=-1)
i_1 = index_along_last_axis(RGB, i_m)
i_2 = index_along_last_axis(chroma, (i_m + 2) % 3)
i_3 = index_along_last_axis(chroma, (i_m + 1) % 3)
indexes = np.stack([i_m, i_1, i_2, i_3], axis=-1)
return self._interpolator(indexes).squeeze()
else:
raise RuntimeError(
"The pre-computed lookup table has not been read or generated!"
)
def RGB_to_sd(
self, RGB: ArrayLike, shape: SpectralShape = SPECTRAL_SHAPE_JAKOB2019
) -> SpectralDistribution:
"""
Look up a given *RGB* colourspace array and return the corresponding
spectral distribution.
Parameters
----------
RGB
*RGB* colourspace array.
shape
Shape used by the spectral distribution.
Returns
-------
:class:`colour.SpectralDistribution`
Spectral distribution corresponding with the RGB* colourspace
array.
Examples
--------
>>> from colour import MSDS_CMFS, SDS_ILLUMINANTS
>>> from colour.models import RGB_COLOURSPACE_sRGB
>>> from colour.utilities import numpy_print_options
>>> cmfs = (
... MSDS_CMFS["CIE 1931 2 Degree Standard Observer"]
... .copy()
... .align(SpectralShape(360, 780, 10))
... )
>>> illuminant = SDS_ILLUMINANTS["D65"].copy().align(cmfs.shape)
>>> LUT = LUT3D_Jakob2019()
>>> LUT.generate(RGB_COLOURSPACE_sRGB, cmfs, illuminant, 3, lambda x: x)
>>> RGB = np.array([0.70573936, 0.19248266, 0.22354169])
>>> with numpy_print_options(suppress=True):
... LUT.RGB_to_sd(RGB, cmfs.shape) # doctest: +ELLIPSIS
SpectralDistribution([[ 360. , 0.7666803...],
[ 370. , 0.6251547...],
[ 380. , 0.4584310...],
[ 390. , 0.3161633...],
[ 400. , 0.2196155...],
[ 410. , 0.1596575...],
[ 420. , 0.1225525...],
[ 430. , 0.0989784...],
[ 440. , 0.0835782...],
[ 450. , 0.0733535...],
[ 460. , 0.0666049...],
[ 470. , 0.0623569...],
[ 480. , 0.06006 ...],
[ 490. , 0.0594383...],
[ 500. , 0.0604201...],
[ 510. , 0.0631195...],
[ 520. , 0.0678648...],
[ 530. , 0.0752834...],
[ 540. , 0.0864790...],
[ 550. , 0.1033773...],
[ 560. , 0.1293883...],
[ 570. , 0.1706018...],
[ 580. , 0.2374178...],
[ 590. , 0.3439472...],
[ 600. , 0.4950548...],
[ 610. , 0.6604253...],
[ 620. , 0.7914669...],
[ 630. , 0.8738724...],
[ 640. , 0.9213216...],
[ 650. , 0.9486880...],
[ 660. , 0.9650550...],
[ 670. , 0.9752838...],
[ 680. , 0.9819499...],
[ 690. , 0.9864585...],
[ 700. , 0.9896073...],
[ 710. , 0.9918680...],
[ 720. , 0.9935302...],
[ 730. , 0.9947778...],
[ 740. , 0.9957312...],
[ 750. , 0.9964714...],
[ 760. , 0.9970543...],
[ 770. , 0.9975190...],
[ 780. , 0.9978936...]],
SpragueInterpolator,
{},
Extrapolator,
{'method': 'Constant', 'left': None, 'right': None})
"""
sd = sd_Jakob2019(self.RGB_to_coefficients(RGB), shape)
sd.name = f"{RGB!r} (RGB) - Jakob (2019)"
return sd
def read(self, path: str | Path) -> LUT3D_Jakob2019:
"""
Load a lookup table from a *\\*.coeff* file.
Parameters
----------
path
Path to the file.
Returns
-------
LUT3D_Jakob2019
*Jakob and Hanika (2019)* lookup table.
Examples
--------
>>> import os
>>> import colour
>>> from colour import MSDS_CMFS, SDS_ILLUMINANTS
>>> from colour.models import RGB_COLOURSPACE_sRGB
>>> from colour.utilities import numpy_print_options
>>> cmfs = (
... MSDS_CMFS["CIE 1931 2 Degree Standard Observer"]
... .copy()
... .align(SpectralShape(360, 780, 10))
... )
>>> illuminant = SDS_ILLUMINANTS["D65"].copy().align(cmfs.shape)
>>> LUT = LUT3D_Jakob2019()
>>> LUT.generate(RGB_COLOURSPACE_sRGB, cmfs, illuminant, 3, lambda x: x)
>>> path = os.path.join(
... colour.__path__[0],
... "recovery",
... "tests",
... "resources",
... "sRGB_Jakob2019.coeff",
... )
>>> LUT.write(path) # doctest: +SKIP
>>> LUT.read(path) # doctest: +SKIP
"""
path = str(path)
with open(path, "rb") as coeff_file:
if coeff_file.read(4).decode("ISO-8859-1") != "SPEC":
raise ValueError(
"Bad magic number, this is likely not the right file type!"
)
self._size = struct.unpack("i", coeff_file.read(4))[0]
self._lightness_scale = np.fromfile(
coeff_file, count=self._size, dtype=np.float32
)
self._coefficients = np.fromfile(
coeff_file, count=3 * (self._size**3) * 3, dtype=np.float32
)
self._coefficients = np.reshape(
self._coefficients, (3, self._size, self._size, self._size, 3)
)
self._create_interpolator()
return self
def write(self, path: str | Path) -> bool:
"""
Write the lookup table to a *\\*.coeff* file.
Parameters
----------
path
Path to the file.
Returns
-------
:class:`bool`
Definition success.
Examples
--------
>>> import os
>>> import colour
>>> from colour import MSDS_CMFS, SDS_ILLUMINANTS
>>> from colour.models import RGB_COLOURSPACE_sRGB
>>> from colour.utilities import numpy_print_options
>>> cmfs = (
... MSDS_CMFS["CIE 1931 2 Degree Standard Observer"]
... .copy()
... .align(SpectralShape(360, 780, 10))
... )
>>> illuminant = SDS_ILLUMINANTS["D65"].copy().align(cmfs.shape)
>>> LUT = LUT3D_Jakob2019()
>>> LUT.generate(RGB_COLOURSPACE_sRGB, cmfs, illuminant, 3, lambda x: x)
>>> path = os.path.join(
... colour.__path__[0],
... "recovery",
... "tests",
... "resources",
... "sRGB_Jakob2019.coeff",
... )
>>> LUT.write(path) # doctest: +SKIP
>>> LUT.read(path) # doctest: +SKIP
"""
path = str(path)
with open(path, "wb") as coeff_file:
coeff_file.write(b"SPEC")
coeff_file.write(struct.pack("i", self._coefficients.shape[1]))
np.float32(self._lightness_scale).tofile(coeff_file)
np.float32(self._coefficients).tofile(coeff_file)
return True