(* Copyright (C) 1992, Digital Equipment Corporation                         *)
(* All rights reserved.                                                      *)
(* See the file COPYRIGHT for a full description.                            *)

(* Created by stolfi on Tue Apr 25 20:50:11 1989               *)
(* Last modified on Tue Feb 11 21:39:39 PST 1992 by muller     *)
(*      modified on Thu Oct 25 10:33:13 PDT 1990 by stolfi     *)
(*      modified on Tue Feb 27 00:00:44 1990 by harrison       *)

MODULE RGBScale;

IMPORT Color, RGB, RGBDist, Intensity, IntensityScale;

PROCEDURE Straight(
    READONLY a, b: RGB.T; 
    READONLY s: IntensityScale.T;
  ): REF T =
  VAR cc: REF T;
  BEGIN
    cc := NEW(REF T, NUMBER(s));
    FOR i := 0 TO LAST(cc^) DO
      IF s[i] = Intensity.Undefined THEN
        cc[i] := RGB.Undefined
      ELSE
        cc[i] := RGB.Mix(a, 1.0 - s[i], b, s[i])
      END;
    END;
    RETURN cc
  END Straight;

PROCEDURE Grey(
    n: CARDINAL; 
    a, b: Intensity.T; 
    ratio: REAL := 1.0
  ): REF T =
  BEGIN
    RETURN Straight(RGB.Black, RGB.White, IntensityScale.New(n, a, b, ratio)^)
  END Grey;

PROCEDURE Grid (READONLY rs, gs, bs: IntensityScale.T): REF T =
  VAR i: CARDINAL; 
      c: RGB.T; 
      cc: REF ARRAY OF RGB.T;
      NR, NG, NB: CARDINAL;
  BEGIN
    NR := NUMBER(rs);
    NG := NUMBER(gs);
    NB := NUMBER(bs);
    cc := NEW(REF T, NR * NG * NB);
    i := 0;
    FOR b := 0 TO NB - 1 DO
      c[2] := bs[b];
      FOR g := 0 TO NG - 1 DO
        c[1] := gs[g];
        FOR r := 0 TO NR - 1 DO 
          c[0] := rs[r]; 
          cc[i] := c;
          INC(i)
        END
      END
    END;
    RETURN cc
  END Grid;

PROCEDURE Path(
    READONLY s: ARRAY OF REAL; 
    READONLY p: ARRAY OF RGB.T;
    dist: Color.DistFn := RGBDist.Perceptual;
    lumWeight: REAL := 1.0;
  ): REF T =
  VAR
    len: REF IntensityScale.T; (* len[i] is the length of path p to p[i] *)
    r: REF T;                  (* r[j] is the jth color of current scale *)
    t: REF IntensityScale.T;   (* t[j] is the position of r[j] along p *)
    rb: REF T;                 (* Best r[j] found so far *)
    tb: REF IntensityScale.T;  (* t of best r found so far *)
    tn: REF IntensityScale.T;  (* t of next r *)
    i, m, n, nit, k: CARDINAL;
    sj, err, errb, leni, ti, bold, c: REAL;
  CONST MaxTries = 3;
  BEGIN
    n := NUMBER(s);
    m := NUMBER(p);
    <* ASSERT m >= 2 *>
    r := NEW(REF T, n);
    t := NEW(REF IntensityScale.T, n);
    IF n = 0 THEN RETURN r END;
    (* Compute cumulative length of given path, normalized to 1.0: *)
    len := NEW(REF IntensityScale.T, m);
    len[0] := 0.0;
    FOR i := 1 TO m - 1 DO
      len[i] := len[i - 1] + dist(p[i], p[i - 1], lumWeight)
    END;
    <* ASSERT len[m - 1] > 0.00001 *>
    FOR i := 1 TO m - 2 DO len[i] := len[i] / len[m - 1] END;
    len[m - 1] := 1.0;
    (* Compute initial guess by looking up s[j] in path length: *)
    i := 0;
    FOR j := 0 TO n - 1 DO
      sj := s[j];
      <* ASSERT (0.0 <= sj) AND (sj <= 1.0) *>
      WHILE sj < len[i] DO DEC(i) END;
      WHILE sj > len[i + 1] DO INC(i) END;
      t[j] := FLOAT(i);
      leni := len[i];
      IF (len[i + 1] - len[i] > 0.00001) THEN
        t[j] := t[j] + (sj - leni) / (len[i + 1] - leni)
      END;
      r[j] := RGB.Interpolate(p[i], sj - len[i], p[i + 1], len[i + 1] - sj)
    END;
    (* Now refine guess: *)
    rb := NIL;
    tb := NIL;
    errb := 1.0e30;
    len := NEW(REF IntensityScale.T, n + 1);
    tn := NEW(REF IntensityScale.T, n);
    nit := 0;
    bold := 0.5;
    LOOP
      (* Compute actual cumulative distance along scale: *)
      len[0] := dist(r[0], p[0], lumWeight);
      FOR j := 1 TO n - 1 DO
        len[j] := len[j - 1] + dist(r[j], r[j - 1], lumWeight)
      END;
      len[n] := len[n - 1] + dist(p[m - 1], r[n - 1], lumWeight);
      <* ASSERT len[n] > 0.00001 *>
      FOR j := 1 TO n - 1 DO len[j] := len[j] / len[n] END;
      len[n] := 1.0;
      (* Compare with desired distances: *)
      err := 0.0;
      FOR j := 0 TO n - 1 DO
        err := MAX(err, ABS(len[j] - s[j]));
      END;
      IF err < errb THEN
        errb := err;
        rb := r;
        tb := t;
        bold := 1.5 * bold
      ELSE
        bold := bold / 2.0
      END;
      IF (errb < 1.0000000e-03) OR (nit >= MaxTries) THEN RETURN rb END;
      (* Compute new positions by interpolating desired values: *)
      i := 0;
      FOR j := 0 TO n - 1 DO
        sj := s[j];
        WHILE (i > 0) AND (sj < len[i]) DO DEC(i) END;
        WHILE sj > len[i + 1] DO INC(i) END;
        IF i < 0 THEN
          leni := 0.0;
          ti := 0.0
        ELSE
          leni := len[i];
          ti := t[i]
        END;
        tn[j] := ti;
        IF (len[i + 1] - leni > 0.00001) THEN
          tn[j] := tn[j] + (t[i + 1] - ti) * (sj - leni) / (len[i + 1] - leni)
        END;
        (* Merge with best solution: *)
        tn[j] := bold * tn[j] + (1.0 - bold) * tb[j];
        k := TRUNC(tn[j]);
        k := MIN(m - 2, MAX(0, k));
        c := tn[j] - FLOAT(k);
        r[j] := RGB.Interpolate(p[k], c, p[k + 1], 1.0 - c);
      END;
      FOR i := 0 TO n - 1 DO t[i] := tn[i] END;
      INC(nit);
    END
  END Path;

BEGIN

END RGBScale.