(* Last edited on 2000-01-05 04:16:43 by stolfi *)
MODULE PGMRemoveGhost EXPORTS Main;

(*
  "Ghost" (bleedthrough) remover for scanned transparent images.
  
  Reads grayscale images "x0,x1,y0,y1" of both sides ("x" = recto, "y"
  = verso) of a translucent folio, scanned with black (0) and white
  (1) backgounds, respectively. The images are assumed to be in
  precise alignment; in particular, the verso images must have been
  flipped left-to-right.
  
  Ouputs images "xc,yc" of the recto and verso, with bleedthrough 
  removed.
*)    

IMPORT Rd, FileRd, Wr, FileWr, Fmt, Thread, OSError, Process, ParseParams;
IMPORT PGM, PGMImageFilter, IA;
FROM Stdio IMPORT stderr;
FROM Math IMPORT sqrt;

<* FATAL Wr.Failure, Thread.Alerted *>

TYPE
  LONG = LONGREAL;
  NAT = CARDINAL;
  BOOL = BOOLEAN;
  INT = INTEGER;
  
  Interval = IA.T;

  Options = RECORD
      inName: TEXT;          (* Input image file name prefix. *)
      outName: TEXT;         (* Output image file name prefix. *)
      showError: BOOL;       (* TRUE outputs error instead of midvalue. *)
      gamma: LONG;           (* Assumes "intensity = (pixelValue/255)^gamma". *)
      paperColor: LONG;      (* Albedo of thick paper. *)
      blackColor: LONG;      (* Actual albedo of "black" background. *)
      whiteColor: LONG;      (* Actual albedo of "white" background. *)
      writePixels: BOOL;     (* TRUE writes pixel triplets to "outName.pxt". *)
      debugX, debugY: NAT;   (* Coordinates of pixel to debug. *)
    END;
    
