nuclear@0: #include <stdio.h>
nuclear@0: #include <stdlib.h>
nuclear@0: #include <math.h>
nuclear@0: #include <assert.h>
nuclear@0: #include <GL/glut.h>
nuclear@0: #include <metasurf.h>
nuclear@0: 
nuclear@0: struct quat {
nuclear@0: 	float x, y, z, w;
nuclear@0: };
nuclear@0: 
nuclear@0: int init(void);
nuclear@0: void cleanup(void);
nuclear@0: void display(void);
nuclear@0: void reshape(int x, int y);
nuclear@0: void keypress(unsigned char key, int x, int y);
nuclear@0: void mouse(int bn, int st, int x, int y);
nuclear@0: void motion(int x, int y);
nuclear@0: 
nuclear@0: float eval(struct metasurface *ms, float x, float y, float z);
nuclear@0: void vertex(struct metasurface *ms, float x, float y, float z);
nuclear@0: 
nuclear@0: int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter);
nuclear@0: float julia_dist(struct quat *z, struct quat *seed, int max_iter);
nuclear@1: void julia_grad(float *grad, struct quat *q, struct quat *seed, int max_iter);
nuclear@0: 
nuclear@0: void show_waitscr(void);
nuclear@0: 
nuclear@0: float cam_theta, cam_phi = 25, cam_dist = 5;
nuclear@1: int grid_size = 300;
nuclear@0: float grid_scale = 1.7;
nuclear@0: struct quat seed = {0.4, 0.0, 0.0, -0.8};
nuclear@1: int max_iter = 10;
nuclear@0: 
nuclear@0: struct metasurface *msurf;
nuclear@0: int dlist;
nuclear@0: FILE *fp;
nuclear@0: 
nuclear@0: int main(int argc, char **argv)
nuclear@0: {
nuclear@1: 	int i;
nuclear@1: 
nuclear@1: 	for(i=1; i<argc; i++) {
nuclear@1: 		if(argv[i][0] == '-') {
nuclear@1: 			if(argv[i][2] != 0) {
nuclear@1: 				fprintf(stderr, "invalid option: %s\n", argv[i]);
nuclear@1: 				return 1;
nuclear@1: 			}
nuclear@1: 			switch(argv[i][1]) {
nuclear@1: 			case 'g':
nuclear@1: 				if((grid_size = atoi(argv[++i])) <= 0) {
nuclear@1: 					fprintf(stderr, "invalid grid size value: %s\n", argv[i]);
nuclear@1: 					return 1;
nuclear@1: 				}
nuclear@1: 				break;
nuclear@1: 
nuclear@1: 			case 'i':
nuclear@1: 				if((max_iter = atoi(argv[++i])) <= 0) {
nuclear@1: 					fprintf(stderr, "invalid iterations count: %s\n", argv[i]);
nuclear@1: 					return 1;
nuclear@1: 				}
nuclear@1: 				break;
nuclear@1: 
nuclear@1: 			default:
nuclear@1: 				fprintf(stderr, "invalid option: %s\n", argv[i]);
nuclear@1: 				return 1;
nuclear@1: 			}
nuclear@1: 		} else {
nuclear@1: 			fprintf(stderr, "unexpected argument: %s\n", argv[i]);
nuclear@1: 			return 1;
nuclear@1: 		}
nuclear@1: 	}
nuclear@1: 	printf("grid size: %d\n", grid_size);
nuclear@1: 	printf("max iter: %d\n", max_iter);
nuclear@1: 
nuclear@0: 	glutInit(&argc, argv);
nuclear@0: 	glutInitWindowSize(800, 600);
nuclear@0: 	glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE);
nuclear@0: 	glutCreateWindow("voxel fractals");
nuclear@0: 
nuclear@0: 	glutDisplayFunc(display);
nuclear@0: 	glutReshapeFunc(reshape);
nuclear@0: 	glutKeyboardFunc(keypress);
nuclear@0: 	glutMouseFunc(mouse);
nuclear@0: 	glutMotionFunc(motion);
nuclear@0: 
nuclear@0: 	if(init() == -1) {
nuclear@0: 		return 1;
nuclear@0: 	}
nuclear@0: 	atexit(cleanup);
nuclear@0: 
nuclear@0: 	glutMainLoop();
nuclear@0: 	return 0;
nuclear@0: }
nuclear@0: 
nuclear@0: int init(void)
nuclear@0: {
nuclear@0: 	glEnable(GL_DEPTH_TEST);
nuclear@0: 	glEnable(GL_CULL_FACE);
nuclear@0: 	glEnable(GL_LIGHTING);
nuclear@0: 	glEnable(GL_LIGHT0);
nuclear@0: 	/*glShadeModel(GL_FLAT);*/
nuclear@0: 
nuclear@0: 	fp = fopen("julia.obj", "wb");
nuclear@0: 
nuclear@0: 	if(!(msurf = msurf_create())) {
nuclear@0: 		return -1;
nuclear@0: 	}
nuclear@0: 	msurf_eval_func(msurf, eval);
nuclear@0: 	msurf_vertex_func(msurf, vertex);
nuclear@0: 	msurf_set_resolution(msurf, grid_size, grid_size, grid_size);
nuclear@0: 	msurf_set_threshold(msurf, 0.0);
nuclear@0: 	msurf_set_inside(msurf, MSURF_LESS);
nuclear@0: 
nuclear@0: 	show_waitscr();
nuclear@0: 
nuclear@0: 	dlist = glGenLists(1);
nuclear@0: 	glNewList(dlist, GL_COMPILE);
nuclear@0: 
nuclear@0: 	glBegin(GL_TRIANGLES);
nuclear@0: 	msurf_polygonize(msurf);
nuclear@0: 	glEnd();
nuclear@0: 
nuclear@0: 	glEndList();
nuclear@0: 
nuclear@0: 	if(fp) {
nuclear@0: 		fclose(fp);
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	glClearColor(0.05, 0.05, 0.05, 1);
nuclear@0: 
nuclear@0: 	return 0;
nuclear@0: }
nuclear@0: 
nuclear@0: void cleanup(void)
nuclear@0: {
nuclear@0: 	msurf_free(msurf);
nuclear@0: }
nuclear@0: 
nuclear@0: void display(void)
nuclear@0: {
nuclear@0: 	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
nuclear@0: 
nuclear@0: 	glMatrixMode(GL_MODELVIEW);
nuclear@0: 	glLoadIdentity();
nuclear@0: 	glTranslatef(0, 0, -cam_dist);
nuclear@0: 	glRotatef(cam_phi, 1, 0, 0);
nuclear@0: 	glRotatef(cam_theta, 0, 1, 0);
nuclear@0: 
nuclear@0: 	glCallList(dlist);
nuclear@0: 
nuclear@0: 	assert(glGetError() == GL_NO_ERROR);
nuclear@0: 	glutSwapBuffers();
nuclear@0: }
nuclear@0: 
nuclear@0: void reshape(int x, int y)
nuclear@0: {
nuclear@0: 	glViewport(0, 0, x, y);
nuclear@0: 	glMatrixMode(GL_PROJECTION);
nuclear@0: 	glLoadIdentity();
nuclear@0: 	gluPerspective(50.0, (float)x / (float)y, 0.5, 500.0);
nuclear@0: }
nuclear@0: 
nuclear@0: void keypress(unsigned char key, int x, int y)
nuclear@0: {
nuclear@1: 	static int flat;
nuclear@1: 
nuclear@0: 	switch(key) {
nuclear@0: 	case 27:
nuclear@0: 		exit(0);
nuclear@1: 
nuclear@1: 	case 's':
nuclear@1: 		flat = !flat;
nuclear@1: 		if(flat) {
nuclear@1: 			glShadeModel(GL_FLAT);
nuclear@1: 		} else {
nuclear@1: 			glShadeModel(GL_SMOOTH);
nuclear@1: 		}
nuclear@1: 		glutPostRedisplay();
nuclear@0: 	}
nuclear@0: }
nuclear@0: 
nuclear@0: int prev_x, prev_y;
nuclear@0: int bnstate[8];
nuclear@0: 
nuclear@0: void mouse(int bn, int st, int x, int y)
nuclear@0: {
nuclear@0: 	bnstate[bn - GLUT_LEFT_BUTTON] = st == GLUT_DOWN;
nuclear@0: 	prev_x = x;
nuclear@0: 	prev_y = y;
nuclear@0: }
nuclear@0: 
nuclear@0: void motion(int x, int y)
nuclear@0: {
nuclear@0: 	int dx = x - prev_x;
nuclear@0: 	int dy = y - prev_y;
nuclear@0: 	prev_x = x;
nuclear@0: 	prev_y = y;
nuclear@0: 
nuclear@0: 	if(!dx && !dy) return;
nuclear@0: 
nuclear@0: 	if(bnstate[0]) {
nuclear@0: 		cam_theta += dx * 0.5;
nuclear@0: 		cam_phi += dy * 0.5;
nuclear@0: 
nuclear@0: 		if(cam_phi < -90) cam_phi = -90;
nuclear@0: 		if(cam_phi > 90) cam_phi = 90;
nuclear@0: 		glutPostRedisplay();
nuclear@0: 	}
nuclear@0: 	if(bnstate[2]) {
nuclear@0: 		cam_dist += dy * 0.1;
nuclear@0: 
nuclear@0: 		if(cam_dist < 0.0) cam_dist = 0.0;
nuclear@0: 		glutPostRedisplay();
nuclear@0: 	}
nuclear@0: }
nuclear@0: 
nuclear@0: float eval(struct metasurface *ms, float x, float y, float z)
nuclear@0: {
nuclear@0: 	struct quat q;
nuclear@0: 
nuclear@0: 	q.w = grid_scale * x;
nuclear@0: 	q.x = grid_scale * y;
nuclear@0: 	q.y = grid_scale * z;
nuclear@0: 	q.z = 0.0f;
nuclear@0: 
nuclear@0: 	return julia_dist(&q, &seed, max_iter);
nuclear@0: }
nuclear@0: 
nuclear@0: void vertex(struct metasurface *ms, float x, float y, float z)
nuclear@0: {
nuclear@0: 	struct quat q;
nuclear@0: 	float norm[3];
nuclear@0: 	float norm_len, nscale = 1.0;
nuclear@0: 
nuclear@0: 	q.w = grid_scale * x;
nuclear@0: 	q.x = grid_scale * y;
nuclear@0: 	q.y = grid_scale * z;
nuclear@0: 	q.z = 0.0f;
nuclear@0: 
nuclear@1: 	julia_grad(norm, &q, &seed, max_iter);
nuclear@0: 	norm_len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
nuclear@0: 
nuclear@1: 	if(norm_len != 0.0) {
nuclear@1: 		nscale = 1.0 / norm_len;
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	glNormal3f(norm[0] * nscale, norm[1] * nscale, norm[2] * nscale);
nuclear@0: 	glVertex3f(x, y, z);
nuclear@0: 
nuclear@0: 	if(fp) {
nuclear@0: 		static int nverts;
nuclear@0: 
nuclear@0: 		fprintf(fp, "v %f %f %f\n", x, y, z);
nuclear@0: 
nuclear@0: 		if((++nverts) % 3 == 0) {
nuclear@0: 			fprintf(fp, "f %d %d %d\n", nverts - 2, nverts - 1, nverts);
nuclear@0: 		}
nuclear@0: 	}
nuclear@0: }
nuclear@0: 
nuclear@0: void quat_mul(struct quat *res, struct quat *a, struct quat *b)
nuclear@0: {
nuclear@0: 	float w = a->w * b->w - (a->x * b->x + a->y * b->y + a->z * b->z);
nuclear@0: 	float x = a->w * b->x + b->w * a->x + (a->y * b->z - a->z * b->y);
nuclear@0: 	float y = a->w * b->y + b->w * a->y + (a->z * b->x - a->x * b->z);
nuclear@0: 	res->z = a->w * b->z + b->w * a->z + (a->x * b->y - a->y * b->x);
nuclear@0: 	res->x = x;
nuclear@0: 	res->y = y;
nuclear@0: 	res->w = w;
nuclear@0: }
nuclear@0: 
nuclear@0: void quat_sq(struct quat *q)
nuclear@0: {
nuclear@0: 	float w = q->w * q->w - (q->x * q->x + q->y * q->y + q->z * q->z);
nuclear@0: 	float x = 2.0 * q->w * q->x;
nuclear@0: 	float y = 2.0 * q->w * q->y;
nuclear@0: 	q->z = 2.0 * q->w * q->z;
nuclear@0: 	q->x = x;
nuclear@0: 	q->y = y;
nuclear@0: 	q->w = w;
nuclear@0: }
nuclear@0: 
nuclear@0: float quat_lensq(struct quat *q)
nuclear@0: {
nuclear@0: 	return q->x * q->x + q->y * q->y + q->z * q->z + q->w * q->w;
nuclear@0: }
nuclear@0: 
nuclear@0: int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter)
nuclear@0: {
nuclear@0: 	int i;
nuclear@0: 
nuclear@0: 	*qres = *q;
nuclear@0: 	qprime->x = qprime->y = qprime->z = 0.0f;
nuclear@0: 	qprime->w = 1.0f;
nuclear@0: 
nuclear@0: 	for(i=0; i<max_iter; i++) {
nuclear@0: 		quat_mul(qprime, qres, qprime);
nuclear@0: 		qprime->x *= 2.0;
nuclear@0: 		qprime->y *= 2.0;
nuclear@0: 		qprime->z *= 2.0;
nuclear@0: 		qprime->w *= 2.0;
nuclear@0: 
nuclear@0: 		quat_sq(qres);
nuclear@0: 		qres->x += seed->x;
nuclear@0: 		qres->y += seed->y;
nuclear@0: 		qres->z += seed->z;
nuclear@0: 		qres->w += seed->w;
nuclear@0: 
nuclear@0: 		if(quat_lensq(qres) > 8.0) {
nuclear@0: 			return 0;
nuclear@0: 		}
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	return 1;
nuclear@0: }
nuclear@0: 
nuclear@0: float julia_dist(struct quat *z, struct quat *seed, int max_iter)
nuclear@0: {
nuclear@0: 	struct quat qprime, q;
nuclear@0: 
nuclear@0: 	/* calc julia at z */
nuclear@0: 	/*int inside = */julia(&q, &qprime, z, seed, max_iter);
nuclear@0: 
nuclear@0: 	float lenq = sqrt(quat_lensq(&q));
nuclear@0: 	float lenqp = sqrt(quat_lensq(&qprime));
nuclear@0: 	return 0.5 * lenq * log(lenq) / lenqp;
nuclear@0: }
nuclear@0: 
nuclear@1: #define OFFS 1e-5
nuclear@1: void julia_grad(float *grad, struct quat *q, struct quat *seed, int max_iter)
nuclear@0: {
nuclear@0: 	struct quat qnext = *q;
nuclear@0: 	struct quat qprev = *q;
nuclear@0: 
nuclear@0: 	qnext.w += OFFS;
nuclear@0: 	qprev.w -= OFFS;
nuclear@0: 	grad[0] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
nuclear@0: 	qnext.w = qprev.w = q->w;
nuclear@0: 
nuclear@0: 	qnext.x += OFFS;
nuclear@0: 	qprev.x -= OFFS;
nuclear@0: 	grad[1] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
nuclear@0: 	qnext.x = qprev.x = q->x;
nuclear@0: 
nuclear@0: 	qnext.y += OFFS;
nuclear@0: 	qprev.y -= OFFS;
nuclear@0: 	grad[2] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
nuclear@0: }
nuclear@0: 
nuclear@0: void show_waitscr(void)
nuclear@0: {
nuclear@0: 	const char *text = "Please wait, generating fractal...";
nuclear@0: 	glClear(GL_COLOR_BUFFER_BIT);
nuclear@0: 
nuclear@0: 	glPushAttrib(GL_ENABLE_BIT);
nuclear@0: 	glDisable(GL_LIGHTING);
nuclear@0: 	glDisable(GL_DEPTH_TEST);
nuclear@0: 
nuclear@0: 	glMatrixMode(GL_PROJECTION);
nuclear@0: 	glPushMatrix();
nuclear@0: 	glLoadIdentity();
nuclear@0: 	glTranslatef(-0.75, 0, 0);
nuclear@0: 	glScalef(0.00075, 0.00075, 0.00075);
nuclear@0: 	glMatrixMode(GL_MODELVIEW);
nuclear@0: 
nuclear@0: 	glEnable(GL_LINE_SMOOTH);
nuclear@0: 	glEnable(GL_BLEND);
nuclear@0: 	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
nuclear@0: 
nuclear@0: 	glLineWidth(1.5);
nuclear@0: 
nuclear@0: 	glColor3f(1, 1, 1);
nuclear@0: 	while(*text) {
nuclear@0: 		glutStrokeCharacter(GLUT_STROKE_ROMAN, *text++);
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	glMatrixMode(GL_PROJECTION);
nuclear@0: 	glPopMatrix();
nuclear@0: 
nuclear@0: 	glPopAttrib();
nuclear@0: 
nuclear@0: 	glutSwapBuffers();
nuclear@0: }