metrics.h

Go to the documentation of this file.

 
#pragma once
 
#include <cassert>
#include <cmath>
#include <qualpal/colors.h>
 
namespace qualpal {
 
namespace metrics {
namespace detail {
inline double
cosd(double degrees)
{
  return std::cos(degrees * M_PI / 180.0);
}
 
inline double
sind(double degrees)
{
  return std::sin(degrees * M_PI / 180.0);
}
 
inline double
tand(double degrees)
{
  return std::tan(degrees * M_PI / 180.0);
}
 
inline double
atan2d(double y, double x)
{
  double deg = std::atan2(y, x) * 180.0 / M_PI;
  return (deg >= 0) ? deg : deg + 360.0;
}
 
inline double
square(double x)
{
  return x * x;
}
} // namespace detail
 
// Forward declarations
class RGB;
class HSL;
class XYZ;
class Lab;
class DIN99d;
 
enum class MetricType
{
  DIN99d,   
  CIE76,    
  CIEDE2000 
};
 
class DIN99d
{
private:
  bool use_power_transform;
  double power;
  double scale;
 
public:
  explicit DIN99d(bool use_power_transform = true,
                  double power = 0.74,
                  double scale = 1.28)
    : use_power_transform(use_power_transform)
    , power(power)
    , scale(scale)
  {
  }
 
  template<typename ColorType1, typename ColorType2>
  double operator()(const ColorType1& c1, const ColorType2& c2) const
  {
    colors::DIN99d d1(c1), d2(c2);
    double d = std::hypot(d1.l() - d2.l(), d1.a() - d2.a(), d1.b() - d2.b());
 
    if (use_power_transform) {
      return std::pow(d, power) * scale;
    } else {
      return d;
    }
  }
};
 
struct CIE76
{
  template<typename ColorType1, typename ColorType2>
  double operator()(const ColorType1& c1, const ColorType2& c2) const
  {
    colors::Lab l1(c1), l2(c2);
    return std::hypot(l1.l() - l2.l(), l1.a() - l2.a(), l1.b() - l2.b());
  }
};
 
class CIEDE2000
{
private:
  double K_L;
  double K_C;
  double K_H;
 
public:
  explicit CIEDE2000(double K_L = 1.0, double K_C = 1.0, double K_H = 1.0)
    : K_L(K_L)
    , K_C(K_C)
    , K_H(K_H)
  {
    assert((K_L > 0 && K_C > 0 && K_H > 0) &&
           "CIEDE2000 weighting factors must be positive");
  }
 
  template<typename ColorType1, typename ColorType2>
  double operator()(const ColorType1& c1, const ColorType2& c2) const
  {
    using namespace detail;
 
    colors::Lab x(c1), y(c2);
 
    double L_hat_prime = (x.l() + y.l()) / 2.0;
    double C1 = std::hypot(x.a(), x.b());
    double C2 = std::hypot(y.a(), y.b());
    double C_hat = (C1 + C2) / 2.0;
    double G = 0.5 * (1 - std::sqrt(std::pow(C_hat, 7) /
                                    (std::pow(C_hat, 7) + std::pow(25.0, 7))));
    double a1_prime = x.a() * (1.0 + G);
    double a2_prime = y.a() * (1.0 + G);
    double C1_prime = std::hypot(a1_prime, x.b());
    double C2_prime = std::hypot(a2_prime, y.b());
    double C_hat_prime = (C1_prime + C2_prime) / 2.0;
 
    double h1_prime = atan2d(x.b(), a1_prime);
 
    if (h1_prime < 0) {
      h1_prime += 360;
    }
 
    double h2_prime = atan2d(y.b(), a2_prime);
 
    if (h2_prime < 0) {
      h2_prime += 360;
    }
 
    double H_hat_prime = std::abs(h1_prime - h2_prime) > 180
                           ? (h1_prime + h2_prime + 360) / 2.0
                           : (h1_prime + h2_prime) / 2.0;
 
    double T = 1.0 - 0.17 * cosd(H_hat_prime - 30) +
               0.24 * cosd(2 * H_hat_prime) + 0.32 * cosd(3 * H_hat_prime + 6) -
               0.20 * cosd(4 * H_hat_prime - 63);
 
    double delta_h_prime = h2_prime - h1_prime;
    if (std::abs(delta_h_prime) > 180) {
      if (h2_prime <= h1_prime)
        delta_h_prime += 360;
      else
        delta_h_prime -= 360;
    }
 
    double delta_L_prime = y.l() - x.l();
    double delta_C_prime = C2_prime - C1_prime;
 
    double delta_H_prime =
      2 * std::sqrt(C1_prime * C2_prime) * sind(delta_h_prime / 2.0);
    double S_L = 1 + (0.015 * std::pow(L_hat_prime - 50, 2)) /
                       std::sqrt(20 + std::pow(L_hat_prime - 50, 2));
 
    double S_C = 1 + 0.045 * C_hat_prime;
    double S_H = 1 + 0.015 * C_hat_prime * T;
 
    double delta_theta = 30 * std::exp(-std::pow((H_hat_prime - 275) / 25, 2));
 
    double R_C = 2 * std::sqrt(std::pow(C_hat_prime, 7) /
                               (std::pow(C_hat_prime, 7) + std::pow(25.0, 7)));
 
    double R_T = -R_C * sind(2 * delta_theta);
    double out = std::sqrt(square(delta_L_prime / (K_L * S_L)) +
                           square(delta_C_prime / (K_C * S_C)) +
                           square(delta_H_prime / (K_H * S_H)) +
                           R_T * (delta_C_prime / (K_C * S_C)) *
                             (delta_H_prime / (K_H * S_H)));
 
    assert(out >= 0 && "CIEDE2000 color difference must be non-negative");
    assert(std::isfinite(out) &&
           "CIEDE2000 color difference must not be finite");
 
    return out;
  }
};
} // namespace metrics
} // namespace qualpal