PROCEDURE Main() =
  <* FATAL Wr.Failure, Thread.Alerted *>
  VAR 
    NX, NY: NAT := LAST(CARDINAL);
    iMax: NAT := LAST(CARDINAL);
    nConflicts: NAT := 0; (* Number of pixels with inconsistens eqs. *)
  CONST
    ZeroOne = Interval{0.0d0, 1.0d0};  (* Normal range of pixel values. *)
    MaxPlotPixels = 10000;
  BEGIN
    WITH 
      o = GetOptions(),

      xBlackRd = OpenReader(o.inName & "-x0", NX, NY, iMax),
      xWhiteRd = OpenReader(o.inName & "-x1", NX, NY, iMax),
      yBlackRd = OpenReader(o.inName & "-y0", NX, NY, iMax),
      yWhiteRd = OpenReader(o.inName & "-y1", NX, NY, iMax),
      
      NXY = NX*NY,
      plotStep = (NXY - 1) DIV MaxPlotPixels + 1,

      xcWr = OpenWriter(o.outName & "-xc", NX, NY, iMax),
      ycWr = OpenWriter(o.outName & "-yc", NX, NY, iMax),
      trWr = OpenWriter(o.outName & "-tr", NX, NY, iMax),

      iRow = NEW(REF ARRAY OF NAT, NX)^,
      
      pixelValue = MakeGammaTable(iMax, o.gamma)^,
      
      xBlackRow = NEW(REF ARRAY OF Interval, NX)^,
      xWhiteRow = NEW(REF ARRAY OF Interval, NX)^,
      yBlackRow = NEW(REF ARRAY OF Interval, NX)^,
      yWhiteRow = NEW(REF ARRAY OF Interval, NX)^,
        (* 
          Invariant: while processing scanline "yp", 
          the corresponding input scanlines are stored in
          these variables. *)
      
      (* For image output: *)
      xcRow = NEW(REF ARRAY OF Interval, NX)^,
      ycRow = NEW(REF ARRAY OF Interval, NX)^,
      trRow = NEW(REF ARRAY OF Interval, NX)^,
      
      oRow = NEW(REF ARRAY OF NAT, NX)^
    DO
      Wr.PutText(stderr, "images: " & Fmt.Int(NX) & " × " & Fmt.Int(NY) & " pixels");
      Wr.PutText(stderr, ", ");
      Wr.PutText(stderr, "maxval = " & Fmt.Int(iMax) & "\n");
      
      PROCEDURE GetRow(rd: PGM.Reader; VAR rRow: ARRAY OF Interval) =
        BEGIN
          PGMImageFilter.ReadPixels(rd, iRow, NIL, LAST(NAT));
          FOR xp := 0 TO NX-1 DO rRow[xp] := pixelValue[iRow[xp]] END;
        END GetRow;

      PROCEDURE PutRow(wr: PGM.Writer; READONLY rRow: ARRAY OF Interval) =
        BEGIN
          IF o.showError THEN
            PGMImageFilter.FixIntervalPixelErrors(rRow, iMax, oRow, LAST(NAT))
          ELSE
            PGMImageFilter.FixIntervalPixels(rRow, iMax, 0.499999d0, oRow, LAST(NAT));
            PGMImageFilter.ReplacePixelValue(LAST(NAT), iMax DIV 2, oRow);
            PGMImageFilter.WritePixels(wr, oRow, mkWr := NIL, noPixel := LAST(NAT))
          END;
        END PutRow;

      VAR pxWr: Wr.T; 
          x2R, y2R, S2: Interval;
      BEGIN
        IF o.writePixels THEN pxWr := OpenPlot(o.outName, o) END;
    
        (* Perform filtering *)
        FOR yp := 0 TO NY-1 DO
          GetRow(xBlackRd, xBlackRow); GetRow(xWhiteRd, xWhiteRow);
          GetRow(yBlackRd, yBlackRow); GetRow(yWhiteRd, yWhiteRow);
          (* Process pixels: *)
          FOR xp := 0 TO NX-1 DO
            WITH
              xyp = NX*yp + xp,
              xBlack = xBlackRow[xp], xWhite = xWhiteRow[xp], xc = xcRow[xp],
              yBlack = yBlackRow[xp], yWhite = yWhiteRow[xp], yc = ycRow[xp],
              tr = trRow[xp],
              xDelta = IA.Sub(xWhite, xBlack),
              yDelta = IA.Sub(yWhite, yBlack),
              consistent = SolveOldEqs(
                xBlack, yBlack, 
                xDelta, yDelta, 
                o.whiteColor, o.blackColor,
                (*OUT*) x2R, y2R, S2,
                tol := 0.001d0,
                debug := (xp = o.debugX AND yp = o.debugY)
              )
            DO
              IF NOT consistent THEN INC(nConflicts) END;
              WITH
                (* Compute "xTy = x*T*y" where T is transmittance: *)
                xTy = IA.Sqrt(S2),
                (* Compute the product "xy = x*y" of the two clean images: *)
                xRy = IA.Sqrt(IA.Mul(x2R, y2R)),
                xy = IA.Clip(IA.Add(xTy, IA.Scale(1.0d0/o.paperColor, xRy)), ZeroOne),
                (* Estimate the reflectance "R" and the transmittance "T": *)
                RRange = Interval{MAX(x2R.lo, y2R.lo), 1.0d0},
                R = IA.Clip(IA.Div(xRy, xy), RRange),
                TRange = Interval{0.0d0, MAX(0.0d0, 1.0d0 - R.hi/o.paperColor)},
                T = IA.Clip(IA.Div(xTy, xy), TRange),
                (* Compute the clean images "xClean = x" and "yClean = y": *)
                xRange = Interval{MAX(sqrt(x2R.lo), MAX(xRy.lo, xTy.lo)), 1.0d0},
                xClean = IA.Div(x2R, R),
                yRange = Interval{MAX(sqrt(y2R.lo), MAX(xRy.lo, xTy.lo)), 1.0d0},
                yClean = IA.Div(y2R, R)
              DO
                xc := IA.Clip(xClean, xRange);
                yc := IA.Clip(yClean, yRange);
                tr := IA.Clip(T, ZeroOne);
                IF o.writePixels AND (xyp + plotStep DIV 2) MOD plotStep = 0 THEN
                  PlotPixel(pxWr,
                    xp, yp, 
                    xBlack, yBlack, 
                    xWhite, yWhite, 
                    x2R, y2R, S2,
                    R, T,
                    xClean, yClean
                  )
                END;
                IF xp = o.debugX AND yp = o.debugY THEN
                  DebugPoint("xp,yp", xp, yp);
                  DebugInterval("xBlack", xBlack);
                  DebugInterval("xWhite", xWhite);
                  DebugInterval("yBlack", yBlack);
                  DebugInterval("yWhite", yWhite);
                  DebugInterval("x2R", x2R);
                  DebugInterval("y2R", y2R);
                  DebugInterval("S2", S2);
                  DebugInterval("xTy", xTy);
                  DebugInterval("xRy", xRy);
                  DebugInterval("xy", xy);
                  DebugInterval("R", R);
                  DebugInterval("T", T);
                  DebugInterval("xClean", xClean);
                  DebugInterval("yClean", yClean);
                END
              END
            END
          END;
          PutRow(xcWr, xcRow); 
          PutRow(ycWr, ycRow);
          PutRow(trWr, trRow);
        END;
        xcWr.finish();
        ycWr.finish();
        IF o.writePixels THEN Wr.Close(pxWr) END;
        IF nConflicts > 0 THEN 
          Wr.PutText(stderr, "**Inconstent solutions in " & FI(nConflicts,1) & " pixels\n");
        END
      END
    END
  END Main;
 
