Arbitrary shapes in the plane

Ivan Cnop
Vrije universiteit Brussel
icnop@vub.ac.be

Trigonometric approximations for starlike shapes

Initialisation

We have to color directions and use coordinates for a value in the plane

<< Graphics`ArgColors`

coords[z_] := {Re[z], Im[z]}

Define  Spirograph

Drawing an approximate square using the spirograph

From Points in the plane to polar coordinates and back

Choosing points

Start with an empty graph and select points by mouseclicks along a curve that does not self-intersect and is starlike around the origin. 
For the first example we use a square.

Plot[0, {x, -2, 2}, PlotRange -> {-2, 2}, 
      AspectRatio -> Automatic, GridLines -> Automatic]

pts = {{0.991495, 0.0548687}, {0.991495, 0.115418}, {0.991495, 
            0.175967}, {0.991495, 0.236517}, {0.991495, 0.297066}, {0.991495, 
            0.357615}, {0.991495, 0.418165}, {1.00663, 0.478714}, {1.00663, 
            0.569538}, {1.00663, 0.630087}, {1.00663, 0.720911}, {1.00663, 
            0.781461}, {1.00663, 0.84201}, {1.00663, 0.902559}, {1.00663, 
            0.963109}, {0.946083, 0.993383}, {0.885534, 0.993383}, {0.809847, 
            0.993383}, {0.719023, 0.993383}, {0.643337, 0.993383}, {0.582788, 
            0.993383}, {0.522238, 0.993383}, {0.446551, 0.993383}, {0.34059, 
            0.978246}, {0.280041, 0.963109}, {0.219491, 0.932834}, {0.158942, 
            0.932834}, {0.0832555, 0.932834}, {0.00756884, 
            0.932834}, {-0.0832553, 0.932834}, {-0.174079, 
            0.932834}, {-0.234629, 0.932834}, {-0.295178, 
            0.932834}, {-0.355727, 0.932834}, {-0.416277, 
            0.932834}, {-0.476826, 0.932834}, {-0.537375, 
            0.932834}, {-0.597925, 0.932834}, {-0.688748, 
            0.947971}, {-0.749298, 0.947971}, {-0.809847, 
            0.963109}, {-0.900671, 0.963109}, {-0.961221, 
            0.963109}, {-0.991495, 0.902559}, {-0.991495, 
            0.84201}, {-0.991495, 0.781461}, {-0.991495, 
            0.720911}, {-0.991495, 0.660362}, {-0.991495, 
            0.599813}, {-0.991495, 0.539263}, {-0.976358, 
            0.478714}, {-0.976358, 0.418165}, {-0.976358, 
            0.357615}, {-1.00663, 0.297066}, {-1.00663, 
            0.236517}, {-1.02177, 0.175967}, {-1.02177, 
            0.115418}, {-1.03691, 
            0.0548687}, {-1.03691, -0.00568056}, {-1.05204, -0.0662301}, {-1.05204, -0.126779}, {-1.05204, -0.202466}, {-1.05204, -0.278153}, {-1.05204, -0.353839}, {-1.05204, -0.429526}, {-1.05204, -0.505213}, {-1.05204, -0.580899}, {-1.05204, -0.641449}, {-1.05204, -0.701998}, {-1.03691, -0.777685}, {-1.00663, -0.838234}, {-0.976358, -0.898783}, {-0.900671, -0.913921}, {-0.840122, -0.929058}, {-0.764435, -0.929058}, {-0.703886, -0.929058}, {-0.613062, -0.929058}, {-0.537375, -0.929058}, {-0.476826, -0.929058}, {-0.416277, -0.929058}, {-0.355727, -0.929058}, {-0.264903, -0.944195}, {-0.204354, -0.959333}, {-0.128667, -0.959333}, {-0.0681179, -0.97447}, {-0.0075686, -0.97447}, {0.0529807, -0.989607}, {0.128667, -0.989607}, {0.189217, -0.989607}, {0.249766, -0.989607}, {0.310315, -0.989607}, {0.370865, -0.989607}, {0.431414, -1.00474}, {0.522238, -1.00474}, {0.613062, -1.01988}, {0.673611, -1.01988}, {0.749298, -1.01988}, {0.809847, -1.01988}, {0.870397, -1.01988}, {0.930946, -0.989607}, {0.991495, -0.944195}, {1.02177, -0.883646}, {1.02177, -0.807959}, {1.02177, -0.717135}, {1.02177, -0.656586}, {0.991495, -0.550625}, {0.976358, -0.490075}, {0.961221, -0.384114}, {0.961221, -0.323565}, {0.946083, -0.247878}, {0.946083, -0.157054}, {0.946083, -0.0965046}};

