MODULE SmoothShape EXPORTS Main;

(*
  Applies smoothing heuristics to a given ".top" file.
  
  Also writes an ".st" file with a movie of the smoothing;
  a ".plot" file with the evolution of various energies during 
  the smoothing; and a ".gnuplot" file with the "gnuplot" commands
  to show the latter.
  
  Created 1994 by Rober M. Rosi and J. Stolfi.
*)

IMPORT Thread, Text, Wr, Fmt, FileWr, OSError, Process, ParseParams;
IMPORT LR3, Triang, Heuristics, Energy, MixedEnergy, ParseEnergyParams;
IMPORT RandomPerm, Random;
FROM Energy IMPORT Gradient;
FROM Triang IMPORT Topology, Coords;
FROM Stdio IMPORT stderr;

TYPE
  Options = RECORD              (* Input file name (minus ".top" extension) *)        
      inFile: TEXT;             (* Output file name prefix *)                         
      outFile: TEXT;            (* Number of configurations to generate *)            
      jitter: LONGREAL;         (* How much to randomize "inFile" *)
      seed: CARDINAL;           (* Random number seed *)
      flatten: CARDINAL;        (* Passes of "flatten" heuristic *)     
      spread: CARDINAL;         (* Passes of "spread" heuristic *)  
      normalize: BOOLEAN;       (* Keep edge lengths normalized *)  
      eFunction: MixedEnergy.T; (* Optional: energy to watch *)
      showEvery: CARDINAL;      (* When to write ".top" file. *)
    END;
    
TYPE
  SmoothProc = PROCEDURE(
      READONLY top: Topology; 
      VAR c: Coords; 
      READONLY variable: ARRAY OF BOOLEAN;
      coins: Random.T;
      jitter: LONGREAL;
    );

PROCEDURE DoIt() =
  VAR flattenC: CARDINAL := 0;
      spreadC: CARDINAL := 0;
      h: SmoothProc;
      e: LONGREAL;
  <* FATAL Thread.Alerted, Wr.Failure, OSError.E *>
  BEGIN
    WITH 
      o = GetOptions(), 
      tc = Triang.Read(o.inFile),
      top = tc.top,
      c = tc.c^,
      NV = top.NV,
      coins = NEW(Random.Default).init(TRUE),
      g = NEW(REF Gradient, NV)^,
      variable = NEW(REF ARRAY OF BOOLEAN, NV)^,
      plotWr = FileWr.Open(o.outFile & ".plot"),
      stWr = FileWr.Open(o.outFile & ".st"),
      NPasses = o.flatten + o.spread
    DO
      (* Initialize the generator: *)
      FOR i := 1 TO o.seed DO EVAL coins.integer() END;
      
      (* Write triangles for animation: *)
      Triang.WriteRLuisSuMa(o.outFile, top, all := FALSE, comments := tc.comments);

      (* Get non-fixed vertices: *)
      Triang.GetVariableVertices(top, variable);
      
      (* Write energy plot header and "gnuplot" commands: *)
      IF o.eFunction # NIL THEN
        o.eFunction.defTop(top);
        o.eFunction.defVar(variable);
        WritePlotComments(plotWr, o.eFunction);
        WritePlotHeader(plotWr, o.eFunction);
        WITH gnuWr = FileWr.Open(o.outFile & ".gnuplot") DO
          WriteGNUPlotCommands(gnuWr, o.outFile, o.eFunction.term^);
          Wr.Close(gnuWr)
        END;
      END;
      
      FOR i := 0 TO NPasses DO 
        (* Choose the heuristic to use next: *)
        WITH
          flattenR = (FLOAT(flattenC)+0.25)/FLOAT(o.flatten),
          spreadR = (FLOAT(spreadC)+0.75)/FLOAT(o.spread),
          minR = MIN(flattenR, spreadR)
        DO
          IF i = 0 THEN
            h := JitterVertices;
          ELSIF minR = flattenR THEN
            h := FlattenVertices; INC(flattenC)
          ELSIF minR = spreadR THEN
            h := SpreadVertices;  INC(spreadC)
          ELSE
            <* ASSERT FALSE *>
          END;
        END;
        (* Apply the heuristic: *)
        h(top, c, variable, coins, o.jitter);
        IF o.normalize THEN 
          Triang.NormalizeVertexNorms(top, c);
        END;
        (* Evaluate the energy: *)
        IF o.eFunction # NIL THEN
          o.eFunction.eval(c, e, FALSE, g);
          PlotEnergy(
            plotWr, i, e, 
            o.eFunction.termValue^, o.eFunction.weight^
          );
        END;
        IF i = 0 THEN
          WITH
            name = o.outFile & "-initial"
          DO
            WriteConfiguration(name, top, c, flattenC, spreadC, tc.comments)
          END
        END;
        IF NiceNumberCoords(i, NPasses, o.showEvery) THEN
          WriteCoords(stWr, c, g, flattenC, spreadC, tc.comments)
        END;
      END;
      WITH
        name = o.outFile & "-final"
      DO
        WriteConfiguration(name, top, c, flattenC, spreadC, tc.comments)
      END;
      Wr.Close(plotWr);
      Wr.Close(stWr)
    END
  END DoIt;
  
