diff --git a/colour/__init__.py b/colour/__init__.py index b3325c5c03..02a2a74bfa 100644 --- a/colour/__init__.py +++ b/colour/__init__.py @@ -438,6 +438,7 @@ is_within_mesh_volume, is_within_pointer_gamut, is_within_visible_spectrum, + macadam_limits, ) from .graph import describe_conversion_path, convert @@ -853,6 +854,7 @@ def __getattr__(self, attribute) -> Any: "RGB_colourspace_volume_MonteCarlo", "RGB_colourspace_volume_coverage_MonteCarlo", "is_within_macadam_limits", + "macadam_limits", "is_within_mesh_volume", "is_within_pointer_gamut", "is_within_visible_spectrum", diff --git a/colour/volume/__init__.py b/colour/volume/__init__.py index 9871de1bec..aa5ac5c7d3 100644 --- a/colour/volume/__init__.py +++ b/colour/volume/__init__.py @@ -1,6 +1,7 @@ from .datasets import * # noqa from . import datasets from .macadam_limits import is_within_macadam_limits +from .macadam_limits import macadam_limits from .mesh import is_within_mesh_volume from .pointer_gamut import is_within_pointer_gamut from .spectrum import ( @@ -22,6 +23,7 @@ __all__ += [ "is_within_macadam_limits", ] +__all__ += ["macadam_limits"] __all__ += [ "is_within_mesh_volume", ] diff --git a/colour/volume/macadam_limits.py b/colour/volume/macadam_limits.py index c7a8d7bf4b..f2d6df93a8 100644 --- a/colour/volume/macadam_limits.py +++ b/colour/volume/macadam_limits.py @@ -1,7 +1,6 @@ """ Optimal Colour Stimuli - MacAdam Limits ======================================= - Defines the objects related to *Optimal Colour Stimuli* computations. """ @@ -9,7 +8,6 @@ import numpy as np from scipy.spatial import Delaunay - from colour.hints import ( ArrayLike, Dict, @@ -19,11 +17,17 @@ Optional, Union, ) +from colour.colorimetry import ( + MSDS_CMFS, + reshape_sd, + SpectralShape, + SDS_ILLUMINANTS, +) from colour.models import xyY_to_XYZ from colour.volume import OPTIMAL_COLOUR_STIMULI_ILLUMINANTS -from colour.utilities import CACHE_REGISTRY, validate_method +from colour.utilities import CACHE_REGISTRY, tsplit, zeros, validate_method -__author__ = "Colour Developers" +__author__ = "Colour Developers", "Christian Greim" __copyright__ = "Copyright 2013 Colour Developers" __license__ = "New BSD License - https://opensource.org/licenses/BSD-3-Clause" __maintainer__ = "Colour Developers" @@ -142,3 +146,202 @@ def is_within_macadam_limits( simplex = np.where(simplex >= 0, True, False) return simplex + + +def macadam_limits( + luminance: Floating = 0.5, + illuminant=SDS_ILLUMINANTS["E"], + spectral_range=SpectralShape(360, 830, 1), + cmfs=MSDS_CMFS["CIE 1931 2 Degree Standard Observer"], +) -> NDArray: + """ + Return an array of CIE -X,Y,Z - Triples containing colour-coordinates + of the MacAdam-limit for the defined luminance for every + whavelength defined in spectral_range. + Target ist a fast running codey, by + not simply testing the possible optimums step by step but + more effectively targeted by steps of power of two. The wavelengths + left and right of a rough optimum are fitted by a rule of proportion, + so that the wished brightness will be reached exactly. + + Parameters + ---------- + luminance + set the wanted luminance + has to be between 0 and 1 + illuminant + Illuminant spectral distribution, default to *CIE Illuminant E* + spectral_range + SpectralShape according to colour.SpectralShape + cmfs + Standard observer colour matching functions, default to the + *CIE 1931 2 Degree Standard Observer*. + + Returns + ------- + :class:`numpy.ndarray` + an array of CIE -X,Y,Z - Triples containing colour-coordinates + of the MacAdam-limit for the definde luminance for every + whavelength defined in spectral_range + array([[ 3.83917134e-01, 5.00000000e-01, 3.55171511e-01], + [ 3.56913361e-01, 5.00000000e-01, 3.55215349e-01], + [ 3.32781985e-01, 5.00000000e-01, 3.55249953e-01], + ... + [ 4.44310989e-01, 5.00000000e-01, 3.55056751e-01], + [ 4.13165551e-01, 5.00000000e-01, 3.55118668e-01]]) + + References + ---------- + - cite: Wyszecki, G., & Stiles, W. S. (2000). + In Color Science: Concepts and Methods, + Quantitative Data and Formulae (pp. 181–184). Wiley. + ISBN:978-0-471-39918-6 + - cite: Francisco Martínez-Verdú, Esther Perales, + Elisabet Chorro, Dolores de Fez, + Valentín Viqueira, and Eduardo Gilabert, "Computation and + visualization of the MacAdam limits + for any lightness, hue angle, and light source," J. + Opt. Soc. Am. A 24, 1501-1515 (2007) + - cite: Kenichiro Masaoka. In OPTICS LETTERS, June 15, 2010 + / Vol. 35, No. 1 (pp. 2031 - 2033) + + Examples + -------- + from matplotlib import pyplot as plt + import numpy as np + import math + fig = plt.figure(figsize=(7,7)) + ax = fig.add_axes([0,0,1,1]) + illuminant = colour.SDS_ILLUMINANTS['D65'] + def plot_Narrowband_Spectra (Yxy_Narrowband_Spectra): + FirstColumn = 0 + SecondColumn = 1 + x = Yxy_Narrowband_Spectra[...,FirstColumn] + y = Yxy_Narrowband_Spectra[...,SecondColumn] + ax.plot(x,y,'orange',label='Spectrum Loci') + x = [Yxy_Narrowband_Spectra[-1][FirstColumn], + Yxy_Narrowband_Spectra[0][FirstColumn]] + y = [Yxy_Narrowband_Spectra[-1][SecondColumn], + Yxy_Narrowband_Spectra[0][SecondColumn]] + ax.plot(x,y,'purple',label='Purple Boundary') + return() + for n in range(1, 20): + Yxy_Narrowband_Spectra = colour.XYZ_to_xy( + colour.macadam_limits(n/20, illuminant)) + plot_Narrowband_Spectra (Yxy_Narrowband_Spectra) + plt.show() + """ + + target_bright = luminance + if target_bright > 1 or target_bright < 0: + raise TypeError( + "brightness of function macadam_limits( )" + "has to be between 0 and 1" + ) + # workarround because illuminant and cmfs are rounded + # in a different way. + illuminant = reshape_sd(illuminant, cmfs.shape) + cmfs = reshape_sd(cmfs, spectral_range) + illuminant = reshape_sd(illuminant, spectral_range) + + # The cmfs are convolved with the given illuminant + X_illuminated, Y_illuminated, Z_illuminated = ( + tsplit(cmfs.values) * illuminant.values + ) + # Generate empty output-array + out_limits = np.zeros_like(cmfs.values) + # For examle a SpectralShape(360, 830, 1) has 471 entries + opti_colour = np.zeros_like(Y_illuminated) + # The array of optimal colours has the same dimensions + # like Y_illuminated, in our example: 471 + colour_range = illuminant.values.shape[0] + a = np.arange(12) + a = np.ceil(colour_range / 2 ** (a + 1)).astype(int) + step_sizes = np.append(np.flip(np.unique(a)), 1) + middle_opti_colour = step_sizes[0] - 1 + # in our example: 235 + width = step_sizes[1] - 1 + # in our example: 117 + step_sizes = np.delete(step_sizes, [0, 1]) + # in our example: np.array([59 30 15 8 4 2 1]) + # The first optimum color has its center initially at zero + maximum_brightness = np.sum(Y_illuminated) + + def optimum_colour(width, center): + opti_colour = zeros(colour_range) + # creates array of 471 zeros and ones which represents optimum-colours + # All values of the opti_colour-array are initially set to zero + half_width = width + center_opti_colour = center + opti_colour[ + middle_opti_colour + - half_width : middle_opti_colour + + 1 + + half_width + ] = 1 + # we start the construction of the optimum color + # at the center of the opti_colour-array + opti_colour = np.roll( + opti_colour, center_opti_colour - middle_opti_colour + ) + # the optimum colour is rolled to the right wavelength + return opti_colour + + def bright_opti_colour(width, center, lightsource): + brightness = ( + np.sum(optimum_colour(width, center) * lightsource) + / maximum_brightness + ) + return brightness + + # here we do some kind of Newton's Method to aproximate the + # wandted illuminance at the whavelengt. + # therefore the numbers 127, 64, 32 and so on + # step_size is in this case np.array([59 30 15 8 4 2 1]) + for wavelength in range(0, colour_range): + for n in step_sizes: + brightness = bright_opti_colour(width, wavelength, Y_illuminated) + if brightness > target_bright or width >= middle_opti_colour: + width -= n + else: + width += n + + brightness = bright_opti_colour(width, wavelength, Y_illuminated) + if brightness < target_bright: + width += 1 + + rough_optimum = optimum_colour(width, wavelength) + brightness = np.sum(rough_optimum * Y_illuminated) / maximum_brightness + + # in the following, the both borders of the found rough_optimum + # are reduced to get more exact results + bright_difference = (brightness - target_bright) * maximum_brightness + # discrimination for single-wavelength-spectra + if width > 0: + opti_colour = zeros(colour_range) + opti_colour[ + middle_opti_colour - width : middle_opti_colour + width + 1 + ] = 1 + # instead rolling forward opti_colour, light is rolled backward + rolled_light = np.roll( + Y_illuminated, middle_opti_colour - wavelength + ) + opti_colour_light = opti_colour * rolled_light + left_opti = opti_colour_light[middle_opti_colour - width] + right_opti = opti_colour_light[middle_opti_colour + width] + interpolation = 1 - (bright_difference / (left_opti + right_opti)) + opti_colour[middle_opti_colour - width] = interpolation + opti_colour[middle_opti_colour + width] = interpolation + # opti_colour is rolled to right position + final_optimum = np.roll( + opti_colour, wavelength - middle_opti_colour + ) + else: + final_optimum = rough_optimum / brightness * target_bright + + out_X = np.sum(final_optimum * X_illuminated) / maximum_brightness + out_Y = target_bright + out_Z = np.sum(final_optimum * Z_illuminated) / maximum_brightness + triple = np.array([out_X, out_Y, out_Z]) + out_limits[wavelength] = triple + return out_limits