from __future__ import division from scipy.interpolate import splev from numpy import (sqrt, linspace, sin, cos, array, hstack, transpose, pi, ones_like, zeros_like, where, r_) import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from matplotlib.tri import Triangulation def splevn(t, tks): """Like splev but control points are vectors.""" return zip(*[splev(t, [tks[0], p, tks[-1]]) for p in zip(*tks[1])]) # Knots, weighted control points, degree: arcNURBS = [[0,0,0,0.25,0.25,0.5,0.5,0.75,0.75,1,1,1], [[1,0,1], [sqrt(0.5), sqrt(0.5), sqrt(0.5)], [0,1,1], [-sqrt(0.5), sqrt(0.5), sqrt(0.5)], [-1,0,1], [-sqrt(0.5), -sqrt(0.5), sqrt(0.5)], [0,-1,1], [sqrt(0.5), -sqrt(0.5), sqrt(0.5)], [1,0,1]], 2] theta = linspace(0, 2*pi, 7) pts = transpose([cos(theta), sin(theta)]) w = 1.0 / array([1.,2,1,2,1,2,1]) arcNURBS = [array([0.,0,0,1,1,2,2,3,3,3]) / 3.0, hstack([pts, w[:,None]]), 2] Pw = array(arcNURBS[1]) t = linspace(arcNURBS[0][0], arcNURBS[0][-1], 6*15+1) x, y, w = transpose(splevn(t, arcNURBS)) fig = plt.figure(figsize=(10, 10)) ax = fig.add_subplot(111, projection='3d') ax.view_init(azim=-50, elev=20) for setRange in (ax.set_xlim, ax.set_ylim): setRange(-1.5, 1.5) ax.set_zlim(-0.5, 1.5) # Shade in the cone: xs = hstack([[0], x/w]) ys = hstack([[0], y/w]) ws = hstack([[0], ones_like(w)]) tris = [] for i in range(len(x)): tris.append([0, (i + 1) % len(xs), (i + 2) % len(xs)]) color = (0.9,)*3 cone = ax.plot_trisurf(xs, ys, tris, ws, color=color, edgecolor=color, alpha=1.0, lw=0) Pw = transpose(arcNURBS[1]) ax.plot(Pw[0]/Pw[2], Pw[1]/Pw[2], Pw[2]/Pw[2], 'k.-', label='NURBS control w/o weights') p = ax.plot(Pw[0], Pw[1], Pw[2], '.-', color='b', label='3D control polygon')[0] p.set_zorder(-1000) ax.plot(x, y, w, 'b', label='3D parabolas') ax.plot(x/w, y/w, ones_like(w), color=(0.85,0.1,0.1), linewidth=2, label='NURBS circle at w = 1') if True: ax.plot(Pw[0], Pw[1], ones_like(Pw[2]), '.-', color='b', alpha=0.2) ax.plot(x, y, ones_like(w), 'b', alpha=0.2) if not True: plots = [] plots += ax.plot(Pw[0]/Pw[2], Pw[1]/Pw[2], 0*Pw[2]/Pw[2], 'k.-', alpha=0.2) plots += ax.plot(Pw[0], Pw[1], 0*Pw[2], '.-', color='b', alpha=0.2) plots += ax.plot(x, y, 0*w, 'b', alpha=0.2) plots += ax.plot(x/w, y/w, zeros_like(w), color=(0.85,0.1,0.1), linewidth=2, alpha=0.2) for i, p in enumerate(plots): p.set_zorder(-1000-i) DRAW_RULES_ON_CONE = False if DRAW_RULES_ON_CONE: for i in r_[:len(t):10]: ax.plot([0, x[i], x[i]/w[i]], [0, y[i], y[i]/w[i]], [0, w[i], 1], 'k.--') for i in where(Pw[-1] != 1)[0]: ax.plot([0, Pw[0,i]/Pw[2,i]], [0, Pw[1,i]/Pw[2,i]], [0, 1], 'k.--', alpha=0.5) plt.xlabel(r"$x$", fontsize=18) plt.ylabel(r"$y$", fontsize=18) ax.set_zlabel(r"$z$", fontsize=18) plt.legend() plt.savefig("NURBS-circle-3D.svg", transparent=True, bbox_inches="tight") plt.show()