#define PROG_NAME "Visibility"
#define PROG_DESC "???"
#define PROG_VERS "1.0"
/* Last edited on 2024-12-21 11:36:43 by stolfi */

#define Visibility_C_COPYRIGHT \
  ""
  
#define PROG_INFO \
  "" \
  " "


  /*
   This program implements the "visibility" of a 3D projection.

   Given a 4D (optimized) configuration, and the viewing parameters of
   the 4D observer (From4, Up4, etc), writes a POVRAY ".inc" file con-
   taining the 3D projection of the 4D configuration as "seen" from that
   viewpoint.

   The projection will use distinct textures (and optionally omit) the
   "front" and "back" elements, as well as, the silhouette walls.

   The back/front test assumes that the tetrahedra constitute the shell
   of a convex polytope and that the observer is outside that shell.
   The test consists in computing the orientation of each tetrahedron of
   the 3D-projected triangulation with reference to the observer and with
   reference to the "barycenter" of the configuration.

   Depending the position of the 4D observer we obtain different orientations.

   Use the option "both" to show both front and back elements. Use the
   option "-fade" to paint the back elements with faded colors. Use "back"
   or "front" to show only those elements.

   If "From4" is not given as an option, the program reads zero or
   more lines from stdin, each containing an index "i", a value of "From4",
   one output option ("front"/"back"/"both") and writes one ".inc" file for
   each such line.
   
*/

#include <stdio.h>
#include <string.h>

#include <argparser.h>
#include <r3.h>
#include <fget.h>
#include <r4.h>

#include <Pov.h>
#include <Octf.h>
#include <Triangulation.h>
#include <Tridimensional.h>

CONST
double Epsilon = 0.0000000001;

TYPE
 Options_t == RECORD
  char *inFileTp;        /* Input file prefix (topology, tabel, materials).*/
  char *inFileSt;        /* Input file prefix (4D geometry). */
  char *outFile;         /* Output file prefix (POVray include files). */
  bool_t wire;  /* Don't draw walls. */
  bool_t detail;  /* Debugging. */
  bool_t silhouette;  /* Draw the  silhouette walls. */
  bool_t perspective;  /* TRUE means perspective projection. */
  bool_t autoProject;  /* The program choses the viewpoint. */
  bool_t filter;  /* TRUE uses filtered colors, FALSE transmit colors.*/
  bool_t fade;  /* Make the back elements fainter. */
  bool_t normalize;  /* Norm. al nodes onto the S^3 with that radius. */
  bool_t multiple;     /* T reads multiple "From4"s and "side"s from stdin.*/
  From4:   r4_t;        /* Viewpoint in 4D. */
  To4:     r4_t;        /* 4D viewing parameters as expected by the  */
  Up4:     r4_t;        /* "Wire4"- Interactive 4D Wireframe Display */
  Over4:   r4_t;        /* Program.                                  */
  Vangle4: double;     /*                                           */
  Side select;       /* Which elements to show in the output.     */
  char *selectTxt;       /* Ditto, in text format.                    */
  bool_t colored;    /* coloring the original walls with a palette*/
 }

    /* [[!!! DOES THIS COMMENT FROM Octf.h APPLY HERE ??? */
   /* A '@{chip->?}' is a group of four @places, consisting of two mutually dual
   wedges. */

FourNodes_t  == RECORD u, v, w, x: uint; }
 Side  == {Back, Front, Both};
 Color == frgb_t;
 PACK  == RECORD bool_t belong; num: INTEGER; }

 <* FATAL Rd.Failure, Wr.Failure, Thread.Alerted, FloatMode.Trap, Lex.Error);

PROCEDURE @{Edge->?}BelongsToOriginalWall(Place_t @p, ElemTableRec_t *top): PACK ==
  /* 
    Returns TRUE if the edge component of "a" belongs to a wall of the 
    original map also accomplish the number of the original wall. */
  VAR
     Place_t b = a;
    PACK pack;
  { if ((top->edge[PEdge(a)->num]->root!=-1)){
      pack.belong = FALSE;
      pack->num = -1;
      return pack;
    }
    /* Check if any wall incident to "a" belongs to a wall of the map. */
    do 
    {
      if ((top->wall[PWall(b)->num]->root!=-1))
      {
        pack.belong = TRUE;
        pack->num = PWall(b)->num;
        return pack;
      }
      b = NextF(b)
    } while (b != a);
    pack.belong = FALSE;
    pack->num = -1;
    return pack;
  } /* END @{Edge->?}BelongsToOriginalWall */

