#! /usr/bin/gawk -f
# Last edited on 2007-01-10 12:17:07 by stolfi

BEGIN {
  abort = -1;
  usage = ( \
    "make-symmetric-palette \\\n" \
    "  -v ini={iR},{iG},{iB} \\\n" \
    "  -v fin={fR},{fG},{fB} \\\n" \
    "  [ -v debug={BOOL} ] \\\n" \
    "  [ -v showlum={BOOL} ] \\\n" \
    "  > OUTFILE.txt" \
  ); 

  # Produces a color palette for signed data that 
  # maps positive colors along a path from {ini} to {fin}, 
  # and negative colors along a path with complementary hues.
  # Zero is mapped to a gray color with same luminosity as {ini}.
  # 
  # The output has one color per line, in the format 
  #    "{V}  {R} {G} {B}"
  # where {V} is a number, between -1 and +1, that
  # should be mapped to the specified color.
  #
  # If {showlum} is nonzero, it also prints the luminosity as a comment
  #   "# Y = {LUM}"
  #
  # The table {t} is such that {t(-V)} and {t(V)} have the same
  # brightness, opposite hues, and the same relative saturation.
  
  if (debug == "") { debug = 0; }
  if (showlum == "") { showlum = 0; }
  if (ini == "") { arg_error("must define \"ini\""); }
  if (fin == "") { arg_error("must define \"fin\""); }
  
  split ("", tmp);

  # Parse initial color, extract initial hue and lum:
  parse_rgb(ini, tmp); iR = tmp[0]; iG = tmp[1]; iB = tmp[2];
  if (debug) { debug_color("ini color", iR, iG, iB); }
  iH = hue(iR,iG,iB);
  if (debug) { debug_scalar("ini hue (iH)", iH); }
  iY = lum(iR,iG,iB);
  if (debug) { debug_scalar("ini lum (iY)", iY); }

  # Parse final color, extract final hue and lum:
  parse_rgb(fin, tmp); fR = tmp[0]; fG = tmp[1]; fB = tmp[2];
  if (debug) { debug_color("fin color", fR, fG, fB); }
  fH = hue(fR,fG,fB);
  if (debug) { debug_scalar("fin hue (fH)", fH); }
  fY = lum(fR,fG,fB);
  if (debug) { debug_scalar("fin lum (fY)", fY); }
  
  # Adjust final hue so that {|ini-fin| <= 1/2}:
  while (fH > iH + 0.5) { fH = fH - 1; }
  while (fH < iH - 0.5) { fH = fH + 1; }
  if (debug) { debug_scalar("adj fin hue (fH)", fH); }
  
  # Interpolate gradually between hues and lums:
  split("", pR); split("", pG); split("", pB); 
  np = 7;  # The palette will have {2*np+1} entries, {pR[-np..+np]}.
  pow2 = 1; # {pow2} is {2^i}.
  for (i = 1; i <= np; i++)
    { # Update {pow2}:
      pow2 = pow2 * 2;
      
      # Fraction of the way along the path:
      t = (i + 0.0)/(np + 0.0);
      if (debug) { debug_scalar("path arg (t)", t); }
      
      # Current hue {tH}:
      tH = (1-t)*iH + t*fH;
      if (debug) { debug_scalar("hue (tH)", tH); }

      # Current luminosity {tY}:
      tY = (1-t)*iY + t*fY;
      if (debug) { debug_scalar("lum (tY)", tY); }
      
      # Compute desired relative saturation {tS}:
      tS = 1.0;
      if (debug) { debug_scalar("sat (tS)", tS); }
      
      # Compute colorwith hue {tH}, lum {tY} and relsat {S}:
      color_from_HYS(tH,tY,tS, tmp);
      
      # Save in table:
      pR[+i] = tmp[0]; pG[+i] = tmp[1]; pB[+i] = tmp[2]; 
      
      if(debug) { printf "\n" > "/dev/stderr"; }

      # Compute supplementary color with hue {tH+0.5}, lum {tY} and relsat {S}:
      color_from_HYS(tH+0.5,tY,tS, tmp);  tR = tmp[0]; tG = tmp[2]; tB = tmp[3];
      
      # Save in table:
      pR[-i] = tmp[0]; pG[-i] = tmp[1]; pB[-i] = tmp[2]; 
      
      if(debug) { printf "\n\n" > "/dev/stderr"; }
    }
    
  # Compute brightness {mY} of middle entry:
  mY = 1.1*iY - 0.1*fY;
  if (mY < 0) { mY = 0; }
  if (mY > 1) { mY = 1; }
  
  # Save {mY}-gray in table:
  pR[0] = mY; pG[0] = mY; pB[0] = mY;
  
  # Dump table:
  for (i = -np; i <= np; i++)
    { # Fraction of the way along the path:
      t = (i + 0.0)/(np + 0.0);
      
      # Print entry: 
      printf "%+5.3f  %5.3f %5.3f %5.3f", t, pR[i], pG[i], pB[i];
      if (showlum) { printf " # Y = %6.4f", lum(pR[i], pG[i], pB[i]); }
      printf "\n"; 
    }
}

