(* Last edited on 1999-03-18 01:32:31 by stolfi *)

MODULE PGMConvol EXPORTS Main;

IMPORT PGM, IA, Math;
FROM IA IMPORT Infinity;
IMPORT Text, Fmt, Rd, FileRd, Wr, Thread, OSError, Process, ParseParams;
FROM Stdio IMPORT stdin, stdout, stderr;

CONST
  MaxMaxVal = 65535;

TYPE
  Options = RECORD
      rdA: Rd.T;         (* The convolution mask *)
      rangeA: Interval;  (* The nominal range of image "A". *)
      rdB: Rd.T;         (* The input image *)
      rangeB: Interval;  (* The nominal range OF image "B". *)
      maxVal: CARDINAL;  (* Max output pixel, or 0 if automatic *)
      quiet: BOOLEAN;    (* TRUE supresses informative messages *)
    END;
    
  LONG = LONGREAL;
  NAT = CARDINAL;
  Interval = IA.T;
    
<* FATAL Wr.Failure, Thread.Alerted, Rd.Failure *>

PROCEDURE Main() =
  VAR NXA, NYA, MaxA, NXB, NYB, MaxB: NAT;
  BEGIN
    WITH 
      o = GetOptions(),
      pgmA = NEW(PGM.Reader).init(o.rdA, NXA, NYA, MaxA),
      pgmB = NEW(PGM.Reader).init(o.rdB, NXB, NYB, MaxB)
    DO
      IF NYA <= NYB THEN 
        DoConvol(o, pgmA, NXA, NYA, MaxA, pgmB, NXB, NYB, MaxB)
      ELSE
        DoConvol(o, pgmB, NXB, NYB, MaxB, pgmA, NXA, NYA, MaxA)
      END
    END
  END Main;
  
PROCEDURE DoConvol(
    READONLY o: Options;
    pgmA: PGM.Reader; NXA, NYA: NAT; MaxA: NAT;
    pgmB: PGM.Reader; NXB, NYB: NAT; MaxB: NAT;
  ) =
  BEGIN
    <* ASSERT NYA <= NYB *>
    WITH
      NXC = NXA + NXB - 1,
      NYC = NYA + NYB - 1,
      MaxC = ComputeMaxC(o.maxVal, MaxA, MaxB, NP := MIN(NXA,NXB)*MIN(NYA,NYB)),
      pgmC = NEW(PGM.Writer).init(stdout, NYC, NXC, MaxC),

      imgA = NEW(REF ARRAY OF ARRAY OF LONG, NYA, NXA)^,

      (* We need only "NYA" lines of the "B" image at a time: *)
      imgB = NEW(REF ARRAY OF ARRAY OF LONG, NYA, NXB)^
    DO
      VAR scaleA := (o.rangeA.hi - o.rangeA.lo)/FLOAT(MaxA, LONG);
          scaleB := (o.rangeB.hi - o.rangeB.lo)/FLOAT(MaxB, LONG);
          scaleC: LONG;
          sum: LONG;
          c: NAT;
          rangeC: Interval := Interval{+Infinity, -Infinity};
      BEGIN
        IF NOT o.quiet THEN
          Wr.PutText(stderr, "PGMConvol: output maxVal = " & Fmt.Int(MaxC) & "\n")
        END;
        
        (* Read mask, normalize to "rangeA", *)
        (*   and compute min/max convolution value: *)
        FOR y := 0 TO NYA-1 DO
          FOR x := 0 TO NXA-1 DO
            pgmA.get(c); 
            WITH 
              v = o.rangeA.lo + scaleA * FLOAT(c, LONG),
              rangeP = IA.Scale(v, o.rangeB)
            DO
              rangeC := IA.Add(rangeC, rangeP);
              imgA[y,x] := v
            END;
          END
        END;
        
        IF rangeC.lo = rangeC.hi THEN 
          (* Mask is all zeros; write gray image out *)
          WITH c = MaxC DIV 2 DO
            FOR y := 0 TO NYC-1 DO
              FOR x := 0 TO NXC-1 DO
                pgmC.put(c)
              END
            END
          END
        ELSE
          (* Perform convolution *)
          (* Invariant: scanline "yb" of input image, *)
          (* if available, is stored in "imgB[yb MOD NYA]" *)
          
          WITH eps = 1.0d-6 * IA.Rad(rangeC) DO
            rangeC := IA.Add(rangeC, Interval{-eps, +eps})
          END;
          scaleC := FLOAT(MaxC, LONG)/(rangeC.hi - rangeC.lo);
          
          FOR yc := 0 TO NYC-1 DO
            (* Compute line "yc" of convolution *)
            
            (* Read line "yc" of image "B", if it exists: *)
            WITH ym = yc MOD NYA DO
              IF yc < NYB THEN
                FOR x := 0 TO NXB-1 DO 
                  pgmB.get(c); 
                  imgB[ym,x] := o.rangeB.lo + scaleB * FLOAT(c, LONG)
                END
              END
            END;

            WITH
              yaLo = MAX(0, NYA - 1 - yc),
              yaHi = MIN(NYA - 1,  NYC - 1 - yc)
            DO
              FOR xc := 0 TO NXC-1 DO
                (* Compute dot product: *)
                sum := 0.0d0;
                WITH 
                  xaLo = MAX(0, NXA - 1 - xc),
                  xaHi = MIN(NXA - 1,  NXC - 1 - xc)
                DO
                  FOR ya := yaLo TO yaHi DO
                    WITH yb = ya + yc - (NYA - 1), ybm = yb MOD NYA DO
                      FOR xa := xaLo TO xaHi DO
                        WITH xb = xa + xc - (NXA - 1) DO
                          sum := sum + imgA[ya,xa] * imgB[ybm,xb]
                        END
                      END
                    END
                  END
                END;
                <* ASSERT sum >= rangeC.lo AND sum <= rangeC.hi *>
                pgmC.put(ROUND(scaleC * (sum - rangeC.lo)))
              END
            END;
          END
        END;
        pgmC.finish()
      END;
    END
  END DoConvol;
  
