2,2,2 smooth curve source code

The main difference between this and Bezier curves or Nurbs is this doesn't use control points to guide the curve, it just goes through the points given directly. The way it does it is given 4 points it finds the best quadratic fit through the first 3 points parameterized from t=0..2 and the best quadractic fit for the last 3 points taken in reverse order parameterized from t=0..2 and makes a new point at the average of both curves at t=1.5. Then it does that for every group of 4 points to get twice the number of points, and recurses until it reaches the required resolution... 
The source code below takes for example these points:
points = [[50,300],[190,210],[344,256],[446,384],[580,294],[510,115]]
and generates a smooth curve through them in constant time with respect to the final number of points in the curve as described in the blog post before this...
here's a plot of the points, y increases as you move down because in image programs that's how they arrange the axis for some reason.

After the first iteration we find points between those given:

second iteration:

third:

After 10 iterations:

**Runtime complexity**
The program runs in constant time vs the final number of points generated.
**Source Code(Python) This source code requires PIL image library for python 2.7**
import Image
def e(coeff, t):
    return coeff[2]*t**2 + coeff[1]*t+coeff[0]

def m(a,b,c,d):
    coefflx = [b[0], -.5*a[0] + .5*c[0], .5*c[0]+.5*a[0]-b[0]]
    coeffly =  [b[1], -.5*a[1] + .5*c[1], .5*c[1]+.5*a[1]-b[1]]
    a,b,c
    d,c,b
    coeffrx = [c[0], -.5*d[0] + .5*b[0], .5*b[0]+.5*d[0]-c[0]]
    coeffry =  [c[1], -.5*d[1] + .5*b[1], .5*b[1]+.5*d[1]-c[1]]
    Xl = [e(coefflx, .5), e(coeffly, .5)]
    Xr = [e(coeffrx, .5), e(coeffry, .5)]
    return [.5*(Xl[0]+Xr[0]), .5*(Xl[1]+Xr[1])]
def n(a,b,c, side):
    coeffx = [b[0], -.5*a[0] + .5*c[0], .5*c[0]+.5*a[0]-b[0]]
    coeffy =  [b[1], -.5*a[1] + .5*c[1], .5*c[1]+.5*a[1]-b[1]]
    if side ==0:
        return [e(coeffx, -.5), e(coeffy, -.5)]
    else:
        return [e(coeffx, .5), e(coeffy, .5)]
def twotwotwo(plot, points, iter):
    points2 = [points[0], n(points[0], points[1], points[2], 0), points[1]]
    for i in range(0, len(points)-3):
        points2.append(m(points[i], points[i+1], points[i+2], points[i+3]))
        points2.append(points[i+2])
    points2.append(n(points[len(points)-3], points[len(points)-2], points[len(points)-1], 1))
    points2.append(points[len(points)-1])
    for i in points2:
        plot.putpixel((int(i[0]), int(i[1])), (255,255,255))
    plot.save("plot"+str(iter)+".png")
    return points2
def main():
    plot = Image.new("RGB", (640,640))
    points = [[50,300],[190,210],[344,256],[446,384],[580,294],[510,115]]
    for iter in range(1, 5):
        points = twotwotwo(plot, points, iter)
    
main()