bool_t Sign(d: double)
  /* Return TRUE iff the longreal value is positive, FALSE c.c. */
  { /*assert(d!=0.0);*/
    if (d < 0.0){ return FALSE }else{ return TRUE; };
  } /* END Sign */

Options_t *GetOptions(int argc, char **argv);

int main(int argc, char **argv)
  Options_t *o = GetOptions(argc, argv);
    char *topo_cmt = jsprintf("Created by %s on %s", PROG_NAME, Today());
    /* Random_t coins = MakeRandomSource(4615); */
    uint i;
  { ??? tc = Triangulation.ReadToMa(o->inFileTp);
    ??? top = tc.top;
    ??? rc4 = Triangulation.ReadState(o->inFileSt);
    ??? rc3 = NEW(REF Tridimensional.Coords3D_t;
    with (top->node.nel),
      c4  == rc4^,
      c3  == rc3^
  ){
      if (o->normalize)
      {
        fprintf(stderr, "projecting nodes onto the unit S^3\n");
        Triangulation.NormalizeNodeDistances(top, c4, o->normalize);
      }
      if (o->autoProject)
      {
        fprintf(stderr, "selecting camera parameters\n");
        SelectProjection(o, c4, top);
      }
      fprintf(stderr, "projecting and converting to POV-Ray\n");
      if (o->multiple)
      {
        while (1)
        {
          Lex.Skip(stdin, Lex.Blanks);
          if ((Rd.EOF(stdin))){ EXIT; }
          i = fget_int32(stdin);
          o->From4[0] = Lex.LongReal(stdin);
          o->From4[1] = Lex.LongReal(stdin);
          o->From4[2] = Lex.LongReal(stdin);
          o->From4[3] = Lex.LongReal(stdin);
          Lex.Skip(stdin, Lex.Blanks);
          o->selectTxt = Lex.Scan(stdin);
          if (0 == strcmp(o->selectTxt, "back")))
          {
            o->select = Side.Back;
          }
          else
           if (0 == strcmp(o->selectTxt, "front"))){
            o->select = Side.Front;
          }
          else
           if (0 == strcmp(o->selectTxt, "both"))){
            o->select = Side.Both;
          }
          else
          {
            fprintf(stderr, "Bad shape \"" &  o->selectTxt & "\"\n");
            Process.Exit¦(1)
          }
          ProjectTo3D(o, c4, c3, top);
          ??? itx = Fmt.Pad(Fmt.Int(i), 4, '0');
    {
            WritePOV(o->outFile & "-" & o->selectTxt & "-" & itx, top, c3, c4, o)
          }
        }
      }
      else
      {
        ProjectTo3D(o, c4, c3, top);
        WritePOV(o->outFile & "-" & o->selectTxt, top, c3, c4, o)
      }
    
    return 0;
  }

void WritePOV(char *name, ElemTableRec_t *top, 
    * c3: Tridimensional.Coords3D_t;
    * c4: Coords_t;
    Options_t *  o;
 )
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E);
  { ??? wr = FileWr.Open(name & ".inc");
    {
      fprintf(wr, "// Include File: <" & name & ".inc>\n");   
      DoWritePOV(wr, top, c3, c4, o);
      fclose(wr)
    }
  } /* END WritePOV */

