From leech@alanine.UUCP Fri Mar 24 18:29:14 1989
Path: leah!csd4.milw.wisc.edu!mailrus!tut.cis.ohio-state.edu!rutgers!mcnc!thorin!alanine!leech
From: leech@alanine.cs.unc.edu (Jonathan Leech)
Newsgroups: comp.graphics
Subject: Re: Polygon Representation of a Sphere's Surface
Message-ID: <7409@thorin.cs.unc.edu>
Date: 24 Mar 89 23:29:14 GMT
References: <270@ai.etl.army.mil>
Sender: news@thorin.cs.unc.edu
Reply-To: leech@alanine.UUCP (Jonathan Leech)
Organization: University Of North Carolina, Chapel Hill
Lines: 339
Summary:
Expires:
Sender:
Followup-To:
Distribution:
Keywords:

In article <270@ai.etl.army.mil> richr@ai.etl.army.mil. (Richard Rosenthal) writes:
>I want to have in 3-D space (x, y, z) a polygon near-representation
>of the surface of a sphere where each polygon is identical
>and regular (if that's the word).
>...
>Any comments, algorithms, or reference suggestions?

    This request comes up (seemingly) once a month or so.  For a
change, here's the code to do it instead of a description.  Note that
the triangles generated are all equilateral, but not all the same
size. This program generates tesselations with

    #triangles(n'th approximation) = 8 * 4^(n-1),

since the starting point is an octahedron. Starting with other
platonic solids would change the coefficient, of course.

    The program retains generated triangles in core. It could be
easily modified to print out results on the fly if memory is limited.

    Code modifications will be required to generate your preferred
database format - see the comments.

P.S. Eugene, or whoever is keeping the commonly asked questions list,
how about adding this?

/*
 * sphere - generate a polygon mesh approximating a sphere by
 *  recursive subdivision. First approximation is an octahedron;
 *  each level of refinement increases the number of polygons by
 *  a factor of 4.
 * Level 3 (128 polygons) is a good tradeoff if gouraud
 *  shading is used to render the database.
 *
 * Usage: sphere [level]
 *	level is an integer >= 1 setting the recursion level (default 1).
 *
 * Notes:
 *
 *  The triangles are generated with vertices in clockwise order as
 *  viewed from the outside in a right-handed coordinate system.
 *  To reverse the order, compile with COUNTERCLOCKWISE defined.
 *
 *  Shared vertices are not retained, so numerical errors may produce
 *  cracks between polygons at high subdivision levels.
 *
 *  The subroutines print_object() and print_triangle() should
 *  be changed to generate whatever the desired database format is.
 *  If UNC is defined, a PPHIGS text archive is generated.
 *
 * Jon Leech 3/24/89
 */
#include 
#include 

typedef struct {
    double  x, y, z;
} point;

typedef struct {
    point     pt[3];	/* Vertices of triangle */
    double    area;	/* Unused; might be used for adaptive subdivision */
} triangle;

typedef struct {
    int       npoly;	/* # of polygons in object */
    triangle *poly;	/* Polygons in no particular order */
} object;

/* Six equidistant points lying on the unit sphere */
#define XPLUS {  1,  0,  0 }	/*  X */
#define XMIN  { -1,  0,  0 }	/* -X */
#define YPLUS {  0,  1,  0 }	/*  Y */
#define YMIN  {  0, -1,  0 }	/* -Y */
#define ZPLUS {  0,  0,  1 }	/*  Z */
#define ZMIN  {  0,  0, -1 }	/* -Z */

/* Vertices of a unit octahedron */
triangle octahedron[] = {
    { XPLUS, ZPLUS, YPLUS }, 0.0,
    { YPLUS, ZPLUS, XMIN  }, 0.0,
    { XMIN , ZPLUS, YMIN  }, 0.0,
    { YMIN , ZPLUS, XPLUS }, 0.0,
    { XPLUS, YPLUS, ZMIN  }, 0.0,
    { YPLUS, XMIN , ZMIN  }, 0.0,
    { XMIN , YMIN , ZMIN  }, 0.0,
    { YMIN , XPLUS, ZMIN  }, 0.0
};

/* An octahedron */
object oct = {
    sizeof(octahedron) / sizeof(octahedron[0]),
    &octahedron[0]
};

/* Forward declarations */
point *normalize(/* point *p */);
point *midpoint(/* point *a, point *b */);
void print_object(/* object *obj, int level */);
void print_triangle(/* triangle *t */);
double sqr(/* double x */);
double area_of_triangle(/* triangle *t */);

extern char *malloc(/* unsigned */);

main(ac, av)
int ac;
char *av[];
{
    object *old, *new;
    int     i, level, maxlevels = 1;

    if (ac > 1)
	if ((maxlevels = atoi(av[1])) < 1) {
	    fprintf(stderr, "%s: # of levels must be >= 1\n", av[0]);
	    exit(1);
	}

#ifdef COUNTERCLOCKWISE
    /* Reverse order of points in each triangle */
    for (i = 0; i < oct.npoly; i++) {
	point tmp;
		      tmp = oct.poly[i].pt[0];
	oct.poly[i].pt[0] = oct.poly[i].pt[2];
	oct.poly[i].pt[2] = tmp;
    }
#endif

    old = &oct;

    /* Subdivide each starting triangle (maxlevels - 1) times */
    for (level = 1; level < maxlevels; level++) {
	/* Allocate a new object */
	new = (object *)malloc(sizeof(object));
	if (new == NULL) {
	    fprintf(stderr, "%s: Out of memory on subdivision level %d\n",
		av[0], level);
	    exit(1);
	}
	new->npoly = old->npoly * 4;

	/* Allocate 4* the number of points in the current approximation */
	new->poly  = (triangle *)malloc(new->npoly * sizeof(triangle));
	if (new->poly == NULL) {
	    fprintf(stderr, "%s: Out of memory on subdivision level %d\n",
		av[0], level);
	    exit(1);
	}

	/* Subdivide each polygon in the old approximation and normalize
	 *  the new points thus generated to lie on the surface of the unit
	 *  sphere.
	 * Each input triangle with vertices labelled [0,1,2] as shown
	 *  below will be turned into four new triangles:
	 *
	 *			Make new points
	 *			    a = (0+2)/2
	 *			    b = (0+1)/2
	 *			    c = (1+2)/2
	 *	  1
	 *	 /\		Normalize a, b, c
	 *	/  \
	 *    b/____\ c		Construct new triangles
	 *    /\    /\		    [0,b,a]
	 *   /	\  /  \		    [b,1,c]
	 *  /____\/____\	    [a,b,c]
	 * 0	  a	2	    [a,c,2]
	 */
	for (i = 0; i < old->npoly; i++) {
	    triangle
		 *oldt = &old->poly[i],
		 *newt = &new->poly[i*4];
	    point a, b, c;

	    a = *normalize(midpoint(&oldt->pt[0], &oldt->pt[2]));
	    b = *normalize(midpoint(&oldt->pt[0], &oldt->pt[1]));
	    c = *normalize(midpoint(&oldt->pt[1], &oldt->pt[2]));

	    newt->pt[0] = oldt->pt[0];
	    newt->pt[1] = b;
	    newt->pt[2] = a;
	    newt++;

	    newt->pt[0] = b;
	    newt->pt[1] = oldt->pt[1];
	    newt->pt[2] = c;
	    newt++;

	    newt->pt[0] = a;
	    newt->pt[1] = b;
	    newt->pt[2] = c;
	    newt++;

	    newt->pt[0] = a;
	    newt->pt[1] = c;
	    newt->pt[2] = oldt->pt[2];
	}

	if (level > 1) {
	    free(old->poly);
	    free(old);
	}

	/* Continue subdividing new triangles */
	old = new;
    }

    /* Print out resulting approximation */
    print_object(old, maxlevels);
}

/* Normalize a point p */
point *normalize(p)
point *p;
{
    static point r;
    double mag;

    r = *p;
    mag = r.x * r.x + r.y * r.y + r.z * r.z;
    if (mag != 0.0) {
	mag = 1.0 / sqrt(mag);
	r.x *= mag;
	r.y *= mag;
	r.z *= mag;
    }

    return &r;
}

/* Return the average of two points */
point *midpoint(a, b)
point *a, *b;
{
    static point r;

    r.x = (a->x + b->x) * 0.5;
    r.y = (a->y + b->y) * 0.5;
    r.z = (a->z + b->z) * 0.5;

    return &r;
}

/* Write out all polygons in an object */
void print_object(obj, level)
object *obj;
int level;
{
    int i;

#ifdef UNC  /* Generate PPHIGS archive format */
    int dx, dy, dz;

    printf("structure sphere%d posted {\n", level);
    printf("\tcolor polygon {\n");
    printf("\t\t200 100  50   0     50 100 200   0\n");
    printf("\t};\n");

    switch (level) {
	case 1:
	    dx = -2000; dy =  2000; dz = 0;
	    break;
	case 2:
	    dx =  2000; dy =  2000; dz = 0;
	    break;
	case 3:
	    dx = -2000; dy = -2000; dz = 0;
	    break;
	case 4:
	    dx =  2000; dy = -2000; dz = 0;
	    break;
	case 5:
	    dx =     0; dy =	 0; dz = 0;
	    break;
	default:
	    dx = dy = dz = 0;
	    break;
    }

    printf("\tmatrix Pre scale 1000 1000 1000;\n");
    printf("\tmatrix Pre translate %d %d %d ;\n", dx, dy, dz);
#endif	/* UNC */

    /* Spit out coordinates for each triangle */
    for (i = 0; i < obj->npoly; i++)
	print_triangle(&obj->poly[i]);

#ifdef UNC
    printf("};\n");
#endif	/* UNC */
}

/* Write a single polygon */
void print_triangle(t)
triangle *t;
{
    int i;

#ifdef UNC
    printf("\tpolygon 3 {\n");
    for (i = 0; i < 3; i++)
	printf("\t\t%g\t%g\t%g %g\t%g\t%g ;\n",
	    t->pt[i].x, t->pt[i].y, t->pt[i].z,     /* Point */
	    t->pt[i].x, t->pt[i].y, t->pt[i].z);    /* Normal */
    printf("\t};\n");
#else
    /* Insert code for your output format here.
     * Triangle vertices are in t->pt[0..2].{x,y,z}
     */
    static int errflag = 0;
    fprintf(stderr, "Error! You must modify print_triangle() to generate triangles in\n");
    fprintf(stderr, "\tthe desired format!\n");
    exit(1);
#endif
}

double sqr(x) double x; { return x * x ; }

/* Returns area of triangle (currently unused) */
double area_of_triangle(t)
triangle *t;
{
    double a;

    a  = sqr(t->pt[0].y * (t->pt[1].z - t->pt[2].z) -
	     t->pt[1].y * (t->pt[0].z - t->pt[2].z) +
	     t->pt[2].y * (t->pt[0].z - t->pt[1].z));
    a += sqr(t->pt[0].z * (t->pt[1].x - t->pt[2].x) -
	     t->pt[1].z * (t->pt[0].x - t->pt[2].x) +
	     t->pt[2].z * (t->pt[0].x - t->pt[1].x));
    a += sqr(t->pt[0].x * (t->pt[1].y - t->pt[2].y) -
	     t->pt[1].x * (t->pt[0].y - t->pt[2].y) +
	     t->pt[2].x * (t->pt[0].y - t->pt[1].y));

    return 0.5 * sqrt(a);
}
--
    Jon Leech (leech@cs.unc.edu)    __@/
    SUSHIDO: the Way of the Tuna

Discuss this article in the forums


Date this article was posted to GameDev.net: 7/16/1999
(Note that this date does not necessarily correspond to the date the article was written)

See Also:
Miscellaneous

© 1999-2011 Gamedev.net. All rights reserved. Terms of Use Privacy Policy
Comments? Questions? Feedback? Click here!