home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / modelers / geomview / source.lha / Geomview / src / lib / gprim / discgrp / dgdirdom.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-08-28  |  10.4 KB  |  366 lines
  1. #include "geom.h"
  2. #include "polylistP.h"
  3. #include "discgrpP.h"
  4. #include "point.h"
  5. #include "winged_edge.h"
  6. #include "transform.h"
  7. #include "math.h"
  8.  
  9.     static WEpolyhedron    *wepoly1, **wepoly2; 
  10.         HPoint3 DGorigin = {0,0,0,1};
  11.         HPoint3 DGrandom = {.1,.2,.3,.4};
  12.  
  13. static int
  14. is_id(Transform t)
  15. {
  16.       register int i,j;
  17.  
  18.       for (i=0; i<4; ++i)
  19.         for (j=0; j<4; ++j)
  20.             if (fabs(t[i][j] - (i == j)) > .0005) return(0);
  21.       return(1);
  22. }
  23.  
  24. /* make sure that the center point attached to the discrete group
  25.    isn't fixed by some generator; if it is; redefine the center 
  26.    point to be the center of gravity of the orbit by the generators 
  27. */
  28. void
  29. DiscGrpCheckCPoint(register DiscGrp *dg)
  30. {
  31.     int i, cnt, fixed;
  32.     HPoint3 tmp, out;
  33.     float d;
  34.     
  35.     if (dg->gens == NULL) {
  36.     return;
  37.     }
  38.     
  39.     /* go through the generators, checking for fixed-pointed-ness */
  40.     for  (i = 0, fixed = 0 ; i < dg->gens->num_el; ++i )    
  41.     {
  42.     HPt3Transform(dg->gens->el_list[i].tform, &dg->cpoint, &tmp);
  43.     d = HPt3SpaceDistance(&dg->cpoint, &tmp, dg->attributes & DG_METRIC_BITS);
  44.     if (fabs(d) < .0005)  {
  45.         fixed = 1;
  46.         break;
  47.         }
  48.     }    
  49.  
  50.     /* no fixed points */
  51.     if (fixed == 0) return;
  52.  
  53.     /* clean out the special bit */
  54.     for  (i = 0 ; i < dg->gens->num_el; ++i )    
  55.     dg->gens->el_list[i].attributes &= ~DG_TMP;
  56.  
  57.     out.x = out.y = out.z = out.w = 0.0;
  58.     /* don't average but one of each generator, inverse pair */
  59.     for  (cnt = 0, i = 0 ; i < dg->gens->num_el; ++i )        {
  60.     if (!(dg->gens->el_list[i].attributes & DG_TMP))    {
  61.         HPt3Transform(dg->gens->el_list[i].tform, &DGrandom, &tmp);
  62.         HPt3Add(&tmp, &out, &out);
  63.         dg->gens->el_list[i].inverse->attributes |= DG_TMP;
  64.         cnt++;
  65.         }
  66.     }
  67.     HPt3Dehomogenize(&out, &out);
  68.     /* return it or set cpoint?? */
  69.     dg->cpoint = out;
  70. }
  71.  
  72. void
  73. DiscGrpAddInverses(register DiscGrp *discgrp)
  74. {
  75.     int i, j, found = 0;
  76.     Transform product, inverse;
  77.     DiscGrpElList *newgens;
  78.     
  79.     /* remove all identity matrices */
  80.     for (j=0, i=0; i<discgrp->gens->num_el; ++i)    {
  81.     if ( !is_id(discgrp->gens->el_list[i].tform) )    {
  82.         /* ought to have a DiscGrpElCopy() */
  83.         discgrp->gens->el_list[j] = discgrp->gens->el_list[i];
  84.         TmCopy(discgrp->gens->el_list[i].tform, 
  85.         discgrp->gens->el_list[j].tform);
  86.         j++;
  87.         }
  88.     }
  89.     /* store the new count */
  90.     discgrp->gens->num_el = j;
  91.  
  92.     for (i=0; i<discgrp->gens->num_el; ++i)    {
  93.       if (discgrp->gens->el_list[i].inverse == NULL)    {
  94.     /* look for inverse among the existing generators */
  95.     for (j=i; j<discgrp->gens->num_el; ++j)      {
  96.         TmConcat(discgrp->gens->el_list[i].tform, discgrp->gens->el_list[j].tform, product);
  97.         if ( is_id(product) )     {
  98.         discgrp->gens->el_list[i].inverse = &discgrp->gens->el_list[j];
  99.         discgrp->gens->el_list[j].inverse = &discgrp->gens->el_list[i];
  100.         found++;
  101.         }
  102.         }
  103.     }
  104.       else found++;
  105.     }
  106.  
  107.     newgens = OOGLNew(DiscGrpElList);
  108.     newgens->num_el = 2 * discgrp->gens->num_el - found;
  109.     newgens->el_list = OOGLNewN(DiscGrpEl, newgens->num_el );
  110.  
  111.     bcopy(discgrp->gens->el_list, newgens->el_list, sizeof(DiscGrpEl) * discgrp->gens->num_el);
  112.  
  113.     /* now go through looking for group elements without inverses */
  114.     {
  115.     char c;
  116.     j = discgrp->gens->num_el;
  117.     for (i=0; i<discgrp->gens->num_el; ++i)    {
  118.     if (newgens->el_list[i].inverse == NULL)    {
  119.         newgens->el_list[j+i] = newgens->el_list[i];
  120.         /* make the symbol of the inverse be the 'other case' */
  121.         c = newgens->el_list[i].word[0];
  122.         if (c < 'a') newgens->el_list[j+i].word[0] = c + 32;
  123.         else newgens->el_list[j+i].word[0] = c - 32;
  124.         TmInvert( newgens->el_list[i].tform, newgens->el_list[j+i].tform);
  125.         newgens->el_list[j+i].inverse = &newgens->el_list[i];
  126.         newgens->el_list[i].inverse = &newgens->el_list[j+i];
  127.         }
  128.     else j--;
  129.     }
  130.     }
  131.  
  132.     DiscGrpElListDelete(discgrp->gens);
  133.     discgrp->gens = newgens;
  134. }
  135.  
  136. void
  137. DiscGrpSetupDirdom(register DiscGrp *discgrp)
  138. {
  139.     WEpolyhedron *dd;
  140.  
  141.     if ( discgrp->nhbr_list )    {
  142.     OOGLFree(discgrp->nhbr_list->el_list);
  143.     OOGLFree(discgrp->nhbr_list);
  144.     }
  145.  
  146.     /* worry about fixed points */
  147.     DiscGrpCheckCPoint(discgrp);
  148.     dd = DiscGrpMakeDirdom(discgrp, &discgrp->cpoint, 0);
  149.     discgrp->nhbr_list = DiscGrpExtractNhbrs(dd);
  150. }
  151.  
  152. /* find the group element whose 'center point' is closest to the point poi */
  153. DiscGrpEl *
  154. DiscGrpClosestGroupEl(register DiscGrp *discgrp, HPoint3 *poi)
  155. {
  156.     int count, i, closeri;
  157.     int metric;
  158.     float min, d;
  159.     HPoint3 pt0, pt1;
  160.     DiscGrpEl *closer, *closest = OOGLNew(DiscGrpEl);
  161.     Transform cinv;
  162.  
  163.     TmIdentity(closest->tform);
  164.  
  165.     if (!discgrp->nhbr_list) DiscGrpSetupDirdom(discgrp);
  166.  
  167.     metric = discgrp->attributes & (DG_METRIC_BITS);
  168.  
  169.     /* iterate until we're in the fundamental domain */
  170.     count = 0;
  171.     closeri = -1;
  172.     pt0 = *poi;
  173.     while (count < 1000 && closeri != 0)    {
  174.       for (i = 0; i<discgrp->nhbr_list->num_el; ++i)    {
  175.     HPt3Transform(discgrp->nhbr_list->el_list[i].tform, &discgrp->cpoint, &pt1);    
  176.      d = HPt3SpaceDistance(&pt0, &pt1, metric);
  177.     if (i==0)     {
  178.         min = d;
  179.         closer = &discgrp->nhbr_list->el_list[i];
  180.         closeri = i;
  181.         }
  182.     else if (d < min)    {
  183.         min = d;
  184.         closer = &discgrp->nhbr_list->el_list[i];
  185.         closeri = i;
  186.         }
  187.     } 
  188.       count++;
  189.       if (closeri)    {
  190.           TmConcat(closer->tform, closest->tform, closest->tform);
  191.           /* move the point of interest by the inverse of the closest nhbr 
  192.           and iterate */
  193.           TmInvert(closest->tform, cinv);
  194.           HPt3Transform(cinv, poi, &pt0);
  195.       }
  196.     }
  197.     return (closest);
  198. }
  199.     static ColorA white = {1,1,1,1};
  200.  
  201. /* return a list of group elements corresponding to the faces of the
  202.    dirichlet domain */
  203. DiscGrpElList *
  204. DiscGrpExtractNhbrs( WEpolyhedron *wepoly )
  205. {
  206.     register int     i,j,k;
  207.     int         flag = 0, wl;
  208.     WEface         *fptr;    
  209.     DiscGrpElList    *mylist;
  210.     ColorA        GetCmapEntry();
  211.     
  212.     if (!wepoly)    return(NULL);
  213.  
  214.     /* should use realloc() here to take care of large groups...*/
  215.     mylist = OOGLNew(DiscGrpElList);
  216.     mylist->el_list = OOGLNewN(DiscGrpEl, wepoly->num_faces + 1);
  217.     mylist->num_el = wepoly->num_faces + 1;
  218.     
  219.     /* include the identity matrix */
  220.     TmIdentity( mylist->el_list[0].tform);
  221.     mylist->el_list[0].color = white;
  222.  
  223.     /* read the matrices corresponding to the faces of dirichlet domain */
  224.     for  (fptr = wepoly->face_list, k = 1; 
  225.     k<=wepoly->num_faces && fptr != NULL; 
  226.     k++, fptr = fptr->next)
  227.       {
  228.       for (i=0; i<4; ++i)
  229.     for (j=0; j<4; ++j)
  230.       /* the group elements stored in fptr are transposed! */
  231.       mylist->el_list[k].tform[j][i] = fptr->group_element[i][j];
  232.       mylist->el_list[k].color = GetCmapEntry(fptr->fill_tone);
  233.       }
  234.     if (mylist->num_el != k) 
  235.     OOGLError(1,"Incorrect number of nhbrs.\n");;
  236.  
  237.     return(mylist);
  238. }
  239.  
  240. /* attempt to create a scaled copy of fundamental domain for spherical
  241.    groups by taking weighted sum of vertices with the distinguished point */
  242. static void
  243. DiscGrpScalePolyList(DiscGrp *dg, PolyList *dirdom,  HPoint3 *pt0, float scale)
  244. {
  245.     PolyList *copy;
  246.     int i, metric;
  247.     HPoint3 tmp1, tmp2, tpt0, tpt1;
  248.     HPt3Copy(pt0, &tpt0);
  249.     metric = dg->attributes & DG_METRIC_BITS;
  250.     if (metric != DG_EUCLIDEAN)     {
  251.     HPt3SpaceNormalize(&tpt0, metric);
  252.         HPt3Scale( 1.0 - scale, &tpt0, &tmp2);
  253.         for (i=0; i<dirdom->n_verts; ++i)    {
  254.         HPt3Copy(&dirdom->vl[i].pt, &tpt1);
  255.             HPt3SpaceNormalize(&tpt1, metric);
  256.             HPt3SpaceNormalize(&tpt1, metric);
  257.         HPt3Scale( scale, &tpt1, &tmp1);
  258.         HPt3Add(&tmp1, &tmp2, &dirdom->vl[i].pt);
  259.         }
  260.     }
  261.     else    {
  262.     Transform T, TT, ITT, tmp;
  263.     static HPoint3 average = {0,0,0,0};
  264.     /* compute average */
  265.         for (i=0; i<dirdom->n_verts; ++i)    
  266.         HPt3Add(&average, &dirdom->vl[i].pt, &average);
  267.     HPt3Dehomogenize(&average, &average);
  268.  
  269.     TmTranslate(TT, average.x, average.y, average.z );
  270.     TmInvert(TT, ITT);
  271.     TmScale(T, scale, scale, scale);
  272.     TmConcat(ITT, T, tmp);
  273.     TmConcat(tmp, TT, tmp);
  274.         for (i=0; i<dirdom->n_verts; ++i)    
  275.         HPt3Transform(tmp, &dirdom->vl[i].pt, &dirdom->vl[i].pt);    
  276.     }
  277. }
  278.  
  279.     Geom *small_dd, *large_dd;
  280. Geom *
  281. DiscGrpDirDom(register DiscGrp *dg)
  282. {
  283.     Geom *oogldirdom;
  284.     WEpolyhedron *dd;
  285.     HPoint3 pt1;
  286.     extern Geom             *WEPolyhedronToPolyList();
  287.     Geom *mylist, *smlist;
  288.  
  289.       if (dg->flag & DG_DDBEAM)    {
  290.         WEpolyhedron *poly = DiscGrpMakeDirdom( dg, &dg->cpoint, 0); 
  291.         Geom *beams;
  292.         beams = WEPolyhedronToBeams(poly, dg->scale);
  293.         return(beams);
  294.     }
  295.       else    {
  296.     float scale;
  297.     /* first a full-size, wireframe version of dd */
  298.         dd = DiscGrpMakeDirdom(dg, &dg->cpoint, 0);
  299.         if (dd) {
  300.         oogldirdom = WEPolyhedronToPolyList(dd);
  301.          scale = 1.0;
  302.         DiscGrpScalePolyList(dg, (PolyList*)oogldirdom, &dg->cpoint, scale);
  303.         large_dd = oogldirdom;
  304.         large_dd->ap = ApCreate(AP_DO, APF_EDGEDRAW, AP_DONT, APF_FACEDRAW, AP_END);
  305.         }
  306.         else return((Geom *) NULL);
  307.     /* next a scaled version with cusps cut off */
  308.         dd = DiscGrpMakeDirdom(dg, &dg->cpoint, 1);
  309.         if (dd) {
  310.         oogldirdom = WEPolyhedronToPolyList(dd);
  311.         DiscGrpScalePolyList(dg, (PolyList*)oogldirdom, &dg->cpoint, dg->scale);
  312.         small_dd = oogldirdom;
  313.         small_dd->ap = ApCreate(AP_DONT, APF_EDGEDRAW, AP_DO, APF_FACEDRAW, AP_END);
  314.         }
  315.         else return(NULL);
  316.     smlist = GeomCreate("list", CR_GEOM, small_dd, CR_END);
  317.     mylist = GeomCreate("list", CR_GEOM, large_dd, CR_CDR, smlist, CR_END);
  318.     return(mylist);
  319.     }
  320. }
  321.  
  322. static WEpolyhedron *
  323. DiscGrpMakeDirdom(DiscGrp *gamma, HPoint3 *poi, int slice)
  324. {
  325.     int i, j, k;
  326.     WEvertex    *vptr;
  327.     Transform tr;
  328.     proj_matrix *gen_list;
  329.         point origin;
  330.     int metric, transp;
  331.  
  332.     transp = gamma->attributes & DG_TRANSPOSED;
  333.     /* transform from floating point to double, essentially */
  334.         gen_list = OOGLNewN(proj_matrix, gamma->gens->num_el);
  335.     /* jeff week's code basically uses transposed matrices, so
  336.     if the transposed flag is set, we do nothing, otherwise we
  337.     have to transpose! */
  338.         for (i=0; i<gamma->gens->num_el; ++i) 
  339.         for (j=0; j<4; ++j)
  340.             for (k=0; k<4; ++k)
  341.             {
  342.             if (transp) gen_list[i][j][k] = gamma->gens->el_list[i].tform[j][k];
  343.             else  gen_list[i][k][j] = gamma->gens->el_list[i].tform[j][k];
  344.             }
  345.         origin[0] = poi->x;
  346.         origin[1] = poi->y;
  347.         origin[2] = poi->z;
  348.         origin[3] = poi->w;
  349.  
  350.     wepoly2 = &wepoly1;
  351.     metric = (gamma->attributes & DG_METRIC_BITS); 
  352.     do_weeks_code(wepoly2, origin, gen_list, gamma->gens->num_el, metric,slice);
  353.  
  354.     free(gen_list);
  355.  
  356.     /* turn off the 'new dirdom' bit */
  357.         gamma->flag &= ~DG_NEWDIRDOM;
  358.     return(*wepoly2);
  359. }
  360.  
  361. Geom *
  362. DiscGrpBeamedDirDom(DiscGrp *dg, float alpha)
  363. {
  364. }
  365.  
  366.