MODULE PZFractalDim EXPORTS Main;

(* Generates a plot to estimate the fractal dimension of a curve. *)
(* Last edited on 1999-08-07 06:34:05 by hcgl *)

IMPORT Wr, FileRd, FileWr, OSError, Thread, Stdio, Fmt, Math, ParseParams, Process;
IMPORT LR3;
IMPORT PZLR3Chain;

FROM Math IMPORT sin, log;
FROM LR3 IMPORT Dist;
FROM Stdio IMPORT stderr;   
FROM PZTypes IMPORT LONG, LONGS, INT, NAT;

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

TYPE
  MismatchOptions = PZMismatch.Options;
  Options = RECORD
      inFile: TEXT;      (* Input float chain file (without ".flc") *)
      outFile: TEXT;     (* Output plot data file (without ".plt") *)
      minDist: LONG;     (* Minimum beeline distance to consider *)
    END;
  
CONST 
  Pi = 3.14159265358979323d0;
  Sqrt2 = 1.4142135623730950488d0;
  Log2 = 0.69314718055994530941d0;

PROCEDURE Main() =      
  BEGIN 
    WITH 
      o = GetOptions(),
      fc = PZLR3Chain.Read(FileRd.Open(o.inFile & ".flc")),
      chain = fc.c^
    DO 
      WriteLengthScalingPlot(FileWr.Open(o.outFile & ".plt"), chain, o.minDist);
    END;
  END Main;


PROCEDURE GetOptions(): Options =
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO
      TRY
        pp.getKeyword("-inFile");
        o.inFile := pp.getNext();
        
        pp.getKeyword("-outFile");
        o.outFile := pp.getNext();
        
        pp.getKeyword("-minDist");
        o.minDist := pp.getNextLongReal(0.1d0, 100000.0d0);
        
        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: PZFilter \\\n");
          Wr.PutText(stderr, "  -inFile NAME \\\n");
          Wr.PutText(stderr, "  -outFile NAME \\\n");
          Wr.PutText(stderr, "  -minDist: NUM\n");
        Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;

  
PROCEDURE WriteLengthScalingPlot(wr: Wr.T; READONLY p: PZLR3Chain.T; minDist: LONG)=
  VAR j, k: INT;
      (* These arrays are indexed by "ROUND(log_2(dij/minDist)/2)": *)
      tSum: ARRAY [0..60] OF LONG;
      dSum: ARRAY [0..60] OF LONG;
      sSum: ARRAY [0..60] OF LONG;
      nSteps: ARRAY [0..60] OF NAT;
  BEGIN
    WITH 
      n = NUMBER(p),
      t = NEW(REF LONGS, n+1) 
    DO 
      Wr.PutText(stderr, "number of samples = " & Fmt.Int(n) & "\n");
      
      (* Store in "t[i]" the length of the curve from "p[0]" to "p[i]". *)
      t[0] := 0.00d0; 
      FOR i:= 1 TO n-1 DO 
        t[i] := t[i-1] + Dist(p[i-1], p[i])
      END;
      t[n] := t[n-1] + Dist(p[n-1], p[0]);
      Wr.PutText(stderr, "Total length = " & Fmt.LongReal(t[n]) & " pixels\n");
      
      (* Clear accumulators: *)
      FOR r := 0 TO LAST(tSum) DO 
        tSum[r] := 0.0d0; dSum[r] := 0.0d0; sSum[r] := 0.0d0;
        nSteps[r] := 0
      END;
      
      (* Now generate a plot of "log(L)/log(D)": *)
      k := 2;
      WHILE k < n DIV 2 DO 
        FOR i := 0 TO n - 1 BY k DIV 2 DO 
          j:= i + k; 
          WITH 
            (* Beeline distance (chord): *)
            dij = Dist(p[i], p[j MOD n]),
            
            (* Length along the curve: *)
            tij = t[j MOD n] + t[n]*FLOAT(j DIV n, LONG) - t[i],
            
            (* Half of angle subtended by "tij" on a circle of length "t[n]": *)
            theta2 = Pi * tij / t[n],
            
            (* Length of arc of circle with same angle and chord "dij": *)
            sij = dij * theta2 / ABS(sin(theta2))
          DO
            IF tij > minDist/Sqrt2 AND dij > minDist/Sqrt2 THEN 
              WITH
                r = MAX(0, ROUND(2.0d0*log(dij/minDist)/Log2))
              DO
                tSum[r] := tSum[r] + tij;
                dSum[r] := dSum[r] + dij;
                sSum[r] := sSum[r] + sij;
                INC(nSteps[r])
              END
            END
          END
        END;
        k := MAX(k+1, ROUND(FLOAT(k, LONG)*Sqrt2))
      END;
      
      (* Now output the plot: *)
      FOR r := 0 TO LAST(tSum) DO 
        IF nSteps[r] >= 2 THEN
          WITH
            ns = FLOAT(nSteps[r], LONG),
            tAvg = tSum[r]/ns,
            dAvg = dSum[r]/ns,
            sAvg = sSum[r]/ns,
            dDim = log(tSum[r])/log(dSum[r]),
            sDim = log(tSum[r])/log(sSum[r])
          DO
            Wr.PutText(wr, FLR(tAvg, 7, 2));
            Wr.PutText(wr, " ");
            Wr.PutText(wr, FLR(dAvg, 7, 2)); 
            Wr.PutText(wr, " ");
            Wr.PutText(wr, FLR(sAvg, 7, 2)); 
            Wr.PutText(wr, " ");
            Wr.PutText(wr, FLR(dDim, 7, 4));
            Wr.PutText(wr, " ");
            Wr.PutText(wr, FLR(sDim, 7, 4));
            Wr.PutText(wr, "\n");
            Wr.Flush(wr);
          END
        END
      END
    END; 
  END WriteLengthScalingPlot;
  
PROCEDURE FLR(x: LONG; w, p: NAT): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, style := Fmt.Style.Fix, prec := p), w)
  END FLR;

BEGIN
  Main()
END PZFractalDim.



(*
  Copyright © 2001 Universidade Estadual de Campinas (UNICAMP).
  Authors: Helena C. G. Leitão and Jorge Stolfi.
  
  This file can be freely distributed, used, and modified, provided
  that this copyright and authorship notice is preserved, and that any
  modified versions are clearly marked as such.
  
  This software has NO WARRANTY of correctness or applicability for
  any purpose. Neither the authors nor their employers chall be held
  responsible for any losses or damages that may result from its use.
*)