PROCEDURE DoWritePOV
  ( FILE *wr,
   *ElemTableRec_t *top;
   * c3: Tridimensional.Coords3D_t;
   Coords_t *  c;
   Options_t *  o;
 ) ==
    VAR
      vseen = NEW(REF ARRAY OF bool_t, top->node.nel); /* stack for nodes */
      eseen = NEW(REF ARRAY OF bool_t, top->NE); /* stack for edges    */
      fseen = NEW(REF ARRAY OF bool_t, top->wall.nel); /* stack for walls    */
      quad  : REF ARRAY OF FourNodes_t;

    bool_t IsFront(q: FourNodes_t)
      {
        ??? u = r4_Sub(c[q.u];
        with (cb),
          v  == r4_Sub(c[q.v],cb),
          w  == r4_Sub(c[q.w],cb),
          x  == r4_Sub(c[q.x],cb), 
        double d1 = r4_det(u,v,w,x);
        double u_ = r4_Sub(c[q.u],o->From4);
        double v_ = r4_Sub(c[q.v],o->From4);
        double w_ = r4_Sub(c[q.w],o->From4);
        double x_ = r4_Sub(c[q.x],o->From4);
        double d2 = r4_det(u_,v_,w_,x_)
        DO
          return (Sign(d1)) && (NOT Sign(d2))) || (( NOT Sign(d1)) && (Sign(d2));
        }
      } IsFront;

    double FindOrientR3(q: FourNodes_t)
      {
        with (
        double a = (r4_t){c3[q.u][0], c3[q.u][1], c3[q.u][2], 1.0};
        double b = (r4_t){c3[q.v][0], c3[q.v][1], c3[q.v][2], 1.0};
        double c = (r4_t){c3[q.w][0], c3[q.w][1], c3[q.w][2], 1.0};
        double d = (r4_t){c3[q.x][0], c3[q.x][1], c3[q.x][2], 1.0}
      ){
          return r4_det(a,b,c,d);
        }   
      } FindOrientR3;

    bool_t IsSilhouette(Place_t @p)
      /* Return TRUE iff the wall associated to the @place "a" is a silhouette
         wall, FALSE c.c.
      */
      {
        ??? t = TetraNegPosNodes(a);
        ??? un = t[0]->num;
        ??? vn = t[1]->num;
        ??? wn = t[2]->num;
        ??? xn = t[3]->num;
        ??? yn = t[4]->num;
        with (
        double t1 = FourNodes_t{un,vn,wn,xn};
        double t2 = FourNodes_t{un,vn,wn,yn};
        double d1 = FindOrientR3(t1);
        double d2 = FindOrientR3(t2)
        DO  
          return (Sign(d1)) && (Sign(d2))) || ((NOT Sign(d1)) && (NOT Sign(d2));
        }
      } IsSilhouette;

  void WritePOVNode(s: Triangulation.Node; fade: REAL; fbTag: char *)
    char **dimTag;
    {
      ??? ff = ((double)fade);
        ??? ts = s.transp;
      ??? ts = (ts[0] + ts[1] + ts[2]) / 3.0 * fade;
      ??? tc = R3.Mix(ff;
      with (s.color, 1.0 - ff, (frgb_t){0.5,0.5,0.5})
    ){
        if ((s->exists) && (NOT vseen[s->num]))
        {
          if (0 == strcmp(s.label, "VV"))){ dimTag = "node"
            }
            else
            { dimTag = "@{edge->?}"
          }
          WritePOVSphere(wr, c3[s->num], s.radius, tc, ts, TRUE, dimTag, fbTag);
          vseen[s->num] = TRUE
        }
      }
    } WritePOVNode;

  void WritePOVWall(t: Triangulation.Wall; fade: REAL; fbTag: char *)
    {
      ??? ff = ((double)fade);
        ??? ts = t.transp;
      ??? ts = (ts[0] + ts[1] + ts[2]) / 3.0 * fade;
      ??? tc = R3.Mix(ff;
      with (t.color, 1.0 - ff, (frgb_t){0.5,0.5,0.5}),
       double a = t.pa;
       double b = NextE(a);
       double c = NextE(b);
      double v0 = OrgV(a)->num;
      double v1 = OrgV(b)->num;
      double v2 = OrgV(c)->num
    ){
        if ((t->exists) && (NOT fseen[t->num]))
        {
          WritePOVTriangle(wr, c3[v0], c3[v1], c3[v2], tc, ts,
            TRUE, "wall", fbTag);
          fseen[t->num] = TRUE
        }
      }
    } WritePOVWall;

  void WritePOVSilhouetteWall(t: Triangulation.Wall)
    {
      ??? ts = 0.99;
        with (
      double tc = (frgb_t){1.00, 0.99, 0.99};
       double a = t.pa;
       double b = NextE(a);
       double c = NextE(b);
      double v0 = OrgV(a)->num;
      double v1 = OrgV(b)->num;
      double v2 = OrgV(c)->num
    ){
        WritePOVTriangle(wr, c3[v0], c3[v1], c3[v2], tc, ts, o->filter,
          "wall", "silhouette");
      }
    } WritePOVSilhouetteWall;
    
  PROCEDURE WritePOV@{Edge->?}(a: Triangulation.edge; fade: REAL;
    Color color; fbTag: char *) ==
    {
      ??? ff = ((double)fade);
        ??? ts = a.transp;
      ??? ts = (ts[0] + ts[1] + ts[2]) / 3.0 * fade;
      ??? tc = R3.Mix(ff;
      with (color, 1.0 - ff, (frgb_t){0.5,0.5,0.5}),
      double v0 = OrgV(a.pa)->num;
      double v1 = OrgV(Clock(a.pa))->num
    ){
        if ((a->exists) && (NOT eseen[a->num]))
        {
          WritePOVCylinder(wr, c3[v0], c3[v1], a.radius, tc, ts, o->filter,
            "@{edge->?}", fbTag);
          eseen[a->num] = TRUE
        };
      }
    } WritePOV@{Edge->?};

  PROCEDURE Write@{Edge->?}Colored(e: Triangulation.edge) ==
    /* draws one colored @{edge->?}. */
    {  
      ??? a = e.pa;
      {
        if ((e->exists) && (@{Edge->?}BelongsToOriginalWall(a,top).belong))
        {
          /* an @{edge->?} on the net belonging to some original wall */
          ??? facc = @{Edge->?}BelongsToOriginalWall(a;
          with (top)->num,
          double froot = top->wall[facc]->root
        ){                
            WritePOV@{Edge->?}(e, fadeFactor, fac[froot], "Net");
          }
        }
        else
         if ((e->exists) && (NOT @{Edge->?}BelongsToOriginalWall(a,top).belong)){ 
          /* an @{edge->?} of the map */
          WritePOV@{Edge->?}(e, fadeFactor, e.color, "Map");
        }
      }
    } Write@{Edge->?}Colored;
    
  VAR
    bool_t isFront;
    REAL fadeFactor;
    char *textureTag;
    cb: r4_t;
    Color palette[5+1];
    max : uint = 0;
    fac : REF ARRAY OF Color;  /* the wall colors to use */
    
  { cb    = Barycenter(top,c,TRUE);
    quad  = NEW(REF ARRAY OF FourNodes_t, top->cell.nel);
    for(i = 0; i < vseen.nel; i++){ vseen[i] = FALSE; }
    for(i = 0; i < fseen.nel; i++){ fseen[i] = FALSE; }
    for(i = 0; i < eseen.nel; i++){ eseen[i] = FALSE; }
    if (o->colored)
    {
      /* computing the original number walls */
      for (i = 0; i < top->wall.nel; i++)
      {
        ??? f = top->wall[i];
    {
          if (f->exists)
          {
            max = MAX(max,f->root);
          }
        }
      }
      fprintf(stderr, "The original number walls is " \
        Fmt.Int(max+1)&"\n");

      fac = NEW(REF ARRAY OF Color,max+1);

      /* defining the palette colors with brightness aprox. 0.7 */
      palette[0] = Color{1.00,0.50,1.00};
      palette[1] = Color{0.00,1.00,1.00};
      palette[2] = Color{0.80,0.80,0.00};
      palette[3] = Color{1.00,0.60,0.60};
      palette[4] = Color{0.25,1.00,0.25};
      palette[5] = Color{0.70,0.70,1.00};

      /* attribution of differents colors for each original wall */
      for (j = 0; j <= (max); j++)
      {
        fac[j] = palette[j % 6];
      }
      for (i = 0; i < top->wall.nel; i++)
      {
        ??? f = top->wall[i];
    {
          if (f->exists)
          {
            ??? froot = f->root;
    { f.color = fac[froot];}
          }
        }
      }
    }
    
    for (i = 0; i < top->cell.nel; i++)
    {
      ??? da = Srot(top.cell[i]);
        Place_t a0 = Tors(da);
      Place_t db = Clock(PrevE(da));
      Place_t b0 = Tors(db);
      Place_t a1 = NextE(a0);
      Place_t a2 = NextE(a1);
      Node_t un = OrgV(a0)->num;
      Node_t vn = OrgV(a1)->num;
      Node_t wn = OrgV(a2)->num;
      Node_t xn = OrgV(PrevE(b0))->num;
      {
        quad[i] = FourNodes_t{un,vn,wn,xn};
        isFront = IsFront(quad[i]);
        if (o->detail)
        {
          fprintf(stderr,"Tetrahedron " & Fmt.Int(i) & " on ");
        }
        if (((isFront) && (o->select!=Side.Back) OR
          (NOT isFront) && (o->select!=Side.Front))){
          if (o->detail)
          {
            fprintf(stderr, "isFront == " & Fmt.Int(ORD(isFront)) & "\n");
          }
          if ((o->fade) && (NOT isFront))
          {
            fadeFactor = 0.5; textureTag = "back" 
          }
          else
          {
            fadeFactor = 1.0; textureTag = "front"
          }
          ??? ri = Srot(top.cell[i]);
            Place_t a = Tors(ri);
          Place_t db = Clock(PrevE(ri));
          Place_t b0 = Tors(db);
          Place_t a1 = NextE(a);
          Place_t a2 = NextE(a1);
          Node_t v0 = OrgV(a)->num;
          ??? v1 = OrgV(a1)->num;
          ??? v2 = OrgV(a2)->num;
          ??? v3 = OrgV(PrevE(b0))->num;
          ??? s0 = top->node[v0];
          ??? s1 = top->node[v1];
          ??? s2 = top->node[v2];
          ??? s3 = top->node[v3];
          with (
            /* walls */
          double walls = Triangulation.TetraWalls(a);
          double f0 = walls[0]->num;
          double f1 = walls[1]->num;
          double f2 = walls[2]->num;
          double f3 = walls[3]->num;
          double t0 = top->wall[f0];
          double t1 = top->wall[f1];
          double t2 = top->wall[f2];
          double t3 = top->wall[f3];
            /* @{edge->?}s */
          double @{edge->?}s = Triangulation.TetraEdges(a);
          double e0 = @{edge->?}s[0]->num;
          double e1 = @{edge->?}s[1]->num;
          double e2 = @{edge->?}s[2]->num;
          double e3 = @{edge->?}s[3]->num;
          double e4 = @{edge->?}s[4]->num;
          double e5 = @{edge->?}s[5]->num;
          double a0 = top->edge[e0];
          double a1 = top->edge[e1];
          double a2 = top->edge[e2];
          double a3 = top->edge[e3];
          double a4 = top->edge[e4];
          double a5 = top->edge[e5];
          double dc = (frgb_t){1.00,1.000,0.500} /* DF color por @{edge->?} net */
        ){
            WritePOVNode(s0, fadeFactor, textureTag);
            WritePOVNode(s1, fadeFactor, textureTag);
            WritePOVNode(s2, fadeFactor, textureTag);
            WritePOVNode(s3, fadeFactor, textureTag);

            /* write walls */
            if (! o->wire)
            {
              WritePOVWall(t0, fadeFactor, textureTag);
              WritePOVWall(t1, fadeFactor, textureTag);
              WritePOVWall(t2, fadeFactor, textureTag);
              WritePOVWall(t3, fadeFactor, textureTag);
            }
            if (! o->colored)
            {
              /* write @{edge->?}s */
              if (a0->root!=-1)
              {
                WritePOV@{Edge->?}(a0, fadeFactor, a0.color, textureTag);
              }
              else
              {
                WritePOV@{Edge->?}(a0, fadeFactor, dc, "Net");
              }
              if (a1->root!=-1)
              {
                WritePOV@{Edge->?}(a1, fadeFactor, a1.color, textureTag);
              }
              else
              {
                WritePOV@{Edge->?}(a1, fadeFactor, dc, "Net");  
              }
              if (a2->root!=-1)
              {
                WritePOV@{Edge->?}(a2, fadeFactor, a2.color, textureTag);
              }
              else
              {
                WritePOV@{Edge->?}(a2, fadeFactor, dc, "Net");
              }
              if (a3->root!=-1)
              {
                WritePOV@{Edge->?}(a3, fadeFactor, a3.color, textureTag);
              }
              else
              {
                WritePOV@{Edge->?}(a3, fadeFactor, dc, "Net");
              }
              if (a4->root!=-1)
              {
                WritePOV@{Edge->?}(a4, fadeFactor, a4.color, textureTag);
              }
              else
              {
                WritePOV@{Edge->?}(a4, fadeFactor, dc, "Net");
              }
              if (a5->root!=-1)
              {
                WritePOV@{Edge->?}(a5, fadeFactor, a5.color, textureTag);
              }
              else
              {
                WritePOV@{Edge->?}(a5, fadeFactor, dc, "Net");
              }
            }
            else
            {
              Write@{Edge->?}Colored(a0);
              Write@{Edge->?}Colored(a1);
              Write@{Edge->?}Colored(a2);
              Write@{Edge->?}Colored(a3);
              Write@{Edge->?}Colored(a4);
              Write@{Edge->?}Colored(a5);
            }          
          }
        }
        else
        {
          if (o->detail)
          {
            fprintf(stderr, "front == " & Fmt.Int(ORD(isFront)) & "\n");
          }
        }        
      }
    }

    if (o->silhouette)
    {
      for (i = 0; i < top->wall.nel; i++)
      {
        ??? f = top->wall[i];
        {
          if (! f->exists)
          {
            with (a  == f.pa)
            {
              if ((IsSilhouette(a)))
              {
                if (o->detail)
                {
                  fprintf(stderr, "Silhouette wall: " & Fmt.Int(i) & "\n");
                }
                WritePOVSilhouetteWall(f)
              }
            }
          }
        }
      }
    }
  } /* END DoWritePOV */    

