Symmetrical image transformation for a kaleidoscope-type image

Question

I've been trying to write a little kaleidoscope-type manipulate, using images rather than graphics, and using ImageTransformation as a way in. However, I'm starting to suspect that my approach isn't going to work. Before I abandon it, though, I thought I'd ask here if it can be rescued:

im = ImageResize[ExampleData[{"TestImage", "Mandrill"}], 100]

Manipulate[
 ImageTransformation[
     ImagePad[im, {{x, -x}, {-x, x}}, Padding -> "Reversed"],
     Function[{pt}, {ArcTan[-#2, #1] & @@ pt, Norm[pt]}], 
   DataRange -> Automatic, PlotRange -> {{-2, 2}, {-2, 2}}, 
   Padding -> "Reversed"], {x, 0, 100, 5}]

While the results are OK, it isn't symmetrical, and I'm not sure that the transform can be corrected to make it so?

(Strictly speaking this is a teleidoscope, rather than a kaleidoscope...)

Answer 1

I think the main problem is that ArcTan runs from $-\pi$ to $\pi$, but the DataRange in the x direction runs from 0 to 1, so you're actually using 6 and a bit copies of the original image to cover the whole circle. This is easily fixed by specifying explicit values for DataRange, i.e.

Manipulate[
 With[{sectors = 8},
  ImageTransformation[
   ImagePad[im, {{-x, x}, {-x, x}}, Padding -> "Reversed"], 
   Function[{pt}, {ArcTan[-#2, #1] & @@ (pt), Norm[pt]}], 250,
   DataRange -> {Pi/sectors {-1, 1}, {0, 1}}, 
   PlotRange -> {{-2, 2}, {-2, 2}}, Padding -> "Reversed"]], {x, 0, 
  100, 5}]

Here, sectors is the number of sectors in the transformed image.

Edit

Note that while ImageTransformation works it isn't very fast. If you are using Mathematica 8 you could use for example ParametricPlot in combination with TextureCoordinateFunction to get a similar but faster result.

First, we need a source image for the Texture. To get the right tiling where two neighbouring images are each other's reflection, I'm using the following:

imReflected = ImageAssemble[{
   {ImageReflect[im, Bottom], ImageRotate[im, Pi]},
   {im, ImageReflect[im, Left]}}]

Next, we need a TextureCoordinateFunction. To get the reflections of the shifted image in the circles r=a and the radial lines t=b where a and b are integers I'm using a triangle wave with period 2 and amplitude 1/2 which is vertically shifted by some offset corresponding to x in the solution above.

func[p_, x_] := x + TriangleWave[{0., 1/2.}, p/2 - 1/4]

The transformed image can then be plotted according to

Manipulate[
 ParametricPlot[{r Cos[2 Pi/sectors t], r Sin[2 Pi/sectors t]}, 
    {r, 0, 3}, {t, 0, sectors},
  Mesh -> None, BoundaryStyle -> None,
  Axes -> False,
  PlotStyle -> {Opacity[1], Texture[imReflected]},
  TextureCoordinateFunction -> ({func[#4, x], func[#3, x]} &),
  TextureCoordinateScaling -> False,
  PerformanceGoal -> "Quality", 
  PlotPoints -> {4, 1 + Round[40, sectors]},
  PlotRange -> {{-2, 2}, {-2, 2}}],
 {x, 0, 1},
 {{sectors, 6}, 4, 10, 2}]

Note that for the best result the number of plot points in the r and t direction should be one plus a multiple of the number of rings and sectors, respectively.

Edit

Inspired by faleichik's answer below, I've decided to implement his solution using textured triangles. There are basically two different triangles in the tessellation, one being the mirror image of the other. The other triangles are rotated and/or translated copies of either of these triangles which can be constructed from the base triangles using Rotate and Translate. Therefore, one way to get the tessellation is as follows:

Manipulate[
 DynamicModule[{texcrds, base, hex, im},
  im = ImageResize[ExampleData[{"TestImage", "Mandrill"}], 100];

  texcrds = offset + {{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}};
  base = {Polygon[{{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}}, 
     VertexTextureCoordinates -> texcrds],
    Polygon[{{0, 0}, {1, 0}, {1/2, -Sqrt[3]/2}}, 
     VertexTextureCoordinates -> texcrds]};
  hex = Rotate[base, #, {0, 0}] & /@ {0, 2 Pi/3, -2 Pi/3};

  Graphics[{Texture[im],
    Translate[Translate[hex, {{0, 0}, {3/2, -Sqrt[3]/2}}], 
     Tuples[{3 Range[0, n/2], Range[0, n] Sqrt[3]}]]},
   PlotRange -> {{-1/2, 3 n/2 + 1/2}, {-Sqrt[3]/2, n Sqrt[3]}}]],
 {{n, 3, "resolution"}, 1, 5, 1},
 {{offset, 0}, 0, 2}]

Answer 2

The topic has intrigued me for a long time, and this question became a final boost for sorting it out. Let me start from the result itself:

To see what is going on we should first analyse the structure of real kaleidoscope tesselations. Take a kaleidoscope picture from the inet and make some marking:

So we have six different types of equilaterial triangles which are the reflections of the original (centered) triangle. Things become easier when we introduce the basis $v=\{v_1,v_2\}$ (red vectors). Note that the marked coordinates on the picture are given in this basis. The periodic structure of tesselation in this basis is simple: its tile is the $[0,3)\times[0,3)$ square (diamond in original orthogonal coordinates).

Each tile consists of 9 diamonds combibed from two triangles. There are three types of diamonds (I, II and III), the corresponding triangle components are denoted t.1 and t.2, where t is the diamond type.

Now we need some basic yet boring linear algebra to start it all up. At the end we get the function mappt which for a given point on the plane determines the type of the triangle it belongs to and maps it to the original triangle to obtain the required value.

v = {v1, v2} = scale {{1, 0}, {1/2, Sqrt[3]/2}} // N;
vI = Inverse@Transpose@v;
patt = a : {a1_, a2_};
a2x[patt] := a1 v1 + a2 v2;
x2a[patt] := vI.a;

tI[patt] := a;
tI[patt /; 1 <= a1 + a2 < 2] := {1 - a2, 1 - a1};
tII[patt] := {1 - a2 - a1, a1};
tII[patt /; 1 <= a1 + a2 < 2] := {a1 + a2 - 1, 1 - a2};
tIII[patt] := {a2, 1 - a2 - a1};
tIII[patt /; 1 <= a1 + a2 < 2] := {1 - a1, a1 + a2 - 1};

tTable = {{tI, tII, tIII}, {tIII, tI, tII}, {tII, tIII, tI}} // 
   Transpose;
gett[patt] := tTable[[Sequence @@ ({1, 1} + Floor@Mod[a, 3])]];
map[patt] := gett[a][(# - Floor@#) &@Mod[a, 3]];

mappt = a2x@map@x2a@# &;

To demonstrate the result it is better to get non-symmetric original image so we have to cut the mandrill.

im = ImageResize[ImageTake[
   ExampleData[{"TestImage", "Mandrill"}], {150, 450}, {200, 500}], 50]

The object to manipulate is the scale which determines the width of the original equilaterial triangle.

Manipulate[
  scale = s; 
  ImageTransformation[ im, mappt, 200,
     DataRange -> Automatic,
     AspectRatio -> Automatic,
     PlotRange -> Sqrt[3]/(2 s) {{-1, 2}, {-1, 2}}
  ], 
  {s, .5, 1, .05},
  SynchronousUpdating -> False
]

That's it. The only issue now is that this implementation is quite slow. One might wish to make a compiled version of mappt, but I like this particular code pretty much.