
Date: Thu, 1 Mar 90 19:07:53 -0600
From: Pat Moran <moran@a.cs.uiuc.edu>
Message-Id: <9003020107.AA11204@p.cs.uiuc.edu>
To: irst!bellutta@a.cs.uiuc.edu, mcw@ukc.ac.uk
Subject: Sphere tesselation
Sender: moran <@p.cs.uiuc.edu:moran@a.cs.uiuc.edu>
Status: R


Included with this message is the algorithm that I mentioned in
comp.graphics.  I hope you find it useful.

Pat Moran
University of Illinois at Urbana-Champaign
moran@cs.uiuc.edu

==========
/*
 * NAME:  solidsample.c
 *
 * DESCRIPTION:
 *	"solidsample" generates the anglular part of vectors in a
 * spherical coordinate system with the goal of having the directions
 * evenly distributed over the coordinate space.  These angles are
 * intended to be used to control the sampling (gradient) directions
 * in a magnetic resonance system.
 *	The algorithm used below starts with a tetrahedron defined
 * by it vertices in spherical space.  Each face of the solid is then
 * recursively subdivided into four smaller triangles
 *	The convention used for the spherical coordinate system is to
 * describe the vector direction by two vectors, theta and phi.  In
 * terms of a Cartesian system, theta describes the angle down from the
 * positive z direction.  Phi gives the counter-clockwise angle in
 * the xy plane starting from the positive x direction.  This program
 * outputs one (theta, phi) pair per line, with the angles expressed in
 * degrees.
 *
 * HISTORY:
 *      90/02/07:  initial version
 */
#include <stdio.h>
#include <math.h>

#define EPSILON 0.000001
#define EQ0(x) (-(EPSILON) < (x)  &&  (x) < EPSILON)
#define ABS(x) ((x) >= 0 ? (x) : -(x))
#define DTOR(x) ((x) * M_PI / 180.0)
#define RTOD(x) ((x) * 180.0 / M_PI)
#define ACOS_1_3 109.4712207			/* Acos(-1/3) in degrees */

typedef struct {
  double theta, phi;
} sphereAngle;

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

typedef struct {
  vector3 v[3];
} face;

sphereAngle tetra[4] = {
  {0.0, 0.0}, {ACOS_1_3, 0.0}, {ACOS_1_3, 120.0}, {ACOS_1_3, 240.0}
};


void usage(name)
char *name;
{
  fprintf(stderr, "\nusage:\t%s n\t", name);
  fprintf(stderr, "(to generate 4**n entries)\n");
}


/* spherical angle assign */
void saAssign(dst, src)
sphereAngle *dst, *src;
{
  dst->theta = src->theta;
  dst->phi = src->phi;
}


/* copy a vector from src to dst (assignment) */
void vAssign(dst, src)
vector3 *dst, *src;
{
  dst->x = src->x;
  dst->y = src->y;
  dst->z = src->z;
}


/*   Normalize a vector, i.e. make its magnitude equal to 1 */
void normalize(v)
vector3 *v;
{
  double mag = sqrt(v->x*v->x + v->y*v->y + v->z*v->z);

  if (! EQ0(mag)) {
    v->x /= mag;
    v->y /= mag;
    v->z /= mag;
  }
}


/*   Given a spherical angle, return the corresponding unit vector */
void saToCart(sa, v)
sphereAngle *sa;
vector3 *v;
{
  double sinTheta = sin(DTOR(sa->theta));

  v->x = sinTheta * cos(DTOR(sa->phi));
  v->y = sinTheta * sin(DTOR(sa->phi));
  v->z = cos(DTOR(sa->theta));
}


/*   Given a vector, return the spherical angle.  Note that
 * this routine assumes that the vector is of unit magnitude.
 */
void cartToSa(sa, v)
sphereAngle *sa;
vector3 *v;
{
  sa->theta = RTOD(acos(v->z));		/* acos returns in range 0 to pi */

  if (EQ0(sa->theta))
    sa->phi = 0.0;
  else if (EQ0(180.0 - sa->theta))
    sa->phi = 0.0;
  else {
    sa->phi = RTOD(atan2(v->y, v->x));	/* atan2 returns in range -pi to pi */
    if (sa->phi < 0.0) sa->phi += 360.0;
  }
}


/*   "vBisector" returns the normalized bisector of two vectors.  It is
 * assumed that the input vectors are of equal magnitude (eg 1)
 */
void vBisector(bis, a, b)
vector3 *bis, *a, *b;
{
  bis->x = a->x + b->x;
  bis->y = a->y + b->y;
  bis->z = a->z + b->z;
  normalize(bis);
}


/* The average of 3 vectors */
void vAve3(ave, a, b, c)
vector3 *ave, *a, *b, *c;
{
  ave->x = a->x + b->x + c->x;
  ave->y = a->y + b->y + c->y;
  ave->z = a->z + b->z + c->z;
  normalize(ave);
}


double dotProd(a, b)
vector3 *a, *b;
{
  return a->x*b->x + a->y*b->y + a->z*b->z;
}


/*   faceNormal returns the normal vector to a face. */
void faceNormal(f, i, v)
face f[];
int i;
vector3 *v;
{
  vector3 ave;
  double x0 = f[i].v[0].x, y0 = f[i].v[0].y, z0 = f[i].v[0].z;
  double x1 = f[i].v[1].x, y1 = f[i].v[1].y, z1 = f[i].v[1].z;
  double x2 = f[i].v[2].x, y2 = f[i].v[2].y, z2 = f[i].v[2].z;

  v->x = y0 * (z1-z2) + y1 * (z2-z0) + y2 * (z0-z1);
  v->y = z0 * (x1-x2) + z1 * (x2-x0) + z2 * (x0-x1);
  v->z = x0 * (y1-y2) + x1 * (y2-y0) + x2 * (y0-y1);

  normalize(v);

  /* since the original points are not necessarily in counter-clockwise
   * direction, v may point 180 degrees in the wrong direction.  Compute
   * the dot product of v with the average of the three original vectors
   * in order to see if v is in the correct direction, and change it if
   * it is not.
   */
  vAve3(&ave, &f[i].v[0], &f[i].v[1], &f[i].v[2]);
  if (dotProd(v, &ave) < 0.0) {
    v->x = -v->x; v->y = -v->y; v->z = -v->z;
  }
}


