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

(* Created by Jorge Stolfi on Sat Oct  7  2:51:33 PDT 1989     *)
(* Last modified on Tue Feb 11 21:39:43 PST 1992 by muller     *)
(*      modified on Mon Nov 11 16:15:12 PST 1991 by steveg     *)
(*      modified on Thu Oct 25  8:58:07 PDT 1990 by stolfi     *)

MODULE LDW;

IMPORT Intensity, Math, RGB;
FROM YPQ IMPORT
  rY, rP, rQ,  gY, gP, gQ,  bY, bP, bQ,
  Yr, Yg, Yb,  Pr, Pg, Pb,  Qr, Qg, Qb;

CONST
  epsY = 0.05;  (* Brightness threshold for log behavior of Lightness *)
  magY = 0.271039152013153; (* 1.0/asinh(1.0/epsY) *)
  offY = 0.10;  (* Brightness threshold for ratio behavior of Chroma *)
  epsC = 0.20;  (* P/Y and Q/Y threshold for log behavior of Chroma *)
  magC = 0.432443962786746; (* 1.0/asinh(1.0/epsC) *)

PROCEDURE YFromL (l: REAL): REAL =
  VAR y: REAL;
  BEGIN
    IF l = Intensity.Undefined THEN RETURN Intensity.Undefined END;
    y := epsY * FLOAT(Math.sinh(FLOAT(l / magY, LONGREAL)));
    RETURN y
  END YFromL;

PROCEDURE LFromY (y: REAL): REAL =
  BEGIN
    IF y = Intensity.Undefined THEN RETURN Intensity.Undefined END;
    RETURN FLOAT(Math.asinh(FLOAT(y / epsY, LONGREAL))) * magY;
  END LFromY;

PROCEDURE LFromRGB (READONLY rgb: RGB.T): REAL =
  VAR y: REAL;
  BEGIN
    IF rgb = RGB.Undefined THEN RETURN Intensity.Undefined END;
    y := rgb[0] * rY + rgb[1] * gY + rgb[2] * bY;
    RETURN LFromY(y);
  END LFromRGB;

PROCEDURE RGBFromL (l: REAL): RGB.T =
  VAR y: REAL;
  BEGIN
    IF l = Intensity.Undefined THEN RETURN RGB.Undefined END;
    y := YFromL(l);
    RETURN T{y, y, y}
  END RGBFromL;

PROCEDURE DFromRGB (READONLY rgb: RGB.T): REAL =
  VAR y, p: REAL;
  BEGIN
    IF rgb = RGB.Undefined THEN RETURN Intensity.Undefined END;
    y := rgb[0] * rY + rgb[1] * gY + rgb[2] * bY;
    p := rgb[0] * rP + rgb[1] * gP + rgb[2] * bP;
    RETURN FLOAT(Math.asinh(FLOAT((p / FLOAT(Math.hypot(FLOAT(y, LONGREAL), FLOAT(offY, LONGREAL)))) / epsC, LONGREAL))) * magC;
  END DFromRGB;

PROCEDURE WFromRGB (READONLY rgb: RGB.T): REAL =
  VAR y, q: REAL;
  BEGIN
    IF rgb = RGB.Undefined THEN RETURN Intensity.Undefined END;
    y := rgb[0] * rY + rgb[1] * gY + rgb[2] * bY;
    q := rgb[0] * rQ + rgb[1] * gQ + rgb[2] * bQ;
    RETURN FLOAT(Math.asinh(FLOAT((q / FLOAT(Math.hypot(FLOAT(y, LONGREAL), FLOAT(offY, LONGREAL)))) / epsC, LONGREAL))) * magC;
  END WFromRGB;

PROCEDURE ChromaFromRGB (READONLY rgb: RGB.T): REAL =
  VAR ldw: T;
  BEGIN
    IF rgb = RGB.Undefined THEN RETURN Intensity.Undefined END;
    ldw := LDWFromRGB (rgb);
    RETURN FLOAT(Math.hypot(FLOAT(ldw[1], LONGREAL), FLOAT(ldw[2], LONGREAL)))
  END ChromaFromRGB;

PROCEDURE LDWFromRGB (READONLY rgb: RGB.T): T =
  VAR y, yc, p, q: REAL;
  BEGIN
    IF rgb = RGB.Undefined THEN RETURN Undefined END;
    y := rgb[0] * rY + rgb[1] * gY + rgb[2] * bY;
    p := rgb[0] * rP + rgb[1] * gP + rgb[2] * bP;
    q := rgb[0] * rQ + rgb[1] * gQ + rgb[2] * bQ;
    yc := FLOAT(Math.hypot(FLOAT(y, LONGREAL), FLOAT(offY, LONGREAL)));
    RETURN T{
      FLOAT(Math.asinh(FLOAT(y / epsY, LONGREAL))) * magY,
      FLOAT(Math.asinh(FLOAT((p / yc) / epsC, LONGREAL))) * magC,
      FLOAT(Math.asinh(FLOAT((q / yc) / epsC, LONGREAL))) * magC
    }
  END LDWFromRGB;

PROCEDURE RGBFromLDW (READONLY ldw: T): RGB.T =
  VAR p, q, y, yc: REAL;
  BEGIN
    IF ldw = Undefined THEN RETURN RGB.Undefined END;
    y := epsY * FLOAT(Math.sinh(FLOAT(ldw[0] / magY, LONGREAL)));
    yc := FLOAT(Math.hypot(FLOAT(y, LONGREAL), FLOAT(offY, LONGREAL)));
    p := yc * epsC * FLOAT(Math.sinh(FLOAT(ldw[1] / magC, LONGREAL)));
    q := yc * epsC * FLOAT(Math.sinh(FLOAT(ldw[2] / magC, LONGREAL)));
    RETURN T{
      y * Yr + p * Pr + q * Qr,
      y * Yg + p * Pg + q * Qg,
      y * Yb + p * Pb + q * Qb
    }
  END RGBFromLDW;

PROCEDURE Dist (READONLY a, b: T; lumWeight: REAL := 1.0): REAL =
  BEGIN
    IF a = Undefined OR b = Undefined THEN RETURN Intensity.Undefined END;
    RETURN FLOAT(Math.hypot(Math.hypot(
        FLOAT(lumWeight * (a[0] - b[0]), LONGREAL), 
        FLOAT(a[1] - b[1], LONGREAL)
      ), 
      FLOAT(a[2] - b[2], LONGREAL)
    ))
  END Dist;

BEGIN
END LDW.