PROCEDURE CombineIntervals(a, b, meet: Interval): Interval =
  BEGIN
    IF meet.lo <= meet.hi THEN
      RETURN meet
    ELSE
      WITH
        hi = MIN(a.hi, b.hi),
        wd = 0.5d0 * MIN(a.hi-a.lo, b.hi-b.lo),
        lo = MAX(0.0d0, hi - wd)
      DO
        RETURN Interval{lo, hi}
      END
    END
  END CombineIntervals;

PROCEDURE SolveOldEqs(
    READONLY xB, yB: Interval; 
    READONLY xDelta, yDelta: Interval;
    a, b: LONG;
    VAR (*OUT*) x2R, y2R: Interval;
    VAR (*OUT*) S2: Interval;
    tol: LONG;
    debug: BOOL;
  ): BOOL =
  (*
    Solves the equations
    |   x2R = xB - S2*b/(1 - b*y2R)
    |   y2R = yB - S2*b/(1 - b*x2R)
    |   S2 = xDelta/(a/(1-a*y2R) - b/(1-b*y2R))
    |   S2 = yDelta/(a/(1-a*x2R) - b/(1-b*x2R))
    for "x2R", "y2R", and "S2".
    
    The given data are "xB,yB", the black background images (recto and
    verso), and "xDelta,yDelta", the differences between white and
    black background images.
    
    The unknowns are defined in terms of the sandwich model
    parameters: the ink transmittances "x" and "y", the the paper's
    reflectance "R", and its transmittance "T".
    
    The procedure actually solves for the auxiliary variables "x2R =
    x^2*R", "y2R = y^2*R", and "S2 = x^2*T^2*y^2". Note that there 
    are two equations defining "S2", i.e. the system is overdetermined.
    We compute both values "S2_x" and "S2_y", and take some average.
    Returns TRUE if in the last iteration the intervals for "S2_x" and "S2_y"
    had non-empty intersection, FALSE if they were disjoint.
  *)
  VAR x2RPrev, y2RPrev: Interval;
      S2_x, S2_y, S2_meet: Interval;
  CONST 
    EpsOne = IA.T{1.0d-5, 1.0d0};
    OneOne = IA.T{1.0d0, 1.0d0};
    ZeroOne = IA.T{0.0d0, 1.0d0};
  BEGIN
    x2R := IA.T{0.0d0, xB.hi};
    y2R := IA.T{0.0d0, yB.hi};
    IF debug THEN
      Wr.PutText(stderr, "--- BEGIN SolveOldEqs ---\n");
      DebugInterval("x2R", x2R);
      DebugInterval("y2R", y2R);
    END;
    LOOP
      WITH
        aYDen = IA.Clip(IA.Sub(OneOne, IA.Scale(a, y2R)), EpsOne),
        aYFrac = IA.Scale(a, IA.Inv(aYDen)),

        bYDen = IA.Clip(IA.Sub(OneOne, IA.Scale(b, y2R)), EpsOne),
        bYFrac = IA.Scale(b, IA.Inv(bYDen)),

        xDen = IA.Sub(aYFrac, bYFrac),

        aXDen = IA.Clip(IA.Sub(OneOne, IA.Scale(a, x2R)), EpsOne),
        aXFrac = IA.Scale(a, IA.Inv(aXDen)),

        bXDen = IA.Clip(IA.Sub(OneOne, IA.Scale(b, x2R)), EpsOne),
        bXFrac = IA.Scale(b, IA.Inv(bXDen)),

        yDen = IA.Sub(aXFrac, bXFrac)
      DO
        S2_x := IA.Clip(IA.Div(xDelta, xDen), ZeroOne);
        S2_y := IA.Clip(IA.Div(yDelta, yDen), ZeroOne);
        S2_meet := IA.Meet(S2_x, S2_y);
        S2 := CombineIntervals(S2_x, S2_y, S2_meet);
        x2RPrev := x2R; y2RPrev := y2R;
        x2R := IA.Clip(IA.Sub(xB, IA.Mul(bYFrac, S2)), x2RPrev);
        y2R := IA.Clip(IA.Sub(yB, IA.Mul(bXFrac, S2)), y2RPrev);
        IF debug THEN
          Wr.PutText(stderr, "\n");
          DebugInterval("xDen", xDen);
          DebugInterval("yDen", yDen);
          DebugInterval("S2_x", S2_x);
          DebugInterval("S2_y", S2_y);
          DebugInterval("S2", S2);
          DebugInterval("x2R", x2R);
          DebugInterval("y2R", y2R);
        END;
        <* ASSERT x2R.lo <= x2R.hi *>
        <* ASSERT y2R.lo <= y2R.hi *>
        IF  ABS(x2R.lo - x2RPrev.lo) + ABS(x2R.hi - x2RPrev.hi) < tol
        AND ABS(y2R.lo - y2RPrev.lo) + ABS(y2R.hi - y2RPrev.hi) < tol THEN
          EXIT
        END
      END
    END;
    IF debug THEN
      Wr.PutText(stderr, "--- END SolveOldEqs ---\n");
    END;
    RETURN S2_meet.lo <= S2_meet.hi
  END SolveOldEqs;
  
