findGamutIntersection function

double findGamutIntersection(double a, double b, double L1, double C1, double L0, LC cusp,)

Implementation

double findGamutIntersection(double a, double b, double L1, double C1, double L0, LC cusp) {
  // Find the intersection for upper and lower half separately
  double t;
  if (((L1 - L0) * cusp.C - (cusp.L - L0) * C1) <= 0) {
    // Lower half
    t = cusp.C * L0 / (C1 * cusp.L + cusp.C * (L0 - L1));
  } else {
    // Upper half

    // First intersect with triangle
    t = cusp.C * (L0 - 1) / (C1 * (cusp.L - 1) + cusp.C * (L0 - L1));

    // Then one step Halley's method
    {
      double dL = L1 - L0;
      double dC = C1;

      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 lDt = dL + dC * kL;
      double mDt = dL + dC * kM;
      double sDt = dL + dC * kS;

      // If higher accuracy is required, 2 or 3 iterations of the following block can be used:
      {
        double L = L0 * (1 - t) + t * L1;
        double C = t * C1;

        double l_ = L + C * kL;
        double m_ = L + C * kM;
        double s_ = L + C * kS;

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

        double ldt = 3 * lDt * l_ * l_;
        double mdt = 3 * mDt * m_ * m_;
        double sdt = 3 * sDt * s_ * s_;

        double ldt2 = 6 * lDt * lDt * l_;
        double mdt2 = 6 * mDt * mDt * m_;
        double sdt2 = 6 * sDt * sDt * s_;

        double r = 4.0767416621 * l - 3.3077115913 * m + 0.2309699292 * s - 1;
        double r1 = 4.0767416621 * ldt - 3.3077115913 * mdt + 0.2309699292 * sdt;
        double r2 = 4.0767416621 * ldt2 - 3.3077115913 * mdt2 + 0.2309699292 * sdt2;

        double uR = r1 / (r1 * r1 - 0.5 * r * r2);
        double tR = -r * uR;

        double g = -1.2684380046 * l + 2.6097574011 * m - 0.3413193965 * s - 1;
        double g1 = -1.2684380046 * ldt + 2.6097574011 * mdt - 0.3413193965 * sdt;
        double g2 = -1.2684380046 * ldt2 + 2.6097574011 * mdt2 - 0.3413193965 * sdt2;

        double uG = g1 / (g1 * g1 - 0.5 * g * g2);
        double tG = -g * uG;

        double b = -0.0041960863 * l - 0.7034186147 * m + 1.7076147010 * s - 1;
        double b1 = -0.0041960863 * ldt - 0.7034186147 * mdt + 1.7076147010 * sdt;
        double b2 = -0.0041960863 * ldt2 - 0.7034186147 * mdt2 + 1.7076147010 * sdt2;

        double uB = b1 / (b1 * b1 - 0.5 * b * b2);
        double tB = -b * uB;

        tR = uR >= 0 ? tR : double.infinity;
        tG = uG >= 0 ? tG : double.infinity;
        tB = uB >= 0 ? tB : double.infinity;

        t += math.min(tR, math.min(tG, tB));
      }
    }
  }

  return t;
}