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

#define Explode_C_COPYRIGHT \
  ""
  
#define PROG_INFO \
  "" \
  " "

/* This program receives as input a single cell (the envelope) and produces
  a include file with the barycentric subdivision of walls (2-skeleton) and
  cells (3-skeleton) but in an exploding way. */

#define Explode_C_author \
  "Modified by L.A.P.Lozada on 2000-05-19."

#include <r3.h>
#include <string.h>
#define _GNU_SOURCE
#include <stdio.h>
#include <r4.h>
#include <Triangulation.h>
#include <Octf.h>
#include <argparser.h>
#include <Pov.h>

TYPE
typedef struct Options_t {                
    char *inFileTp;           /* Input file name (minus ".tp" extension) */
    char *inFileSt;           /* Input file name (minus ".st" extension) */
    char *outFile;            /* Output file name prefix */  
    char *projectionName;
    bool_t autoProject;
    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 explode;      /* For explode the configuration             */
    Element element;      /* Choose the element to explode (wall or cell)  */
    char *eleName;
  }                        

double Element = {Wall, Cell};
    
CONST
double Epsilon = 0.0000000001;  
  
Options_t *GetOptions(int argc, char **argv);

int main(int argc, char **argv)
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E);
  uint NNE = 0;
    NNV:  uint = 0;
    gi :  REF ARRAY OF Place_vec_t;
    ElemTableRec_t ntop;
    vf :  REF Node_vec_t;
    o   : Options_t = GetOptions();
    inc : FileWr.T;
  { ??? tc = Triangulation.ReadToMa(o->inFileTp);
    ??? top = tc.top;
    ??? rc = Triangulation.ReadState(o->inFileSt);
    ??? c =  rc^;
    ??? half = Place_vec_new(2*top->wall.nelE)^;
    ??? vnew = Node_vec_new(top->node.nel)^;
    {
      if (o->explode)
      {
        if (o->element == Element.Wall)
        {
          inc = FileWr.Open(o->outFile & "-Ex-Wall.inc");
        }
        else
         if ( o->element == Element.Cell){  
          inc = FileWr.Open(o->outFile & "-Ex-Cell.inc");
        }
      }            
    
      fprintf(stderr, "Subdividing from: " & o->inFileTp & ".tp\n");
      gi = NEW(REF ARRAY OF Place_vec_t, top->cell.nel, 4);
      Node_vec_t vf = Node_vec_new(top->wall.nel);

      Place_t Half(Place_t @p)
        {
          ??? na = PWedge(a)->num;
            ??? oa = OrientationBit(a);
          ??? sa = SpinBit(a);
          ??? h = half[2*na + oa];
          {
            if (sa == 0){ return h }else{ return Spin(h); }
          }
        } Half;
  
      bool_t CreateNode(Place_t @p)
        Place_t t = a;
        VAR {
          ??? fe = NARROW(PWedge(a)-> TriangulationWedge_t);
    {
            if (! fe.mark)
            {
              fe.mark = TRUE;
              do 
              {
                ??? fe = NARROW(PWedge(t), TriangulationWedge_t);
    { 
                  fe.mark = TRUE;  
                }
                t = NextF(t)
              } while (t != a);
              return TRUE;
            }
            return FALSE;
          }
        } CreateNode;

      void SetAllVhnum(Place_t @p, Node_t v)

        void SetVhnum(Place_t @p)
          {
            ??? fe = NARROW(PWedge(a)-> TriangulationWedge_t);
    { 
              fe.vh = v;  
            }
          } SetVhnum;

        Place_t t = a;
        VAR {
          do 
          {
            SetVhnum(t);
            t = NextF(t);
          } while (t != a);
        } SetAllVhnum;
        
      Node Vhnum(Place_t @p)
        {
          ??? fe = NARROW(PWedge(a)-> TriangulationWedge_t);
    {
            return fe.vh;
          }
        } Vhnum; 
            
      {
        /* Copy original nodes, save correspondence in "vnew" array: */
        for (iu = 0; iu < top->node.nel; iu++)
        {
          ??? u = top->node[iu];
            Node_t v = MakeNode(u->num);
          {
            v->exists = u->exists;
            v.fixed = u.fixed;
            v.color = u.color;
            v.radius = u.radius;
            vnew[iu] = v
          }
        }

        /* 
          Create two new wedges for each original wedge "fe".
          The new wedge corresponding to the origin half of "fe",
          with same spin and orientation, will be
          | Half(a) == Spin^s(half[2*PWedge(a)->num + s])
          where s == SpinBit(a), s == OrientationBit(a)
        */
        
      REF Node_vec_t  ve = Node_vec_new(top->NE);
        uint i = 0;
      {
        for (ie = 0; ie < top->wall.nelE; ie++)
        {
          ??? a = top->wedge[ie];
            ??? oa = OrientationBit(a);
          ??? sa = SpinBit(a);
          ??? ho = half[2*ie + oa];
          ??? hd = half[2*ie + 1 - oa];
          ??? fe = NARROW(PWedge(a)->
          with ( TriangulationWedge_t),
          double b = CreateNode(a)
        ){
            if (b)
            {
              ve[i] = MakeNode(top->node.nel + NNV); INC (NNV);
              ve[i]->exists = top->edge[fe.edge->num]->exists;
              ve[i].fixed = FALSE;
              ve[i].radius = fe.edge.radius;
              ve[i].color = fe.edge.color;
              SetAllVhnum(a, ve[i]);
              i++;
            }
            ho = MakeWedge();
            ??? hoe = NARROW(PWedge(ho)-> TriangulationWedge_t);
    {
              hoe->num = NNE; NNE++;
              hoe.edge->exists = fe.edge->exists;
              hoe.edge.radius = fe.edge.radius;
              hoe.edge.color = fe.edge.color;
            }
            
            if (sa == 1){ ho = Spin(ho); }
            hd = MakeWedge();
            ??? hde = NARROW(PWedge(hd)-> TriangulationWedge_t);
    {
              hde->num = NNE; NNE++;
              hde.edge->exists = fe.edge->exists;
              hde.edge.radius = fe.edge.radius;
              hde.edge.color = fe.edge.color;
            }
            
            if (sa == 1){ ho = Spin(ho); }
            SpliceEdges(ho, Clock(hd));
 
            ??? m = Vhnum(a);
    {
              SetOrg(Clock(ho),ve[m->num-top->node.nel]);
              SetOrg(Clock(hd),ve[m->num-top->node.nel]);
            }

            SetOrg(ho, vnew[OrgV(a)->num]);
            SetOrg(hd, vnew[OrgV(Clock(a))->num]);            
          }
        }
      }
                
      /* Connect the half-wedges as in the original triangulation */

      for (ie = 0; ie < top->wall.nelE; ie++)
      {
        ??? a = top->wedge[ie];
          Place_t b = NextF(a);
        Place_t c = NextE(a);
        ??? ha = Half(a);
        ??? hac = Half(Clock(a));
        ??? hb = Half(b);
        ??? hc = Half(c);
        {
          if ((b!=a) && (NextF(ha)!=hb))
          {
            SetNextF(ha, hb);             /* so, NextF(ha) == hb */
            SetEdgeInfo(hb, PEdge(ha));
          }

          if ((c!=a) && (NextE(Clock(hac))!=hc))
          {
            SetNextE(Clock(hac), hc);     /* so, NextE(Clock(hac) == hc */
          }
        }

        Place_t a = Clock(top->wedge[ie]);
          Place_t b = NextF(a);
        Place_t c = NextE(a);
        ??? ha = Half(a);
        ??? hac = Half(Clock(a));
        ??? hb = Half(b);
        ??? hc = Half(c);
        {
          if ((b!=a) && (NextF(ha)!=hb))
          {
            SetNextF(ha, hb);             /* so, NextF(ha) == hb */
            SetEdgeInfo(hb, PEdge(ha));
          }

          if ((NextE(Clock(hac))!=hc))
          {
            SetNextE(Clock(hac), hc);     /* so, NextE(Clock(hac) == hc */
          }
        }
      }

      for (i = 0; i < top->wall.nel; i++)
      {
        vf[i] = MakeNode(top->node.nel + NNV); INC (NNV);
      }

      for (i = 0; i < top->wall.nel; i++)
      {
        Place_t *a,b,c,d,aa;
        {
          ??? p = top->wall[i].pa;
    {
            aa = p;
            for (j = 0 TO Octf.DegreeOfEdge(p)-1)
            {
              ??? ha = Half(aa);
                ??? hac = Half(Clock(aa));
              {
                a = MakeWedge();
                PWedge(a)->num = NNE; NNE++;
                b = MakeWedge();
                PWedge(b)->num = NNE; NNE++;
                c = MakeWedge();
                PWedge(c)->num = NNE; NNE++;
                d = MakeWedge();
                PWedge(d)->num = NNE; NNE++;
              
                SetOrg(a, OrgV(Clock(ha)));
                SetOrg(Clock(a),vf[i]);
                SetOrg(b,vf[i]);
                SetOrg(Clock(b), OrgV(ha));
                SetNextE(ha,a);
                SetNextE(a,b);
                SetNextE(b,ha);
                SetRingWallInfo(a,PWall(a));

                SetOrg(c, OrgV(Clock(hac)));
                SetOrg(Clock(c),vf[i]);
                SetOrg(d,vf[i]);
                SetOrg(Clock(d), OrgV(hac));
                SetNextE(hac,c);
                SetNextE(c,d);
                SetNextE(d,hac);
                SetRingWallInfo(c,PWall(c));
                SetNextF(a,c);
                SetNextF(c,a);
                SetOrgAll(a, OrgV(a));
              }
              aa = NextE(aa);
            }
          }
        }
      }

      for (i = 0; i < top->wall.nel; i++)
      {
        Place_t *aa;
        {
          ??? p = top->wall[i].pa;
          {
            aa = p;
            for (j = 0 TO Octf.DegreeOfEdge(p)-1)
            {
              ??? hpc = Half(Clock(aa));
                ??? hq = Half(NextE(aa));
              {
                SetNextF(PrevE(hq), PrevE(hpc));
                SetRingEdgeInfo(PrevE(hq),PEdge(PrevE(hq)));
              }
              aa = NextE(aa);
            }
            SetOrgAll(PrevE(Half(p)),vf[i]);
          }
        }
      }

      for (i = 0; i < top->wall.nelE; i++)
      {
        ??? a = top->wedge[i];
        ??? ha = Half(a);
        ??? hac = Half(Clock(a))
        DO
          SetRingEdgeInfo(ha;
        with (PEdge(ha));
          SetRingEdgeInfo(hac, PEdge(hac));
          SetOrgAll(hac, OrgV(hac));
          SetOrgAll(ha, OrgV(ha));
          SetOrgAll(Clock(hac), OrgV(Clock(hac)));
          SetOrgAll(Clock(ha), OrgV(Clock(ha)));
        }
      } 
        
      for (i = 0; i < top->cell.nel; i++)
      {
        ??? v = Srot(top.cell[i]), a == Tors(v), af == PrevF(a),
           double ae = PrevE(a), aee == NextE(a)
;
    {
          assert(PnegP(a)->num == i);
          assert(OrgV(v)->num == i);
          if ((PposP(af)!=NULL)){ assert(PnegP(a) == PposP (af)); }
            
          gi[i,0] = Clock(NextE(PrevF(NextE(Half(a)))));
          gi[i,1] = Clock(NextE(NextF(NextE(Half(af)))));
          gi[i,2] = Clock(NextE(NextF(NextE(Half(PrevF(ae))))));
          gi[i,3] = NextF(NextE(Half(PrevF(aee))));
        }
      }

      ntop = MakeElemTable(Half(top->wedge[0]));        
        
      ??? nc = NEW(REF Coords_t;
        with (ntop.node.nel)^
    ){
        double Compute4DRadius()
          /* compute Data4Radius */
          double Data4Radius = 0.0;
          {
            for (i = 0; i < ntop.node.nel; i++)
            {
              ??? v = ntop.node[i];
              ??? Temp4 = r4_Sub(nc[v->num];
              with (o->To4),
              double dist = r4_dot(Temp4, Temp4)
            ){
                if (dist > Data4Radius)
                {
                  Data4Radius = dist;
                }
              }
            }
            return Data4Radius;
          }  Compute4DRadius;   

        void CalcV4Matrix(VAR WA,WB,WC,WD: r4_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.
          */
          double norm ;
          { 
            /* Calculate Wd, the 4th coordinate 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;  
    
        r3_t ProjectTo3D(p: r4_t)
          r3_t c3;
          VAR {
            ??? TempV = r4_Sub(p;
            with (o->From4),
            double rtemp = 1.0 / data4
          ){  
              c3[0] = rtemp * r4_dot(TempV, Wa);
              c3[1] = rtemp * r4_dot(TempV, Wb);
              c3[2] = rtemp * r4_dot(TempV, Wc); 
              return c3;
            }
          } ProjectTo3D;
        
        void ExplodeWall(un,vn,wn: uint)
          r4_t B = (r4_t){0.0, ..};
          VAR {
            ??? cun = nc[un];
            ??? cvn = nc[vn];
            ??? cwn = nc[wn];
            ??? us = r4_Scale(4.00;
            with (cun),
            double vs = r4_Scale(4.00,cvn);
            double ws = r4_Scale(4.00,cwn);
              b  == r4_Add(r4_Add(r4_Add(B, cun),cvn),cwn),
            double ba = r4_Scale(1.0/3.0, b);
            double ue = r4_Add(us, ba);
            double ve = r4_Add(vs, ba);
            double we = r4_Add(ws, ba);
            double ue3 = ProjectTo3D(ue);
            double ve3 = ProjectTo3D(ve);
            double we3 = ProjectTo3D(we);
              ra  == 0.008,
              cf  == (frgb_t){1.000, 0.745, 0.745},
              ce  == (frgb_t){0.0,0.0,0.0},
              tr  == 0.0
          ){
              /* Draw cylinders */
              WritePOVCylinder(inc,ue3,ve3,ra,ce,tr,TRUE);
              WritePOVCylinder(inc,ve3,we3,ra,ce,tr,TRUE);
              WritePOVCylinder(inc,ue3,we3,ra,ce,tr,TRUE);
              /* Draw wall */
              WritePOVTriangle(inc,ue3,ve3,we3,cf,tr,TRUE);
            }
          } ExplodeWall;

          void ExplodeTetrahedron(un,vn,wn: uint)
            r4_t B = (r4_t){0.0, ..};
            VAR {
              ??? cun = nc[un];
              ??? cvn = nc[vn];
              ??? cwn = nc[wn];
              ??? cxn = bary;
              ??? us = r4_Scale(4.00;
              with (cun),
              double vs = r4_Scale(4.00,cvn);
              double ws = r4_Scale(4.00,cwn);
              double xs = r4_Scale(4.00,cxn);
                b  == r4_Add(r4_Add(r4_Add(r4_Add(B, cun),cvn),cwn),cxn),
              double ba = r4_Scale(1.0/4.0, b);
              double ue = r4_Add(us, ba);
              double ve = r4_Add(vs, ba);
              double we = r4_Add(ws, ba);
              double xe = r4_Add(xs, ba);
              double ue3 = ProjectTo3D(ue);
              double ve3 = ProjectTo3D(ve);
              double we3 = ProjectTo3D(we);
              double xe3 = ProjectTo3D(xe);
                ra  == 0.008,
                cf  == (frgb_t){1.000, 0.745, 0.745},
                ce  == (frgb_t){0.0,0.0,0.0},
                tr  == 0.0
            ){
                /* Draw cylinders */
                WritePOVCylinder(inc,ue3,ve3,ra,ce,tr,TRUE);
                WritePOVCylinder(inc,ve3,we3,ra,ce,tr,TRUE);
                WritePOVCylinder(inc,ue3,we3,ra,ce,tr,TRUE);            
                WritePOVCylinder(inc,ue3,xe3,ra,ce,tr,TRUE);
                WritePOVCylinder(inc,ve3,xe3,ra,ce,tr,TRUE);
                WritePOVCylinder(inc,we3,xe3,ra,ce,tr,TRUE);

                /* Draw triangles  */
                WritePOVTriangle(inc,ue3,ve3,we3,cf,tr,TRUE);
                WritePOVTriangle(inc,ue3,ve3,xe3,cf,tr,TRUE);
                WritePOVTriangle(inc,ve3,we3,xe3,cf,tr,TRUE);
                WritePOVTriangle(inc,ue3,we3,xe3,cf,tr,TRUE);
              }
            } ExplodeTetrahedron;

        VAR
          Wa,Wb,Wc,Wd: r4_t;
          data4 : double;
          bary  : r4_t;
        {
          for (j = 0; j < top->wall.nelE; j++)
          {
            ??? a = top->wedge[j];
              Node_t ou = OrgV(a)->num;
            Node_t ov = OrgV(Clock(a))->num;
            Node_t nu = OrgV(Half(a))->num;
            Node_t nv = OrgV(Half(Clock(a)))->num;
            Node_t nx = OrgV(Clock(Half(a)))->num;
            Node_t ny = OrgV(Clock(Half(Clock(a))))->num;
            {
              assert(nx == ny);
              nc[nu] = c[ou];
              nc[nv] = c[ov];
              nc[ny] = r4_Scale(0.5, r4_Add(c[ou], c[ov]));
              nc[nx] = nc[ny];
            }
          }

          for (j = 0; j < top->wall.nel; j++)
          {
            ??? a = top->wall[j].pa;
              ??? vj = vf[j];
            ??? vjn = vj->num;
            {
              nc[vjn] = Triangulation.WallBarycenter(a,c);
            }          
          }
          
          /* compute the 4D viewing Matrix */
          CalcV4Matrix(Wa,Wb,Wc,Wd);

          /* compute the radius of the 4D configuration */
          data4 = Compute4DRadius();

          /* compute the barycenter of the 4D configuration */
          bary = Triangulation.Barycenter(ntop,nc,TRUE);
          
          for (j = 0; j < ntop.wall.nel; j++)
          {
            ??? a = ntop.wall[j].pa;
            Node_t un = OrgV(a)->num;
            Node_t vn = OrgV(NextE(a))->num;
            Node_t wn = OrgV(PrevE(a))->num;
            {
              if (o->element == Element.Cell)
              {
                ExplodeTetrahedron(un,vn,wn);
              }
              else
               if (o->element == Element.Wall){
                ExplodeWall(un,vn,wn);
              }                
            }          
          }
        }
      }        
    }
  }
} DoIt;

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);
    {                         
      ¦TRY

        argparser_get_keyword(pp, "-inFileTp");                               
        o->inFileTp = argparser_get_next(pp);
        argparser_get_keyword(pp, "-inFileSt");                               
        o->inFileSt = argparser_get_next(pp);                         
        argparser_get_keyword(pp, "-outFile");                               
        o->outFile = argparser_get_next(pp);
                argparser_get_keyword(pp, "-projection"); 
        o->projectionName =  argparser_get_next(pp);
        
        if ((! 0 == strcmp(o->projectionName, "Parallel") OR
                Text.Equal(o->projectionName, "Perspective"))
      )){
          argparser_error(pp, "Bad projection \"" & argparser_get_next(pp) & "\"\n")
        }

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

        if ((argparser_keyword_present(pp, "-explode")))
        {
          o->explode = TRUE;
          o->eleName = argparser_get_next(pp);
          if (0 == strcmp(o->eleName,"Wall")))
          {
            o->element = Element.Wall;
          }
          else
           if (0 == strcmp(o->eleName,"Cell"))){
            o->element = Element.Cell;
          }
          else
          {
            argparser_error(pp, "Bad element \"" & argparser_get_next(pp) & "\"\n")
          }
        };
                
        argparser_finish(pp);                                       
      
  -----------------------------------
#define _HELP \
  fprintf(stderr, "Usage: Subdivide \\\n" \
   "  -inFileTp <name> \\\n" \
   "  -inFileSt <name> \\\n" \
   "  -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" \
   "   ] [ -explode { Wall | Cell } ] \n");
          
END¦                                              
      }
    }
    return o;
  } /* END GetOptions */

{
  DoIt()
} Explode. 



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