PROCEDURE OpenReader(name: TEXT; VAR NX, NY: NAT; VAR iMax: NAT): PGM.Reader =
  (*
    Opens a PGM image reader. If "iMax" is defined, checks 
    whether the values of "NX", "NY", "iMax" in the file
    are consistent with the values in the variables; else
    sets the latter from the image file.
  *)
  <* FATAL Rd.Failure *>
  VAR NXT, NYT, iMaxT: NAT; 
  BEGIN
    TRY 
      WITH 
        rd = FileRd.Open(name & ".pgm"),
        rdPGM = NEW(PGM.Reader).init(rd, NXT, NYT, iMaxT)
      DO
        IF iMax = LAST(CARDINAL) THEN
          NX := NXT; NY := NYT; iMax := iMaxT 
        ELSE
          <* ASSERT NX = NXT *> 
          <* ASSERT NY = NYT *> 
          <* ASSERT iMax = iMaxT *> 
        END;
        RETURN rdPGM
      END
    EXCEPT
    | OSError.E => 
        Wr.PutText(stderr, "can't open input file \"" & name & "\"");
        Process.Exit(1);
        RETURN NIL
    END;
  END OpenReader;
  
PROCEDURE OpenWriter(name: TEXT; NX, NY: NAT; iMax: NAT): PGM.Writer =
  BEGIN
    TRY 
      WITH 
        wr = FileWr.Open(name & ".pgm"),
        wrPGM = NEW(PGM.Writer).init(wr, NX, NY, iMax)
      DO
        RETURN wrPGM
      END
    EXCEPT
    | OSError.E => 
        Wr.PutText(stderr, "can't open output file \"" & name & "\"");
        Process.Exit(1);
        RETURN NIL
    END
  END OpenWriter;
  
