# HG changeset patch
# User John Tsiombikas <nuclear@member.fsf.org>
# Date 1497452637 -10800
# Node ID 8171de5a000b084e9ba7a40b64f8b3463429deb9

initial commit

diff -r 000000000000 -r 8171de5a000b Makefile
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Makefile	Wed Jun 14 18:03:57 2017 +0300
@@ -0,0 +1,13 @@
+src = $(wildcard src/*.c)
+obj = $(src:.c=.o)
+bin = voxfract
+
+CFLAGS = -pedantic -Wall -g
+LDFLAGS = -lGL -lGLU -lglut -lmetasurf -lm
+
+$(bin): $(obj)
+	$(CC) -o $@ $(obj) $(LDFLAGS)
+
+.PHONY: clean
+clean:
+	rm -f $(obj) $(bin)
diff -r 000000000000 -r 8171de5a000b README
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/README	Wed Jun 14 18:03:57 2017 +0300
@@ -0,0 +1,5 @@
+Voxel 3D fractal generator and polygonizer
+------------------------------------------
+
+Author: John Tsiombikas <nuclear@member.fsf.org>
+This program is public domain. Feel free to use it any way you like.
diff -r 000000000000 -r 8171de5a000b src/main.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/main.c	Wed Jun 14 18:03:57 2017 +0300
@@ -0,0 +1,341 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <assert.h>
+#include <GL/glut.h>
+#include <metasurf.h>
+
+struct quat {
+	float x, y, z, w;
+};
+
+int init(void);
+void cleanup(void);
+void display(void);
+void reshape(int x, int y);
+void keypress(unsigned char key, int x, int y);
+void mouse(int bn, int st, int x, int y);
+void motion(int x, int y);
+
+float eval(struct metasurface *ms, float x, float y, float z);
+void vertex(struct metasurface *ms, float x, float y, float z);
+
+int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter);
+float julia_dist(struct quat *z, struct quat *seed, int max_iter);
+void julia_grad(float *grad, float dist, struct quat *q, struct quat *seed, int max_iter);
+
+void show_waitscr(void);
+
+float cam_theta, cam_phi = 25, cam_dist = 5;
+int grid_size = 350;
+float grid_scale = 1.7;
+struct quat seed = {0.4, 0.0, 0.0, -0.8};
+int max_iter = 9;
+
+struct metasurface *msurf;
+int dlist;
+FILE *fp;
+
+int main(int argc, char **argv)
+{
+	glutInit(&argc, argv);
+	glutInitWindowSize(800, 600);
+	glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE);
+	glutCreateWindow("voxel fractals");
+
+	glutDisplayFunc(display);
+	glutReshapeFunc(reshape);
+	glutKeyboardFunc(keypress);
+	glutMouseFunc(mouse);
+	glutMotionFunc(motion);
+
+	if(init() == -1) {
+		return 1;
+	}
+	atexit(cleanup);
+
+	glutMainLoop();
+	return 0;
+}
+
+int init(void)
+{
+	glEnable(GL_DEPTH_TEST);
+	glEnable(GL_CULL_FACE);
+	glEnable(GL_LIGHTING);
+	glEnable(GL_LIGHT0);
+	/*glShadeModel(GL_FLAT);*/
+
+	fp = fopen("julia.obj", "wb");
+
+	if(!(msurf = msurf_create())) {
+		return -1;
+	}
+	msurf_eval_func(msurf, eval);
+	msurf_vertex_func(msurf, vertex);
+	msurf_set_resolution(msurf, grid_size, grid_size, grid_size);
+	msurf_set_threshold(msurf, 0.0);
+	msurf_set_inside(msurf, MSURF_LESS);
+
+	show_waitscr();
+
+	dlist = glGenLists(1);
+	glNewList(dlist, GL_COMPILE);
+
+	glBegin(GL_TRIANGLES);
+	msurf_polygonize(msurf);
+	glEnd();
+
+	glEndList();
+
+	if(fp) {
+		fclose(fp);
+	}
+
+	glClearColor(0.05, 0.05, 0.05, 1);
+
+	return 0;
+}
+
+void cleanup(void)
+{
+	msurf_free(msurf);
+}
+
+void display(void)
+{
+	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+	glMatrixMode(GL_MODELVIEW);
+	glLoadIdentity();
+	glTranslatef(0, 0, -cam_dist);
+	glRotatef(cam_phi, 1, 0, 0);
+	glRotatef(cam_theta, 0, 1, 0);
+
+	glCallList(dlist);
+
+	assert(glGetError() == GL_NO_ERROR);
+	glutSwapBuffers();
+}
+
+void reshape(int x, int y)
+{
+	glViewport(0, 0, x, y);
+	glMatrixMode(GL_PROJECTION);
+	glLoadIdentity();
+	gluPerspective(50.0, (float)x / (float)y, 0.5, 500.0);
+}
+
+void keypress(unsigned char key, int x, int y)
+{
+	switch(key) {
+	case 27:
+		exit(0);
+	}
+}
+
+int prev_x, prev_y;
+int bnstate[8];
+
+void mouse(int bn, int st, int x, int y)
+{
+	bnstate[bn - GLUT_LEFT_BUTTON] = st == GLUT_DOWN;
+	prev_x = x;
+	prev_y = y;
+}
+
+void motion(int x, int y)
+{
+	int dx = x - prev_x;
+	int dy = y - prev_y;
+	prev_x = x;
+	prev_y = y;
+
+	if(!dx && !dy) return;
+
+	if(bnstate[0]) {
+		cam_theta += dx * 0.5;
+		cam_phi += dy * 0.5;
+
+		if(cam_phi < -90) cam_phi = -90;
+		if(cam_phi > 90) cam_phi = 90;
+		glutPostRedisplay();
+	}
+	if(bnstate[2]) {
+		cam_dist += dy * 0.1;
+
+		if(cam_dist < 0.0) cam_dist = 0.0;
+		glutPostRedisplay();
+	}
+}
+
+float eval(struct metasurface *ms, float x, float y, float z)
+{
+	struct quat q;
+
+	q.w = grid_scale * x;
+	q.x = grid_scale * y;
+	q.y = grid_scale * z;
+	q.z = 0.0f;
+
+	return julia_dist(&q, &seed, max_iter);
+}
+
+void vertex(struct metasurface *ms, float x, float y, float z)
+{
+	struct quat q;
+	float dist;
+	float norm[3];
+	float norm_len, nscale = 1.0;
+
+	q.w = grid_scale * x;
+	q.x = grid_scale * y;
+	q.y = grid_scale * z;
+	q.z = 0.0f;
+
+	dist = julia_dist(&q, &seed, max_iter);
+
+	julia_grad(norm, dist, &q, &seed, max_iter);
+	norm_len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
+
+	if(norm_len != 0.0f) {
+		nscale = 1.0f / norm_len;
+	}
+
+	glNormal3f(norm[0] * nscale, norm[1] * nscale, norm[2] * nscale);
+	glVertex3f(x, y, z);
+
+	if(fp) {
+		static int nverts;
+
+		fprintf(fp, "v %f %f %f\n", x, y, z);
+
+		if((++nverts) % 3 == 0) {
+			fprintf(fp, "f %d %d %d\n", nverts - 2, nverts - 1, nverts);
+		}
+	}
+}
+
+void quat_mul(struct quat *res, struct quat *a, struct quat *b)
+{
+	float w = a->w * b->w - (a->x * b->x + a->y * b->y + a->z * b->z);
+	float x = a->w * b->x + b->w * a->x + (a->y * b->z - a->z * b->y);
+	float y = a->w * b->y + b->w * a->y + (a->z * b->x - a->x * b->z);
+	res->z = a->w * b->z + b->w * a->z + (a->x * b->y - a->y * b->x);
+	res->x = x;
+	res->y = y;
+	res->w = w;
+}
+
+void quat_sq(struct quat *q)
+{
+	float w = q->w * q->w - (q->x * q->x + q->y * q->y + q->z * q->z);
+	float x = 2.0 * q->w * q->x;
+	float y = 2.0 * q->w * q->y;
+	q->z = 2.0 * q->w * q->z;
+	q->x = x;
+	q->y = y;
+	q->w = w;
+}
+
+float quat_lensq(struct quat *q)
+{
+	return q->x * q->x + q->y * q->y + q->z * q->z + q->w * q->w;
+}
+
+int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter)
+{
+	int i;
+
+	*qres = *q;
+	qprime->x = qprime->y = qprime->z = 0.0f;
+	qprime->w = 1.0f;
+
+	for(i=0; i<max_iter; i++) {
+		quat_mul(qprime, qres, qprime);
+		qprime->x *= 2.0;
+		qprime->y *= 2.0;
+		qprime->z *= 2.0;
+		qprime->w *= 2.0;
+
+		quat_sq(qres);
+		qres->x += seed->x;
+		qres->y += seed->y;
+		qres->z += seed->z;
+		qres->w += seed->w;
+
+		if(quat_lensq(qres) > 8.0) {
+			return 0;
+		}
+	}
+
+	return 1;
+}
+
+float julia_dist(struct quat *z, struct quat *seed, int max_iter)
+{
+	struct quat qprime, q;
+
+	/* calc julia at z */
+	/*int inside = */julia(&q, &qprime, z, seed, max_iter);
+
+	float lenq = sqrt(quat_lensq(&q));
+	float lenqp = sqrt(quat_lensq(&qprime));
+	return 0.5 * lenq * log(lenq) / lenqp;
+}
+
+#define OFFS 1e-4
+void julia_grad(float *grad, float dist, struct quat *q, struct quat *seed, int max_iter)
+{
+	struct quat qnext = *q;
+	struct quat qprev = *q;
+
+	qnext.w += OFFS;
+	qprev.w -= OFFS;
+	grad[0] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
+	qnext.w = qprev.w = q->w;
+
+	qnext.x += OFFS;
+	qprev.x -= OFFS;
+	grad[1] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
+	qnext.x = qprev.x = q->x;
+
+	qnext.y += OFFS;
+	qprev.y -= OFFS;
+	grad[2] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
+}
+
+void show_waitscr(void)
+{
+	const char *text = "Please wait, generating fractal...";
+	glClear(GL_COLOR_BUFFER_BIT);
+
+	glPushAttrib(GL_ENABLE_BIT);
+	glDisable(GL_LIGHTING);
+	glDisable(GL_DEPTH_TEST);
+
+	glMatrixMode(GL_PROJECTION);
+	glPushMatrix();
+	glLoadIdentity();
+	glTranslatef(-0.75, 0, 0);
+	glScalef(0.00075, 0.00075, 0.00075);
+	glMatrixMode(GL_MODELVIEW);
+
+	glEnable(GL_LINE_SMOOTH);
+	glEnable(GL_BLEND);
+	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+
+	glLineWidth(1.5);
+
+	glColor3f(1, 1, 1);
+	while(*text) {
+		glutStrokeCharacter(GLUT_STROKE_ROMAN, *text++);
+	}
+
+	glMatrixMode(GL_PROJECTION);
+	glPopMatrix();
+
+	glPopAttrib();
+
+	glutSwapBuffers();
+}