/*   The function solidAngle returns the solid angle between two unit
 * vectors using equations derived from the dot product.
 */
double solidAngle(a, b)
vector3 *a, *b;
{
  return RTOD(acos(a->x*b->x + a->y*b->y + a->z*b->z));
}


/*   The vertices of the initial faces can be described as all the 
 * all the combinations from the 4 tetrahedron vertices, choosing 3.
 * Generates faces (t0,t1,t2), (t1,t2,t3), (t2,t3,t0) and (t3,t0,t1)
 * where tx is tetrahedron angle x.
 */
void initFaces(f)
face f[];
{
  int i, j;
  vector3 v;

  for (i=0; i < 4; i++)
    for (j=0; j < 3; j++) {
      saToCart(&tetra[(i+j)%4], &v);
      vAssign(&f[i].v[j], &v);
    }
}


void readArgs(argc, argv, n)
int argc;
char *argv[];
int *n;
{
  if (argc != 2) {
    usage(argv[0]);
    exit(1);
  }

  if (! sscanf(argv[1], "%d", n)) {
    (void) fprintf(stderr, "%s: Bad argument \"%d\".\n", argv[0], *n);
    usage(argv[0]);
    exit(1);
  }

  *n = ABS(*n);
  if (*n <= 0  ||  *n > 10) {
    (void) fprintf(stderr, "%s: Argument not in expected range.\n", argv[0]);
    usage(argv[0]);
    exit(1);
  }
}


/* power4 returns 4^n */
int power4(n)
{
  int result = 1;

  if (n < 0) return 0;
  while(n--) result <<= 2;
  return result;
}


/*   Subdivide takes as arguments the array of faces along with two
 * indices: the first specifying the old face and the second indicating
 * where the 3 newly generated faces will be stored.  
 *   The approach that this routine uses can be thought of as follows:
 *
 *		a0 *			  *
 *		  / \			 /u\
 *	      m2 +   + m1		+---+
 *		/     \		       /v\x/w\
 *	    a1 *---+---* a2	      *---+---+
 *		   m0
 * Points a0, a1 and a2 are the originals given for the
 * face at index oldf.  The new points m0, m1 and m2 are the midpoints
 * between the originals.  The new faces given by the solid angle triplets
 * (a0, m1, m2) (a1, m2, m0) and (a2, m0, m1) are stored in the face
 * array at newfs, newfs+1 and newfs+2 respectively.  In the diagram
 * above the faces correspond to u, v and w.  The face
 * specified by the triplet (m0, m1, m2) (face x in the diagram above)
 * is stored at oldf, where the original face used to be.  
 */
void subdivide(f, oldf, newfs)
face f[];
int oldf, newfs;
{
  int i;
  vector3 m[3];

  /* calculate the midpoints */
  for (i=0; i < 3; i++)
    vBisector(&m[i], &f[oldf].v[(i+1)%3], &f[oldf].v[(i+2)%3]);

  /* generate the three new outer faces */
  for (i=0; i < 3; i++) {
    vAssign(&f[newfs+i].v[0], &f[oldf].v[i]);
    vAssign(&f[newfs+i].v[1], &m[(i+1)%3]);
    vAssign(&f[newfs+i].v[2], &m[(i+2)%3]);
  }

  /* save the new middle face (m0, m1, m2) where the original used to be */
  for (i=0; i < 3; i++)
    vAssign(&f[oldf].v[i], &m[i]);
}
  

main(argc, argv)
int argc;
char *argv[];
{
  int i, n, totalFaces, pass, nfacesOld;
  face *f;
  vector3 v;
  sphereAngle norm;
  extern char *calloc();

  readArgs(argc, argv, &n);

  /* totalFaces = 4^n */
  totalFaces = power4(n);

  /* allocate space for table of faces */
  f = (face *) calloc((unsigned) totalFaces, sizeof(face));
  if (f == NULL) {
    (void) fprintf(stderr, "%s: Could not calloc enough space.\n", argv[0]);
    exit(1);
  }

  /* initialize array with tetrahedron faces */
  initFaces(f);

  /* subdivide faces */
  for (pass=1; pass < n; pass++) {
    nfacesOld = power4(pass);
    for (i=0; i < nfacesOld; i++)
      subdivide(f, i, nfacesOld + 3*i);
  }

  /* output the Cartesian coordinates for the corners of each face */
  /*for (i=0; i < totalFaces; i++) {
    printf("%d: (%6.3lf,%6.3lf,%6.3lf), ",
	   i, f[i].v[0].x, f[i].v[0].y, f[i].v[0].z);
    printf("(%6.3lf,%6.3lf,%6.3lf), ", f[i].v[1].x, f[i].v[1].y, f[i].v[1].z);
    printf("(%6.3lf,%6.3lf,%6.3lf)\n", f[i].v[2].x, f[i].v[2].y, f[i].v[2].z);
  } */

  /* output the spherical angles */
  for (i=0; i < totalFaces; i++) {
    faceNormal(f, i, &v);
    cartToSa(&norm, &v);
    printf("%6.2lf %6.2lf\n", norm.theta, norm.phi);
  }
}