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