# How to achieve faster performance on plotting complex valued functions

I'm very interested in visualizing complex valued functions $$f:\mathbb{C} \rightarrow \mathbb{C}$$. Especially interesting to me are Möbius transformations.

I found a nice way of visualizing these kind of functions at

How do I put an image on the complex plane

I do use the code provided by J. M. very often.

Here is an example:

c=1/2;
image=


Using the Möbius transformation $$\dfrac{z-c}{\overline{c}z-1}$$ which maps the unit circle onto the unit circle itself, we get:

Image[ImageForwardTransformation[image, Through[{Re, Im}[((#[[1]] + I #[[2]])- c)/(Conjugate[c]*(#[[1]] + I #[[2]]) - 1)]] &, Background -> 1, DataRange -> {{-1, 1}, {-1, 1}}, PlotRange -> {{-1, 1}, {-1, 1}}]]


As you can see, the image quality is a little bit on the lower side at some spots, which is totally expected using this plot. I do get better results using input images with a higher pixel count. I generally like to use 4096×4096 px, which takes some time to plot. I'm asking if there is any way of speeding up this plotting method using bigger input images?

The plot takes about 2-3 minutes. Also Mathematica does not make use of the 12 logical CPU cores my machine has. Is there maybe a way of splitting this up for multicore performance?

• You might be able to squeeze out some speed by explicitly determining the expressions for the real and imaginary parts, and then using those in Image(Forward)Transformation[]. – J. M.'s ennui Oct 21 '18 at 17:09
• I have the same problem and haven't found a good answer. The best case scenario is if you can shoehorn your problem to fit as a texture. ImageTransformation and ImageForwardTransformation are incredibly slow, but Wolfram believes they are fast enough and won't do anything about it. – Robert Jacobson Oct 21 '18 at 17:20
• My approach would be to construct a GraphicsComplex for the initial image and apply the transformation to the packed coordinate array of the complex. – Michael E2 Oct 21 '18 at 17:50
• See also my related answer here which gives an example of the inverse transformation trick which reduces computation and gives antialiased images. Essentially do an inverse Moebius transformation on the complex coordinates of your result pixels and these will be the (interpolated) coordinates for the texture lookup. – Thies Heidecke Oct 21 '18 at 21:46

image=Import["https://i.stack.imgur.com/GZcUT.jpg"]
c=1/2;
ImageForwardTransformation[image,
Through[{Re,Im}[((#[[1]]+I #[[2]])-c)/(Conjugate[c]*(#[[1]]+I #[[2]])-1)]]&,
Background->1,DataRange->{{-1,1},{-1,1}},PlotRange->{{-1,1},{-1,1}}]//AbsoluteTiming


Compile transformation function

expr=ReIm[((#[[1]]+I #[[2]])-c)/(Conjugate[c]*(#[[1]]+I #[[2]])-1)]&[{x,y}]//
ComplexExpand//Simplify

cf=With[{expr=expr},Compile[{{v,_Real,1}},Block[{x=v[[1]],y=v[[2]]},expr]]];

ImageForwardTransformation[image,cf,
Background->1,DataRange->{{-1,1},{-1,1}},PlotRange->{{-1,1},{-1,1}}]//AbsoluteTiming


Using ParametricPlot with Texture

ParametricPlot[expr,{x,-1,1},{y,-1,1},
PlotStyle->{Opacity[1],Texture[image]},ImageSize->ImageDimensions[image],
BoundaryStyle->None,Axes->False,Frame->False,PlotRange->1]//AbsoluteTiming


Here's a quick idea of what I had in mind in my comment:

plot = ParametricPlot[r {Cos[t], Sin[t]}, {r, 0, 1}, {t, 0, 2 Pi},
PlotPoints -> {35, 201}, Mesh -> 11, Axes -> False, Frame -> False,
ColorFunction -> (Hue[#4] &)];

{gc} = Cases[plot, _GraphicsComplex];

xf[c_] := ReIm@((# - c)/(Conjugate[c] # - 1) &)@(#.{1, I}) &;

Manipulate[
Graphics@MapAt[xf[c.{1, I}], gc, 1],
{{c, {0, 0}}, {-1, -1}, {1, 1},
TrackingFunction -> ((c = #/Max[1, 1.1 Norm[#]]) &)}
]


It could be made even quicker with a more efficient initial image plot.

Update: Here's a more efficient gc, using single polygons for each block:

ClearAll[polarRect];
polarRect[{r1_, r2_}, {t1_, t2_}, n_: 60] := Polygon@Join[
Table[
r2 {Cos[t], Sin[t]}, {t, Subdivide[t1, t2, 2 + Round[n*r2/(t2 - t1)]]}],
Rest@Table[r {Cos[t2], Sin[t2]}, {r, r2, r1, (r1 - r2)/Round[n/4]}],
{r1 {Cos[t2], Sin[t2]}},
If[r1 == 0, {},
Table[
r1 {Cos[t], Sin[t]}, {t, Subdivide[t2, t1, 2 + Round[n*r2/(t2 - t1)]]}]],
Rest@Table[r {Cos[t1], Sin[t1]}, {r, r1, r2, (r2 - r1)/Round[n/4]}]
];

With[{dr = 1/10, dt = Pi/4},
base = {EdgeForm[Directive[Thin, Black]],
Table[
{Hue[t/2/Pi],
N@polarRect[{r, r + dr}, {t, t + dt}]},
{r, 0, 1 - dr, dr}, {t, 0, 2 Pi - dt, dt}
]
}
];
coords = DeleteDuplicates@
DeveloperToPackedArray[
Flatten[Cases[base, Polygon[p_] :> p, Infinity], 1], Real];
nf = Nearest[coords -> "Index"];
gc = GraphicsComplex[coords, base /. Polygon -> Polygon@*Flatten@*nf];
`