From 77221e174cdab64de60450646de77f849b043062 Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Fri, 19 Jun 2020 23:42:16 +0200 Subject: Major cleanup - The interpolator code is no longer here (but in the Colour fork) - Refactoring made the code shorter and cleaner - All tests work properly again --- gsoc_common.py | 120 ++++++++++++++++++++++++++++++++++++++++++++++++ jakob_hanika.py | 93 ------------------------------------- test_bad_convergence.py | 16 ++----- test_coeff.py | 68 --------------------------- test_colorchecker.py | 49 ++++---------------- test_coverage.py | 15 +++--- test_example.py | 24 ++++------ test_integration.py | 11 ++--- test_interpolator.py | 27 +++++++++++ 9 files changed, 178 insertions(+), 245 deletions(-) create mode 100644 gsoc_common.py delete mode 100644 jakob_hanika.py delete mode 100644 test_coeff.py create mode 100644 test_interpolator.py diff --git a/gsoc_common.py b/gsoc_common.py new file mode 100644 index 0000000..b028e45 --- /dev/null +++ b/gsoc_common.py @@ -0,0 +1,120 @@ +import numpy as np, scipy.optimize as optimize +from colour import * +from colour.difference import delta_E_CIE1976 +from colour.colorimetry import * +from colour.plotting import * +from matplotlib import pyplot as plt + +D65_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"]["D65"] +D65 = SpectralDistribution(ILLUMINANT_SDS["D65"]) + +# The same wavelength grid is used throughout +wvl = np.arange(360, 830, 5) +wvlp = (wvl - 360) / (830 - 360) + +# This is the model of spectral reflectivity described in the article. +def model(wvlp, ccp): + yy = ccp[0] * wvlp**2 + ccp[1] * wvlp + ccp[2] + return 1 / 2 + yy / (2 * np.sqrt(1 + yy ** 2)) + +# Create a SpectralDistribution using given coefficients +def model_sd(ccp, primed=True): + # FIXME: don't hardcode the wavelength grid; there should be a way + # of creating a SpectralDistribution from the function alone + grid = wvlp if primed else wvl + return SpectralDistribution(model(grid, ccp), wvl, name="Model") + +# The goal is to minimize the color difference between a given distrbution +# and the one computed from the model above. +def error_function(ccp, target, ill_sd, ill_xy): + # One thing they didn't mention in the text is that their polynomial + # is nondimensionalized (mapped from 360--800 nm to 0--1). + sd = model_sd(ccp) + Lab = XYZ_to_Lab(sd_to_XYZ(sd, illuminant=ill_sd) / 100, ill_xy) + return delta_E_CIE1976(target, Lab) + + + +# Finds the parameters for Jakob and Hanika's model +def jakob_hanika(target_XYZ, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True, try_hard=True): + def do_optimize(XYZ, ccp0): + Lab = XYZ_to_Lab(XYZ, ill_xy) + opt = optimize.minimize( + error_function, ccp0, (Lab, ill_sd, ill_xy), + method="Nelder-Mead", options={"disp": verbose} + ) + if verbose: + print(opt) + return opt + + # A special case that's hard to solve numerically + if np.allclose(target_XYZ, [0, 0, 0]): + return np.array([0, 0, -1e+9]), 0 # FIXME: dtype? + + if verbose: + print("Trying the target directly, XYZ=%s" % target_XYZ) + opt = do_optimize(target_XYZ, ccp0) + if opt.fun < 0.1 or not try_hard: + return opt.x, opt.fun + + good_XYZ = (1/3, 1/3, 1/3) + good_ccp = (2.1276356, -1.07293026, -0.29583292) # FIXME: valid only for D65 + + divisions = 3 + while divisions < 10: + if verbose: + print("Trying with %d divisions" % divisions) + + keep_divisions = False + ref_XYZ = good_XYZ + ref_ccp = good_ccp + + ccp0 = ref_ccp + for i in range(1, divisions): + intermediate_XYZ = ref_XYZ + (target_XYZ - ref_XYZ) * i / (divisions - 1) + if verbose: + print("Intermediate step %d/%d, XYZ=%s with ccp0=%s" % + (i + 1, divisions, intermediate_XYZ, ccp0)) + opt = do_optimize(intermediate_XYZ, ccp0) + if opt.fun > 0.1: + if verbose: + print("WARNING: intermediate optimization failed") + break + else: + good_XYZ = intermediate_XYZ + good_ccp = opt.x + keep_divisions = True + + ccp0 = opt.x + else: + return opt.x, opt.fun + + if not keep_divisions: + divisions += 2 + + raise Exception("optimization failed for target_XYZ=%s, ccp0=%s" \ + % (target_XYZ, ccp0)) + +# Makes a comparison plot with SDs and swatches +def plot_comparison(target, matched_sd, label, error, ill_sd, show=True): + if type(target) is SpectralDistribution: + target_XYZ = sd_to_XYZ(target, illuminant=ill_sd) / 100 + else: + target_XYZ = target + target_RGB = np.clip(XYZ_to_sRGB(target_XYZ), 0, 1) + target_swatch = ColourSwatch(label, target_RGB) + matched_XYZ = sd_to_XYZ(matched_sd, illuminant=ill_sd) / 100 + matched_RGB = np.clip(XYZ_to_sRGB(matched_XYZ), 0, 1) + matched_swatch = ColourSwatch("Model", matched_RGB) + + axes = plt.subplot(2, 1, 1) + plt.title(label) + if type(target) is SpectralDistribution: + plot_multi_sds([target, matched_sd], axes=axes, standalone=False) + else: + plot_single_sd(matched_sd, axes=axes, standalone=False) + + axes = plt.subplot(2, 1, 2) + plt.title("ΔE = %g" % error) + plot_multi_colour_swatches([target_swatch, matched_swatch], + standalone=show, axes=axes) diff --git a/jakob_hanika.py b/jakob_hanika.py deleted file mode 100644 index 22b87fc..0000000 --- a/jakob_hanika.py +++ /dev/null @@ -1,93 +0,0 @@ -import numpy as np, scipy.optimize as optimize -from colour import * -from colour.difference import delta_E_CIE1976 -from colour.colorimetry import * - - - -# The same wavelength grid is used throughout -wvl = np.arange(360, 830, 5) -wvlp = (wvl - 360) / (830 - 360) - -# This is the model of spectral reflectivity described in the article. -def model(wvlp, ccp): - yy = ccp[0] * wvlp**2 + ccp[1] * wvlp + ccp[2] - return 1 / 2 + yy / (2 * np.sqrt(1 + yy ** 2)) - -# Create a SpectralDistribution using given coefficients -def model_sd(ccp, primed=True): - # FIXME: don't hardcode the wavelength grid; there should be a way - # of creating a SpectralDistribution from the function alone - grid = wvlp if primed else wvl - return SpectralDistribution(model(grid, ccp), wvl, name="Model") - -# The goal is to minimize the color difference between a given distrbution -# and the one computed from the model above. -def error_function(ccp, target, ill_sd, ill_xy): - # One thing they didn't mention in the text is that their polynomial - # is nondimensionalized (mapped from 360--800 nm to 0--1). - sd = model_sd(ccp) - Lab = XYZ_to_Lab(sd_to_XYZ(sd, illuminant=ill_sd) / 100, ill_xy) - return delta_E_CIE1976(target, Lab) - - - -# Finds the parameters for Jakob and Hanika's model -def jakob_hanika(target_XYZ, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True, try_hard=True): - def do_optimize(XYZ, ccp0): - Lab = XYZ_to_Lab(XYZ, ill_xy) - opt = optimize.minimize( - error_function, ccp0, (Lab, ill_sd, ill_xy), - method="Nelder-Mead", options={"disp": verbose} - ) - if verbose: - print(opt) - return opt - - # A special case that's hard to solve numerically - if np.allclose(target_XYZ, [0, 0, 0]): - return np.array([0, 0, -1e+9]), 0 # FIXME: dtype? - - if verbose: - print("Trying the target directly, XYZ=%s" % target_XYZ) - opt = do_optimize(target_XYZ, ccp0) - if opt.fun < 0.1 or not try_hard: - return opt.x, opt.fun - - good_XYZ = (1/3, 1/3, 1/3) - good_ccp = (2.1276356, -1.07293026, -0.29583292) # FIXME: valid only for D65 - - divisions = 3 - while divisions < 10: - if verbose: - print("Trying with %d divisions" % divisions) - - keep_divisions = False - ref_XYZ = good_XYZ - ref_ccp = good_ccp - - ccp0 = ref_ccp - for i in range(1, divisions): - intermediate_XYZ = ref_XYZ + (target_XYZ - ref_XYZ) * i / (divisions - 1) - if verbose: - print("Intermediate step %d/%d, XYZ=%s with ccp0=%s" % - (i + 1, divisions, intermediate_XYZ, ccp0)) - opt = do_optimize(intermediate_XYZ, ccp0) - if opt.fun > 0.1: - if verbose: - print("WARNING: intermediate optimization failed") - break - else: - good_XYZ = intermediate_XYZ - good_ccp = opt.x - keep_divisions = True - - ccp0 = opt.x - else: - return opt.x, opt.fun - - if not keep_divisions: - divisions += 2 - - raise Exception("optimization failed for target_XYZ=%s, ccp0=%s" \ - % (target_XYZ, ccp0)) diff --git a/test_bad_convergence.py b/test_bad_convergence.py index dadb21e..a47af5a 100644 --- a/test_bad_convergence.py +++ b/test_bad_convergence.py @@ -1,15 +1,6 @@ import numpy as np from colour import * -from jakob_hanika import jakob_hanika, model_sd -from test_colorchecker import plot_comparison - - - -illuminant = "D65" -ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant] -ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant]) - - +from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison # A set of inputs that make the Nelder-Mead solver work really hard. # This might not be interesting in future versions of jakob_hanika.py, once @@ -19,6 +10,5 @@ if __name__ == "__main__": XYZ = np.array([0.00263241, 0.00174905, 0.00092602]) ccp0 = np.array([5.11263599, -2.65344447, -2.37301856]) - ccp, error = jakob_hanika(XYZ, ill_sd, ill_xy, ccp0) - if error > 0.1: - plot_comparison(XYZ, model_sd(ccp), "Target", error, ill_sd) + ccp, error = jakob_hanika(XYZ, D65, D65_xy, ccp0) + plot_comparison(XYZ, model_sd(ccp), "Target", error, D65) diff --git a/test_coeff.py b/test_coeff.py deleted file mode 100644 index c2a8970..0000000 --- a/test_coeff.py +++ /dev/null @@ -1,68 +0,0 @@ -import numpy as np, struct -from scipy.interpolate import RegularGridInterpolator -from colour import * -from colour.plotting import * -from colour.models import eotf_inverse_sRGB -from colour.difference import delta_E_CIE1976 -from jakob_hanika import jakob_hanika, model_sd -from test_colorchecker import plot_comparison - - - -illuminant = "D65" -ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant] -ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant]) - - - -class JakobHanikaInterpolator: - def __init__(self): - pass - - def read_file(self, path): - with open(path, "rb") as fd: - if fd.read(4).decode("ISO-8859-1") != "SPEC": - raise ValueError("bad magic number") - - self.res = struct.unpack("i", fd.read(4))[0] - self.scale = np.fromfile(fd, count=self.res, dtype=np.float32) # FIXME: default dtype - coeffs = np.fromfile(fd, count=3*self.res**3*3, dtype=np.float32) # FIXME: default dtype - coeffs = coeffs.reshape(3, self.res, self.res, self.res, 3) - - t = np.linspace(0, 1, self.res) - axes = [[0, 1, 2], self.scale, t, t] - self.cubes = RegularGridInterpolator(axes, coeffs[:, :, :, :, :], bounds_error=False) - - def __call__(self, RGB): - RGB = np.asarray(RGB, dtype=np.float) # FIXME: default dtype - vmax = np.max(RGB, axis=-1) - imax = np.argmax(RGB, axis=-1) - chroma = RGB / (np.expand_dims(vmax, -1) + 1e-10) # Avoid division by zero - vmax = np.max(RGB, axis=-1) - v2 = np.take_along_axis(chroma, np.expand_dims((imax + 2) % 3, axis=-1), axis=-1).squeeze(axis=-1) - v3 = np.take_along_axis(chroma, np.expand_dims((imax + 1) % 3, axis=-1), axis=-1).squeeze(axis=-1) - coords = np.stack([imax, vmax, v2, v3], axis=-1) - return self.cubes(coords).squeeze() - - - -if __name__ == "__main__": - jh = JakobHanikaInterpolator() - jh.read_file("data/srgb.coeff") - - RGBs = np.random.random((7, 6, 5, 4, 3)) - ccs = jh(RGBs) - - RGB = eotf_inverse_sRGB(RGBs[0, 0, 0, 0, :]) - XYZ = sRGB_to_XYZ(RGB) - Lab = XYZ_to_Lab(XYZ, ill_xy) - cc = ccs[0, 0, 0, 0, :] - - matched_sd = model_sd(cc, primed=False) - matched_XYZ = sd_to_XYZ(matched_sd, illuminant=ill_sd) - matched_Lab = XYZ_to_Lab(matched_XYZ, ill_xy) - error = delta_E_CIE1976(Lab, matched_Lab) - - print(matched_sd) - - plot_comparison(XYZ, matched_sd, "Model", error, ill_sd, ill_xy) diff --git a/test_colorchecker.py b/test_colorchecker.py index b0c9688..fc54d84 100644 --- a/test_colorchecker.py +++ b/test_colorchecker.py @@ -1,48 +1,15 @@ import numpy as np from colour import * -from colour.plotting import * -from matplotlib import pyplot as plt -from jakob_hanika import jakob_hanika, model_sd - - - -illuminant = "D65" -ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant] -ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant]) - - - -# Makes a comparison plot with SDs and swatches -def plot_comparison(target, matched_sd, label, error, ill_sd, show=True): - if type(target) is SpectralDistribution: - target_XYZ = sd_to_XYZ(target, illuminant=ill_sd) / 100 - else: - target_XYZ = target - target_RGB = np.clip(XYZ_to_sRGB(target_XYZ), 0, 1) - target_swatch = ColourSwatch(label, target_RGB) - matched_XYZ = sd_to_XYZ(matched_sd, illuminant=ill_sd) / 100 - matched_RGB = np.clip(XYZ_to_sRGB(matched_XYZ), 0, 1) - matched_swatch = ColourSwatch("Model", matched_RGB) - - axes = plt.subplot(2, 1, 1) - plt.title(label) - if type(target) is SpectralDistribution: - plot_multi_sds([target, matched_sd], axes=axes, standalone=False) - else: - plot_single_sd(matched_sd, axes=axes, standalone=False) - - axes = plt.subplot(2, 1, 2) - plt.title("ΔE = %g" % error) - plot_multi_colour_swatches([target_swatch, matched_swatch], - standalone=show, axes=axes) +from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison +# This demo goes through SDs in a color checker if __name__ == "__main__": - # This demo goes through SDs in a color checker - for name, sd in COLOURCHECKERS_SDS['ColorChecker N Ohta'].items(): - XYZ = sd_to_XYZ(sd, illuminant=ill_sd) / 100 + for name, sd in COLOURCHECKER_SDS['ColorChecker N Ohta'].items(): + XYZ = sd_to_XYZ(sd, illuminant=D65) / 100 - print("Color checker: The target is '%s' with X=%g, Y=%g, Z=%g" % (name, *XYZ)) - ccp, error = jakob_hanika(XYZ, ill_sd, ill_xy) + print("Color checker: The target is '%s' with X=%g, Y=%g, Z=%g" + % (name, *XYZ)) + ccp, error = jakob_hanika(XYZ, D65, D65_xy) matched_sd = model_sd(ccp) - plot_comparison(sd, matched_sd, name, error, ill_sd) + plot_comparison(sd, matched_sd, name, error, D65) diff --git a/test_coverage.py b/test_coverage.py index e73e015..90e3c79 100644 --- a/test_coverage.py +++ b/test_coverage.py @@ -1,14 +1,10 @@ import numpy as np, os, multiprocessing from colour import * from colour.models import eotf_inverse_sRGB -from jakob_hanika import jakob_hanika, model_sd +from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison -illuminant = "D65" -ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant] -ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant]) # Note: changed in dev (to ILLUMINANT_SD) - # Resolution of the discretization cubes # Both should be 64, but this code is a bit too slow at the moment CHROMA_STEPS = 8 @@ -19,10 +15,10 @@ LIGHTNESS_STEPS = 8 # Solve for a specific RGB color def optimize_RGB(linear_RGB, ccp0): RGB = eotf_inverse_sRGB(linear_RGB) - target = sRGB_to_XYZ(RGB, ill_xy) - ccp, error = jakob_hanika(target, ill_sd, ill_xy, ccp0, verbose=False) + target = sRGB_to_XYZ(RGB, D65_xy) + ccp, error = jakob_hanika(target, D65, D65_xy, ccp0, verbose=False) if error > 0.1: - print("WARNING: convergence failed with L=%s, starting from %s" \ + print("WARNING: convergence failed with L=%s, starting from %s" % (linear_RGB, ccp0)) return ccp, error @@ -90,7 +86,8 @@ if __name__ == "__main__": dirname = "out/R%.3fG%.3fB%.3f" % (*RGB,) total += 1 - if os.path.exists(dirname) and len(os.listdir(dirname)) == LIGHTNESS_STEPS: + if os.path.exists(dirname) and \ + len(os.listdir(dirname)) == LIGHTNESS_STEPS: done += 1 continue diff --git a/test_example.py b/test_example.py index f2e5345..6be982d 100644 --- a/test_example.py +++ b/test_example.py @@ -1,15 +1,11 @@ import numpy as np from colour import * from colour.difference import delta_E_CIE1976 -from jakob_hanika import jakob_hanika, model_sd -from test_colorchecker import plot_comparison +from colour.models import eotf_inverse_sRGB +from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison -illuminant = "D65" -ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant] -ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant]) - # These numbers are taken from Jakob and Hanika's Jupyter notebook. RGB_ref = np.array([0.79264853, 0.4, 0.63703843]) # *linear* sRGB cc_ref = np.array([ 18.70184886, -13.19804478, 2.12180137]) @@ -17,21 +13,21 @@ cc_ref = np.array([ 18.70184886, -13.19804478, 2.12180137]) if __name__ == "__main__": - XYZ = sRGB_to_XYZ(1.055 * RGB_ref ** (1/2.4) - 0.055, ill_xy) - Lab = XYZ_to_Lab(XYZ, ill_xy) + XYZ = sRGB_to_XYZ(eotf_inverse_sRGB(RGB_ref), D65_xy) + Lab = XYZ_to_Lab(XYZ, D65_xy) print("Target: X=%g, Y=%g, Z=%g, L=%g, a=%g, b=%g" % (*XYZ, *Lab)) - ccp, error = jakob_hanika(XYZ, ill_sd, ill_xy) + ccp, error = jakob_hanika(XYZ, D65, D65_xy) reference_sd = model_sd(cc_ref) - reference_XYZ = sd_to_XYZ(reference_sd, illuminant=ill_sd) - reference_Lab = XYZ_to_Lab(reference_XYZ, ill_xy) + reference_XYZ = sd_to_XYZ(reference_sd, illuminant=D65) + reference_Lab = XYZ_to_Lab(reference_XYZ, D65_xy) found_sd = model_sd(ccp) - found_XYZ = sd_to_XYZ(found_sd, illuminant=ill_sd) - found_Lab = XYZ_to_Lab(found_XYZ, ill_xy) + found_XYZ = sd_to_XYZ(found_sd, illuminant=D65) + found_Lab = XYZ_to_Lab(found_XYZ, D65_xy) print("Our results differ from the reference by ΔE = %g" \ % delta_E_CIE1976(reference_Lab, found_Lab)) - plot_comparison(XYZ, found_sd, "Reference", error, ill_sd) + plot_comparison(XYZ, found_sd, "Reference", error, D65) diff --git a/test_integration.py b/test_integration.py index 35c19ed..e6aa819 100644 --- a/test_integration.py +++ b/test_integration.py @@ -1,19 +1,16 @@ import numpy as np from colour import * +from colour.models import eotf_inverse_sRGB from colour.recovery import Jakob2019Interpolator from colour.plotting import plot_multi_sds +from gsoc_common import D65_xy -illuminant = "D65" -ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant] -ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant]) - - - +# This is a basic test of the feature/jakob2019 branch if __name__ == "__main__": RGB = np.array([0.79264853, 0.4, 0.63703843]) - XYZ = sRGB_to_XYZ(1.055 * RGB ** (1/2.4) - 0.055, ill_xy) + XYZ = sRGB_to_XYZ(eotf_inverse_sRGB(RGB), D65_xy) sd = XYZ_to_sd(XYZ, method="Jakob 2019", verbose=True) interp = Jakob2019Interpolator() diff --git a/test_interpolator.py b/test_interpolator.py new file mode 100644 index 0000000..21a50c5 --- /dev/null +++ b/test_interpolator.py @@ -0,0 +1,27 @@ +import numpy as np, struct +from colour import * +from colour.difference import delta_E_CIE1976 +from colour.models import eotf_inverse_sRGB +from colour.recovery import Jakob2019Interpolator +from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison + +# This script tests if the interpolator correctly handles multi-dimensional +# inputs. +if __name__ == "__main__": + interp = Jakob2019Interpolator() + interp.read_file("data/srgb.coeff") + + RGBs = np.random.random((7, 6, 5, 4, 3)) + ccs = interp.coeffs(RGBs) + + RGB = eotf_inverse_sRGB(RGBs[0, 0, 0, 0, :]) + XYZ = sRGB_to_XYZ(RGB) + Lab = XYZ_to_Lab(XYZ, D65_xy) + cc = ccs[0, 0, 0, 0, :] + + matched_sd = model_sd(cc, primed=False) + matched_XYZ = sd_to_XYZ(matched_sd, illuminant=D65) + matched_Lab = XYZ_to_Lab(matched_XYZ, D65_xy) + error = delta_E_CIE1976(Lab, matched_Lab) + + plot_comparison(XYZ, matched_sd, "Model", error, D65) -- cgit