home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Assembler / dse-src3.dms / in.adf / Source / Vectors / Sphere.lha / sphere.c
Encoding:
C/C++ Source or Header  |  1993-01-15  |  7.1 KB  |  307 lines
  1. /*% cc -g sphere.c -o sphere -lm
  2.  *
  3.  * sphere - generate a triangle mesh approximating a sphere by
  4.  *  recursive subdivision. First approximation is an octahedron;
  5.  *  each level of refinement increases the number of triangles by
  6.  *  a factor of 4.
  7.  * Level 3 (128 triangles) is a good tradeoff if gouraud
  8.  *  shading is used to render the database.
  9.  *
  10.  * Usage: sphere [level] [-p] [-c]
  11.  *    level is an integer >= 1 setting the recursion level (default 1).
  12.  *    -p causes generation of a PPHIGS format ASCII archive
  13.  *        instead of the default generic output format.
  14.  *    -c causes triangles to be generated with vertices in counterclockwise
  15.  *        order as viewed from the outside in a RHS coordinate system.
  16.  *        The default is clockwise order.
  17.  *
  18.  *  The subroutines print_object() and print_triangle() should
  19.  *  be changed to generate whatever the desired database format is.
  20.  *
  21.  * Jon Leech (leech@cs.unc.edu) 3/24/89
  22.  */
  23. #include <stdio.h>
  24. #include <math.h>
  25.  
  26. typedef struct {
  27.     double  x, y, z;
  28. } point;
  29.  
  30. typedef struct {
  31.     point     pt[3];    /* Vertices of triangle */
  32.     double    area;    /* Unused; might be used for adaptive subdivision */
  33. } triangle;
  34.  
  35. typedef struct {
  36.     int       npoly;    /* # of triangles in object */
  37.     triangle *poly;    /* Triangles */
  38. } object;
  39.  
  40. /* Six equidistant points lying on the unit sphere */
  41. #define XPLUS {  1,  0,  0 }    /*  X */
  42. #define XMIN  { -1,  0,  0 }    /* -X */
  43. #define YPLUS {  0,  1,  0 }    /*  Y */
  44. #define YMIN  {  0, -1,  0 }    /* -Y */
  45. #define ZPLUS {  0,  0,  1 }    /*  Z */
  46. #define ZMIN  {  0,  0, -1 }    /* -Z */
  47.  
  48. /* Vertices of a unit octahedron */
  49. triangle octahedron[] = {
  50.     { { XPLUS, ZPLUS, YPLUS }, 0.0 },
  51.     { { YPLUS, ZPLUS, XMIN  }, 0.0 },
  52.     { { XMIN , ZPLUS, YMIN  }, 0.0 },
  53.     { { YMIN , ZPLUS, XPLUS }, 0.0 },
  54.     { { XPLUS, YPLUS, ZMIN  }, 0.0 },
  55.     { { YPLUS, XMIN , ZMIN  }, 0.0 },
  56.     { { XMIN , YMIN , ZMIN  }, 0.0 },
  57.     { { YMIN , XPLUS, ZMIN  }, 0.0 }
  58. };
  59.  
  60. /* A unit octahedron */
  61. object oct = {
  62.     sizeof(octahedron) / sizeof(octahedron[0]),
  63.     &octahedron[0]
  64. };
  65.  
  66. int PPHIGSflag = 0;
  67.  
  68. /* Forward declarations */
  69. point *normalize(/* point *p */);
  70. point *midpoint(/* point *a, point *b */);
  71. void print_object(/* object *obj, int level */);
  72. void print_triangle(/* triangle *t */);
  73. void pphigs_header(/* int */);
  74. void pphigs_trailer();
  75.  
  76. extern char *malloc(/* unsigned */);
  77.  
  78. main(ac, av)
  79. int ac;
  80. char *av[];
  81. {
  82.     object *old,
  83.        *new;
  84.     int     ccwflag = 0,    /* Reverse vertex order if true */
  85.         i,
  86.         level,        /* Current subdivision level */
  87.         maxlevel = 1;    /* Maximum subdivision level */
  88.  
  89.     /* Parse arguments */
  90.     for (i = 1; i < ac; i++) {
  91.     if (!strcmp(av[i], "-p"))
  92.         PPHIGSflag = 1;
  93.     else if (!strcmp(av[i], "-c"))
  94.         ccwflag = 1;
  95.     else if ((maxlevel = atoi(av[i])) < 1) {
  96.         fprintf(stderr, "%s: # of levels must be >= 1\n", av[0]);
  97.         exit(1);
  98.     }
  99.     }
  100.  
  101.     if (ccwflag) {
  102.     /* Reverse order of points in each triangle */
  103.     for (i = 0; i < oct.npoly; i++) {
  104.         point tmp;
  105.               tmp = oct.poly[i].pt[0];
  106.         oct.poly[i].pt[0] = oct.poly[i].pt[2];
  107.         oct.poly[i].pt[2] = tmp;
  108.     }
  109.     }
  110.  
  111.     old = &oct;
  112.  
  113.     /* Subdivide each starting triangle (maxlevel - 1) times */
  114.     for (level = 1; level < maxlevel; level++) {
  115.     /* Allocate a new object */
  116.     new = (object *)malloc(sizeof(object));
  117.     if (new == NULL) {
  118.         fprintf(stderr, "%s: Out of memory on subdivision level %d\n",
  119.         av[0], level);
  120.         exit(1);
  121.     }
  122.     new->npoly = old->npoly * 4;
  123.  
  124.     /* Allocate 4* the number of points in the current approximation */
  125.     new->poly  = (triangle *)malloc(new->npoly * sizeof(triangle));
  126.     if (new->poly == NULL) {
  127.         fprintf(stderr, "%s: Out of memory on subdivision level %d\n",
  128.         av[0], level);
  129.         exit(1);
  130.     }
  131.  
  132.     /* Subdivide each triangle in the old approximation and normalize
  133.      *  the new points thus generated to lie on the surface of the unit
  134.      *  sphere.
  135.      * Each input triangle with vertices labelled [0,1,2] as shown
  136.      *  below will be turned into four new triangles:
  137.      *
  138.      *            Make new points
  139.      *                a = (0+2)/2
  140.      *                b = (0+1)/2
  141.      *                c = (1+2)/2
  142.      *      1
  143.      *     /\        Normalize a, b, c
  144.      *    /  \
  145.      *    b/____\ c        Construct new triangles
  146.      *    /\    /\            [0,b,a]
  147.      *   /    \  /  \            [b,1,c]
  148.      *  /____\/____\        [a,b,c]
  149.      * 0      a    2        [a,c,2]
  150.      */
  151.     for (i = 0; i < old->npoly; i++) {
  152.         triangle
  153.          *oldt = &old->poly[i],
  154.          *newt = &new->poly[i*4];
  155.         point a, b, c;
  156.  
  157.         a = *normalize(midpoint(&oldt->pt[0], &oldt->pt[2]));
  158.         b = *normalize(midpoint(&oldt->pt[0], &oldt->pt[1]));
  159.         c = *normalize(midpoint(&oldt->pt[1], &oldt->pt[2]));
  160.  
  161.         newt->pt[0] = oldt->pt[0];
  162.         newt->pt[1] = b;
  163.         newt->pt[2] = a;
  164.         newt++;
  165.  
  166.         newt->pt[0] = b;
  167.         newt->pt[1] = oldt->pt[1];
  168.         newt->pt[2] = c;
  169.         newt++;
  170.  
  171.         newt->pt[0] = a;
  172.         newt->pt[1] = b;
  173.         newt->pt[2] = c;
  174.         newt++;
  175.  
  176.         newt->pt[0] = a;
  177.         newt->pt[1] = c;
  178.         newt->pt[2] = oldt->pt[2];
  179.     }
  180.  
  181.     if (level > 1) {
  182.         free(old->poly);
  183.         free(old);
  184.     }
  185.  
  186.     /* Continue subdividing new triangles */
  187.     old = new;
  188.     }
  189.  
  190.     /* Print out resulting approximation */
  191.     print_object(old, maxlevel);
  192. }
  193.  
  194. /* Normalize a point p */
  195. point *normalize(p)
  196. point *p;
  197. {
  198.     static point r;
  199.     double mag;
  200.  
  201.     r = *p;
  202.     mag = r.x * r.x + r.y * r.y + r.z * r.z;
  203.     if (mag != 0.0) {
  204.     mag = 1.0 / sqrt(mag);
  205.     r.x *= mag;
  206.     r.y *= mag;
  207.     r.z *= mag;
  208.     }
  209.  
  210.     return &r;
  211. }
  212.  
  213. /* Return the midpoint on the line between two points */
  214. point *midpoint(a, b)
  215. point *a, *b;
  216. {
  217.     static point r;
  218.  
  219.     r.x = (a->x + b->x) * 0.5;
  220.     r.y = (a->y + b->y) * 0.5;
  221.     r.z = (a->z + b->z) * 0.5;
  222.  
  223.     return &r;
  224. }
  225.  
  226. /* Write out all triangles in an object */
  227. void print_object(obj, level)
  228. object *obj;
  229. int level;
  230. {
  231.     int i;
  232.  
  233.     if (PPHIGSflag)
  234.     pphigs_header(level);
  235.  
  236.     /* Spit out coordinates for each triangle */
  237.     for (i = 0; i < obj->npoly; i++)
  238.     print_triangle(&obj->poly[i]);
  239.  
  240.     if (PPHIGSflag)
  241.     pphigs_trailer();
  242. }
  243.  
  244. /* Output a triangle */
  245. void print_triangle(t)
  246. triangle *t;
  247. {
  248.     int i;
  249.  
  250.     if (PPHIGSflag) {
  251.     printf("\tpolygon 3 {\n");
  252.     for (i = 0; i < 3; i++)
  253.         printf("\t\t%g\t%g\t%g %g\t%g\t%g ;\n",
  254.         t->pt[i].x, t->pt[i].y, t->pt[i].z,    /* Point */
  255.         t->pt[i].x, t->pt[i].y, t->pt[i].z);    /* Normal */
  256.     printf("\t};\n");
  257.     } else {
  258.     /* Modify this to generate your favorite output format
  259.      * Triangle vertices are in t->pt[0..2].{x,y,z}
  260.      * A generic format is provided.
  261.      */
  262.     printf("triangle\n");
  263.     for (i = 0; i < 3; i++)
  264.         printf("\t%g %g %g\n", t->pt[i].x, t->pt[i].y, t->pt[i].z);
  265.     }
  266. }
  267.  
  268. /* Generate header for a PPHIGS ASCII archive */
  269. void pphigs_header(level)
  270. int level;
  271. {
  272.     int dx, dy, dz;
  273.  
  274.     printf("structure sphere%d posted {\n", level);
  275.     printf("\tcolor polygon {\n");
  276.     printf("\t\t200 100  50   0     50 100 200   0\n");
  277.     printf("\t};\n");
  278.  
  279.     switch (level) {
  280.     case 1:
  281.         dx = -2000; dy =  2000; dz = 0;
  282.         break;
  283.     case 2:
  284.         dx =  2000; dy =  2000; dz = 0;
  285.         break;
  286.     case 3:
  287.         dx = -2000; dy = -2000; dz = 0;
  288.         break;
  289.     case 4:
  290.         dx =  2000; dy = -2000; dz = 0;
  291.         break;
  292.     case 5:
  293.         dx =     0; dy =     0; dz = 0;
  294.         break;
  295.     default:
  296.         dx = dy = dz = 0;
  297.         break;
  298.     }
  299.  
  300.     printf("\tmatrix Pre scale 1000 1000 1000;\n");
  301.     printf("\tmatrix Pre translate %d %d %d ;\n", dx, dy, dz);
  302. }
  303.  
  304. void pphigs_trailer() {
  305.     printf("};\n");
  306. }
  307.