merge

.gitignore

@@ -1,3 +1,4 @@
 env/
 .pytest_cache/
-**/__pycache__/
+**/__pycache__/
+example/

README.md

@@ -1 +1,63 @@
-# shape-similarity-
+<h2 align="center">
+  SHAPE SIMILARITY
+</h2>
+
+## :bulb: About
+The package allows you to check the similarity between two shapes/curves, using Frechet distance together with Procrustes analysis.
+Internally, `shape_similarity` works by first normalizing the curves using Procrustes analysis and then calculating Fréchet distance between the curves.
+
+## :page_facing_up: Content
+* Frechet Distance
+  * In mathematics, the Fréchet distance is a measure of similarity between curves that takes into account the location and ordering of the points along the curves. Imagine a person traversing a finite curved path while walking their dog on a leash, with the dog traversing a separate finite curved path. Each can vary their speed to keep slack in the leash, but neither can move backwards. The Fréchet distance between the two curves is the length of the shortest leash sufficient for both to traverse their separate paths from start to finish. Note that the definition is symmetric with respect to the two curves—the Fréchet distance would be the same if the dog were walking its owner.
+
+* Procrustes Analysis
+  * In statistics, Procrustes analysis is a form of statistical shape analysis used to analyse the distribution of a set of shapes. To compare the shapes of two or more objects, the objects must be first optimally "superimposed". Procrustes superimposition (PS) is performed by optimally translating, rotating and uniformly scaling the objects. In other words, both the placement in space and the size of the objects are freely adjusted. The aim is to obtain a similar placement and size, by minimizing a measure of shape difference called the Procrustes distance between the objects.  
+
+## :rocket: Technologies
+* [Python3](https://www.python.org/)
+
+## :package: Installation
+1. Install with pip
+```shell
+$ pip3 install shapesimilarity
+```
+2. Install from source
+```shell
+$ git clone https://github.com/nelsonwenner/shape-similarity.git
+
+$ pip3 install ./shape-similarity
+```
+
+## :information_source: Example useage
+```python
+from shapesimilarity import shape_similarity
+import matplotlib.pyplot as plt
+import numpy as np
+
+x = np.linspace(1, -1, num=200)
+
+y1 = 2*x**2 + 1
+y2 = 2*x**2 + 2
+
+shape1 = np.column_stack((x, y1))
+shape2 = np.column_stack((x, y2))
+
+similarity = shape_similarity(shape1, shape2)
+
+plt.plot(shape1[:,0], shape1[:,1], linewidth=2.0)
+plt.plot(shape2[:,0], shape2[:,1], linewidth=2.0)
+
+plt.title(f'Shape similarity is: {similarity}', fontsize=14, fontweight='bold')
+plt.show()
+```
+
+## :chart_with_downwards_trend: Results
+![export](https://user-images.githubusercontent.com/40550247/126214358-6aa995aa-15b1-4c60-9f0e-34bbef91a99b.png)
+![export](https://user-images.githubusercontent.com/40550247/126214579-302d9220-98ed-4823-992b-d4439145bc5a.png)
+## :pushpin: Referencies
+* [Frechet distance](https://en.wikipedia.org/wiki/Fr%C3%A9chet_distance)
+* [Procrustes analysis](https://en.wikipedia.org/wiki/Procrustes_analysis)
+* [Curve Matcher](https://github.com/chanind/curve-matcher)
+
+---
+Made with :hearts: by Nelson Wenner :wave: [Get in touch!](https://www.linkedin.com/in/nelsonwenner/)

src/__init__.py → shapesimilarity/__init__.py

@@ -1,3 +1,4 @@
 from .procrustesanalysis import *
+from .shapesimilarity import *
 from .frechetdistance import *
 from .geometry import *

shapesimilarity/frechetdistance.py

@@ -0,0 +1,48 @@
+from .geometry import euclidean_distance
+
+'''
+Discrete Frechet distance between 2 curves
+based on http://www.kr.tuwien.ac.at/staff/eiter/et-archive/cdtr9464.pdf
+modified to be iterative and have better memory usage
+'''
+
+def frechet_distance(curve1, curve2):
+  '''
+  Args:
+    polyP: polynomial representing curve 1
+    polyQ: polynomial representing curve 2
+  Returns:
+    Frechet distance between two curves
+  Descriptions:
+    Calculate Frechet distance between two curves
+  '''
+
+  longcalcurve = curve1 if len(curve1) >= len(curve2) else curve2
+  shortcalcurve = curve2 if len(curve1) >= len(curve2) else curve1
+
+  prev_resultscalcol = []
+  for i in range(0, len(longcalcurve)):
+    current_resultscalcol = []
+    for j in range(0, len(shortcalcurve)):
+      current_resultscalcol.append(
+        calc_value(
+          i, j, prev_resultscalcol, 
+          current_resultscalcol, 
+          longcalcurve, shortcalcurve
+        )
+      )
+    prev_resultscalcol = current_resultscalcol
+  return prev_resultscalcol[len(shortcalcurve) - 1]
+
+def calc_value(i, j, prevResultsCol, currentResultsCol, longCurve, shortCurve):
+  if i == 0 and j == 0:
+    return euclidean_distance(longCurve[0], shortCurve[0])
+  if i > 0 and j == 0:
+    return max(prevResultsCol[0], euclidean_distance(longCurve[i], shortCurve[0]))
+  last_result = currentResultsCol[len(currentResultsCol) - 1]
+  if i == 0 and j > 0:
+    return max(last_result, euclidean_distance(longCurve[0], shortCurve[j]))
+  return max(
+    min(prevResultsCol[j], prevResultsCol[j - 1], last_result),
+    euclidean_distance(longCurve[i], shortCurve[j])
+  )

shapesimilarity/geometry.py

@@ -0,0 +1,63 @@
+import math
+
+def euclidean_distance(point1, point2):
+  '''
+  Args:
+    point1: type array two values [x, y]
+    point2: type array two values [x, y]
+  Returns:
+    Distance of two points
+  Descriptions:
+    Calculate Euclidian distance of two points in Euclidian space
+  '''
+
+  return round(math.sqrt((point1[0]-point2[0])**2 + (point1[1]-point2[1])**2), 2)
+
+def curve_length(curve):
+  '''
+  Args:
+    points: type arrays two values [[x, y], [x, y]]
+  Returns:
+    acc_length: curve length
+  Descriptions:
+    Calculate the length of the curve
+  '''
+
+  acc_length = 0
+  for i in range(0, len(curve)-1):
+    acc_length += euclidean_distance(curve[i], curve[i+1])
+  return acc_length
+
+def extend_point_on_line(point1, point2, distance):
+  '''
+  Args:
+    point1: type array two values [x, y]
+    point2: type array two values [x, y]
+    distance: type float
+  Returns:
+    A new point, point3, which is on the same
+    line generated by point1 and point2
+  '''
+
+  norm = round(distance / euclidean_distance(point1, point2), 2)
+  new_point_x = point2[0] + norm * (point1[0] - point2[0])
+  new_point_y = point2[1] + norm * (point1[1] - point2[1])
+  return [new_point_x, new_point_y]
+
+def rotate_curve(curve, thetaRad):
+  '''
+  Args:
+    curve: original curve, type array two values [[x, y], [x, y]]
+    thetaRad: rotation angle type float
+  Returns:
+    rot_curve: rotated curve
+  Descriptions:
+    Rotate the curve around the origin
+  '''
+
+  rot_curve = []
+  for i in range(len(curve)):
+    x_cord = math.cos(-1 * thetaRad) * curve[i][0] - math.sin(-1 * thetaRad) * curve[i][1]
+    y_cord = math.sin(-1 * thetaRad) * curve[i][0] + math.cos(-1 * thetaRad) * curve[i][1]
+    rot_curve.append([x_cord, y_cord])
+  return rot_curve

src/procrustesanalysis.py → shapesimilarity/procrustesanalysis.py

@@ -5,28 +5,28 @@
 '''
 from https://en.wikipedia.org/wiki/Procrustes_analysis
 '''
-def procrustes_normalize_curve(curve, rebalance=True, estimationPoints=50):
+
+def procrustes_normalize_curve(curve):
   '''
   Args:
     curve: type array [[x, y]], [x, y]].
-    rebalance: Optionally runs default True
-    estimationPoints: Optionally runs default 50
   Returns:
     procrustes_normalize_curve: procrustes_normalize_curve
   Descriptions:
     Translate and scale curve by Procrustes Analysis
   '''
-  balanced_curve = rebalance_curve(curve, estimationPoints) if rebalance else curve
-  mean_x = array_average(list(map(lambda item: item[0], balanced_curve)))
-  mean_y = array_average(list(map(lambda item: item[1], balanced_curve)))
-  mean = [mean_x, mean_y]
-  translated_curve = list(map(lambda point: substract(point, mean), balanced_curve))
-  scale = math.sqrt(array_average(list(map(lambda item: item[0]*item[0] + item[1]*item[1], translated_curve))))
-  return list(map(lambda item: [item[0] / scale, item[1] / scale], translated_curve))
+
+  curve_length = len(curve)
+  mean_x = array_average(list(map(lambda item: item[0], curve)))
+  mean_y = array_average(list(map(lambda item: item[1], curve)))
+  curve = list(map(lambda item: [item[0]-mean_x, item[1]-mean_y], curve))
+
+  squared_sum = 0
+  for i in range(curve_length):
+    squared_sum += (curve[i][0])**2 + (curve[i][1])**2
+  scale = round(squared_sum / curve_length, 2)
+  return list(map(lambda item: [item[0]/scale, item[1]/scale], curve))
 
-'''
-from https://en.wikipedia.org/wiki/Procrustes_analysis
-'''
 def find_procrustes_rotation_angle(curve, relativeCurve):
   '''
   Args:
@@ -41,16 +41,14 @@ def find_procrustes_rotation_angle(curve, relativeCurve):
     Find the angle to rotate `curve` to match the rotation
     of `relativeCurve` using procrustes analysis
   '''
+
   assert len(curve) == len(relativeCurve), 'curve and relativeCurve must have the same length'
   numerator, denominator = 0, 0
-  for i in range(0, len(curve)):
-    numerator += curve[i][1] * relativeCurve[i][0] - curve[i][0] * relativeCurve[i][1]
-    denominator += curve[i][0] * relativeCurve[i][0] + curve[i][1] * relativeCurve[i][1]
+  for i in range(len(curve)):
+    numerator += relativeCurve[i][0]*curve[i][1] - relativeCurve[i][1]*curve[i][0]
+    denominator += relativeCurve[i][0]*curve[i][0] + relativeCurve[i][1]*curve[i][1]
   return math.atan2(numerator, denominator)
 
-'''
-from https://en.wikipedia.org/wiki/Procrustes_analysis
-'''
 def procrustes_normalize_rotation(curve, relativeCurve):
   '''
   Args:
@@ -65,5 +63,6 @@ def procrustes_normalize_rotation(curve, relativeCurve):
     Rotate `curve` to match the rotation of 
     `relativeCurve` using procrustes analysis
   '''
+
   angle = find_procrustes_rotation_angle(curve, relativeCurve)
   return rotate_curve(curve, angle)

shapesimilarity/shapesimilarity.py

@@ -0,0 +1,41 @@
+from .procrustesanalysis import *
+from .frechetdistance import *
+from .geometry import *
+import math
+
+def shape_similarity(shape1, shape2, rotations=10, checkRotation=True):
+  procrustes_normalized_curve1 = procrustes_normalize_curve(shape1)
+  procrustes_normalized_curve2 = procrustes_normalize_curve(shape2)
+  geo_avg_curve_len = math.sqrt(
+    curve_length(procrustes_normalized_curve1) * 
+    curve_length(procrustes_normalized_curve2)
+  )
+
+  thetas_to_check = [0]
+  if checkRotation:
+    procrustes_theta = find_procrustes_rotation_angle(
+      procrustes_normalized_curve1, 
+      procrustes_normalized_curve2
+    )
+    # use a negative rotation rather than a large positive rotation
+    if procrustes_theta > math.pi:
+      procrustes_theta = procrustes_theta - 2 * math.pi
+    if procrustes_theta != 0 and abs(procrustes_theta) < math.pi:
+      thetas_to_check.append(procrustes_theta)
+    for i in range(0, rotations):
+      theta = -1 * math.pi + (2 * i * math.pi) / (rotations - 1)
+      # 0 and Math.PI are already being checked, no need to check twice
+      if theta != 0 and theta != math.pi:
+        thetas_to_check.append(theta)
+
+  # Using Frechet distance to check the similarity level 
+  min_frechet_distance = float('inf')
+  # check some other thetas here just in case the procrustes theta isn't the best rotation
+  for theta in thetas_to_check:
+    rotated_curve1 = rotate_curve(procrustes_normalized_curve1, theta)
+    frechet_dist = frechet_distance(rotated_curve1, procrustes_normalized_curve2)
+    if frechet_dist < min_frechet_distance:
+      min_frechet_distance = frechet_dist
+  # divide by Math.sqrt(2) to try to get the low results closer to 
+  result = max(1 - min_frechet_distance / (geo_avg_curve_len / math.sqrt(2)), 0)
+  return round(result, 4)