/*% cc -g sphere.c -o sphere -lm
 *
 * sphere - generate a triangle mesh approximating a sphere by
 *  recursive subdivision. First approximation is an platonic
 *  solid; each level of refinement increases the number of
 *  triangles by a factor of 4.
 *
 * Level 3 (128 triangles for an octahedron) is a good tradeoff if
 *  gouraud shading is used to render the database.
 *
 * Usage: sphere [level] [-p] [-c] [-f] [-t] [-i]
 *	level is an integer >= 1 setting the recursion level (default 1).
 *	-p causes generation of a PPHIGS format ASCII archive
 *	    instead of the default generic output format.
 *	-c causes triangles to be generated with vertices in counterclockwise
 *	    order as viewed from the outside in a RHS coordinate system.
 *	    The default is clockwise order.
 *	-f generates triangle without per-vertex normals (PPHIGS only)
 *	-t starts with a tetrahedron instead of an octahedron
 *	-i starts with a icosahedron instead of an octahedron
 *
 *  The subroutines print_object() and print_triangle() should
 *  be changed to generate whatever the desired database format is.
 *
 * Jon Leech (leech@cs.unc.edu) 3/24/89
 * icosahedral code added by Jim Buddenhagen (jb1556@daditz.sbc.com) 5/93
 */
#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include <string.h>

#define PI	3.1415926535897932384626433832795
#define radtodeg	180/PI

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 triangles in object */
    triangle *poly;	/* Triangles */
} 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 }
};

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

/* Vertices of a tetrahedron */
#define sqrt_3 0.5773502692
#define PPP {  sqrt_3,	sqrt_3,  sqrt_3 }   /* +X, +Y, +Z */
#define MMP { -sqrt_3, -sqrt_3,  sqrt_3 }   /* -X, -Y, +Z */
#define MPM { -sqrt_3,	sqrt_3, -sqrt_3 }   /* -X, +Y, -Z */
#define PMM {  sqrt_3, -sqrt_3, -sqrt_3 }   /* +X, -Y, -Z */

/* Structure describing a tetrahedron */
triangle tetrahedron[] = {
    {{ PPP, MMP, MPM }, 0.0},
    {{ PPP, PMM, MMP }, 0.0},
    {{ MPM, MMP, PMM }, 0.0},
    {{ PMM, PPP, MPM }, 0.0}
};

object tet = {
    sizeof(tetrahedron) / sizeof(tetrahedron[0]),
    &tetrahedron[0]
};

/* Twelve vertices of icosahedron on unit sphere */
#define tau 0.8506508084      /* t=(1+sqrt(5))/2, tau=t/sqrt(1+t^2)  */
#define one 0.5257311121      /* one=1/sqrt(1+t^2) , unit sphere     */
#define ZA {  tau,  one,    0 }
#define ZB { -tau,  one,    0 }
#define ZC { -tau, -one,    0 }
#define ZD {  tau, -one,    0 }
#define YA {  one,   0 ,  tau }
#define YB {  one,   0 , -tau }
#define YC { -one,   0 , -tau }
#define YD { -one,   0 ,  tau }
#define XA {   0 ,  tau,  one }
#define XB {   0 , -tau,  one }
#define XC {   0 , -tau, -one }
#define XD {   0 ,  tau, -one }

/* Structure for unit icosahedron */
triangle icosahedron[] = {
    { { YA, XA, YD }, 0.0 },
    { { YA, YD, XB }, 0.0 },
    { { YB, YC, XD }, 0.0 },
    { { YB, XC, YC }, 0.0 },
    { { ZA, YA, ZD }, 0.0 },
    { { ZA, ZD, YB }, 0.0 },
    { { ZC, YD, ZB }, 0.0 },
    { { ZC, ZB, YC }, 0.0 },
    { { XA, ZA, XD }, 0.0 },
    { { XA, XD, ZB }, 0.0 },
    { { XB, XC, ZD }, 0.0 },
    { { XB, ZC, XC }, 0.0 },
    { { XA, YA, ZA }, 0.0 },
    { { XD, ZA, YB }, 0.0 },
    { { YA, XB, ZD }, 0.0 },
    { { YB, ZD, XC }, 0.0 },
    { { YD, XA, ZB }, 0.0 },
    { { YC, ZB, XD }, 0.0 },
    { { YD, ZC, XB }, 0.0 },
    { { YC, XC, ZC }, 0.0 }
};