function color_from_HYS(H,Y,S,vec,  z,q,R,G,B,L,r0,r1)
{
  # Split {H} into sextant index {q} and fraction {z}:
  z = H*6;
  q = int(z);
  z = z - q;
  while (z < 0) { z = z + 1; q = q - 1; }
  while (z >= 1) { z = z - 1; q = q + 1; }
  q = q % 6;
  while (q < 0) { q += 6; }
  while (q >= 6) { q -= 6; }
  
  # Compute ``spectral'' color of hue {H}:
  if (q == 0) { R = 1;   G = z;   B = 0;   }
  if (q == 1) { R = 1-z; G = 1;   B = 0;   }
  if (q == 2) { R = 0;   G = 1;   B = z;   }
  if (q == 3) { R = 0;   G = 1-z; B = 1;   }
  if (q == 4) { R = z;   G = 0;   B = 1;   }
  if (q == 5) { R = 1;   G = 0;   B = 1-z; }
  if (debug) { debug_color("spectral color", R, G, B); }

  # Compute its luminosity:
  L = lum(R,G,B);
  
  # Mix {R,G,B} with white or black so that ist lum is {Y}:
  if (L < Y)
    { # Mix with white: 
      r0 = (Y - L)/(1 - L); r1 = 1-r0;
      R = r1*R + r0; 
      G = r1*G + r0; 
      B = r1*B + r0; 
    }
  else
    { # Mix with black: 
      r1 = Y/L;
      R = r1*R; 
      G = r1*G; 
      B = r1*B;
    }
  if (debug) { debug_color("lum-adj color", R, G, B); }

  # Adjust {R,G,B} to have relative saturation {S}:
  R = S*R + (1-S)*Y;
  G = S*G + (1-S)*Y;
  B = S*B + (1-S)*Y;
  
  # Clip to unit cube:
  R = clip_value(R);
  G = clip_value(G);
  B = clip_value(B);
  if (debug) { debug_color("table color", R, G, B); }

  # Return in {vec}:
  vec[0] = R; vec[1] = G; vec[2] = B;
}      

function parse_rgb(str,vec,  n,fld)
{
  # Parses the string {str} as a color value
  # (three numbers separated by commas or spaces,
  # possibly delimited by parentheses, bracket, braces, etc..
  # Returns the result in {vec[0]} to {vec[2]}.
  
  # Remove delimiters:
  gsub(/[][(){}<>,\011\015]/, " ", str);
  
  # Check validity of format:
  
  # Split at blanks: 
  n = split(str,fld);
  if (n != 3) { data_error(("bad color = \"" $0 "\"\n")); }
  
  # Parse RGB components:
  vec[0] = check_value(fld[1]);
  vec[1] = check_value(fld[2]); 
  vec[2] = check_value(fld[3]);
}
  
