So I discovered some links created by others about a 3D kind of fractal that is sufficiently interesting (a 3D mandelbrot rather than a Julia).
Then this link came after I had done this work! Slashdotted
So I went ahead an implemented myself just so that I could explore it.
Haskell Voxel Implementation
I felt like playing around with parallelism in haskell so I wrote this and stole the mapP function (seriously though I was going to write exactly that anyways). I wanted to play with parallelism and see it. So I decided I'd rather generate voxels and plot those. The end result was this source code:
I had a heck of a time in Paraview partially because the concepts of the datafiles were no clear. I had started by printing out points where the object existed. The format for this is a little complicated, but I couldn't get volume rendering! I was able to convert points to glyphs though:
Eventually I figured out that it was best to say that my data is voxels and the points that aren't important but here are the scalar values at the points. This was using structured_points vtk data type where the grids is assumed and then you stream out all of the scalars associated with each point. But at least it worked!
An example of volume rendering. See more at http://softwareprocess.es/x/x/mandelbulb/showoff1/
Here are the VTK files I made BTW since the haskell is horribly slow and I didn't parallelize it right (the scheduling overheard was too great and the parallelism was too granular, IMHO).
- 100 x 100 x 100 voxel of 1/8th of the mandelbulb 95kb http://softwareprocess.es/x/x/mandelbulb/bigoutx100x100x100.vtk.bz2
- 300 x 300 x 300 voxels 2MB http://softwareprocess.es/x/x/mandelbulb/bigoutx300x300x300.vtk.bz2
- 1000 x 1000 x 1000 voxels 50MB http://softwareprocess.es/x/x/mandelbulb/bigoutx1000x1000x1000.vtk.bz2
- Here's the haskell function I used.
-- iterate over this! mandebulbf n (x,y,z) (xo,yo,zo) = (nxc,nyc,nzc) -- (newx,newy,newz) where r = sqrt ( x^2 + y^2 + z^2 ) theta = atan2 (sqrt (x^2 + y^2)) z phi = atan2 y x rn = r**n thetan = theta * n phin = phi * n sinthetan = sin thetan newx = rn * sinthetan * (cos phin) newy = rn * sinthetan * (sin phin) newz = rn * (cos thetan) nxc = newx + xo nyc = newy + yo nzc = newz + zo
Here's a video of using Paraview if you want to see more angles: http://softwareprocess.es/x/x/mandelbulb/paraview-mandelbulb-demo.ogv
I decided raytracing was necessary (in fact a friend kept harassing me about it). So I modified Povray because I have previous POV experience, I like POVray, and it is a quality render.
I found this tutorial and made my own changes: http://www.abx.art.pl/pov/patches/f_triangle.php
I found that values less than the threshold were drawn by the isosurface (so 0.0 or less). Which I implemented myself to solve the Mandelbulb problem.
- Modified povray 3.6.1 source/fnintern.cpp file:
- Patch http://softwareprocess.es/x/x/mandelbulb/fnintern.cpp.patch
- My modified copy of the file http://softwareprocess.es/x/x/mandelbulb/fnintern.cpp
- Povray 3.7-beta34 Patch (Created and Optimized by Folkert van Heusden)
- If you need to be stingy about it, I put my mandelbulb changes under the WTFPL V2 with boilerplate, this means POVRAY can integrate with it just fine.
- Note due to the distribution rules of POVRAY you must apply this patch yourself
- Povray scene used to model the mandelbulb with an ISOSurface http://softwareprocess.es/x/x/mandelbulb/mandelbulb.pov
First run without Antialiasing. 1920x1200 version
With all of POVRAY's bell's and whistles on. 1920x1200 version
Later I considered that if you were only going to do so many iterations that povray's own internal language would've sufficed as it'd just be loop unrolling.
I really suggest looking closely at the hi-res anti-aliased image. Investigate the shadows. You'll see that solid structures are actually quite hollow and lacey as shown in their shadows.
I'd like to improve the rendering or animate but I simply don't have time.
Perl Mongers Talk
- During the Perl monger talk I explained some of the math behind mandelbrot fractals by showing the actual code: