def N(self,U,i,p,u): if p == 0: if U[i] <= u and u < U[i+1]: return 1 else: return 0 else: v1 = 0 n1 = u - U[i] d1 = U[i+p] - U[i] if d1 != 0: v1 = 1.*n1/d1*self.N(U,i,p-1,u) v2 = 0 n2 = U[i+p+1] - u d2 = U[i+p+1] - U[i+1] if d2 != 0: v2 = 1.*n2/d2*self.N(U,i+1,p-1,u) return v1 + v2 def V(self,u,v,index): num = 0.0 den = 0.0 for i in xrange(0,self.m): for j in xrange(0,self.n): num += self.N(self.uknot,i,self.p,u)*self.N(self.vknot,j,self.q,v)*self.P[i][j][3]*self.P[i][j][index] den += self.N(self.uknot,i,self.p,u)*self.N(self.vknot,j,self.q,v)*self.P[i][j][3] if den == 0: return 0 else: return num/den