GamutCompress.blink

kernel GamutCompression : public ImageComputationKernel<ePixelWise> {
  Image<eRead, eAccessPoint, eEdgeClamped> src;
  Image<eWrite> dst;

  param:
    float3 threshold;
    float p;
    float cyan;
    float magenta;
    float yellow;
    float p_Height;
    bool invert;
    bool overlay;

  local:
  float3 thr;
  float3 lim;

  void init() {
    // thr is the percentage of the core gamut to preserve
    thr = float3(min(0.9999f, threshold.x), min(0.9999f, threshold.y), min(0.9999f, threshold.z));

    // lim is the max distance from the gamut boundary that will be compressed
    // 0 is a no-op, 1 will compress colors from a distance of 2.0 from achromatic to the gamut boundary
    lim = float3(cyan+1.0f, magenta+1.0f, yellow+1.0f);
  }

  // calculate compressed distance
  float compress(float x, float l, float t) {
    float cx, s;
    if (x < t) {
      cx = x;
    } else {
      // power(p) compression function plot https://www.desmos.com/calculator/54aytu7hek
      if (l < 1.0001) {
        return x; // disable compression, avoid nan
      }
      s = (l-t)/pow(pow((1.0f-t)/(l-t),-p)-1.0f,1.0f/p); // calc y=1 intersect
      if (invert == 0) {
        cx = t+s*((x-t)/s)/(pow(1.0f+pow((x-t)/s,p),1.0f/p)); // compress
      } else {
        if (x > (t + s)) {
          cx = x; // avoid singularity
        } else {
          cx = t+s*pow(-(pow((x-t)/s,p)/(pow((x-t)/s,p)-1.0f)),1.0f/p); // uncompress
        }
      }
    }
    return cx;
  }


  void process(int2 pos) {
    // source pixels
    SampleType(src) rgba = src();
    float3 rgb = float3(rgba.x, rgba.y, rgba.z);

    // normalised pixel coordinates
    float2 fpos = float2((float)pos.x / dst.bounds.width(), (float)pos.y / p_Height);

    // achromatic axis 
    float ach = max(rgb.x, max(rgb.y, rgb.z));

    // distance from the achromatic axis for each color component aka inverse rgb ratios
    // distance is normalized by achromatic, so that 1.0f is at gamut boundary, avoid 0 div
    float3 dist = ach == 0 ? float3(0.0f, 0.0f, 0.0f) : (ach-rgb)/fabs(ach);

    // compress distance with user controlled parameterized shaper function
    float sat;
    float3 csat, cdist;
    cdist = float3(
      compress(dist.x, lim.x, thr.x),
      compress(dist.y, lim.y, thr.y),
      compress(dist.z, lim.z, thr.z));

    // recalculate rgb from compressed distance and achromatic
    // effectively this scales each color component relative to achromatic axis by the compressed distance
    float3 crgb = ach-cdist*fabs(ach);

    // Graph overlay method based on one by Paul Dore
    // https://github.com/baldavenger/DCTLs/tree/master/ACES%20TOOLS
    if (overlay) {
      float3 cramp = float3(
        compress(2.0f * fpos.x, lim.x, thr.x),
        compress(2.0f * fpos.x, lim.y, thr.y),
        compress(2.0f * fpos.x, lim.z, thr.z));
      bool overlay_r = fabs(2.0f * fpos.y - cramp.x) < 0.004f || fabs(fpos.y - 0.5f) < 0.0005f ? true : false;
      bool overlay_g = fabs(2.0f * fpos.y - cramp.y) < 0.004f || fabs(fpos.y - 0.5f) < 0.0005f ? true : false;
      bool overlay_b = fabs(2.0f * fpos.y - cramp.z) < 0.004f || fabs(fpos.y - 0.5f) < 0.0005f ? true : false;
      crgb.x = overlay_g || overlay_b ? 1.0f : crgb.x;
      crgb.y = overlay_b || overlay_r ? 1.0f : crgb.y;
      crgb.z = overlay_r || overlay_g ? 1.0f : crgb.z;
    }

    // write to output
    dst() = float4(crgb.x, crgb.y, crgb.z, rgba.w);
  }
};