#define PROG_NAME "ProjectTo3D"
#define PROG_DESC "???"
#define PROG_VERS "1.0"

#define ProjectTo3D_C_COPYRIGHT \
  ""
  
#define PROG_INFO \
  "" \
  " "


/* Projects a 4D configuration to 3D using the 4D viewing
   parameters as expected by the "Wire4" - Interactive 4D Wireframe Display
   Program.*/

#define ProjectTo3D_C_author \
  "Modified by L.A.P.Lozada on 2000-02-25."

#include <Stat.h>
#include <Tridimensional.h>
#include <math.h>
#include <Triangulation.h>
#include <argparser.h>
#define _GNU_SOURCE
#include <stdio.h>
#include <string.h>
#include <Mis.h>
#include <r4.h>
#include <Octf.h>

CONST
double Epsilon = 0.0000000001;

TYPE
double Coords4D = Triangulation.Coords_t;
double Coords3D_t = Tridimensional.Coords3D_t;

typedef struct Options_t {
      char *inFileTp;
      char *inFileSt;
      char *outFile;
      bool_t oblique;       /* TRUE == sets From4 at specified angle from mean normal. */
      double obliqueAngle; /* (if "oblique") Angle between From4 and mean normal, or . */
      bool_t adjustCenter;  /* TRUE == translates the model to the origin before projection. */
      bool_t adjustScale;   /* TRUE == scales model so that it fits in a sphere of radius 1. */
      bool_t perspective;   /* TRUE == uses perspective projection. */
      bool_t statistics;    /* TRUE == shows model statistics. */
      bool_t printDepth;    /* TRUE == prints 4D depth of every node. */
      /* 4D viewing parameters */
      From4:   r4_t;         /* Observer position in R^4. */
      To4:     r4_t;         /* Center of attention in R^4. */
      Up4:     r4_t;         /* Z coordinate of projection space. */
      Over4:   r4_t;         /* Y coordinate of projection space. */
      Vangle4: double;      /* Viewing angle. */
  }



PROCEDURE WL(double x) : char *==
  { return(Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, prec = 4), 8));
  } /* END WL */

PROCEDURE WL4(*x: r4_t) : char *==
  { return WL(x[0]) & " " & WL(x[1]) & " " & WL(x[2]) & " "  & WL(x[3]);
  } /* END WL4 */

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); */
    char *newcmt;
  { fprintf(stderr, "--------------------------------------------------\n");
    ??? tc = Triangulation.ReadToMa(o->inFileTp);
    ??? top = tc.top;
    ??? rc = Triangulation.ReadState(o->inFileSt);
    with ( c == rc^,
    double rc3 = NEW(REF Tridimensional.Coords3D_t, top->node.nel), c3 == rc3^;
    double cmt = tc.cmt & "\nProjectTo3D on " & Today() & "\noptions: " \
        "\n  -outFile " & o->outFile \
        "\n  -projection " & ARRAY bool_t OF char *{"parallel", "perspective"}[o->perspective],
    double bar = Triangulation.Barycenter(top, c, all = FALSE);
    double norm = MeanNorm(top, c);
    double meanRadius = MeanRadius(top, c, bar);
    double meanThickness = MeanThickness(top, c, bar, norm)
  ){
      fprintf(stderr, "Model shape parameters (before normalization and projection):\n");
      fprintf(stderr, "  normal ==          ( " & WL4(norm) & ")\n");
      fprintf(stderr, "  mean radius ==    " & WL(meanRadius) & "\n");
      fprintf(stderr, "  mean thickness == " & WL(meanThickness) & "\n");
      fprintf(stderr, "  rel. thickness == " & WL(meanThickness/meanRadius) & "\n");
      
      if (o->oblique)
      {
        /* Places "From4" so that "From4-To4" is oblique to the mean normal */
        SelectObliqueProjection(o, norm, o->obliqueAngle);
        newcmt = cmt & "\n   -oblique " & WL(o->obliqueAngle);
      }
      if (o->adjustCenter)
      {
        /* Move figure so that barycenter is at origin */
        for (i = 0 TO (c.nel-1)){ c[i] = r4_Sub(c[i], bar); }
        newcmt = cmt & "\n  -adjustCenter";
      }
      if (o->adjustScale)
      {
        /* Adjust scale of figure so that it fits in a sphere of radius 1 centered at To4 */
        ??? maxRadius = MaxRadius(top, c, o->To4);
    {
          for (i = 0 TO (c.nel-1)){ c[i] = r4_Scale(1.0/maxRadius, c[i]);}
        }
        newcmt = cmt & "\n  -adjustScale";
      }
      newcmt = cmt & "\nparameters: " \
        "\n  -From4 " & WL4(o->From4) \
        "\n  -To4   " & WL4(o->To4) \
        "\n  -Up4   " & WL4(o->Up4) \
        "\n  -Over4 " & WL4(o->Over4) \
        "\n  -Vangle4 " & WL(o->Vangle4);
      ProjectTo3D(o, c, c3, top);
      if (o->statistics)
      {
        ShowStatistics(top, c, c3)
      }
      WriteState3D(o->outFile, top, c3, newcmt)
    }
    fprintf(stderr, "--------------------------------------------------\n")
  } /* END DoIt */

