Smooth curve generator implemented in python

Given a list of ordered points as in the bottom of this image, my program produces the top part of the image:

For more about the theory behind it see my last blog post.
The run time for the program was on the order of a hundredth of a second but I was plotting the pixels one by one I think there are faster ways to plot the output.

The code;s not commented extensively but short and straightforward. if there are any questions let me know... it requires the Image module from PIL to make the image shown above (the only difference is I made the dots bigger in the image above with GIMP to make them stand out....

import Image
import math
import time

def plot(curve, blank):
    for point in curve:
        blank.putpixel(point, (255,255,255))
def bplot(curve, blank):
    for point in curve:
        blank.putpixel((point[0], point[1]+250), (255, 0, 0))
def circ(point1, point2, point3):
    ax = point1[0]
    ay = point1[1]
    ax = point1[0]
    ay = point1[1]
    bx = point2[0]
    by = point2[1]
    cx = point3[0]
    cy = point3[1]
    '''find the x,y and radius for the circle through the 3 points'''
    if (ax*by-ax*cy-cx*by+cy*bx-bx*ay+cx*ay) != 0:
        x=.5*(-pow(ay, 2)*cy+pow(ay, 2)*by-ay*pow(bx, 2)\
        -ay*pow(by, 2)+ay*pow(cy, 2)+ay*pow(cx, 2)-\
        pow(cx, 2)*by+pow(ax, 2)*by+pow(bx, 2)*\
        cy-pow(ax, 2)*cy-pow(cy, 2)*by+cy*pow(by, 2))\
        /(ax*by-ax*cy-cx*by+cy*bx-bx*ay+cx*ay)
        y=-.5*(-pow(ax, 2)*cx+pow(ax, 2)*bx-ax*pow(by, 2)\
        -ax*pow(bx, 2)+ax*pow(cx, 2)+ax*pow(cy, 2)-\
        pow(cy, 2)*bx+pow(ay, 2)*bx+pow(by, 2)*cx\
        -pow(ay, 2)*cx-pow(cx, 2)*bx+cx*pow(bx, 2))\
        /(ax*by-ax*cy-cx*by+cy*bx-bx*ay+cx*ay)
    else:
        return False
    
    r=pow(pow(x-ax, 2)+pow(y-ay, 2), .5)
    
    return x, y, r
def findpoint(eq1, eq2, point1, point2):
    '''find the centroid of the overlapping part of two circles
    from their equations'''
    thetabeg = math.acos((point1[0]-eq1[0])/eq1[2])
    thetaend = math.acos((point2[0]-eq1[0])/eq1[2])
    mid1x = eq1[2]*math.cos((thetabeg+thetaend)/2)+eq1[0]
    thetaybeg = math.asin((point1[1]-eq1[1])/eq1[2])
    thetayend = math.asin((point2[1]-eq1[1])/eq1[2])
    mid1y = eq1[2]*math.sin((thetaybeg+thetayend)/2)+eq1[1]

    thetabeg2 = math.acos((point1[0]-eq2[0])/eq2[2])
    thetaend2 = math.acos((point2[0]-eq2[0])/eq2[2])
    mid2x = eq2[2]*math.cos((thetabeg2+thetaend2)/2)+eq2[0]
    thetaybeg2 = math.asin((point1[1]-eq2[1])/eq2[2])
    thetayend2 = math.asin((point2[1]-eq2[1])/eq2[2])
    mid2y = eq2[2]*math.sin((thetaybeg2+thetayend2)/2)+eq2[1]
    return [(mid2x+mid1x)/2, (mid2y+mid1y)/2]
def main():
    starttime = time.time()
    ''' Ordered list of points, the first and last affect the shape
    of the curve but are not connected though drawn'''
    curve = [[23,150],[43,127],[100,49],[179,83],[253,155],[362,215],[445,122],[551,57],[560,23], [540, 10], [520,20],[500,50],[450,70]]
    blank = Image.new("RGB", (600, 600))
    bplot(curve, blank)
    
    for j in range(0, 6):
        newpoints = []
        for i in range(0, len(curve)-3):
            eq = circ(curve[i], curve[i+1], curve[i+2])
            eq2 = circ(curve[i+1], curve[i+2], curve[i+3])
            if eq == False or eq2 == False:
                newpoints.append([(curve[i+1][0]+curve[i+2][0])/2, (curve[i+1][1]+curve[i+2][1])/2])
            else:    
                newpoints.append(findpoint(eq, eq2, curve[i+1], curve[i+2]))
        for point in newpoints:
            point[0] = int(round(point[0]))
            point[1] = int(round(point[1]))
        for m in range(0, len(newpoints)):
            curve.insert(2*m+2, newpoints[m])
        
    plot(curve, blank)
    print time.time()-starttime
    blank.save("smoothcurve.png")
    
main()

**The above code runs on Windows 7 Home, Python 2.7.3, PIL 32 bit

I think it's nice to be able to do like the following, usually selection in image programs follows polygonal lines but with this algorithm you can do this instead:

And it's fast enough that you could make adjustments to add new nodes or move individual nodes and it recalculate in real time, which I think would be handy....
***LOL I know it looks like I just hand drew that outline around the cloud but I really gave it the coordinates of the red dots and it did all the rest
**Comment from Ofnuts from the gimp mailing list:
Ofnuts: The curve doesn't include the two end points and it seems it doesn't even get through one of the points (next before last at top right):

Me: Right I should have made it more clear, I left it as a comment in the code, but the first and last point figure into the calculation and do have some effect but are not directly connected to the rest of the curve. They could be left undrawn or automatically be selected as the top left of the bounding rectangle and the bottom right or as I was saying put anywhere to give the beginning and end of the curve the "feel" that it was coming from a certain direction before and after  the points that actually are all connected.

The second issue is that I was placing the red dots in gimp and reading off their coordinates also with gimp to put in the program, so it wasn't exact, but the coordinates actually put in the program are all plotted as pixels so in the end are drawn as part of the curve... I just checked and made sure 520,20 given in the program is where it should be, the red dots were wrong...