PROCEDURE ProjectTo3D(
    Options_t *o = (Options_t *)malloc(sizeof(Options_t));
    Coords_t *c;
    VAR c3: Tridimensional.Coords3D_t;
    ElemTableRec_t *top
   )
  

  void CalcV4Matrix()
    /* This procedure computes the four basis vectors for the 4D viewing
       matrix, Wa,Wb,Wc, and Wd. Note that the Up vector transforms to Wb,
       the Over vector transforms to Wc, and the line of sight transforms to
       Wd.  The Wa vector is then computed from Wb,Wc and Wd.
    */
    
    { 
      /* Calculate Wd, the 4th coordinate basis vector and line-of-sight. */
      Wd = r4_Sub(o->To4,o->From4);
      norm = r4_Norm(Wd);  
   
      if (norm < Epsilon)
      {
        fprintf(stderr,"4D To Point and From Point are the same\n");
        Process.Exit¦(1);
      }      
      Wd = r4_Scale(1.0/norm, Wd);

      /* Calculate Wa, the X-axis basis vector. */

      Wa = r4_cross(o->Up4,o->Over4,Wd);
      norm = r4_Norm(Wa);
      if (norm < Epsilon)
      {
        fprintf(stderr,
                   "4D up,over and view vectors are not perpendicular\n");
        Process.Exit¦(1);
      }      
      Wa = r4_Scale(1.0/norm, Wa);

      /* Calculate Wb, the perpendicularized Up vector. */

      Wb = r4_cross(o->Over4,Wd,Wa);
      norm = r4_Norm(Wb);
      if (norm < Epsilon)
      {
         fprintf(stderr,"Invalid 4D over vector\n");
         Process.Exit¦(1);
      }       
      Wb = r4_Scale(1.0/norm, Wb);

      /* Calculate Wc, the perpendicularized Over vector.  Note that the
         resulting vector is already normalized, since Wa, Wb and Wd are all
         unit vectors.
      */

      Wc = r4_cross(Wd,Wa,Wb);

    } CalcV4Matrix;  

  VAR
    Tan2Vangle4, Data4Radius, pconst, rtemp, depth: double;
    TempV : r4_t;
    Wa,Wb,Wc,Wd : r4_t;
    double norm;
  { ??? angle = o->Vangle4/2.0;
    ??? angler = (FLOAT(Math.Pi, double)*angle)/180.0);
  with(  )
  {
      Tan2Vangle4 = Math.tan(angler);
    }

    /* Find the radius of the 4D data.  The radius of the 4D data is the
       radius of the smallest enclosing sphere, centered at the To point.
       Note that during the loop through the nodes, Data4Radius holds
       the squared radius value.
    */

    Data4Radius = 0.0;
    for (i = 0; i < top->node.nel; i++)
    {
      Node_t v = OrgV(top->node.el[i]);
        ??? Temp4 = r4_Sub(c[v->num];
      with (o->To4),
      double dist = r4_dot(Temp4,Temp4)
    ){
        if (dist > Data4Radius)
        {
          Data4Radius = dist;
        }
      }
    }          

    Data4Radius = sqrt(Data4Radius);

    fprintf(stderr,"Data4Radius: "&
       Fmt.Pad(Fmt.LongReal(Data4Radius, Fmt.Style.Fix, prec = 4),8)&"\n");
    CalcV4Matrix();
    if (o->perspective)
    {
      pconst = 1.0 / Tan2Vangle4;  
    }
    else
    {
      rtemp = 1.0 / Data4Radius;
    }

    for (i = 0; i < top->node.nel; i++)
    {
      /* Transform the nodes from 4d World coordinates to 4D eye
         coordinates.
      */
      Node_t v = OrgV(top->node.el[i]);
      {
        TempV = r4_Sub(c[v->num],o->From4);
        depth = r4_dot(TempV,Wd);
        if (o->perspective)
        {
          rtemp = pconst / depth;
        }  
        c3[v->num][0] = rtemp * r4_dot(TempV, Wa);
        c3[v->num][1] = rtemp * r4_dot(TempV, Wb);
        c3[v->num][2] = rtemp * r4_dot(TempV, Wc);     
      }
    }           
  } /* END ProjectTo3D */
   