PROCEDURE NiceNumberCoords(
    pass: CARDINAL; 
    NPasses: CARDINAL;
    every: CARDINAL;
  ): BOOLEAN =
  BEGIN
    RETURN (pass MOD every = 0)
      OR pass <= 10
      OR pass = 10 OR pass = 20 OR pass = 50
      OR (pass MOD 100 = 0 AND pass <= 5000) 
      OR (pass MOD 500 = 0)
      OR (pass = NPasses)
  END NiceNumberCoords;
  
PROCEDURE WriteConfiguration(
    name: TEXT;
    READONLY top: Topology;
    READONLY c: Coords;
    flattenC, spreadC: CARDINAL;
    comments: TEXT;
  ) =
  BEGIN
    WITH
      cmts = 
        "Output of SmoothShape after " &
        Fmt.Int(flattenC) & " \"flatten\" passes" & " and " &
        Fmt.Int(spreadC)  & " \"spread\" passes" & "\n" &
        comments
    DO
      Triang.Write(name, top, c, cmts)
    END
  END WriteConfiguration;
  
PROCEDURE WriteCoords(
    wr: Wr.T;
    READONLY c: Coords;
    READONLY g: Gradient;
    flattenC, spreadC: CARDINAL;
    comments: TEXT;
  ) =
  BEGIN
    WITH
      cmts = 
        "Output of SmoothShape after " &
        Fmt.Int(flattenC) & " \"flatten\" passes" & " and " &
        Fmt.Int(spreadC)  & " \"spread\" passes" & "\n" &
        comments
    DO
      Triang.WriteRLuisSt(wr, c, g, comments := cmts)
    END
  END WriteCoords;
  
PROCEDURE PlotEnergy(
    wr: Wr.T;
    iteration: CARDINAL;
    e: LONGREAL;
    READONLY termValue: ARRAY OF LONGREAL;
    READONLY weight: ARRAY OF REAL;
  ) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad(Fmt.Int(iteration), 8));
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad(Fmt.LongReal(e, Fmt.Style.Fix, 5), 10));
    FOR i := 0 TO LAST(termValue) DO 
      Wr.PutText(wr, " ");
      WITH
        t = termValue[i] * FLOAT(weight[i], LONGREAL)
      DO
        Wr.PutText(wr, Fmt.Pad(Fmt.LongReal(t, Fmt.Style.Fix, 5), 10))
      END;
    END;
    Wr.PutText(wr, "\n");
    Wr.Flush(wr);
  END PlotEnergy;
  
PROCEDURE WritePlotComments(wr: Wr.T; e: MixedEnergy.T) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    Wr.PutText(wr, "#");
    Wr.PutText(wr, "energy function: " & e.name());
    Wr.PutText(wr, "\n");

    Wr.Flush(wr);
  END WritePlotComments;

PROCEDURE WritePlotHeader(wr: Wr.T; eFunction: MixedEnergy.T) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    Wr.PutText(wr, "#");
    Wr.PutText(wr, Fmt.Pad("iter", 8));
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad("energy", 10));
    FOR i := 0 TO LAST(eFunction.term^) DO 
      Wr.PutText(wr, " ");
      Wr.PutText(wr, Fmt.Pad(EnergyTag(eFunction.term[i].name()), 10));
    END;
    Wr.PutText(wr, "\n");
    Wr.Flush(wr);
  END WritePlotHeader;

PROCEDURE WriteGNUPlotCommands(
    wr: Wr.T; 
    outFile: TEXT;
    READONLY term: ARRAY OF Energy.T;
  ) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    Wr.PutText(wr, "set terminal X11\n");
    Wr.PutText(wr, "set title \"" & outFile & "\"\n");
    Wr.PutText(wr, 
      "plot \"" & outFile & ".plot\" using 1:2 title \"etot\" with lines, \\\n"
    );
    FOR i := 0 TO LAST(term) DO 
      WITH col = i + 3 DO
        Wr.PutText(wr, 
          "     \"" & outFile & ".plot\" using 1:" & 
          Fmt.Int(col) & " title \"" &
          EnergyTag(term[i].name()) & 
          "\" with linespoints"
        );
        IF i < LAST(term) THEN Wr.PutText(wr, ", \\") END;
        Wr.PutText(wr, "\n");
      END
    END;
    Wr.PutText(wr, "pause 300\n");
    Wr.PutText(wr, "quit\n");
    Wr.Flush(wr);
  END WriteGNUPlotCommands;
  
