rev |
line source |
nuclear@0
|
1 #include <stdio.h>
|
nuclear@0
|
2 #include <stdlib.h>
|
nuclear@0
|
3 #include <math.h>
|
nuclear@0
|
4 #include <assert.h>
|
nuclear@0
|
5 #include <GL/glut.h>
|
nuclear@0
|
6 #include <metasurf.h>
|
nuclear@0
|
7
|
nuclear@0
|
8 struct quat {
|
nuclear@0
|
9 float x, y, z, w;
|
nuclear@0
|
10 };
|
nuclear@0
|
11
|
nuclear@0
|
12 int init(void);
|
nuclear@0
|
13 void cleanup(void);
|
nuclear@0
|
14 void display(void);
|
nuclear@0
|
15 void reshape(int x, int y);
|
nuclear@0
|
16 void keypress(unsigned char key, int x, int y);
|
nuclear@0
|
17 void mouse(int bn, int st, int x, int y);
|
nuclear@0
|
18 void motion(int x, int y);
|
nuclear@0
|
19
|
nuclear@0
|
20 float eval(struct metasurface *ms, float x, float y, float z);
|
nuclear@0
|
21 void vertex(struct metasurface *ms, float x, float y, float z);
|
nuclear@0
|
22
|
nuclear@0
|
23 int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter);
|
nuclear@0
|
24 float julia_dist(struct quat *z, struct quat *seed, int max_iter);
|
nuclear@1
|
25 void julia_grad(float *grad, struct quat *q, struct quat *seed, int max_iter);
|
nuclear@0
|
26
|
nuclear@0
|
27 void show_waitscr(void);
|
nuclear@0
|
28
|
nuclear@0
|
29 float cam_theta, cam_phi = 25, cam_dist = 5;
|
nuclear@1
|
30 int grid_size = 300;
|
nuclear@0
|
31 float grid_scale = 1.7;
|
nuclear@0
|
32 struct quat seed = {0.4, 0.0, 0.0, -0.8};
|
nuclear@1
|
33 int max_iter = 10;
|
nuclear@0
|
34
|
nuclear@0
|
35 struct metasurface *msurf;
|
nuclear@0
|
36 int dlist;
|
nuclear@0
|
37 FILE *fp;
|
nuclear@0
|
38
|
nuclear@0
|
39 int main(int argc, char **argv)
|
nuclear@0
|
40 {
|
nuclear@1
|
41 int i;
|
nuclear@1
|
42
|
nuclear@1
|
43 for(i=1; i<argc; i++) {
|
nuclear@1
|
44 if(argv[i][0] == '-') {
|
nuclear@1
|
45 if(argv[i][2] != 0) {
|
nuclear@1
|
46 fprintf(stderr, "invalid option: %s\n", argv[i]);
|
nuclear@1
|
47 return 1;
|
nuclear@1
|
48 }
|
nuclear@1
|
49 switch(argv[i][1]) {
|
nuclear@1
|
50 case 'g':
|
nuclear@1
|
51 if((grid_size = atoi(argv[++i])) <= 0) {
|
nuclear@1
|
52 fprintf(stderr, "invalid grid size value: %s\n", argv[i]);
|
nuclear@1
|
53 return 1;
|
nuclear@1
|
54 }
|
nuclear@1
|
55 break;
|
nuclear@1
|
56
|
nuclear@1
|
57 case 'i':
|
nuclear@1
|
58 if((max_iter = atoi(argv[++i])) <= 0) {
|
nuclear@1
|
59 fprintf(stderr, "invalid iterations count: %s\n", argv[i]);
|
nuclear@1
|
60 return 1;
|
nuclear@1
|
61 }
|
nuclear@1
|
62 break;
|
nuclear@1
|
63
|
nuclear@1
|
64 default:
|
nuclear@1
|
65 fprintf(stderr, "invalid option: %s\n", argv[i]);
|
nuclear@1
|
66 return 1;
|
nuclear@1
|
67 }
|
nuclear@1
|
68 } else {
|
nuclear@1
|
69 fprintf(stderr, "unexpected argument: %s\n", argv[i]);
|
nuclear@1
|
70 return 1;
|
nuclear@1
|
71 }
|
nuclear@1
|
72 }
|
nuclear@1
|
73 printf("grid size: %d\n", grid_size);
|
nuclear@1
|
74 printf("max iter: %d\n", max_iter);
|
nuclear@1
|
75
|
nuclear@0
|
76 glutInit(&argc, argv);
|
nuclear@0
|
77 glutInitWindowSize(800, 600);
|
nuclear@0
|
78 glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE);
|
nuclear@0
|
79 glutCreateWindow("voxel fractals");
|
nuclear@0
|
80
|
nuclear@0
|
81 glutDisplayFunc(display);
|
nuclear@0
|
82 glutReshapeFunc(reshape);
|
nuclear@0
|
83 glutKeyboardFunc(keypress);
|
nuclear@0
|
84 glutMouseFunc(mouse);
|
nuclear@0
|
85 glutMotionFunc(motion);
|
nuclear@0
|
86
|
nuclear@0
|
87 if(init() == -1) {
|
nuclear@0
|
88 return 1;
|
nuclear@0
|
89 }
|
nuclear@0
|
90 atexit(cleanup);
|
nuclear@0
|
91
|
nuclear@0
|
92 glutMainLoop();
|
nuclear@0
|
93 return 0;
|
nuclear@0
|
94 }
|
nuclear@0
|
95
|
nuclear@0
|
96 int init(void)
|
nuclear@0
|
97 {
|
nuclear@0
|
98 glEnable(GL_DEPTH_TEST);
|
nuclear@0
|
99 glEnable(GL_CULL_FACE);
|
nuclear@0
|
100 glEnable(GL_LIGHTING);
|
nuclear@0
|
101 glEnable(GL_LIGHT0);
|
nuclear@0
|
102 /*glShadeModel(GL_FLAT);*/
|
nuclear@0
|
103
|
nuclear@0
|
104 fp = fopen("julia.obj", "wb");
|
nuclear@0
|
105
|
nuclear@0
|
106 if(!(msurf = msurf_create())) {
|
nuclear@0
|
107 return -1;
|
nuclear@0
|
108 }
|
nuclear@0
|
109 msurf_eval_func(msurf, eval);
|
nuclear@0
|
110 msurf_vertex_func(msurf, vertex);
|
nuclear@0
|
111 msurf_set_resolution(msurf, grid_size, grid_size, grid_size);
|
nuclear@0
|
112 msurf_set_threshold(msurf, 0.0);
|
nuclear@0
|
113 msurf_set_inside(msurf, MSURF_LESS);
|
nuclear@0
|
114
|
nuclear@0
|
115 show_waitscr();
|
nuclear@0
|
116
|
nuclear@0
|
117 dlist = glGenLists(1);
|
nuclear@0
|
118 glNewList(dlist, GL_COMPILE);
|
nuclear@0
|
119
|
nuclear@0
|
120 glBegin(GL_TRIANGLES);
|
nuclear@0
|
121 msurf_polygonize(msurf);
|
nuclear@0
|
122 glEnd();
|
nuclear@0
|
123
|
nuclear@0
|
124 glEndList();
|
nuclear@0
|
125
|
nuclear@0
|
126 if(fp) {
|
nuclear@0
|
127 fclose(fp);
|
nuclear@0
|
128 }
|
nuclear@0
|
129
|
nuclear@0
|
130 glClearColor(0.05, 0.05, 0.05, 1);
|
nuclear@0
|
131
|
nuclear@0
|
132 return 0;
|
nuclear@0
|
133 }
|
nuclear@0
|
134
|
nuclear@0
|
135 void cleanup(void)
|
nuclear@0
|
136 {
|
nuclear@0
|
137 msurf_free(msurf);
|
nuclear@0
|
138 }
|
nuclear@0
|
139
|
nuclear@0
|
140 void display(void)
|
nuclear@0
|
141 {
|
nuclear@0
|
142 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
|
nuclear@0
|
143
|
nuclear@0
|
144 glMatrixMode(GL_MODELVIEW);
|
nuclear@0
|
145 glLoadIdentity();
|
nuclear@0
|
146 glTranslatef(0, 0, -cam_dist);
|
nuclear@0
|
147 glRotatef(cam_phi, 1, 0, 0);
|
nuclear@0
|
148 glRotatef(cam_theta, 0, 1, 0);
|
nuclear@0
|
149
|
nuclear@0
|
150 glCallList(dlist);
|
nuclear@0
|
151
|
nuclear@0
|
152 assert(glGetError() == GL_NO_ERROR);
|
nuclear@0
|
153 glutSwapBuffers();
|
nuclear@0
|
154 }
|
nuclear@0
|
155
|
nuclear@0
|
156 void reshape(int x, int y)
|
nuclear@0
|
157 {
|
nuclear@0
|
158 glViewport(0, 0, x, y);
|
nuclear@0
|
159 glMatrixMode(GL_PROJECTION);
|
nuclear@0
|
160 glLoadIdentity();
|
nuclear@0
|
161 gluPerspective(50.0, (float)x / (float)y, 0.5, 500.0);
|
nuclear@0
|
162 }
|
nuclear@0
|
163
|
nuclear@0
|
164 void keypress(unsigned char key, int x, int y)
|
nuclear@0
|
165 {
|
nuclear@1
|
166 static int flat;
|
nuclear@1
|
167
|
nuclear@0
|
168 switch(key) {
|
nuclear@0
|
169 case 27:
|
nuclear@0
|
170 exit(0);
|
nuclear@1
|
171
|
nuclear@1
|
172 case 's':
|
nuclear@1
|
173 flat = !flat;
|
nuclear@1
|
174 if(flat) {
|
nuclear@1
|
175 glShadeModel(GL_FLAT);
|
nuclear@1
|
176 } else {
|
nuclear@1
|
177 glShadeModel(GL_SMOOTH);
|
nuclear@1
|
178 }
|
nuclear@1
|
179 glutPostRedisplay();
|
nuclear@0
|
180 }
|
nuclear@0
|
181 }
|
nuclear@0
|
182
|
nuclear@0
|
183 int prev_x, prev_y;
|
nuclear@0
|
184 int bnstate[8];
|
nuclear@0
|
185
|
nuclear@0
|
186 void mouse(int bn, int st, int x, int y)
|
nuclear@0
|
187 {
|
nuclear@0
|
188 bnstate[bn - GLUT_LEFT_BUTTON] = st == GLUT_DOWN;
|
nuclear@0
|
189 prev_x = x;
|
nuclear@0
|
190 prev_y = y;
|
nuclear@0
|
191 }
|
nuclear@0
|
192
|
nuclear@0
|
193 void motion(int x, int y)
|
nuclear@0
|
194 {
|
nuclear@0
|
195 int dx = x - prev_x;
|
nuclear@0
|
196 int dy = y - prev_y;
|
nuclear@0
|
197 prev_x = x;
|
nuclear@0
|
198 prev_y = y;
|
nuclear@0
|
199
|
nuclear@0
|
200 if(!dx && !dy) return;
|
nuclear@0
|
201
|
nuclear@0
|
202 if(bnstate[0]) {
|
nuclear@0
|
203 cam_theta += dx * 0.5;
|
nuclear@0
|
204 cam_phi += dy * 0.5;
|
nuclear@0
|
205
|
nuclear@0
|
206 if(cam_phi < -90) cam_phi = -90;
|
nuclear@0
|
207 if(cam_phi > 90) cam_phi = 90;
|
nuclear@0
|
208 glutPostRedisplay();
|
nuclear@0
|
209 }
|
nuclear@0
|
210 if(bnstate[2]) {
|
nuclear@0
|
211 cam_dist += dy * 0.1;
|
nuclear@0
|
212
|
nuclear@0
|
213 if(cam_dist < 0.0) cam_dist = 0.0;
|
nuclear@0
|
214 glutPostRedisplay();
|
nuclear@0
|
215 }
|
nuclear@0
|
216 }
|
nuclear@0
|
217
|
nuclear@0
|
218 float eval(struct metasurface *ms, float x, float y, float z)
|
nuclear@0
|
219 {
|
nuclear@0
|
220 struct quat q;
|
nuclear@0
|
221
|
nuclear@0
|
222 q.w = grid_scale * x;
|
nuclear@0
|
223 q.x = grid_scale * y;
|
nuclear@0
|
224 q.y = grid_scale * z;
|
nuclear@0
|
225 q.z = 0.0f;
|
nuclear@0
|
226
|
nuclear@0
|
227 return julia_dist(&q, &seed, max_iter);
|
nuclear@0
|
228 }
|
nuclear@0
|
229
|
nuclear@0
|
230 void vertex(struct metasurface *ms, float x, float y, float z)
|
nuclear@0
|
231 {
|
nuclear@0
|
232 struct quat q;
|
nuclear@0
|
233 float norm[3];
|
nuclear@0
|
234 float norm_len, nscale = 1.0;
|
nuclear@0
|
235
|
nuclear@0
|
236 q.w = grid_scale * x;
|
nuclear@0
|
237 q.x = grid_scale * y;
|
nuclear@0
|
238 q.y = grid_scale * z;
|
nuclear@0
|
239 q.z = 0.0f;
|
nuclear@0
|
240
|
nuclear@1
|
241 julia_grad(norm, &q, &seed, max_iter);
|
nuclear@0
|
242 norm_len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
|
nuclear@0
|
243
|
nuclear@1
|
244 if(norm_len != 0.0) {
|
nuclear@1
|
245 nscale = 1.0 / norm_len;
|
nuclear@0
|
246 }
|
nuclear@0
|
247
|
nuclear@0
|
248 glNormal3f(norm[0] * nscale, norm[1] * nscale, norm[2] * nscale);
|
nuclear@0
|
249 glVertex3f(x, y, z);
|
nuclear@0
|
250
|
nuclear@0
|
251 if(fp) {
|
nuclear@0
|
252 static int nverts;
|
nuclear@0
|
253
|
nuclear@0
|
254 fprintf(fp, "v %f %f %f\n", x, y, z);
|
nuclear@0
|
255
|
nuclear@0
|
256 if((++nverts) % 3 == 0) {
|
nuclear@0
|
257 fprintf(fp, "f %d %d %d\n", nverts - 2, nverts - 1, nverts);
|
nuclear@0
|
258 }
|
nuclear@0
|
259 }
|
nuclear@0
|
260 }
|
nuclear@0
|
261
|
nuclear@0
|
262 void quat_mul(struct quat *res, struct quat *a, struct quat *b)
|
nuclear@0
|
263 {
|
nuclear@0
|
264 float w = a->w * b->w - (a->x * b->x + a->y * b->y + a->z * b->z);
|
nuclear@0
|
265 float x = a->w * b->x + b->w * a->x + (a->y * b->z - a->z * b->y);
|
nuclear@0
|
266 float y = a->w * b->y + b->w * a->y + (a->z * b->x - a->x * b->z);
|
nuclear@0
|
267 res->z = a->w * b->z + b->w * a->z + (a->x * b->y - a->y * b->x);
|
nuclear@0
|
268 res->x = x;
|
nuclear@0
|
269 res->y = y;
|
nuclear@0
|
270 res->w = w;
|
nuclear@0
|
271 }
|
nuclear@0
|
272
|
nuclear@0
|
273 void quat_sq(struct quat *q)
|
nuclear@0
|
274 {
|
nuclear@0
|
275 float w = q->w * q->w - (q->x * q->x + q->y * q->y + q->z * q->z);
|
nuclear@0
|
276 float x = 2.0 * q->w * q->x;
|
nuclear@0
|
277 float y = 2.0 * q->w * q->y;
|
nuclear@0
|
278 q->z = 2.0 * q->w * q->z;
|
nuclear@0
|
279 q->x = x;
|
nuclear@0
|
280 q->y = y;
|
nuclear@0
|
281 q->w = w;
|
nuclear@0
|
282 }
|
nuclear@0
|
283
|
nuclear@0
|
284 float quat_lensq(struct quat *q)
|
nuclear@0
|
285 {
|
nuclear@0
|
286 return q->x * q->x + q->y * q->y + q->z * q->z + q->w * q->w;
|
nuclear@0
|
287 }
|
nuclear@0
|
288
|
nuclear@0
|
289 int julia(struct quat *qres, struct quat *qprime, struct quat *q, struct quat *seed, int max_iter)
|
nuclear@0
|
290 {
|
nuclear@0
|
291 int i;
|
nuclear@0
|
292
|
nuclear@0
|
293 *qres = *q;
|
nuclear@0
|
294 qprime->x = qprime->y = qprime->z = 0.0f;
|
nuclear@0
|
295 qprime->w = 1.0f;
|
nuclear@0
|
296
|
nuclear@0
|
297 for(i=0; i<max_iter; i++) {
|
nuclear@0
|
298 quat_mul(qprime, qres, qprime);
|
nuclear@0
|
299 qprime->x *= 2.0;
|
nuclear@0
|
300 qprime->y *= 2.0;
|
nuclear@0
|
301 qprime->z *= 2.0;
|
nuclear@0
|
302 qprime->w *= 2.0;
|
nuclear@0
|
303
|
nuclear@0
|
304 quat_sq(qres);
|
nuclear@0
|
305 qres->x += seed->x;
|
nuclear@0
|
306 qres->y += seed->y;
|
nuclear@0
|
307 qres->z += seed->z;
|
nuclear@0
|
308 qres->w += seed->w;
|
nuclear@0
|
309
|
nuclear@0
|
310 if(quat_lensq(qres) > 8.0) {
|
nuclear@0
|
311 return 0;
|
nuclear@0
|
312 }
|
nuclear@0
|
313 }
|
nuclear@0
|
314
|
nuclear@0
|
315 return 1;
|
nuclear@0
|
316 }
|
nuclear@0
|
317
|
nuclear@0
|
318 float julia_dist(struct quat *z, struct quat *seed, int max_iter)
|
nuclear@0
|
319 {
|
nuclear@0
|
320 struct quat qprime, q;
|
nuclear@0
|
321
|
nuclear@0
|
322 /* calc julia at z */
|
nuclear@0
|
323 /*int inside = */julia(&q, &qprime, z, seed, max_iter);
|
nuclear@0
|
324
|
nuclear@0
|
325 float lenq = sqrt(quat_lensq(&q));
|
nuclear@0
|
326 float lenqp = sqrt(quat_lensq(&qprime));
|
nuclear@0
|
327 return 0.5 * lenq * log(lenq) / lenqp;
|
nuclear@0
|
328 }
|
nuclear@0
|
329
|
nuclear@1
|
330 #define OFFS 1e-5
|
nuclear@1
|
331 void julia_grad(float *grad, struct quat *q, struct quat *seed, int max_iter)
|
nuclear@0
|
332 {
|
nuclear@0
|
333 struct quat qnext = *q;
|
nuclear@0
|
334 struct quat qprev = *q;
|
nuclear@0
|
335
|
nuclear@0
|
336 qnext.w += OFFS;
|
nuclear@0
|
337 qprev.w -= OFFS;
|
nuclear@0
|
338 grad[0] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
|
nuclear@0
|
339 qnext.w = qprev.w = q->w;
|
nuclear@0
|
340
|
nuclear@0
|
341 qnext.x += OFFS;
|
nuclear@0
|
342 qprev.x -= OFFS;
|
nuclear@0
|
343 grad[1] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
|
nuclear@0
|
344 qnext.x = qprev.x = q->x;
|
nuclear@0
|
345
|
nuclear@0
|
346 qnext.y += OFFS;
|
nuclear@0
|
347 qprev.y -= OFFS;
|
nuclear@0
|
348 grad[2] = julia_dist(&qnext, seed, max_iter) - julia_dist(&qprev, seed, max_iter);
|
nuclear@0
|
349 }
|
nuclear@0
|
350
|
nuclear@0
|
351 void show_waitscr(void)
|
nuclear@0
|
352 {
|
nuclear@0
|
353 const char *text = "Please wait, generating fractal...";
|
nuclear@0
|
354 glClear(GL_COLOR_BUFFER_BIT);
|
nuclear@0
|
355
|
nuclear@0
|
356 glPushAttrib(GL_ENABLE_BIT);
|
nuclear@0
|
357 glDisable(GL_LIGHTING);
|
nuclear@0
|
358 glDisable(GL_DEPTH_TEST);
|
nuclear@0
|
359
|
nuclear@0
|
360 glMatrixMode(GL_PROJECTION);
|
nuclear@0
|
361 glPushMatrix();
|
nuclear@0
|
362 glLoadIdentity();
|
nuclear@0
|
363 glTranslatef(-0.75, 0, 0);
|
nuclear@0
|
364 glScalef(0.00075, 0.00075, 0.00075);
|
nuclear@0
|
365 glMatrixMode(GL_MODELVIEW);
|
nuclear@0
|
366
|
nuclear@0
|
367 glEnable(GL_LINE_SMOOTH);
|
nuclear@0
|
368 glEnable(GL_BLEND);
|
nuclear@0
|
369 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
|
nuclear@0
|
370
|
nuclear@0
|
371 glLineWidth(1.5);
|
nuclear@0
|
372
|
nuclear@0
|
373 glColor3f(1, 1, 1);
|
nuclear@0
|
374 while(*text) {
|
nuclear@0
|
375 glutStrokeCharacter(GLUT_STROKE_ROMAN, *text++);
|
nuclear@0
|
376 }
|
nuclear@0
|
377
|
nuclear@0
|
378 glMatrixMode(GL_PROJECTION);
|
nuclear@0
|
379 glPopMatrix();
|
nuclear@0
|
380
|
nuclear@0
|
381 glPopAttrib();
|
nuclear@0
|
382
|
nuclear@0
|
383 glutSwapBuffers();
|
nuclear@0
|
384 }
|