1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/Makefile	Wed Jun 14 18:03:57 2017 +0300
     1.3 @@ -0,0 +1,13 @@
     1.4 +src = $(wildcard src/*.c)
     1.5 +obj = $(src:.c=.o)
     1.6 +bin = voxfract
     1.7 +
     1.8 +CFLAGS = -pedantic -Wall -g
     1.9 +LDFLAGS = -lGL -lGLU -lglut -lmetasurf -lm
    1.10 +
    1.11 +$(bin): $(obj)
    1.12 +	$(CC) -o $@ $(obj) $(LDFLAGS)
    1.13 +
    1.14 +.PHONY: clean
    1.15 +clean:
    1.16 +	rm -f $(obj) $(bin)

     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/README	Wed Jun 14 18:03:57 2017 +0300
     2.3 @@ -0,0 +1,5 @@
     2.4 +Voxel 3D fractal generator and polygonizer
     2.5 +------------------------------------------
     2.6 +
     2.7 +Author: John Tsiombikas <nuclear@member.fsf.org>
     2.8 +This program is public domain. Feel free to use it any way you like.

     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/src/main.c	Wed Jun 14 18:03:57 2017 +0300
     3.3 @@ -0,0 +1,341 @@
     3.4 +#include <stdio.h>
     3.5 +#include <stdlib.h>
     3.6 +#include <math.h>
     3.7 +#include <assert.h>
     3.8 +#include <GL/glut.h>
     3.9 +#include <metasurf.h>
    3.10 +
    3.11 +struct quat {
    3.12 +	float x, y, z, w;
    3.13 +};
    3.14 +
    3.15 +int init(void);
    3.16 +void cleanup(void);
    3.17 +void display(void);
    3.18 +void reshape(int x, int y);
    3.19 +void keypress(unsigned char key, int x, int y);
    3.20 +void mouse(int bn, int st, int x, int y);
    3.21 +void motion(int x, int y);
    3.22 +
    3.23 +float eval(struct metasurface *ms, float x, float y, float z);
    3.24 +void vertex(struct metasurface *ms, float x, float y, float z);
    3.25 +
    3.26 +int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter);
    3.27 +float julia_dist(struct quat *z, struct quat *seed, int max_iter);
    3.28 +void julia_grad(float *grad, float dist, struct quat *q, struct quat *seed, int max_iter);
    3.29 +
    3.30 +void show_waitscr(void);
    3.31 +
    3.32 +float cam_theta, cam_phi = 25, cam_dist = 5;
    3.33 +int grid_size = 350;
    3.34 +float grid_scale = 1.7;
    3.35 +struct quat seed = {0.4, 0.0, 0.0, -0.8};
    3.36 +int max_iter = 9;
    3.37 +
    3.38 +struct metasurface *msurf;
    3.39 +int dlist;
    3.40 +FILE *fp;
    3.41 +
    3.42 +int main(int argc, char **argv)
    3.43 +{
    3.44 +	glutInit(&argc, argv);
    3.45 +	glutInitWindowSize(800, 600);
    3.46 +	glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE);
    3.47 +	glutCreateWindow("voxel fractals");
    3.48 +
    3.49 +	glutDisplayFunc(display);
    3.50 +	glutReshapeFunc(reshape);
    3.51 +	glutKeyboardFunc(keypress);
    3.52 +	glutMouseFunc(mouse);
    3.53 +	glutMotionFunc(motion);
    3.54 +
    3.55 +	if(init() == -1) {
    3.56 +		return 1;
    3.57 +	}
    3.58 +	atexit(cleanup);
    3.59 +
    3.60 +	glutMainLoop();
    3.61 +	return 0;
    3.62 +}
    3.63 +
    3.64 +int init(void)
    3.65 +{
    3.66 +	glEnable(GL_DEPTH_TEST);
    3.67 +	glEnable(GL_CULL_FACE);
    3.68 +	glEnable(GL_LIGHTING);
    3.69 +	glEnable(GL_LIGHT0);
    3.70 +	/*glShadeModel(GL_FLAT);*/
    3.71 +
    3.72 +	fp = fopen("julia.obj", "wb");
    3.73 +
    3.74 +	if(!(msurf = msurf_create())) {
    3.75 +		return -1;
    3.76 +	}
    3.77 +	msurf_eval_func(msurf, eval);
    3.78 +	msurf_vertex_func(msurf, vertex);
    3.79 +	msurf_set_resolution(msurf, grid_size, grid_size, grid_size);
    3.80 +	msurf_set_threshold(msurf, 0.0);
    3.81 +	msurf_set_inside(msurf, MSURF_LESS);
    3.82 +
    3.83 +	show_waitscr();
    3.84 +
    3.85 +	dlist = glGenLists(1);
    3.86 +	glNewList(dlist, GL_COMPILE);
    3.87 +
    3.88 +	glBegin(GL_TRIANGLES);
    3.89 +	msurf_polygonize(msurf);
    3.90 +	glEnd();
    3.91 +
    3.92 +	glEndList();
    3.93 +
    3.94 +	if(fp) {
    3.95 +		fclose(fp);
    3.96 +	}
    3.97 +
    3.98 +	glClearColor(0.05, 0.05, 0.05, 1);
    3.99 +
   3.100 +	return 0;
   3.101 +}
   3.102 +
   3.103 +void cleanup(void)
   3.104 +{
   3.105 +	msurf_free(msurf);
   3.106 +}
   3.107 +
   3.108 +void display(void)
   3.109 +{
   3.110 +	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
   3.111 +
   3.112 +	glMatrixMode(GL_MODELVIEW);
   3.113 +	glLoadIdentity();
   3.114 +	glTranslatef(0, 0, -cam_dist);
   3.115 +	glRotatef(cam_phi, 1, 0, 0);
   3.116 +	glRotatef(cam_theta, 0, 1, 0);
   3.117 +
   3.118 +	glCallList(dlist);
   3.119 +
   3.120 +	assert(glGetError() == GL_NO_ERROR);
   3.121 +	glutSwapBuffers();
   3.122 +}
   3.123 +
   3.124 +void reshape(int x, int y)
   3.125 +{
   3.126 +	glViewport(0, 0, x, y);
   3.127 +	glMatrixMode(GL_PROJECTION);
   3.128 +	glLoadIdentity();
   3.129 +	gluPerspective(50.0, (float)x / (float)y, 0.5, 500.0);
   3.130 +}
   3.131 +
   3.132 +void keypress(unsigned char key, int x, int y)
   3.133 +{
   3.134 +	switch(key) {
   3.135 +	case 27:
   3.136 +		exit(0);
   3.137 +	}
   3.138 +}
   3.139 +
   3.140 +int prev_x, prev_y;
   3.141 +int bnstate[8];
   3.142 +
   3.143 +void mouse(int bn, int st, int x, int y)
   3.144 +{
   3.145 +	bnstate[bn - GLUT_LEFT_BUTTON] = st == GLUT_DOWN;
   3.146 +	prev_x = x;
   3.147 +	prev_y = y;
   3.148 +}
   3.149 +
   3.150 +void motion(int x, int y)
   3.151 +{
   3.152 +	int dx = x - prev_x;
   3.153 +	int dy = y - prev_y;
   3.154 +	prev_x = x;
   3.155 +	prev_y = y;
   3.156 +
   3.157 +	if(!dx && !dy) return;
   3.158 +
   3.159 +	if(bnstate[0]) {
   3.160 +		cam_theta += dx * 0.5;
   3.161 +		cam_phi += dy * 0.5;
   3.162 +
   3.163 +		if(cam_phi < -90) cam_phi = -90;
   3.164 +		if(cam_phi > 90) cam_phi = 90;
   3.165 +		glutPostRedisplay();
   3.166 +	}
   3.167 +	if(bnstate[2]) {
   3.168 +		cam_dist += dy * 0.1;
   3.169 +
   3.170 +		if(cam_dist < 0.0) cam_dist = 0.0;
   3.171 +		glutPostRedisplay();
   3.172 +	}
   3.173 +}
   3.174 +
   3.175 +float eval(struct metasurface *ms, float x, float y, float z)
   3.176 +{
   3.177 +	struct quat q;
   3.178 +
   3.179 +	q.w = grid_scale * x;
   3.180 +	q.x = grid_scale * y;
   3.181 +	q.y = grid_scale * z;
   3.182 +	q.z = 0.0f;
   3.183 +
   3.184 +	return julia_dist(&q, &seed, max_iter);
   3.185 +}
   3.186 +
   3.187 +void vertex(struct metasurface *ms, float x, float y, float z)
   3.188 +{
   3.189 +	struct quat q;
   3.190 +	float dist;
   3.191 +	float norm[3];
   3.192 +	float norm_len, nscale = 1.0;
   3.193 +
   3.194 +	q.w = grid_scale * x;
   3.195 +	q.x = grid_scale * y;
   3.196 +	q.y = grid_scale * z;
   3.197 +	q.z = 0.0f;
   3.198 +
   3.199 +	dist = julia_dist(&q, &seed, max_iter);
   3.200 +
   3.201 +	julia_grad(norm, dist, &q, &seed, max_iter);
   3.202 +	norm_len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
   3.203 +
   3.204 +	if(norm_len != 0.0f) {
   3.205 +		nscale = 1.0f / norm_len;
   3.206 +	}
   3.207 +
   3.208 +	glNormal3f(norm[0] * nscale, norm[1] * nscale, norm[2] * nscale);
   3.209 +	glVertex3f(x, y, z);
   3.210 +
   3.211 +	if(fp) {
   3.212 +		static int nverts;
   3.213 +
   3.214 +		fprintf(fp, "v %f %f %f\n", x, y, z);
   3.215 +
   3.216 +		if((++nverts) % 3 == 0) {
   3.217 +			fprintf(fp, "f %d %d %d\n", nverts - 2, nverts - 1, nverts);
   3.218 +		}
   3.219 +	}
   3.220 +}
   3.221 +
   3.222 +void quat_mul(struct quat *res, struct quat *a, struct quat *b)
   3.223 +{
   3.224 +	float w = a->w * b->w - (a->x * b->x + a->y * b->y + a->z * b->z);
   3.225 +	float x = a->w * b->x + b->w * a->x + (a->y * b->z - a->z * b->y);
   3.226 +	float y = a->w * b->y + b->w * a->y + (a->z * b->x - a->x * b->z);
   3.227 +	res->z = a->w * b->z + b->w * a->z + (a->x * b->y - a->y * b->x);
   3.228 +	res->x = x;
   3.229 +	res->y = y;
   3.230 +	res->w = w;
   3.231 +}
   3.232 +
   3.233 +void quat_sq(struct quat *q)
   3.234 +{
   3.235 +	float w = q->w * q->w - (q->x * q->x + q->y * q->y + q->z * q->z);
   3.236 +	float x = 2.0 * q->w * q->x;
   3.237 +	float y = 2.0 * q->w * q->y;
   3.238 +	q->z = 2.0 * q->w * q->z;
   3.239 +	q->x = x;
   3.240 +	q->y = y;
   3.241 +	q->w = w;
   3.242 +}
   3.243 +
   3.244 +float quat_lensq(struct quat *q)
   3.245 +{
   3.246 +	return q->x * q->x + q->y * q->y + q->z * q->z + q->w * q->w;
   3.247 +}
   3.248 +
   3.249 +int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter)
   3.250 +{
   3.251 +	int i;
   3.252 +
   3.253 +	*qres = *q;
   3.254 +	qprime->x = qprime->y = qprime->z = 0.0f;
   3.255 +	qprime->w = 1.0f;
   3.256 +
   3.257 +	for(i=0; i<max_iter; i++) {
   3.258 +		quat_mul(qprime, qres, qprime);
   3.259 +		qprime->x *= 2.0;
   3.260 +		qprime->y *= 2.0;
   3.261 +		qprime->z *= 2.0;
   3.262 +		qprime->w *= 2.0;
   3.263 +
   3.264 +		quat_sq(qres);
   3.265 +		qres->x += seed->x;
   3.266 +		qres->y += seed->y;
   3.267 +		qres->z += seed->z;
   3.268 +		qres->w += seed->w;
   3.269 +
   3.270 +		if(quat_lensq(qres) > 8.0) {
   3.271 +			return 0;
   3.272 +		}
   3.273 +	}
   3.274 +
   3.275 +	return 1;
   3.276 +}
   3.277 +
   3.278 +float julia_dist(struct quat *z, struct quat *seed, int max_iter)
   3.279 +{
   3.280 +	struct quat qprime, q;
   3.281 +
   3.282 +	/* calc julia at z */
   3.283 +	/*int inside = */julia(&q, &qprime, z, seed, max_iter);
   3.284 +
   3.285 +	float lenq = sqrt(quat_lensq(&q));
   3.286 +	float lenqp = sqrt(quat_lensq(&qprime));
   3.287 +	return 0.5 * lenq * log(lenq) / lenqp;
   3.288 +}
   3.289 +
   3.290 +#define OFFS 1e-4
   3.291 +void julia_grad(float *grad, float dist, struct quat *q, struct quat *seed, int max_iter)
   3.292 +{
   3.293 +	struct quat qnext = *q;
   3.294 +	struct quat qprev = *q;
   3.295 +
   3.296 +	qnext.w += OFFS;
   3.297 +	qprev.w -= OFFS;
   3.298 +	grad[0] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
   3.299 +	qnext.w = qprev.w = q->w;
   3.300 +
   3.301 +	qnext.x += OFFS;
   3.302 +	qprev.x -= OFFS;
   3.303 +	grad[1] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
   3.304 +	qnext.x = qprev.x = q->x;
   3.305 +
   3.306 +	qnext.y += OFFS;
   3.307 +	qprev.y -= OFFS;
   3.308 +	grad[2] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
   3.309 +}
   3.310 +
   3.311 +void show_waitscr(void)
   3.312 +{
   3.313 +	const char *text = "Please wait, generating fractal...";
   3.314 +	glClear(GL_COLOR_BUFFER_BIT);
   3.315 +
   3.316 +	glPushAttrib(GL_ENABLE_BIT);
   3.317 +	glDisable(GL_LIGHTING);
   3.318 +	glDisable(GL_DEPTH_TEST);
   3.319 +
   3.320 +	glMatrixMode(GL_PROJECTION);
   3.321 +	glPushMatrix();
   3.322 +	glLoadIdentity();
   3.323 +	glTranslatef(-0.75, 0, 0);
   3.324 +	glScalef(0.00075, 0.00075, 0.00075);
   3.325 +	glMatrixMode(GL_MODELVIEW);
   3.326 +
   3.327 +	glEnable(GL_LINE_SMOOTH);
   3.328 +	glEnable(GL_BLEND);
   3.329 +	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
   3.330 +
   3.331 +	glLineWidth(1.5);
   3.332 +
   3.333 +	glColor3f(1, 1, 1);
   3.334 +	while(*text) {
   3.335 +		glutStrokeCharacter(GLUT_STROKE_ROMAN, *text++);
   3.336 +	}
   3.337 +
   3.338 +	glMatrixMode(GL_PROJECTION);
   3.339 +	glPopMatrix();
   3.340 +
   3.341 +	glPopAttrib();
   3.342 +
   3.343 +	glutSwapBuffers();
   3.344 +}