/*
 * 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 EXP_FACTOR	0.1
#define SPIN		60
#define WOBBLE_FACTOR	0.15
#define EXP_WOBBLE	0.2
#define FLY_OFF		6
#define FLY_AMP		8.0

#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]);
}
  
float
real_rand(scale)
float	scale;
{
	return (1.0 + ((float)(random() % 1000000) - 500000) / 500000.0 * scale);
}

main(argc, argv)
int argc;
char *argv[];
{
	int i, j, n, totalFaces, pass, nfacesOld;
	int plane_no;
	int sphere_no;
	float tot_x = 0, tot_y = 0, tot_z = 0;
	int bodge;

	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);
	  } */

	plane_no = 0;
	sphere_no = 0;

	srandom(totalFaces);

	for (i = 0; i < totalFaces; i++) {
		int	rev_mode;
		float	mid_x,
			mid_y,
			mid_z,
			x_dist,
			y_dist,
			z_dist,
			radius,
			bsphere_radius,
			exp_factor,
			spin_scale;

		tot_x = tot_y = tot_z = 0.0;
		for (j = 0; j < 3; j++) {
			tot_x += f[i].v[j].x;
			tot_y += f[i].v[j].y;
			tot_z += f[i].v[j].z;
		}
		mid_x = tot_x / 3.0;
		mid_y = tot_y / 3.0;
		mid_z = tot_z / 3.0;
		

		bsphere_radius = 0.0;
		for (j = 0; j < 3; j++) {
			x_dist = f[i].v[j].x - mid_x;
			y_dist = f[i].v[j].y - mid_y;
			z_dist = f[i].v[j].z - mid_z;
			radius = sqrt(x_dist * x_dist +
				      y_dist * y_dist +
				      z_dist * z_dist);
			if (bsphere_radius < radius) {
				bsphere_radius = radius;
			}
		}

		for (j = 0; j < 3; j++) {
			f[i].v[j].x *= real_rand(bsphere_radius * WOBBLE_FACTOR);
			f[i].v[j].y *= real_rand(bsphere_radius * WOBBLE_FACTOR);
			f[i].v[j].z *= real_rand(bsphere_radius * WOBBLE_FACTOR);
		}

		bsphere_radius = 0.0;
		for (j = 0; j < 3; j++) {
			x_dist = f[i].v[j].x - mid_x;
			y_dist = f[i].v[j].y - mid_y;
			z_dist = f[i].v[j].z - mid_z;
			radius = sqrt(x_dist * x_dist +
				      y_dist * y_dist +
				      z_dist * z_dist);
			if (bsphere_radius < radius) {
				bsphere_radius = radius;
			}
		}

		exp_factor = 1.0;
		spin_scale = 1.0;

		if (random() % FLY_OFF == 0) {
			float	extra =
				 ((float)(random() % 100000) / 100000) * FLY_AMP;

			exp_factor *= extra;
			spin_scale *= extra;
		}

		tot_x = mid_x * EXP_FACTOR * fabs(real_rand(EXP_WOBBLE)) *
				exp_factor;
		tot_y = mid_y * EXP_FACTOR * fabs(real_rand(EXP_WOBBLE)) *
				exp_factor;
		tot_z = mid_z * EXP_FACTOR * fabs(real_rand(EXP_WOBBLE)) *
				exp_factor;

		rev_mode = (i % 8) >= 4;

		printf("SOLID o%d {\n", i);
		printf("MEET_LIST { sphere\n");
		for (j = 0; j < 3; j++) {
			int k;

			printf("  { PLANE4 0.0 0.0 0.0 ");
			for (k = 0; k < 2; k++) {
				printf("%1.5f %1.5f %1.5f", f[i].v[(j+k)%3].x, f[i].v[(j+k)%3].y, f[i].v[(j+k)%3].z);
				putchar(' ');
			}
			printf("%1.5f %1.5f %1.5f }\n", f[i].v[(j+k)%3].x, f[i].v[(j+k)%3].y, f[i].v[(j+k)%3].z);
		}
		printf("  { SPHERE SCALE %1.5f %1.5f %1.5f TRANSLATE %1.5f %1.5f %1.5f }\n", bsphere_radius, bsphere_radius, bsphere_radius, mid_x, mid_y, mid_z);
		printf("} TRANSLATE %1.5f %1.5f %1.5f ROTATE X_AXIS %1.5f DEG ROTATE Y_AXIS %1.5f DEG\n  TRANSLATE %1.5f %1.5f %1.5f SCALE 2 2 2 BBOX }\n\n",
			-mid_x, -mid_y, -mid_z,
			(float) (random() % SPIN - (SPIN / 2)) * spin_scale, (float)(random() % SPIN - (SPIN / 2)) * spin_scale,
			mid_x + tot_x, mid_y + tot_y, mid_z + tot_z);
	}

	printf("SOLID thing {\n  JOIN_LIST { ");
	for (i = 0; i < totalFaces; i++) {
		printf("o%d ", i);
	}
	puts("} SURFACE thing_col TRANSLATE 0 8 0\n}");
}