PROCEDURE SelectProjection(
    Options_t *o = (Options_t *)malloc(sizeof(Options_t));
    Coords_t *c, ElemTableRec_t *top)
      norm: r4_t = (r4_t){0.0, ..};
  { for (i = 0; i < top->cell.nel; i++){
      ??? f = Srot(top.cell[i]);
        ??? k = Triangulation.TetraNegNodes(f);
      ??? p = c[k[0]->num];
      ??? q = c[k[1]->num];
      ??? r = c[k[2]->num];
      ??? s = c[k[3]->num];
      ??? pq = r4_Sub(q;
      with ( p),
      double pr = r4_Sub(r, p);
      double ps = r4_Sub(s, p);
      double v = r4_cross(pq, pr, ps);
      double n = r4_Dir(v)
    ){
        norm = r4_Add(norm, n)
      }
    }
    ??? m = r4_Norm(norm);
    {
      if (m < 1.0d-20)
      {
        norm = (r4_t){1.0, 0.0, ..}
      }
      else
      {
        norm = r4_Scale(1.0/m, norm)
      }
    }
    ??? bar = Barycenter(top;
    with (c,TRUE)
  ){
      o->To4   = bar;
      o->From4 = r4_Add(o->To4, norm);
      SelectTwoIndepDirs(norm, o->Up4, o->Over4);
      o->Vangle4 = 120.0;
    }
  } /* END SelectProjection */

