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

#define Depthcueing4D_C_COPYRIGHT \
  ""
  
#define PROG_INFO \
  "" \
  " "

/*
   This program performs the 4D depth-cueing by POVray tracer.
   The include file must be process with the "povray-3.1a-1" version
   in ~stolfi/pkg/povray-3.1a-1/PUB/sun4-SunOS-5.5

.

   Revisions:

   09-05-2000: Added in the rendering the nodes (spheres)
   20-05-2000: Added the silhouette walls.
   19-07-2000: Added the original walls.
*/

#include <r3.h>
#include <r3.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>
#include <Pov.h>

r4_t bari ;
  cmap : REF ARRAY OF Row3I;
  Wa,Wb,Wc,Wd : r4_t;
  vstack = NEW(REF uint_vec_t , INIT_STACK_SIZE); /* stack for nodes */
  vtop : uint = 0;                               /* top for the above stack */

CONST
double Epsilon = 0.0000000001;
double INIT_STACK_SIZE = 100000;

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

TYPE
double Row3I = ARRAY[0..2] OF INTEGER;
double FourNodes_t = RECORD u, v, w, x: uint; }
typedef struct Options_t {
    char *inFileTp;
    char *inFileSt;
    char *outFile;
    char *projectionName;
    bool_t autoProject;
    double normalize;   /* Normalize al nodes onto the S^3 with that radius. */
    From4:   r4_t;        /*                                           */
    To4:     r4_t;        /* 4D viewing parameters as expected by the  */
    Up4:     r4_t;        /* "Wire4"- Interactive 4D Wireframe Display */
    Over4:   r4_t;        /* Program.                                  */
    Vangle4: double;     /*                                           */   
    bool_t printDepth;
    bool_t colormap;
    dc_ColorFar:  Row3I;  
    dc_ColorNear: Row3I;
    double dc_DepthFar;
    double dc_DepthNear;
    INTEGER dc_Levels;
    bool_t silhouette;  /* draws the silhouette walls   */
    color      : frgb_t;     /* attributes of color and opacity for silhouette   */
    REAL opacity;     /* walls; Transparent (opacity==1) Opaque (opacity==0)*/
    bool_t walls;  /* draws the original walls     */
  }

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

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); */
    char *newcmt;
  { ??? tc = Triangulation.ReadToMa(o->inFileTp);
    ??? top = tc.top;
    ??? rc = Triangulation.ReadState(o->inFileSt);
    ??? c = rc^;
    ??? rc3 = NEW(REF Tridimensional.Coord3D;
    with ( top->node.nel),
    double c3 = rc3^;
    double rdepth =  double_vec_new(top->node.nel);
    double depth = rdepth^;
    double cmt = tc.cmt & "\nProjected in " & o->projectionName \
                 ": " & o->outFile & ".st3  on " & Today()
  ){
      cmap = NEW(REF ARRAY OF Row3I, o->dc_Levels);
      SetColorMap(o);
      if (o->colormap){ PrintColorMap(o); }
      if (o->normalize > 0.0)
      {
        fprintf(stderr, "projecting nodes onto the unit S^3\n");
        NormalizeNodeCoords(c, o->normalize);
      }
      if (o->autoProject
      ){
          SelectProjection(o, c, top);
          newcmt = cmt & "\nWith AutoProject";
        }
        else
        {
          newcmt = cmt & "\nParameters: " \
                         "\nFrom4: " & WL(o->From4[0]) & WL(o->From4[1]) 
                                     & WL(o->From4[2]) & WL(o->From4[3]) \
                         "\nTo4:   " & WL(o->To4[0]) & WL(o->To4[1]) 
                                     & WL(o->To4[2]) & WL(o->To4[3]) \
                         "\nUp4:   " & WL(o->Up4[0]) & WL(o->Up4[1]) 
                                     & WL(o->Up4[2]) & WL(o->Up4[3]) \
                         "\nOver4: " & WL(o->Over4[0]) & WL(o->Over4[1]) 
                                     & WL(o->Over4[2]) & WL(o->Over4[3]) \
                         "\nVangle4: " & WL(o->Vangle4);
      }
      ProjectTo3D(o, c, c3, depth, top);
      if (0 == strcmp(o->projectionName, "Parallel")))
      {
        WritePOVFile(o, depth, o->outFile & "-Par", top, c3);
      }
      else
       if (0 == strcmp(o->projectionName, "Perspective"))){
        WritePOVFile(o, depth, o->outFile & "-Per", top, c3);
      }
      bari = Triangulation.Barycenter(top, c);  
    
    return 0;
  }

  void NormalizeNodeCoords(VAR c: Coords_t; newR: double)
    {
      ??? b = Barycenter(c);
    {
        for (i = 0 TO (c.nel-1))
        {
          ??? p = c[i], q == r4_Sub(p, b), r == r4_Norm(q);
    {
            if (r > 0.0){ p = r4_Scale(newR/r, q);}
          }
        }
      }
    } NormalizeNodeCoords;

  r4_t Barycenter(*c: Coords_t)
    r4_t B = (r4_t){0.0, ..};
    VAR {
      for (i = 0 TO (c.nel-1)){ B = r4_Add(B, c[i]); }
      return r4_Scale(1.0/FLOAT((c.nel), double), B);
    } Barycenter;

  void CalcV4Matrix(*o: Options_t)
    /* 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.
    */
    
    VAR
      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;  

