# Drawing Clifford Attractors in Mathematica

I came upon the "Clifford Attractor" today, which are equations of the form: $$x_{n+1}=\sin(a\ y_n) + c \cos(x_n)$$ $$y_{n+1}=\sin(b\ x_n) + d \cos(y_n)$$ They create some incredible images by choosing different values of $a,b,c,d$:

With a = -1.24458; b = -1.25191; c = -1.815908; d = -1.90866;, I've been trying to reproduce these in Mathematica using the following code:

tbl = RecurrenceTable[{x[n + 1] == Sin[a*y[n]] + c*Cos[a*x[n]],
y[n + 1] == Sin[b*x[n]] + d*Cos[a*y[n]], x[0] == 0,
y[0] == 0}, {x, y}, {n, 1, 25000}];
ListPlot[tbl, PlotRange -> All, Axes -> False, PlotStyle -> Black]


Which produces:

Obviously, there is a difference in both of these images. Mostly, the fact that in the first image, the color of each pixel is determined by how many times that point appears. I don't care that the scale of the image is different, but how can I at least approximate the "style" of the first image in Mathematica?

Parameters:

a = -1.24458; b = -1.25191; c = -1.815908; d = -1.90866;


Compiled function used for iteration:

cf = Compile[{{pt, _Real, 1}},
{Sin[a pt[[2]]] + c Cos[pt[[1]]],
Sin[b pt[[1]]] + d Cos[pt[[2]]]},
CompilationOptions -> {"InlineExternalDefinitions" -> True}
];


Iterate, rescale the result to the box {{0,1},{0,1}}, histogram it, convert to an image, and finally apply a gamma adjustment (^0.5) for better visibility.

im = ImageAdjust@Image@BinCounts[
Rescale@NestList[cf, {0., 0.}, 1000000],
1/500., 1/500.
];
ColorNegate[im^0.5]


You may want to throw in a Colorize too.

With 10 million points and $\gamma = 1/3$, I get

Update: I realized that by trying to show off too many features of LTemplate, I made an overengineered mess that will sooner deter people from LTemplate then attract them. Here's a single function solution, which is 95% as good as the complicated one below, but much shorter.

template =
LClass["CliffordAttractor2",
{LFun["compute", {{Real, 1}, Integer (* iterations *), Integer (*
image width *)}, Image]}
];

code = "
#include <cmath>

using namespace std;
using namespace mma;

struct CliffordAttractor2 {

GenericImageRef compute(RealTensorRef param, mint n, mint size) {
massert(param.size() == 4 && size > 0);

double a, b, c, d, h, w;
a = param[0]; b = param[1]; c = param[2]; d = param[3];
w = 1 + abs(c);
h = 1 + abs(d);

auto image = makeImage<float>(size, size * (h/w));
std::fill(image.begin(), image.end(), 0.0f);

double x = 0, y = 0;
for (mint i=0; i < n; ++i) {
double newx, newy;
newx = sin(a*y) + c*cos(a*x);
newy = sin(b*x) + d*cos(b*y);
x = newx; y = newy;
image( (image.rows()-1) * (y+h)/(2*h), (image.cols()-1) * (x+w)/(2*w) ) += 1;
}

return image;
}
};
";
Export["CliffordAttractor2.h", code, "String"];

CompileTemplate[template, "CompileOptions" -> {"-ffast-math"}]

cliff = Make[CliffordAttractor2];

im = cliff@"compute"[{-1.7, 1.3, -0.1, -1.2}, 100000000, 400];

ColorFunction -> (Blend[{White, RGBColor[
0.97, 0.9500000000000001, 0.79], RGBColor[
0.97, 0.25, 0.04]}, #] &)]


Here's a LibraryLink implementation with LTemplate.

Create class and initialize with the the desired $a,b,c,d$ parameters and image width.

clifford = Make[CliffordAttractor];
clifford@"init"[{-1.8, -2.0, -0.5, -0.9}, 600]


Compute 50 million iterations. If the image is not of sufficient quality, more iterations can be computed without losing the old data.

clifford@"compute"[50000000] // AbsoluteTiming
(* {3.1645, Null} *)


Visualize with a custom colour function:

Colorize[
ColorFunction -> (Blend[{White, RGBColor[0.87, 0.94, 1], RGBColor[0.48, 0.33333, 0.66667], Red}, #] &)
]


Check the current $(x,y)$ value:

clifford@"state"[]
(* {-0.773525, 1.1536} *)


The library code follows. Admittedly, this could have been done with a single function that takes the parameters and returns an image. To refine the result, we could have averaged multiple returned images.

Here I wanted to demonstrate how to maintain a state within the library and update it or retrieve information about it as needed.

Needs["LTemplate"]
SetDirectory[\$TemporaryDirectory];

template =
LClass["CliffordAttractor",
{LFun["init", {{Real, 1, "Constant"} (* {a,b,c,d} *), Integer (* image width *)}, "Void"],
LFun["setState", {{Real, 1, "Constant"} (* {x,y} *)}, "Void"],
LFun["state", {}, {Real, 1}], (* get {x,y} *)

LFun["compute", {Integer (* iterations *)}, "Void"],
LFun["image", {}, Image]}
];

code = "
using namespace mma;

class CliffordAttractor {
double a = 0.0, b = 0.0, c = 0.0, d = 0.0;
double x = 0.0, y = 0.0;

ImageRef<float> *im = nullptr;

double w, h; // image half-width and half-height in real coordinates

void free() {
if (im) {
im->free();
delete im;
}
}

public:
~CliffordAttractor() { free(); }

void init(RealTensorRef param, mint size) {
massert(param.size() == 4);

a = param[0]; b = param[1]; c = param[2]; d = param[3];
w = 1+std::abs(c);
h = 1+std::abs(d);

free();
im = new ImageRef<float>(makeImage<float>(size, std::ceil(size * (h/w))));
std::fill(im->begin(), im->end(), 0.0);
}

void setState(mma::RealTensorRef state) {
massert(state.size() == 2);
x = state[0]; y = state[1];
}

mma::RealTensorRef state() const { return mma::makeVector<double>({x,y}); }

void compute(mint n) {
massert(im);
for (mint i=0; i < n; ++i) {
double newx, newy;
newx = std::sin(a*y) + c*std::cos(a*x);
newy = std::sin(b*x) + d*std::cos(b*y);
x = newx; y = newy;
(*im)( (im->rows()-1) * (y+h)/(2*h), (im->cols()-1) * (x+w)/(2*w) ) += 1;
}
}

GenericImageRef image() const { massert(im); return im->clone(); }
};
";
Export["CliffordAttractor.h", code, "String"];

CompileTemplate[template, "CompileOptions" -> {"-O3 -ffast-math"}]


• Nice. I can't find another (web) reference with the OPs definition for the equations. I find Sin[a y] + c Cos[a x] and Sin[b x] + d Cos[b y]` (e.g. here.) Obviously doesn't change your code, but I wonder if one form or the other is 'correct'. (Although your last code block uses this definition.) – bobthechemist Sep 22 '18 at 17:38