RecursiveBezierishCurve: Try it

This does not seem to make any difference.  How interesting.  I
suspect that in fact this construction arranges, somehow, to get the
same cubics for the two halves as for the original.  But I CBA to do
the algebra to check.

Ah well, we will revert this and try something else.

Signed-off-by: Ian Jackson <ijackson@chiark.greenend.org.uk>

diff --git a/moebius.py b/moebius.py
index 34111728052e03c6b6dbdcceaf4a88f03d7f0a60..851a8a8029fabbed929648e9089a434fb3af07c9 100644 (file)
--- a/moebius.py
+++ b/moebius.py
@@ -5,6 +5,7 @@ import numpy as np
 from numpy import cos, sin
 
 from bezier import BezierSegment
+from recursivebezier import RecursiveBezierishCurve
 from moenp import *
 
 import sys
@@ -81,7 +82,7 @@ class MoebiusHalf:
     m._thetas = [ u * tau for u in np.linspace(0, 1, nu+1) ]
     m._cp2b = BezierSegment([ (c,) for c in [0.33,0.33, 1.50]])
     m._beziers = [ m._bezier(theta) for theta in m._thetas ]
-  def _bezier(m, theta, constructor=DoubleCubicBezier):
+  def _bezier(m, theta, constructor=RecursiveBezierishCurve):
     cp = [None] * 4
     cp[0] =               m.edge   .point(theta)
     cp[3] =               m.midline.point(theta*2)

diff --git a/recursivebezier.py b/recursivebezier.py
new file mode 100644 (file)
index 0000000..9ec7f27
--- /dev/null
+++ b/recursivebezier.py
@@ -0,0 +1,36 @@
+
+from __future__ import print_function
+
+import numpy as np
+
+from bezier import BezierSegment
+from moenp import *
+
+class RecursiveBezierishCurve():
+  def __init__(db, cp, depth=10):
+    single = BezierSegment(cp)
+    midpoint = np.array(single.point_at_t(0.5))
+    mid_dirn = single.point_at_t(0.5 + 0.001) - midpoint
+    mid_dirn /= np.linalg.norm(mid_dirn)
+    ocp_factor = 0.5
+    mid_scale = ocp_factor * 0.5 * (np.linalg.norm(cp[1] - cp[0]) +
+                                    np.linalg.norm(cp[3] - cp[2]))
+    if depth > 0:
+      constructor = lambda ccp: RecursiveBezierishCurve(ccp, depth-1)
+    else:
+      constructor = BezierSegment
+
+    db.b0 = constructor([ cp[0],
+                          cp[1] * ocp_factor + cp[0] * (1-ocp_factor),
+                          midpoint - mid_dirn * mid_scale,
+                          midpoint ])
+    db.b1 = constructor([ midpoint,
+                          midpoint + mid_dirn * mid_scale,
+                          cp[2] * ocp_factor + cp[3] * (1-ocp_factor),
+                          cp[3] ])
+  def point_at_t(db, t):
+    if t < 0.5:
+      return db.b0.point_at_t(t*2)
+    else:
+      return db.b1.point_at_t(t*2 - 1)
+