PROCEDURE ProjectTo3D(
    Options_t *o = (Options_t *)malloc(sizeof(Options_t));
    Coords_t *c;
    VAR c3: Tridimensional.Coord3D;
    double_vec_t depth, ElemTableRec_t *top)
  
  VAR
    Tan2Vangle4, Data4Radius, pconst, rtemp: double;
    TempV : r4_t;
  { ??? 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\n");
    CalcV4Matrix(o);
    if (0 == strcmp(o->projectionName, "Parallel")))
    {
      rtemp = 1.0 / Data4Radius;
    }
    else
    {
      pconst = 1.0 / Tan2Vangle4;  
    }

    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[v->num] = r4_dot(TempV,Wd);
        if (0 == strcmp(o->projectionName,  "Perspective")))
        {
          rtemp = pconst / depth[v->num];
        };  
        if (o->printDepth)
        {
          fprintf(stderr, "v[" & Fmt.Pad(Fmt.Int(v->num),4) & "] depth == " \
           Fmt.Pad(Fmt.LongReal(depth[v->num],Fmt.Style.Fix,prec=4),8) & "\n");
        }
        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 TetrahedronNodes(f:TriangulationPlace_t): ARRAY [0..3] OF uint
  { Place_t g = Tors(f);
    Place_t h = Tors(Clock(PrevE(f)));
    Node_t p = OrgV(g)->num;
    Node_t q = OrgV(NextE(g))->num;
    Node_t r = OrgV(PrevE(g))->num;
    Node_t s = OrgV(PrevE(h))->num;
    {
      return ARRAY [0..3] OF uint {p, q, r, s};
    }
  } /* END TetrahedronNodes */
  
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 = 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)
      {
        norm = (r4_t){1.0, 0.0, ..}
      }
      else
      {
        norm = r4_Scale(1.0/m, norm)
      }
    }
    ??? bar = Triangulation.Barycenter(top;
    with ( c)
  ){
      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 */

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

void WritePOV(
    Options_t *op;
    *depth: double_vec_t;
    wr: Wr.T;
    *ElemTableRec_t *top; 
    *c3: Coord3D;
 )
  

  void WritePOVCylinder4DC(
      *o,d: r3_t;
      REAL radius;
      *cr: ARRAY [0..1] OF R3.T;
   )
    /*
      Defines a finity length cylinder without parallel end caps with
      4D depthcueing. This effect is obtained by the use of color map
      that modified the pigments colors through the length cylinder.
    */
    {
      fprintf(wr,"  cylinder {\n");
      fprintf(wr,"    "); WritePOVPoint(wr,o); fprintf(wr, ",\n");
      fprintf(wr,"    "); WritePOVPoint(wr,d); fprintf(wr, ",\n");
      fprintf(wr,"    " & Fmt.Real(radius) & "\n");
      fprintf(wr,"    open\n");
      fprintf(wr,"    pigment {\n");
      fprintf(wr,"      gradient ");  WritePOVPoint(wr,r3_sub(d,o));
                                         fprintf(wr,"\n");
      fprintf(wr,"      color_map{\n");
      fprintf(wr,"        [0.0 color "); WritePOVColor(wr,cr[0]);
                                            fprintf(wr,"]\n");
      fprintf(wr,"        [1.0 color "); WritePOVColor(wr,cr[1]);
                                            fprintf(wr,"]\n");
      fprintf(wr,"      }\n");
      fprintf(wr,"      translate "); WritePOVPoint(wr,o);
                                         fprintf(wr,"\n");
      fprintf(wr,"    }\n");
      fprintf(wr,"  }\n\n");
      fflush(wr);
    } WritePOVCylinder4DC;

  double FindOriR3(q: FourNodes_t)
    /*
      For each tetrahedron with extremus nodes numbers u,v,w,x compute us
      its orientation in R^{3} through the 4x4 determinant:
              _                                          _
             | c3[q.u][0]  c3[q.u][1]   c3[q.u][2]  1.0 |
       double B = | c3[q.v][0]  c3[q.v][1]   c3[q.v][2]  1.0 |
             | c3[q.w][0]  c3[q.w][1]   c3[q.w][2]  1.0 |
             | c3[q.x][0]  c3[q.x][1]   c3[q.x][2]  1.0 |
              -                                          -
    */
    {
      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);
      }   
    } FindOriR3;

  bool_t SilhouetteWalls(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 = FindOriR3(t1);
      double d2 = FindOriR3(t2)
    ){  
        if (((Sign(d1)) && (Sign(d2))) || (((NOT Sign(d1))) && ((NOT Sign(d2)))))
        {
          return TRUE;
        }
        else
        { return FALSE;
        }
      }
    } SilhouetteWalls;

      frgb_t cr[1+1]; 
  { for (i = 0; i < top->NE; i++){ 
      ??? e = top->edge[i];
      { 
        if (e->exists)
        {
          ??? a = e.pa;
             Node_t ov = OrgV(a);
          Node_t dv = OrgV(Clock(a));
          ??? or = ov->num;
          ??? de = dv->num;
          ??? cm = LineDepthCue(op;
          with (ov,dv,depth)
        ){
            for (i = 0; i < 2; i++)
            {
              for (j = 0; j < 3; j++)
              {
                cr[i][j] = FLOAT(cm[i][j],REAL)/255.0;  
             }                                         
            }                                         
            WritePOVCylinder4DC(c3[or], c3[de], e.radius, cr);
            ??? ov = (Node_t)(ov);
              {
            double dv = (Node_t)(dv)
          ){
              if ((! Present(vstack,vtop,or)))
              {
                WritePOVSphere(wr,c3[or], ov.radius, cr[0], 0.0, TRUE);
                Save(vstack,vtop,ov->num);
              }
              if ((! Present(vstack,vtop,de)))
              {
                WritePOVSphere(wr,c3[de], dv.radius, cr[1], 0.0, TRUE);
                Save(vstack,vtop,dv->num);
              }
            }
          }
        }
      }
    }

    if (op.silhouette)
    {
      for (i = 0; i < top->wall.nel; i++)
      {
        ??? f = top->wall[i];
        ??? a = f.pa;
        Node_t un = OrgV(a)->num;
        Node_t vn = OrgV(NextE(a))->num;
        Node_t wn = OrgV(PrevE(a))->num;
        ??? c = op.color;
        ??? t = op.opacity;
        {
          if ((SilhouetteWalls(a)) && (NOT f->exists))
          {
            WritePOVTriangle(wr,c3[un], c3[vn], c3[wn],c,t, TRUE);
          }
        }
      }
    }

    if (op.walls)
    {
      for (i = 0; i < top->wall.nel; i++)
      {
        ??? f = top->wall[i];
        ??? a = f.pa;
        Node_t un = OrgV(a)->num;
        Node_t vn = OrgV(NextE(a))->num;
        Node_t wn = OrgV(PrevE(a))->num;
        ??? t3 = f.transp;
        ??? tp = (t3[0] + t3[1] + t3[2]) / 3.0;
        ??? co = f.color;
        {
          if (f->exists)
          {
            WritePOVTriangle(wr,c3[un], c3[vn], c3[wn], co, tp, TRUE);
          }
        }
      }
    }
    fflush(wr);
  } /* END WritePOV */

