home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch4_5 / ellipsoid.c < prev    next >
C/C++ Source or Header  |  1995-04-04  |  10KB  |  291 lines
  1. /* ellipsoid.c */
  2.  
  3. #include <stdio.h>
  4. #include <math.h>
  5. #include <malloc.h>
  6. #include "ellipsoid.h"
  7.  
  8. #ifndef M_PI_2
  9. #define    M_PI_2    1.57079632679489661923
  10. #endif
  11.  
  12. typedef struct slot { float cos, sin; enum { None, Only, Done } flag; } slot;
  13.  
  14. static int n_max = 0;         /* current maximum degree of subdivision */
  15. static slot *table = NULL;    /* an array of slots */
  16. static vertex *octant = NULL; /* the base octant of the ellipsoid */
  17.  
  18. #define SetP(p,px,py,pz) (p).x=(px), (p).y=(py), (p).z=(pz)
  19. #define SetV(v,px,py,pz,nx,ny,nz) SetP((v)->p,px,py,pz), SetP((v)->n,nx,ny,nz)
  20. #define SetF(f,i0,i1,i2) (f)->v0 = i0, (f)->v1 = i1, (f)->v2 = i2
  21.  
  22. /*
  23. // Compute the necessary cosine and sine values for generating ellipsoids
  24. // with the degree of subdivision n, and initialize the array table[].
  25. // The largest n becomes n_max, and calls with n <= n_max return immediately.
  26. // The memory for the base octant is allocated to cope with any n <= n_max.
  27. */
  28. void ellipsoid_init (int n)
  29. {
  30.    int n_table, i, j, k, l, m, h, d;
  31.    slot *t0, *t1, *t2;
  32.    float theta;
  33.  
  34.    if (n > n_max) {
  35.       n_max = n;
  36.       if (table) free (table);
  37.       if ((n_table = ((n-1)*n)/2) == 0) table = NULL;
  38.       else table = (slot *) malloc (n_table * sizeof(slot));
  39.       if (octant) free (octant);
  40.       octant = (vertex *) malloc (((n+1)*(n+2))/2 * sizeof(vertex));
  41.  
  42.       for (t0 = table, k = n_table; k > 0; k--, t0++) t0->flag = None;
  43.       for (t0 = table, k = 0, l = 1, m = 3, i = 2; i <= n_max; i++) {
  44.          l += m, m += 2, h = n_max / i - 1;
  45.          for (t1 = t0+i - 2, j = 1; j < i; j++, k++, t0++, t1--) {
  46.             if (t0->flag == None) {
  47.                theta = (M_PI_2 * j) / i;
  48.                t0->cos = t1->sin = cos (theta);
  49.                t0->sin = t1->cos = sin (theta);
  50.                t0->flag = t1->flag = Only;
  51.             }
  52.             if (t0->flag == Only) {
  53.                t0->flag = Done;
  54.                for (d = k+l, t2 = t0; h > 0; h--) {
  55.                   t2 += d, d += l;
  56.                   t2->cos = t0->cos;
  57.                   t2->sin = t0->sin;
  58.                   t2->flag = Done;
  59.                }
  60.             }
  61.          }
  62.       }
  63.    }
  64. }
  65.  
  66. /*
  67. // Construct the base octant of the ellipsoid whose parameters are a, b, and c,
  68. // with the degree of subdivision n using the cosine and sine values in table[].
  69. // It is assumed that n <= n_max.
  70. */
  71. static void ellipsoid_octant (int n, float a, float b, float c)
  72. {
  73.    int i, j;
  74.    float a_1, b_1, c_1;
  75.    float cos_ph, sin_ph, px, py, pz, nx, ny, nz, nznz, rnorm, tmp;
  76.    vertex *o = octant;
  77.    slot *table_th, *table_ph;
  78.  
  79.    a_1 = 1.0 / a; b_1 = 1.0 / b; c_1 = 1.0 / c;
  80.    o = octant;
  81.    table_th = table;
  82.    table_ph = table + ((n-1)*(n-2))/2;
  83.  
  84.    SetV (o, 0.0, 0.0, c, 0.0, 0.0, 1.0), o++;         /* i = 0, j = 0 */
  85.    for (i = 1; i < n; i++, table_ph++) {
  86.       cos_ph = table_ph->cos;
  87.       sin_ph = table_ph->sin;
  88.       pz = cos_ph * c;
  89.       nz = cos_ph * c_1;
  90.       nznz = nz * nz;
  91.  
  92.       px = sin_ph * a;
  93.       nx = sin_ph * a_1;
  94.       rnorm = 1.0 / sqrt (nx*nx + nznz);              /* 0 < i < n, j = 0 */
  95.       SetV (o, px, 0.0, pz, nx*rnorm, 0.0, nz*rnorm), o++;
  96.       for (j = i; --j > 0; table_th++) {
  97.          tmp = table_th->cos * sin_ph;
  98.          px = tmp * a;
  99.          nx = tmp * a_1;
  100.          tmp = table_th->sin * sin_ph;
  101.          py = tmp * b;
  102.          ny = tmp * b_1;
  103.          rnorm = 1.0 / sqrt (nx*nx + ny*ny + nznz);   /* 0 < i < n, 0 < j < i */
  104.          SetV (o, px, py, pz, nx*rnorm, ny*rnorm, nz*rnorm), o++;
  105.       }
  106.       py = sin_ph * b;
  107.       ny = sin_ph * b_1;
  108.       rnorm = 1.0 / sqrt (ny*ny + nznz);              /* 0 < i < n, j = i */
  109.       SetV (o, 0.0, py, pz, 0.0, ny*rnorm, nz*rnorm), o++;
  110.    }
  111.    SetV (o, a, 0.0, 0.0, 1.0, 0.0, 0.0), o++;         /* i = n, j = 0 */
  112.    for (j = i; --j > 0; table_th++) {
  113.       tmp = table_th->cos;
  114.       px = tmp * a;
  115.       nx = tmp * a_1;
  116.       tmp = table_th->sin;
  117.       py = tmp * b;
  118.       ny = tmp * b_1;
  119.       rnorm = 1.0 / sqrt (nx*nx + ny*ny);             /* i = n, 0 < j < i */
  120.       SetV (o, px, py, 0.0, nx*rnorm, ny*rnorm, 0.0), o++;
  121.    }
  122.    SetV (o, 0.0, b, 0.0, 0.0, 1.0, 0.0);              /* i = n, j = i */
  123. }
  124.  
  125. /*
  126. // Note the following conventions in ellipsoid_seq() and ellipsoid_par():
  127. // the north pole:        th =   0,      ph  =   0,
  128. // the 1st octant:   0 <= th <  90,  0 < ph <=  90,
  129. // the 2nd octant:  90 <= th < 180,  0 < ph <=  90,
  130. // the 3rd octant: 180 <= th < 270,  0 < ph <=  90,
  131. // the 4th octant: 270 <= th < 360,  0 < ph <=  90,
  132. // the 5th octant:   0 <= th <  90, 90 < ph <= 180,
  133. // the 6th octant:  90 <= th < 180, 90 < ph <= 180,
  134. // the 7th octant: 180 <= th < 270, 90 < ph <= 180,
  135. // the 8th octant: 270 <= th < 360, 90 < ph <= 180, and
  136. // the south pole:        th =   0,      ph  = 180.
  137. */
  138.  
  139. /*
  140. // Generate the vertices fot the ellipsoid with parameters a, b, and c
  141. // with the degree of subdivision n, by reflecting the base octant.
  142. // Also generate triangular faces of the ellipsoid with vertices ordered
  143. // counterclockwise when viewed from the outside.
  144. */
  145.  
  146. /* sequential version */
  147. void ellipsoid_seq (object *ellipsoid, int n, float a, float b, float c)
  148. {
  149.    vertex *v, *o;
  150.    face *f;
  151.    int i, j, ko, kv, kw, kv0, kw0;
  152.  
  153.    /* Check parameters for validity. */
  154.    if (n <= 0 || n_max < n || a <= 0.0 || b <= 0.0 || c <= 0.0) {
  155.       ellipsoid->nv = 0; ellipsoid->v = NULL;
  156.       ellipsoid->nf = 0; ellipsoid->f = NULL;
  157.       return;
  158.    }
  159.  
  160.    /* Initialize the base octant. */
  161.    ellipsoid_octant (n, a, b, c);
  162.  
  163.    /* Allocate memories for vertices and faces. */
  164.    ellipsoid->nv = 4*n*n + 2;
  165.    ellipsoid->nf = 8*n*n;
  166.    ellipsoid->v = (vertex *) malloc (ellipsoid->nv * sizeof(vertex));
  167.    ellipsoid->f = (face *) malloc (ellipsoid->nf * sizeof(face));
  168.  
  169.    /* Generate vertices of the ellipsoid from octant[]. */
  170.    v = ellipsoid->v;
  171.    o = octant;
  172. #define op o->p
  173. #define on o->n
  174.    SetV (v,  op.x,  op.y,  op.z,  on.x,  on.y,  on.z), v++; /* the north pole */
  175.    for (i = 0; ++i <= n;) {
  176.       o += i;
  177.       for (j = i; --j >= 0; o++, v++)                       /* 1st octant */
  178.          SetV (v,  op.x,  op.y,  op.z,  on.x,  on.y,  on.z);
  179.       for (j = i; --j >= 0; o--, v++)                       /* 2nd octant */
  180.          SetV (v, -op.x,  op.y,  op.z, -on.x,  on.y,  on.z);
  181.       for (j = i; --j >= 0; o++, v++)                       /* 3rd octant */
  182.          SetV (v, -op.x, -op.y,  op.z, -on.x, -on.y,  on.z);
  183.       for (j = i; --j >= 0; o--, v++)                       /* 4th octant */
  184.          SetV (v,  op.x, -op.y,  op.z,  on.x, -on.y,  on.z);
  185.    }
  186.    for (; --i > 1;) {
  187.       o -= i;
  188.       for (j = i; --j > 0; o++, v++)                        /* 5th octant */
  189.          SetV (v,  op.x,  op.y, -op.z,  on.x,  on.y, -on.z);
  190.       for (j = i; --j > 0; o--, v++)                        /* 6th octant */
  191.          SetV (v, -op.x,  op.y, -op.z, -on.x,  on.y, -on.z);
  192.       for (j = i; --j > 0; o++, v++)                        /* 7th octant */
  193.          SetV (v, -op.x, -op.y, -op.z, -on.x, -on.y, -on.z);
  194.       for (j = i; --j > 0; o--, v++)                        /* 8th octant */
  195.          SetV (v,  op.x, -op.y, -op.z,  on.x, -on.y, -on.z);
  196.    }
  197.    o--, SetV (v, -op.x, -op.y, -op.z, -on.x, -on.y, -on.z); /* the south pole */
  198. #undef op
  199. #undef on
  200.  
  201.    /* Generate triangular faces of the ellipsoid. */
  202.    f = ellipsoid->f;
  203.    kv = 0, kw = 1;
  204.    for (i = 0; i < n; i++) {
  205.       kv0 = kv, kw0 = kw;
  206.       for (ko = 1; ko <= 3; ko++)                /* the 1st, 2nd, 3rd octants */
  207.          for (j = i;; j--) {
  208.             SetF (f, kv, kw, ++kw), f++;
  209.             if (j == 0) break;
  210.             SetF (f, kv, kw, ++kv), f++;
  211.          }
  212.       for (j = i;; j--) {                        /* the 4th octant */
  213.          if (j == 0) { SetF (f, kv0, kw, kw0), kv++, kw++, f++; break; }
  214.          SetF (f, kv, kw, ++kw), f++;
  215.          if (j == 1) SetF (f, kv, kw, kv0), f++;
  216.          else SetF (f, kv, kw, ++kv), f++;
  217.       }
  218.    }
  219.    for (; --i >= 0;) {
  220.       kv0 = kv, kw0 = kw;
  221.       for (ko = 5; ko <= 7; ko++)                /* the 5th, 6th, 7th octants */
  222.          for (j = i;; j--) {
  223.             SetF (f, kv, kw, ++kv), f++;
  224.             if (j == 0) break;
  225.             SetF (f, kv, kw, ++kw), f++;
  226.          }
  227.       for (j = i;; j--) {                        /* the 8th octant */
  228.          if (j == 0) { SetF (f, kv, kw0, kv0), kv++, kw++, f++; break; }
  229.          SetF (f, kv, kw, ++kv), f++;
  230.          if (j == 1) SetF (f, kv, kw, kw0), f++;
  231.          else SetF (f, kv, kw, ++kw), f++;
  232.       }
  233.    }
  234. }
  235.  
  236. /* parallel version */
  237. void ellipsoid_par (object *ellipsoid, int n, float a, float b, float c)
  238. {
  239.    int nv, nf;
  240.    vertex *ev;
  241.    face *ef;
  242.  
  243.    /* Check parameters for validity. */
  244.    if (n <= 0 || n_max < n || a <= 0.0 || b <= 0.0 || c <= 0.0) {
  245.       ellipsoid->nv = 0; ellipsoid->v = NULL;
  246.       ellipsoid->nf = 0; ellipsoid->f = NULL;
  247.       return;
  248.    }
  249.  
  250.    /* Initialize the base octant. */
  251.    ellipsoid_octant (n, a, b, c);
  252.  
  253.    /* Allocate memories for vertices and faces. */
  254.    nv = ellipsoid->nv = 4*n*n + 2;
  255.    nf = ellipsoid->nf = 8*n*n;
  256.    ev = ellipsoid->v = (vertex *) malloc (nv * sizeof(vertex));
  257.    ef = ellipsoid->f = (face *) malloc (nf * sizeof(face));
  258.  
  259.    ellipsoid_par_help (n, nv, nf, ev, ef, octant);
  260. }
  261.  
  262. /* A utility function that simply outputs the ellipsoid data. */
  263. void ellipsoid_print (object *ellipsoid)
  264. {
  265.    int i;
  266.    vertex *v;
  267.    face *f;
  268.  
  269.    for (v = ellipsoid->v, i = 0; i < ellipsoid->nv; i++, v++)
  270.       printf ("v %g %g %g %g %g %g\n",
  271.          v->p.x, v->p.y, v->p.z, v->n.x, v->n.y, v->n.z);
  272.  
  273.    for (f = ellipsoid->f, i = 0; i < ellipsoid->nf; i++, f++)
  274.       printf ("f %d %d %d\n", f->v0, f->v1, f->v2);
  275. }
  276.  
  277. /* undo of ellipsoid_init () */
  278. void ellipsoid_done (void)
  279. {
  280.    n_max = 0;
  281.    if (table) free (table), table = NULL;
  282.    if (octant) free (octant), octant = NULL;
  283. }
  284.  
  285. /* undo of ellipsoid_seq () */
  286. void ellipsoid_free (object *ellipsoid)
  287. {
  288.    free (ellipsoid->v), ellipsoid->nv = 0, ellipsoid->v = NULL;
  289.    free (ellipsoid->f), ellipsoid->nf = 0, ellipsoid->f = NULL;
  290. }
  291.