helixish: attempt at the whole thing

[moebius3.git] / moebius.py

diff --git a/moebius.py b/moebius.py
index c6041244c3985b0342aa155e6117d763addd4463..e5a81acad9d0d458bbc823769b5cea55176e0801 100644 (file)
--- a/moebius.py
+++ b/moebius.py
@@ -5,15 +5,33 @@ import numpy as np
 from numpy import cos, sin
 
 from bezier import BezierSegment
 from numpy import cos, sin
 
 from bezier import BezierSegment

+from helixish import HelixishCurve
+from moenp import *

 
 import sys
 
 
 import sys
 

-tau = np.pi * 2
-
-origin = np.array((0,0,0))
-unit_x = np.array((1,0,0))
-unit_y = np.array((0,1,0))
-unit_z = np.array((0,0,1))
+class DoubleCubicBezier():
+  def __init__(db, cp):
+    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]))
+    db.b0 = BezierSegment([ cp[0],
+                            cp[1] * ocp_factor + cp[0] * (1-ocp_factor),
+                            midpoint - mid_dirn * mid_scale,
+                            midpoint ])
+    db.b1 = BezierSegment([ 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)

 
 class ParametricCircle:
   def __init__(pc, c, r0, r1):
 
 class ParametricCircle:
   def __init__(pc, c, r0, r1):


@@ -34,31 +52,133 @@ class Twirler(ParametricCircle):
         but with an orientation at each point, orthogonal to
           the circle
         the orientation circles round cycles times during the
         but with an orientation at each point, orthogonal to
           the circle
         the orientation circles round cycles times during the

-          whole cycle
+          whole cycle (if cycles<1, to get a whole circling of
+          the dirn around the circle, must pass some theta>tau)

         begin_zeta is the angle from outwards at theta==0
         begin_zeta is the angle from outwards at theta==0

-          positive meaning in the direction of r0 x r1
+          positive meaning in the direction of r0 x r1 from r0

     '''
     ParametricCircle.__init__(tw, c, r0, r1)
     tw._cycles = cycles
     tw._begin_zeta = begin_zeta
     tw._axis = np.cross(r0, r1)
     '''
     ParametricCircle.__init__(tw, c, r0, r1)
     tw._cycles = cycles
     tw._begin_zeta = begin_zeta
     tw._axis = np.cross(r0, r1)

-  def dirn(tw, theta):
-    zeta = tw._begin_zeta + theta * tw._cycles
+  def dirn(tw, theta, extra_zeta=0):
+    ''' unit vector for dirn at theta,
+         but rotated extra_theta more around the circle
+    '''
+    zeta = tw._begin_zeta + theta * tw._cycles + extra_zeta

     r = tw.radius(theta)
     return cos(zeta) * r + sin(zeta) * tw._axis
 
     r = tw.radius(theta)
     return cos(zeta) * r + sin(zeta) * tw._axis
 

-class Moebius:
-  def __init__(m, n_u): # ix_u will be in [0, n_u] for [0, 1]
+class MoebiusHalf:
+  def __init__(m, nu):
+    '''
+     MoebiusHalf().edge is a Twirler for the edge,
+      expecting theta = u * tau (see MoebiusHalf().point)
+      with dirn pointing into the surface
+    '''

     m.edge    = Twirler(origin,  unit_z, unit_x, -2, tau/2)
     m.midline = Twirler(-unit_z, unit_z, unit_y, -0.5, 0)
     m.edge    = Twirler(origin,  unit_z, unit_x, -2, tau/2)
     m.midline = Twirler(-unit_z, unit_z, unit_y, -0.5, 0)

-    m._beziers = [ m._bezier(u) for u in np.linspace(0, 1, n_u+1) ]
-  def _bezier(m,u):
-    theta = u * tau
+    m.nu      = nu
+    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, HelixishCurve) for theta in m._thetas ]
+    #check = int(nu/3)
+    #m._beziers[check] = m._bezier(m._thetas[check], HelixishCurve)
+  def _bezier(m, theta, constructor=DoubleCubicBezier):

     cp = [None] * 4
     cp[0] =               m.edge   .point(theta)
     cp = [None] * 4
     cp[0] =               m.edge   .point(theta)