/* A unit icosahedron */
object ico = {
    sizeof(icosahedron) / sizeof(icosahedron[0]),
    &icosahedron[0]
};

int PPHIGSflag = 0; /* Don't generate PPHIGS format output */
int Flatflag = 0;   /* Don't generate per-vertex normals */
FILE* fh;

/* Forward declarations */
point *normalize(/* point *p */);
point *midpoint(/* point *a, point *b */);
void flip_object(/* object *obj */);
void print_object(/* object *obj, int level */);
void print_triangle(/* triangle *t */);
void pphigs_header(/* int */);
void pphigs_trailer();
void print_script();

main(ac, av)
int ac;
char *av[];
{
    object *old = &oct,		/* Default is octahedron */
	   *new;
    int     ccwflag = 0,	/* Reverse vertex order if true */
	    i,
	    level,		/* Current subdivision level */
	    maxlevel = 1;	/* Maximum subdivision level */

	char fname[255];
	char outf[255];

    /* Parse arguments */
    for (i = 1; i < ac; i++) {
	if (!strcmp(av[i], "-p"))
	    PPHIGSflag = 1;
	else if (!strcmp(av[i], "-c"))
	    ccwflag = 1;
	else if (!strcmp(av[i], "-f"))
	    Flatflag = 1;
	else if (!strcmp(av[i], "-t"))
	    old = &tet;
	else if (!strcmp(av[i], "-i"))
	    old = &ico;
	else if (isdigit(av[i][0])) 
	{
	    if ((maxlevel = atoi(av[i])) < 1) 
		{
			fprintf(stderr, "%s: # of levels must be >= 1\n", av[0]);
			exit(1);
	    }
		
		i++;
		if( (fh = fopen(av[i],"w")) == NULL)
		{
			fprintf(stderr, "Cannot create output file %s\n", av[i]);
			exit(1);
		}

		i++;	
		strcpy(fname, av[i]);

		i++;	
		strcpy(outf, av[i]);
	} 
	else {
	    fprintf(stderr, "Usage: %s [-p] [-c] [-f] [-t] [-i] [n]\n", av[0]);
	    exit(1);
	}
    }

    if (ccwflag)
	flip_object(old);

    /* Subdivide each starting triangle (maxlevel - 1) times */
    for (level = 1; level < maxlevel; 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 triangle 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_script(fname);
    /* Print out resulting approximation */
    print_object(old, maxlevel);
	fprintf(fh, "render framerange:#active outputwidth:640 outputheight:480 superblack:false renderatmosphericeffects:false outputfile: \"%s\" renderer:#draft shadows:false antiAliasing:false autoreflect:false objectmotionblur:false",outf);
}

/* 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 midpoint on the line between 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;
}

/* Reverse order of points in each triangle */
void flip_object(obj)
object *obj;
{
    int i;
    for (i = 0; i < obj->npoly; i++) {
	point tmp;
		       tmp = obj->poly[i].pt[0];
	obj->poly[i].pt[0] = obj->poly[i].pt[2];
	obj->poly[i].pt[2] = tmp;
    }
}

void print_script(name)
char* name;
{
	int i;
	int noq = 1;

	i = strlen(name) - 1;
	if(name[i]=='"')
	{ 
		i--; noq = 0;
	}

	if(name[i]=='s' || name[i]=='S')
	{
		fprintf(fh,"importFile ");
		if(noq) fprintf(fh, "\"");
		fprintf(fh,"%s",name);
		if(noq) fprintf(fh, "\"\n");
	}
	else
	{
		fprintf(fh, "loadMaxFile ", name);
		if(noq) fprintf(fh, "\"");
		fprintf(fh,"%s",name);
		if(noq) fprintf(fh, "\"\n");
	}

	fprintf(fh, "select geometry\ngroup selection name:\"Group01\"\nset coordsys world\n");
	fprintf(fh, "ambientColor = color 0 0 0\nlightTintColor = color 0 0 0\nbackgroundColor = color 255 255 255\n");
	fprintf(fh, "setUseEnvironmentMap false\nselect lights\ndeselect $Group01\ndelete selection\n");
	fprintf(fh, "$Group01.pos=[0,0,0]\ndist=(length($Group01.max - $Group01.min))*2.0\n");
	fprintf(fh, "c = targetCamera pos:[0,0,dist] target:(targetObject pos:[0,0,0])\n");
	fprintf(fh, "viewport.setCamera c\n");
}

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

    if (PPHIGSflag)
	pphigs_header(level);

	fprintf(fh,"animationRange = interval 1 %i\n", obj->npoly);
	fprintf(fh,"c.pos.controller = Position_XYZ ()\n");

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

    if (PPHIGSflag)
	pphigs_trailer();
}

