moebius.py

   1
   2 from __future__ import print_function
   3
   4 import numpy as np
   5 from numpy import cos, sin
   6
   7 from bezier import BezierSegment
   8 from helixish import HelixishCurve
   9 from moenp import *
  10
  11 import sys
  12
  13 class DoubleCubicBezier():
  14   def __init__(db, cp):
  15     single = BezierSegment(cp)
  16     midpoint = np.array(single.point_at_t(0.5))
  17     mid_dirn = single.point_at_t(0.5 + 0.001) - midpoint
  18     mid_dirn /= np.linalg.norm(mid_dirn)
  19     ocp_factor = 0.5
  20     mid_scale = ocp_factor * 0.5 * (np.linalg.norm(cp[1] - cp[0]) +
  21                                     np.linalg.norm(cp[3] - cp[2]))
  22     db.b0 = BezierSegment([ cp[0],
  23                             cp[1] * ocp_factor + cp[0] * (1-ocp_factor),
  24                             midpoint - mid_dirn * mid_scale,
  25                             midpoint ])
  26     db.b1 = BezierSegment([ midpoint,
  27                             midpoint + mid_dirn * mid_scale,
  28                             cp[2] * ocp_factor + cp[3] * (1-ocp_factor),
  29                             cp[3] ])
  30   def point_at_t(db, t):
  31     if t < 0.5:
  32       return db.b0.point_at_t(t*2)
  33     else:
  34       return db.b1.point_at_t(t*2 - 1)
  35
  36 class ParametricCircle:
  37   def __init__(pc, c, r0, r1):
  38     ''' circle centred on c
  39         with theta=0 point at c+r0
  40         and with theta=tau/4 point at c+r1 '''
  41     pc._c  = c
  42     pc._r0 = r0
  43     pc._r1 = r1
  44   def radius(pc, theta):
  45     return pc._r0 * cos(theta) + pc._r1 * sin(theta)
  46   def point(pc, theta):
  47     return pc._c + pc.radius(theta)
  48
  49 class Twirler(ParametricCircle):
  50   def __init__(tw, c, r0, r1, cycles, begin_zeta):
  51     ''' circle centred on c, etc.
  52         but with an orientation at each point, orthogonal to
  53           the circle
  54         the orientation circles round cycles times during the
  55           whole cycle (if cycles<1, to get a whole circling of
  56           the dirn around the circle, must pass some theta>tau)
  57         begin_zeta is the angle from outwards at theta==0
  58           positive meaning in the direction of r0 x r1 from r0
  59     '''
  60     ParametricCircle.__init__(tw, c, r0, r1)
  61     tw._cycles = cycles
  62     tw._begin_zeta = begin_zeta
  63     tw._axis = np.cross(r0, r1)
  64   def dirn(tw, theta, extra_zeta=0):
  65     ''' unit vector for dirn at theta,
  66          but rotated extra_theta more around the circle
  67     '''
  68     zeta = tw._begin_zeta + theta * tw._cycles + extra_zeta
  69     r = tw.radius(theta)
  70     return cos(zeta) * r + sin(zeta) * tw._axis
  71
  72 class MoebiusHalf:
  73   def __init__(m, nu):
  74     '''
  75      MoebiusHalf().edge is a Twirler for the edge,
  76       expecting theta = u * tau (see MoebiusHalf().point)
  77       with dirn pointing into the surface
  78     '''
  79     m.edge    = Twirler(origin,  unit_z, unit_x, -2, tau/2)
  80     m.midline = Twirler(-unit_z, unit_z, unit_y, -0.5, 0)
  81     m.nu      = nu
  82     m._thetas = [ u * tau for u in np.linspace(0, 1, nu+1) ]
  83     m._cp2b = BezierSegment([ (c,) for c in [0.33,0.33, 1.50]])
  84     m._beziers = [ m._bezier(theta) for theta in m._thetas ]
  85     check = int(nu/3)-1
  86     m._bezier(m._thetas[check], HelixishCurve)
  87     m._bezier(m._thetas[check+1], HelixishCurve)
  88     m._bezier(m._thetas[check+2], HelixishCurve)
  89   def _bezier(m, theta, constructor=DoubleCubicBezier):
  90     cp = [None] * 4
  91     cp[0] =               m.edge   .point(theta)
  92     cp[3] =               m.midline.point(theta*2)
  93     ncp3 = np.linalg.norm(cp[3])
  94     cpt = ncp3 * ncp3 / 4
  95     cp2scale = m._cp2b.point_at_t(ncp3/2)
  96     cp1scale = 1.0 * cpt + 0.33 * (1-cpt)
  97     #print('u=%d ncp3=%f cp2scale=%f' % (u, ncp3, cp2scale),
  98     #      file=sys.stderr)
  99     cp[1] = cp[0] + cp1scale * m.edge   .dirn (theta)
 100     cp[2] = cp[3] + cp2scale * m.midline.dirn (theta*2)
 101     return constructor(cp)
 102   def point(m, iu, t):
 103     '''
 104     0 <= iu <= nu     meaning 0 <= u <= 1
 105                       along the extent (well, along the edge)
 106                        0 and 1 are both the top half of the flat traverse
 107                        0.5        is the bottom half of the flat traverse
 108     0 <=  t   <= 1    across the half-traverse
 109                        0 is the edge, 1 is the midline
 110     '''
 111     return np.array(m._beziers[iu].point_at_t(t))
 112
 113   def details(m, iu, t):
 114     '''
 115     returns tuple of 4 vectors:
 116            - point location
 117            - normal (+ve is in the +ve y direction at iu=t=0) unit vector
 118            - along extent   (towrds +ve iu)                   unit vector
 119            - along traverse (towards +ve t)                   unit vector
 120     '''
 121     if iu == m.nu:
 122       return m.details(0, t)
 123     p = m.point(iu, t)
 124     vec_t = unit_v( m.point(iu, t + 0.01) - p )
 125     if t == 0:
 126       normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
 127       vec_u = np.cross(vec_t, normal)
 128     else:
 129       vec_u = unit_v( m.point(iu+1, t) - p )
 130       normal = np.cross(vec_u, vec_t)
 131     return p, normal, vec_u, vec_t
 132
 133   def point_offset(m, iu, t, offset):
 134     '''
 135     offset by offset perpendicular to the surface
 136      at the top (iu=t=0), +ve offset is in the +ve y direction
 137     '''
 138     p, normal, dummy, dummy = m.details(iu, t)
 139     return p + offset * normal
 140
 141 class Moebius():
 142   def __init__(m, nv, nw):
 143     '''
 144       0 <= v <= nw    along the extent,    v=0 is the flat traverse
 145       0 <= w <= nv    across the traverse  nw must be even
 146       the top is both   v=0, w=0  v=nv, w=nw
 147     '''
 148     assert(nw % 1 == 0)
 149     m.nv = nv
 150     m.nw = nw
 151     m.nt = nw/2
 152     m._t_vals = np.linspace(0, 1, m.nt+1)
 153     m.h = MoebiusHalf(nu=nv*2)
 154
 155   def _vw2tiu_kw(m, v, w):
 156     if w <= m.nt:
 157       it = w
 158       iu = v
 159     else:
 160       it = m.nw - w
 161       iu = m.nv + v
 162     #print('v,w=%d,%d => it,iu=%d,%d' % (v,w,it,iu),
 163     #      file=sys.stderr)
 164     return { 't': m._t_vals[it], 'iu': iu }
 165
 166   def point(m, v, w):
 167     return m.h.point(**m._vw2tiu_kw(v,w))
 168
 169   def point_offset(m, v, w, offset):
 170     return m.h.point_offset(offset=
 171                             offset if w <= m.nt else -offset,
 172                             **m._vw2tiu_kw(v,w))
 173
 174   def details(m, v, w):
 175     '''
 176     returns tuple of 4 vectors:
 177            - point location
 178            - normal (+ve is in the +ve y direction at iu=t=0) unit vector
 179            - along extent   (towrds +ve v)                    unit vector
 180            - along traverse (towards +ve w)                   unit vector
 181     '''
 182     p, normal, vec_v, vec_w = m.h.details(**m._vw2tiu_kw(v,w))
 183     if w > m.nt:
 184       normal = -normal
 185       vec_w  = -vec_w
 186     return p, normal, vec_v, vec_w