MODULE PGMPlanetMask EXPORTS Main;
(*
  This program generates a greyscale image that is useful for masking
  and weighting planetary images. The image shows a matte white sphere
  illuminated by a distant point source.  This image is optionally
  multiplied by a user-specified mask.
  
  The sphere's center and apparent radius in the image are either
  given directly by the user, or computed from user-given pixel
  coordinates of three points along its limb.
  
  The direction "sunDir" of the light source user may be specified by the
  user, either directly or by giving two points on the terminator
  (counerclockwise as seen from the lit hemisphere).  If "sunDir" not
  specified, the sphere is assumed to be illuminated frontally,
  i.e. "sunDir=(0,0,1)".
  
  The sphere's apparent relative brightness at a point where the
  outward normal is "N" is computed as
  
    "MAX(0, (light-dark) * (bias + (1-bias)*cos)) + dark"
    
  where "cos" is "dot(sunDir, N)", and "bias" is a user-specified parameter
  that determines the displacement of the termiantor relative to the
  great circle orthogonal to "sunDir". ("bias=0" puts the terminator at the
  great circle, "bias=1" says the whole sphere is illuminated,
  "bias=-1" makes the whole sphere dark except directly under the
  light.)
  
  Pixels that cover some part of the sphere (even its dark side)
  will have values ranging from 1 to the given "maxVal" parameter.  
  Pixels entirely in the background have value 0.
  
  All coordinates are in pixels, with origin at the top left corner, X
  pointing right, Y pointing down, and Z pointing out.
*)

IMPORT PGM, Math, LR3, LR3Extras;
IMPORT Wr, Rd, Text, Thread, Process, FileRd, OSError, ParseParams;
FROM Stdio IMPORT stdin, stdout, stderr;

TYPE
  NAT = CARDINAL;
  LONG = LONGREAL;

TYPE
  Options = RECORD
      NX, NY: NAT;            (* Image size *)
      mkRd: Rd.T;             (* The input mask image file, or NIL *)
      maxVal: NAT;            (* Max pixel value *)
      centerX, centerY: LONG; (* Center of sphere, in screen coords *)
      radius: LONG;           (* Apparent radius of sphere in pixels *)
      sun: LR3.T;             (* Sun direction, screen coords *)
      light: LONG;            (* Illumination intensity *)
      dark: LONG;             (* Dark side light level *)
      bias: LONG;             (* Termiantor displacement bias *)
    END;
    
PROCEDURE Main() =
  <* FATAL Wr.Failure, Thread.Alerted, Rd.Failure *>
  VAR totLight: LONG; 
      nHits: NAT;
      mkPGM: PGM.Reader;   (* Mask pixel reader *)
      mkMax: NAT;          (* Max pixel value in mask *)
      mkStep: LONG;        (* 1/mkMax *)
      mkC, otC: NAT;
  BEGIN
    WITH 
      o = GetOptions(),
      NX = o.NX,
      NY = o.NY,
      otPGM = NEW(PGM.Writer).init(stdout, NX, NY, o.maxVal),
      otRange = FLOAT(o.maxVal-1, LONG)
    DO
      IF o.mkRd # NIL THEN 
        VAR  mkNX, mkNY: NAT;
        BEGIN
          mkPGM := NEW(PGM.Reader).init(o.mkRd, mkNX, mkNY, mkMax);
          <* ASSERT mkNX = NX AND mkNY = NY *>
          mkStep := 1.0d0/FLOAT(mkMax, LONG);
        END
      ELSE
        mkPGM := NIL ; mkMax := 1;
        mkStep := 1.0d0
      END;

      (* Generate image: *)
      FOR iY := 0 TO o.NY-1 DO
        FOR iX := 0 TO o.NX-1 DO
          IF mkPGM # NIL THEN mkPGM.get(mkC) ELSE mkC := mkMax END;
          IF mkC = 0 THEN
            otC := 0
          ELSE
            WITH mkF = mkStep * FLOAT(mkC, LONG) DO
              nHits := 0; totLight := 0.0d0;
              ColorSample(iX, iY,  +0.25d0, +0.25d0, o, nHits, totLight);
              ColorSample(iX, iY,  -0.25d0, +0.25d0, o, nHits, totLight);
              ColorSample(iX, iY,  -0.25d0, -0.25d0, o, nHits, totLight);
              ColorSample(iX, iY,  +0.25d0, -0.25d0, o, nHits, totLight);
              IF nHits = 0 THEN
                otC := 0
              ELSE
                otC := 1 + ROUND(otRange * mkF * totLight/4.0d0)
              END
            END
          END;
          otPGM.put(otC)
        END
      END;
      otPGM.finish()
    END
  END Main;
  
