# 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? :^) Apr 27, 2013 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 Apr 27, 2013 at 7:33
• OH My eyes...$\phantom{}$ Apr 27, 2013 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 Apr 27, 2013 at 9:42
• Also this site about all kinds of optical illusions, especially this example. Apr 27, 2013 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, {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)]^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)]^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? Apr 27, 2013 at 10:38
• +1 Would like to see the explanation very much! Apr 27, 2013 at 11:40
• So it's the Log! Looks like I should not limit my mind with polynomial.. Apr 27, 2013 at 11:53
• Maybe you can use Colorize[] instead of ImageApply[] for the coloring... Apr 27, 2013 at 14:21
• @J.M. Yes, this is definitely faster. Apr 27, 2013 at 14:24

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[#[], #[]] &,
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][]

Export["LPTzoom.gif",
Table[
ImageResize[
ImageTransformation[
ImageTake[
img,
{1, 14*d/16}, {1 + (2 - 2 t)*d/32, (32 - 2 t)*d/32}],
LogPolar[#[], #[]] &, 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! May 4, 2013 at 18:16
• Good answer! And thanks for the references! May 4, 2013 at 21:41
• Really love the linked blog! Thanks! May 4, 2013 at 22:05
• That last animation could be modified easily to simulate what the scintillation from a visual migraine looks like. Dec 20, 2020 at 18:15

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, 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 Cos[# π/2 + (3 π)/4] + 3/2 & /@
Range[2 #],
Flatten[{#, #}\[Transpose]] &@Range[#]}\[Transpose] &[
divsetlen];
MapIndexed[{EdgeForm[
ColorData["Rainbow"][Rescale[#2[], {1, divsetlen}]]],
FaceForm[
ColorData["Rainbow"][Rescale[#2[], {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. :) Apr 27, 2013 at 14:30
• @MarkMcClure Thanks! I tried so hard just can't find the right equation to make the picture "twinkle".. Apr 27, 2013 at 16:31
• "Twinkle" is one way to put it, but "sea sickness" is another - but in +1 seasick kind of way :-) May 2, 2013 at 14:04

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}] • nice attention to detail! May 4, 2013 at 18:04
• This one has the strongest effect for me. It appears the colors are important. +1
– shrx
Mar 6, 2014 at 23:51
• I think this is the best one. (And yes, most of the other presented solutions are great too.) Jun 10, 2016 at 12:01
• finalImg = ImageCompose[rays, {RemoveAlphaChannel[colormask /. _EdgeForm -> EdgeForm[Opacity[0.]]], 0.35}] works for me. Nov 21, 2020 at 8:44

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",
Exclusions -> None, Frame -> False, PlotPoints -> 505]] 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",
Exclusions -> None, Frame -> False, PlotPoints -> 505]] LogPolar[x_, y_] := {Log[Sqrt[x^2 + y^2]], ArcTan[x, y]}

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


Juts for clarification.

I am allmost shure that the above doese the inverse Log Polar Transform. To actualy do the transform one would have to use ImageForwardTransform[]. Alternatively

ImageTransformation[img, Through[{Cos, Sin}[#[]]]*Exp[#[]] &,
DataRange -> {{-1, 1}, {-1, 1}},
PlotRange -> {Pi {-1, 1}, Automatic}]


Of course this will not give the desired star/flower like pattern.

By the way

The transform functions for ImageTransformation[] can be provided more elegantly, by using Through[] to evaluate two functions on one argument and Apply[] to evaluate a (non list compatible) function on arguments given as a list

Log Polar:

Apply[{ϕ, r} \[Function] Exp@r*Through[{Cos, Sin}@ϕ]]


Inverse Log Polar:

(p \[Function] Through[{Log@*Norm, Apply@ArcTan}@p])


The expressions need to be copied to a MMA notebook to see the beauty :-)

Robert

• But Robert, how does one generate img to begin with? Otherwise your answer looks rather diminished compared to the other answers here. Jun 10, 2016 at 9:30
• Hi J. M. my answer was not intended to create the pattern/image. It was to clarify how to do a Log Polar Transform with MMA. I should have post this as a comment, but stackexchange did not allow me to comment - so I posted it as an answer. Robert Jun 10, 2016 at 10:09
• Since this doesn't attempt to answer the question, it probably shouldn't be posted here as an answer. Perhaps it would be better as a separate self-answered question. Jun 10, 2016 at 14:26
• Also for elegance consider ReIm@Exp[#.{1,I}]& and ReIm@Log[#.{1,I}]& Jun 10, 2016 at 14:29
• Apply[{\[Phi], r} \[Function] Through[{Re, Im}[Exp[r + I \[Phi]]]]] Apply[{x, y} \[Function] Through[{Re, Im}[Log[x + I y]]]] Hi @SimonWoods, live perhaps is a hunt for elegance, I got your elegant idea and gave it an elegant appearance. Now the Transform and its inverse are almost perfect symmetric. For reasons of elegance I try to avoid Slot[] in favour of the arrow by \ [Function]. By the way no ReIm[] on MMA 10. Jun 10, 2016 at 23:04