# How to make a Nebulabrot?

A Nebulabrot is a generalization of the Buddhabrot, a fractal rendering technique related to the Mandelbrot set that sort of looks like a meditating buddha. The Buddhabrot rendering technique was discovered and later described in a 1993 Usenet post to sci.fractals by Melinda Green.

I'm trying to generate this high-res image with Mathematica. Step 1: Implement a Buddhabrot.

The map is created by counting the number of times in the iterative creation algorithm a point is visited, it can look like this: Here is a great description of the algorithm.

Here's one old scrap of old implementation I found from Paul Nylander:

Mandelbrot[zc_] := Length[NestWhileList[#^2 + zc &, 0, Abs[#] < 2 &, 1, 1000]];
n = 275; image = Table[0, {n}, {n}];
Do[If[Mandelbrot[xc + I yc] < 1000, Module[{z = 0, iter = 0}, While[Abs[z] < 2, z = z^2 + xc + I yc; {i, j} = Floor[n{Im[z] + 1.5, Re[z] + 2}/3]; If[0 < i <= n && 0 < j <= n, image[[i, j]] += xc + I yc]; iter++]]], {xc, -2.0, 1.0, 0.01}, {yc, -1.5, 1.5, 0.01}];
Show[Graphics[RasterArray[Map[Hue[Arg[#]/Pi, 1, Min[Abs[#]/18, 1]] &, image, {2}]], AspectRatio -> 1]];


But this doesn't really do the trick, it generates this: Step 3: Use Image3D

A true three-dimensional Nebulabrot, like this is what I'm after: • related: Smooth Peter de Jong attractor
– Kuba
Jul 29 '15 at 18:32
• To get rid of the grid, you have to discard the fast-converging orbits. In the article that you linked, Michal Dichter says: "I found 297 points with iteration depths between 500k and one million." which suggests he only uses one half of all orbits. You can achieve this by adding the lower bound condition in the Do loop: If[ 500 <= Mandelbrot[xc + I yc] < 1000, ...
– shrx
Jul 31 '15 at 19:26
• Also, as Mathematica will warn you, the visualization code is deprecated and slow. The last line of Paul Nylander's code should be replaced with ColorConvert[ Image[Map[{Arg[#]/Pi + 1, Min[Abs[#]/18, 1]} &, image, {2}], ColorSpace -> "HSB"], "RGB"].
– shrx
Jul 31 '15 at 19:30
• I have a usable implementation written in C on Github. It produces the following output: ![nebulabrot](i.imgur.com/SAfD6.jpeg) Nov 14 at 22:44

Here is some code I used, which may partially answer your question. The key is to calculate efficiently points near the border of the Mandelbrot set. These points have large iteration counts which, when iterated, produce the Buddhabrot form.

The algorithm linked to in the question uses an adaptive mesh of squares to locate border points. The alternate code below was developed from pages 896 to 901 of Chaos and Fractals by Peitgen, Jürgens, and Saupe. Their technique follows directly the boundary of the Mandelbrot set with triangles of a given size, without adapting a mesh covering the entire region of interest. I am not sure of the comparitive efficiency of the two algorithms, but I like the chances of the triangulation approach.

A non-detailed outline of the algorithm follows.

Triangle vertices are labelled with 0 if they are outside the set, and with 1 if they are inside the set. Hence, each vertex is stored as {x,y,m}, where x and y are the real and imaginary coordinates of the vertex, and m=1 if the vertex is inside the set, and m=0 otherwise.

Define an initial triangle located on the horizontal x axis with sides dx. Increase the initial x=0.2 by dx until the point {x,0} is outside the set. Complete the initial triangle with {x,dx} labelled with its membership value.

InitialTriangle[dx_] :=
Block[{x = 0.2},
While[Quiet[MandelbrotSetMemberQ[x]], x += dx];
{{x - dx, 0.0, 1},
{x, 0.0, 0},
{x, dx, Boole[Quiet[MandelbrotSetMemberQ[x + I*dx]]]}}]


The new point defining the new triangle is znew = zinside + zoutside - zflip.

The function BorderTriangles[dx] constructs a chain of triangles, stopping when it reaches the starting triangle again.

BorderTriangles[dx_] :=
Block[{t, t1, zflip, zin, zout, znew},
t = InitialTriangle[dx];
t1 = Sort[Round[t, 10.^-10]];
If[t[[3, 3]] == 1,
zflip = t[]; zin = t[]; zout = t[],
zflip = t[]; zin = t[]; zout = t[]
];
t = {t}; (* enclose in braces for Join below *)
While[
znew = zin + zout - zflip;
If[Quiet[MandelbrotSetMemberQ[Most[znew].{1, I}]],
zflip = zin; zin = ReplacePart[znew, 3 -> 1], (* same zout *)
zflip = zout; zout = ReplacePart[znew, 3 -> 0] (* same zin *)
];
Sort[Round[{zflip, zin, zout}, 10.^-10]] =!= t1,
t = Join[t, {{zflip, zin, zout}}]
];
t]


Plot the line connecting all triangle centres. These centroids are good candidates for points with large iteration counts.

Graphics[Line[Map[Most[Mean[#]] &, BorderTriangles[0.01]]]] The function DeepPoints[dx,itmin,itmax] finds iteration points near the border with iteration counts between itmin and itmax.

DeepPoints[dx_, itmin_, itmax_] :=
With[{p = Map[Most[Mean[#]] &, BorderTriangles[dx]]},
Cases[
Transpose[{p,
MandelbrotSetIterationCount[p.{1, I}, MaxIterations->itmax]}],
_?(itmin <= Last[#] < itmax &)]
]


Store about 5000 points in the global variable borderpoints.

borderpoints =
Reverse[SortBy[Map[Flatten, N[DeepPoints[0.001, 1000, 500000]]], Last]];


Given a border point c with a large iteration count, the function iterates[c] calculates the list of iterates of this initial point.

iterates = Compile[{{c, _Real, 1}},
Map[{Re[#], Im[#]} &,
NestList[#^2 + {1, I}.c[[{1, 2}]] &, 0.0 + 0.0*I, Floor[c[] - 10]]],
CompilationTarget -> "C", Parallelization -> True,
RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed"];


Calculate an array of bin counts.

Block[{bp, bins, d, p, len, xmin = -1.5, xmax = 1.0, ymin = -1.1, ymax = 1.1,
delta = 0.005, tx, ty},
bp = borderpoints[[1;;-1;;5]];
bins = ConstantArray[0, Floor[{ymax - ymin, xmax - xmin}/delta] + {1, 1}];
d = Dimensions[bins];
Do[
p = Drop[iterates[bp[[j]]], 100];
len = Length[p];
Do[
tx = Floor[(p[[i, 1]] - xmin)/delta] + 1;
ty = Floor[(p[[i, 2]] + ymax)/delta] + 1;
If[tx >= 1 && tx <= d[] && ty >= 1 && ty <= d[], bins[[ty, tx]] += 1],
{i, 1, len}],
{j, 1, Length[bp]}];
Export["bins.dat", Transpose[bins]]]


Note the exponent of 0.2 which collapses the range of bin counts so that the colour function displays more detail.

ArrayPlot[
Import["bins.dat"],
ColorFunction -> (GrayLevel[#^0.2] &), ImageSize -> 1500] EDIT

M.R. wanted a colour example, so here is a Buddha Face from the central region. Take the log of the counts in the bins and using exponent 0.4 help make subtle changes visible.

With[{bins = Import["Bins1000_0010.dat"]},
ArrayPlot[Log[Log[bins + 1.] + 1.],
ColorFunction -> (ColorData["DarkRainbow", #^0.4] &),
ColorFunctionScaling -> True, ImageSize -> 800,
PlotRange -> {{700, 1500}, {400, 1800}}]] • +1 but ... :D Aug 4 '15 at 23:57
• The final gamma correction was a nice touch. :) Aug 5 '15 at 1:30
• ok, but show me the colors :)
– M.R.
Jul 12 '18 at 4:15
• @KennyColnago thanks for the update with color, but is there some way of getting the colors to have that glowing rainbow look like in the OP?
– M.R.
Jul 24 '18 at 1:24
• There’s no glow here, must be missing something big...?
– M.R.
Feb 4 '19 at 6:30