This post shows, how to get some NURBS editing capabilities in Python by making use of a mature and sophisticated, yet open source NURBS library called SISL.
Basically, the SISL library lets you do everything you would expect from a CAD system. I was especially interested in:
- find intersection of curves/surfaces
- find closest points on curve/surface
- calculate curvature and derivatives
The following image shows a B-Spline curve with the control points and curvature. Three closest points to the curve are evaluated.
 |
| B-Spline with curvature and closest points on curve |
Roughtly the necessary steps can be described as:
- compile SISL library
- since SISL compiles into a .lib (windows file associations are assumed here) static libary, write wrapper to pack needed functions into a .dll dynamic library
- load .dll into Python, access it with help of cTypes
- write wrapper/class in Python to easily access SISL functions
Get and compile SISL library
The library can be downloaded from here. There is available the so called GoTools, which is a huge software package, but we only need SISL. As a result you get the sisl.lib file.Pack needed functions into .dll dynamic library
Since in Python, as far as I know, the standard way to access libraries is via a .dll dynamic library, we need to write one that exposes the needed functions.This minimal example creates a SISL curve and calculates the closest points to the curve.
A more general description about how to compile a .dll dynamic library to access from Python is given by the Scipy Cookbook. The following code is mostly taken from the Cookbook and adapted accordingly. For more information about how to call functions from the SISL library, see the SISL manual.
SISLlib.cpp:
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <string>
#include <stdexcept>
#include <cstdlib>
#include "sisl.h"
using namespace std;
#ifdef SISLlib_DLL
#ifdef SISLlib_EXPORTS
#define SISLlib_API __declspec(dllexport)
#else
#define SISLlib_API __declspec(dllimport)
#endif /* SISLlib_EXPORTS */
#else
#define SISLlib_API extern /* XXX confirm this */
#endif /* SISLlib_DLL */
#ifdef __cplusplus
extern "C" {
#endif
extern SISLlib_API int SISLclosestPointToCurve(int kind, double *knots1, double *coef1, int number1, int order1, int dim1, double *PointsIn, int nPointsIn, double *pointparout) {
/* BSpline Kurve erzeugen */
SISLCurve* curve1 = newCurve(number1, // number of control points
order1, // order of spline curve (degree + 1)
knots1, // pointer to knot vector (parametrization)
coef1, // pointer to coefficient vector (control points)
kind, // kind, 1 = polynomial B-spline curve, 2 = Rational B-spline (nurbs) curve
dim1, // dimension
1); // 0=no copying of information, 'borrow' arrays
if (!curve1) {
throw runtime_error("Error occured while generating curve1.");
}
/* search closest points to surface */
double epsco = 1.0e-8; /*computational resolution*/
double epsge = 1.0e-8;
double point[2];
int numintpt;
double *intpar;
int numintcu=0;
SISLIntcurve **intcurve=NULL;
int jstat;
for (int i = 0; i < nPointsIn; i++){
point[0]=PointsIn[0+i*2];
point[1]=PointsIn[1+i*2];
s1953(curve1, point, dim1, epsco, epsge, &numintpt, &intpar, &numintcu, &intcurve, &jstat);
if (jstat < 0){
throw runtime_error("Error occured inside call to SISL routine s1953.");
} else if (jstat > 0){
cerr << "WARNING: warning occured inside call to SISL routine s1953.\n" << endl;
}
if(numintpt >= 1){
*(pointparout+i) = intpar[0];
}
}
// cleaning up
freeCurve(curve1);
freeIntcrvlist(intcurve,numintcu);
return 0;
}
To compile, use e.g. this makefile, adjust the paths (libpath) accordingly:
CXX = cl.exe
LINK = link.exe
CPPFLAGS = -D_WIN32 -D_USRDLL -DSISLlib_DLL -DSISLlib_EXPORTS -Iinclude -nologo -GS -W3 -Wp64 -MD -O2 -EHsc
CXXFLAGSALL = -EHsc -nologo -GS -W3 -Wp64 $(CPPFLAGS)
CXXFLAGSDBG = -MDd -Od -Z7 -RTCcsu
CXXFLAGSOPT = -MD -O2
#CXXFLAGS = $(CXXFLAGSALL) $(CXXFLAGSDBG)
CXXFLAGS = $(CXXFLAGSALL) $(CXXFLAGSOPT)
LINKFLAGSALL = /nologo /DLL /libpath:lib sisl.lib
LINKFLAGSDBG = /DEBUG
LINKFLAGSOPT =
#LINKFLAGS = $(LINKFLAGSALL) $(LINKFLAGSDBG)
LINKFLAGS = $(LINKFLAGSALL) $(LINKFLAGSOPT)
all: SISLlib.dll
SISLlib.dll: SISLlib.obj
$(LINK) $(LINKFLAGS) SISLlib.obj /OUT:SISLlib.dll
svm.obj: svm.cpp svm.h
$(CXX) $(CXXFLAGS) -c SISLlib.cpp
clean:
-erase /Q *.obj *.dll *.exp *.lib
With Visual Studio this compiles to SISLlib.dll with the nmake command from the Visual Studio Command Prompt. The SISLlib.dll can then be loaded into Python.
Load SISLlib.dll into Python
When the .dll compiles correctly the hardest work is already done. To load the library into Python following command can be used: import numpy as np
SISLlib_dll = np.ctypeslib.load_library("SISLlib",r"pathToLibrary\SISLlib")
Python Class to easily access SISL functions
In Python it is convenient to write a Curve class, which is the interface to the SISL functions.When a curve object is initialized, it takes as arguments the knot vector, the control points, information about the order of the curve and possibly weights of the coefficients (needed if curve should be a NURBS). It must declare, which data types are expected by a function in the SISLlib.dll as input and what is the ouput data type. It also should provide a Pythonic way to call functions to e.g. calculate the closest point to the curve.
# create SISL Curve Object
curve = SISLcurve(knots,coefs,order,weights)
# get closest points to curve
parOut = curve.SISLclosestPointToCurve(points)
This is a proposal of how a curve class could look like.
pySISL.py:
from __future__ import division
import ctypes
import numpy as np
import matplotlib.pyplot as plt
class SISLcurve():
def __init__(self, knots, coefs, order, weights=[]):
if weights == []:
self.kind = 1 # polynomial B-Spline
self.weights = weights
else:
self.kind = 2 # rational B-Spline (NURBS)
self.weights = weights.astype("double")
self.knots = knots.astype("double")
self.coefs = coefs.astype("double")
self.order = order
self.dim = coefs.shape[0]
# import SISL libaray
self.SISLlib_dll = np.ctypeslib.load_library("SISLlib",r"G:\ProgrammierProjekte\SISL\SISLlib")
# declare SISLclosestPointToCurve function
self.SISLlib_dll.SISLclosestPointToCurve.restype = ctypes.c_int
self.SISLlib_dll.SISLclosestPointToCurve.argtypes = [ctypes.c_int, #kind 1: polynomial B-Spline, 2: rational B-Spline (NURBS)
ctypes.POINTER(ctypes.c_double), #knots
ctypes.POINTER(ctypes.c_double), #coefs
ctypes.c_int, ctypes.c_int, ctypes.c_int, #nCoefs, order, dim
ctypes.POINTER(ctypes.c_double), #points
ctypes.c_int, #nPoints
ctypes.POINTER(ctypes.c_double)] #pointsParOut
def prepareSISLdata(self):
kind = self.kind
knots = self.knots
number = self.coefs.shape[1]
coefs = self.coefs.reshape((self.coefs.size),order="F")
order = self.order
dim = self.dim
# if curve is rational, coefs need to be multiplied by weight
if self.weights != []:
weights = self.weights
tmp_w = np.tile(weights,(dim,1))
tmp_c = self.coefs * tmp_w
tmp_cw = np.vstack((tmp_c,weights))
coefs = tmp_cw.reshape((tmp_cw.size),order="F")
return kind, knots, number, coefs, order, dim
def SISLclosestPointToCurve(self, points):
kind, knots, number, coefs, order, dim = self.prepareSISLdata()
nPoints = points.shape[1]
points = points.reshape((points.size),order="F")
points = points.astype("double")
pointsParOut = np.zeros(nPoints)
self.SISLlib_dll.SISLclosestPointToCurve(kind,
knots.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
coefs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
number, order, dim,
points.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
nPoints,
pointsParOut.ctypes.data_as(ctypes.POINTER(ctypes.c_double)))
return pointsParOut
if __name__ == "__main__":
# Data for 2D Spline
knots = np.array([0,0,0,1,1,2,2,3,3,4,4,4])
coefs_x = np.array([0,1,1,1,0,-1,-1,-1,0])
coefs_y = np.array([-1,-1,0,1,1,1,0,-1,-1])
coefs_z = np.array([0,0,0,0,0,0,0,0,0])
coefs = np.vstack((coefs_x,coefs_y,coefs_z))
weights = np.array([1,(2**0.5)/2,1,(2**0.5)/2,1,(2**0.5)/2,1,(2**0.5)/2,1])
order = 3
# Points
points_x = np.array([0.950000000000000,0.948365660353047,0.946681081750349,0.944627123146684,0.942188773156426,0.939399032892644,0.936232939571863,0.932690108317814,0.928770254558779,0.924462521735862,0.919770399050986,0.914760631089697,0.909376263946973,0.903585812191018,0.897468630037332,0.891015499686028,0.884148637040933,0.876982392633972,0.869461355406687,0.861554859384639,0.853377314168268,0.844768886370890,0.835927093263629,0.826659503750130,0.817152570478120,0.807224602039578,0.797065796079186,0.786515110030005,0.775693455728189,0.764537838907564,0.753050219497601,0.741289980627127,0.729171580374359,0.716773946925715,0.704059298956755,0.690978248179852,0.677612415560273,0.663917247813943,0.649864255028625,0.635503203034109,0.620821802447561,0.605786959796453,0.590392566215039,0.574694342989329,0.558630565252025,0.542185206713670,0.525373448459848,0.508204975104239,0.490632813730356,0.472635874021836,0.454236471406481,0.435439193363185,0.416218675261615,0.396570222858047,0.376542216747913,0.356146885282818,0.335382062721527,0.314254523947087,0.292855861004967,0.271191942243154,0.249289199749707,0.227213115084334,0.205041874526635,0.182805870132726,0.160565960001584,0.138438450747921,0.116457821240756,0.0947018141041910,0.0732783327465884,0.0522156119082409,0.0316246164284243,0.0115659850707168,-0.00791202104385392,-0.0267042000226111,-0.0448139470270519,-0.0621263608066405,-0.0786633164090686,-0.0943337919206955,-0.109168987800449,-0.123103860278432,-0.136168875212454,-0.148344178963388,-0.159644690042500,-0.170098972799865,-0.179716486713481,-0.188530723275019,-0.196584499430869,-0.203908474072773,-0.210535813047408,-0.216518891589666,-0.221902400415273,-0.226721422700915,-0.231011507186882,-0.234811457180879,-0.238156960106410,-0.241084055321751,-0.243621812078548,-0.245789804776811,-0.247616108931611,-0.249497141016798])
points_y = np.array([-0.0994321452870676,-0.0994226258843406,-0.0994011421318266,-0.0993638635104812,-0.0993061039921876,-0.0992191204394156,-0.0991094962412971,-0.0989758995773281,-0.0988172225349698,-0.0986324951964819,-0.0984017275113659,-0.0981271747430788,-0.0978224338994053,-0.0974815622227530,-0.0970794433470467,-0.0966373768472060,-0.0961540453186260,-0.0956029034674430,-0.0950148959841584,-0.0943597983531194,-0.0936499939404355,-0.0928880723801365,-0.0920550818411189,-0.0911715420258190,-0.0902160232481631,-0.0892055634568082,-0.0881255694188480,-0.0869769598211418,-0.0857658423138637,-0.0844835330960565,-0.0831380138003898,-0.0817249600103063,-0.0802376046712566,-0.0786796175308041,-0.0770530700400894,-0.0753426682406730,-0.0735450754025746,-0.0716739843232092,-0.0697241104204034,-0.0676828864957862,-0.0655604886165625,-0.0633668842861504,-0.0610870449049815,-0.0587123402033290,-0.0562669315097688,-0.0537535218622797,-0.0511523935680138,-0.0484660107488078,-0.0457142971790514,-0.0429012017883088,-0.0400093517599968,-0.0370374806273727,-0.0340077065131031,-0.0309256570274795,-0.0277765350027134,-0.0245578071914544,-0.0212958999960175,-0.0179995362192438,-0.0146458041074342,-0.0112460425963343,-0.00782424211239759,-0.00437524956478380,-0.000893563813683097,0.00260435993597310,0.00610447748869498,0.00962021509300744,0.0131358483933153,0.0166339093545291,0.0201111949029309,0.0235539076034093,0.0269612433102411,0.0303173926105221,0.0336211753342074,0.0368681605591590,0.0400406671704278,0.0431536061729911,0.0461705095054608,0.0491158099116442,0.0519509427343211,0.0547034649982273,0.0573375222960913,0.0598641536018524,0.0622741490265492,0.0645495916327474,0.0667100720625331,0.0687237608042033,0.0705943638721636,0.0723321403710434,0.0739092963536412,0.0753254443305403,0.0765908300800888,0.0776974052822079,0.0786305979877449,0.0793888141372688,0.0799705067334637,0.0803762588498084,0.0806083507407432,0.0806585219718005,0.0805152222178871,0.0800000000000000])
points = np.vstack((points_x,points_y))
# create SISL Curve Object
curve = SISLcurve(knots,coefs,order,weights)
# get closest points to curve
parOut = curve.SISLclosestPointToCurve(points)
Now to actually see whats going on we need to calculate the points in euclidean space from the parameters in parOut. This can be done using SISL function s1542 (evaluate points on curve by parameter) but that would already go beyond this rough guide. The result however could look something like that:
 |
| Circle described as NURBS curve, red points on the curve are the closest points on the circle to the green points in the middle |
Keine Kommentare:
Kommentar veröffentlichen