home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / renaisnc / lib.lha / Lib / meshsub.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-01-04  |  11.4 KB  |  445 lines
  1. /* File: meshsub.c
  2. ** Author: Charles Loop
  3. ** Purpose: Contains routines needed for the creation of a subdivided mesh.
  4. ** Last Modified: 2 July, 1987
  5. ** Reason : Creation.
  6. ** 
  7. **  new version
  8. */
  9.  
  10. #include <stdio.h>
  11. #include "Math.h"
  12. #include "4D.h"
  13. #include "mesh.h"
  14.  
  15.  
  16. static double a[maxN+1],a_[maxN+1],
  17.               b[maxN+1],b_[maxN+1];
  18.  
  19. static int FlatShading;
  20. static int Converge;
  21.  
  22. /*****************************************************************************/
  23.  
  24. set_a(l)
  25.      int l;
  26. {
  27.   double r;
  28.   int i;
  29.  
  30.   for  (i = 3; i <= maxN; i++){
  31.     if (l) {
  32.       r = 0.375 + 0.25*cos(TwoPi/i);
  33.       a[i] = r*r + 0.375;
  34.       a_[i] = (1.0 - a[i])/i;
  35.       b[i] = 3.0/(8.0*(1.0 - r*r));
  36.       b_[i] = (1.0 - b[i])/i;
  37.     }
  38.     else {
  39.       b[i] = 1.0;
  40.       b_[i] = 0.0;
  41.     }
  42.   }
  43. }
  44.  
  45. /*****************************************************************************/
  46.  
  47. unsigned OrdAdjList(P,e)
  48.      PointType4D P[];
  49.      Edge *e;
  50. {
  51.   Edge *ep;
  52.   unsigned N = 1;
  53.  
  54.   if (e) {
  55.     if (e->e) {
  56.       P[0] = *e->next->p;
  57.       for (ep = e->e; ep != e->next && ep != (Edge *) NULL; ep = ep->next->next->e)
  58.     P[N++] = *ep->p;
  59.  
  60.       if (ep)
  61.     return N;
  62.       else
  63.     return 0;
  64.     }
  65.     return 0;
  66.   }
  67.   return 0;
  68. }
  69.  
  70. /*****************************************************************************/
  71.  
  72. Subdivide(Pr,Ps,Pt,pr,ps,pt,qr,qs,qt,Nr,Ns,Nt)
  73. PointType4D Pr[],Ps[],Pt[],pr[],ps[],pt[],qr[],qs[],qt[];
  74. unsigned Nr,Ns,Nt;
  75. {
  76.   PointType4D Vr,Vs,Vt,vr,vs,vt;
  77.   unsigned i;
  78.  
  79.   Vr = Pt[0]; Vs = Pr[0]; Vt = Ps[0];
  80.  
  81.   vr.P[X] = 0.0;  vr.P[Y] = 0.0;  vr.P[Z] = 0.0;
  82.   vs.P[X] = 0.0;  vs.P[Y] = 0.0;  vs.P[Z] = 0.0;
  83.   vt.P[X] = 0.0;  vt.P[Y] = 0.0;  vt.P[Z] = 0.0;
  84.  
  85.   for (i = 0; i < Nr; i++){
  86.     pr[i].P[X] = 0.125*(Pr[(i+Nr-1)%Nr].P[X]+Pr[(i+1)%Nr].P[X])+0.375*(Pr[i].P[X] + Vr.P[X]);
  87.     pr[i].P[Y] = 0.125*(Pr[(i+Nr-1)%Nr].P[Y]+Pr[(i+1)%Nr].P[Y])+0.375*(Pr[i].P[Y] + Vr.P[Y]);
  88.     pr[i].P[Z] = 0.125*(Pr[(i+Nr-1)%Nr].P[Z]+Pr[(i+1)%Nr].P[Z])+0.375*(Pr[i].P[Z] + Vr.P[Z]);
  89.     vr.P[X] = vr.P[X] + Pr[i].P[X];
  90.     vr.P[Y] = vr.P[Y] + Pr[i].P[Y];
  91.     vr.P[Z] = vr.P[Z] + Pr[i].P[Z];
  92.   }
  93.   vr.P[X] = a[Nr]*Vr.P[X] + a_[Nr]*vr.P[X];
  94.   vr.P[Y] = a[Nr]*Vr.P[Y] + a_[Nr]*vr.P[Y];
  95.   vr.P[Z] = a[Nr]*Vr.P[Z] + a_[Nr]*vr.P[Z];
  96.  
  97.   for (i = 0; i < Ns; i++){
  98.     ps[i].P[X] = 0.125*(Ps[(i+Ns-1)%Ns].P[X]+Ps[(i+1)%Ns].P[X])+0.375*(Ps[i].P[X] + Vs.P[X]);
  99.     ps[i].P[Y] = 0.125*(Ps[(i+Ns-1)%Ns].P[Y]+Ps[(i+1)%Ns].P[Y])+0.375*(Ps[i].P[Y] + Vs.P[Y]);
  100.     ps[i].P[Z] = 0.125*(Ps[(i+Ns-1)%Ns].P[Z]+Ps[(i+1)%Ns].P[Z])+0.375*(Ps[i].P[Z] + Vs.P[Z]);
  101.     vs.P[X] = vs.P[X] + Ps[i].P[X];
  102.     vs.P[Y] = vs.P[Y] + Ps[i].P[Y];
  103.     vs.P[Z] = vs.P[Z] + Ps[i].P[Z];
  104.   }
  105.   vs.P[X] = a[Ns]*Vs.P[X] + a_[Ns]*vs.P[X];
  106.   vs.P[Y] = a[Ns]*Vs.P[Y] + a_[Ns]*vs.P[Y];
  107.   vs.P[Z] = a[Ns]*Vs.P[Z] + a_[Ns]*vs.P[Z];
  108.  
  109.   for (i = 0; i < Nt; i++){
  110.     pt[i].P[X] = 0.125*(Pt[(i+Nt-1)%Nt].P[X]+Pt[(i+1)%Nt].P[X])+0.375*(Pt[i].P[X] + Vt.P[X]);
  111.     pt[i].P[Y] = 0.125*(Pt[(i+Nt-1)%Nt].P[Y]+Pt[(i+1)%Nt].P[Y])+0.375*(Pt[i].P[Y] + Vt.P[Y]);
  112.     pt[i].P[Z] = 0.125*(Pt[(i+Nt-1)%Nt].P[Z]+Pt[(i+1)%Nt].P[Z])+0.375*(Pt[i].P[Z] + Vt.P[Z]);
  113.     vt.P[X] = vt.P[X] + Pt[i].P[X];
  114.     vt.P[Y] = vt.P[Y] + Pt[i].P[Y];
  115.     vt.P[Z] = vt.P[Z] + Pt[i].P[Z];
  116.   }
  117.   vt.P[X] = a[Nt]*Vt.P[X] + a_[Nt]*vt.P[X];
  118.   vt.P[Y] = a[Nt]*Vt.P[Y] + a_[Nt]*vt.P[Y];
  119.   vt.P[Z] = a[Nt]*Vt.P[Z] + a_[Nt]*vt.P[Z];
  120.  
  121.   qr[0] = vt;        qs[0] = vr;        qt[0] = vs;
  122.   qr[1] = pt[0];     qs[1] = pr[0];     qt[1] = ps[0];
  123.   qr[2] = pr[0];     qs[2] = ps[0];     qt[2] = pt[0];
  124.   qr[3] = vs;        qs[3] = vt;        qt[3] = vr;
  125.   qr[4] = ps[Ns-1];  qs[4] = pt[Nt-1];  qt[4] = pr[Nr-1];
  126.   qr[5] = pt[2];     qs[5] = pr[2];     qt[5] = ps[2];
  127.   qr[6] = qr[0];     qs[6] = qs[0];     qt[6] = qt[0];
  128.   qr[7] = qr[1];     qs[7] = qs[1];     qt[7] = qt[1];
  129.  
  130. }
  131.  
  132. /*****************************************************************************/
  133.  
  134. Plane *PlaneEquation(Vr,Vs,Vt)
  135.      Vertex *Vr,*Vs,*Vt;
  136. {
  137.   Plane *P;
  138.  
  139.   P = (Plane *)malloc(sizeof(Plane));
  140.  
  141.   P->a = (Vs->y - Vr->y)*(Vt->z - Vr->z) - (Vt->y - Vr->y)*(Vs->z - Vr->z);
  142.   P->b = (Vt->x - Vr->x)*(Vs->z - Vr->z) - (Vs->x - Vr->x)*(Vt->z - Vr->z);
  143.   P->c = (Vs->x - Vr->x)*(Vt->y - Vr->y) - (Vt->x - Vr->x)*(Vs->y - Vr->y);
  144.   P->d = -(P->a*Vr->x + P->b*Vr->y + P->c*Vr->z);
  145.  
  146.   if (fabs(P->a) > fabs(P->b))
  147.     if (fabs(P->a) > fabs(P->c))
  148.       P->maxcomponent = 0;
  149.     else
  150.       P->maxcomponent = 2;
  151.   else if (fabs(P->b) > fabs(P->c))
  152.     P->maxcomponent = 1;
  153.   else
  154.     P->maxcomponent = 2;
  155.  
  156.   return(P);
  157.  
  158. }
  159.  
  160. /****************************************************************************/
  161.  
  162. Vertex *SetVertex(V,P,N)
  163.      PointType4D *V,*P;
  164.      unsigned N;
  165. {
  166.   Vertex *Q = (Vertex *)malloc(sizeof(Vertex));
  167.   PointType4D Fr0, Fr1, Fi1, Fr2, Fi2;
  168.   double C1, S1, C2, S2, Theta, TwoTheta, tmp;
  169.   unsigned i;
  170.  
  171.   Fr0.P[X] = 0.0;  Fr0.P[Y] = 0.0;  Fr0.P[Z] = 0.0;
  172.   Fr1.P[X] = 0.0;  Fr1.P[Y] = 0.0;  Fr1.P[Z] = 0.0;
  173.   Fi1.P[X] = 0.0;  Fi1.P[Y] = 0.0;  Fi1.P[Z] = 0.0;
  174.   Fr2.P[X] = 0.0;  Fr2.P[Y] = 0.0;  Fr2.P[Z] = 0.0;
  175.   Fi2.P[X] = 0.0;  Fi2.P[Y] = 0.0;  Fi2.P[Z] = 0.0;
  176.  
  177.   Theta = TwoPi/N;
  178.   TwoTheta = 2.0*Theta;
  179.  
  180.   for (i = 0; i < N; i++) {
  181.     C1 = cos(Theta*i);  S1 = sin(Theta*i);
  182.     C2 = cos(TwoTheta*i);  S2 = sin(TwoTheta*i);
  183.  
  184.     Fr0.P[X] += P[i].P[X];
  185.     Fr0.P[Y] += P[i].P[Y];
  186.     Fr0.P[Z] += P[i].P[Z];
  187.  
  188.     Fr1.P[X] += C1*P[i].P[X];
  189.     Fr1.P[Y] += C1*P[i].P[Y];
  190.     Fr1.P[Z] += C1*P[i].P[Z];
  191.     Fi1.P[X] += S1*P[i].P[X];
  192.     Fi1.P[Y] += S1*P[i].P[Y];
  193.     Fi1.P[Z] += S1*P[i].P[Z];
  194.  
  195.     Fr2.P[X] += C2*P[i].P[X];
  196.     Fr2.P[Y] += C2*P[i].P[Y];
  197.     Fr2.P[Z] += C2*P[i].P[Z];
  198.     Fi2.P[X] += S2*P[i].P[X];
  199.     Fi2.P[Y] += S2*P[i].P[Y];
  200.     Fi2.P[Z] += S2*P[i].P[Z];
  201.   }
  202.  
  203.   if (Converge) {
  204.     Q->x = b[N]*V->P[X] + b_[N]*Fr0.P[X];
  205.     Q->y = b[N]*V->P[Y] + b_[N]*Fr0.P[Y];
  206.     Q->z = b[N]*V->P[Z] + b_[N]*Fr0.P[Z];
  207.   } else {
  208.     Q->x = V->P[X];
  209.     Q->y = V->P[Y];
  210.     Q->z = V->P[Z];
  211.   }
  212.  
  213.   Q->a = Fr1.P[Y]*Fi1.P[Z] - Fi1.P[Y]*Fr1.P[Z];
  214.   Q->b = Fr1.P[Z]*Fi1.P[X] - Fi1.P[Z]*Fr1.P[X];
  215.   Q->c = Fr1.P[X]*Fi1.P[Y] - Fi1.P[X]*Fr1.P[Y];
  216.  
  217.   tmp = sqrt(Q->a*Q->a + Q->b*Q->b + Q->c*Q->c);
  218.  
  219.   Q->a = Q->a/tmp;
  220.   Q->b = Q->b/tmp;
  221.   Q->c = Q->c/tmp;
  222.  
  223. /* compute coefficients of first fundamental form */
  224. /*
  225.   e = Fr1.P[X]*Fr1.P[X] + Fr1.P[Y]*Fr1.P[Y] + Fr1.P[Z]*Fr1.P[Z];
  226.   f = Fr1.P[X]*Fi1.P[X] + Fr1.P[Y]*Fi1.P[Y] + Fr1.P[Z]*Fi1.P[Z];
  227.   g = Fi1.P[X]*Fi1.P[X] + Fi1.P[Y]*Fi1.P[Y] + Fi1.P[Z]*Fi1.P[Z];
  228. */
  229. /* compute coefficients of second fundamental form */
  230. /*
  231.   l = Q->a*(0.5*(Fr0.P[X] - N*V->P[X]) + Fr2.P[X])
  232.     + Q->b*(0.5*(Fr0.P[Y] - N*V->P[Y]) + Fr2.P[Y])
  233.     + Q->c*(0.5*(Fr0.P[Z] - N*V->P[Z]) + Fr2.P[Z]);
  234.  
  235.   m = Q->a*Fi2.P[X] + Q->b*Fi2.P[Y] + Q->c*Fi2.P[Z];
  236.  
  237.   n = Q->a*(0.5*(Fr0.P[X] - N*V->P[X]) - Fr2.P[X])
  238.     + Q->b*(0.5*(Fr0.P[Y] - N*V->P[Y]) - Fr2.P[Y])
  239.     + Q->c*(0.5*(Fr0.P[Z] - N*V->P[Z]) - Fr2.P[Z]);
  240. */
  241. /* compute Gaussian curvature */
  242. /*
  243.   Q->k = N*N*(l*n - m*m)/(e*g - f*f);
  244. */
  245.   return(Q);
  246.  
  247. }
  248.  
  249. /****************************************************************************/
  250.  
  251. minmax(r,s,t,m,M)
  252. double r,s,t,*m,*M;
  253. {
  254.   if (r > s)
  255.     if (s > t)
  256.       {*M = r; *m = t;}
  257.     else if (r > t)
  258.       {*M = r; *m = s;}
  259.     else
  260.       {*M = t; *m = s;}
  261.   else if (s > t)
  262.     if (r > t)
  263.       {*M = s; *m = t;}
  264.     else
  265.       {*M = s; *m = r;}
  266.   else
  267.     {*M = t; *m = r;}
  268. }
  269.  
  270. /****************************************************************************/
  271.  
  272. Min(a,b,c,d,m)
  273. double a,b,c,d,*m;
  274. {
  275.   *m = (a < b) ? a : b;
  276.   *m = (c < *m) ? c : *m;
  277.   *m = (d < *m) ? d : *m;
  278. }
  279.  
  280. /****************************************************************************/
  281.  
  282. Max(a,b,c,d,M)
  283. double a,b,c,d,*M;
  284. {
  285.   *M = (a > b) ? a : b;
  286.   *M = (c > *M) ? c : *M;
  287.   *M = (d > *M) ? d : *M;
  288. }
  289.  
  290. /****************************************************************************/
  291.  
  292. CombindBoxes(box,s0,s1,s2,s3)
  293. double *box;
  294. Mesh *s0,*s1,*s2,*s3;
  295. {
  296.  
  297.     Min(s0->box[X][MIN],s1->box[X][MIN],s2->box[X][MIN],s3->box[X][MIN],box);
  298.     Min(s0->box[Y][MIN],s1->box[Y][MIN],s2->box[Y][MIN],s3->box[Y][MIN],box+2);
  299.     Min(s0->box[Z][MIN],s1->box[Z][MIN],s2->box[Z][MIN],s3->box[Z][MIN],box+4);
  300.  
  301.     Max(s0->box[X][MAX],s1->box[X][MAX],s2->box[X][MAX],s3->box[X][MAX],box+1);
  302.     Max(s0->box[Y][MAX],s1->box[Y][MAX],s2->box[Y][MAX],s3->box[Y][MAX],box+3);
  303.     Max(s0->box[Z][MAX],s1->box[Z][MAX],s2->box[Z][MAX],s3->box[Z][MAX],box+5);
  304.  
  305. }
  306.  
  307. /****************************************************************************/
  308.  
  309. Generate(node, Pr, Ps, Pt, Nr, Ns, Nt, l)
  310. Mesh *node;
  311. PointType4D Pr[], Ps[], Pt[];
  312. int Nr, Ns, Nt, l;
  313. {
  314.   Mesh *submesh;
  315.   PointType4D pr[maxN], ps[maxN], pt[maxN], qr[8], qs[8], qt[8];
  316.   int i;
  317.  
  318.   if (l) {
  319.     node->type = SUBMESH;
  320.  
  321.     submesh = (Mesh *)malloc(4*sizeof(Mesh));
  322.  
  323.     Subdivide(Pr, Ps, Pt, pr, ps, pt, qr, qs, qt, Nr, Ns, Nt);
  324.  
  325.     Generate(submesh  , qr+1, qs+1, qt+1,  6,  6,  6, l-1);
  326.     Generate(submesh+1,   pr, qt+2,   qs, Nr,  6,  6, l-1);
  327.     Generate(submesh+2,   qt,   ps, qr+2,  6, Ns,  6, l-1);
  328.     Generate(submesh+3, qs+2,   qr,   pt,  6,  6, Nt, l-1);
  329.  
  330.     CombindBoxes(node->box, submesh, submesh+1, submesh+2, submesh+3);
  331.  
  332.     for (i = 0; i < 3; i++)
  333.       submesh[i].next = submesh+i+1;
  334.  
  335.     node->sub.m = submesh;
  336.  
  337.   }
  338.   else {
  339.     node->type = TRIANGLE;
  340.  
  341.     node->sub.t = (Triangle *)malloc(sizeof(Triangle));
  342.  
  343.     node->sub.t->V[0] = SetVertex(Pt, Pr, Nr);
  344.     node->sub.t->V[1] = SetVertex(Pr, Ps, Ns);
  345.     node->sub.t->V[2] = SetVertex(Ps, Pt, Nt);
  346.  
  347.     minmax(node->sub.t->V[0]->x,node->sub.t->V[1]->x,node->sub.t->V[2]->x,
  348.        &node->box[X][MIN],&node->box[X][MAX]);
  349.     minmax(node->sub.t->V[0]->y,node->sub.t->V[1]->y,node->sub.t->V[2]->y,
  350.        &node->box[Y][MIN],&node->box[Y][MAX]);
  351.     minmax(node->sub.t->V[0]->z,node->sub.t->V[1]->z,node->sub.t->V[2]->z,
  352.        &node->box[Z][MIN],&node->box[Z][MAX]);
  353.  
  354.     node->sub.t->pl = PlaneEquation(node->sub.t->V[0],
  355.                   node->sub.t->V[1],
  356.                   node->sub.t->V[2]);
  357.     if (FlatShading) {
  358.       node->sub.t->V[0]->a =
  359.       node->sub.t->V[1]->a =
  360.       node->sub.t->V[2]->a = node->sub.t->pl->a;
  361.       node->sub.t->V[0]->b =
  362.       node->sub.t->V[1]->b =
  363.       node->sub.t->V[2]->b = node->sub.t->pl->b;
  364.       node->sub.t->V[0]->c =
  365.       node->sub.t->V[1]->c =
  366.       node->sub.t->V[2]->c = node->sub.t->pl->c;
  367.     }
  368.   }
  369. }
  370.  
  371. /****************************************************************************/
  372.  
  373.  
  374. Traverse(m, level)
  375.      Mesh *m;
  376.      int level;
  377. {
  378.   Mesh *mp;
  379.   PointType4D Pr[maxN], Ps[maxN], Pt[maxN];
  380.   int Nr, Ns, Nt, i;
  381.  
  382.   if (m) {
  383.     if (m->type == MESH) {
  384.       if (m->sub.m) {
  385.     for (mp = m->sub.m; mp != (Mesh *) 0; mp = mp->next)
  386.       Traverse(mp, level);
  387.  
  388.     for (i = 0; i < 3; i++) {
  389.       m->box[i][MIN] = m->sub.m->box[i][MIN];
  390.       m->box[i][MAX] = m->sub.m->box[i][MAX];
  391.     }
  392.  
  393.     for (mp = m->sub.m->next; mp != (Mesh *) 0; mp = mp->next)
  394.       for (i = 0; i < 3; i++) {
  395.         if (mp->box[i][MIN] < m->box[i][MIN])
  396.           m->box[i][MIN] = mp->box[i][MIN];
  397.         if (mp->box[i][MAX] > m->box[i][MAX])
  398.           m->box[i][MAX] = mp->box[i][MAX];
  399.       }
  400.       }
  401.     }
  402.     else if (m->type == FACE) {
  403.       if (m->sub.e) {
  404.  
  405.     Nr = OrdAdjList(Pr, m->sub.e);
  406.     Ns = OrdAdjList(Ps, m->sub.e->next);
  407.     Nt = OrdAdjList(Pt, m->sub.e->next->next);
  408.  
  409.     if (Nr && Ns && Nt)
  410.       Generate(m, Pr, Ps, Pt, Nr, Ns, Nt, level);
  411.     else
  412.       {
  413.         minmax(Pr[0].P[X], Ps[0].P[X], Pt[0].P[X],
  414.            &m->box[X][MIN],&m->box[X][MAX]);
  415.         minmax(Pr[0].P[Y], Ps[0].P[Y], Pt[0].P[Y],
  416.            &m->box[Y][MIN],&m->box[Y][MAX]);
  417.         minmax(Pr[0].P[Z], Ps[0].P[Z], Pt[0].P[Z],
  418.            &m->box[Z][MIN],&m->box[Z][MAX]);
  419.       }
  420.       }
  421.     }
  422.   }
  423. }
  424.  
  425. MeshSub(m, phong, converge, level)
  426.      Mesh *m;
  427.      int phong, converge, level;
  428. {
  429.  
  430.   if (phong)
  431.     FlatShading = 0;
  432.   else
  433.     FlatShading = 1;
  434.  
  435.   Converge = converge;
  436.  
  437.   set_a(level);
  438.   
  439.  
  440.   MeshPartition(m);
  441.   Traverse(m, level);
  442.  
  443. }
  444.  
  445.