PythonCAD/Generic/nurbs.py

   1 #
   2 # Copyright (c) 2003 Art Haas
   3 #
   4 # This file is part of PythonCAD.
   5 #
   6 # PythonCAD is free software; you can redistribute it and/or modify
   7 # it under the terms of the GNU General Public License as published by
   8 # the Free Software Foundation; either version 2 of the License, or
   9 # (at your option) any later version.
  10 #
  11 # PythonCAD is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15 #
  16 # You should have received a copy of the GNU General Public License
  17 # along with PythonCAD; if not, write to the Free Software
  18 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  19 #
  20 # NURBS curves
  21
  22
  23 import array
  24
  25 class Nurb(object):
  26     def __init__(self, ctrlpts, knots, order):
  27         if not isinstance(ctrlpts, list):
  28             raise TypeError, "Invalid control point list: " + str(ctrlpts)
  29         if not isinstance(knots, list):
  30             raise TypeError, "Invalid knot list: " + str(knots)
  31         if not isinstance(order, int):
  32             raise TypeError, "Invalid order; " + str(order)
  33         if order < 2 or order > 16: # what is a good max value?
  34             raise ValueError, "Invalid order: %d" % order
  35         _ctrlpts = []
  36         for _pt in ctrlpts:
  37             if not isinstance(_pt, tuple):
  38                 raise TypeError, "Invalid control point: " + str(_pt)
  39             _len = len(_pt)
  40             if not (1 < _len < 4):
  41                 raise ValueError, "Invalid tuple length: " + str(_pt)
  42             if _len == 2:
  43                 _x, _y = _pt
  44                 _w = 1.0
  45             else:
  46                 _x, _y, _w = _pt
  47             if not isinstance(_x, float):
  48                 _x = float(_x)
  49             if not isinstance(_y, float):
  50                 _y = float(_y)
  51             if not isinstance(_w, float):
  52                 _w = float(_w)
  53             if not (_w > 0.0):
  54                 raise ValueError, "Invalid weight: %g" % _w
  55             _ctrlpts.append((_x, _y, _w))
  56         _knots = []
  57         for _knot in knots:
  58             if not isinstance(_knot, float):
  59                 _knot = float(_knot)
  60             if (_knot < 0.0 or _knot > 1.0):
  61                 raise ValueError, "Invalid knot value: %g" % _knot
  62             for _val in _knots:
  63                 if _knot < (_val - 1e-10):
  64                     raise (ValueError, "Invalid decreasing knot: %g < %g" %
  65                            (_knot, _val))
  66             _knots.append(_knot)
  67         print "knots: " + str(_knots)
  68         print "ctrl: " + str(_ctrlpts)
  69         print "order: %d" % order
  70         _clen = len(_ctrlpts)
  71         if _clen < order:
  72             raise ValueError, "Order greater than number of control points."
  73         if len(_knots) != (_clen + order):
  74             raise ValueError, "Knot/Control Point/Order number error."
  75         self.__ctrlpts = _ctrlpts
  76         self.__knots = _knots
  77         self.__order = order
  78
  79     def getControlPoints(self):
  80         return self.__ctrlpts[:]
  81
  82     def getKnots(self):
  83         return self.__knots[:]
  84
  85     def getOrder(self):
  86         return self.__order
  87
  88     def calculate(self, count):
  89         if not isinstance(count, int):
  90             raise TypeError, "Invalid count: " + str(count)
  91         _cpts = self.__ctrlpts
  92         _knots = self.__knots
  93         _dt = 1.0/float(count)
  94         _p = self.__order - 1
  95         _pts = []
  96         for _c in range(count):
  97             _t = _c * _dt
  98             # print "time: %g" % _t
  99             _nx = _ny = _nw = 0.0
 100             for _i in range(len(_cpts)):
 101                 # print "using cpt %d" % _i
 102                 _x, _y, _w = _cpts[_i]
 103                 _Ni = self._N(_i, _p, _t)
 104                 _nx = _nx + (_Ni * _x)
 105                 _ny = _ny + (_Ni * _y)
 106                 _nw = _nw + (_Ni * _w)
 107             # print "nw: %.3f" % _nw
 108             # print "nx: %.3f" % _nx
 109             # print "ny: %.3f" % _ny
 110             if abs(_nw) > 1e-10:
 111                 _pts.append((_nx/_nw, _ny/_nw))
 112             else:
 113                 print "zero weight: %f, %f" % (_nx, _ny)
 114
 115         return _pts
 116
 117     def _N(self, i, p, t):
 118         # print "_N() ..."
 119         _flag = False
 120         if abs(t - 1.0) < 1e-10 and False:
 121             _flag = True
 122         if _flag:
 123             print "i: %d" % i
 124             print "p: %d" % p
 125             print "t: %.3f" % t
 126         _knots = self.__knots
 127         _ki = _knots[i]
 128         _kin = _knots[i + 1]
 129         if _flag:
 130             print "ki: %.3f" % _ki
 131             print "kin: %.3f" % _kin
 132         if p == 0:
 133             if ((_ki - 1e-10) < t < _kin):
 134                 _val = 1.0
 135             else:
 136                 _val = 0.0
 137         else:
 138             _kip = _knots[i + p]
 139             _kipn = _knots[i + p + 1]
 140             if _flag:
 141                 print "kip: %.3f" % _kip
 142                 print "kipn: %.3f" % _kipn
 143             _t1 = 0.0
 144             _v1 = _kip - _ki
 145             if abs(_v1) > 1e-10:
 146                 _v2 = t - _ki
 147                 if abs(_v2) > 1e-10:
 148                     _t1 = (_v2/_v1) * self._N(i, (p - 1), t)
 149             _t2 = 0.0
 150             _v1 = _kipn - _kin
 151             if abs(_v1) > 1e-10:
 152                 _v2 = _kipn - t
 153                 if abs(_v2) > 1e-10:
 154                     _t2 = (_v2/_v1) * self._N((i + 1), (p - 1), t)
 155             _val = _t1 + _t2
 156         if _flag:
 157             print "val: %f" % _val
 158         return _val
 159
 160     def writedata(self, count, fname):
 161         if not isinstance(count, int):
 162             raise TypeError, "Invalid count: " + str(count)
 163         _f = file(fname, "w")
 164         for _pt in self.calculate(count):
 165             _x, _y = _pt
 166             _f.write("%f %f\n" % (_x, _y))
 167         _f.close()
 168         _f = file('control_points', "w")
 169         for _pt in self.getControlPoints():
 170             _x, _y, _w = _pt # ignore weight
 171             _f.write("%f %f\n" % (_x, _y))
 172         _f.close()
 173         _f = file('knots', "w")
 174         for _knot in self.getKnots():
 175             _f.write("%f 0.0\n" % _knot)
 176         _f.close()
 177
 178     def _NN(self, i, p, t):
 179         _cpts = self.__ctrlpts
 180         _cl = len(_cpts)
 181         _knots = self.__knots
 182         _kl = len(_knots) - 1
 183         _val = _kl * [0.0]
 184         #
 185         # calculate values for 0
 186         #
 187         for _i in range(_kl):
 188             if ((_knots[i] - 1e-10) < t < _knots[i + 1]):
 189                 _val[_i] = 1.0
 190         #
 191         # calculate values up to the degree
 192         #
 193         for _j in range(1, (p + 1)):
 194             for _i in range(_kl - _j):
 195                 _ki = _knots[_i]
 196                 _kin = _knots[_i + 1]
 197                 _kip = _knots[_i + p]
 198                 _kipn = _knots[_i + p + 1]
 199                 _t1 = 0.0
 200                 _n = _val[_i]
 201                 if abs(_n) > 1e-10:
 202                     _v1 = _kip - _ki
 203                     if abs(_v1) > 1e-10:
 204                         _v2 = t - _ki
 205                         if abs(_v2) > 1e-10:
 206                             _t1 = (_v2/_v1) * _n
 207                 _t2 = 0.0
 208                 _n = _val[_i + 1]
 209                 if abs(_n) > 1e-10:
 210                     _v1 = _kipn - _kin
 211                     if abs(_v1) > 1e-10:
 212                         _v2 = _kipn - t
 213                         if abs(_v2) > 1e-10:
 214                             _t2 = (_v2/_v1) * _n
 215                 _val[_i] = _t1 + _t2