PROCEDURE LineDepthCue(
  Options_t *o = (Options_t *)malloc(sizeof(Options_t));
  Node P,Q;
  *depth: double_vec_t;
 ) : ARRAY[0..1] OF Row3I ==
  /*
   This function computes the colors for 4D depthcueing.  The main approach used is
   to subdivide the pigment into dc_Levels subsegments, that is realized by the
   "gradient <vector>" implemented by J. Stolfi in the "povray-3.1a-1" version,
   where <vector> is an arbitrary vector. */
      Row3I cd[1+1];
  { ??? depthP = depth[P->num];
    ??? depthQ = depth[Q->num];
    ??? dc_DFar_DNear = o->dc_DepthFar - o->dc_DepthNear;
    ??? dscale = FLOAT(o->dc_Levels, double) / dc_DFar_DNear);
  with(  )
  {
       if (depthP < o->dc_DepthNear
       ){
          cd[0] = Row3I{cmap[0][0],cmap[0][1],cmap[0][2]};
       }
       else
        if (depthP > o->dc_DepthFar){
         cd[0] = Row3I{cmap[o->dc_Levels-1][0],
                        cmap[o->dc_Levels-1][1],
                        cmap[o->dc_Levels-1][2]};
       }
       else
       {
          cd[0] = Row3I{cmap[FLOOR(dscale*(depthP-o->dc_DepthNear))][0],
                         cmap[FLOOR(dscale*(depthP-o->dc_DepthNear))][1],
                         cmap[FLOOR(dscale*(depthP-o->dc_DepthNear))][2]};
       }        

       if (depthQ < o->dc_DepthNear
       ){ 
           cd[1] = Row3I{cmap[0][0],cmap[0][1],cmap[0][2]};
       }
       else
        if (depthQ > o->dc_DepthFar){
         cd[1] = Row3I{cmap[o->dc_Levels-1][0],
                        cmap[o->dc_Levels-1][1],
                        cmap[o->dc_Levels-1][2]};
       }
       else
       {
         cd[1] = Row3I{cmap[FLOOR(dscale*(depthQ-o->dc_DepthNear))][0],
                        cmap[FLOOR(dscale*(depthQ-o->dc_DepthNear))][1],
                        cmap[FLOOR(dscale*(depthQ-o->dc_DepthNear))][2]}; 
       }
       return cd;
     }  
   } LineDepthCue;

