# How can this image (optical illusion) be created with Mathematica?

I came across this image the other day:

and liked the sensation of it pulsing. I was wondering if anyone would know how to create something similar with Mathematica (without the Pink Floyd Dark Side of the Moon logo).

Edit

Very nice work from @halirutan, @Silvia, and @J.M -- you've given me headaches -- literally :)

-
They did it to me: what have you tried? :^) –  Mr.Wizard Apr 27 '13 at 7:16
@Mr.Wizard actually I am breaking protocols. Haven't been here for quite a while and saw this and wondered if it could be done -- or similar -- with Mma but had no idea about how I might start. Given I have not made any attempt it is ok to remove it. It is mainly a curiosity. Mike –  Mike Honeychurch Apr 27 '13 at 7:33
OH My eyes...$\phantom{}$ –  Silvia Apr 27 '13 at 7:39
Mike, I'm not going to remove it! You've contributed far more than enough to justify some no-effort-shown questions IMHO, and these questions can turn out to be some of the most interesting. I don't have an answer for you but you might get inspiration form reading this: intothecontinuum.tumblr.com –  Mr.Wizard Apr 27 '13 at 9:42
Also this site about all kinds of optical illusions, especially this example. –  Silvia Apr 27 '13 at 11:43

## How to make your eyes hurt

Mike asked whether it is possible to recreate the image he posted in his question. Although I haven't searched the web whether the equations for the above image are published somewhere, I will show how you can create such kind of image by pure inspection.

By inspecting Mike's original image, one recognizes the following things:

• the pattern is rotationally symmetric, which means to me that it's probably easier to recreate the pattern along a radius and then transform it to polar coordinates.
• the pattern is some kind of wave with sharp peaks and smooth bottom.
• in addition to the (repeated) coloring of the pattern itself, we see a change of color when going radially outwards.
• when you go radially outward, you see that the repetition of the pattern slows down.

Using this information, you could first try to reproduce the wave pattern. Here, I started with a squared sine function and used Log to make sharp peaks:

With[{n = 3},
Plot[-Log[Sin[n x]^2 + 1/100]/Log[100], {x, 0, 2 Pi}]
]


This is the basic idea. What's left is to transform this into heights in polar space and the inclusion of the other phenomena. This is basically playing with sine functions of different frequencies. The two most important parameters are probably the number of divisions and repetitions.

Manipulate[
Plot3D[
With[{r = Norm[{x, y}], phi = ArcTan[x, y]},
Sin[division2*
Exp[-r/50] (r + (1/2 + 1/4 r) Log[Sin[division1 phi]^2 + 1/100]/
Log[100])]^2]
, {x, -10, 10}, {y, -10, 10},
PlotStyle -> ControlActive[None, Automatic],
Mesh -> ControlActive[Full, Automatic], PlotPoints -> 50,
ColorFunction -> "Rainbow"],
{{division1, 3}, 2, 4},
{{division2, 1}, 0.1, 2}
]


The full code which rasterizes the plane and creates the color image is the following:

f = Compile[{{x, _Real, 0}, {y, _Real, 0}},
Module[{phi, r = Norm[{x, y}]},
If[r == 0, 0,
phi = ArcTan[x, y];
1/10 Sin[r] + Sin[
3 Exp[-r/50] (r + (1/2 + 1/4 r) Log[Sin[15 phi]^2 + 1/100]/Log[100])]^2]
],
CompilationTarget -> "C", Parallelization -> True,
RuntimeAttributes -> {Listable}
];

With[{imageSize = 512},
pts = Table[p, {p, -20, 20, 40/(imageSize - 1)}];
]

img = Image[Rescale@Outer[f, pts, pts]]


The coloring can be done with

Colorize[img, ColorFunction -> "SunsetColors"]