Length[pts]

Here they are again:

ListPlot[pts, PlotRange -> {{-2, 2}, {-2, 2}}, 
      AspectRatio -> Automatic]

Compute direction angles:

angles[{a_, b_}] := ArcTan[b/a]

Map[angles, pts]

REMARK: this way we only obtain directions in the right half plane. 
      We will improve the definition below to cover the full circle of directions. 
      Otherwise we would have to take negative distances for points in the left half plane. 
      This is not the usual defintion of polar coordinates.

Compute distances from the origin

d[{a_, b_}] := Sqrt[(a^2) + (b^2)]

A better definition of angles

angle[{a_, b_}] := 
      If[a ≥ 0, t /. First[NSolve[Tan[t] == b/a, {t}]], 
        If[b ≥ 0, π - t /. 
            First[NSolve[Tan[t] == -b/a, {t}]], -π + t /. 
            First[NSolve[Tan[t] == b/a, {t}]]]]

Map[d, pts];

angles = Flatten[Map[angle, pts]]

comp = 
        Map[d, pts] (E^(I angles));

Show[
      Graphics[Transpose[{Map[ArgColor, comp], 
            Map[Point, Map[coords, comp]]}], 
        PlotRange -> {{-2, 2}, {-2, 2}}, 
        AspectRatio -> Automatic]]

Coloring values according to directions (optional)

The values in the upper half plane can be colored. This is done to illustrate the text-editing capabilities of   Mathematica   .

Take[
        Table[StyleBox[ToBoxes[d[pts[[i]]]], 
            Background -> Map[ArgColor, comp][[i]]], {i, 
            Length[pts]}], Floor[Length[pts]/2]] // DisplayForm

Coloring (ArgColor command) has to be modified to match arguments in the lower half plane.

Polynomial approximation

To approximate a function we have to use both cos and sin functions: cosine for the even part, sine for the odd part. Later we will use the identities (known from taking coordinates for directions)

2 Cos[
          t]  = (E^(I t))  +  (E^(-I t))
2 I Sin[
          t]  = (E^(I t))  -  (E^(-I t))

\(to\ change\ these\ into\ a\ \((rational)\)\ expression\ of\ the\ directions\ \ \ \ \ \(\(e\^\(i\ t\)\)\(\ \ \ \)\(.\)\)\)

standard = Flatten[{Cos[Range[0, 6] x], Sin[Range[6] x]}]

f[x_] = 
      Fit[Transpose[{Flatten[Map[angle, pts]], Map[d, pts]}], standard, 
        x]

ListPlot[Transpose[{Flatten[Map[angle, pts]], Map[d, pts]}], 
      PlotRange -> {0, 2}]

Remark: modify PlotRange to see more detail

Show[{%, Plot[f[x], {x, -π, π}]}, 
      PlotRange -> {0.8, 1.5}]

Conversion to the plane

\(A\ usual\ fit\ will\ not\ work\ since\ {f . 
                cos \((θ)\), \ 
            f . sin \((θ)\)}\ is\ not\ a\ polynomial\ in\ \ \ z\  = e\^\(i \((θ)\)\) . 
Since\ the\ function\ under\ review\ is\ periodic\), it is never possible to approximate these by polynomials in the variable theta .  
We first make a trigonometric approximation and then turn this in a complex shifted polynomial . 
The effect of reversing polar coordinates is simply multiplication by  z  . 
Indeed, again writing   f  for  f (θ):

x  =  f . cos (θ)
y  =  f . sin (θ)
x + I y  = f .  (E^(I (θ))) =  f . z

Making complex this real distance function

f[x] /. {Cos[x] ->  .5 (z + 1/z), 
        Sin[ x] ->  .5/I (z - 1/z), 
        Cos[n_ x] ->  .5 ((z^n) + (z^(-n))), 
        Sin[n_ x] ->  .5/I ((z^n) - 
                (z^(-n)))}

Expand[z %]

CoefficientList[% (z^7), z] // TableForm

Apply Spirograph

spirograph[%% /. z -> (E^(I t)), 
      t]

Show[%, PlotRange -> {{-2, 2}, {-2, 2}}, 
      AspectRatio -> Automatic]

This matches nicely with the given points. This is clear if we move the graphic and align the graphics frames.
Increasing the periods of cos and sin gives a higher degree polynomial and even better approximations.

Arbitrary curves

Example

Plot[0, {x, -2, 2}, PlotRange -> {-2, 2}, 
      AspectRatio -> Automatic, GridLines -> Automatic]

