home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / modelers / geomview / source.lha / Geomview / src / lib / geometry / cmodel / cm_geometry.c next >
Encoding:
C/C++ Source or Header  |  1993-06-09  |  3.9 KB  |  159 lines
  1. /* do the actual geometric computation parts of the "cmodel" program */
  2.  
  3. #include <math.h>
  4. #include "cmodelP.h"
  5.  
  6. void projective_to_conformal(int curv, HPoint3 *proj, Transform T, Point3 *conf)
  7. {
  8.    double norm, scale;
  9.    HPoint3 pt;
  10.  
  11.    HPt3Transform(T, proj, &pt);
  12.    norm = pt.x*pt.x + pt.y*pt.y + pt.z*pt.z;
  13.    if (curv)
  14.    {
  15.        norm = curv*norm + pt.w*pt.w;
  16.        norm = (norm < 0)? 0 : sqrt(norm);
  17.        scale = pt.w-curv*norm;
  18.    }
  19.    else
  20.        scale = -norm/pt.w;
  21.    Pt3Mul(1 / scale, (Point3 *)&pt, conf);
  22.    return;
  23.    }
  24.  
  25. void TgtTransform(Transform T, HPoint3 *p, Point3 *v, HPoint3 *tp, Point3 *tv)
  26. {
  27.    HPoint3 hv, thp, thv;
  28.    
  29.    Pt3Copy(v, (Point3 *)&hv);
  30.    hv.w = 0;
  31.    
  32.    HPt3Transform(T, p, tp);
  33.    HPt3Transform(T, &hv, &thv);
  34.    
  35.    Pt3Comb(1/tp->w, (Point3 *)&thv, -thv.w/tp->w/tp->w, (Point3 *)tp, tv);
  36.    
  37.    return;
  38.    }
  39.    
  40.  
  41. /* map a tangent vector from Projective model to Conformal model */
  42. void projective_vector_to_conformal(int curv, HPoint3 *pt,  Point3 *v,
  43.                     Transform T, Point3 *ppt, Point3 *pv)
  44. {
  45.    HPoint3 tp;
  46.    Point3 tv;
  47.    double norm,scale;
  48.  
  49.    /* transform point */
  50.    TgtTransform(T, pt, v, &tp, &tv);
  51.  
  52.    norm = tp.x*tp.x + tp.y*tp.y + tp.z*tp.z;
  53.    if (curv)
  54.    {
  55.        norm = curv*norm + tp.w*tp.w;
  56.        norm = (norm < 0)? 0 : sqrt(norm);
  57.        scale = tp.w-curv*norm;
  58.    }
  59.    else
  60.        scale = -norm/tp.w;
  61.  
  62.    Pt3Mul(1 / scale, (Point3 *)&tp, ppt);
  63.    if (curv)
  64.        Pt3Comb(norm/scale, &tv, Pt3Dot(ppt, &tv), ppt, pv);
  65.    else
  66.        Pt3Comb(tp.w/scale, &tv, 2*Pt3Dot(ppt, &tv), ppt, pv);
  67.  
  68.    Pt3Unit(pv);
  69.    return; 
  70.    }
  71.  
  72. /* given three vertices of a triangle in the conformal model this 
  73.    computes the center of the sphere on which the triangle lies. */
  74.  
  75. void triangle_polar_point(int curv,
  76.               const Point3 *a, const Point3 *b, const Point3 *c, 
  77.               HPoint3 *p)
  78. {
  79.    Point3 ab, bc, ca;
  80.  
  81.    Pt3Cross(a, b, &ab);
  82.    Pt3Cross(b, c, &bc);
  83.    Pt3Cross(c, a, &ca);
  84.    p->w = 2*Pt3Dot(a, &bc);
  85.    
  86.    Pt3Mul(Pt3Dot(a, a) - curv, &bc, &bc);
  87.    Pt3Mul(Pt3Dot(b, b) - curv, &ca, &ca);
  88.    Pt3Mul(Pt3Dot(c, c) - curv, &ab, &ab);
  89.  
  90.    Pt3Add(&ab, &bc, (Point3 *)p);
  91.    Pt3Add(&ca, (Point3 *)p, (Point3 *)p);
  92.  /*fprintf(stderr,"In triangle_polar got %f %f %f %f\n",p->x,p->y,p->z,p->w);*/
  93.  
  94.    return;
  95.    }
  96.  
  97. /* given two points on a geodesic in the conformal model this 
  98.    computes the centre of the euclidean circle formed by the geodesic */
  99.  
  100. void edge_polar_point(int curv, const Point3 *a, const Point3 *b, HPoint3 *p)
  101. {
  102.    double aa, ab, bb, ca, cb;
  103.    
  104.    aa = Pt3Dot(a, a);
  105.    ab = Pt3Dot(a, b);
  106.    bb = Pt3Dot(b, b);
  107.  
  108.    ca = (aa-ab)*bb + curv*(ab - bb);
  109.    cb = (bb-ab)*aa + curv*(ab - aa);
  110.    
  111.    Pt3Comb(ca, a, cb, b, (Point3 *)p);
  112.    p->w = 2 * (aa * bb - ab * ab);
  113.  
  114.    return;
  115.    }
  116.  
  117. struct vertex *edge_split(struct edge *e, double cosmaxbend)
  118. {
  119.    double cosbend, w;
  120.    Point3 m, mp, x, y, *a, *b, p;
  121.    double aa,ab,bb,ma,mb;
  122.  
  123.    a = (Point3 *)&e->v1->V.pt;
  124.    b = (Point3 *)&e->v2->V.pt;
  125.    w = e->polar.w;
  126.  
  127.  /*fprintf(stderr,"In edge_split\n");*/
  128.    if (w < .001) return NULL;
  129.    
  130.    Pt3Mul(1./w, (Point3 *)&e->polar, &p);
  131.    
  132.    Pt3Sub(a, &p, &x);
  133.    Pt3Sub(b, &p, &y);
  134.  
  135.    cosbend = Pt3Dot(&x,&y)/sqrt(Pt3Dot(&x,&x) * Pt3Dot(&y,&y));
  136.    if (cosmaxbend < cosbend)
  137.       return NULL;
  138.  
  139.  /*fprintf(stderr,"end pts %f %f %f, %f %f %f, polar %f %f %f %f\n",
  140.              a->x,a->y,a->z, b->x,b->y,b->z,
  141.              e->polar.x,e->polar.y,e->polar.z,e->polar.w);*/
  142.    Pt3Add(&x, &y, &m);
  143.    Pt3Mul(sqrt(Pt3Dot(&x, &x) / Pt3Dot(&m, &m)), 
  144.       &m, &m);
  145.    Pt3Add(&p, &m, &mp);
  146.    aa = Pt3Dot(a,a); ab = Pt3Dot(a,b); bb = Pt3Dot(b,b);
  147.    ma = Pt3Dot(a,&mp); mb = Pt3Dot(b,&mp);
  148.    if (aa*mb < ab*ma  ||  bb*ma < ab*mb)
  149.       Pt3Sub(&p,&m,&mp);
  150.  
  151.    /* to check we're doing the right thing we could do:
  152.    ma = Pt3Dot(a,&mp); mb = Pt3Dot(b,&mp);
  153.    if (aa*mb < ab*ma  ||  bb*ma < ab*mb)
  154.       fprintf(stderr,"Can't find right subdivision\n");
  155.    */
  156.  
  157.    return new_vertex(&mp, e->v1, e->v2);
  158.    }
  159.