# Smooth Peter de Jong attractor

Today I was playing with Peter de Jong attractor.

At the bottom of the page I've linked there are beautiful examples like:

My attempts are not so great:

It is around 10^5 points. For more than 5*10^5 my 4GB RAM gives up. How can I achieve such smooth result on standard/oldish pc?

Here is the piece of code to play with:

Points generator:

fr2 = Compile[{{p, _Real, 1}, {n, _Integer}, {a, _Real, 1}, {b, _Real, 1},
{c, _Real, 1}},
NestList[
{ Sin[a[[1]] #[[2]]] - Cos[a[[2]] #[[1]]],
Sin[b[[1]] #[[1]]] - Cos[b[[2]] #[[2]]],
Sin[c[[1]] #[[1]]] - Cos[c[[2]] #[[1]]]
} &, p, n]];


Interactive toy:

For printed color image I've used approach introduced by Szabolcs in antialiasing3D.

DynamicModule[{a, b, c, n, type, fig, controls, print},
Deploy@Dynamic[Refresh[
Panel@Grid[{{fig, controls}}, Spacings -> 2],
None]],

Initialization :> (
type = 1;
n = 100000;
{a, b, c} = {{1.4, -2.3}, {2.4, -2.1}, {2.41, 1.64}};

controls = Column[{
Slider2D[Dynamic@a, {-Pi, Pi, .01}], Dynamic@a,
Slider2D[Dynamic@b, {-Pi, Pi, .01}], Dynamic@b,
Slider2D[Dynamic@c, {-Pi, Pi, .01}], Dynamic@c,
Button["Print", print[], Method -> "Queued"] }];

fig = Graphics[{White, AbsolutePointSize@1,
Dynamic[
Point@fr2[{.0, .0, .0}, ControlActive[2 10^4, 10^5], a, b, c][[All , ;; 2]]]
}, ImageSize -> {1, 1} 500, Background -> Black, AspectRatio -> Automatic];

print[] := With[{t = 3, pointsize = 1, pts = 10^5, res = 72},
Composition[
CreateDocument,
ImageResize[Rasterize[#, "Image", ImageResolution -> t res], Scaled[1/t]] &,
Graphics[{AbsolutePointSize@pointsize,
Riffle[Hue@Rescale[#, {-2, 2}, {0, 1}] & /@ #[[;; , 3]],
Point /@ #[[;; , ;; 2]]] &@#
}, ImageSize -> 800, Background -> Black] &
][fr2[{.0, .0, .0}, pts, a, b, c]]];
)]


-
Just related, here is a demonstration project from @vitaliy-kaurov – Murta Jun 15 '14 at 14:12
The first image you showed is made with 12 million points. They do not keep all the points in memory; whenever each point is computed they just set the corresponding pixel location in the image and move on. This way they don't need gigabytes of storage. – Rahul Jun 15 '14 at 14:19
Google "Flame fractals". I create this type of fractals a lot; 10⁵ is not nearly enough iterations. You would like about 2000 iterations PER PIXEL in the final image to get a decent result. – Per Alexandersson Jun 16 '14 at 6:55
@Paxinum Well, I understand :) Wanna share with us some fancy code/results? :) – Kuba Jun 16 '14 at 7:15
@Kuba: The implementation I have is written in Java, and is open-sourced here: sourceforge.net/projects/flamethyst Some images created can be found here: people.su.se/~peal0658/index.php?page=fractals – Per Alexandersson Jun 16 '14 at 8:02

UPDATE I thought it would be neat to try and animate the thing, so I let the $a$ parameter run between $-\pi$ and $\pi$. I generated 600 images and put them together using ffmpeg. Check it out on youtube.

It might not be in the spirit of Mathematica Stack Exchange, but allow me an objection - stuff that is slow in Mathematica should be kept out of it. To wit consider how this little C++ nugget does the grunt work:

#include <stdio.h>
#include <cmath>
#include <omp.h>

int main()
{
const int dim = 4096;
const float a = 1.4f, b = -2.3f, c = 2.4f, d = -2.1f;

int size = dim*dim;
float *image = new float[size];
for (int i = 0; i < size; ++i) image[i] = 1;

#pragma omp parallel
{
float x = omp_get_thread_num(), y = 0;
for (int i = 0; i < 10000000; ++i)
{
float xn = sin(a * y) - cos(b * x);
y = sin(c * x) - cos(d * y);
x = xn;

auto xp = ((dim - 1) * (1 + x * 0.43) * 0.5);
auto yp = (int)((dim - 1) * (1 - y * 0.43) * 0.5);
image[(int)((yp * dim + xp))] *= 0.99f;
}
}

FILE *file = fopen("image.bin", "wb");
fwrite(image, sizeof(float), size, file);
fclose(file);
delete[] image;
return 0;
}


This should be compiled with fast math for better performance. I've also included omp in there, because my system has 12 cores and if I don't use them for this - there's no justification for me buying it.

And this Mathematica code makes an image and colorizes it from the produced data:

buffer = BinaryReadList["image.bin", "Real32"];
dim = Sqrt[Length@buffer];
bigimg = Image[Partition[buffer, dim], "Real32"];
Colorize[Rasterize[bigimg, ImageSize -> dim/4],
ColorFunction -> ColorData["SunsetColors"]]


Note that I'm rendering the original on the 4096x4096 canvas, and then down-sampling it in Mathematica - which I find produces a more pleasant aesthetic. This stuff would have taken several days to do proper in C++, and while I'm sure it's possible to write a fast iterator in Mathematica - it would probably take a long time as well.

Final image:

-
Nice! BTW, this approach (using Java/C/C++ with Mathematica) is quite common on this site and even encouraged when appropriate (as it is in this case). Welcome to the site (even though you've had an account for ~2 years) :) – R. M. Jun 17 '14 at 17:24
@rm-rf Thanks. I used to be quite a Mathematica buff back in the day when there was no separate site. ;) – Gleno Jun 17 '14 at 17:38
Öska, that error message means that you are not running the latest compiler; try adding -std=c++0x compilation flag. Then, since I was greedy and added openmp, you should add -fopenmp flag as well (assuming GCC). I compiled and ran this on windows under latest visual studio. – Gleno Jun 17 '14 at 19:46
That's great you've decided to post that answer :) – Kuba Jun 17 '14 at 20:02
@Gleno Works like a charm with g++ -std=c++0x -fopenmp file.cpp indeed, thanks. – Öskå Jun 18 '14 at 9:21

As @RahulNarain says, forming the image point by point saves significant memory because the number of image pixels is typically much smaller than the hundreds of millions of iterations that compose it. Therefore, iterate the attractor equations, and for each point generated, find its location within the image matrix. Colour coding of the number of hits in each pixel makes the image.

The following is a rough example.

{a, b, c} = {{1.4, -2.3}, {2.4, -2.1}, {2.41, 1.64}};
Block[{xmin=-2., xmax=2., ymin=-2., ymax=2., delta=0.01, bins, d,
itmax=10^5, x, y, tx, ty},
bins = ConstantArray[0, Floor[{xmax-xmin, ymax-ymin}/delta] + {1, 1}];
d = Dimensions[bins];
{x, y} = {0., 0.};
Do[
{x, y} = {Sin[a[[1]] y] - Cos[a[[2]] x],
Sin[b[[1]] x] - Cos[b[[2]] y]};
tx = Floor[(x - xmin)/delta] + 1;
ty = Floor[(y - ymin)/delta] + 1;
If[tx >= 1 && tx <= d[[1]] && ty >= 1 && ty <= d[[2]],
bins[[tx, ty]] += 1],
{i, 1, itmax}];
ArrayPlot[Log[bins+1], ColorFunction->(ColorData["DarkRainbow",#^0.4]&)]]


This block may be compiled in the usual way. Other plotting routines such as MatrixPlot or ListDensityPlot may be used. TheLog[bins+1]function and the exponent within theColorFunctionare example methods to adjust the dynamic range of pixel counts to bring out subtle features.

COMPILED CODE

attractor =
Compile[{{xmin,_Real}, {xmax,_Real}, {ymin,_Real}, {ymax,_Real},
{delta,_Real}, {itmax,_Integer}, {a,_Real,1}, {b,_Real,1}},
Block[{bins, d, x, y, tx, ty},
bins = ConstantArray[0, Floor[{xmax-xmin, ymax-ymin}/delta] + {1,1}];
d = Dimensions[bins];
{x, y} = {0., 0.};
Do[
{x, y} = {Sin[a[[1]] y] - Cos[a[[2]] x],
Sin[b[[1]] x] - Cos[b[[2]] y]};
tx = Floor[(x - xmin)/delta] + 1;
ty = Floor[(y - ymin)/delta] + 1;
If[tx >= 1 && tx <= d[[1]] && ty >= 1 && ty <= d[[2]],
bins[[tx, ty]] += 1],
{i, 1, itmax}];
bins],
CompilationTarget :> "C"];


Example run:

AbsoluteTiming[bins = N[attractor[-2., 2., -2.2, 2.2, 0.005, 5*10^6,
{1.2, -2.1}, {2.4, -2.1}]];]
ArrayPlot[Log[bins+1], ColorFunction->(ColorData["FallColors",#^0.4] &)]


EDIT

I compiled the above code for speed. I also wrote a Mathematica interface to Fortran code which allows me to rapidly calculate hundreds of millions of iterations. That way aManipulatecan quickly find the parameters giving an interesting form. The resulting matrix of number of hits in each pixel may be converted to an image via the built-in colour functions, opacity, etc.

This image is formed with random affine transforms in the hyperbolic plane.

The following is a very simple 5-fold icon from the Symmetry In Chaos reference.

But the really interesting (for me) images are made by converting the number of hits in each pixel to spheres with radii (non-linearly) proportional to the number of hits. The POVRay ray tracer (or other similar tool) "blobs" these spheres together, then textures may be applied. I have not as much experience with the Mathematica textures.

My Stack Exchange icon is derived from the Symmetry In Chaos equation, then inverted in its bounding circle to have matching inner and outer loops.

Other manipulations include a Moebius transform to open up the generally circular symmetry into a half plane. Cuttle fish eyes added for effect...

-
– R. M. Jun 15 '14 at 21:17
Here's an animation where $a_1 \in [-2,2]$. – wxffles Jun 16 '14 at 0:55
@KennyColnago Great reference, thanks :) p.s. don't you know that without pictures posts get half the upvotes they can? :p – Kuba Jun 16 '14 at 6:00
I think it'll be better if you add the compiled code since the existed piece of code can't be compiled directly. Compile isn't easy for everyone :) – xzczd Jun 18 '14 at 7:46
compiled code example has been included – KennyColnago Jun 19 '14 at 15:09

I revisited this problem - this time in pure Mathematica. The trick to any kind of performance is the Compile[] function, which in itself can be a bit moody - so you need to set global options to warn you when it refuses compilation and work around that. The performance I'm seeing is on the order of magnitude slower than that I get from C++, and two orders of magnitude slower than I get with C++ AMP (GPGPU on Windows; I used AMP to generate images for this animation).

Another point of note, is that I've discovered that the algorithm starting points actually matter! If you take a uniform sample in the region $[-2, 2]$ as starting points, instead of $(x_0, y_0) = (0, 0)$ you get a slightly different image, typically with more artifacts. As I've tried to animate the $a$ parameter of the equations going from $-\pi$ to $\pi$, starting with $(x_0, y_0) = (0, 0)$ produces random blank frames for some values of $a$, while selecting points uniformly produces a smooth animation.

Anyway, the Mathematica code for the $(x_0, y_0) = (0, 0)$ case is as follows:

func = Compile[{{dim, _Integer}, {params, _Real,
1}, {iters, _Integer}},
Module[{matrix, d1, d2, a, b, c, d, x, y, iter, xp, xn, yp, yn,
max},
matrix = Table[0.0, {dim dim}];
{a, b, c, d, x, y, iter} = params~Join~{0, 0, 1};
{d1, d2} = {1 + 0.5 (dim - 1), 0.25 (dim - 1)};
While[++iter < iters,
xn = Sin[a y] - Cos[b x]; y = Sin[c x] - Cos[d y]; x = xn;
xp = Floor[d1 + d2 x]; yp = Floor[d1 - d2 y];
++matrix[[yp dim + xp]];
];
max = Max@matrix;
Partition[(Sqrt[#]/max &) /@ matrix, dim]
], CompilationTarget -> "C", RuntimeOptions -> "Speed"];

dejong[dim_, params_, iters_] := Module[{matrix},
matrix = func[dim, params, iters];
Colorize[ArrayPlot[matrix, ImageSize -> dim, Frame -> False],
ColorFunction -> ColorData["SunsetColors"]]];


You can get a value matrix by calling func, or an image by calling dejong. Example would be:

dejong[1024, {1.4, -2.3, 2.4, -2.1}, 20000000]


Which in roughly 2 seconds produces:

Faster output can be generated in black and white, if you forgo the colorization step, which takes roughly half the time.

-
Nice to mention the artifacts, I've enconutered them as well when I wrote my code to generate such images and was surprised when I couldn't see them in other images. – shrx Jul 6 '14 at 21:06
I'd like to mention that, the "CompileOptions" set in \$CCompiler will significantly influence the speed of the code in this case. For GCC, Mathematica uses "CompileOptions"->"-O2" by default, but when I changed it to "CompileOptions"->"-Ofast", I got a 2X speed-up! BTW, with my dual core old machine, your C++ code is even a little slower than the Mathematica code with -Ofast. – xzczd Oct 15 '14 at 10:08