pts = {{0.426816, 1.15392}, {0.366382, 1.13881}, {0.260621, 
            1.07838}, {0.185077, 1.04816}, {0.094425, 1.00283}, {0.00377274, 
            0.972616}, {-0.0566621, 0.942398}, {-0.147314, 
            0.897072}, {-0.222858, 0.866855}, {-0.31351, 
            0.836637}, {-0.419271, 0.776203}, {-0.555249, 
            0.700659}, {-0.676119, 0.610007}, {-0.812097, 
            0.534463}, {-0.932967, 0.45892}, {-1.02362, 
            0.398485}, {-1.14449, 0.33805}, {-1.20492, 
            0.292724}, {-1.28047, 0.247398}, {-1.3409, 
            0.202072}, {-1.40134, 0.156746}, {-1.49199, 
            0.0812025}, {-1.55242, -0.00944972}, {-1.56753, -0.0849932}, {-1.62797, -0.190754}, {-1.64308, -0.296515}, {-1.6884, -0.387167}, {-1.70351, -0.538254}, {-1.74884, -0.659124}, {-1.74884, -0.840428}, {-1.74884, -0.946189}, {-1.73373, -1.03684}, {-1.6884, -1.15771}, {-1.65818, -1.21815}, {-1.62797, -1.27858}, {-1.53731, -1.38434}, {-1.46177, -1.44478}, {-1.40134, -1.4901}, {-1.3409, -1.53543}, {-1.28047, -1.56565}, {-1.17471, -1.61097}, {-1.06894, -1.6563}, {-0.993401, -1.6563}, {-0.88764, -1.70162}, {-0.766771, -1.71673}, {-0.676119, -1.71673}, {-0.600575, -1.6563}, {-0.54014, -1.62608}, {-0.449488, -1.59586}, {-0.358836, -1.52032}, {-0.283293, -1.45988}, {-0.19264, -1.36923}, {-0.147314, -1.3088}, {-0.0868795, -1.24836}, {0.00377274, -1.15771}, {0.0642076, -1.08217}, {0.124642, -1.00662}, {0.200186, -0.946189}, {0.275729, -0.900863}, {0.351273, -0.855537}, {0.411708, -0.810211}, {0.472142, -0.795102}, {0.562795, -0.764885}, {0.62323, -0.734667}, {0.683664, -0.719559}, {0.759208, -0.70445}, {0.819643, -0.689341}, {0.880077, -0.674232}, {0.955621, -0.674232}, {1.01606, -0.644015}, {1.15203, -0.58358}, {1.22758, -0.553363}, {1.28801, -0.508037}, {1.34845, -0.477819}, {1.40888, -0.432493}, {1.46932, -0.387167}, {1.49953, -0.326732}, {1.54486, -0.266298}, {1.55997, -0.205863}, {1.55997, -0.145428}, {1.55997, -0.0698845}, {1.55997, -0.00944972}, {1.54486, 
            0.0509851}, {1.52975, 0.126529}, {1.51464, 0.186963}, {1.45421, 
            0.247398}, {1.39377, 0.292724}, {1.31823, 0.307833}, {1.25779, 
            0.322942}, {1.19736, 0.33805}, {1.13693, 0.353159}, {1.07649, 
            0.383376}, {1.01606, 0.443811}, {0.985838, 0.504246}, {0.955621, 
            0.564681}, {0.910295, 0.625116}, {0.880077, 0.68555}, {0.84986, 
            0.745985}, {0.834751, 0.80642}, {0.819643, 0.866855}, {0.774317, 
            0.92729}, {0.72899, 0.987725}, {0.668555, 1.03305}, {0.62323, 
            1.09349}, {0.562795, 1.10859}};

Length[pts]

angles[{a_, b_}] := ArcTan[b/a]

Map[angles, pts];

ListPlot[pts, PlotRange -> {{-2, 2}, {-2, 2}}, 
      AspectRatio -> Automatic]

d[{a_, b_}] := Sqrt[(a^2) + (b^2)]

angle[{a_, b_}] := 
      If[a ≥ 0, t /. First[NSolve[Tan[t] == b/a, {t}]], 
        If[b ≥ 0, π - t /. 
            First[NSolve[Tan[t] == -b/a, {t}]], -π + t /. 
            First[NSolve[Tan[t] == b/a, {t}]]]]

Map[d, pts];

angles = Flatten[Map[angle, pts]];

comp = 
        Map[d, pts] (E^(I angles));