point getmidpoint(triangle* t)
{
	point p;

	p.x = (t->pt[0].x + t->pt[1].x)/2;
	p.y = (t->pt[0].y + t->pt[1].y)/2;
	p.z = (t->pt[0].z + t->pt[1].z)/2;

	p.x = (t->pt[2].x - p.x)/3 + p.x;
	p.y = (t->pt[2].y - p.y)/3 + p.y;
	p.z = (t->pt[2].z - p.z)/3 + p.z;

	return p;
}

/* Output a triangle */
void print_triangle(t,count)
triangle *t;
int count;
{
    int i;
	double alpha;
	double betta;
	double theta;
	double leng;
	point p;

    if (PPHIGSflag) {
	fprintf(fh,"\tpolygon 3 {\n");
	for (i = 0; i < 3; i++)
	    if (Flatflag) {
		fprintf(fh,"\t\t%g\t%g\t%g ;\n",
		    t->pt[i].x, t->pt[i].y, t->pt[i].z);    /* Point */
	    } else {
		fprintf(fh,"\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 */
	    }
	fprintf(fh,"\t};\n");
    } else {
	/* Modify this to generate your favorite output format
	 * Triangle vertices are in t->pt[0..2].{x,y,z}
	 * A generic format is provided.
	 */

	  p = getmidpoint(t);
	  leng = sqrt(p.x*p.x + p.y*p.y + p.z*p.z);
	  alpha = p.x/leng;
	  betta = p.y/leng;
	  theta = p.z/leng;

	  fprintf(fh,"addNewKey c.pos.controller %if\n",count);
	  fprintf(fh,"c.pos.controller.X_Position.controller.keys[%i].value = dist*%0.3f\n",count,alpha);
	  fprintf(fh,"c.pos.controller.Y_Position.controller.keys[%i].value = dist*%0.3f\n",count,betta);
	  fprintf(fh,"c.pos.controller.Z_Position.controller.keys[%i].value = dist*%0.3f\n",count,theta);
	}
}

/* Generate header for a PPHIGS ASCII archive */
void pphigs_header(level)
int level;
{
    int dx, dy, dz;

    if (Flatflag)
	fprintf(fh,"structure flat%d posted {\n", level);
    else
	fprintf(fh,"structure sphere%d posted {\n", level);
    fprintf(fh,"\tcolor polygon {\n");
    fprintf(fh,"\t\t200 100  50   0     50 100 200   0\n");
    fprintf(fh,"\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;
    }

    fprintf(fh,"\tmatrix Pre scale 1000 1000 1000;\n");
    fprintf(fh,"\tmatrix Pre translate %d %d %d ;\n", dx, dy, dz);
}

void pphigs_trailer() {
    fprintf(fh,"};\n");
}
