|
|
waggy wrote:
> Oh, and here is the same view with the isosurface function limited to eight
> iterations.
>
> I just found the function also works well, and is far, far quicker, as a media
> density.
>
> Good times.
>
> ~David
Well, I'm seeing that the actual surface generation and rendering is
probably best left to people like you who seem to get good results, so I
figured I might as well go another direction with it.
Instead of going for the full 3d, as many of you are, I decided to make
a way to get 2d 'slices' of this beast that allow for nice smooth
coloring based on divergence rate, like most of the pictures for the
original Mandelblot, instead of only edge finding like everyone seems to
be doing with this so far. I wrote some code that will generate
grayscale images based on number of iterations according to various
parameters.
Here is one slice, generated with parameters mirroring povray:
center = <0.7275,0.272,0.5>
right = <0.001,0,0>
up = <0,0.001,0>
n = 8
max_iterations = 200
This image is NOT pov-created, but can be used for patterns and the
like. It was originally generated at 5000x5000 (.png) which took about a
minute, then downsized for uploading here.
C++ source available on request, the only dependencies are ImageMagick
for creating the image file itself, and a compiler supporting OpenMP.
Post a reply to this message
Attachments:
Download 'mb_slice_demo.jpg' (270 KB)
Preview of image 'mb_slice_demo.jpg'
|
|
|
|
"scott" <sco### [at] scottcom> wrote:
> Hehe I just did the same :-) Mine was 1024x1024x1024 df3 file which took
> about half an hour multithreaded under c#.net (I decided the speedup I could
> get by easy multithreading in c# outweighed the speedup I could get with C++
> over C#). The render only took about 2 minutes, it's just an isosurface and
> a few lights.
Nice work, Scott. And others!
I have a slow old machine, with only 768MB of RAM, so my renders take a bit
I'm especially jealous of that special version of povray that can do custom
fractal functions. :) But I suppose it's probably slower than using a density
file. (I haven't played with them much before, and I'm frankly amazed at how
fast those things parse, even fairly huge ones.)
I *think* I figured a way to do this Mandelbulb function using standard SDL. We
don't need the ability to do arbitrarily deep iterations, since 6 or
7 are quite adequate. So we could unroll the iteration loop... into a veritable
rats' nest of functions, which would probably involve lots of duplicated
evaluations & be incredibly slow. :)
> > I haven't optimized the code too far yet, it's still using the trig
> > functions, but it can be used for any value N.
>
> Yeh mine too, it looks like there should be a way to lose the trig even with
> N=8, but it might get a bit scary...
You can lose the trig functions by using the formulae of sin(2t) & cos(2t) in
terms of sin(t) & cos(t). We want the 8th power, so we just double the angle 3
times. My first explorations of the Mandelbulb used Python to generate the voxel
data, but I translated my code to C today, and it should compile with no
problems on any system (if it supports doubles :)).
To avoid clutter, I'll post my C source in a separate message in this thread.
> The issue is that the isosurface is still made up of millions of little cube
> shaped voxels, I think the only way to improve that is to write a
> fractal-ray tracer too! Or maybe patch POV...
Your renders look a fair bit smoother than mine. What accuracy are you using on
the isosurface, and what sort of max_gradient for that size grid? The one below
uses max_gradient 1072, accuracy 0.002.
The pigment on the attached image is just function{sqrt(x*x+z*z)*4.5}, with an
old favourite colour map I created many years ago when I was writing 2D
Mandelbrot programs on the Amiga.
Post a reply to this message
Attachments:
Download 'mandelbulbg8sa90.jpg' (160 KB)
Preview of image 'mandelbulbg8sa90.jpg'
|
|
|
|
Simple "Mandelbulb" C source code.
Compile with
gcc -o Mandelbulb Mandelbulb.c -lm -O3
/*
* Mandelbulb
*
* Written by PM 2Ring November 2009.
* Converted from Python to C, 2009.11.21
*
*/
/*
From http://www.skytopia.com/project/fractal/mandelbulb.html
Similar to the original 2D Mandelbrot , the 3D formula is defined by:
z -> z^n + c
....but where 'z' and 'c' are hypercomplex ('triplex') numbers, representing
Cartesian x, y, and z coordinates. The exponentiation term is defined by:
{x,y,z}^n = r^n {sin(theta*n) * cos(phi*n), sin(theta*n) * sin(phi*n),
cos(theta*n)}
where
r = sqrt(x^2 + y^2 + z^2)
theta = atan2(sqrt(x^2+y^2), z)
phi = atan2(y, x)
n is the order of the 3D Mandelbulb. Use n=8 to find the exact object in this
article.
The addition term in z -> z^n + c is similar to standard complex addition,
and is simply defined by:
{x,y,z}+{a,b,c} = {x+a, y+b, z+c}
The rest of the algorithm is similar to the 2D Mandelbrot!
*/
#include <math.h>
#include <stdio.h>
#include <time.h>
//Compute sin(8t) & cos(8t) given r*sin(t), r*cos(t) & r
void trig8(double *ps, double *pc, double r)
{
//int i;
double c = *pc, s = *ps, t;
if (r > 0.0)
{
s /= r;
c /= r;
//Iterate formulae for sin(2t) & cos(2t) 3 times
//to get sin(8t) & cos(8t)
//for (i=0; i<3; i++)
//{
t = 2.0*s*c; c = c*c - s*s; s = t;
t = 2.0*s*c; c = c*c - s*s; s = t;
t = 2.0*s*c; c = c*c - s*s; s = t;
//}
*ps = s;
*pc = c;
}
}
//Determine if voxel (cx, cy, cz) is in the Mandelbulb
int Mandel3D(double cx, double cy, double cz, int iters)
{
double x = cx, y = cy, z = cz, w, w2, r, r2, r8,
sth8, cth8, sph8, cph8;
int i;
for (i=0; i<iters; i++)
{
//Convert (x, y, z) coordinates to polar (r, th, ph)
w2 = x*x + y*y;
w = sqrt(w2);
r2 = w2 + z*z;
r = sqrt(r2);
if (r >= 2.0)
return 0;
//Raise (r, th, ph) to the 8th power
sth8 = w; cth8 = z; trig8(&sth8, &cth8, r);
sph8 = y; cph8 = x; trig8(&sph8, &cph8, w);
//(r, th, ph) <- c + (r, th, ph)^8
r8 = r2 * r2 * r2 * r2;
x = cx + r8 * sth8 * cph8;
y = cy + r8 * sth8 * sph8;
z = cz + r8 * cth8;
}
return 1;
}
int Mandelbulb(double vx, double vy, double vz, double scale, int side, int
iters, char *fname)
{
FILE *f;
unsigned char *row;
int i, j, k;
time_t itime;
double x, y, z;
if(!(f = fopen(fname, "wb")))
return 1;
if(!(row = (unsigned char *) malloc(side)))
return 1;
//Write header: voxel grid dimensions in little endian
{
unsigned char lo = side & 0xff, hi = side >> 8;
fprintf(f, "%c%c%c%c%c%c", hi, lo, hi, lo, hi, lo);
}
itime = time(NULL);
for (k=0; k<side; k++)
{
//Planes
z = vz + scale * k;
for (j=0; j<side; j++)
{
//Rows
y = vy + scale * j;
for (i=0; i<side; i++)
//Voxels
row[i] = Mandel3D(vx + scale * i, y, z, iters) ? 255 : 0;
fwrite(row, 1, side, f);
}
if (k%10 == 0)
{
int n = (int) difftime(time(NULL), itime);
//fputc('.', stderr);
fprintf(stderr, "\r %2d:%02d:%02d Plane %d of %d ",
n/3600, (n/60)%60, n%60, k, side);
}
}
free(row);
fclose(f);
return 0;
}
int main(int argc, char *argv[])
{
int side, iters, rc;
char *fname;
double r, scale;
side = argc > 1 ? atoi(argv[1]) : 100;
fname = argc > 2 ? argv[2] : "Mandelbulb.df3";
iters = argc > 3 ? atoi(argv[3]) : 7;
r = 1.1; //Mandelbulb Set "radius"
//r = 1.05 * (1.0 + 2.0 / side);
scale = 2.0*r / (side - 1);
rc = Mandelbulb(-r, -r, -r, scale, side, iters, fname);
fputc('\n', stderr);
return rc;
}
Post a reply to this message
|
|