True 3D mandelbrot type fractal

[David Makin]

I'm interested to know if this 3D fractal figure could be output as a CAD file like a .dxf, .stl or .aoi format? 

& would anyone share the complete code for this with me? I could write the file format output part myself if I had that.

What I would like to do is have one of these rendered as a an actual object by laser sintering, fdm, cnc or some similar method.

It'd look really cool in brass or cut from silicon

See:

http://www.fractalforums.com/3d-fractal-generation/finding-the-edges/msg8626/#msg8626

[David Makin]

<snip>
In the meanwhile I think I should try David-Makin´s way of DE, which only uses 2 points (mine uses 4).
<snip>
JosLeys, so in the end which of the formulas is the correct one for DE?

The key part of my delta DE is the storing/testing of the maxdist values at a given iteration for a given ray, for each new step distance calculated I use:

              if j<=miter
                if step<dists[j]
                  repeat
                    dists[j] = step
                    j = j + 1
                  until step>=dists[j] || j>miter
                  if j>miter
                    mdist = step
                  endif
                else
                  step = dists[j]
                endif
              else;if j>miter
                miter = miter + 1
                if step>mdist
                  step = mdist
                endif
                while miter<j
                  dists[(miter=miter+1)] = mdist
                endwhile
                dists[j] = step
              endif

Where mdist is the absolute minimum found on the ray and miter is the maximum iterations used on the ray.

With respect to the scale factor for DE I'm just about to do some testing - I think I have an older version of my formula that calculates in the global section so I can check values the way Jos did but for an entire render by setting the threshold very small indeed and assuming the final value is "correct", this should give a good idea whether values at 100* greater distance are accurate or not.

