-
-
Notifications
You must be signed in to change notification settings - Fork 266
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
PR: Implement support for "MacAdam Limits" computation. #768
base: develop
Are you sure you want to change the base?
Changes from all commits
16ae787
47883a8
9e27d71
f204568
99872c3
308712e
6d34507
5b8593a
afc9946
170eb03
5cbd91d
d6fe7ee
b03eaf6
f3bc4fb
7296c92
f5ee2a4
388c0f6
91c45f7
6eb76d5
41a277b
977be68
9dc7145
6f306dc
ac12a6b
75ade26
463a1b2
92a1cd1
78141ce
a190c00
c67984d
a603c18
0659c32
f91801f
95a446e
38b727a
fa6485b
c1cd25a
9f28353
e58c32b
fd3a52f
5d817da
c87523b
31f2dbd
5ed84a0
cbdb8b4
d1eca3b
b5f06c8
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change | ||
---|---|---|---|---|
@@ -1,15 +1,13 @@ | ||||
""" | ||||
Optimal Colour Stimuli - MacAdam Limits | ||||
======================================= | ||||
|
||||
Defines the objects related to *Optimal Colour Stimuli* computations. | ||||
""" | ||||
|
||||
from __future__ import annotations | ||||
|
||||
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 | ||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Looking at this implementation, I'm left wondering if we should not simply output the XYZ values from this definition which does a similar work but using Trimesh and geometrical intersection of the Rösch-MacAdam colour solid with a plane: colour/colour/plotting/section.py Line 479 in a412b67
https://colour.readthedocs.io/en/develop/_images/Plotting_Plot_Visible_Spectrum_Section.png There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I was thinking the same thing after looking at the spectrum.py code in #994 There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. At a first glimps it seems to be a good idea. Some reasons for my code:
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Y
is luminance, not brightness.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Changed in actual version.