1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/main.c	Wed Jun 14 18:03:57 2017 +0300
     1.3 @@ -0,0 +1,341 @@
     1.4 +#include <stdio.h>
     1.5 +#include <stdlib.h>
     1.6 +#include <math.h>
     1.7 +#include <assert.h>
     1.8 +#include <GL/glut.h>
     1.9 +#include <metasurf.h>
    1.10 +
    1.11 +struct quat {
    1.12 +	float x, y, z, w;
    1.13 +};
    1.14 +
    1.15 +int init(void);
    1.16 +void cleanup(void);
    1.17 +void display(void);
    1.18 +void reshape(int x, int y);
    1.19 +void keypress(unsigned char key, int x, int y);
    1.20 +void mouse(int bn, int st, int x, int y);
    1.21 +void motion(int x, int y);
    1.22 +
    1.23 +float eval(struct metasurface *ms, float x, float y, float z);
    1.24 +void vertex(struct metasurface *ms, float x, float y, float z);
    1.25 +
    1.26 +int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter);
    1.27 +float julia_dist(struct quat *z, struct quat *seed, int max_iter);
    1.28 +void julia_grad(float *grad, float dist, struct quat *q, struct quat *seed, int max_iter);
    1.29 +
    1.30 +void show_waitscr(void);
    1.31 +
    1.32 +float cam_theta, cam_phi = 25, cam_dist = 5;
    1.33 +int grid_size = 350;
    1.34 +float grid_scale = 1.7;
    1.35 +struct quat seed = {0.4, 0.0, 0.0, -0.8};
    1.36 +int max_iter = 9;
    1.37 +
    1.38 +struct metasurface *msurf;
    1.39 +int dlist;
    1.40 +FILE *fp;
    1.41 +
    1.42 +int main(int argc, char **argv)
    1.43 +{
    1.44 +	glutInit(&argc, argv);
    1.45 +	glutInitWindowSize(800, 600);
    1.46 +	glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE);
    1.47 +	glutCreateWindow("voxel fractals");
    1.48 +
    1.49 +	glutDisplayFunc(display);
    1.50 +	glutReshapeFunc(reshape);
    1.51 +	glutKeyboardFunc(keypress);
    1.52 +	glutMouseFunc(mouse);
    1.53 +	glutMotionFunc(motion);
    1.54 +
    1.55 +	if(init() == -1) {
    1.56 +		return 1;
    1.57 +	}
    1.58 +	atexit(cleanup);
    1.59 +
    1.60 +	glutMainLoop();
    1.61 +	return 0;
    1.62 +}
    1.63 +
    1.64 +int init(void)
    1.65 +{
    1.66 +	glEnable(GL_DEPTH_TEST);
    1.67 +	glEnable(GL_CULL_FACE);
    1.68 +	glEnable(GL_LIGHTING);
    1.69 +	glEnable(GL_LIGHT0);
    1.70 +	/*glShadeModel(GL_FLAT);*/
    1.71 +
    1.72 +	fp = fopen("julia.obj", "wb");
    1.73 +
    1.74 +	if(!(msurf = msurf_create())) {
    1.75 +		return -1;
    1.76 +	}
    1.77 +	msurf_eval_func(msurf, eval);
    1.78 +	msurf_vertex_func(msurf, vertex);
    1.79 +	msurf_set_resolution(msurf, grid_size, grid_size, grid_size);
    1.80 +	msurf_set_threshold(msurf, 0.0);
    1.81 +	msurf_set_inside(msurf, MSURF_LESS);
    1.82 +
    1.83 +	show_waitscr();
    1.84 +
    1.85 +	dlist = glGenLists(1);
    1.86 +	glNewList(dlist, GL_COMPILE);
    1.87 +
    1.88 +	glBegin(GL_TRIANGLES);
    1.89 +	msurf_polygonize(msurf);
    1.90 +	glEnd();
    1.91 +
    1.92 +	glEndList();
    1.93 +
    1.94 +	if(fp) {
    1.95 +		fclose(fp);
    1.96 +	}
    1.97 +
    1.98 +	glClearColor(0.05, 0.05, 0.05, 1);
    1.99 +
   1.100 +	return 0;
   1.101 +}
   1.102 +
   1.103 +void cleanup(void)
   1.104 +{
   1.105 +	msurf_free(msurf);
   1.106 +}
   1.107 +
   1.108 +void display(void)
   1.109 +{
   1.110 +	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
   1.111 +
   1.112 +	glMatrixMode(GL_MODELVIEW);
   1.113 +	glLoadIdentity();
   1.114 +	glTranslatef(0, 0, -cam_dist);
   1.115 +	glRotatef(cam_phi, 1, 0, 0);
   1.116 +	glRotatef(cam_theta, 0, 1, 0);
   1.117 +
   1.118 +	glCallList(dlist);
   1.119 +
   1.120 +	assert(glGetError() == GL_NO_ERROR);
   1.121 +	glutSwapBuffers();
   1.122 +}
   1.123 +
   1.124 +void reshape(int x, int y)
   1.125 +{
   1.126 +	glViewport(0, 0, x, y);
   1.127 +	glMatrixMode(GL_PROJECTION);
   1.128 +	glLoadIdentity();
   1.129 +	gluPerspective(50.0, (float)x / (float)y, 0.5, 500.0);
   1.130 +}
   1.131 +
   1.132 +void keypress(unsigned char key, int x, int y)
   1.133 +{
   1.134 +	switch(key) {
   1.135 +	case 27:
   1.136 +		exit(0);
   1.137 +	}
   1.138 +}
   1.139 +
   1.140 +int prev_x, prev_y;
   1.141 +int bnstate[8];
   1.142 +
   1.143 +void mouse(int bn, int st, int x, int y)
   1.144 +{
   1.145 +	bnstate[bn - GLUT_LEFT_BUTTON] = st == GLUT_DOWN;
   1.146 +	prev_x = x;
   1.147 +	prev_y = y;
   1.148 +}
   1.149 +
   1.150 +void motion(int x, int y)
   1.151 +{
   1.152 +	int dx = x - prev_x;
   1.153 +	int dy = y - prev_y;
   1.154 +	prev_x = x;
   1.155 +	prev_y = y;
   1.156 +
   1.157 +	if(!dx && !dy) return;
   1.158 +
   1.159 +	if(bnstate[0]) {
   1.160 +		cam_theta += dx * 0.5;
   1.161 +		cam_phi += dy * 0.5;
   1.162 +
   1.163 +		if(cam_phi < -90) cam_phi = -90;
   1.164 +		if(cam_phi > 90) cam_phi = 90;
   1.165 +		glutPostRedisplay();
   1.166 +	}
   1.167 +	if(bnstate[2]) {
   1.168 +		cam_dist += dy * 0.1;
   1.169 +
   1.170 +		if(cam_dist < 0.0) cam_dist = 0.0;
   1.171 +		glutPostRedisplay();
   1.172 +	}
   1.173 +}
   1.174 +
   1.175 +float eval(struct metasurface *ms, float x, float y, float z)
   1.176 +{
   1.177 +	struct quat q;
   1.178 +
   1.179 +	q.w = grid_scale * x;
   1.180 +	q.x = grid_scale * y;
   1.181 +	q.y = grid_scale * z;
   1.182 +	q.z = 0.0f;
   1.183 +
   1.184 +	return julia_dist(&q, &seed, max_iter);
   1.185 +}
   1.186 +
   1.187 +void vertex(struct metasurface *ms, float x, float y, float z)
   1.188 +{
   1.189 +	struct quat q;
   1.190 +	float dist;
   1.191 +	float norm[3];
   1.192 +	float norm_len, nscale = 1.0;
   1.193 +
   1.194 +	q.w = grid_scale * x;
   1.195 +	q.x = grid_scale * y;
   1.196 +	q.y = grid_scale * z;
   1.197 +	q.z = 0.0f;
   1.198 +
   1.199 +	dist = julia_dist(&q, &seed, max_iter);
   1.200 +
   1.201 +	julia_grad(norm, dist, &q, &seed, max_iter);
   1.202 +	norm_len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
   1.203 +
   1.204 +	if(norm_len != 0.0f) {
   1.205 +		nscale = 1.0f / norm_len;
   1.206 +	}
   1.207 +
   1.208 +	glNormal3f(norm[0] * nscale, norm[1] * nscale, norm[2] * nscale);
   1.209 +	glVertex3f(x, y, z);
   1.210 +
   1.211 +	if(fp) {
   1.212 +		static int nverts;
   1.213 +
   1.214 +		fprintf(fp, "v %f %f %f\n", x, y, z);
   1.215 +
   1.216 +		if((++nverts) % 3 == 0) {
   1.217 +			fprintf(fp, "f %d %d %d\n", nverts - 2, nverts - 1, nverts);
   1.218 +		}
   1.219 +	}
   1.220 +}
   1.221 +
   1.222 +void quat_mul(struct quat *res, struct quat *a, struct quat *b)
   1.223 +{
   1.224 +	float w = a->w * b->w - (a->x * b->x + a->y * b->y + a->z * b->z);
   1.225 +	float x = a->w * b->x + b->w * a->x + (a->y * b->z - a->z * b->y);
   1.226 +	float y = a->w * b->y + b->w * a->y + (a->z * b->x - a->x * b->z);
   1.227 +	res->z = a->w * b->z + b->w * a->z + (a->x * b->y - a->y * b->x);
   1.228 +	res->x = x;
   1.229 +	res->y = y;
   1.230 +	res->w = w;
   1.231 +}
   1.232 +
   1.233 +void quat_sq(struct quat *q)
   1.234 +{
   1.235 +	float w = q->w * q->w - (q->x * q->x + q->y * q->y + q->z * q->z);
   1.236 +	float x = 2.0 * q->w * q->x;
   1.237 +	float y = 2.0 * q->w * q->y;
   1.238 +	q->z = 2.0 * q->w * q->z;
   1.239 +	q->x = x;
   1.240 +	q->y = y;
   1.241 +	q->w = w;
   1.242 +}
   1.243 +
   1.244 +float quat_lensq(struct quat *q)
   1.245 +{
   1.246 +	return q->x * q->x + q->y * q->y + q->z * q->z + q->w * q->w;
   1.247 +}
   1.248 +
   1.249 +int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter)
   1.250 +{
   1.251 +	int i;
   1.252 +
   1.253 +	*qres = *q;
   1.254 +	qprime->x = qprime->y = qprime->z = 0.0f;
   1.255 +	qprime->w = 1.0f;
   1.256 +
   1.257 +	for(i=0; i<max_iter; i++) {
   1.258 +		quat_mul(qprime, qres, qprime);
   1.259 +		qprime->x *= 2.0;
   1.260 +		qprime->y *= 2.0;
   1.261 +		qprime->z *= 2.0;
   1.262 +		qprime->w *= 2.0;
   1.263 +
   1.264 +		quat_sq(qres);
   1.265 +		qres->x += seed->x;
   1.266 +		qres->y += seed->y;
   1.267 +		qres->z += seed->z;
   1.268 +		qres->w += seed->w;
   1.269 +
   1.270 +		if(quat_lensq(qres) > 8.0) {
   1.271 +			return 0;
   1.272 +		}
   1.273 +	}
   1.274 +
   1.275 +	return 1;
   1.276 +}
   1.277 +
   1.278 +float julia_dist(struct quat *z, struct quat *seed, int max_iter)
   1.279 +{
   1.280 +	struct quat qprime, q;
   1.281 +
   1.282 +	/* calc julia at z */
   1.283 +	/*int inside = */julia(&q, &qprime, z, seed, max_iter);
   1.284 +
   1.285 +	float lenq = sqrt(quat_lensq(&q));
   1.286 +	float lenqp = sqrt(quat_lensq(&qprime));
   1.287 +	return 0.5 * lenq * log(lenq) / lenqp;
   1.288 +}
   1.289 +
   1.290 +#define OFFS 1e-4
   1.291 +void julia_grad(float *grad, float dist, struct quat *q, struct quat *seed, int max_iter)
   1.292 +{
   1.293 +	struct quat qnext = *q;
   1.294 +	struct quat qprev = *q;
   1.295 +
   1.296 +	qnext.w += OFFS;
   1.297 +	qprev.w -= OFFS;
   1.298 +	grad[0] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
   1.299 +	qnext.w = qprev.w = q->w;
   1.300 +
   1.301 +	qnext.x += OFFS;
   1.302 +	qprev.x -= OFFS;
   1.303 +	grad[1] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
   1.304 +	qnext.x = qprev.x = q->x;
   1.305 +
   1.306 +	qnext.y += OFFS;
   1.307 +	qprev.y -= OFFS;
   1.308 +	grad[2] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
   1.309 +}
   1.310 +
   1.311 +void show_waitscr(void)
   1.312 +{
   1.313 +	const char *text = "Please wait, generating fractal...";
   1.314 +	glClear(GL_COLOR_BUFFER_BIT);
   1.315 +
   1.316 +	glPushAttrib(GL_ENABLE_BIT);
   1.317 +	glDisable(GL_LIGHTING);
   1.318 +	glDisable(GL_DEPTH_TEST);
   1.319 +
   1.320 +	glMatrixMode(GL_PROJECTION);
   1.321 +	glPushMatrix();
   1.322 +	glLoadIdentity();
   1.323 +	glTranslatef(-0.75, 0, 0);
   1.324 +	glScalef(0.00075, 0.00075, 0.00075);
   1.325 +	glMatrixMode(GL_MODELVIEW);
   1.326 +
   1.327 +	glEnable(GL_LINE_SMOOTH);
   1.328 +	glEnable(GL_BLEND);
   1.329 +	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
   1.330 +
   1.331 +	glLineWidth(1.5);
   1.332 +
   1.333 +	glColor3f(1, 1, 1);
   1.334 +	while(*text) {
   1.335 +		glutStrokeCharacter(GLUT_STROKE_ROMAN, *text++);
   1.336 +	}
   1.337 +
   1.338 +	glMatrixMode(GL_PROJECTION);
   1.339 +	glPopMatrix();
   1.340 +
   1.341 +	glPopAttrib();
   1.342 +
   1.343 +	glutSwapBuffers();
   1.344 +}