PROCEDURE EnergyTag(name: TEXT): TEXT =
  BEGIN
    WITH
      n = Text.FindChar(name, '(')
    DO
      IF n = -1 THEN 
        RETURN Text.Sub(name, 0, 8)
      ELSE
        RETURN Text.Sub(name, 0, MIN(n, 8))
      END
    END
  END EnergyTag;

PROCEDURE JitterVertices(
    READONLY top: Topology;
    VAR c: Coords; 
    READONLY variable: ARRAY OF BOOLEAN;
    coins: Random.T; 
    jitter: LONGREAL;
  ) =
  VAR R: LONGREAL := Triang.MeanVertexNorm(top, c);
      p: LR3.T;
  BEGIN
    IF jitter = 0.0d0 THEN RETURN END;
    FOR i := 0 TO LAST(c) DO
      IF variable[i] THEN
        REPEAT
          p := LR3.T{
            coins.longreal(-1.0d0, +1.0d0),
            coins.longreal(-1.0d0, +1.0d0),
            coins.longreal(-1.0d0, +1.0d0)
          }
        UNTIL LR3.Norm(p) <= 1.0d0;
        c[i] := LR3.Mix(1.0d0 - jitter, c[i], jitter * R, p)
      END
    END
  END JitterVertices;

PROCEDURE FlattenVertices(
    READONLY top: Topology; 
    VAR c: Coords; 
    READONLY variable: ARRAY OF BOOLEAN;
    coins: Random.T;
    <*UNUSED*> jitter: LONGREAL;
  ) =
  <* FATAL RandomPerm.Exhausted *>
  BEGIN
    WITH perm = NEW(RandomPerm.LowQuality).init(top.NV, coins) DO
      FOR i := 0 TO top.NV-1 DO 
        WITH v = perm.next() DO
          IF variable[v] THEN
            Heuristics.FlattenVertex(top.out[v], c, variable, coins)
          END
        END
      END;
    END;
  END FlattenVertices;
  
PROCEDURE SpreadVertices(
    READONLY top: Topology; 
    VAR c: Coords; 
    READONLY variable: ARRAY OF BOOLEAN;
    coins: Random.T;
    <*UNUSED*> jitter: LONGREAL;
  ) =
  <* FATAL RandomPerm.Exhausted *>
  BEGIN
    WITH perm = NEW(RandomPerm.LowQuality).init(top.NV, coins) DO
      FOR i := 0 TO top.NV-1 DO 
        WITH v = perm.next() DO
          IF variable[v] THEN
            Heuristics.SpreadVertex(top.out[v], c, variable, coins)
          END;
        END
      END;
    END;
  END SpreadVertices;

PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  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();

        IF pp.keywordPresent("-jitter") THEN
          o.jitter := pp.getNextLongReal(0.0d0, 1.0d0);
          IF pp.keywordPresent("-seed") THEN
            o.seed := pp.getNextInt(0, 999)
          ELSE
            o.seed := 0
          END
        ELSE
          o.jitter := 0.0d0
        END;
        
        IF pp.keywordPresent("-flatten") THEN
          o.flatten := pp.getNextInt(0, 1000000)
        ELSE
          o.flatten := 0
        END;
        
        IF pp.keywordPresent("-spread") THEN
          o.spread := pp.getNextInt(0, 1000000)
        ELSE
          o.spread := 0
        END;
        
        IF o.flatten + o.spread = 0 THEN
          pp.error("no smoothing heuristic specified");   
        END;
        
        o.normalize := pp.keywordPresent("-normalize");
        
        o.eFunction := ParseEnergyParams.Parse(pp);
        
        IF pp.keywordPresent("-showEvery") THEN
          o.showEvery := pp.getNextInt(0, 1000000)
        ELSIF pp.keywordPresent("-showAll") THEN
          o.showEvery := 1
        ELSE
          o.showEvery := LAST(CARDINAL)
        END;

        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: SmoothShape \\\n");
          Wr.PutText(stderr, "  -inFile <name>   -outFile <name> \\\n");
          Wr.PutText(stderr, "  [ -jitter <num> [ -seed <num> ] ] \\\n");
          Wr.PutText(stderr, "  [ -flatten <num> ] \\\n");
          Wr.PutText(stderr, "  [ -spread <num> ] \\\n");
          Wr.PutText(stderr, "  [ -normalize ] \\\n");
          Wr.PutText(stderr, "  [ -showEvery <num> | -showAll ] \\\n");
          Wr.PutText(stderr, ParseEnergyParams.Help);
          Wr.PutText(stderr, "\n");
          Process.Exit (1);                                              
      END
    END;
    RETURN o
  END GetOptions; 
  
BEGIN
  DoIt()
END SmoothShape.