function round(x  )
{
  # Rounds {x} to integer:
  return int(x + (x < 0 ? -0.5 : +0.5));
}
  
function hue(R,G,B)
{
  # If {R==G==B}, the hue is undefined (grey):
  if ((R == G) && (G == B)) { return 0; }
  
  # Othwerwise these computations should work:
  if ((R >= G) && (G >= B)) { return (0 + (G-B)/(R-B))/6; }
  if ((G >= R) && (R >= B)) { return (1 + (G-R)/(G-B))/6; }
  
  if ((G >= B) && (B >= R)) { return (2 + (B-R)/(G-R))/6; }
  if ((B >= G) && (G >= R)) { return (3 + (B-G)/(B-R))/6; }
  
  if ((B >= R) && (R >= G)) { return (4 + (R-G)/(B-G))/6; }
  if ((R >= B) && (B >= G)) { return (5 + (R-B)/(R-G))/6; }
  
  # We should never get here:
  assert(0);
  return 0;
}
  
function lum(R,G,B)
{
  # Luminosity of {R,G,B} (European TV formula):
  return 0.298911*R + 0.586611*G + 0.114478*B;
}

function sat(R,G,B,  Y,S,sR,sG,sB)
{
  Y = lum(R,G,B);
  S = 0;
  sR = sat1(R-Y,Y); if (sR > S) { S = sR; }
  sG = sat1(G-Y,Y); if (sG > S) { S = sG; }
  sB = sat1(B-Y,Y); if (sB > S) { S = sB; }
  return S;
}

function sat1(C,Y  )
{
  # Saturation of {C} relative to the range {[-Y _ 1-Y]}
  return (C < 0 ? -C/Y : C/(1-Y));
}

function check_value(x,  v)
{
  # Checks format of {x}.
  if (x !~ /^[ ]*[-+]?[0-9]*([.][0-9]*|[0-9])[ ]*$/)
    { data_error(("bad number = \"" x "\"\n")); }
  v = x + 0.000;
  return v;
}

function clip_value(V,  eps)
{
  # Returns V clipped to range [0_1]; warns if outside.
  eps = 1.0e-5; # Fudge factor to allow for roundoff errors.
  if ((V < -eps) || (V > 1+eps))
    { data_warning(("value " V " not in unit cube")); }
  if (V < 0) { V = 0; } else if (V > 1) { V = 1; }
  return V;
}

function assert(cond)
{
  if (! cond) { prog_error("bug"); }
}

function data_warning(msg)
{
  printf "%s:%s: warning - %s\n", FILENAME, FNR, msg > "/dev/stderr";
}

function arg_error(msg)
{
  printf "** %s\n", msg > "/dev/stderr";
  printf "usage: %s\n", usage;
  abort = 1;
  exit abort;
}

function data_error(msg)
{
  printf "%s:%s: ** %s\n", FILENAME, FNR, msg > "/dev/stderr";
  abort = 1;
  exit abort;
}

function prog_error(msg)
{
  printf "** %s\n", msg > "/dev/stderr";
  abort = 1;
  exit abort;
}

function debug_color(lab,R,G,B)
{
  printf "%-20s = ( %+6.3f %+6.3f %+6.3f )", lab, R, G, B > "/dev/stderr";
  printf " hue = %5.3f", hue(R,G,B) > "/dev/stderr";
  printf " lum = %5.3f", lum(R,G,B) > "/dev/stderr";
  printf " sat = %5.3f", sat(R,G,B) > "/dev/stderr";
  printf "\n" > "/dev/stderr";
}

function debug_scalar(lab,fmt,V)
{
  printf "%-20s = ", lab > "/dev/stderr";
  printf fmt, V > "/dev/stderr";
  printf "\n" > "/dev/stderr";
}