PROCEDURE SelectTwoIndepDirs(
    *u: r4_t;
    VAR v, w: r4_t;
 )
  /*
    Selects two vectors "v", "w", independent of each other
    and of the given vector "u". */
      m: uint = 0;
  { /* Find the largest coordinate of "u": */
    for (i = 1; i < 4; i++){ if ((fabs(u[i]) > fabs(u[m]))){ m = i; } }
    for (i = 0; i < 4; i++){ v[i] = 0.0; w[i] = 0.0; }
    v[(m+1) % 4] = 1.0;
    w[(m+2) % 4] = 1.0;
  } /* END SelectTwoIndepDirs */

Options_t GetOptions ()
  { Options_t *o = (Options_t *)malloc(sizeof(Options_t));
    argparser_t *pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);
    
    argparser_get_keyword(pp, "-inFileTp");                               
        o->inFileTp = argparser_get_next(pp);

        argparser_get_keyword(pp, "-inFileSt");                               
        o->inFileSt = argparser_get_next(pp);  
        
        if ((argparser_keyword_present(pp, "-outFile")))
        {
          o->outFile =  argparser_get_next(pp)
        }
        else
        {
          o->outFile =  o->inFileTp
        }

        o->fade   = argparser_keyword_present(pp, "-fade");
        o->wire   = argparser_keyword_present(pp, "-wire");
        o->detail = argparser_keyword_present(pp, "-detail");
        o->filter = argparser_keyword_present(pp, "-filter");
        
        if ((argparser_keyword_present(pp, "-silhouette")))
        {
          o->silhouette = TRUE;
        }

        argparser_get_keyword(pp, "-projection"); 
        ??? pname =  argparser_get_next(pp);
    {
          if (0 == strcmp(pname, "parallel")))
          {
            o->perspective = FALSE;
          }
          else
           if (0 == strcmp(pname, "perspective"))){
            o->perspective = TRUE;
          }
          else
          {
            argparser_error(pp, "Bad projection \"" & pname & "\"\n");
          }
        }

        if ((argparser_keyword_present(pp, "-normalize")))
        {
          o->normalize = TRUE;
        }
        else
        {
          o->normalize = FALSE;
        };

        if ((argparser_keyword_present(pp, "-autoProject")))
        {
          /* Program chooses From4: */
          o->autoProject = TRUE;
          o->multiple = FALSE;
        }
        else
        {
          /* User may give "From4" */
          o->autoProject = FALSE;
          o->multiple = argparser_keyword_present(pp, "-multiple");
          if (o->multiple)
          {
            /* Multiple From4"s will be read from stdin */
            o->From4 = (r4_t){0.0,0.0,0.0,0.0}
          }
          else
          {
            /* User may give single From4 option, or take default */
            if ((argparser_keyword_present(pp, "-From4")))
            {
              for (j = 0; j < 4; j++)
              {
                o->From4[j] = pp.getNextLongReal(-100.0, 100.0);
              }
            }
            else
            {
              o->From4 = (r4_t){100.0,50.0,20.0,0.0}
            }
          }
        }

        if (! o->multiple)
        {
          argparser_get_keyword(pp, "-select");
          o->selectTxt = argparser_get_next(pp);
  
          if (0 == strcmp(o->selectTxt, "back")))
          {
            o->select = Side.Back;
          }
          else
           if (0 == strcmp(o->selectTxt, "front"))){
            o->select = Side.Front;
          }
          else
           if (0 == strcmp(o->selectTxt, "both"))){
            o->select = Side.Both;
          }
          else
          {
            argparser_error(pp, "Bad shape \"" & o->selectTxt & "\"\n");
          }
        };
   
        if ((argparser_keyword_present(pp, "-To4")))
        {
          for (j = 0; j < 4; j++)
          {
            o->To4[j] = pp.getNextLongReal(-100.0, 100.0);
          }
        }
        else
        {
          o->To4 = (r4_t){0.0,0.0,0.0,0.0};
        }
   
        if ((argparser_keyword_present(pp, "-Up4")))
        {
          for (j = 0; j < 4; j++)
          {
            o->Up4[j] = pp.getNextLongReal(-100.0, 100.0);
          }
        }
        else
        {
          o->Up4 = (r4_t){0.0,1.0,0.0,0.0};
        }
        
        if ((argparser_keyword_present(pp, "-Over4")))
        {
          for (j = 0; j < 4; j++)
          {
            o->Over4[j] = pp.getNextLongReal(-100.0, 100.0);
          }
        }
        else
        {
          o->Over4 = (r4_t){0.0,0.0,1.0,0.0};
        }
        
        if ((argparser_keyword_present(pp, "-Vangle4")))
        {
          o->Vangle4 = pp.getNextLongReal(1.0, 179.0);
        }
        else
        {
          o->Vangle4 = 25.0;
        };
        if ((argparser_keyword_present(pp, "-colored")))
        {
          o->colored = TRUE;
        }
        else
        {
          o->colored = FALSE;
        }
        argparser_finish(pp);                                       
      
  -----------------------------------
#define _HELP \
  fprintf(stderr, "Usage: Visibility \\\n" \
   "        -inFileTp <name> \\\n" \
   "        -inFileSt <name> -outFile <name> \\\n" \
   "        -select [back | front | both] [-fade]\\\n" \
   "        [ -wire ] [ -detail ] [-filter ] \\\n" \
   "        [ -silhouette ] [-normalize <num> ] \\\n" \
   "        -projection [ perspective | parallel ]\\\n" \
   "        [ [ -autoProject ] | \\\n" \
   "          [ -From4 <num> <num> <num> <num> ] \\\n" \
   "          [ -To4   <num> <num> <num> <num> ] \\\n" \
   "          [ -Up4   <num> <num> <num> <num> ] \\\n" \
   "          [ -Over4 <num> <num> <num> <num> ] \\\n" \
   "        ] \\\n" \
   "        [ -Vangle4 <num> ] [ -colored ]\n");
          
END¦
      }
    }
    return o;
  } /* END GetOptions */

nil
{
  DoIt();
} Visibility.


  /* Copyright © 2001     Universidade Estadual de Campinas (UNICAMP)     */