-
How did you do that? –  bill s Apr 27 '13 at 10:38
+1 Would like to see the explanation very much! –  Silvia Apr 27 '13 at 11:40
So it's the Log! Looks like I should not limit my mind with polynomial.. –  Silvia Apr 27 '13 at 11:53
Maybe you can use Colorize[] instead of ImageApply[] for the coloring... –  Guess who it is. Apr 27 '13 at 14:21
@J.M. Yes, this is definitely faster. –  halirutan Apr 27 '13 at 14:24

The best I can do..

Well I think I missed some important properties on optical illusion that is presented in OP and halirutan's answer. I would very much like to wait for halirutan's explanation on it (if he is willing :)

Here is how I did it.

The outline shape is governed by equation 50 (1 + r/2) (Abs[Mod[θ, (2 π)/50] - π/50]^2 + r^0.1 10^-2) with 0 < θ < 2 π and 0.1 < r < 14, and then render them and add foreground objects.

lineFunc = Compile[{{z, _Complex}},
Module[{r, θ},
{r, θ} = {Abs[z], Arg[z]};
50 (1 + r/2) (Abs[Mod[θ, (2 π)/50] - π/50]^2 + r^.1  10^-2)
],
CompilationTarget -> "C",
CompilationOptions -> {"ExpressionOptimization" -> True},
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"]

θrange = Range[0, 2 π, .01];
rrange = Reverse@Range[.1, 14, .02];

ρset = lineFunc[Flatten[ Table[r Exp[I θ], {r, rrange}, {θ, θrange}]]];

lineSet = Function[θrange,
Function[{ρrange}, (ρrange # & /@ Through[{Cos, Sin}[θrange]])\[Transpose]] /@
Partition[ρset, Length[θrange]]]@θrange;

basegraph = Graphics[{
MapIndexed[
{ColorData["SunsetColors"][
Rescale[(Mod[Rescale[#2[[1]], {1, Length[rrange]}, {500, 0}],
30]/30)^1.3, {0, 1}, {0, 1}]
], Polygon[#1]} &,
lineSet]
}, PlotRange -> 3 {{-1, 1}, {-1, 1}}]

foreground = Graphics[{
EdgeForm[{Lighter[Yellow, .8], Thin}],
FaceForm[Black], Polygon[.5 {{-2, -1}, {2, -1}, {0, 1.6}}],
Lighter[Yellow, .8], Thickness[.005],
Line[{{-3, -0.5311}, {-0.4623, 0.1967}}],
EdgeForm[],
Polygon[{{-0.4623, 0.1967}, {0.5508, 0.06885}, {0.3148, 0.3738}},
VertexColors -> {Lighter[Yellow, .8], Black, Black}],
Function[{pt00, pt01, pt10, pt11, n},
Module[{divset, pts0, pts1, divsetlen, resortset},
divset = Range[0, 1, 1/n];
divsetlen = Length[divset];
pts0 = # pt00 + (1 - #) pt01 & /@ divset;
pts1 = # pt10 + (1 - #) pt11 & /@ divset;
resortset = {1/Sqrt[2] Cos[# π/2 + (3 π)/4] + 3/2 & /@
Range[2 #],
Flatten[{#, #}\[Transpose]] &@Range[#]}\[Transpose] &[
divsetlen];
MapIndexed[{EdgeForm[
ColorData["Rainbow"][Rescale[#2[[1]], {1, divsetlen}]]],
FaceForm[
ColorData["Rainbow"][Rescale[#2[[1]], {1, divsetlen}]]],
Polygon[#1]} &,
Partition[{pts0, pts1}[[##]] & @@ # & /@ resortset, 4, 2]]
]][{0.5508, 0.06885}, {0.3148, 0.3738}, {3, -.3}, {3, .3}, 10]
}]

Show[basegraph, foreground]

-
If that's the best you can do, then I guess I'll just have to vote for that. :) –  Mark McClure Apr 27 '13 at 14:30
@MarkMcClure Thanks! I tried so hard just can't find the right equation to make the picture "twinkle".. –  Silvia Apr 27 '13 at 16:31
"Twinkle" is one way to put it, but "sea sickness" is another - but in +1 seasick kind of way :-) –  Yves Klett May 2 '13 at 14:04

After trying things out, it would seem that I had ended up with a solution qualitatively similar to halirutan's; that is, wrap a sharply-peaked wave around a circle, and then have the "radius" vary. Nevertheless, my choice of the wave used is different, so I thought that I might as well post my variation:

With[{b = 10, f = 30, h = 2, w = 9},
DensityPlot[Sin[w (Norm[{u, v}] - (h +
Cos[f ArcTan[u, v]]^2/(1 + b Sin[f ArcTan[u, v]]^2)))],
{u, -2 π, 2 π}, {v, -2 π, 2 π}, ColorFunction -> "Heat", Frame -> False,
MaxRecursion -> 1, PlotPoints -> 405]]


Due to the limitations of the machine I am now using, I had to limit the settings of MaxRecursion and PlotPoints; this shows in the "jitter" that shows up near the image's center. A higher setting of these ought to give a better image.

Here's a slightly simpler variation:

With[{f = 30, h = 3, w = 9},
DensityPlot[Sin[w (Norm[{u, v}] - (h - Abs[Sin[f ArcTan[u, v]]]))],
{u, -2 π, 2 π}, {v, -2 π, 2 π}, ColorFunction -> "Heat",
Frame -> False, PlotPoints -> 305]]


-

The only significant difference is compositing rings of colors over the wavy pattern.

  (* colors for color mask *)
colFn = Blend[{{0, RGBColor[1, 0, 0.6]}, {0.03,
RGBColor[1, 0, 0.6]}, {0.04, RGBColor[0.4, 0, 0.8]}, {0.06,
RGBColor[0.5, 1, 0.8]}, {0.08, Lighter[Yellow, 0.6]}, {0.14,
Lighter[Yellow, 0.6]}, {0.2, RGBColor[0.2, 0, 0.9]}, {0.26,
RGBColor[0.8, 0, 0.2]}, {0.32, Lighter[Yellow, 0.5]}, {0.4,
RGBColor[0.8, 0, 0.2]}, {0.46, RGBColor[0.2, 0, 0.8]}, {0.5,
RGBColor[0.8, 0.8, 0.5]}, {0.54, RGBColor[0.8, 0, 0.2]}, {0.58,
Darker@Red}, {0.62, RGBColor[0.2, 0, 0.8]}, {0.66,
RGBColor[0.8, 0.8, 0.5]}, {0.7, Darker@Red}, {0.8,
Darker@Blue}, {0.9, Yellow}, {1.0, Darker@Red}}, #] &;
ContourPlot[Identity[x^2 + y^2], {x, -3, 3}, {y, -3, 3},
Contours -> 60, ContourStyle -> None,
ColorFunction -> (Directive[Opacity[0.75],
colFn[Mean[{#, Round[#, 1/50]}]]] &), Frame -> None,

(* for the wavy pattern *)
f = Compile[{{x, _Real}, {y, _Real}, {n, _Integer}},
Module[{r, t, dx, dy},
dx = (x - 1)/(n - 1) - 1/2;
dy = (y - 1)/(n - 1) - 1/2;
r = Sqrt[dx^2 + dy^2];
t = If[dx == 0 && dy == 0, 0, ArcTan[dx, dy]];
2 r^(1/4) - Sqrt[r] (6 Mod[t, 2 π/49, -π/49])^2]];

(* colors for color function of wavy pattern *)
bl = Darker[Blue, 0.7]; rd = Darker[Red, 0.2]; yl = Lighter[Yellow, 0.7];
rays = Image@Array[
List @@ Blend[{bl, bl, bl, bl, rd, rd, Orange, yl, yl, yl, bl},
Mod[20 #, 1]] &@f[##, 600] &, {600, 600}];

finalImg = ImageCompose[rays, {colormask /. _EdgeForm -> EdgeForm[Opacity[0.]], 0.35}]


-
+1 for the fine color tuning. –  Silvia Apr 29 '13 at 12:51
nice attention to detail! –  Sum it May 4 '13 at 18:04
This one has the strongest effect for me. It appears the colors are important. +1 –  shrx Mar 6 '14 at 23:51

I think all the other answers do a better job at exactly replicating the original image than what I am going to share, but my main intention here is to provide some exposition and show the utility in a particular coordinate transformation that naturally results in graphics having similar properties to the original image. I will refer to this transformation as a log-polar transform (also referred to as log-polar coordinates) for reasons which will become clear after its definition is given below.

Interestingly enough, what makes this transformation "natural" and yield psychedelic characteristics is its relationship with the anatomical properties of the human eye and its neurological basis in explaining the various form constants perceived during visual hallucinations. To the best of my knowledge, the earliest account of such a mathematical modelling in the literature seems to be the 1979 paper, by J. D. Cowan and G. B, Ermentrout, “A Mathematical theory of Visual Hallucinations”.

For some motivation, consider some image in the complex plane with coordinates given in polar form as:

$z=re^{i\theta}$.

Taking the complex natural logarithm of $z$ gives:

$\ln(z)=\ln(r)+i\theta$,

which is now expressed in standard form. Here, the real part is the logarithm of the radial component of $z$ and the imaginary part is just the angular component of $z$.

The log-polar transform, in the context of the complex plane, is just the mapping which results from taking the complex logarithm of each of the points in the plane.

Instead, in the context of the Cartesian plane, the log-polar transform can be though of as the mapping which takes points $(x,y)=(r\cos(\theta),r\sin(\theta))$ to the points $(x',y')=(\ln(r),\theta))$, or more explicitly as:

$(x',y')=(\ln(\sqrt{x^2+y^2}),\mathrm{atan2}(y,x))$.

Some properties of this conformal mapping:

• vertical lines turn into circles (constant radius)
• horizontal lines turn into radial rays (constant angle)
• lines at other angles spiral out from the origin

As an illustration, consider the following periodic density plot:

    img =
ImageCrop@DensityPlot[
Sin[2 x - 20 Log[2 (Sin[y]^2 + 1), 2]],
{x, 0, 16 Pi}, {y, 0, 32 Pi},
PlotPoints -> 250, ColorFunction -> "SunsetColors",
Frame -> False, ImageSize -> 600]


To apply the log-polar transform to this image, first define the map:

LogPolar[x_, y_] := {Log[Sqrt[x^2 + y^2]], ArcTan[x, y]}


Then use Mathematica's ImageTransformation command on the original image:

 ImageTransformation[img, LogPolar[#[[1]], #[[2]]] &,
DataRange -> {{-Pi, Pi}, {-Pi, Pi}}]


Note: in order for the transformed image to appear seamless at the angle corresponding to $\pi$ radians, the top and bottom edges of the original image should appear seamless if joined together.

We can exploit the translational symmetry of the original plot to create a zooming animation after the log-polar transform has been applied. Instead of having to recompute the plot for every frame of the animation, lets use ImageTake to crop a portion of the original image and then shift this crop vertically by an amount that corresponds to the periodicity of the plot:

d = ImageDimensions[img][[1]]

Export["LPTzoom.gif",
Table[
ImageResize[
ImageTransformation[
ImageTake[
img,
{1, 14*d/16}, {1 + (2 - 2 t)*d/32, (32 - 2 t)*d/32}],
LogPolar[#[[1]], #[[2]]] &, DataRange -> {{-Pi, Pi}, {-Pi, Pi}}],
500],
{t, 0, 6/7, 1/7}]
]


Similarly, translating the original image horizontally would produce a spinning animations instead of a zooming one. For good measure, combining both of these two directions of motion results in a spiraling animation:

ProTip: try looking at the still image after staring at the animation for a little motion aftereffect.

The interested viewer is invited to explore log-polar transforms of various images in excess at this link.

-
Neat paper, that one. Thanks for linking to it! –  Guess who it is. May 4 '13 at 18:16
Good answer! And thanks for the references! –  Silvia May 4 '13 at 21:41
Really love the linked blog! Thanks! –  Silvia May 4 '13 at 22:05