weavingspace.symmetry

  1"""Classes for working with symmetries."""
  2
  3from __future__ import annotations
  4
  5import inspect
  6from dataclasses import dataclass
  7from typing import TYPE_CHECKING
  8
  9import geopandas as gpd
 10import numpy as np
 11import shapely
 12import shapely.affinity as affine
 13import shapely.geometry as geom
 14from matplotlib import pyplot as plt
 15
 16from weavingspace import tiling_utils
 17
 18if TYPE_CHECKING:
 19  from collections.abc import Sequence
 20
 21
 22@dataclass(slots=True)
 23class KMP_Matcher:  # noqa: N801
 24  """Class to find matching subsequences in a sequence."""
 25
 26  sequence:Sequence
 27  """Iterable in which subsequences are to be found."""
 28
 29  def find_matches(
 30      self,
 31      pattern:Sequence[tuple[float,...]]) -> list[int]:
 32    """Find matches in supplied pattern using Knuth-Morris-Pratt algorithm.
 33
 34    See:
 35    https://en.wikipedia.org/wiki/Knuth-Morris-Pratt_algorithm which provides
 36    detailed pseudo-code on which this code is directly based. See also:
 37
 38    Knuth DE, JH Morris Jr, and VR Pratt. 1977. Fast pattern  matching in
 39    strings. SIAM Journal on Computing 6(2): 323-350. doi: 10.1137/0206024.
 40
 41    This implementation expects sequences of tuples of floats, although in
 42    principle any objects could be contained in the Sequence.
 43
 44    Args:
 45      pattern (Sequence[tuple[float,...]]): The sequence to match.
 46
 47    Returns:
 48      Sequence[int]: the index positions at which matches are found.
 49
 50    """
 51    j, k, = 0, 0
 52    finds = []
 53    table = self._get_table(pattern)
 54    while j < len(self.sequence):
 55      if self._equal_tuples(pattern[k], self.sequence[j]):
 56        j = j + 1
 57        k = k + 1
 58        if k == len(pattern):
 59          finds.append(j - k)
 60          k = table[k]
 61      else:
 62        k = table[k]
 63        if k < 0:
 64          j = j + 1
 65          k = k + 1
 66    return finds
 67
 68  def _get_table(
 69      self,
 70      pattern:Sequence) -> Sequence[int]:
 71    """Return 'offsets' table as required by self.find_matches.
 72
 73    Note that the terse variable names are those used in the KMP 1977 paper.
 74
 75    Args:
 76      pattern (Iterable): The pattern to set up the table for.
 77      equals (function, optional): A function to use to test for equality of
 78        elements in patterns. Defaults to _equal_tuples().
 79
 80    Returns:
 81      Iterable[int]: _description_
 82
 83    """
 84    pos = 1
 85    cnd = 0
 86    T = {0: -1}
 87    while pos < len(pattern):
 88      if np.allclose(pattern[pos], pattern[cnd]):
 89        T[pos] = T[cnd]
 90      else:
 91        T[pos] = cnd
 92        while cnd >= 0 and not self._equal_tuples(pattern[pos], pattern[cnd]):
 93          cnd = T[cnd]
 94      pos = pos + 1
 95      cnd = cnd + 1
 96    T[pos] = cnd
 97    return tuple(T.values())
 98
 99  def _equal_tuples(
100      self,
101      t1:Sequence[float],
102      t2:Sequence[float]) -> bool:
103    """Test for near equality of two iterables of floats.
104
105    We use numpy.allclose to avoid issues with floating point exact equality
106    tests. Wrapped to allow find_matches function to take alternative equality
107    tests as input.
108
109    Args:
110      t1 (Sequence[float]): First set of floats.
111      t2 (Sequence[float]): Second set of floats.
112
113    Returns:
114      bool: True if t1 == t2 for all corresponding items in each sequence.
115
116    """
117    return np.allclose(t1, t2, atol = tiling_utils.RESOLUTION * 10)
118
119  def _round_tuple(
120      self,
121      t:Sequence[float],
122      digits:int = 3) -> Sequence[float]:
123    """Round all floats in a sequence.
124
125    Args:
126      t (Sequence[float]): the floats to round.
127      digits (int): number of digits to round to.
128
129    Returns:
130      Sequence[float]: the rounded Sequence.
131
132    """
133    return tuple([np.round(x, digits) for x in t])
134
135
136@dataclass(slots=True, repr=False)
137class Transform:
138  """Class to store details of a transform and draw it."""
139
140  transform_type:str
141  """Type of transform 'rotation', 'reflection', 'translation' or 'identity'."""
142  angle:float
143  """Angle of rotation (degrees)."""
144  centre:geom.Point
145  """Centre of the transformation."""
146  translation:tuple[float,...]
147  """X and Y coordinates shifts of a translation transform. A glide reflection
148  may also include this."""
149  transform:tuple[float,...]
150  """Six element tuple for the transform in shapely.transform format. See
151  https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
152  and methods in `weavingspace.tiling_utils`."""
153
154  def __str__(self) -> str:
155    if self.transform_type == "rotation":
156      return (
157        f"{self.transform_type:<11}{self.angle:>6.1f}° "
158        f"POINT ({self.centre.x:6.1f} {self.centre.y:6.1f}) "
159        f"{tuple([float(np.round(x, 3)) for x in self.transform])}"
160      )
161    if self.transform_type == "translation":
162      return (
163        f"{self.transform_type:<11}{'':>14}"
164        f"({self.translation[0]:>6.1f},{self.translation[0]:>6.1f}) "
165        f"{tuple([float(np.round(x, 3)) for x in self.transform])}"
166      )
167    if self.transform_type == "reflection":
168      return (
169        f"{self.transform_type:<11}{self.angle:>6.1f}° "
170        f"POINT ({self.centre.x:6.1f} {self.centre.y:6.1f}) "
171        f"{tuple([float(np.round(x, 3)) for x in self.transform])}"
172      )
173    return (
174      f"{self.transform_type:<11}{'':>30}"
175      f"{tuple([float(np.round(x, 3)) for x in self.transform])}"
176    )
177
178
179  def __repr__(self) -> str:
180    return str(self)
181
182
183  def apply(
184      self,
185      geometry:shapely.Geometry) -> shapely.Geometry:
186    """Apply this transform to supplied shapely geometry and return result.
187
188    Args:
189      geometry (shapely.Geometry): a shapely geometry to transform.
190
191    Returns:
192      shapely.Geometry: the resulting transformed shapely geometry.
193
194    """
195    return affine.affine_transform(geometry, self.transform)
196
197
198  def draw(
199      self,
200      ax:plt.Axes,
201      **kwargs:str,
202    ) -> plt.Axes:
203    """Draw this transform on the supplied axes.
204
205    Arguments specific to each transform type are supplied as **kwargs and are
206    documented in `draw_rotation`, `draw_reflection`, and `draw_translation`.
207
208    Args:
209      ax (pyplot.axes): the axes on which to draw the transform.
210      kwargs: passed on to specific draw method.
211
212    Returns:
213      pyplot.Axes: the axes with the rendering of this transform added.
214
215    """
216    match self.transform_type:
217      case "rotation":
218        args = inspect.signature(self.draw_rotation).parameters
219        rotn_args = {k: kwargs.pop(k) for k in dict(kwargs) if k in args}
220        return self.draw_rotation(ax = ax, **rotn_args)
221      case "reflection":
222        args = inspect.signature(self.draw_reflection).parameters
223        refn_args = {k: kwargs.pop(k) for k in dict(kwargs) if k in args}
224        return self.draw_reflection(ax = ax, **refn_args)
225      case "translation":
226        args = inspect.signature(self.draw_translation).parameters
227        trans_args = {k: kwargs.pop(k) for k in dict(kwargs) if k in args}
228        return self.draw_translation(ax = ax, **trans_args)
229      case "identity":
230        w = ax.get_window_extent().width
231        return ax.set_title(f"{self.transform_type}", fontsize = w / 20)
232
233
234  def draw_rotation(
235      self,
236      ax:plt.Axes,
237      radius:float = 200,
238      add_title:bool = True,
239    ) -> plt.Axes:
240    """Draw a rotation transform on the supplied Axes.
241
242    Args:
243      ax (pyplot.Axes): axes on which to draw the representation.
244      radius (float): radius in object units to draw arcs representing the
245        angle. Defaults to 200.
246      add_title (bool): whether or not to add a title to the drawing.
247        Defaults to True.
248
249    Returns:
250      the Axes with the rendering of this transform added.
251
252    """
253    x, y = self.centre.x, self.centre.y
254    axis = geom.LineString([(x, y), (x + radius * 1.25, y)])
255    arc = geom.LineString([
256      geom.Point(x + radius * np.cos(np.radians(a)),
257                 y + radius * np.sin(np.radians(a)))
258      for a in np.linspace(0, self.angle, 50)])
259    gpd.GeoSeries(self.centre).plot(
260      ax = ax, color = "k", markersize = 4, zorder = 5)
261    gpd.GeoSeries([axis, arc]).plot(ax = ax, color = "k", lw = .5)
262    if add_title:
263      w = ax.get_window_extent().width
264      ax.set_title(
265        f"{self.transform_type} {round(self.angle, 2)}", fontsize = w / 20)
266    return ax
267
268
269  def draw_reflection(
270      self,
271      ax:plt.Axes,
272      w:float = 5,
273      mirror_length:float = 100,
274      add_title:bool = True,
275    ) -> plt.Axes:
276    """Draw a reflection transform on the supplied Axes.
277
278    Args:
279      ax (pyplot.Axes): axes on which to draw the representation.
280      w (float): linewidth for the mirror line. Defaults to 5.
281      mirror_length (float): length of the mirror line in units of the
282        objects being reflected. Defaults to 100.
283      add_title (bool): whether or not to add a title to the drawing.
284        Defaults to True.
285
286    Returns:
287      the Axes with the rendering of this transform added.
288
289    """
290    x, y = self.centre.x, self.centre.y
291    dx, dy = self.translation
292    r = np.sqrt(self.translation[0] ** 2 + self.translation[1] ** 2)
293    mirror = geom.LineString([
294      (x - mirror_length / 2 * np.cos(np.radians(self.angle)),
295       y - mirror_length / 2 * np.sin(np.radians(self.angle))),
296      (x + mirror_length / 2 * np.cos(np.radians(self.angle)),
297       y + mirror_length / 2 * np.sin(np.radians(self.angle)))])
298    gpd.GeoSeries([mirror]).plot(
299      ax = ax, color = "k", lw = 1, ls = "dashdot", zorder = 5)
300    no_slide = np.isclose(r, 0, rtol = 1e-6, atol = 1e-6)
301    if not no_slide:
302      ax.arrow(
303        x - dx / 2, y - dy / 2, dx, dy, length_includes_head = True,
304        width = w, fc = "k", ec = None, head_width = w * 6, zorder = 5)
305    if add_title:
306      w = ax.get_window_extent().width
307      ax.set_title(
308        f"{self.transform_type} {round(self.angle, 2)}", fontsize = w / 20)
309    return ax
310
311
312  def draw_translation(
313      self,
314      ax:plt.Axes,
315      c:geom.Point,
316      w:float = 5,
317      add_title:bool = True,
318    ) -> plt.Axes:
319    """Draw a translation transform on the supplied Axes.
320
321    Args:
322      ax (pyplot.Axes): axes on which to draw the representation.
323      c (geom.Point): an origin point in the coordinate space of the objects
324        from which to start the arrow representing the translation.
325      w (float): linewidth for the arrow. Defaults to 5.
326      add_title (bool): whether or not to add a title to the drawing.
327        Defaults to True.
328
329    Returns:
330      the Axes with the rendering of this transform added.
331
332    """
333    gpd.GeoSeries([c]).plot(ax = ax, color = "k")
334    ax.arrow(c.x, c.y, self.translation[0], self.translation[1],
335             length_includes_head = True, lw = 0.5, width = w,
336             fc = "k", ec = None, head_width = w * 6, zorder = 5)
337    if add_title:
338      w = ax.get_window_extent().width
339      ax.set_title(
340        (f"{self.transform_type} "
341         f"({self.translation[0]:.1f}, {self.translation[1]:.1f})"),
342        fontsize = w / 20)
343    return ax
344
345
346@dataclass(slots=True, init=False)
347class Symmetries:
348  """Class to identify and store symmetries of a shapely.Polygon."""
349
350  polygon:geom.Polygon
351  """Polygon from which these symmetries are derived."""
352  matcher:KMP_Matcher
353  """the subsequence matcher used in symmetry detection."""
354  n:int
355  """number of vertices of the polygon."""
356  poly_code:list[tuple[float,...]]
357  """the encoding of the polygon as sequence of length, angle pairs."""
358  poly_code_r:list[tuple[float,...]]
359  """the reversed encoding used to detect reflection symmetries."""
360  rotation_shifts:list[int]
361  """list of number of 2pi/n rotation symmetries."""
362  reflection_shifts:list[int]
363  """list of pi/n relection angle symmetries."""
364  symmetries:list[Transform]
365  """list of Transform objects with more complete information."""
366  symmetry_group:str
367  """the code denoting the symmetry group"""
368
369  def __init__(self, polygon:geom.Polygon) -> None:
370    """Construct a Symmetries instance."""
371    self.polygon = polygon
372    self.n = shapely.count_coordinates(self.polygon) - 1
373    self.poly_code = self._get_polygon_code(self.polygon)
374    self.poly_code_r = self._get_polygon_code(self.polygon, mirrored = True)
375    self.matcher = KMP_Matcher(self.poly_code + self.poly_code[:-1])
376    self.symmetries = self.get_symmetries()
377    self.symmetry_group = self.get_symmetry_group_code()
378
379
380  def _get_polygon_code(
381      self,
382      polygon:geom.Polygon,
383      mirrored:bool = False,
384    ) -> list[tuple[float,...]]:
385    r"""Return list of length-angle pairs to encode a polygon shape.
386
387    This encoding is unique up to scale. The alignment of lengths and angles is
388    important and symmetry detection is sensitively dependent on it...
389
390    The forward (i.e. unmirrored) polygon pairs are of edge lengths and the
391    angle between each edge and its successor going CW around the polygon.
392
393           A[i] ----L[i]---- A (i+1)
394          /                   \
395         /                     \
396
397    i.e. edge i is between angles i and i + 1, and angle i is between edge
398    i-1 and i. The mirrored encoding pairs length[i] with angle[i+1].
399
400    The mirrored encoding pairs matched indexes of lengths and angles in the
401    original polygon. This means the angle is still the one between and edge
402    and its successor but proceeding CCW around the polygon.
403
404    See in particular for clarification:
405
406      Eades P. 1988. Symmetry Finding Algorithms. In Machine Intelligence and
407      Pattern Recognition, ed. GT Toussaint, 6:41-51. Computational Morphology.
408      North-Holland. doi: 10.1016/B978-0-444-70467-2.50009-6.
409
410    Args:
411      polygon (geom.Polygon): the polygon to encode
412      mirrored (bool, optional): if true encoding will be in CCW order.
413        Defaults to False.
414
415    Returns:
416      list[tuple[float]]: a list of side-length angle pairs.
417
418    """
419    # get edge lengths and angles
420    lengths = tiling_utils.get_side_lengths(polygon)
421    raw_angles = tiling_utils.get_interior_angles(polygon)
422    if mirrored:
423      lengths_r = lengths[::-1]
424      angles_r = raw_angles[::-1]
425      return list(zip(lengths_r, angles_r, strict = True))
426    angles = raw_angles[1:] + raw_angles[:1]
427    return list(zip(lengths, angles, strict = True))
428
429
430  def get_symmetries(self) -> list[Transform]:
431    """Find rotation and reflection symmetries of a polygon.
432
433    The work is delegated to `get_rotations()` and `get_reflections()`.
434
435    Based on
436
437      Eades P. 1988. Symmetry Finding Algorithms. In Machine Intelligence and
438      Pattern Recognition, ed. GT Toussaint, 6:41-51. Computational Morphology.
439      North-Holland. doi: 10.1016/B978-0-444-70467-2.50009-6.
440
441    and also
442
443      Wolter JD, TC Woo, and RA Volz. 1985. Optimal algorithms for symmetry
444      detection in two and three dimensions. The Visual Computer 1(1): 37-48.
445      doi: 10.1007/BF01901268.
446
447    Details in these papers are not unambiguous. This implementation was
448    developed based on them, but with some trial and error to get the index
449    offsets and (especially) retrieval of the reflection axes angles to work.
450
451    Returns:
452      list[Transform]: a list of Transform objects representing the polygon
453        symmetries.
454
455    """
456    return self.get_rotations() + self.get_reflections()
457
458
459  def get_rotations(
460      self,
461      offsets:None|list[int] = None,
462    ) -> list[Transform]:
463    """Get the rotations associated with this collection of symmetries.
464
465    Returns:
466      list[Transform]: A list of the rotation symmetries associated with this
467        Symmetries instance's polygon.
468
469    """
470    if offsets is None:
471      self.rotation_shifts = self.matcher.find_matches(self.poly_code)
472      offsets = self.rotation_shifts
473    c = self.polygon.centroid
474    rot_angles = [np.round(i * 360 / self.n, 6) for i in offsets]
475    rotation_transforms = [tuple(
476      np.round(tiling_utils.get_rotation_transform(a, (c.x, c.y)), 6))
477      for a in rot_angles]
478    return [
479      Transform("rotation", angle, c, (0, 0), transform) for
480      transform, angle in zip(rotation_transforms, rot_angles, strict = True)]
481
482
483  def get_reflections(
484      self,
485      offsets:None|list[int] = None,
486    ) -> list[Transform]:
487    """Get the reflections associated with this collection of symmetries.
488
489    Args:
490      offsets (list[int]): list of integer offsets i where i / n * 360 where n
491        is the number of sides of this Symmetries instance's polygons is the
492        angle of the mirror in which a reflection symmetry has been identified.
493
494    Returns:
495      list[Transform]: A list of the reflection symmetries associated with this
496        Symmetries instance's polygon.
497
498    """
499    # collate all possible reflection axes - ordering of these is important
500    # for correct picking - so don't mess with this code!
501    bearings = tiling_utils.get_side_bearings(self.polygon)
502    reflection_axes = []
503    interior_angles = tiling_utils.get_interior_angles(self.polygon)
504    if offsets is None:
505      self.reflection_shifts = self.matcher.find_matches(self.poly_code_r)
506      offsets = self.reflection_shifts
507    for i in range(self.n):
508      # the angle bisectors
509      reflection_axes.append(bearings[i] - interior_angles[i] / 2)
510      # the edge_perpendicular bisectors
511      reflection_axes.append(bearings[i] - 90)  # normal to the edge
512    # pick out those where matches occurred using the offsets
513    ref_angles = [np.round(reflection_axes[i], 6) for i in offsets]
514    c = self.polygon.centroid
515    reflection_transforms = [tuple(
516      np.round(tiling_utils.get_reflection_transform(angle, (c.x, c.y)), 6))
517      for angle in ref_angles]
518    return [
519      Transform("reflection", angle, c, (0, 0), transform) for
520      transform, angle in zip(reflection_transforms, ref_angles, strict = True)]
521
522  def get_symmetry_group_code(self) -> str:
523    """Return symmetry group code for this instance's polygon.
524
525    See https://en.wikipedia.org/wiki/Symmetry_group#Two_dimensions.
526
527    Returns:
528        str: code such as C2 or D4 for cyclic or dihedral symmetry groups.
529
530    """
531    if len(self.reflection_shifts) == 0:
532      # no reflection symmetries, so it is a cyclic group Cn
533      return f"C{len(self.rotation_shifts)}"
534    # must be the dihedral group Dn
535    return f"D{len(self.rotation_shifts)}"
536
537
538  def get_corner_offset(
539      self,
540      poly2:geom.Polygon) -> int|None:
541    """Find offset between corners of this polygon and another one.
542
543    The offset is determined under whatever transformation will match one on
544    to the other. Used by the `Topology` class. The offset is expressed as the
545    number of corners clockwise by which the supplied polygon is offset from
546    this one.
547
548    Args:
549      poly2 (geom.Polygon): polygon for which the offset is desired.
550
551    Returns:
552      int|None: integer number of corners clockwise by which the polygon
553        is offset from this one, or None if they don't match.
554
555    """
556    s2 = Symmetries(poly2)
557    matches = self.matcher.find_matches(s2.poly_code)
558    if len(matches) > 0:
559      return matches[0]
560    matches = self.matcher.find_matches(s2.poly_code_r)
561    if len(matches) > 0:
562      return -matches[0]
563    return None
564
565
566@dataclass(slots=True, init=False)
567class ShapeMatcher:
568  """Class to manage finding transforms that map one polygon on to another.
569
570  This was previously handled 'internally' by an earlier version of the
571  Symmetries code, but that quickly became ungainly and difficult to maintain.
572  This approach makes for cleaner code, and does the additional work of finding
573  e.g. centres of rotations (which may lie outside either or both polygons).
574  """
575
576  shape:geom.Polygon
577  """The target polygon on to which matching transforms are required."""
578  s1:Symmetries
579  """The Symmetries of the target polygon."""
580
581  def __init__(  # noqa: D107
582      self,
583      shape: geom.Polygon) -> None:
584    self.shape = shape
585    self.s1 = Symmetries(shape)
586
587
588  def get_polygon_matches(
589      self,
590      shape2:geom.Polygon) -> list[Transform]|None:
591    """Return all transforms that match shape2 onto this instance's polygon.
592
593    Args:
594      shape2 (geom.Polygon): polygon for which matching transforms are
595        requested.
596
597    Returns:
598      list[Transform]|None: a list of Transforms that will match shape2
599        onto this instance's shape. None is returned if no matching transforms
600        are found.
601
602    """
603    s2 = Symmetries(shape2)
604    if self.s1.symmetry_group != s2.symmetry_group:
605      return None
606    match_rot, rots = self._get_rotation_matches(s2)
607    match_ref, refs = self._get_reflection_matches(s2)
608    if match_rot and match_ref:
609      return rots + refs
610    if match_rot:
611      return rots
612    if match_ref:
613      return refs
614    return None
615
616
617  def _get_rotation_matches(
618      self,
619      s2:Symmetries) -> tuple[bool,list[Transform]]:
620    """Get rotations to match s2's polygon onto this instance's polygon.
621
622    Args:
623      s2 (Symmetries): Symmetries associated with the polygon of interest.
624
625    Returns:
626      tuple[bool,list[Transform]]: boolean is True if there are any matches
627        when matching transforms will be included in the associated list. If
628        there are no matches (False, []) will be returned.
629
630    """
631    # need this locally to check for equality of all members of a list
632    def all_equal(lst:list[float]) -> bool:
633      return all(np.isclose(lst[0], x, 1e-3, 1e-3) for x in lst)
634
635    matches = self.s1.matcher.find_matches(s2.poly_code)
636    if len(matches) == 0:
637      return False, None
638    transforms = []
639    # get lists of polygon corners aligned correctly to measure the angle
640    p1_corners = tiling_utils.get_corners(self.shape, repeat_first = False)
641    p2_corners = tiling_utils.get_corners(s2.polygon, repeat_first = False)
642    corner_distances = self._get_corner_distances_between_polys(
643      p1_corners, p2_corners, matches)
644    equalities = [all_equal(d) for d in corner_distances]
645    if all(equalities):
646      # the two shapes are identical and in the same spot
647      # so rotations symmetries are those of the target shape
648      transforms.extend(self.s1.get_rotations(matches))
649      return True, transforms
650    if any(equalities):
651      # the one where all are equal is a translation so we'll skip it below
652      # and add it here as a translation transform
653      m = equalities.index(True)
654      del matches[m]
655      dx, dy = (p1_corners[0].x - p2_corners[m].x,
656                p1_corners[0].y - p2_corners[m].y)
657      transforms.append(
658        Transform("translation", None, None, (dx, dy), (1, 0, 0, 1, dx, dy)))
659    for m in matches:
660      a, c = self._get_angle_between_polygons(p1_corners, p2_corners, m)
661      if a is None:
662        continue
663      transforms.append(Transform("rotation", a, c, (0, 0),
664        tiling_utils.get_rotation_transform(a, (c.x, c.y))))
665    return True, transforms
666
667
668  def _get_corner_distances_between_polys(
669      self,
670      corners1:list[geom.Point],
671      corners2:list[geom.Point],
672      offsets:list[int],
673    ) -> list[list[float]]:
674    return [
675      [p1.distance(p2) for p1, p2 in
676       zip(corners1, corners2[i:] + corners2[:i], strict = True)]
677      for i in offsets]
678
679
680  def _get_angle_between_polygons(
681      self,
682      corners1:list[geom.Point],
683      corners2:list[geom.Point],
684      offset:int,
685    ) -> tuple[None,...]|tuple[str,tuple[float,...]]|tuple[float,geom.Point]:
686    """Return angle by which lists of corners are rotationally offset.
687
688    If lists are at the same orientation then the translation required to match
689    them (if any) will be returned.
690
691    Args:
692      corners1 (list[geom.Point]): list of corners of the first polygon.
693      corners2 (list[geom.Point]): list of corners of the second polygon.
694      offset (int): number of clockwise steps through the corners by which
695        corners2 are offset from corners1 as previously determined.
696
697    Returns:
698      tuple[float,geom.Point]|tuple[str,tuple[float]]: one of
699      ('identity' (0, 0)), ('translation', (x, y)), or (angle, geom.Point),
700        where the last of these is an angle and a centre of rotation.
701
702    """
703    corners2 = corners2[offset:] + corners2[:offset]
704    dists = [
705      p1.distance(p2) for p1, p2 in zip(corners1, corners2, strict = True)]
706    if all(np.isclose(dists[0], d, atol = 1e-3, rtol = 1e-3) for d in dists):
707      # polygons already aligned similarly so we just need to find translation
708      if np.isclose(dists[0], 0, atol = 1e-3, rtol  = 1e-3):
709        return "identity", (0, 0)
710      if offset == 0:
711        return "translation", (corners1[0].x - corners2[0].x,
712                               corners1[0].y - corners2[0].y)
713    # Polygons not lined up. Order them by distances between them
714    ordered_dists = sorted([(i, d) for i, d in enumerate(dists)],
715                           key = lambda x: x[1], reverse = True)
716    # the two furthest will give the most reliable numbers
717    a, b = [i for i, d in ordered_dists[:2]]
718    p_1a, p_1b = corners1[a], corners1[b]
719    p_2a, p_2b = corners2[a], corners2[b]
720    # get the perpendiculars of the lines joining A and B in each polygon
721    perp_1a_2a = tiling_utils.get_straight_line(
722      p_1a, p_2a, perpendicular = True)
723    perp_1b_2b = tiling_utils.get_straight_line(
724      p_1b, p_2b, perpendicular = True)
725    # and intersect them
726    centre = tiling_utils.get_intersection(perp_1a_2a, perp_1b_2b)
727    if centre is None:
728      return None, None
729    angle = -tiling_utils.get_inner_angle(p_1a, centre, p_2a)
730    if angle < -179:
731      angle = angle + 360
732    elif angle > 180:
733      angle = angle - 360
734    return angle, centre
735
736
737  def _get_reflection_matches(
738      self,
739      s2:Symmetries,
740    ) -> tuple[bool,list[Transform]|None]:
741    """Get reflections that match s2's polygon to this instance's polygon.
742
743    Args:
744      s2 (Symmetries): Symmetries associated with the polygon of interest.
745
746    Returns:
747      tuple[bool,list[Transform]]: boolean is True if there are any matches
748        when matching transforms will be included in the associated list. If
749        there are no matches (False, []) will be returned.
750
751    """
752    matches = self.s1.matcher.find_matches(s2.poly_code_r)
753    if len(matches) == 0:
754      return False, None
755    # get mean centroid - all reflection matches will pass through this
756    c1 = self.s1.polygon.centroid
757    c2 = s2.polygon.centroid
758    c = ((c1.x + c2.x) / 2, (c1.y + c2.y) / 2)
759    # get each polygon's reflection symmetries
760    reflections1 = self.s1.get_reflections(matches)
761    reflections2 = s2.get_reflections(matches)
762    # their mean angle will be a line of reflection through c
763    angles = [(ref1.angle + ref2.angle) / 2
764              for ref1, ref2 in zip(reflections1, reflections2, strict = True)]
765    transforms = [tiling_utils.get_reflection_transform(angle, c)
766                  for angle in angles]
767    # there may also be a 'glide' component additional to the reflection
768    # to find it reflect the centre of polygon 2 by each transform and
769    # then find translation necessary to map it onto the centre of polygon 1
770    ctrs2r = [affine.affine_transform(c2, transform)
771              for transform in transforms]
772    dxdys = [(c1.x - ctr2r.x, c1.y - ctr2r.y) for ctr2r in ctrs2r]
773    return (True,
774      [Transform(
775        "reflection", angle, geom.Point(c), dxdy,
776        tiling_utils.combine_transforms([
777          tiling_utils.get_reflection_transform(angle, c),
778          tiling_utils.get_translation_transform(dxdy[0], dxdy[1])]))
779          for angle, dxdy in zip(angles, dxdys, strict = True)])