computeMaxSaturation function

double computeMaxSaturation(double a, double b)

Implementation

double computeMaxSaturation(double a, double b) {
  // Max saturation will be when one of r, g or b goes below zero.

  // Select different coefficients depending on which component goes below zero first
  double k0, k1, k2, k3, k4, wl, wm, ws;

  if (-1.88170328 * a - 0.80936493 * b > 1) {
    // Red component
    k0 = 1.19086277;
    k1 = 1.76576728;
    k2 = 0.59662641;
    k3 = 0.75515197;
    k4 = 0.56771245;
    wl = 4.0767416621;
    wm = -3.3077115913;
    ws = 0.2309699292;
  } else if (1.81444104 * a - 1.19445276 * b > 1) {
    // Green component
    k0 = 0.73956515;
    k1 = -0.45954404;
    k2 = 0.08285427;
    k3 = 0.12541070;
    k4 = 0.14503204;
    wl = -1.2684380046;
    wm = 2.6097574011;
    ws = -0.3413193965;
  } else {
    // Blue component
    k0 = 1.35733652;
    k1 = -0.00915799;
    k2 = -1.15130210;
    k3 = -0.50559606;
    k4 = 0.00692167;
    wl = -0.0041960863;
    wm = -0.7034186147;
    ws = 1.7076147010;
  }

  // Approximate max saturation using a polynomial:
  double S = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b;

  // Do one step Halley's method to get closer
  // this gives an error less than 10e6, except for some blue hues where the dS/dh is close to infinite
  // this should be sufficient for most applications, otherwise do two/three steps

  double kL = 0.3963377774 * a + 0.2158037573 * b;
  double kM = -0.1055613458 * a - 0.0638541728 * b;
  double kS = -0.0894841775 * a - 1.2914855480 * b;

  {
    double l_ = 1 + S * kL;
    double m_ = 1 + S * kM;
    double s_ = 1 + S * kS;

    double l = l_ * l_ * l_;
    double m = m_ * m_ * m_;
    double s = s_ * s_ * s_;

    double lDs = 3 * kL * l_ * l_;
    double mDs = 3 * kM * m_ * m_;
    double sDs = 3 * kS * s_ * s_;

    double lDs2 = 6 * kL * kL * l_;
    double mDs2 = 6 * kM * kM * m_;
    double sDs2 = 6 * kS * kS * s_;

    double f = wl * l + wm * m + ws * s;
    double f1 = wl * lDs + wm * mDs + ws * sDs;
    double f2 = wl * lDs2 + wm * mDs2 + ws * sDs2;

    S = S - f * f1 / (f1 * f1 - 0.5 * f * f2);
  }

  return S;
}