double MaxRadius(ElemTableRec_t *top, *c: Coords4D; *ctr: r4_t)
       radius = 0.0;
  { for (i = 0; i < top->node.nel; i++){
      assert(OrgV(top->node[i])->num == i);
      ??? d = r4_Dist(ctr;
        with ( c[i])
    ){
        if (d > radius){ radius = d;}
      }
    }            
    return radius;
  } /* END MaxRadius */

double MeanRadius(ElemTableRec_t *top, *c: Coords4D; *bar: r4_t)
      sum2 : double = 0.0;
  { for (i = 0; i < top->node.nel; i++){
      ??? di = r4_Dist(c[i];
        with (bar)
    ){
        sum2 = sum2 + di*di;
      }
    }
    return sqrt(sum2/FLOAT(top->node.nel-1,double));
  } /* END MeanRadius */

double MeanThickness(ElemTableRec_t *top, *c: Coords4D; *bar, norm: r4_t)
      sum2 : double = 0.0;
  { for (i = 0; i < top->node.nel; i++){
      ??? di = r4_dot(r4_Sub(c[i];
        with (bar),norm)
    ){
        sum2 = sum2 + di*di;
      }
    }
    return sqrt(sum2/FLOAT(top->node.nel-1,double));
  } /* END MeanThickness */

r4_t MeanNorm(ElemTableRec_t *top, *c: Coords4D)
      norm: r4_t = (r4_t){0.0, ..};
  { for (i = 0; i < top->cell.nel; i++){
      ??? f = Srot(top.cell[i]);
        ??? k = TetrahedronNodes(f);
      ??? p = c[k[0]];
      ??? q = c[k[1]];
      ??? r = c[k[2]];
      ??? s = c[k[3]];
      ??? 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)
      {
        return (r4_t){1.0, 0.0, ..};
      }
      else
      {
        return r4_Scale(1.0/m, norm);
      }
    }
  } /* END MeanNorm */

