bezier update

klayout_dot_config/python/SiEPIC/utils/geometry.py

@@ -61,8 +61,12 @@ def __str__(self):
         return str(Point({self.x}, {self.y}))
 
     def norm(self):
+        '''Euclidean length'''
         return sqrt(self.x**2 + self.y**2)
 
+    def long_edge_length(self):
+        '''return the longest segment of a Manhattan distance'''
+        return max(self.x, self.y)
 
 class Line(Point):
     """ Defines a line """
@@ -150,6 +154,13 @@ def max_curvature(P0, P1, P2, P3):
     max_curv = np.max(np.abs(curv.flatten()))
     return max_curv
 
+def min_curvature(P0, P1, P2, P3):
+    """Gets the minimum curvature of Bezier curve"""
+    t = np.linspace(0, 1, 300)
+    curv = curvature_bezier(P0, P1, P2, P3)(t)
+    min_curv = np.min(np.abs(curv.flatten()))
+    return min_curv
+
 
 def _curvature_penalty(P0, P1, P2, P3):
     """Penalty on the curvyness of Bezier curve"""
@@ -377,27 +388,31 @@ def bezier_optimal(P0, P3, *args, **kwargs):
 except ImportError:
     pass
 
-def bezier_cubic(P0, P3, angle0, angle3, a, b, accuracy = 0.001, *args, **kwargs):
+def bezier_cubic(P0, P3, angle0, angle3, a, b, accuracy = 0.001, verbose=False, plot=False, *args, **kwargs):
     '''
     Calculate a cubic Bezier curve between Points P0 and P3, 
     where the control point positions P1 and P2 are determined by
     the angles at P0 (angle0) and P3 (angle3), 
     at a distance of a * scale from P0, and b * scale from P3,
-    where scale is the distance between P0 and P3.
+    where scale is the longest segment in a Manhattan route between P0 and P3.
 
     Args:
         P0, P3: pya.DPoint (in microns)
         angle0, angle3: radians
-        a, b: <float> greater than 0
+        a, b: <float>. 0 corresponds to P1=P0, and 1 corresponds to P1 at the corner of a 90º bend
         accuracy: 0.001 = 1 nm
 
     Returns:
         list of pya.DPoint
+
+    Example:
+        Bezier curve can approximate a 1/4 circle (arc) for a=b=0.553
+            # https://stackoverflow.com/questions/1734745/how-to-create-circle-with-b%C3%A9zier-curves
     '''
 
     P0 = Point(P0.x, P0.y)
     P3 = Point(P3.x, P3.y)
-    scale = (P3 - P0).norm()  # distance between the two end points
+    scale = (P3 - P0).long_edge_length()  # longest distance between the two end points
     P1 = a * scale * Point(np.cos(angle0), np.sin(angle0)) + P0
     P2 = P3 - b * scale * Point(np.cos(angle3), np.sin(angle3))
     new_bezier_line = bezier_line(P0, P1, P2, P3)