Note that using a maxdist array in this manner also cures the problem of singularities in the derivative (where it's zero or near zero) when using the analytical DE - such problems occur especially on Julia Sets.

[bugman]

I'm interested to know if this 3D fractal figure could be output as a CAD file like a .dxf, .stl or .aoi format? 

& would anyone share the complete code for this with me? I could write the file format output part myself if I had that.

What I would like to do is have one of these rendered as a an actual object by laser sintering, fdm, cnc or some similar method.

It'd look really cool in brass or cut from silicon

Yes. That would be neat. I also have a contact who has offered to create a rapid prototype of this for me. I don't think triangle formats would be very practical for this, but I think voxel data would be practical.

Also, a laser etched crystal could be created from a point cloud, which could be calculated using inverse methods if you are creating a 3D Julia set. For example, look at this beautiful one by Bathsheba Grossman:
http://www.bathsheba.com/crystal/julia/

She explains the process here:
http://www.bathsheba.com/crystal/process/

[lycium]

Can anyone please send me £2000 so I can get a new system ? (this P4HT is a dinosaur - only an ATI X600 too)

i really don't want to take this offtopic thread here, but hopefully someone else can convince you that buying a dual-processor system is going to cause more money problems than it solves computing problems. why is it overpriced dual processor apple or bust?


anyway back on topic, those are really cool renders iq! i especially like the neat shading trick, hope to give that a try soon. i've also started work on a gpu implementation, though i doubt my crusty old 8800 gts will get 2 seconds per frame at that resolution. do you think it's the DE that's giving that kind of performance? my theory is that the DE regularises the computation, so neighbouring pixels will have less branch divergence.

right now i'm fighting with the watchdog timer on windows 7, which is resetting my display driver after 2 seconds of computation *grumblegrumblegrumble*

[David Makin]

After Jos Leys tests using analytical DE on the zero seed Julia I thought I'd try investigation of analytical DE values on a quaternionic Mandelbrot.

I quickly found that indeed the "0.5" is correct which means that the actual DE values in my wip formula for Ultra Fractal are 2* bigger than they should be - though my step distances are OK since the default is the DE divided by two This of course results in a given DE threshold taking us to a point twice as close to "inside" as the specified threshold distance.

Here are the test results I did using the correctly scaled DE value.

These tests were carried out @640*480 with a magnification of 1.6 and a bailout of x^2+y^2>=1024.
Note that in the text below "worst" means "largest" rather than "most incorrect".

Quaternion z^2+c side view in parallel projection
DE Threshold 1e-6
Worst cases for entire fractal

"True" -> calculated

1.0 -> 0.828693280840216
0.1 -> 0.0756111410237615
0.01 -> 0.00826817743483479
1e-3 -> 0.00623750530621662
1e-4 -> 0.00616464060211196
1e-5 -> 0.00593556668345077

Quaternion z^2+c side view in parallel projection
DE Threshold 1e-6
Average cases for entire fractal

"True" -> calculated

1.0 -> 0.470217800965237
0.1 -> 0.0277512230487669
0.01 -> 0.0027385391736115
1e-3 -> 0.000537393167849202
1e-4 -> 0.000298006036418662
1e-5 -> 0.000183765097401958

Quaternion z^2+c pointed end view in parallel projection
DE Threshold 1e-6
Worst cases for entire fractal

"True" -> calculated

1.0 -> 0.97783685049528
0.1 -> 0.0971841767897501
0.01 -> 0.0099424008139091
1e-3 -> 0.00291127344517454
1e-4 -> 0.0028424928302291
1e-5 -> 0.00283555485322532

Quaternion z^2+c pointed end view in parallel projection
DE Threshold 1e-6
Average cases for entire fractal

"True" -> calculated

1.0 -> 0.268576868120558
0.1 -> 0.0330877709877468
0.01 -> 0.00333153377206449
1e-3 -> 0.000601747860965648
1e-4 -> 0.000318767428809674
1e-5 -> 0.000195952986975526

Quaternion z^2+c pointed end view in parallel projection
DE Threshold 1e-6
Values for the centre ray (i.e. hitting the point on the real axis)

"True" -> calculated

1.0 -> 0.97783685049528
0.1 -> 0.0971841767897501
0.01 -> 0.0099424008139091
1e-3 -> 0.000999129722927703
1e-4 -> 0.00010003123154646
1e-5 -> 9.99071625909761E-006

Quaternion z^2+c side view in parallel projection
DE Threshold 1e-16
Worst cases for entire fractal

"True" -> calculated

1.0 -> 0.828694191861107
0.1 -> 0.0756118350632922
0.01 -> 0.00826812262793161
1e-3 -> 0.00623734563239439
1e-4 -> 0.00616447654586974
1e-5 -> 0.00593548110247801
1e-6 -> 0.002967740551239
1e-7 -> 0.0014838702756195
1e-8 -> 0.000741935137809751
1e-9 -> 0.000370967568904876
1e-10 -> 0.000185483784452438
1e-11 -> 9.27418922262189E-005
1e-12 -> 4.63709461131094E-005

Quaternion z^2+c side view in parallel projection
DE Threshold 1e-16
Average cases for entire fractal

"True" -> calculated

1.0 -> 0.470218381043924
0.1 -> 0.0277514611184807
0.01 -> 0.00273875075822565
1e-3 -> 0.000537604762798528
1e-4 -> 0.000298158908926989
1e-5 -> 0.000185850698115385
1e-6 -> 9.97458957492159E-005
1e-7 -> 5.0596384017043E-005
1e-8 -> 2.53716863807056E-005
1e-9 -> 1.26919229340975E-005
1e-10 -> 6.34662886913257E-006
1e-11 -> 3.17331443456628E-006
1e-12 -> 1.58665721728314E-006

Quaternion z^2+c pointed end view in parallel projection
DE Threshold 1e-16
Worst cases for entire fractal

"True" -> calculated

1.0 -> 0.977836860473067
0.1 -> 0.097184186106234
0.01 -> 0.00994241047465512
1e-3 -> 0.00291124435480221
1e-4 -> 0.00284246358773596
1e-5 -> 0.00283552550759178
1e-6 -> 0.0019832407359585
1e-7 -> 0.00182386184009324
1e-8 -> 0.000911930920046621
1e-9 -> 0.000786028563546067
1e-10 -> 0.000786028308654553
1e-11 -> 0.00078602828316538
1e-12 -> 0.000786028280616523

Quaternion z^2+c pointed end view in parallel projection
DE Threshold 1e-16
Average cases for entire fractal

"True" -> calculated

1.0 -> 0.268571283431845
0.1 -> 0.0330876101698883
0.01 -> 0.00333120301202227
1e-3 -> 0.000602036669078832
1e-4 -> 0.000319200775195032
1e-5 -> 0.000198471399420637
1e-6 -> 0.000106430696869362
1e-7 -> 5.41934392969572E-005
1e-8 -> 2.71886867735502E-005
1e-9 -> 1.36275997817489E-005
1e-10 -> 6.84249791092663E-006
1e-11 -> 3.44994698460904E-006
1e-12 -> 1.75367152223619E-006

Quaternion z^2+c pointed end view in parallel projection
DE Threshold 1e-16
Values for the centre ray (i.e. hitting the point of the Mandy on the real axis)

"True" -> calculated

1.0 -> 0.977836860473067
0.1 -> 0.097184186106234
0.01 -> 0.00994241047465512
1e-3 -> 0.000999139471789839
1e-4 -> 0.000100040977667748
1e-5 -> 1.00004784004224E-005
1e-6 -> 1.00006323931994E-006
1e-7 -> 1.001088304125E-007
1e-8 -> 1.00620844487879E-008
1e-9 -> 1.0638804422352E-009
1e-10 -> 1.61612459695359E-010
1e-11 -> 7.17222631161621E-011
1e-12 -> 3.5861131558081E-011

I think the above shows that the problem of areas with small gradient causing erroneously high DE values is more prevalent than I previously thought even in a Mandelbrot As I've mentioned previously a potential cure for this is develop a DE algorithm that incorporates the second derivative in some way. For instance the way I derived the standard DE formula was based on using the Newton, a similar derivation could be made say using Householder's formula.

In the last set of figures I think the 1e-10, 1e-11 and 1e-12 results may be affected by the likelyhood that y and z in quaternionic (x,y,z,w) weren't exactly zero due to slight floating point errors in the 90 degree rotation - i.e. we didn't exactly hit the point.

I intend to redo these tests with say a degree 16 quaternion to see if the values confirm Jos Leys' suggestion that no scaling is required based on the power or to show if such scaling is necessary - and if necessary then to try and find what it should be.

[iq]

interesting thoughts about the second derivative. For regular complex M-sets that should give these iteration rules:

ddz = p·(p-1)·z^(p-2)·dz² + p·z^(p-1)·ddz
dz = p·z^(p-1)·dz + 1
z = z^p + c

for p=2:

ddz = 2·dz² + 2·z·ddz
dz = 2·z·dz + 1
z = z² + c

I suppose for the first iteration z = 0, dz = 1, ddz = 0

The first Taylor expansion of G(c+eps) is G(c) + eps·G'(c) which leads to eps = distance = G/|G'| = the famous |z| · log|z| / |dz|

So the second order expansion should be something like G(c+eps) = G(c) + eps·G'(c) + eps²/2 · G''(c), so eps = distance = (solving the quadratic) =  (-|G'| ± sqrt( |G'|² - 2·|G''|·|G|))/|G''| I guess the solution we want is the one with the negative square root... As usual,

|G| = (1/2^n) · log |z|
|G'| = (1/2^n) · |dz| / |z|
|G''| = (1/2^n) · (|ddz|·|z|-|dz|²) / |z|²

which means that we have a huge monster:

k = |dz|² - |z|·|ddz|

distance = ( -|z|·|dz| ± sqrt(|dz|² - 2·k·log|z|) )/k

which is quite more complex that the first order distance estimator |z|·log|z|/|dz|  I might very well have made mistakes in this deductive process. Most likely with the signs...

I'm curious to see if it works, at least in the old-good-2d M-Set. Can somebody try?

[eternaleye]

I find myself wondering if there could be a 'lazy' way to voxel-map this shape. Bear with me a sec while I figure this out.

Okay, so what you guys are doing now is creating a volume of x*y*z voxels, and iterating over those voxels checking whether there is or is not a portion of the shape in that voxel; am I correct?

If I am, there is a much more efficient way to plot it.

The trick is not to directly iterate over the voxels; rather, one subdivides the entire volume, ignoring sectors where no voxels contain parts of the shape, and recursing into sectors where some or all of the voxels _do_.

Let's say we had an 8x8x8 volume we're plotting it into.

First, we divide it into 8 4x4x4 volumes, and check whether any of these volumes is completely empty. This can be done identically to mapping the shape into a 2x2x2 volume. We then never visit such empty subvolumes again.

For each subvolume that _does_ contain part of the shape, we again split it and perform the same action. Any sub-subvolumes that are empty are discarded, and the voxels in non-empty sub-subvolumes are iterated over.

A somewhat inconvenient downside is the requirement of x==y==z && x == 2 ** n, where n is any integer.
A rather convenient upside is that, if one has exported the volumetric bitmap that results, it is trivial to import it and use it to double the (linear) resolution without recalculating the previous work.

Here's some rough C++ish pseudocode:

Code:

#include <queue>
// I'm going to assume that the function to check whether a voxel is occupied is called 'bool is_v_occ( int x, int y, int z, int size )',
// where 'size' is the length of each side of the volume the voxel is in.
// I will also assume that a 3d C++ vector of bool has already been allocated under the name 'volume'.

struct SectorSpec {
int stride, x, y, z;
SectorSpec( int str, int xnew, int ynew, int znew) {
stride = str;
x = xnew;
y = ynew;
z = znew;
}
}

void map_volume( SectorSpec* sector, Volume& vol, std::queue<SectorSpec*>& work ) {
if ( sector->stride == 1 ) {
// A sector is always 2x2x2 strides
for ( int x = sector->x; x < sector->x + 2; ++x ) {
for ( int y = sector->y; y < sector->y + 2; ++y ) {
for ( int z = sector->z; z < sector->z + 2; ++z ) {
vol[x][y][z] = is_v_occ( x, y, z, vol.size() );
}
}
}
} else {
int sz = vol.size() / sector->stride;
for ( int x = sector->x / sector->stride; x < ( sector->x / sector->stride ) + 2; ++x ) {
for ( int y = sector->y / sector->stride; y < ( sector->y / sector->stride ) + 2; ++y ) {
for ( int z = sector->z / sector->stride; z < ( sector->z / sector->stride ) + 2; ++z ) {
if( is_v_occ( x, y, z, sz ) ) {
SectorSpec* newsec = new SectorSpec(
sector->stride / 2,
x * sector->stride,
y * sector->stride,
z * sector->stride
);
work.push( newspec )
} else {
for ( int xinner = sector->x; x < sector->x + ( 2 * sector->stride ); ++x ) {
for ( int yinner = sector->y; y < sector->y + ( 2 * sector->stride ); ++y ) {
for ( int zinner = sector->z; z < sector->z + ( 2 * sector->stride ); ++z ) {
vol[xinner][yinner][zinner] = false;
}
}
}
}
}
}
}
}
}

int main() {
std::queue<SectorSpec*> sectors;
SectorSpec* seed = new SectorSpec( volume.size, 0, 0, 0 );
map_volume( seed, volume, sectors );
for ( SectorSpec* nextSec = sectors.pop(); !sectors.empty(); nextSec = sectors.pop() ) {
map_volume( nextSec, volume, sectors )
}
// volume now contains a 3d bitmap of occupied voxels. Congrats! Add your code here to do as you wish.
return( 0 );
}

EDIT: I just realized that the 'inconvenient downside' above only holds true for a naive implementation such as the one I have here. A more sophisticated implementation could eliminate that by some boundschecking and dimension-dependent stride values.

[David Makin]

I find myself wondering if there could be a 'lazy' way to voxel-map this shape. Bear with me a sec while I figure this out.

Okay, so what you guys are doing now is creating a volume of x*y*z voxels, and iterating over those voxels checking whether there is or is not a portion of the shape in that voxel; am I correct?

If I am, there is a much more efficient way to plot it.

The trick is not to directly iterate over the voxels; rather, one subdivides the entire volume, ignoring sectors where no voxels contain parts of the shape, and recursing into sectors where some or all of the voxels _do_.

Let's say we had an 8x8x8 volume we're plotting it into.

First, we divide it into 8 4x4x4 volumes, and check whether any of these volumes is completely empty. This can be done identically to mapping the shape into a 2x2x2 volume. We then never visit such empty subvolumes again.

For each subvolume that _does_ contain part of the shape, we again split it and perform the same action. Any sub-subvolumes that are empty are discarded, and the voxels in non-empty sub-subvolumes are iterated over.

A somewhat inconvenient downside is the requirement of x==y==z && x == 2 ** n, where n is any integer.
A rather convenient upside is that, if one has exported the volumetric bitmap that results, it is trivial to import it and use it to double the (linear) resolution without recalculating the previous work.

Here's some rough C++ish pseudocode:

Code:

#include <queue>
// I'm going to assume that the function to check whether a voxel is occupied is called 'bool is_v_occ( int x, int y, int z, int size )',
// where 'size' is the length of each side of the volume the voxel is in.
// I will also assume that a 3d C++ vector of bool has already been allocated under the name 'volume'.

struct SectorSpec {
int stride, x, y, z;
SectorSpec( int str, int xnew, int ynew, int znew) {
stride = str;
x = xnew;
y = ynew;
z = znew;
}
}

void map_volume( SectorSpec* sector, Volume& vol, std::queue<SectorSpec*>& work ) {
if ( sector->stride == 1 ) {
// A sector is always 2x2x2 strides
for ( int x = sector->x; x < sector->x + 2; ++x ) {
for ( int y = sector->y; y < sector->y + 2; ++y ) {
for ( int z = sector->z; z < sector->z + 2; ++z ) {
vol[x][y][z] = is_v_occ( x, y, z, vol.size() );
}
}
}
} else {
int sz = vol.size() / sector->stride;
for ( int x = sector->x / sector->stride; x < ( sector->x / sector->stride ) + 2; ++x ) {
for ( int y = sector->y / sector->stride; y < ( sector->y / sector->stride ) + 2; ++y ) {
for ( int z = sector->z / sector->stride; z < ( sector->z / sector->stride ) + 2; ++z ) {
if( is_v_occ( x, y, z, sz ) ) {
SectorSpec* newsec = new SectorSpec(
sector->stride / 2,
x * sector->stride,
y * sector->stride,
z * sector->stride
);
work.push( newspec )
} else {
for ( int xinner = sector->x; x < sector->x + ( 2 * sector->stride ); ++x ) {
for ( int yinner = sector->y; y < sector->y + ( 2 * sector->stride ); ++y ) {
for ( int zinner = sector->z; z < sector->z + ( 2 * sector->stride ); ++z ) {
vol[xinner][yinner][zinner] = false;
}
}
}
}
}
}
}
}
}

int main() {
std::queue<SectorSpec*> sectors;
SectorSpec* seed = new SectorSpec( volume.size, 0, 0, 0 );
map_volume( seed, volume, sectors );
for ( SectorSpec* nextSec = sectors.pop(); !sectors.empty(); nextSec = sectors.pop() ) {
map_volume( nextSec, volume, sectors )
}
// volume now contains a 3d bitmap of occupied voxels. Congrats! Add your code here to do as you wish.
return( 0 );
}

I know what you're saying and you're essentially correct, however applying such an algorithm to fractals is similar to optimisations you can make in 2D (e.g. "guessing" mode in UF) and tends to result in missed solid especially if the fractal is not connected, having said that I just realised that using analytical distance estimator for the "empty voxel" test woiuld probably largely avoid the problems.

[iq]

The technical name for what you are describing is "octree". I'm nos sure how many people are using voxels to do their renderings, tho, I though most are doing direct raymarching?

bib, very nice video! I like the zooming.

I have done a little video. Well, in fact I had to stop it in half rendering because after a few minutes my GPU was so hot that my computer started to show garbage pixels in the screen, and then it was all broken and unreadable!!! So I reimplemented the thing in the CPU and started over again. For 1280x720 I'm getting around 4 seconds rendering in the CPU, with a 8 core CPU (so it should be like half minute in a monocore).

Anyway, this is the half video I rendered. I will re-render it in the CPU (don't wanna risk my HW again, and I don't want to start putting Sleep() every few frames to let the GPU take some breath) and upload to the same place, so you can try the link again by tomorrow or so for the full video.

I personally think that just rotation around the set is not as interesting as zomming in as in bib's video. But well, it's a start.

www.iquilezles.org/trastero/f8p5.avi

[David Makin]

interesting thoughts about the second derivative. For regular complex M-sets that should give these iteration rules:

ddz = p·(p-1)·z^(p-2)·dz² + p·z^(p-1)·ddz
dz = p·z^(p-1)·dz + 1
z = z^p + c

for p=2:

ddz = 2·dz² + 2·z·ddz
dz = 2·z·dz + 1
z = z² + c

I suppose for the first iteration z = 0, dz = 1, ddz = 0

The first Taylor expansion of G(c+eps) is G(c) + eps·G'(c) which leads to eps = distance = G/|G'| = the famous |z| · log|z| / |dz|

So the second order expansion should be something like G(c+eps) = G(c) + eps·G'(c) + eps²/2 · G''(c), so eps = distance = (solving the quadratic) =  (-|G'| ± sqrt( |G'|² - 2·|G''|·|G|))/|G''| I guess the solution we want is the one with the negative square root... As usual,

|G| = (1/2^n) · log |z|
|G'| = (1/2^n) · |dz| / |z|
|G''| = (1/2^n) · (|ddz|·|z|-|dz|²) / |z|²

which means that we have a huge monster:

k = |dz|² - |z|·|ddz|

distance = ( -|z|·|dz| ± sqrt(|dz|² - 2·k·log|z|) )/k

which is quite more complex that the first order distance estimator |z|·log|z|/|dz|  I might very well have made mistakes in this deductive process. Most likely with the signs...

I'm curious to see if it works, at least in the old-good-2d M-Set. Can somebody try?

Will try this evening...

[iq]

fantastic David! I'm done with fractals for today I think, it's the turn of my CPUs to work now. As my last contribution of the day to the thread, I link to something funny I found out there:



Definitely humans are designed to find patterns.

[David Makin]

The technical name for what you are describing is "octree". I'm nos sure how many people are using voxels to do their renderings, tho, I though most are doing direct raymarching?

bib, very nice video! I like the zooming.

I have done a little video. Well, in fact I had to stop it in half rendering because after a few minutes my GPU was so hot that my computer started to show garbage pixels in the screen, and then it was all broken and unreadable!!! So I reimplemented the thing in the CPU and started over again. For 1280x720 I'm getting around 4 seconds rendering in the CPU, with a 8 core CPU (so it should be like half minute in a monocore).

Anyway, this is the half video I rendered. I will re-render it in the CPU (don't wanna risk my HW again, and I don't want to start putting Sleep() every few frames to let the GPU take some breath) and upload to the same place, so you can try the link again by tomorrow or so for the full video.

I personally think that just rotation around the set is not as interesting as zomming in as in bib's video. But well, it's a start.

www.iquilezles.org/trastero/f8p5.avi

Excellent (though short) rotation anim !
As to using voxels I don't think anyone's tried it yet, I think other than raymarching the only other method tried was Buddhi who rendered many 2D slices producing full high resolution 3D data that could then be rendered quickly (to the given detail level) even with volumetric effects based on iteration values at any point.

[lycium]

zomg, that's quite a remarkable resemblance!

i'm a little bit concerned about what you said with gpu overheating; which gpu are you using? i haven't left my gpu rendering for extended periods because i was running into the watchdog timer limit of 2 seconds, but i'm about to break the computation into smaller pieces so that i can do arbitrarily long renders... actually i'm not too worried about breaking my 8800 gts - looking for an excuse to replace it - but if it takes my other system components with it, that would suck

[iq]

8600 gts here (a really low end gfx card) I think I should break the image into tiles too, otherwise I will have the same problem for huge posters like 100k x 100k pixels that will take 8 hours or so. Don't want kill my PC.

[twinbee]

Wow, it's all there - the eyes, nose, beard and hair!

Hi Jules, have you tried rendering large/raytraced images of the 3D set you found, and maybe higher powers too? If you know the boolean function (accepting x,y,z, and outputting true or false - i.e. in the set or not), it'd be interesting to see the output.