movie of the convex hull added on top of the contour

convergence_extension/perimeter_long_axis.py

@@ -4,30 +4,36 @@ import pandas as pd
 from matplotlib.widgets import Slider
 from scipy.spatial import ConvexHull
 
-def perimeter_area(X, Y):
+def perimeter_area_majoraxis(X, Y):
+
     points = np.column_stack((X, Y))
     hull = ConvexHull(points)
+
+    # perimeter
     N = len(hull.vertices)
     perimeter_value = 0
     for i in range(N):
         perimeter_value += np.linalg.norm(points[hull.vertices[i]] - points[hull.vertices[(i+1)%N]]) 
         # The (i+1)%N part ensures that after the last vertex, we go back to the first vertex, creating a closed loop
-    return (perimeter_value, hull.volume)
     
-def long_axis(x, y):
+    # area
+    area = hull.volume
+
+    # long axis
+    x = hull.points[hull.vertices,0]  # hull.vertices gives the index of points on the hull
+    y = hull.points[hull.vertices,1]
     x_diff = x[:, None] - x
     y_diff = y[:, None] - y
-    distances = np.sqrt(x_diff**2 + y_diff**2)
-    return np.max(distances)
+    max_distance = np.max(np.sqrt(x_diff**2 + y_diff**2))
+    return (perimeter_value, area, max_distance, hull)
 
-# NOTE: check perimeter and long_axis functions on known shapes
+# NOTE: check on known shapes
 # x = np.array([1, -1, -1, 1])
 # y = np.array([1, 1, -1, -1])
 
 # points = np.column_stack((x, y))
 # hull = ConvexHull(points)
-# print(perimeter_area(x, y))
-# print(long_axis(x, y))
+# print(perimeter_area_majoraxis(x, y))
 # raise SystemExit
 
 file = np.loadtxt('xy_coordinates_movie1.txt')
@@ -52,19 +58,32 @@ y = np.array(y.interpolate(method='linear', axis=0, limit_direction='both'))
 x = x.T
 y = y.T
 
-# make movie of the contours
+# make movie of the contours and convex hull
 fig, ax = plt.subplots(1, 1, sharex=True, figsize=(8,8))
 ax.set_xlabel(r'$x$')
 ax.set_xlim(np.min(x), np.max(x))
 ax.set_ylim(np.min(y), np.max(y))
+
 ax.set_ylabel(r'$y$')
 
-outlineplot, = ax.plot(x[0], y[0], marker = 'o', ms=3, linestyle='None')
+outlineplot, = ax.plot(x[0], y[0], marker = '.', linestyle='None')
+hull = perimeter_area_majoraxis(x[0], y[0])[-1]
+hull_lines = [ax.plot(x[0][simplex], y[0][simplex], 'r-')[0] for simplex in hull.simplices]
 
 def update(value):
     ti = np.abs(time_array - value).argmin()
     outlineplot.set_ydata(y[ti])
     outlineplot.set_xdata(x[ti])
+    
+    if hasattr(update, 'hull_lines'):
+        for line in update.hull_lines:
+            line.remove()
+    else:  
+        for line in hull_lines:
+            line.remove()
+
+    hull = perimeter_area_majoraxis(x[ti], y[ti])[-1]
+    update.hull_lines = [ax.plot(x[ti][simplex], y[ti][simplex], 'r-')[0] for simplex in hull.simplices]
     plt.draw()
         
 sax = plt.axes([0.1, 0.92, 0.7, 0.02])
@@ -75,28 +94,23 @@ slider.drawon = False
 slider.on_changed(update)
 plt.show()
 
-# find the long axis
-long_axis_array = np.zeros(int(time_points/2))
-for k in range(0, int(time_points/2)):
-    long_axis_array[k] = long_axis(x[k], y[k])
-plt.plot(time_array, long_axis_array, marker='o', linestyle='None')
-plt.xlabel('Time')
-plt.ylabel('Length of the long axis')
-plt.show()
-
 # find the perimeter and area
 perimeter_array = np.zeros(int(time_points/2))
 area_array = np.zeros(int(time_points/2))
+long_axis_array = np.zeros(int(time_points/2))
 
 for k in range(0,  int(time_points/2)):
-    perimeter_array[k], area_array[k] = perimeter_area(x[k], y[k])   
+    perimeter_array[k], area_array[k], long_axis_array[k], _ = perimeter_area_majoraxis(x[k], y[k])   
 
 plt.plot(time_array, perimeter_array, marker='o', linestyle='None')     
 plt.xlabel('Time (min)')
-plt.ylabel('Perimeter ($\mu m$)')
+plt.ylabel('Perimeter of the convex hull ($\mu m$)')
+plt.show()
+plt.plot(time_array, long_axis_array, marker='o', linestyle='None')     
+plt.xlabel('Time (min)')
+plt.ylabel('Long axis of the convex hull ($\mu m$)')
 plt.show()
-
 plt.plot(time_array, area_array, marker='o', linestyle='None')     
 plt.xlabel('Time (min)')
-plt.ylabel('Area ($\mu m^2$)')
+plt.ylabel('Area of the convex hull ($\mu m^2$)')
 plt.show()