weavingspace.symmetry

  1#!/usr/bin/env python
  2# coding: utf-8
  3
  4import inspect
  5from typing import Iterable
  6from typing import Any
  7from typing import Union
  8from collections import namedtuple
  9from dataclasses import dataclass
 10import string
 11import copy
 12
 13import numpy as np
 14import matplotlib.pyplot as pyplot
 15import io
 16import PIL
 17
 18import shapely.geometry as geom
 19import shapely.affinity as affine
 20import shapely.ops
 21import geopandas as gpd
 22
 23import weavingspace.tiling_utils as tiling_utils
 24
 25
 26class KMP_Matcher:
 27  """Class to find matching subsequences in a sequence."""
 28
 29  sequence: Iterable
 30  """Iterable in which subsequences are to be found."""
 31
 32  def __init__(self, sequence:Iterable):
 33    self.sequence = sequence
 34
 35  def find_matches(self, pat:Iterable[tuple[float]]) -> list[int]:
 36    """ Implements Knuth-Morris-Pratt string pattern matching algorithm. See:
 37    https://en.wikipedia.org/wiki/Knuth–Morris–Pratt_algorithm which provides
 38    detailed pseudo-code on which this code is directly based. See also:
 39    
 40    Knuth DE, JH Morris Jr, and VR Pratt. 1977. Fast pattern  matching in 
 41    strings. SIAM Journal on Computing 6(2): 323–350. doi: 10.1137/0206024.
 42
 43    This implementation expects sequences of tuples of floats, although in 
 44    principle any objects could be contained in the Iterables.
 45
 46    Args:
 47      pat (Iterable[tuple[float]]): The sequence to match.
 48
 49    Returns:
 50      Iterable[int]: _description_
 51    """
 52    j, k, = 0, 0
 53    finds = []
 54    table = self._get_table(pat)
 55    while j < len(self.sequence):
 56      if self._equal_tuples(pat[k], self.sequence[j]):
 57        j = j + 1
 58        k = k + 1
 59        if k == len(pat):
 60          finds.append(j - k)
 61          k = table[k]
 62      else:
 63        k = table[k]
 64        if k < 0:
 65          j = j + 1
 66          k = k + 1
 67    return finds
 68
 69  def _get_table(self, pattern:Iterable) -> Iterable[int]:
 70    """Returns the 'offsets' table used by KMP pattern matching algorithm as
 71    required by self.find_matches, based on the supplied pattern.
 72
 73    Args:
 74      pattern (Iterable): The pattern to set up the table for.
 75      equals (function, optional): A function to use to test for equality of
 76        elements in patterns. Defaults to _equal_tuples().
 77
 78    Returns:
 79        Iterable[int]: _description_
 80    """
 81    pos = 1
 82    cnd = 0
 83    T = {0: -1}
 84    while pos < len(pattern):
 85      if np.allclose(pattern[pos], pattern[cnd]):
 86        T[pos] = T[cnd]
 87      else:
 88        T[pos] = cnd
 89        while cnd >= 0 and not self._equal_tuples(pattern[pos], pattern[cnd]):
 90          cnd = T[cnd]
 91      pos = pos + 1
 92      cnd = cnd + 1
 93    T[pos] = cnd
 94    return tuple(T.values())
 95
 96  def _equal_tuples(self, t1:Iterable[float], t2:Iterable[float]) -> bool:
 97    """Tests for near equality of two iterables of floats using numpy
 98    allclose. Wrapped so as to allow find_matches function to take 
 99    alternative equality tests as input.
100    """
101    return np.allclose(t1, t2, atol = tiling_utils.RESOLUTION * 10)
102
103  def _round_tuple(self, t:Iterable[float], digits:int = 3) -> Iterable[float]:
104    """Convenience function to round all members of an interable. Used for 
105    display purposes."""
106    return tuple([np.round(x, digits) for x in t])
107
108class Transform:
109  """Class to store details of a transform and draw it."""
110  transform_type: str
111  """Type of transform, 'rotation', 'reflection', 'translation' or 'identity'."""
112  angle: float
113  """Angle of rotation (degrees)."""
114  centre: geom.Point
115  """Centre of the transformation."""
116  translation: tuple[float] 
117  """X and Y coordinates shifts of a translation transform. A glide reflection
118  may also include this."""
119  transform: tuple[float]
120  """Six element tuple for the transform in shapely.transform format. See
121  https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations and
122  methods in `weavingspace.tiling_utils`."""
123  offset: int
124
125  def __init__(self, transform_type:str, angle:float, centre:geom.Point, 
126               translation: tuple[float], transform:tuple[float]):
127    self.transform_type = transform_type
128    self.angle = angle
129    self.centre = centre
130    self.translation = translation
131    self.transform = transform
132
133  def __str__(self) -> str:
134    if self.transform_type == "rotation":
135      return f"{self.transform_type} {self.angle:.1f}° POINT ({self.centre.x:.1f} {self.centre.y:.1f}) {tuple([np.round(x, 3) for x in self.transform])}"
136    elif self.transform_type == "reflection":
137      return f"{self.transform_type} {self.angle:.1f}° {tuple([np.round(x, 3) for x in self.transform])}"
138    else:
139      return f"{self.transform_type} {tuple([np.round(x, 3) for x in self.transform])}"
140
141  def __repr__(self) -> str:
142    return str(self)
143
144  def apply(self, geometry:Any) -> Any:
145    """Applies this transform to supplied shapely geometry and returns result.
146
147    Args:
148      geometry (Any): a shapely geometry to transform.
149    
150    Returns:
151      Any: the resulting transformed shapely geometry.
152    """
153    return affine.affine_transform(geometry, self.transform)
154
155  def draw(self, ax:pyplot.axes, **kwargs) -> pyplot.Axes:
156    """Draws this transform on the supplied axes. Arguments specific to each
157    transform type are supplied as **kwargs and are documented in 
158    `draw_rotation`, `draw_reflection`, and `draw_translation`.
159
160    Args:
161        ax (pyplot.axes): the axes on which to draw the transform.
162
163    Returns:
164        pyplot.Axes: the axes with the rendering of this transform added.
165    """
166    if self.transform_type == "rotation":
167      rotn_args = list(inspect.signature(self.draw_rotation).parameters)
168      rotn_dict = {k: kwargs.pop(k) for k in dict(kwargs) if k in rotn_args}
169      return self.draw_rotation(ax = ax, **rotn_dict)
170    elif self.transform_type == "reflection":
171      refn_args = list(inspect.signature(self.draw_reflection).parameters)
172      refn_dict = {k: kwargs.pop(k) for k in dict(kwargs) if k in refn_args}
173      return self.draw_reflection(ax = ax, **refn_dict)
174    if self.transform_type == "translation":
175      trans_args = list(inspect.signature(self.draw_translation).parameters)
176      trans_dict = {k: kwargs.pop(k) for k in dict(kwargs) if k in trans_args}
177      return self.draw_translation(ax = ax, **trans_dict)
178    elif self.transform_type == "identity":
179      w = ax.get_window_extent().width
180      return ax.set_title(f"{self.transform_type}", fontsize = w / 20)
181
182  def draw_rotation(self, ax:pyplot.Axes, 
183                    radius = 200, add_title = True) -> pyplot.Axes:
184    x, y = self.centre.x, self.centre.y
185    axis = geom.LineString([(x, y), (x + radius * 1.25, y)])
186    arc = geom.LineString([
187      geom.Point(x + radius * np.cos(np.radians(a)),
188                  y + radius * np.sin(np.radians(a)))
189      for a in np.linspace(0, self.angle, 50)])
190    gpd.GeoSeries([self.centre]).plot(
191      ax = ax, color = "r", markersize = 4, zorder = 5)
192    gpd.GeoSeries([axis, arc]).plot(ax = ax, color = "r", lw = .5)
193    if add_title:
194      w = ax.get_window_extent().width
195      ax.set_title(f"{self.transform_type} {np.round(self.angle, 2)}",
196                   fontsize = w / 20)
197    return ax
198
199  def draw_reflection(self, ax:pyplot.Axes, w = 5, 
200                      mirror_length = 100, add_title = True) -> pyplot.Axes:
201    x, y = self.centre.x, self.centre.y
202    dx, dy = self.translation
203    r = np.sqrt(self.translation[0] ** 2 + self.translation[1] ** 2) 
204    mirror = geom.LineString([
205      (x - mirror_length / 2 * np.cos(np.radians(self.angle)), 
206        y - mirror_length / 2 * np.sin(np.radians(self.angle))),
207      (x + mirror_length / 2 * np.cos(np.radians(self.angle)), 
208        y + mirror_length / 2 * np.sin(np.radians(self.angle)))])
209    gpd.GeoSeries([mirror]).plot(
210      ax = ax, color = "r", lw = 1, ls = "dashdot", zorder = 5)
211    no_slide = np.isclose(r, 0, rtol = 1e-6, atol = 1e-6)
212    if not no_slide:
213      pyplot.arrow(
214        x - dx / 2, y - dy / 2, dx, dy, length_includes_head = True,
215        width = w, fc = "k", ec = None, head_width = w * 6, zorder = 5)
216    if add_title:
217      w = ax.get_window_extent().width
218      ax.set_title(f"{self.transform_type} {np.round(self.angle, 2)}",
219                   fontsize = w / 20)
220    return ax
221
222  def draw_translation(self, ax:pyplot.Axes, c:geom.Point, 
223                       w:float = 5, add_title = True) -> pyplot.Axes:
224    gpd.GeoSeries([c]).plot(ax = ax, color = "b")
225    pyplot.arrow(c.x, c.y, self.translation[0], self.translation[1], lw = 0.5,
226                width = w, fc = "b", ec = None, head_width = w * 6, zorder = 5)
227    if add_title:
228      w = ax.get_window_extent().width
229      rounded_tr = tuple([np.round(x, 1) for x in self.translation])
230      ax.set_title(f"{self.transform_type} {rounded_tr}",
231                   fontsize = w / 20)
232    return ax
233
234
235@dataclass
236class Symmetries():
237  """Class to identify and store the symmetries of a supplied shapely.Polygon
238  as a list of `Transform` objects.
239  """
240  polygon:geom.Polygon = None
241  """Polygon from which these symmetries are derived."""
242  matcher:KMP_Matcher = None
243  """the subsequence matcher used in symmetry detection."""
244  n:int = None
245  """number of vertices of the polygon."""
246  p_code:list[tuple[float]] = None
247  """the encoding of the polygon as sequence of length, angle pairs."""
248  p_code_r:list[tuple[float]] = None
249  """the reversed encoding used to detect reflection symmetries."""
250  rotation_shifts:list[int] = None
251  """list of number of 2pi/n rotation symmetries."""
252  reflection_shifts:list[int] = None
253  """list of pi/n relection angle symmetries."""
254  symmetries:list[Transform] = None
255  """list of Transform objects with more complete information."""
256  symmetry_group:str = None
257  """the code denoting the symmetry group"""
258
259  def __init__(self, polygon:geom.Polygon):
260    self.polygon = polygon
261    self.n = shapely.count_coordinates(self.polygon) - 1
262    self.p_code = self._get_polygon_code(self.polygon)
263    self.p_code_r = self._get_polygon_code(self.polygon, mirrored = True)
264    self.matcher = KMP_Matcher(self.p_code + self.p_code[:-1])
265    self.symmetries = self.get_symmetries()
266    self.symmetry_group = self.get_symmetry_group_code()
267
268  def _get_polygon_code(self, p:geom.Polygon, 
269                        mirrored = False) -> list[tuple[float]]:
270    r"""Returns a list of length-angle pairs to uniquely encode a polygon
271    shape (up to scale). The alignment of lengths and angles is important
272    and symmetry detection is sensitively dependent on it...
273
274    The forward (i.e. unmirrored) polygon pairs are of edge lengths and the 
275    angle between each edge and its successor going CW around the polygon. 
276    
277           A[i] ----L[i]---- A (i+1)
278          /                   \
279         /                     \
280                             
281    i.e. edge i is between angles i and i + 1, and angle i is between edge
282    i-1 and i. The unmirrored encoding pairs length[i] with angle[i+1].
283    
284    The mirrored encoding pairs matched indexes of lengths and angles in the
285    original polygon. This means the angle is still the one between and edge 
286    and its successor but proceeding CCW around the polygon.
287    
288    See in particular for clarification.
289    
290    Eades P. 1988. Symmetry Finding Algorithms. In Machine Intelligence and 
291    Pattern Recognition, ed. GT Toussaint, 6:41-51. Computational Morphology. 
292    North-Holland. doi: 10.1016/B978-0-444-70467-2.50009-6.
293    
294    Args:
295      p (geom.Polygon): the polygon to encode
296      mirrored (bool, optional): if true encoding will be in CCC order. 
297        Defaults to False.
298
299    Returns:
300      list[tuple[float]]: _description_
301    """
302    # get edge lengths and angles
303    lengths = tiling_utils.get_side_lengths(p)
304    raw_angles = tiling_utils.get_interior_angles(p)
305    if mirrored:
306      lengths_r = list(reversed(lengths))
307      angles_r = list(reversed(raw_angles))
308      return list(zip(lengths_r, angles_r))
309    else:
310      angles = raw_angles[1:] + raw_angles[:1]
311      return list(zip(lengths, angles))
312
313  def get_symmetries(self) -> list[Transform]:
314    """Finds rotation and reflection symmetries of the supplied polygon.
315    Based on
316    
317      Eades P. 1988. Symmetry Finding Algorithms. In Machine Intelligence and 
318      Pattern Recognition, ed. GT Toussaint, 6:41-51. Computational Morphology. 
319      North-Holland. doi: 10.1016/B978-0-444-70467-2.50009-6.
320    
321    and also
322    
323      Wolter JD, TC Woo, and RA Volz. 1985. Optimal algorithms for symmetry 
324      detection in two and three dimensions. The Visual Computer 1(1): 37-48. 
325      doi: 10.1007/BF01901268.
326    
327    Details in these papers are not unambiguous. This implementation was
328    developed based on them, but with a lot of trial and error to get the index
329    offsets and (especially) retrieval of the reflection axes angles to work.
330
331    Returns:
332      list[Transform]: a list of Transform objects representing the polygon
333        symmetries.
334    """
335    # ==== Rotations ====
336    self.rotation_shifts = self.matcher.find_matches(self.p_code)
337    # ==== Reflections ====
338    self.reflection_shifts = self.matcher.find_matches(self.p_code_r)
339    return self.get_rotations(self.rotation_shifts) + \
340           self.get_reflections(self.reflection_shifts)
341
342  def get_rotations(self, offsets:list[int]) -> list[Transform]:
343    """Gets the rotations associated with this collection of symmetries.
344
345    Returns:
346      list[Transform]: A list of the rotation symmetries associated with this
347        polygon.
348    """
349    c = self.polygon.centroid
350    rot_angles = [np.round(i * 360 / self.n, 6) for i in offsets] 
351    rotation_transforms = [
352      np.round(tiling_utils.get_rotation_transform(a, (c.x, c.y)), 6) 
353      for a in rot_angles]
354    return [Transform("rotation", angle, c, (0, 0), transform)
355            for transform, angle in zip(rotation_transforms, rot_angles)]
356
357  def get_reflections(self, offsets:list[int]) -> list[Transform]:
358    """Gets the reflections associated with this collection of symmetries.
359
360    Returns:
361      list[Transform]: A list of the reflection symmetries associated with this
362        polygon.
363    """
364    # calculate all possible reflection axes - ordering of these is important
365    # for correct picking - don't mess with this code!
366    c = self.polygon.centroid
367    bearings = tiling_utils.get_side_bearings(self.polygon)
368    reflection_axes = []
369    interior_angles = tiling_utils.get_interior_angles(self.polygon)
370    for i in range(self.n):
371      # angle bisectors
372      reflection_axes.append(bearings[i] - interior_angles[i] / 2)
373      # edge_perpendicular bisectors
374      reflection_axes.append(bearings[i] - 90)  # normal to the edge
375    # pick out those where matches occurred
376    ref_angles = [np.round(reflection_axes[i], 6) for i in offsets]
377    reflection_transforms = [
378      np.round(tiling_utils.get_reflection_transform(a, (c.x, c.y)), 6)
379      for a in ref_angles]
380    return [Transform("reflection", angle, c, (0, 0), transform)
381            for transform, angle in zip(reflection_transforms, ref_angles)]
382
383  def get_symmetry_group_code(self):
384    if len(self.reflection_shifts) == 0:
385      return f"C{len(self.rotation_shifts)}"
386    else:
387      return f"D{len(self.rotation_shifts)}"
388
389  def get_corner_offset(self, poly2:geom.Polygon):
390    s2 = Symmetries(poly2)
391    matches = self.matcher.find_matches(s2.p_code)
392    if len(matches) > 0:
393      return matches[0]
394    matches = self.matcher.find_matches(s2.p_code_r)
395    if len(matches) > 0:
396      return -matches[0]
397    return None
398
399  def get_corner_labels(self) -> dict[str, list[str]]:
400    """Returns all the reorderings of vertex labels corresponding to each
401    symmetry.
402
403    Returns:
404      dict[str, list[str]]: A dictionary with two entries. "rotations" is
405        a list of labels under the rotation symmetries, "reflections" those
406        under the reflection symmetries.
407    """
408    labels = list(string.ascii_letters.upper())[:self.n]
409    return self._get_labels_under_symmetries(labels)
410
411  def get_unique_labels(self, offset:int = 0) -> list[str]:
412    labellings = self.get_corner_labels()
413    labellings = labellings["rotations"] + labellings["reflections"]
414    labellings = ["".join(sorted(x)) for x in zip(*labellings)]
415    n_new_labels = len(set(labellings))
416    letters = list(string.ascii_letters.upper())[offset:offset + n_new_labels]
417    mapping = dict()
418    i = 0
419    for label in labellings:
420      if not label in mapping:
421        mapping[label] = letters[i]
422        i = i + 1
423    return self._get_labels_under_symmetries([mapping[x] for x in labellings])
424  
425  def _get_labels_under_symmetries(
426      self, labels:list[str]) -> dict[str, list[str]]:
427    cycle = labels * 2
428    under_rotation = ["".join(cycle[i:self.n + i]) 
429                      for i in self.rotation_shifts]
430    under_reflection = ["".join(cycle[self.n + i:i:-1]) 
431                        for i in self.reflection_shifts]
432    return {"rotations": under_rotation,
433            "reflections": under_reflection}
434
435  def plot(self, as_image:bool = False, title:str = ""):
436    fig = pyplot.figure()
437    fig.suptitle(title)
438
439    n_subplots = len(self.symmetries)
440    nr = int(np.ceil(np.sqrt(n_subplots)))
441    nc = int(np.ceil(n_subplots / nr))
442    n_plots = 0
443
444    for s in self.symmetries:
445      n_plots = n_plots + 1
446      ax = fig.add_subplot(nr, nc, n_plots)
447      s.plot(ax, self.polygon)
448    
449    if as_image:
450      buf = io.BytesIO()
451      fig.savefig(buf)
452      buf.seek(0)
453      return PIL.Image.open(buf)
454    return ax
455
456
457StraightLine = namedtuple("StraightLine", "A B C")
458"""Named tuple to represent equation of a straight line in standard 
459Ax + By + C = 0 form."""
460
461
462class Shape_Matcher:
463  shape: geom.Polygon
464  s1: Symmetries
465  other: geom.Polygon
466  centre: geom.Point
467  translation: tuple[float]
468  matches: list[Transform]
469  identity_transform: tuple[float] = (1, 0, 0, 1, 0, 0)
470  
471  def __init__(self, shape: geom.Polygon):
472    self.shape = shape
473    self.s1 = Symmetries(shape)
474
475  def get_polygon_matches(self, shape2: geom.Polygon):
476    s2 = Symmetries(shape2)
477    if self.s1.symmetry_group != s2.symmetry_group:
478      # print("No matches")
479      return None
480    else:
481      match_rot, rots = self._get_rotation_matches(s2)
482      match_ref, refs = self._get_reflection_matches(s2)
483      if match_rot and match_ref:
484        # print(f"Rotation and reflection matches found")
485        return rots + refs
486      elif match_rot:
487        # print(f"Only rotation matches found")
488        return rots
489      elif match_ref:
490        # print(f"Only reflection matches found")
491        return refs
492      else:
493        # print (f"No matches found")
494        return None
495
496  def _get_rotation_matches(self, s2:Symmetries):
497    matches = self.s1.matcher.find_matches(s2.p_code)
498    if len(matches) == 0:
499      return False, None
500    else:
501      transforms = []
502      # get lists of polygon corners aligned correctly to measure the angle
503      p1_corners = tiling_utils.get_corners(self.shape, repeat_first = False)
504      p2_corners = tiling_utils.get_corners(s2.polygon, repeat_first = False)
505      for m in matches:
506        a, c = self.get_angle_between_polygons(p1_corners, p2_corners, m)
507        if a == "translation":
508          # print("One of the rotation matches is a translation.")
509          transforms.append(
510            Transform("translation", None, None, c, (1, 0, 0, 1, c[0], c[1])))
511        elif a == "identity":
512          transforms.append(
513            Transform("identity", None, None, c, (1, 0, 0, 1, 0, 0)))
514        elif a is None:
515          continue
516        else:
517          transforms.append(Transform("rotation", a, c, (0, 0),
518            tiling_utils.get_rotation_transform(a, (c.x, c.y))))
519      return True, transforms
520
521  def get_angle_between_polygons(self, corners1:list[geom.Point], 
522      corners2:list[geom.Point], offset:int) -> tuple[Union[float,geom.Point]]:
523    corners2 = corners2[offset:] + corners2[:offset]
524    dists = [p1.distance(p2) for p1, p2 in zip(corners1, corners2)]
525    if all([np.isclose(dists[0], d, atol = 1e-3, rtol = 1e-3) for d in dists]):
526      if np.isclose(dists[0], 0, atol = 1e-3, rtol  = 1e-3):
527        return "identity", (0, 0)
528      else:
529        return "translation", (corners1[0].x - corners2[0].x,
530                               corners1[0].y - corners2[0].y)
531    ordered_dists = sorted([(i, d) for i, d in enumerate(dists)], 
532                           key = lambda x: x[1])
533    AB = ordered_dists[-2:]
534    p1A, p1B = corners1[AB[0][0]], corners1[AB[1][0]]
535    p2A, p2B = corners2[AB[0][0]], corners2[AB[1][0]]
536    perpAA = self.get_straight_line(p1A, p2A, True)
537    perpBB = self.get_straight_line(p1B, p2B, True)
538    centre = self.get_intersection(perpAA, perpBB)
539    if centre is None:
540      angle = None
541    else:
542      angle = -tiling_utils.get_inner_angle(p1A, centre, p2A)
543      if angle < -179:
544        angle = angle + 360
545      elif angle > 180:
546        angle = angle - 360
547    return angle, centre
548
549  def get_straight_line(self, p1:geom.Point, p2:geom.Point, 
550                        perpendicular:bool = False) -> StraightLine:
551    if perpendicular:
552      ls = affine.rotate(geom.LineString([p1, p2]), 90)
553      pts = [p for p in ls.coords]
554      x1, y1 = pts[0]
555      x2, y2 = pts[1]
556    else:
557      x1, y1 = p1.x, p1.y
558      x2, y2 = p2.x, p2.y
559    return StraightLine(y1 - y2, x2 - x1, x1 * y2 - x2 * y1)
560
561  def get_intersection(self, line1:StraightLine, 
562                       line2:StraightLine) -> geom.Point:
563    x_set, y_set = False, False
564    denominator = line1.A * line2.B - line2.A * line1.B
565    if np.isclose(line1.A, 0, atol = 1e-4, rtol = 1e-4):
566      y = -line1.C / line1.B
567      y_set = True
568    elif np.isclose(line2.A, 0, atol = 1e-4, rtol = 1e-4):
569      y = -line2.C / line2.B
570      y_set = True
571    if np.isclose(line1.B, 0, atol = 1e-4, rtol = 1e-4):
572      x = -line1.C / line1.A
573      x_set = True
574    elif np.isclose(line2.B, 0, atol = 1e-4, rtol = 1e-4):
575      x = -line2.C / line2.A
576      x_set = True
577    if np.isclose(denominator, 0, atol = 1e-4, rtol = 1e-4):
578      return None
579    x = x if x_set else (line1.B * line2.C - line2.B * line1.C) / denominator
580    y = y if y_set else (line1.C * line2.A - line2.C * line1.A) / denominator
581    return geom.Point(x, y)
582
583  def _get_reflection_matches(self, s2:Symmetries):
584    ctr1 = self.s1.polygon.centroid
585    ctr2 = s2.polygon.centroid
586    c = ((ctr1.x + ctr2.x) / 2, (ctr1.y + ctr2.y) / 2)
587    matches = self.s1.matcher.find_matches(s2.p_code_r)
588    if len(matches) == 0:
589      return False, None
590    reflections1 = self.s1.get_reflections(matches)
591    reflections2 = s2.get_reflections(matches)
592    angles = [(ref1.angle + ref2.angle) / 2
593              for ref1, ref2 in zip(reflections1, reflections2)]
594    trs = [tiling_utils.get_reflection_transform(angle, c)
595           for angle in angles]
596    ctrs2r = [affine.affine_transform(ctr2, tr) for tr in trs]
597    dxdys = [(ctr1.x - ctr2r.x, ctr1.y - ctr2r.y) for ctr2r in ctrs2r]
598    return True, [Transform(
599      "reflection", angle, geom.Point(c), dxdy,
600      tiling_utils.combine_transforms(
601        [tiling_utils.get_reflection_transform(angle, c), 
602         (1, 0, 0, 1, dxdy[0], dxdy[1])])) 
603      for angle, dxdy in zip(angles, dxdys)]