PROCEDURE MakeGammaTable(iMax: NAT; gamma: LONG): REF ARRAY OF Interval =
  BEGIN
    WITH
      rp = NEW(REF ARRAY OF NAT, iMax+1), tp = rp^,
      rv = NEW(REF ARRAY OF Interval, iMax+1), tv = rv^
    DO
      FOR i := 0 TO iMax DO tp[i] := i END;
      PGMImageFilter.FloatPixels(tp, iMax, tv, LAST(NAT));
      FOR i := 0 TO iMax DO tv[i] := IA.Pow(tv[i], gamma) END;
      RETURN rv
    END
  END MakeGammaTable;
  
PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO
      TRY
        pp.getKeyword("-inName");
        o.inName := pp.getNext();
        
        pp.getKeyword("-outName");
        o.outName := pp.getNext();
        
        IF pp.keywordPresent("-gamma") THEN
          o.gamma := pp.getNextLongReal(0.01d0, 100.0d0);
        ELSE
          o.gamma := 1.0d0;
        END;
        
        pp.getKeyword("-paperColor");
        o.paperColor := pp.getNextLongReal(0.0d0, 1.0d0);
        
        IF pp.keywordPresent("-blackColor") THEN
          o.blackColor := pp.getNextLongReal(0.0d0, 1.0d0);
        ELSE
          o.blackColor := 0.0d0;
        END;
        
        IF pp.keywordPresent("-whiteColor") THEN
          o.whiteColor := pp.getNextLongReal(0.0d0, 1.0d0);
        ELSE
          o.whiteColor := 1.0d0;
        END;
        
        o.showError := pp.keywordPresent("-showError");
        
        o.writePixels := pp.keywordPresent("-writePixels");
        
        IF pp.keywordPresent("-debug") THEN
          o.debugX := pp.getNextInt(0, 10000);
          o.debugY := pp.getNextInt(0, 10000);
        ELSE
          o.debugX := LAST(NAT); o.debugY := LAST(NAT);
        END;
        
        pp.skipParsed();
        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: PGMRemoveGhost \\\n");
          Wr.PutText(stderr, "  -paperColor NUM [ -gamma NUM ] \\\n");
          Wr.PutText(stderr, "  [ -blackColor NUM ] [ -whiteColor NUM ] \\\n");
          Wr.PutText(stderr, "  [ -writePixels ] [ -debug X Y ] \\\n");
          Wr.PutText(stderr, "  -inName NAME \\\n");
          Wr.PutText(stderr, "  -outName NAME \n");
          Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;

PROCEDURE OpenPlot(name: TEXT; READONLY o: Options): Wr.T =
  <* FATAL OSError.E *>
  BEGIN
    WITH 
      wr = FileWr.Open(name & ".ppt")
    DO
      Wr.PutText(wr, "# blackColor = " & FLR(o.blackColor, 6,4) & "\n");
      Wr.PutText(wr, "# whiteColor = " & FLR(o.whiteColor, 6,4) & "\n");
      Wr.PutText(wr, "# paperColor = " & FLR(o.paperColor, 6,4) & "\n");
      Wr.PutText(wr, "#   xp    yp  xBlack yBlack  xwhite yWhite     x2R    y2R     S2       R      T  xClean xError  yClean yError\n");
      Wr.PutText(wr, "# ----  ----  ------ ------  ------ ------  ------ ------ ------  ------ ------  ------ ------  ------ ------\n");
      RETURN wr
    END
  END OpenPlot;
  
