GamutCompress.blink

/*  GamutCompress v0.7
    Written by Jed Smith with lots of help from the ACES Gamut Mapping VWG
    https://github.com/jedypod
    https://community.acescentral.com/t/rgb-saturation-gamut-mapping-approach-and-a-comp-vfx-perspective
    https://community.acescentral.com/c/aces-development-acesnext/vwg-aces-gamut-mapping-working-group
*/
kernel GamutCompression : public ImageComputationKernel<ePixelWise> {
  Image<eRead, eAccessPoint, eEdgeClamped> src;
  Image<eWrite> dst;

  param:
    float3 th;
    float3 dl;
    bool invert;

  local:
    float3 s;

  void init() {
    // Pre-calculate scale so compression function passes through distance limit: (x=dl, y=1)
    s = ((float3)(1.0f)-th)/sqrt(max((float3)(1.001f),dl)-(float3)(1.0f));
  }

  void process() {
    float3 rgb = float3(src().x, src().y, src().z);

    float ac = max(rgb.x, max(rgb.y, rgb.z)); // Achromatic axis
    // Inverse RGB Ratios: distance from achromatic axis
    float3 d = ac == 0.0f ? (float3)(0.0f) : (ac-rgb)/fabs(ac);

    float3 cd; // Compressed distance
    // Parabolic compression function: https://www.desmos.com/calculator/nvhp63hmtj
    if (!invert) {
      cd.x = d.x<th.x?d.x:s.x*sqrt(d.x-th.x+s.x*s.x/4.0f)-s.x*sqrt(s.x*s.x/4.0f)+th.x;
      cd.y = d.y<th.y?d.y:s.y*sqrt(d.y-th.y+s.y*s.y/4.0f)-s.y*sqrt(s.y*s.y/4.0f)+th.y;
      cd.z = d.z<th.z?d.z:s.z*sqrt(d.z-th.z+s.z*s.z/4.0f)-s.z*sqrt(s.z*s.z/4.0f)+th.z;
    } else {
      cd.x = d.x<th.x?d.x:pow(d.x-th.x+s.x*sqrt(s.x*s.x/4.0f),2.0f)/(s.x*s.x)-s.x*s.x/4.0f+th.x;
      cd.y = d.y<th.y?d.y:pow(d.y-th.y+s.y*sqrt(s.y*s.y/4.0f),2.0f)/(s.y*s.y)-s.y*s.y/4.0f+th.y;
      cd.z = d.z<th.z?d.z:pow(d.z-th.z+s.z*sqrt(s.z*s.z/4.0f),2.0f)/(s.z*s.z)-s.z*s.z/4.0f+th.z;
    }

    rgb = ac-cd*fabs(ac);
    dst() = float4(rgb.x, rgb.y, rgb.z, src().w);
  }
};