PROCEDURE ColorSample(
    iX, iY: NAT; 
    dX, dY: LONG; 
    READONLY o: Options;
    VAR nHits: NAT; 
    VAR totLight: LONG;
  ) =
  (*
    Lighting calculations for the image point "(iX+dX, iY+dY)".  
    Increments "totLight" by the point's color.  Also
    increments "nHits" if the points is on the sphere. *)
  BEGIN
    WITH
      x = FLOAT(iX, LONG) + dX - o.centerX, 
      y = FLOAT(iY, LONG) + dY - o.centerY,
      x2y2 = x*x + y*y,
      r2 = o.radius * o.radius
    DO
      IF x2y2 <= 0.999999d0 * r2 THEN
        INC(nHits);
        WITH 
          z = Math.sqrt(r2 - x2y2),
          cos = (o.sun[0] * x + o.sun[1] * y + o.sun[2] * z)/o.radius,
          L = MAX(0.0d0, (o.light-o.dark)*(o.bias + (1.0d0 - o.bias)*cos)) + o.dark
        DO
          <* ASSERT cos <= +1.0d0 AND cos >= -1.0d0 *>
          totLight := totLight + L
        END
      END;
    END
  END ColorSample;

PROCEDURE FindCircle(
    READONLY x, y: ARRAY [0..2] OF LONG; 
    VAR centerX, centerY, radius: LONG;
  ) =
  BEGIN
    WITH
      s0 = x[0]*x[0] + y[0]*y[0],
      s1 = x[1]*x[1] + y[1]*y[1],
      s2 = x[2]*x[2] + y[2]*y[2],
      
      A = + LR3Extras.Det(
        LR3.T{x[0], y[0], 1.0d0},
        LR3.T{x[1], y[1], 1.0d0},
        LR3.T{x[2], y[2], 1.0d0}
      ),
      
      D = - LR3Extras.Det(
        LR3.T{s0, y[0], 1.0d0},
        LR3.T{s1, y[1], 1.0d0},
        LR3.T{s2, y[2], 1.0d0}
      ),
      E = + LR3Extras.Det(
        LR3.T{s0, x[0], 1.0d0},
        LR3.T{s1, x[1], 1.0d0},
        LR3.T{s2, x[2], 1.0d0}
      ),
      F = - LR3Extras.Det(
        LR3.T{s0, x[0], y[0]},
        LR3.T{s1, x[1], y[1]},
        LR3.T{s2, x[2], y[2]}
      )
    DO
      centerX := - 0.5d0 * D/A;
      centerY := - 0.5d0 * E/A;
      radius := Math.sqrt(MAX(0.0d0, centerX*centerX + centerY*centerY - F/A))
    END
  END FindCircle;

PROCEDURE FindSun(
    READONLY x, y: ARRAY [0..1] OF LONG; 
    centerX, centerY, radius: LONG;
    VAR sun: LR3.T;
  ) =
  (*
    Computes a sun direction vector given the 
    sphere center and radius, and two points "(x[i],y[i])" on the terminator.
  *)
  BEGIN
    WITH
      r2 = radius * radius,
      
      uX = x[0] - centerX,
      uY = y[0] - centerY,
      uZ = Math.sqrt(MAX(0.0d0, r2 - uX*uX - uY*uY)),
      
      vX = x[1] - centerX,
      vY = y[1] - centerY,
      vZ = Math.sqrt(MAX(0.0d0, r2 - vX*vX - vY*vY))
    DO
      (* Note that the screen coordinate system is left-handed... *)
      sun := LR3.Dir(LR3Extras.Cross(LR3.T{vX, vY, vZ}, LR3.T{uX, uY, uZ}))
    END
  END FindSun;

PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>

  CONST
    MaxSize: NAT = 4096;
    MaxMaxVal: NAT = 256*256-1;
    MaxCoord: LONG = 10000.0d0;

  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO

      TRY
        pp.getKeyword("-size");
        o.NX := pp.getNextInt(1, MaxSize);
        o.NY := pp.getNextInt(1, MaxSize);
        
        IF pp.keywordPresent("-maxVal") THEN
          o.maxVal := pp.getNextInt(1, MaxMaxVal)
        ELSE
          o.maxVal := 255
        END;
        
        WITH
          minX = -0.5d0 - MaxCoord,
          maxX = FLOAT(o.NX, LONG) - 0.5d0 + MaxCoord,
          minY = -0.5d0 - MaxCoord,
          maxY = FLOAT(o.NY, LONG) - 0.5d0 + MaxCoord,
          maxR = MaxCoord + FLOAT(MAX(o.NX, o.NY), LONG)
        DO
          IF pp.keywordPresent("-limb") THEN
            VAR x, y: ARRAY [0..2] OF LONG;
            BEGIN
              FOR i := 0 TO 2 DO 
                x[i]:= pp.getNextLongReal(minX, maxX);
                y[i]:= pp.getNextLongReal(minY, maxY)
              END;
              FindCircle(x, y, o.centerX, o.centerY, o.radius)
            END;
          ELSIF pp.keywordPresent("-disk") THEN
            o.centerX := pp.getNextLongReal(minX, maxX);
            o.centerY := pp.getNextLongReal(minY, maxY);
            o.radius := pp.getNextLongReal(0.0d0, maxR)
          ELSE
            pp.error("please specify either \"-limb\" or \"-disk\"")
          END
        END;
        
        IF pp.keywordPresent("-sun") THEN
          FOR i := 0 TO 2 DO o.sun[i] := pp.getNextLongReal(-1.0d0, 1.0d0) END;
          o.sun := LR3.Dir(o.sun);
        ELSIF pp.keywordPresent("-terminator") THEN
          VAR x, y: ARRAY [0..1] OF LONG;
          BEGIN
            FOR i := 0 TO 1 DO 
              x[i]:= pp.getNextLongReal(-0.5d0, FLOAT(o.NX, LONG)-0.5d0);
              y[i]:= pp.getNextLongReal(-0.5d0, FLOAT(o.NY, LONG)-0.5d0)
            END;
            FindSun(x, y, o.centerX, o.centerY, o.radius, o.sun)
          END;
        ELSE
          o.sun := LR3.T{0.0d0, 0.0d0, 1.0d0}
        END;
        
        IF pp.keywordPresent("-light") THEN
          WITH 
            maxLit = pp.getNextInt(1, o.maxVal)
          DO
            o.light := FLOAT(maxLit-1, LONG)/FLOAT(MAX(1, o.maxVal-1), LONG)
          END
        ELSE
          o.light := 1.0d0
        END;
        
        IF pp.keywordPresent("-dark") THEN
          WITH 
            maxLit = 1 + ROUND(o.light * FLOAT(o.maxVal-1, LONG)),
            minLit = pp.getNextInt(1, maxLit)
          DO
            o.dark := FLOAT(minLit-1, LONG)/FLOAT(MAX(1, o.maxVal-1), LONG)
          END
        ELSE
          o.dark := 1.0d0
        END;
        
        IF pp.keywordPresent("-bias") THEN
          o.bias := pp.getNextLongReal(-1.0d0, 1.0d0)
        ELSE
          o.bias := 0.0d0
        END;
        
        IF pp.keywordPresent("-mask") THEN
          o.mkRd := GetFileOption(pp);
        ELSE
          o.mkRd := NIL;
        END;
        

        pp.skipParsed();
        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: PGMPlanetMask \\\n");
          Wr.PutText(stderr, "  -size <NX> <NY> [ -maxVal <M> ] \\\n");
          Wr.PutText(stderr, "  [ -mask <pgmfile> ] \\\n");
          Wr.PutText(stderr, "  [ -limb <X1> <Y1> <X2> <Y2> <X3> <Y3> | -disk <cX> <cY> <R> ] \\\n");
          Wr.PutText(stderr, "  [ -sun <X> <Y> <Z> | -terminator <X1> <Y1> <X2> <Y2> ] \\\n");
          Wr.PutText(stderr, "  [ -light <L> ] [-dark <K>] [ -bias <B> ]\n");
          Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;
  
PROCEDURE GetFileOption(pp: ParseParams.T): Rd.T RAISES {ParseParams.Error} =
  BEGIN
    WITH name = pp.getNext() DO
      IF Text.Equal(name, "-") THEN
        RETURN stdin
      ELSE
        TRY
          RETURN FileRd.Open(name)
        EXCEPT
        | OSError.E => 
            pp.error("can't open file \"" & name & "\"");
            RETURN NIL
        END
      END
    END
  END GetFileOption;

BEGIN 
  Main();
END PGMPlanetMask.