@@ -415,6 +430,25 @@ def bezier_cubic(P0, P3, angle0, angle3, a, b, accuracy = 0.001, *args, **kwargs
         np.insert(bezier_point_coordinates_sampled, -1, bezier_point_coordinates(1 - accuracy / scale),
                     axis=1)  # add a point right before the last one
 
+    if verbose:
+        # print the minimum/maximum curvature
+        print ('SiEPIC.utils.geometry.bezier_cubic: minimum radius of curvature = %0.3g' % (1/max_curvature(P0, P1, P2, P3)))
+        print ('SiEPIC.utils.geometry.bezier_cubic: maximum radius of curvature = %0.3g' % (1/min_curvature(P0, P1, P2, P3)))
+    if plot:
+        t = np.linspace(0, 1, 300)
+        curv = curvature_bezier(P0, P1, P2, P3)(t)
+        rc = 1./curv.flatten()
+        import matplotlib.pyplot as plt
+        plt.plot(t, rc, '--pb', label='a=%3g, b=%3g' % (a,b), linewidth=1.5)
+        SizeFont = 19
+        plt.xlabel('Position along path (t)', fontsize=SizeFont)
+        plt.ylabel('Radius of curvature (microns)', fontsize=SizeFont)
+        plt.legend(fontsize=SizeFont)
+        plt.xticks(fontsize=SizeFont)
+        plt.ylim(bottom=0)
+        plt.yticks(fontsize=SizeFont)
+        plt.show()
+
     return [pya.DPoint(x, y) for (x, y) in zip(*(bezier_point_coordinates_sampled))]
 
 

klayout_dot_config/tech/GSiP/pymacros/tests/test_bezier.py

@@ -60,14 +60,31 @@ def test_bezier_bends():
     h = 1
     width = 0.5
     layer = ly.layer(ly.TECHNOLOGY['Waveguide'])
+    layer2 = ly.layer(ly.TECHNOLOGY['FloorPlan'])
+    print(' Layers: %s, %s' % (ly.TECHNOLOGY['Waveguide'], ly.TECHNOLOGY['FloorPlan']))
 
-    # Two S-bends
+
+    # Two S-bends, small ∆h
     wg_Dpts = bezier_parallel(pya.DPoint(0, 0), pya.DPoint(w, h), 0)
+    wg_polygon = pya.DPolygon(translate_from_normal(wg_Dpts, width/2) +
+                             translate_from_normal(wg_Dpts, -width/2)[::-1]
+                             )
+    cell.shapes(layer2).insert(wg_polygon)
+    a = 0.45
+    wg_Dpts = bezier_cubic(pya.DPoint(0, 0), pya.DPoint(w, h), 0, 0, a, a, accuracy = accuracy)
     wg_polygon = pya.DPolygon(translate_from_normal(wg_Dpts, width/2) +
                              translate_from_normal(wg_Dpts, -width/2)[::-1]
                              )
     cell.shapes(layer).insert(wg_polygon)
-    a = 0.5
+
+    # Two S-bends, large ∆h
+    h = 9
+    wg_Dpts = bezier_parallel(pya.DPoint(0, 0), pya.DPoint(w, h), 0)
+    wg_polygon = pya.DPolygon(translate_from_normal(wg_Dpts, width/2) +
+                             translate_from_normal(wg_Dpts, -width/2)[::-1]
+                             )
+    cell.shapes(layer2).insert(wg_polygon)
+    a = 0.75
     wg_Dpts = bezier_cubic(pya.DPoint(0, 0), pya.DPoint(w, h), 0, 0, a, a, accuracy = accuracy)
     wg_polygon = pya.DPolygon(translate_from_normal(wg_Dpts, width/2) +
                              translate_from_normal(wg_Dpts, -width/2)[::-1]
@@ -76,20 +93,28 @@ def test_bezier_bends():
 
     # 90º bend
     r = 5
-    bezier = 0.3
-    a = (1-bezier)/np.sqrt(2)
-    wg_Dpts = bezier_cubic(pya.DPoint(0, 0), pya.DPoint(r, r), 0, 90/180*np.pi, a, a, accuracy = accuracy)
+    bezier = 0.2
+    a = (1-bezier) # /np.sqrt(2)
+    wg_Dpts = bezier_cubic(pya.DPoint(0, 0), pya.DPoint(r, r), 0, 90/180*np.pi, a, a, accuracy = accuracy, verbose=True)
     wg_polygon = pya.DPolygon(translate_from_normal(wg_Dpts, width/2) +
                              translate_from_normal(wg_Dpts, -width/2)[::-1]
                              )
     cell.shapes(layer).insert(wg_polygon)
     wg_pts = arc_bezier(r/dbu, 0, 90, float(bezier))
     wg_Dpts = [p.to_dtype(dbu) for p in wg_pts]
-                    # wg_pts += Path(arc_xy(-pt_radius, pt_radius, pt_radius, 270, 270 + inner_angle_b_vectors(pts[i-1]-pts[i], pts[i+1]-pts[i]), 
-                    #     DevRec='DevRec' in layers[lr], dbu=dbu), 0).transformed(Trans(angle, turn < 0, pts[i])).get_points()
     wg_polygon = pya.DPolygon(translate_from_normal(wg_Dpts, width/2) +
                              translate_from_normal(wg_Dpts, -width/2)[::-1]
                              ).transformed(pya.Trans(r,0))
+    cell.shapes(layer2).insert(wg_polygon)
+
+    # U-turn bend
+    r = 5
+    h = 1.6 *r 
+    a = 2 / (h/r)
+    wg_Dpts = bezier_cubic(pya.DPoint(0, 0), pya.DPoint(0, h), 0, 180/180*np.pi, a, a, accuracy = accuracy, verbose=True)
+    wg_polygon = pya.DPolygon(translate_from_normal(wg_Dpts, width/2) +
+                             translate_from_normal(wg_Dpts, -width/2)[::-1]
+                             )
     cell.shapes(layer).insert(wg_polygon)
 
 
@@ -204,6 +229,6 @@ def test_bezier_tapers():
 
 
 if __name__ == "__main__":
-    test_bezier_bends()
+    # test_bezier_bends()
     test_bezier_tapers()