-    cp[1] = cp[0] + 0.5 * m.edge   .dirn (theta)

     cp[3] =               m.midline.point(theta*2)
     cp[3] =               m.midline.point(theta*2)

-    cp[2] = cp[3] + 0.5 * m.midline.dirn (theta*2)
-    return BezierSegment(cp)
-  def point(m, ix_u, t):
-    return m._beziers[ix_u].point_at_t(t)
+    ncp3 = np.linalg.norm(cp[3])
+    cpt = ncp3 * ncp3 / 4
+    cp2scale = m._cp2b.point_at_t(ncp3/2)
+    cp1scale = 1.0 * cpt + 0.33 * (1-cpt)
+    #print('u=%d ncp3=%f cp2scale=%f' % (u, ncp3, cp2scale),
+    #      file=sys.stderr)
+    cp[1] = cp[0] + cp1scale * m.edge   .dirn (theta)
+    cp[2] = cp[3] + cp2scale * m.midline.dirn (theta*2)
+    return constructor(cp)
+  def point(m, iu, t):
+    '''
+    0 <= iu <= nu     meaning 0 <= u <= 1
+                      along the extent (well, along the edge)
+                       0 and 1 are both the top half of the flat traverse
+                       0.5        is the bottom half of the flat traverse
+    0 <=  t   <= 1    across the half-traverse
+                       0 is the edge, 1 is the midline
+    '''
+    return np.array(m._beziers[iu].point_at_t(t))
+
+  def details(m, iu, t):
+    '''
+    returns tuple of 4 vectors:
+           - point location
+           - normal (+ve is in the +ve y direction at iu=t=0) unit vector
+           - along extent   (towrds +ve iu)                   unit vector
+           - along traverse (towards +ve t)                   unit vector
+    '''
+    if iu == m.nu:
+      return m.details(0, t)
+    p = m.point(iu, t)
+    vec_t = unit_v( m.point(iu, t + 0.01) - p )
+    if t == 0:
+      normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
+      vec_u = np.cross(vec_t, normal)
+    else:
+      vec_u = unit_v( m.point(iu+1, t) - p )
+      normal = np.cross(vec_u, vec_t)
+    return p, normal, vec_u, vec_t
+
+  def point_offset(m, iu, t, offset):
+    '''
+    offset by offset perpendicular to the surface
+     at the top (iu=t=0), +ve offset is in the +ve y direction
+    '''
+    p, normal, dummy, dummy = m.details(iu, t)
+    return p + offset * normal
+
+class Moebius():
+  def __init__(m, nv, nw):
+    '''
+      0 <= v <= nw    along the extent,    v=0 is the flat traverse
+      0 <= w <= nv    across the traverse  nw must be even
+      the top is both   v=0, w=0  v=nv, w=nw
+    '''
+    assert(nw % 1 == 0)
+    m.nv = nv
+    m.nw = nw
+    m.nt = nw/2
+    m._t_vals = np.linspace(0, 1, m.nt+1)
+    m.h = MoebiusHalf(nu=nv*2)
+
+  def _vw2tiu_kw(m, v, w):
+    if w <= m.nt:
+      it = w
+      iu = v
+    else:
+      it = m.nw - w
+      iu = m.nv + v
+    #print('v,w=%d,%d => it,iu=%d,%d' % (v,w,it,iu),
+    #      file=sys.stderr)
+    return { 't': m._t_vals[it], 'iu': iu }
+
+  def point(m, v, w):
+    return m.h.point(**m._vw2tiu_kw(v,w))
+
+  def point_offset(m, v, w, offset):
+    return m.h.point_offset(offset=
+                            offset if w <= m.nt else -offset,
+                            **m._vw2tiu_kw(v,w))
+
+  def details(m, v, w):
+    '''
+    returns tuple of 4 vectors:
+           - point location
+           - normal (+ve is in the +ve y direction at iu=t=0) unit vector
+           - along extent   (towrds +ve v)                    unit vector
+           - along traverse (towards +ve w)                   unit vector
+    '''
+    p, normal, vec_v, vec_w = m.h.details(**m._vw2tiu_kw(v,w))
+    if w > m.nt:
+      normal = -normal
+      vec_w  = -vec_w
+    return p, normal, vec_v, vec_w