PROCEDURE ComputeMaxC(
    maxVal: NAT; (* Requested max output pixel *)
    MaxA: NAT;   (* Max mask pixel *)
    MaxB: NAT;   (* Max input pixel *)
    NP: NAT;     (* MAX number of overlapping pixels *)
  ): NAT (* Chosen max output pixel *) =
  CONST MaxDefaultMaxVal = MaxMaxVal - (MaxMaxVal MOD 2);
  BEGIN
    IF maxVal > 0 THEN
      RETURN maxVal
    ELSE
      WITH
        np = FLOAT(NP, LONG),
        fa = 1.0d0/FLOAT(MaxA, LONG),
        fb = 1.0d0/FLOAT(MaxB, LONG),
        maxC = Math.sqrt(np / (fa*fa + fb*fb))
      DO
        <* ASSERT MaxDefaultMaxVal MOD 2 = 0 *>
        IF maxC >= FLOAT(MaxDefaultMaxVal, LONG) THEN
          RETURN MaxDefaultMaxVal
        ELSE
          WITH tmp = CEILING(maxC), M = tmp + (tmp MOD 2) DO
            <* ASSERT M <= MaxDefaultMaxVal *>
            RETURN M
          END
        END
      END
    END
  END ComputeMaxC;

PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO

      TRY
        IF pp.keywordPresent("-maxVal") THEN
          o.maxVal := pp.getNextInt(1, MaxMaxVal)
        ELSE
          o.maxVal := 0 (* Compute it later from input "maxVal"s *)
        END;
        
        o.quiet := pp.keywordPresent("-quiet");

        pp.skipParsed();

        (* process mask file name: *)
        IF pp.next > LAST(pp.arg^) THEN 
          pp.error("missing kernel file name")
        ELSE
          WITH name = pp.getNext() DO
            IF Text.Equal(name, "-") THEN
              o.rdA := stdin
            ELSE
              TRY
                o.rdA := FileRd.Open(name)
              EXCEPT
              | OSError.E => pp.error("can't open file \"" & name & "\"")
              END
            END
          END
        END;

        (* Process input file name: *)
        IF pp.next > LAST(pp.arg^) THEN 
          o.rdB := stdin
        ELSE
          WITH name = pp.getNext() DO
            IF Text.Equal(name, "-") THEN
              o.rdB := stdin
            ELSE
              TRY
                o.rdB := FileRd.Open(name)
              EXCEPT
              | OSError.E => pp.error("can't open file \"" & name & "\"")
              END
            END
          END
        END;
        IF o.rdA = stdin AND o.rdB = stdin THEN 
          pp.error("cannot use stdinn for both kernel and input")
        END;

        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: PGMConvol [ -quiet ] \\\n");
          Wr.PutText(stderr, "  [ -maxVal <maxoutpixel> ] \\\n");
          Wr.PutText(stderr, "  <kernelfile> [ <infile> ]\n");
          Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;
  
BEGIN 
  Main();
END PGMConvol.