PROCEDURE PlotPixel(
    wr: Wr.T;
    xp, yp: NAT;
    READONLY xBlack, yBlack: Interval;
    READONLY xWhite, yWhite: Interval;
    READONLY x2R, y2R, S2: Interval; 
    READONLY R: Interval;
    READONLY T: Interval;
    READONLY xClean, yClean: Interval;
  ) =
  
  BEGIN
    Wr.PutChar(wr, ' ');
    Wr.PutText(wr, FI(xp, 5));                 (* 01 *)  
    Wr.PutChar(wr, ' ');  
    Wr.PutText(wr, FI(yp, 5));                 (* 02 *)  
    Wr.PutChar(wr, ' ');  
    Wr.PutChar(wr, ' ');  
    Wr.PutText(wr, FLR(IA.Mid(xBlack), 6,4));  (* 03 *)  
    Wr.PutChar(wr, ' ');  
    Wr.PutText(wr, FLR(IA.Mid(yBlack), 6,4));  (* 04 *)  
    Wr.PutChar(wr, ' ');  
    Wr.PutChar(wr, ' ');  
    Wr.PutText(wr, FLR(IA.Mid(xWhite), 6,4));  (* 05 *)  
    Wr.PutChar(wr, ' ');  
    Wr.PutText(wr, FLR(IA.Mid(yWhite), 6,4));  (* 06 *)  
    Wr.PutChar(wr, ' ');  
    Wr.PutChar(wr, ' ');  
    Wr.PutText(wr, FLR(IA.Mid(x2R), 6,4));     (* 07 *)  
    Wr.PutChar(wr, ' ');  
    Wr.PutText(wr, FLR(IA.Mid(y2R), 6,4));     (* 08 *)  
    Wr.PutChar(wr, ' ');  
    Wr.PutText(wr, FLR(IA.Mid(S2),  6,4));     (* 09 *)
    Wr.PutChar(wr, ' ');                               
    Wr.PutChar(wr, ' ');                       
    Wr.PutText(wr, FLR(IA.Mid(R), 6,4));       (* 10 *)
    Wr.PutChar(wr, ' ');                               
    Wr.PutText(wr, FLR(IA.Mid(T), 6,4));       (* 11 *)
    Wr.PutChar(wr, ' ');                               
    Wr.PutChar(wr, ' ');                               
    Wr.PutText(wr, FLR(IA.Mid(xClean), 6,4));  (* 12 *)
    Wr.PutChar(wr, ' ');                               
    Wr.PutText(wr, FLR(IA.Rad(xClean), 6,4));  (* 13 *)
    Wr.PutChar(wr, ' ');                               
    Wr.PutChar(wr, ' ');                               
    Wr.PutText(wr, FLR(IA.Mid(yClean), 6,4));  (* 14 *)
    Wr.PutChar(wr, ' ');                               
    Wr.PutText(wr, FLR(IA.Rad(yClean), 6,4));  (* 15 *)
    Wr.PutChar(wr, '\n');
  END PlotPixel;

PROCEDURE FI(x: INT; w: CARDINAL): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Int(x), w)
  END FI;

PROCEDURE FLR(x: LONGREAL; w, p: CARDINAL): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, p), w)
  END FLR;

PROCEDURE DebugPoint(msg: TEXT; xp, yp: INT) =
  BEGIN
    Wr.PutText(stderr, Fmt.Pad(msg, 25) & " = ");
    Wr.PutText(stderr, Fmt.Int(xp));
    Wr.PutText(stderr, "  ");
    Wr.PutText(stderr, Fmt.Int(yp));
    Wr.PutText(stderr, "\n")
  END DebugPoint;

PROCEDURE DebugInterval(msg: TEXT; READONLY z: Interval) = 
  BEGIN
    Wr.PutText(stderr, Fmt.Pad(msg, 25) & " = ");
    Wr.PutText(stderr, FLR(z.lo, 9, 6));
    Wr.PutText(stderr, " _ ");
    Wr.PutText(stderr, FLR(z.hi, 9, 6));
    Wr.PutText(stderr, "\n")
  END DebugInterval;

<*UNUSED*>
PROCEDURE DebugRow(READONLY z: ARRAY OF Interval) = 
  BEGIN
    Wr.PutText(stderr, "\n");
    FOR x := 0 TO LAST(z) DO
      Wr.PutText(stderr, FLR(z[x].lo, 9, 6));
      Wr.PutText(stderr, " _ ");
      Wr.PutText(stderr, FLR(z[x].hi, 9, 6));
      Wr.PutText(stderr, "\n")
    END;
    Wr.PutText(stderr, "\n")
  END DebugRow;

<*UNUSED*>
PROCEDURE FIA(x: Interval; w, p: CARDINAL): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(IA.Mid(x), Fmt.Style.Fix, p), w)
  END FIA;

BEGIN 
  Main();
END PGMRemoveGhost.