Domain Coloring

Question

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];
  Darker[color, shade/4]]

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?

Answer 1

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 @@ %

Answer 2

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];
  Darker[color, shade/4]]

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];
  Darker[colour, shade/4]]
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, 
    PlotRangePadding -> 0,
    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}}]]

Answer 3

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
    ]
]

Answer 4

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[Developer`ToPackedArray[
      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 = Developer`ToPackedArray[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]]