Filling Space with Pursuit Polygons

Question

I want to make a program which can fill a 2D space with "pursuit polygons".
The following picture will help you understand better what I mean.
You can also look up "pursuit curves" or mice problem or watch this gif

What I've tried so far:

First I tried to produce these polygons by rotation of polygons
the square for example

data = {{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}};
Graphics[{Table[{Scale[Rotate[Line[data], 3*x Degree], {x, x}]}, {x, 
0, 20}]}]

then I decided to use spirals
I realised that instead of using n spirals for every n-polygon,
I can use only one(!) and produce the effect I want with "fine tuning"
here are the "pursuit polygons" that I made

Triangle:

a = 41.9;
b = 100;
W = Table[{t^2*Cos[a*t], t^2*Sin[a*t]}, {t, 0, b}];
AppendTo[W, W[[-4]]];
ListPlot[W]
Graphics[Line[W]]

Square:

S = Table[{t^2*Cos[42.4*t], t^2*Sin[42.4*t]}, {t, 0, 160}];
AppendTo[S, S[[-4]]];
Graphics[Line[S]]

Pentagon

P = Table[{t^2*Cos[45.251*t], t^2*Sin[45.251*t]}, {t, 0, 220}];
AppendTo[P, P[[-5]]];
Graphics[Line[P]]

