10 return rand() / (RAND_MAX + 1.0);
13 typedef void (*mapfunc_t) (double x_in, double y_in, double *x_out, double *y_out, double *z_out);
14 typedef void (*colorfunc_t) (double x, double y, double z, double *r, double *g, double *b);
16 void color_test(double x, double y, double z, double *r, double *g, double *b)
18 // put in a nice function here
24 double mandelbrot_iter(double zx, double zy, double cx, double cy, int maxiter)
33 for(i = 1; i < maxiter; ++i)
36 zx = zx * zx - zy * zy + cx;
37 zy = 2 * tmp * zy + cy;
39 f = zx * zx + zy * zy;
48 // f: the greater, the more in 0 direction
49 // the smaller, the more in 1 direction
51 // the greater, the more in 0 direction
52 return i + 1 / (f - 4 + 1); // f = 16: + 0, f = 4: + 1
54 // i.e. 0 for i=1, 1 for i=maxiter
57 double mandelbrot_range(double zx, double zy, double cx, double cy, int maxiter, double offset)
59 double i = mandelbrot_iter(zx, zy, cx, cy, maxiter);
60 // map linearily 1/(offset + iter) so that:
64 // f(iter) = A/(offset + iter) + B
66 // f(0) = A/offset + B = 0
67 // f(maxiter) = A/(offset + maxiter) + B = 1
69 // 1/(1/(offset + maxiter) - 1/offset) = A
70 // B = 1 + offset / maxiter
71 // A = -offset (1 + offset / maxiter)
73 // f(iter) = -offset (1 + offset / maxiter) / (offset + iter) + 1 + offset / maxiter
74 // = -offset (1 + offset / maxiter) / (offset + iter) + 1 + offset / maxiter
75 // = iter (offset + maxiter) / maxiter (offset + iter)
76 return (i * (offset + maxiter)) / ((i + offset) * maxiter);
79 double color_mandelbrot_parms[13];
80 double mandelbrot_miniter = -1;
83 double iter_mandelbrot_raw(double x, double y, double z)
85 z -= color_mandelbrot_parms[6];
90 return mandelbrot_range(color_mandelbrot_parms[4], color_mandelbrot_parms[5], color_mandelbrot_parms[0] + x * color_mandelbrot_parms[2], color_mandelbrot_parms[1] + y * color_mandelbrot_parms[3], MAXITER, color_mandelbrot_parms[9]);
95 void iter_mandelbrot_raw_initialize_min()
97 if(mandelbrot_miniter >= 0)
99 // randomly sample 256 points
101 // set that as miniter
104 mandelbrot_miniter = MAXITER;
105 for(i = 0; i < 8192; ++i)
110 double f = sqrt(x*x + y*y + z*z);
114 double a = (z - color_mandelbrot_parms[6]) / (color_mandelbrot_parms[7] - color_mandelbrot_parms[6]);
115 a = (a - color_mandelbrot_parms[8]) / (1 - color_mandelbrot_parms[8]);
118 double iterations = iter_mandelbrot_raw(x, y, z);
121 if(iterations < mandelbrot_miniter)
122 mandelbrot_miniter = iterations;
126 void color_mandelbrot(double x, double y, double z, double *r, double *g, double *b)
128 iter_mandelbrot_raw_initialize_min();
130 double iterations = iter_mandelbrot_raw(x, y, z);
131 //printf("iter = %f\n", iterations);
132 double a = (z - color_mandelbrot_parms[6]) / (color_mandelbrot_parms[7] - color_mandelbrot_parms[6]);
133 a = (a - color_mandelbrot_parms[8]) / (1 - color_mandelbrot_parms[8]);
138 iterations = iterations * a + mandelbrot_miniter * (1-a);
139 *r = pow(iterations, color_mandelbrot_parms[10]);
140 *g = pow(iterations, color_mandelbrot_parms[11]);
141 *b = pow(iterations, color_mandelbrot_parms[12]);
155 color_starfield_parms;
161 starfield_t *starfield = NULL;
162 int starfield_cmp(const void *a_, const void *b_)
164 const starfield_t *a = a_;
165 const starfield_t *b = b_;
172 void color_starfield(double x, double y, double z, double *r, double *g, double *b)
179 fprintf(stderr, "Initializing starfield...\n");
180 starfield = malloc(sizeof(*starfield) * color_starfield_parms.n);
181 for(i = 0; i < color_starfield_parms.n; ++i)
186 starfield[i].x = rnd() * 2 - 1;
187 starfield[i].y = rnd() * 2 - 1;
188 starfield[i].z = rnd() * 2 - 1;
189 r = starfield[i].x * starfield[i].x
190 + starfield[i].y * starfield[i].y
191 + starfield[i].z * starfield[i].z;
195 starfield[i].R = rnd();
196 starfield[i].G = rnd();
197 starfield[i].B = rnd();
198 f = starfield[i].R * 0.299 + starfield[i].G * 0.587 + starfield[i].B * 0.114;
202 starfield[i].A = rnd();
203 starfield[i].e = color_starfield_parms.density * pow(rnd(), -color_starfield_parms.power);
205 fprintf(stderr, "Gravitating starfield...\n");
206 for(k = 0; k < color_starfield_parms.g; ++k)
208 i = rand() % color_starfield_parms.n;
209 j = rand() % color_starfield_parms.n;
210 f = pow(rnd(), color_starfield_parms.gpower);
211 f = f * color_starfield_parms.gfactor;
212 d = pow(starfield[j].x - starfield[i].x, 2)
213 + pow(starfield[j].y - starfield[i].y, 2)
214 + pow(starfield[j].z - starfield[i].z, 2);
215 f *= pow(1 / (1 + d), color_starfield_parms.gdpower);
216 starfield[i].x += f * (starfield[j].x - starfield[i].x);
217 starfield[i].y += f * (starfield[j].y - starfield[i].y);
218 starfield[i].z += f * (starfield[j].z - starfield[i].z);
219 double r = starfield[i].x * starfield[i].x + starfield[i].y * starfield[i].y + starfield[i].z * starfield[i].z;
225 fprintf(stderr, "Sorting starfield...\n");
226 qsort(starfield, sizeof(*starfield), color_starfield_parms.n, starfield_cmp);
227 fprintf(stderr, "Done.\n");
231 mindot = pow(1/256.0, 1.0/color_starfield_parms.density);
232 for(i = 0; i < color_starfield_parms.n; ++i)
234 double dot = x * starfield[i].x + y * starfield[i].y + z * starfield[i].z;
237 double f = pow(dot, starfield[i].e) * starfield[i].A;
238 *r += starfield[i].R * f;
239 *g += starfield[i].G * f;
240 *b += starfield[i].B * f;
243 *r = *r / (color_starfield_parms.darkness + *r);
244 *g = *g / (color_starfield_parms.darkness + *g);
245 *b = *b / (color_starfield_parms.darkness + *b);
248 void map_back(double x_in, double y_in, double *x_out, double *y_out, double *z_out)
250 *x_out = 2 * x_in - 1;
252 *z_out = 1 - 2 * y_in;
255 void map_right(double x_in, double y_in, double *x_out, double *y_out, double *z_out)
258 *y_out = 1 - 2 * x_in;
259 *z_out = 1 - 2 * y_in;
262 void map_front(double x_in, double y_in, double *x_out, double *y_out, double *z_out)
264 *x_out = 1 - 2 * x_in;
266 *z_out = 1 - 2 * y_in;
269 void map_left(double x_in, double y_in, double *x_out, double *y_out, double *z_out)
272 *y_out = 2 * x_in - 1;
273 *z_out = 1 - 2 * y_in;
276 void map_up(double x_in, double y_in, double *x_out, double *y_out, double *z_out)
278 *x_out = 2 * y_in - 1;
279 *y_out = 1 - 2 * x_in;
283 void map_down(double x_in, double y_in, double *x_out, double *y_out, double *z_out)
285 *x_out = 1 - 2 * y_in;
286 *y_out = 1 - 2 * x_in;
290 void writepic(colorfunc_t f, mapfunc_t m, const char *fn, int width, int height)
296 FILE *file = fopen(fn, "wb");
298 err(1, "fopen %s", fn);
300 memset(tga, 0, sizeof(tga));
301 tga[2] = 2; // uncompressed type
302 tga[12] = (width >> 0) & 0xFF;
303 tga[13] = (width >> 8) & 0xFF;
304 tga[14] = (height >> 0) & 0xFF;
305 tga[15] = (height >> 8) & 0xFF;
306 tga[16] = 24; // pixel size
308 fwrite(&tga, sizeof(tga), 1, file);
310 for(y = height-1; y >= 0; --y)
312 for(x = 0; x < width; ++x)
317 double xxx, yyy, zzz;
319 xx = (x + 0.5) / width;
320 yy = (y + 0.5) / height;
321 m(xx, yy, &xxx, &yyy, &zzz);
322 r = sqrt(xxx*xxx + yyy*yyy + zzz*zzz);
326 f(xxx, yyy, zzz, &rr, &gg, &bb);
327 rgb[2] = floor(rnd() + rr * 255);
328 rgb[1] = floor(rnd() + gg * 255);
329 rgb[0] = floor(rnd() + bb * 255);
330 fwrite(rgb, sizeof(rgb), 1, file);
332 p = (100 * (height - y)) / height;
336 fprintf(stderr, "%d%%\r", percent);
339 fprintf(stderr, "\n");
344 void map_all(const char *fn, colorfunc_t f, int width, int height)
347 fprintf(stderr, "%s_bk.tga\n", fn);
348 snprintf(buf, sizeof(buf), "%s_bk.tga", fn); buf[sizeof(buf) - 1] = 0; writepic(f, map_back, buf, width, height);
349 fprintf(stderr, "%s_ft.tga\n", fn);
350 snprintf(buf, sizeof(buf), "%s_ft.tga", fn); buf[sizeof(buf) - 1] = 0; writepic(f, map_front, buf, width, height);
351 fprintf(stderr, "%s_rt.tga\n", fn);
352 snprintf(buf, sizeof(buf), "%s_rt.tga", fn); buf[sizeof(buf) - 1] = 0; writepic(f, map_right, buf, width, height);
353 fprintf(stderr, "%s_lf.tga\n", fn);
354 snprintf(buf, sizeof(buf), "%s_lf.tga", fn); buf[sizeof(buf) - 1] = 0; writepic(f, map_left, buf, width, height);
355 fprintf(stderr, "%s_up.tga\n", fn);
356 snprintf(buf, sizeof(buf), "%s_up.tga", fn); buf[sizeof(buf) - 1] = 0; writepic(f, map_up, buf, width, height);
357 fprintf(stderr, "%s_dn.tga\n", fn);
358 snprintf(buf, sizeof(buf), "%s_dn.tga", fn); buf[sizeof(buf) - 1] = 0; writepic(f, map_down, buf, width, height);
361 int main(int argc, char **argv)
365 errx(1, "usage: %s filename res func parms...", *argv);
366 int res = atoi(argv[2]);
367 if(!strcmp(argv[3], "mandel"))
369 f = color_mandelbrot;
370 color_mandelbrot_parms[0] = argc<= 4 ? -0.740 : atof(argv[4]); // shift xy
371 color_mandelbrot_parms[1] = argc<= 5 ? -0.314 : atof(argv[5]);
372 color_mandelbrot_parms[2] = argc<= 6 ? -0.003 : atof(argv[6]); // mul xy
373 color_mandelbrot_parms[3] = argc<= 7 ? -0.003 : atof(argv[7]);
374 color_mandelbrot_parms[4] = argc<= 8 ? 0.420 : atof(argv[8]); // shift z
375 color_mandelbrot_parms[5] = argc<= 9 ? 0.000 : atof(argv[9]);
376 color_mandelbrot_parms[6] = argc<=10 ? -0.8 : atof(argv[10]); // horizon
377 color_mandelbrot_parms[7] = argc<=11 ? -0.7 : atof(argv[11]);
378 color_mandelbrot_parms[8] = argc<=12 ? 0.5 : atof(argv[12]);
379 color_mandelbrot_parms[9] = argc<=13 ? 400 : atof(argv[13]); // coloring
380 color_mandelbrot_parms[10] = argc<=14 ? 0.6 : atof(argv[14]);
381 color_mandelbrot_parms[11] = argc<=15 ? 0.5 : atof(argv[15]);
382 color_mandelbrot_parms[12] = argc<=16 ? 0.2 : atof(argv[16]);
384 else if(!strcmp(argv[3], "starfield"))
387 color_starfield_parms.n = argc<= 4 ? 8192 : atoi(argv[4]);
388 color_starfield_parms.darkness = argc<= 5 ? 0.4 : atof(argv[5]);
389 color_starfield_parms.power = argc<= 6 ? 2.5 : atof(argv[6]);
390 color_starfield_parms.density = argc<= 7 ? 60000 : atof(argv[7]);
391 color_starfield_parms.g = argc<= 8 ? 10000000 : atoi(argv[8]);
392 color_starfield_parms.gpower = argc<= 9 ? 3 : atoi(argv[9]);
393 color_starfield_parms.gfactor = argc<= 10 ? 0.9 : atof(argv[10]);
394 color_starfield_parms.gdpower = argc<= 11 ? 15 : atof(argv[11]);
400 map_all(argv[1], f, res, res);