Slow but dense poisson disc distribution sampling

The idea for this algorithm is to somewhat randomly place many dots that are all at least a certain distance apart, without placing so many that it becomes hexagonal close packed packing which is the densest possible. It ends up looking a lot like how the light sensitive cells are arranged on the retina.
This algorithm is slow but has a good end result here it is for 200x200 no two closer than 5 units apart:

The algorithm makes a list of every single possible pixel value shuffled in random order F, then starts with a field with only one random point on it p of list P, pops the first value v off of F and adds it to P if the distance from v to any point p in P is greater than d, where d is the minimum distance between points. The algorithm stops when there are no longer any values in F. In the above at the end there were 2337 points in P.
**improvements?**
I could add logic so the points are on a toroidal surface with the top and bottom attached and the sides attached, also if all the points p in P were stored in an octree type structure I think there could be some good performance gains.

**Python Source Code**
import random
from random import shuffle
import Image
def fill(grid, sizegrid):
    for i in range(0, sizegrid[0]):
        for j in range(0, sizegrid[1]):
            grid.append([i,j])
    shuffle(grid)
def distance(pointa, pointb, mindistance):
    if abs(pointa[0]-pointb[0]) >= mindistance:
        return mindistance+1
    if abs(pointa[1]-pointb[1]) >= mindistance:
        return mindistance+1
    return ((pointa[0]-pointb[0])**2 + (pointa[1]-pointb[1])**2)**.5 
def plot(points, sizegrid):
    p = Image.new("RGB", sizegrid)
    for i in points:
        if i[0] >= 0 and i[1] >= 0:
            p.putpixel(i, (255,0,0))
    p.save("out.png")
def add(a, points,mindistance):
    for point in points:
        if distance(point, a, mindistance)<mindistance:
            return False
    points.append(a)
    return True
def main():
    sizegrid = (200, 200)
    points = [[random.randint(0, sizegrid[0]-1), random.randint(0, sizegrid[1]-1)]]
    fullgrid = []
    fill(fullgrid, sizegrid)
    mindistance = 5
    while(len(fullgrid) > 0):
          f = fullgrid.pop()
          if add(f,points, mindistance):
               points.append(f)
    plot(points, sizegrid)
    print(len(points))
main()