PROCEDURE SetColorMap(*o: Options_t)
  /* Set up the depthcue color map. */
  INTEGER *range;
  { for (rgb = 0; rgb < 3; rgb++){
      range = o->dc_ColorFar[rgb] - o->dc_ColorNear[rgb];
       
      if (range < 0)
      {
        range = range - 1;
      }
      else
      {
        range = range + 1;
      }        
      
      for (i = 0; i < o->dc_Levels; i++)
      {
       cmap[i][rgb]= o->dc_ColorNear[rgb]+ (i*range) DIV o->dc_Levels;
      }
    }                                 
  }  SetColorMap;

PROCEDURE PrintColorMap(*o: Options_t)
  /* Print the depthcue color map. */
  
  { for (i = 0; i < o->dc_Levels; i++){
      fprintf(stderr," cmap["&Fmt.Pad(Fmt.Int(i),2)&"] is");
      for (rgb = 0; rgb < 3; rgb++)
      {
        fprintf(stderr, Fmt.Pad(Fmt.Int(cmap[i][rgb]), 4));
      }
      fprintf(stderr,"\n");
    }                                 
  } /* END PrintColorMap */

PROCEDURE Save(
      Stack: REF uint_vec_t;
  uint *top,element;
 )
  /* Saves the "element" on the stack "Stack". */
  { Stack.El[top] = element;
    top++;
  } /* END Save */