Then I tried to put all these together by rotating and scaling...
(but I'm not satisfied with the result)
notice that for every polygon I use also the "anti-clockwise" version of the polygon, which produces interesting results

T = Table[{t^2*Cos[41.9*t], t^2*Sin[41.9*t]}, {t, 0, 100}];
AppendTo[T, T[[-4]]]; (*Triangle*)
S = Table[{t^2*Cos[42.4*t], t^2*Sin[42.4*t]}, {t, 0, 160}];
AppendTo[S, S[[-4]]];   (*square*)
S1 = Table[{t^2*Cos[-42.4*t], t^2*Sin[-42.4*t]}, {t, 0, 160}];
AppendTo[S1, S1[[-4]]];    (*Anti-clockwise Square*)
P = Table[{t^2*Cos[45.251*t], t^2*Sin[45.251*t]}, {t, 0, 220}];
AppendTo[P, P[[-5]]];      (*Pentagon*)
P1 = Table[{t^2*Cos[-45.251*t], t^2*Sin[-45.251*t]}, {t, 0, 220}];
AppendTo[P1, P1[[-5]]];     (*Anti-clockwise Pentagon*)
Graphics[{Translate[
Rotate[Scale[Line[(T)], {2.09, 2.09}], -67 Degree], {29500, 3800}],
Rotate[Scale[Line[S + 0], {1, 1}], -30 Degree], 
Translate[
Rotate[Scale[Line[(S1)], {.98, .98}], -87.5 Degree], {41000, 
24700}], 
Translate[
Rotate[Scale[Line[P1], {.605, .665}], -20.5 Degree], {76500, 
6500}], Scale[
Translate[
Rotate[Scale[
  Line[(T)], {1.2, 1.7}], -108 Degree], {50900, -12500}], {1.91, 
1.31}], 
Translate[
Rotate[Scale[Line[(P)], {.525, .665}], -64 Degree], {3000, 36000}],
Translate[
Rotate[Scale[Line[(P)], {.61, .61}], -87 Degree], {59000, 58000}], 
Translate[
Rotate[Scale[Line[(T)], {2.1, 1.8}], 12 Degree], {82500, 39500}], 
Translate[
Rotate[Scale[Line[(T)], {2.3, 1.9}], -133 Degree], {32000, 
67000}]}]

Can you find a way to divide and fill any given space with pursuit polygons?
The result would look better if this could work with ANY convex polygon and not only the regular polygons that I used...

Answer 1

The reason that these are called "pursuit polygons" is because they are formed from a dynamical system in which different agents pursue each other. Example:

In this image, one agent starts in each corner of the triangle. The agent starting in the lower right corner pursues the agent starting in the top corner, the agent in the top corner pursues the agent in the lower left corner, and the agent in the lower left corner pursues the agent in the lower right corner.

Drawing lines in between the agents yields a series of triangles which shrink as the agents get closer to each other, and also rotate:

The corresponding differential equations are: $$\begin{align*} \dot{\mathbf{x}}_i &= \mathbf{x}_{i+1} - \mathbf{x}_i,\ i\in\{1,\dots,N-1\}\\ \dot{\mathbf{x}}_i &= \mathbf{x}_1 - \mathbf{x}_N,\ i=N, \end{align*}$$ where $N$ is the number of agents. The agents start in corners for the visualizations we are making, but in general they don't have to.

This code solves the differential equations using Euler integration and plots the lines between the agents:

integrate[pts_, h_] := NestList[Nest[step[#, h] &, #, 10] &, pts, 100]
step[pts_, h_] := pts - h (pts - RotateRight[pts])

wrapAround[pts_] := Append[pts, First[pts]]

PursuitPolygon[pts_, h_: 0.01] := Graphics[{
   Line[wrapAround[pts]],
   Line[wrapAround /@ integrate[pts, h]]
   }]

For regular polygons, this results in logarithmic spirals, as you have pointed out:

Grid@Partition[Table[PursuitPolygon[N@CirclePoints[n]], {n, 2, 10}], 2]

But this approach is more general and works for any convex polygons, as you required. The following is an example of how we can fill in a triangulated region with spirals using this system:

reg = DiscretizeRegion[Rectangle[], MaxCellMeasure -> 0.1]

Show[PursuitPolygon @@@ MeshPrimitives[reg, 2]]

Wikipedia's entry on the mice problem suggests that the mice move at unit speed, i.e. like this: $$\begin{align*} \dot{\mathbf{x}}_i &= \frac{\mathbf{x}_{i+1} - \mathbf{x}_i}{\|\mathbf{x}_{i+1} - \mathbf{x}_i\|},\ i\in\{1,\ldots,N-1\}\\ \dot{\mathbf{x}}_i &= \frac{\mathbf{x}_1 - \mathbf{x}_N}{\|\mathbf{x}_1 - \mathbf{x}_N\|},\ i=N, \end{align*}$$ This turns out to be difficult to implement because the denominator grows large as the points approach each other. I spoke to halirutan and MichaelE2 about this. Here is MichaelE2's solution, which I implemented for triangles but it could be implemented for any polygon.

PursuitPolygon[{pt1_, pt2_, pt3_}] := Module[{},

  {sol1x, sol1y, sol2x, sol2y, sol3x, sol3y} = NDSolveValue[{
     {x1'[t], y1'[t]} == Normalize[{x2[t], y2[t]} - {x1[t], y1[t]}],
     {x2'[t], y2'[t]} == Normalize[{x3[t], y3[t]} - {x2[t], y2[t]}],
     {x3'[t], y3'[t]} == Normalize[{x1[t], y1[t]} - {x3[t], y3[t]}],
     {x1[0], y1[0]} == pt1,
     {x2[0], y2[0]} == pt2,
     {x3[0], y3[0]} == pt3,
     WhenEvent[Norm[{x2[t], y2[t]} - {x1[t], y1[t]}] < 1*^-8, "StopIntegration"]},
    {x1, y1, x2, y2, x3, y3},
    {t, 0, 10}
    ];

  {tmin, tmax} = MinMax@sol1x["Grid"];

  Graphics@Table[Line[{{sol1x[t], sol1y[t]}, {sol2x[t], sol2y[t]}, {sol3x[t], sol3y[t]}, {sol1x[t], sol1y[t]}}], {t, tmin, tmax, (tmax - tmin)/15}]
  ]

Show[PursuitPolygon @@@ MeshPrimitives[reg, 2]]

The solution is qualitatively different when the agents move at unit speed because they no longer meet at the centroid. Note that the difference in the picture is also that the density of the lines is different, this density can be adjusted in both solutions but changing some of the fixed numbers I put in. 15 in (tmax - tmin)/15 controls the number of lines in the last solution.

Answer 2

Here's an extension of @C.E.'s code to non-triangular polygons - polygons taken from a Voronoi mesh of roughly equidistributed (https://mathematica.stackexchange.com/a/141215/3056) 31 points over a square. Code is awful, I didn't really have time to think.

ClearAll@PursuitPolygon;
PursuitPolygon[pts_] := Module[{vars, sols},
  vars = Table[Unique[], Length@pts, 2]; 
  sols = NDSolveValue[{Sequence @@ (((#'[t] & /@ #1) == 
           Normalize[Subtract @@ Map[#[t] &, {#2, #1}, {2}]]) & @@@ 
        Partition[Join[vars, {First@vars}], 2, 1]), 
     Sequence @@ MapThread[(#[0] & /@ #1) == #2 &, {vars, pts}], 
     WhenEvent[Norm[x] < 1*^-8, "StopIntegration"] /. 
      x -> Subtract @@ Map[#[t] &, Take[vars, 2], {2}]}, 
    Flatten@vars, {t, 0, 10}];
  {tmin, tmax} = MinMax@First[sols]["Grid"];
  Graphics@
   Table[Line[
     Map[#[t] &, Join[Partition[sols, 2], {Take[sols, 2]}], {2}]], {t,
      tmin, tmax, (tmax - tmin)/15}]]

ClearAll@reg;
reg = 
 With[{reg = Rectangle[]}, 
  With[{points = 31, samples = 4000, iterations = 100}, 
    Nest[With[{randoms = Join[#, RandomPoint[reg, samples]]}, 
       RegionNearest[reg][
        Mean@randoms[[#]] & /@ 
         Values@PositionIndex@Nearest[#, randoms]]] &, 
     RandomPoint[reg, points], iterations]] // 
   VoronoiMesh[#, {{0, 1}, {0, 1}}] &];

Show[PursuitPolygon /@ MeshPrimitives[reg, 2][[All, 1]]]