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):
  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._thetas = [ u * tau for u in np.linspace(0, 1, nu+1) ]
  82     m._cp2b = BezierSegment([ (c,) for c in [0.33,0.33, 1.50]])
  83     m._beziers = [ m._bezier(theta) for theta in m._thetas ]
  84   def _bezier(m, theta):
  85     cp = [None] * 4
  86     cp[0] =               m.edge   .point(theta)
  87     cp[3] =               m.midline.point(theta*2)
  88     ncp3 = np.linalg.norm(cp[3])
  89     cpt = ncp3 * ncp3 / 4
  90     cp2scale = m._cp2b.point_at_t(ncp3/2)
  91     cp1scale = 1.0 * cpt + 0.33 * (1-cpt)
  92     #print('u=%d ncp3=%f cp2scale=%f' % (u, ncp3, cp2scale),
  93     #      file=sys.stderr)
  94     cp[1] = cp[0] + cp1scale * m.edge   .dirn (theta)
  95     cp[2] = cp[3] + cp2scale * m.midline.dirn (theta*2)
  96     return DoubleCubicBezier(cp)
  97   def point(m, iu, t):
  98     '''
  99     0 <= iu <= nu     meaning 0 <= u <= 1
 100                       along the extent (well, along the edge)
 101                        0 and 1 are both the top half of the flat traverse
 102                        0.5        is the bottom half of the flat traverse
 103     0 <=  t   <= 1    across the half-traverse
 104                        0 is the edge, 1 is the midline
 105     '''
 106     return np.array(m._beziers[iu].point_at_t(t))
 107
 108   def details(m, iu, t):
 109     '''
 110     returns tuple of 4 vectors:
 111            - point location
 112            - normal (+ve is in the +ve y direction at iu=t=0) unit vector
 113            - along extent   (towrds +ve iu)                   unit vector
 114            - along traverse (towards +ve t)                   unit vector
 115     '''
 116     if iu == m.nu:
 117       return m.details(0, t)
 118     p = m.point(iu, t)
 119     vec_t = unit_v( m.point(iu, t + 0.01) - p )
 120     if t == 0:
 121       normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
 122       vec_u = np.cross(vec_t, normal)
 123     else:
 124       vec_u = unit_v( m.point(iu+1, t) - p )
 125       normal = np.cross(vec_u, vec_t)
 126     return p, normal, vec_u, vec_t
 127
 128   def point_offset(m, iu, t, offset):
 129     '''
 130     offset by offset perpendicular to the surface
 131      at the top (iu=t=0), +ve offset is in the +ve y direction
 132     '''
 133     p, normal, dummy, dummy = m.details(iu, t)
 134     return p + offset * normal
 135
 136 class Moebius():
 137   def __init__(m, nv, nw):
 138     '''
 139       0 <= v <= nw    along the extent,    v=0 is the flat traverse
 140       0 <= w <= nv    across the traverse  nw must be even
 141       the top is both   v=0, w=0  v=nv, w=nw
 142     '''
 143     assert(nw % 1 == 0)
 144     m.nv = nv
 145     m.nw = nw
 146     m.nt = nw/2
 147     m._t_vals = np.linspace(0, 1, m.nt+1)
 148     m.h = MoebiusHalf(nu=nv*2)
 149
 150   def _vw2tiu_kw(m, v, w):
 151     if w <= m.nt:
 152       it = w
 153       iu = v
 154     else:
 155       it = m.nw - w
 156       iu = m.nv + v
 157     #print('v,w=%d,%d => it,iu=%d,%d' % (v,w,it,iu),
 158     #      file=sys.stderr)
 159     return { 't': m._t_vals[it], 'iu': iu }
 160
 161   def point(m, v, w):
 162     return m.h.point(**m._vw2tiu_kw(v,w))
 163
 164   def point_offset(m, v, w, offset):
 165     return m.h.point_offset(offset=
 166                             offset if w <= m.nt else -offset,
 167                             **m._vw2tiu_kw(v,w))
 168
 169   def details(m, v, w):
 170     '''
 171     returns tuple of 4 vectors:
 172            - point location
 173            - normal (+ve is in the +ve y direction at iu=t=0) unit vector
 174            - along extent   (towrds +ve v)                    unit vector
 175            - along traverse (towards +ve w)                   unit vector
 176     '''
 177     p, normal, vec_v, vec_w = m.h.details(**m._vw2tiu_kw(v,w))
 178     if w > m.nt:
 179       normal = -normal
 180       vec_w  = -vec_w
 181     return p, normal, vec_v, vec_w