# ported to Python 26-APR-2003
#include <math.h>
#include <stdio.h>
#
#// Copyright YON - Jan C. Hardenbergh. Permission to copy is
#// granted as long as this copyright and this URL are maintained:
#// http://www.jch.com/~jch/NURBS/
#
#// Primary reference Les Peigl - On NURBS: A Survey published
#// in IEEE Computer Graphics and Applications (CG&A) January 1991.
#
#// http://www.cs.wpi.edu/~matt/courses/cs563/talks/nurbs.html
#
#// this is intended to explore NURBS - how they are evaluated.
#// It is explicitly unoptimized.
#
#// the NURBCurve used as an example here is a circle inscribed
#// in an equalateral triangle. I found it on page 374, 
#// Mathematical Elements for Computer Graphics, 
#// 2nd ed. David F. Rogers and J Alan Adams. McGraw Hill, 1990.

class Point:
    def __init__(self, value):
        self.data = [value[0], value[1], value[2]]
    def __getitem__(self, index):
        return self.data[index]
    def __add__(self, other):
        return  Point([self.data[0]+other[0], self.data[1]+other[1], self.data[2]+other[2]])
    def __mul__(self, other):
        return Point([self.data[0]*other, self.data[1]*other, self.data[2]*other])
    def Print(self):
        print "Point : %g, %g, %g" % (self.data[0], self.data[1], self.data[2])

def B(i,k,t,knots):
    ret = 0
    if k>0:
        n1 = (t-knots[i])*B(i,k-1,t,knots)
        d1 = knots[i+k] - knots[i]
        n2 = (knots[i+k+1] - t) * B(i+1,k-1,t,knots)
        d2 = knots[i+k+1] - knots[i+1]
        if d1 > 0.0001 or d1 < -0.0001:
            a = n1 / d1
        else:
            a = 0
	if d2 > 0.0001 or d2 < -0.0001:
            b = n2 / d2
        else:
            b = 0
        ret = a + b
        #print "B i = %d, k = %d, ret = %g, a = %g, b = %g\n"%(i,k,ret,a,b)
    else:
        if knots[i] <= t and t <= knots[i+1]:
            ret = 1
        else:
            ret = 0
    return ret

def C(t, order, points, weights, knots):
    c = Point([0,0,0])
    rational = 0
    i = 0
    while i < len(points):
	b = B(i, order, t, knots)
        p = points[i] * (b * weights[i])
        c = c + p
        rational = rational + b*weights[i]
        i = i + 1
    return c * (1.0/rational)

# main ***************************************************************
######################################################################
rad3 = 1.732

def main():
    controlPoints = [ Point([0,0,0]), Point([1,0,0.1]), Point([.5,.5*rad3,0.2]),
                      Point([0,rad3,0.3]), Point([-.5,.5*rad3,0.4]),
                      Point([-1,0,0.5]), Point([0,0,0.6])]
    knots = [0,0,0,1,1,2,2,3,3,3]
    weights = [1,.5,1,.5,1,.5,1]
    order = 2
    steps = 15
    i = 0
    while i < 15:
        t = i*(3.0/15)
        c = C(t, order, controlPoints, weights, knots)
        c.Print()
        i = i + 1
        
if __name__ == '__main__':
    main()