Show[
      Graphics[Transpose[{Map[ArgColor, comp], 
            Map[Point, Map[coords, comp]]}], 
        PlotRange -> {{-2, 2}, {-2, 2}}, 
        AspectRatio -> Automatic]]

Coloring values according to directions (optional)

Take[
        Table[StyleBox[ToBoxes[d[pts[[i]]]], 
            Background -> Map[ArgColor, comp][[i]]], {i, 
            Length[pts]}], Floor[Length[pts]/2]] // DisplayForm

Polynomial approximation

standard = Flatten[{Cos[Range[0, 6] x], Sin[Range[6] x]}]

f[x_] = 
      Fit[Transpose[{Flatten[Map[angle, pts]], Map[d, pts]}], standard, 
        x]

ListPlot[Transpose[{Flatten[Map[angle, pts]], Map[d, pts]}], 
      PlotRange -> {0, 2}]

Remark: modify PlotRange to see more detail

Show[{%, Plot[f[x], {x, -π, π}]}, 
      PlotRange -> {0.8, 1.5}]

Conversion to the plane

Make it complex

f[x] /. {Cos[x] ->  .5 (z + 1/z), 
          Sin[ x] ->  .5/I (z - 1/z), 
          Cos[n_ x] ->  .5 ((z^n) + (z^(-n))), 
          Sin[n_ x] ->  .5/I ((z^n) - 
                  (z^(-n)))};

Expand[z %];

CoefficientList[% (z^7), z] // TableForm

Apply Spirograph

spirograph[%% /. z -> (E^(I t)), 
      t]

Show[%, PlotRange -> {{-2, 2}, {-2, 2}}, 
      AspectRatio -> Automatic]

Animate this result over a copy of the data:

Correcting the mean

Shifting the picture

For statistical purposes it is often convenient to have a drawing that is balanced with respect to the origin.
One can hope approximation is better if the initial data do not protrude more in some directions. Therefore we shift the origin to the mean of the data.

Length[pts]

mean[pts_] := Plus @@ pts/Length[pts]

pts - Table[mean[pts], {Length[pts]}];

Take these as new data points

pts = %;

d[{a_, b_}] := Sqrt[(a^2) + (b^2)]

angle[{a_, b_}] := 
      If[a ≥ 0, t /. First[NSolve[Tan[t] == b/a, {t}]], 
        If[b ≥ 0, π - t /. 
            First[NSolve[Tan[t] == -b/a, {t}]], -π + t /. 
            First[NSolve[Tan[t] == b/a, {t}]]]]

Map[d, pts];

angles = Flatten[Map[angle, pts]];

comp = 
        Map[d, pts] (E^(I angles));

Show[
      Graphics[Transpose[{Map[ArgColor, comp], 
            Map[Point, Map[coords, comp]]}], 
        PlotRange -> {{-2, 2}, {-2, 2}}, 
        AspectRatio -> Automatic]]

Polynomial approximation

standard = Flatten[{Cos[Range[0, 6] x], Sin[Range[6] x]}]

f[x_] = 
      Fit[Transpose[{Flatten[Map[angle, pts]], Map[d, pts]}], standard, 
        x]

ListPlot[Transpose[{Flatten[Map[angle, pts]], Map[d, pts]}], 
      PlotRange -> {0, 2}]

Remark: modify PlotRange to see more detail

Show[{%, Plot[f[x], {x, -π, π}]}, 
      PlotRange -> {0.8, 1.5}]

Conversion to the plane

Making complex this real distance function

f[x] /. {Cos[x] ->  .5 (z + 1/z), 
        Sin[ x] ->  .5/I (z - 1/z), 
        Cos[n_ x] ->  .5 ((z^n) + (z^(-n))), 
        Sin[n_ x] ->  .5/I ((z^n) - 
                (z^(-n)))}

Expand[z %]

CoefficientList[% (z^7), z] // TableForm

Apply Spirograph

spirograph[%% /. z -> (E^(I t)), 
      t]

Show[%, PlotRange -> {{-2, 2}, {-2, 2}}, 
      AspectRatio -> Automatic]

Animate this result over a copy of the initial data:

Conclusion

This method can be used for compression of periodic images and smoothing of values of periodic functions. 
It is possible to define shapes by indicating a few points on the shape, or by the coefficients in the polynomial.

The method provides an alternative to the discrete Fourier transform in the case the initial data are not evenly spaced around the unit circle.
This may be useful in signal analysis, in generating sounds or in image compression.

Further applications are the use of conformal maps near the circle image to analyse stress and reinforcements in concrete using these shapes.

Statistical applications in two-dimensional data sets are still under investigation.