bool_t Present(
  *Stack: REF uint_vec_t;
  uint top,element;
 )
  /* Return TRUE if "element" its on the stack, FALSE c.c. */
      nstack1: uint = top;
  { while (nstack1 > 0){
      nstack1 = nstack1 - 1;
      if (Stack.El[nstack1] == element){ return TRUE; }
    }
  return FALSE;
} Present;

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"); 
        o->projectionName =  argparser_get_next(pp);
        
        if ((! 0 == strcmp(o->projectionName, "Parallel") OR
                Text.Equal(o->projectionName, "Perspective"))
        THEN
          argparser_error(pp, "Bad projection \"" & argparser_get_next(pp) & "\"\n")
        }

        if ((argparser_keyword_present(pp, "-normalize")))
        {
          /* Desired radius of model in R^4 */
          o->normalize = pp.getNextLongReal(1.0d-10,1.0d+10);
        }
        else
        {
          o->normalize = 0.0 /* No normalization */
        }

        if ((argparser_keyword_present(pp, "-autoProject")))
        {
          o->autoProject = TRUE
        }
        else
        {
          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->printDepth = argparser_keyword_present(pp, "-printDepth");

        o->colormap   = argparser_keyword_present(pp, "-colormap");

        if ((argparser_keyword_present(pp, "-ColorFar")))
        {
          for (j = 0; j < 3; j++)
          {
            o->dc_ColorFar[j] = argparser_get_next_int(pp, 0,255);
          }
        }
        else
        {
          o->dc_ColorFar= Row3I{255, 255, 255};
        };

        if ((argparser_keyword_present(pp, "-ColorNear")))
        {
          for (j = 0; j < 3; j++)
          {
            o->dc_ColorNear[j] = argparser_get_next_int(pp, 0,255);
          }
        }
        else
        {
          o->dc_ColorNear= Row3I{0, 0, 0};
        };     
        
        if ((argparser_keyword_present(pp, "-DepthNear")))
        {
          o->dc_DepthNear = pp.getNextLongReal(-1000.0, 1000.0);
        }
        else
        {
          o->dc_DepthNear= 2.0;
        };     
        
        if ((argparser_keyword_present(pp, "-DepthFar")))
        {
          o->dc_DepthFar = pp.getNextLongReal(-1000.0, 1000.0);
        }
        else
        {
          o->dc_DepthFar = 6.0;
        };     
        
        if ((argparser_keyword_present(pp, "-DepthCueLevels")))
        {
          o->dc_Levels = argparser_get_next_int(pp, 0, 10000);
        }
        else
        {
          o->dc_Levels = 6;
        };     
        
        if ((argparser_keyword_present(pp, "-silhouette")))
        {
          o->silhouette = TRUE;
          if ((argparser_keyword_present(pp, "-color")))
          {
            for (j = 0; j < 3; j++)
            {
              o->color[j] = pp.getNextReal(0.0,1.0);
            }
          }
          else
          {
            o->color = (frgb_t){1.0,1.0,0.75};
          }
          if ((argparser_keyword_present(pp, "-opacity")))
          {
            o->opacity = pp.getNextReal(0.0,1.0);
          }
          else
          {
            o->opacity = 0.85;
          }
        }

        if ((argparser_keyword_present(pp, "-walls")))
        {
          o->walls = TRUE;
        }

        argparser_finish(pp);                                       
      
  -----------------------------------
#define _HELP \
  fprintf(stderr, "Usage:\n" \
   "  Depthcueing4D -inFileTp <name> \\\n" \
   "    -inFileSt <name> [ -outFile <name> ] \\\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" \
   "    ] [ -Vangle4 <num> ] \\\n" \
   "    [ -printDepth ] [ -normalize <num> ] [ -colormap ] \\\n" \
   "    [ -ColorFar  <num> <num> <num> ] \\\n" \
   "    [ -ColorNear <num> <num> <num> ] \\\n" \
   "    [ -DepthFar  <num> ] [ -DepthNear <num> ] \\\n" \
   "    [ -DepthCueLevels <num> ] \\\n" \
   "    [ -silhouette [ -color | -opacity ] ] \\\n" \
   "    [ -walls ]\n");
          
END¦
      }
    }
    return o;
  } /* END GetOptions */

/* end Depthcueing4D */


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