# Domain Coloring

My next project I'd like to start working on is Domain Coloring. I am aware of the beautiful discussion at:

How can I generate this "domain coloring" plot?

And I am studying it. However, a lot of the articles on domain coloring refer back to Hans Lundmark's page at:

http://www.mai.liu.se/~halun/complex/domain_coloring-unicode.html

So, I would like to begin my work by using Mathematica to draw these three images based on Hans' notes. I'd appreciate if anyone can provide some code that will produce these images, as I could use it to start my study of the rest of Hans' page.

A very small adjustment. Still learning.

g[{f_, cf_}] :=
DensityPlot[f, {x, -1, 1}, {y, -1, 1}, PlotPoints -> 51,
ColorFunction -> cf, Frame -> False];

g /@ {{Arg[-(x + I y)], "SolarColors"},
{Mod[Log[2, Abs[x + I y]], 1], GrayLevel}}

ImageMultiply @@ %


Not sure where to put my current question, so I'll update here. Just came back to visit and discovered some wonderful answers at the bottom of this list. I do understand the opening code:

f[z_] := (z + 2)^2*(z - 1 - 2 I)*(z + I)

paint[z_] :=
Module[{x = Re[z], y = Im[z]},
color = Blend[{Black, Red, Orange, Yellow},
Rescale[ArcTan[-x, -y], {-Pi, Pi}]];
shade = Mod[Log[2, Abs[x + I y]], 1];


But then I encounter difficulty with the following code:

ParametricPlot[{x, y}, {x, -3, 3}, {y, -3, 3},
ColorFunctionScaling -> False,
ColorFunction -> Function[{x, y}, paint[f[x + y I]]], Frame -> False,
Axes -> False, MaxRecursion -> 1, PlotPoints -> 50, Mesh -> 400,
PlotRangePadding -> 0, MeshStyle -> None, ImageSize -> 300]


I'm good with the first few lines. Looks like ParametricPlot is plotting points, where x and y both range from -3 to 3 (correct me if I am wrong). I also understand the ColorFunctionScaling and the ColorFunction lines. I understand Axes, PlotRangePadding, MeshStyle, and ImageSize. Where I am having trouble is with what PlotPoints->50 and Mesh->400 are doing.

First of all, my image size is 300. What does PlotPoints->50 mean? Does that mean it will sample and array of 50x50 points out of 300x300 and scale the results to fit in the domain [-3,3]x[-3,3]? My next question is, then those points get colored? And if so, how are the remainder of the points in the image colored? For example, I tried:

Table[ParametricPlot[{x, y}, {x, -3, 3}, {y, -3, 3},
ColorFunctionScaling -> False,
ColorFunction -> Function[{x, y}, paint[f[x + y I]]],
PlotPoints -> n, MeshStyle -> None], {n, 10, 50, 10}]


And the images got a little sharper as the PointPlots->n increased.

Here's another question. What does Mesh->400 do in this situation. For example, I tried lowering the mesh number:

ParametricPlot[{x, y}, {x, -3, 3}, {y, -3, 3},
ColorFunctionScaling -> False,
ColorFunction -> Function[{x, y}, paint[f[x + y I]]], Frame -> False,
Axes -> False, MaxRecursion -> 1, PlotPoints -> 50, Mesh -> 100,
PlotRangePadding -> 0, MeshStyle -> None, ImageSize -> 300]


And was completely surprised that it had an effect on the image, particularly when MeshStyle->None. Here's the image I get:

Why does setting Mesh->100 decrease the sharpness of the image?

One final question I have regards adding the mesh lines. Simon suggested
For the mesh you could do something like Mesh->{Range[-5,5],Range[-5,5]}, MeshStyle->Opacity[0.5], MeshFunctions->{(Re@f[#1+I #2]&),(Im@f[#1+I #2]&)} and cormullion added them to produce a beautiful result, but I tried this:

ParametricPlot[{x, y}, {x, -3, 3}, {y, -3, 3},
ColorFunctionScaling -> False,
ColorFunction -> Function[{x, y}, paint[f[x + y I]]], Frame -> False,
Axes -> False, MaxRecursion -> 1, PlotPoints -> 50,
Mesh -> {Range[-5, 5], Range[-5, 5]}, PlotRangePadding -> 0,
MeshStyle -> Opacity[0.5],
MeshFunctions -> {(Re@f[#1 + I #2] &), (Im@f[#1 + I #2] &)},
ImageSize -> 300]


And got this resulting image.

So I am clearly missing something. Maybe someone could post the code that gives cormullion's last image?

OK, just purchased and installed Presentations package. Tried this:

With[{f = Function[z, (z + 2)^2 (z - 1 - 2 I) (z + I)],
zmin = -2 - 2 I, zmax = 2 + 2 I,
colorFunction = Function[arg, HotColor[Rescale[arg, {-Pi, Pi}]]],
imgSize = 400},
Draw2D[{ComplexDensityDraw[Arg[f[z]], {z, zmin, zmax},
ColorFunction -> colorFunction, ColorFunctionScaling -> False,
Mesh -> 50, MeshFunctions -> {Function[{x, y}, Abs[f[x + I y]]]},
PlotPoints -> {50, 50}]}, Frame -> True, FrameLabel -> {Re, Im},
PlotLabel -> Row[{"Arg coloring and Abs mesh of ", f[z]}],
RotateLabel -> False, BaseStyle -> 12, ImageSize -> imgSize]]


But got this colorless image.

Any thoughts on how to fix this?

-
David, if you want to get the exact image you may want to look at this Wikipedia page. – Mr.Wizard Jan 2 '13 at 15:18
[1/2] David, seeing as you're a new user, I thought I'd like to point out how the reception to the question changed over the past week — when you first asked this, it was just a blind request question and you got two comments (one from me) saying that this wasn't a coding site. Subsequently, several others downvoted it and voted to close. It was at -5 at one point and 4 close votes (it stopped only after Mr.Wizard requested that it not be closed). However, after you edited it on Jan 2 to show some effort (even though it was mostly belisarius' answer), ... – R. M. Jan 8 '13 at 22:07
[2/2] ... you got a couple more very nice answers and got several upvotes which made up for all the initial downvotes. Please don't take the initial response as hostility towards newcomers. Rather, view it as a learning experience for what's expected on this site. People here are all willing to help, but they don't want to just be doing free work for someone who makes no effort themselves. Almost always, you'll get better responses if you show some enthusiasm and effort in getting an answer to your own question. Hope you find this site useful :) – R. M. Jan 8 '13 at 22:08
Re: your update, almost all of it is answered by the documentation for PlotPoints and Mesh. As for the code that generates cormullion's last plot, it would be better to ping them directly with a comment under their answer. – R. M. Jan 15 '13 at 22:02
At the risk of redundancy but in case somebody seeing this doesn't read beyond the first answer... As I said in a response to your comment upon my own Answer, below, you have to define HotColor as it was done in the first Answer, by @belisarius. – murray Jan 17 '13 at 22:00

Here is something quickly made, and similar to what you are after. You'll have to work out the details, though:

HotColor[ z_ ] :=
Which[
0 <= z <= 3/8,       RGBColor[z 8/3, 0, 0],
3/8 <= z <= 6/8,     RGBColor[1, (z - 3/8) 8/3, 0],
True,                RGBColor[1, 1, (z - 6/8) 8/6]
];
g[{f_, cf_}] := DensityPlot[f, {x, -1, 1}, {y, -1, 1}, PlotPoints -> 50,
ColorFunction -> cf, Frame -> False];
g /@ {{ArcTan[-x, -y], HotColor}, {SawtoothWave[3 Norm[{x, y}]],  GrayLevel}}
ImageMultiply @@ %


-
I hope you don't mind the edit. (+1 btw) – Mr.Wizard Jan 1 '13 at 14:28
@Mr.Wizard Not at all. It's really a pity to have to restrain oneself from giving a full answer and get the very image in the question. But I think doing that will spoil a student :( – Dr. belisarius Jan 1 '13 at 17:34
Interesting, as I don't know the blending used to get the very image. Do you? (Perhaps it is on the linked page but didn't look extensively.) EDIT: never mind, I got it. – Mr.Wizard Jan 2 '13 at 15:10
@Mr.Wizard As a last resort you may get the RGB components along a circle in the original image and make an ad-hoc gradient :) – Dr. belisarius Jan 2 '13 at 16:13
I meant the blending method, e.g. ImageMultiply, not the gradient function. It's "Overlay" in the linked Wikipedia article, FWIW. – Mr.Wizard Jan 2 '13 at 17:14

I was hoping that this question would get some good answers, but it must have been asked at a time when everyone was feeling a bit curmudgeonly after Christmas... :)

My understanding of Hans Lundmark's linked page is that there's basically a painting method that colors each point of the final plot depending on the value of a function which is applied to xy coordinates. So you don't have to start off with these images, you just specify how the points are going to be painted. (Of course, you can start off with these images and transform them using ImageApply, but that's not what he had in mind, I think.)

Unfortunately I couldn't find many ColoringFunctions that accepted x and y coordinates - some of the plot functions accept only a single value for coloring. But apparently ParametricPlot can take a ColorFunction that processes x and y coordinates. So here's a function, a painting method, and a ParametricPlot that draws something similar to what the linked page is doing:

f[z_] := z ^3;

paint[z_] :=
Module[{x = Re[z], y = Im[z]},
color =
Blend[{Black, Red, Orange, Yellow},
Rescale[ArcTan[-x , -y], {-Pi, Pi}]];
shade = Mod[Log[2, Abs[x + I y]], 1];

ParametricPlot[
{x, y}, {x, -2, 2}, {y, -2, 2},
ColorFunctionScaling -> False,
ColorFunction -> Function[{x, y}, paint[f[x + y I]]],
Frame -> False, Axes -> False, MaxRecursion -> 1, PlotPoints -> 50,
Mesh -> 400, PlotRangePadding -> 0,
MeshStyle -> None, ImageSize -> 300
]


which doesn't look too different from the original (given some fiddling and fudging, perhaps):

and his fourth degree polynomial (whatever that might be):

f(z) = (z + 2 )2 (z − 1 − 2i) (z + i),

converts to this in Mathematica:

f[z_] := (z + 2)^2 (z - 1 - 2 I) (z + I);


and looks like this:

compared with his original:

Now, I know what you're all going to say — where are the mesh lines? Well, I've been looking at the mesh lines documentation for some minutes now, and it makes little sense to me yet. That's something which you can perhaps look into - an exercise for the questioner...? Or perhaps a more experienced answerer.

Update: Simon's mesh functions get us closer to the original:

My intuitive understanding of the PlotPoints and Mesh options is that they determine the resolution or quality of the plot. You can use Mathematica to help you explore different combinations of settings (this takes a minute or so on my machine):

f[z_] := (z + 2)^2 (z - 1 - 2 I) (z + I);
paint[z_] :=
Module[{x = Re[z], y = Im[z]},
colour =
Blend[{Black, Red, Orange, Yellow},
Rescale[ArcTan[-x , -y + 0.0001], {-Pi, Pi}]];
shade = Mod[Log[2, Abs[x + I y]], 1];
g = Grid[Table[
ParametricPlot[
{x, y}, {x, -3., 3.}, {y, -3., 3.},
ColorFunctionScaling -> False,
ColorFunction -> Function[{x, y}, paint[f[x + y I]]],
Frame -> False, Axes -> False, MaxRecursion -> 3,
PlotPoints -> pp,
Mesh -> mesh,
PlotLegends ->
Placed[{StringJoin["☝", "Plot Points: ", ToString[pp],
", Mesh: ", ToString[mesh]]}, Below],
MeshStyle -> Directive[Opacity[0.3], Yellow],
MeshFunctions -> {(Re@f[#1 + I #2] &), (Im@f[#1 + I #2] &)},
ImageSize -> 500
],
{pp, {25, 75, 150}},
{mesh, {25, 150, 400}}]]

-
For the mesh you could do something like Mesh->{Range[-5,5],Range[-5,5]}, MeshStyle->Opacity[0.5], MeshFunctions->{(Re@f[#1+I #2]&),(Im@f[#1+I #2]&)} – Simon Woods Jan 8 '13 at 14:47
@SimonWoods Cool - that works! – cormullion Jan 8 '13 at 15:13
I tried the MeshFunctions (see my example above), but I don't believe I did it correctly. Can anyone adjust my code to get this image? – David Jan 16 '13 at 17:47
@david you need to adjust the PlotPoints and Mesh values, and also set MeshStyle -> Directive[Opacity[0.3], Yellow]. Oh, and the range is -3 to 3. I'll add the code tomorrow, perhaps... – cormullion Jan 16 '13 at 19:20
Thanks, looking forward to the code. Can't tell you how helpful you've been. Awesome work. – David Jan 17 '13 at 7:01

There are many approaches to domain coloring. The approach in @cormullion's answer is to use Re and Im as the mesh functions. Another way is to color points just by $\text{Arg}(f(z))$ and then superimpose contours for $\text{Abs}(f(z))$.

While that may readily be done in "pure Mathematica", by pulling out real and imaginary parts of complex numbers, it seems a bit silly to do so given that complex numbers are already built in to Mathematica (and, in a sense, Mathematica prefers working with complex numbers to real numbers). Here, David Park's Presentations add-on (http://home.comcast.net/~djmpark/DrawGraphicsPage.html) may be applied. For example:

With[{f = Function[z, (z + 2)^2 (z - 1 - 2 I) (z + I)],
zmin = -2 - 2 I, zmax = 2 + 2 I,
colorFunction = Function[arg, HotColor[Rescale[arg, {-Pi, Pi}]]],
imgSize = 400},
Draw2D[{ComplexDensityDraw[Arg[f[z]], {z, zmin, zmax},
ColorFunction -> colorFunction, ColorFunctionScaling -> False,
Mesh -> 50, MeshFunctions -> {Function[{x, y}, Abs[f[x + I y]]]},
PlotPoints -> {50, 50}]
},
Frame -> True, FrameLabel -> {Re, Im},
PlotLabel -> Row[{"Arg coloring and Abs mesh of ", f[z]}], RotateLabel -> False,
BaseStyle -> 12, ImageSize -> imgSize
]
]


-
As you can see in my example up above, I didn't get these colors when I tried your code. Any thoughts? – David Jan 16 '13 at 7:04
@David: HotColor is not a function either built into Mathematica or provided by Park's Presentations add-on. It's as defined in the 1st Answer, by belisarius. I didn't repeat the definition in my answer since it was already posted by him. Sorry for any confusion. – murray Jan 16 '13 at 16:35
(OT: @David: Welcome to "Presentations". Aside from the docs that come with it, much linked by the Presentations palette, be sure to see the interesting examples David Park has posted at st-andrews.ac.uk/~pl10/c/djmpark .) – murray Jan 16 '13 at 16:37
My bad, I shouldn't have used uppercase – Dr. belisarius Jan 17 '13 at 11:02

I might as well put in my own take. I tried a bit to match the coloring used by Hans in his pics, and here's the result:

SetAttributes[colorize, Listable];
colorize[z_] :=
Darker[Blend[{{0, Black}, {1/2, Red}, {3/4, Orange}, {1, Yellow}},
Mod[Arg[z]/(2 π), 1]], If[z == 0, 0, 1 - Mod[Log2[Abs[z]], 1]]]


Now, we can do any of two things with this function. The first possibility I shall present is similar to the approach Mark used in this answer for generating fractals; that is, directly generate an Image[] object using the coloring function:

Image[DeveloperToPackedArray[
colorize[N[Table[x + I y, {y, 10, -10, -1/20}, {x, -10, 10, 1/20}]]]
/. RGBColor -> List]]


This, of course, is non-adaptive, and does not allow for the superimposition of meshlines. We thus present the second bit, which is more or less equivalent to the approach proposed in other answers: the use of ParametricPlot[]:

ParametricPlot[{x, y}, {x, -15, 15}, {y, -15, 15}, Axes -> None,
ColorFunction -> Function[{x, y, u, v}, colorize[u + I v]],
ColorFunctionScaling -> False,
MeshStyle -> Opacity[2/5, GrayLevel[2/5]], PlotPoints -> 350]


Probably the only annoying thing about this procedure is that one needs a relatively large amount of plot points just for the colors to come out right. The good thing about this compared to the approach of using DensityPlot[] is that the color function is allowed to depend on the parameters themselves, in contrast to DensityPlot[] where the color function is dependent only on the heights.

Here's Hans's polynomial example:

With[{f = Function[z, (z + 2)^2 (z - 1 - 2 I) (z + I)]},
ParametricPlot[{x, y}, {x, -4, 4}, {y, -4, 4}, Axes -> None,
ColorFunction -> Function[{x, y, u, v}, colorize[f[u + I v]]],
ColorFunctionScaling -> False,
MeshFunctions -> {Function[{x, y, u, v}, Re[f[u + I v]]],
Function[{x, y, u, v}, Im[f[u + I v]]]},
MeshStyle -> Opacity[2/5, GrayLevel[2/5]], PlotPoints -> 400]]


Just for fun, here's a domain-colored map of Hans's polynomial example to the Riemann sphere, using a slightly modified coloring scheme:

SetAttributes[colorize, Listable];
colorize[z_] := Blend[{
Blend[{{0, Black}, {1/2, Red}, {3/4, Orange}, {1, Yellow}}, Mod[Arg[z]/(2 π), 1]],
ColorData[{"GrayTones", "Reverse"}][If[z == 0, 1, Mod[Log2[Abs[z]], 1]]]}]

With[{f = Function[z, (z + 2)^2 (z - 1 - 2 I) (z + I)], h = π/360},
tex = DeveloperToPackedArray[colorize[
f[Table[Cot[φ/2] Exp[I θ], {φ, h, π - h, h}, {θ, -π, π, h}] // N]] /.
RGBColor -> List];]

ParametricPlot3D[{Cos[θ] Sin[φ], Sin[θ] Sin[φ], Cos[φ]}, {θ, -π, π}, {φ, 0, π},
Axes -> False, Boxed -> False, Lighting -> "Neutral",
Mesh -> None, PlotPoints -> 95, PlotStyle -> Texture[tex]]


-