Revert "mesh: introduce is_extt (nfc)"

[moebius3.git] / moebius.py

diff --git a/moebius.py b/moebius.py
index 2bb740c4072ec18deb6b6b4c59b226eaf009b2fd..bd98a825c6b1dac6f6eb1325ba5e06a005d8d580 100644 (file)
--- a/moebius.py
+++ b/moebius.py
@@ -15,6 +15,32 @@ unit_x = np.array((1,0,0))
 unit_y = np.array((0,1,0))
 unit_z = np.array((0,0,1))
 
 unit_y = np.array((0,1,0))
 unit_z = np.array((0,0,1))
 

+def unit_v(v):
+  return v / np.linalg.norm(v)
+
+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):
     ''' circle centred on c
 class ParametricCircle:
   def __init__(pc, c, r0, r1):
     ''' circle centred on c


@@ -34,31 +60,131 @@ 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 (maybe)
+          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) for theta in m._thetas ]
+  def _bezier(m, theta):

     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 DoubleCubicBezier(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