From 1934d1ce12200e26e9e76da9cd77b4593bc0a0c2 Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Sat, 20 Jun 2020 22:48:47 +0200 Subject: Rewrite the error function w/ analytic derivatives This greatly speeds up the code and improves convergence. Something's slightly off with the code, though, and the output colors are slightly skewed. --- gsoc_common.py | 88 ++++++++++++++++++++++++++++++++++++++++++++++++--------- test_diff.py | 51 +++++++++++++++++++++++++++++++++ test_example.py | 3 +- 3 files changed, 127 insertions(+), 15 deletions(-) create mode 100644 test_diff.py diff --git a/gsoc_common.py b/gsoc_common.py index b028e45..5d4aea0 100644 --- a/gsoc_common.py +++ b/gsoc_common.py @@ -9,8 +9,8 @@ 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) +wvl = SpectralShape(360, 830, 5) +wvlp = (wvl.range() - 360) / (830 - 360) # This is the model of spectral reflectivity described in the article. def model(wvlp, ccp): @@ -21,27 +21,87 @@ def model(wvlp, ccp): 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") + grid = wvlp if primed else wvl.range() + return SpectralDistribution(model(grid, ccp), wvl.range(), 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) +# This function also calculates the first derivatives with respect to c's. +def error_function(c, target, wvl, cmfs, illuminant, illuminant_XYZ): + U = c[0] * wvl**2 + c[1] * wvl + c[2] + t1 = np.sqrt(1 + U**2) + R = 1 / 2 + U / (2 * t1) + t2 = 1 / (2 * t1) - U**2 / (2 * t1**3) + dR_dc = [wvl**2 * t2, wvl * t2, t2] + E = illuminant.values * R / 100 + dE_dc = illuminant.values * dR_dc / 100 + + XYZ = np.empty(3) + dXYZ_dc = np.empty((3, 3)) + + dlambda = cmfs.wavelengths[1] - cmfs.wavelengths[0] + for i in range(3): + XYZ[i] = np.dot(E, cmfs.values[:, i]) * dlambda + for j in range(3): + dXYZ_dc[i, j] = np.dot(dE_dc[j], cmfs.values[:, i]) * dlambda + + # FIXME: this isn't the full CIE 1976 lightness function + f = (XYZ / illuminant_XYZ)**(1/3) + + # FIXME: this can be vectorized + df_dc = np.empty((3, 3)) + for i in range(3): + for j in range(3): + df_dc[i, j] = 1 / (3 * illuminant_XYZ[i]**(1/3) + * XYZ[i]**(2/3)) * dXYZ_dc[i, j] + + Lab = np.array([ + 116 * f[1] - 16, + 500 * (f[0] - f[1]), + 200 * (f[1] - f[2]) + ]) + + dLab_dc = np.array([ + 116 * df_dc[1], + 500 * (df_dc[0] - df_dc[1]), + 200 * (df_dc[1] - df_dc[2]) + ]) + + error = np.sqrt(np.sum((Lab - target)**2)) + + derror_dc = np.zeros(3) + for i in range(3): + for j in range(3): + derror_dc[i] += dLab_dc[j, i] * (Lab[j] - target[j]) + derror_dc[i] /= error + + #print("c=[%.3g, %.3g, %.3g], XYZ=[%.3g, %.3g, %.3g], Lab=[%.3g, %.3g, %.3g], error=%g" % (*c, *XYZ, *Lab, error)) + + return error, derror_dc # 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): + ill_sd = ill_sd.align(wvl) + ill_XYZ = sd_to_XYZ(ill_sd) / 100 + def do_optimize(XYZ, ccp0): Lab = XYZ_to_Lab(XYZ, ill_xy) + + def fun(*args): + error, derror_dc = error_function(*args) + return error + + def jac(*args): + error, derror_dc = error_function(*args) + return derror_dc + + cmfs = STANDARD_OBSERVER_CMFS["CIE 1931 2 Degree Standard Observer"].align(wvl) + opt = optimize.minimize( - error_function, ccp0, (Lab, ill_sd, ill_xy), - method="Nelder-Mead", options={"disp": verbose} + fun, ccp0, (Lab, wvlp, cmfs, ill_sd, ill_XYZ), jac=jac, + method="L-BFGS-B", options={"disp": verbose} ) if verbose: print(opt) @@ -61,7 +121,7 @@ def jakob_hanika(target_XYZ, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True, try_h good_ccp = (2.1276356, -1.07293026, -0.29583292) # FIXME: valid only for D65 divisions = 3 - while divisions < 10: + while divisions < 30: if verbose: print("Trying with %d divisions" % divisions) @@ -76,7 +136,7 @@ def jakob_hanika(target_XYZ, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True, try_h 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 opt.fun > 1e-3: if verbose: print("WARNING: intermediate optimization failed") break diff --git a/test_diff.py b/test_diff.py new file mode 100644 index 0000000..b634ca6 --- /dev/null +++ b/test_diff.py @@ -0,0 +1,51 @@ +import numpy as np +from scipy.optimize import minimize +from colour import * +from colour.colorimetry import STANDARD_OBSERVER_CMFS, ILLUMINANT_SDS +from colour.models import eotf_inverse_sRGB, sRGB_to_XYZ +from matplotlib import pyplot as plt +from gsoc_common import plot_comparison, error_function, model_sd, D65_xy + +shape = SpectralShape(360, 830, 1) +cmfs = STANDARD_OBSERVER_CMFS["CIE 1931 2 Degree Standard Observer"].align(shape) + +illuminant = SpectralDistribution(ILLUMINANT_SDS["D65"]).align(shape) +illuminant_XYZ = sd_to_XYZ(illuminant) / 100 +wvl = np.linspace(0, 1, len(shape.range())) + +target = np.array([50, -20, 30]) # Some arbitrary Lab coordinates +xs = np.linspace(-10, 10, 500) +h = xs[1] - xs[0] + +# This test checks if derivatives are calculated correctly by comparing them +# to finite differences. +for c_index in range(3): + errors = np.empty(len(xs)) + derrors = np.empty(len(xs)) + + for i, x in enumerate(xs): + c = np.array([1.0, 1, 1]) + c[c_index] = x + + error, derror_dc = error_function( + c, target, wvl, cmfs, illuminant, illuminant_XYZ + ) + + errors[i] = error + derrors[i] = derror_dc[c_index] + + + plt.subplot(2, 3, 1 + c_index) + plt.xlabel("c%d" % c_index) + plt.ylabel("ΔE") + plt.plot(xs, errors) + + plt.subplot(2, 3, 4 + c_index) + plt.xlabel("c%d" % c_index) + plt.ylabel("dΔE/dc%d" % c_index) + + plt.plot(xs, derrors, "k-") + plt.plot(xs[:-1] + h / 2, np.diff(errors) / h, "r:") + + +plt.show() diff --git a/test_example.py b/test_example.py index 6be982d..80b07c7 100644 --- a/test_example.py +++ b/test_example.py @@ -13,7 +13,8 @@ cc_ref = np.array([ 18.70184886, -13.19804478, 2.12180137]) if __name__ == "__main__": - XYZ = sRGB_to_XYZ(eotf_inverse_sRGB(RGB_ref), D65_xy) + #XYZ = sRGB_to_XYZ(eotf_inverse_sRGB(RGB_ref), D65_xy) + XYZ = np.array([1/3, 1/4, 1/3]) Lab = XYZ_to_Lab(XYZ, D65_xy) print("Target: X=%g, Y=%g, Z=%g, L=%g, a=%g, b=%g" % (*XYZ, *Lab)) -- cgit