# # Une surface NURBS # Vincent Legat - 2018 # Ecole Polytechnique de Louvain # from numpy import * import matplotlib import matplotlib.pyplot as plt import mpl_toolkits.mplot3d # ========================================================================= def b(t,T,i,p): if p == 0: return (T[i] <= t)*(t < T[i+1]) else: u = 0.0 if T[i+p ] == T[i] else (t-T[i])/(T[i+p]- T[i]) * b(t,T,i,p-1) u += 0.0 if T[i+p+1] == T[i+1] else (T[i+p+1]-t)/(T[i+p+1]-T[i+1]) * b(t,T,i+1,p-1) return u # ============================= mainProgram =============================== X = [[1,2,3],[1,2,3],[1,2,3]] Y = [[1,1,1],[2,2,2],[3,3,3]] Z = [[1,2,3],[2,3,1],[3,1,4]] T = [0.1,0,0,1,1,1.1] S = [0.1,0,0,1,1,1.1] p = 2; n = 5 t = arange(T[p],T[n-p]+0.05,0.1) s = arange(S[p],S[n-p]+0.05,0.1) Bt = zeros((n-p,len(t))) Bs = zeros((n-p,len(s))) for i in range(0,n-p): Bt[i,:] = b(t,T,i,p) Bs[i,:] = b(s,S,i,p) x = Bs.T @ X @ Bt y = Bs.T @ Y @ Bt z = Bs.T @ Z @ Bt matplotlib.rcParams['toolbar'] = 'None' myColorMap = matplotlib.cm.jet plt.figure("Surface NURBS") ax = plt.axes(projection='3d',proj_type='ortho',azim=235) ax.plot_surface(x,y,z,alpha=1,cmap=myColorMap,linewidth=0.5,edgecolors='k') ax.set_xticks([1,1.5,2,2.5,3]) ax.set_yticks([1,1.5,2,2.5,3]) ax.set_zticks([1,1.5,2,2.5,3,3.5,4]) ax.margins(x=0,y=0,z=0) plt.tight_layout() plt.show() # =========================================================================