Bézier curve fitting

bezier.py

import numpy as np
import matplotlib.pyplot as plt
import random

def Bx(t, px, n):
    if n == 1:
        return px[0]
    else:
        return (1-t)*Bx(t, px[:n-1],n-1) + t*Bx(t, px[1:],n-1)
   
def By(t, py, n):
    if n == 1:
        return py[0]
    else:
        return (1-t)*By(t, py[:n-1],n-1) + t*By(t, py[1:],n-1)

t = np.arange(0.0,1.0, 0.005)
px = []
py = []

n = int(raw_input())
for i in range(n):
    px.append(random.randint(-20,20))
    py.append(random.randint(-20,20))


px.sort()
for i in range(n):
    plt.annotate(str(i+1), xy=(px[i],py[i]))
for i in range(0, (n-1)/3) :
    plt.plot(Bx(t,px[i*3:(i+1)*3+1],4),By(t,py[i*3:(i+1)*3+1],4), 'b--')

if (n-1)% 3 != 0:
    p = (n-1)%3 + 1
    plt.plot(Bx(t, px[n-p:n],p), By(t,py[n-p:n],p), 'b--')

plt.plot(px,py,'ro')
plt.axis([-25,25,-25,25])
plt.show()