PROCEDURE ProjectTo3D(
    Options_t *o = (Options_t *)malloc(sizeof(Options_t));
    *c: Coords4D;
    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.
    */

    double *norm;
    { 
      /* 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
    double rtemp;
    Wa,Wb,Wc,Wd : r4_t;
  { ??? angle = o->Vangle4/2.0;
    ??? angler = (FLOAT(Math.Pi, double)*angle)/180.0);
  with(
    double tan2vangle4 = Math.tan(angler);
    double pconst = 1.0 / tan2vangle4
  ){  

      CalcV4Matrix();
  
      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[i]);
        ??? t = r4_Sub(c[v->num];
        with (o->From4),
        double depth = r4_dot(t,Wd)
        DO
          if (o->printDepth)
          {
             fprintf(stderr, "v[" & Fmt.Pad(Fmt.Int(v->num),4) & "] depth == " \
                  Fmt.Pad(Fmt.LongReal(depth, Fmt.Style.Fix, prec = 4),8) & "\n"
            )
          }
          if (o->perspective){ rtemp = pconst / depth }else{ rtemp = 1.0; }
          c3[v->num][0] = rtemp * r4_dot(t, Wa);
          c3[v->num][1] = rtemp * r4_dot(t, Wb);
          c3[v->num][2] = rtemp * r4_dot(t, Wc); 
        }
      }          
    }
  } /* END ProjectTo3D */

PROCEDURE SelectObliqueProjection(
    Options_t *o = (Options_t *)malloc(sizeof(Options_t));
    *norm: r4_t;
    double angle;
 )
  CONST Pi == 3.14159265358979323844;
  { ??? rad = Pi/180.0 * angle;
    ??? c = Math.cos(rad);
    ??? s = Math.sin(rad);
    ??? dir = r4_Dir(r4_Mix(c;
    with ( norm, s, Perp(norm)))
  ){
      o->From4 = r4_Add(o->To4, r4_Scale(6.0, dir))
    }
    SelectTwoIndepDirs(norm, o->Up4, o->Over4);
  } /* END SelectObliqueProjection */

r4_t Perp(v: r4_t)
  /* Returns a vector perpendicular to "v" */
  uint *r; double s = 0.0;
  VAR { /* Find largest coordinate */
    r = 0;
    for (i = 1; i < 4; i++){ if ((fabs(v[i]) > fabs(v[r]))){ r = i; } }
    if ((fabs(v[r]) > 0.0))
    {
      /* Compute new value for it */
      for (i = 0; i < 4; i++){ if (i != r){ s = s + v[i]*v[i]; } }
      v[r] = - s / v[r];
    }
    return v;
  } /* END Perp */

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 */

PROCEDURE ShowStatistics(ElemTableRec_t *top; *c: Coords4D; *c3: Coords3D_t)
  { Analyze@{Edge->?}s(top, c);
    AnalyzeWalls(top, c, c3);
    AnalyzeCells(top, c, c3);
  } /* END ShowStatistics */

PROCEDURE Analyze@{Edge->?}s(ElemTableRec_t *top, *c: Coords4D) ==
  
      stl: Stat.T;
  { fprintf(stderr, "Statistics for existing @{edge->?}s:\n");
    Stat.Init(stl);  
    for (i = 0; i < top->NE; i++)
    {
      with (e  == top->edge[i])
      {
        if (e->exists)
        {
          ??? p = e.pa;
            Node_t un = OrgV(p)->num;
          Node_t vn = OrgV(Clock(p))->num;
          ??? l = r4_Dist(c[un];
          with (c[vn]),
          double ll = FLOAT(l, REAL)
        ){
            Stat.Accum(stl, ll)
          }
        }
      }
    }
    fprintf(stderr, "  lengths in R^4: ");
    Stat.Print(stderr, stl);
    fprintf(stderr, "\n");
  } /* END Analyze@{Edge->?}s */

PROCEDURE AnalyzeWalls(ElemTableRec_t *top, *c: Coords4D; *c3: Coords3D_t)
  
      sta: Stat.T;
  { fprintf(stderr, "Statistics for existing walls:\n");
    Stat.Init(sta);
    for (i = 0; i < top->wall.nel; i++)
    {
      ??? f = top->wall[i];
    {
        if (f->exists)
        {
          ??? p = f.pa,
          double un = OrgV(p)->num;
          double vn = OrgV(NextE(p))->num;
          double wn = OrgV(PrevE(p))->num;
          double a = r4_Dist(c[vn],c[un]);
          double b = r4_Dist(c[wn],c[un]);
          double c = r4_Dist(c[vn],c[wn]);
          double p = (a+b+c);
          double s = 0.5 * p;
          double sa = s - a;
          double sb = s - b;
          double sc = s - c;
          double ar = sqrt(MAX(s*sa*sb*sc, 0.0));
          double arr = FLOAT(ar, REAL)
;
    {
            Stat.Accum(sta, arr)
          }
        }
      }
    }
    fprintf(stderr, "  areas in R^4: ");
    Stat.Print(stderr, sta);
    fprintf(stderr, "\n");
    ??? ns = CountSilhouetteWalls(top;
    with (c3),
    double nb = CountBorderWalls(top);
    double nn = top->wall.nel-ns-nb
  ){
      fprintf(stderr, "  projection to R^3:" \
        " silhouette == " & Fmt.Int(ns) \
        " border == " & Fmt.Int(nb) \
        " normal == " &  Fmt.Int(nn) & "\n");
    }
  } /* END AnalyzeWalls */ 

PROCEDURE AnalyzeCells(ElemTableRec_t *top, *c: Coords4D; *c3: Coords3D_t)
      stv: Stat.T;
  
  { fprintf(stderr, "Statistics for existing tetrahedra:\n");
    Stat.Init(stv); 
    for (i = 0; i < top->cell.nel; i++)
    {
      ??? r = Srot(top.cell[i]), t == PnegP(Tors(r));
    {
        if (t->exists)
        {
          ??? vt = TetrahedronNodes(r);
            ??? un = vt[0];
          ??? vn = vt[1];
          ??? wn = vt[2];
          ??? xn = vt[3];
          ??? uv = r4_Sub(c[vn];
          with (c[un]),
          double uw = r4_Sub(c[wn],c[un]);
          double ux = r4_Sub(c[xn],c[un]);
          double n = r4_cross(uv, uw, ux);
          double v = 1.0/6.0 * r4_Norm(n);
          double vv = FLOAT(v, REAL)
        ){
            Stat.Accum(stv, vv)
          }
        }
      }
    }
      
    fprintf(stderr, "  volume in R^4: ");
    Stat.Print(stderr, stv);
    fprintf(stderr, "\n");
    ??? nn = CountNegativeTetrahedra(top;
    with (c3),
    double np = top->cell.nel - nn
  ){      
      fprintf(stderr, "  projection to R^3:" \
        " negative == " & Fmt.Int(nn) \
        " positive == " &  Fmt.Int(np) & "\n");
    }
  } /* END AnalyzeCells */
  /*
PROCEDURE CountBorder@{Edge->?}s(ElemTableRec_t *top) : uint ==
      count : uint = 0;
  { for (i = 0; i < top.edge.nel; i++){
      if ((Triangulation.EdgeIsBorder(top.edge[i].pa))){ count++; }
    }
    return count;
  }  CountBorder@{Edge->?}s;

PROCEDURE CountMiswound@{Edge->?}s(ElemTableRec_t *top; *c3: Coords3D_t) : uint ==
      count: uint = 0;
  { for (i = 0; i < top.edge.nel; i++){
      ??? a = top.edge[i].pa;
    {
        if ((! Triangulation.EdgeIsBorder(a)))
        {
          ??? wnd = Tridimensional.edgeWindingNumber(a, c3);
    {
            if ((fabs(wnd)!=1))
            {
              fprintf(stderr, "  @{edge->?} " & Fmt.Int(i)); 
              fprintf(stderr, " winding number == " & Fmt.Int(wnd) & "\n"); 
              count++
            }
          }
        }
      }
    }
    return count;
  }  CountMiswound@{Edge->?}s;
*/
uint CountSilhouetteWalls(ElemTableRec_t *top; *c3: Coords3D_t)

  bool_t WallIsSilhouette(Place_t @p)
    /* Return TRUE iff the wall associated to the @place "a" is a silhouette wall */
    {
      if ((PposP(a) == NULL) || (PnegP(a) == NULL)){ return FALSE; }
      ??? t = TetraNegPosNodes(a);
      ??? un = t[0]->num;
      ??? vn = t[1]->num;
      ??? wn = t[2]->num;
      ??? xn = t[3]->num;
      ??? yn = t[4]->num;
      ??? d1 = TetraDet3D(un;
      with (vn,wn,xn, c3),
      double d2 = TetraDet3D(un,vn,wn,yn, c3)
    ){  
        return d1*d2 >= 0.0; 
      }
    } WallIsSilhouette;
    
      count : uint = 0;
  { for (i = 0; i < top.wall.nel; i++){
      if ((WallIsSilhouette(top.wall[i].pa))){ count++; }
    }
    return count;
  }  CountSilhouetteWalls;

uint CountBorderWalls(ElemTableRec_t *top)

  bool_t WallIsBorder(Place_t @p)
    /* Return TRUE iff the wall associated to the @place "a" is a border 
       wall, FALSE c.c.
    */
    {
      return PposP(a) == NULL) || (PnegP(a) == NULL;
    } WallIsBorder;
    
      count : uint = 0;
  { for (i = 0; i < top.wall.nel; i++){
      if ((WallIsBorder(top.wall[i].pa))){ count++; }
    }
    return count;
  }  CountBorderWalls;

uint CountNegativeTetrahedra(ElemTableRec_t *top; *c3: Coords3D_t)

      count : uint = 0;
  { ??? t = Triangulation.CollectTetrahedra(top)^;
    {
      for (i = 0; i < top.cell.nel; i++)
      {
        ??? a = t[i];
        Node_t un = OrgV(a)->num;
        Node_t vn = OrgV(NextE(a))->num;
        Node_t wn = OrgV(PrevE(a))->num;
        Node_t xn = OrgV(PrevE(PrevF(a)))->num
        DO     
          if ((TetraDet3D(un;
        with (vn,wn,xn, c3) <= 0.0)){ count++; }
        }
      }
    }  
    return count;
  }  CountNegativeTetrahedra; 

PROCEDURE TetrahedronNodes(f:TriangulationPlace_t): ARRAY [0..3] OF uint ==
  /* Returns the nodes of a cell, given a dual wedge "f"
     with origin at the cell's center.  Assumes that the cell is
     a tetrahedron. */
  { Place_t g = Tors(f);
    Place_t h = Tors(Clock(PrevE(f)));
    ??? p = OrgV(g)->num;
    ??? q = OrgV(NextE(g))->num;
    ??? r = OrgV(PrevE(g))->num;
    ??? s = OrgV(PrevE(h))->num;
    {
      return ARRAY [0..3] OF uint {p, q, r, s};
    }
  } /* END TetrahedronNodes */

double TetraDet3D(uint u, uint v, uint w, uint x, Coords3D_t *c3)
  /*
    For each tetrahedron with nodes numbers u,v,w,x computes
    its orientation in R^{3} through the 4x4 determinant:
            _                                          _
           | c3[u][0]  c3[u][1]   c3[u][2]  1.0 |
     double B = | c3[v][0]  c3[v][1]   c3[v][2]  1.0 |
           | c3[w][0]  c3[w][1]   c3[w][2]  1.0 |
           | c3[x][0]  c3[x][1]   c3[x][2]  1.0 |
            -                                          - */
  { with (
    double a = (r4_t){c3[u][0], c3[u][1], c3[u][2], 1.0};
    double b = (r4_t){c3[v][0], c3[v][1], c3[v][2], 1.0};
    double c = (r4_t){c3[w][0], c3[w][1], c3[w][2], 1.0};
    double d = (r4_t){c3[x][0], c3[x][1], c3[x][2], 1.0}
  ){
      return r4_det(a,b,c,d);
    }   
  } /* END TetraDet3D */

Options_t *GetOptions(int argc, char **argv)
  { 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
        }
        
        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 \"" & argparser_get_next(pp) & "\"\n")
          }
        }

        o->adjustCenter = argparser_keyword_present(pp, "-adjustCenter");
        o->adjustScale = argparser_keyword_present(pp, "-adjustScale");

        if ((argparser_keyword_present(pp, "-oblique")))
        {
          o->oblique = TRUE;
          o->obliqueAngle = pp.getNextLongReal(-180.0, +180.0);
        }
        else
        {
          o->oblique = FALSE;
          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){0.0,0.0,0.0,-3.0};
          }
  
        }

        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;
        }

        o->statistics = argparser_keyword_present(pp, "-statistics");
        o->printDepth = argparser_keyword_present(pp, "-printDepth");

        argparser_finish(pp);                                       
      
  -----------------------------------
#define _HELP \
  fprintf(stderr, "Usage:\n" \
   "  ProjectTo3D -inFileTp <name> \\\n" \
   "    -inFileSt <name> [ -outFile <name> ] \\\n" \
   "    -projection [ Perspective | Parallel ]\\\n" \
   "    [ -adjustCenter ] [ -adjustScale ] \\\n" \
   "    [ -oblique ANGLE | -From4 <num> <num> <num> <num> ] \\\n" \
   "    [ -To4 <num> <num> <num> <num> ] \\\n" \
   "    [ -Up4 <num> <num> <num> <num> ] \\\n" \
   "    [ -Over4 <num> <num> <num> <num> ] \\\n" \
   "    [ -Vangle4 <num> ] \\\n" \
   "    [ -statistics ] \\\n" \
   "    [ -printDepth ] \n");
          
END¦
      }
    }
    return o;
  } /* END GetOptions */

/* end ProjectTo3D */


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