I am trying to classify an image by selecting a pixel at random, then finding all pixels in the image that are a certain euclidian distance in color space from that original pixel. My current script takes a prohibitively long time. I wonder if I am able to use this equation to generate a boolean matrix that will allow quicker manipulation of the image.
( x-cx ) ^2 + (y-cy) ^2 + (z-cz) ^ 2 < r^2 
Here is the code I am using now:
import PIL, glob, numpy, random, math, time
def zone_map(picture, threshold):
    im = PIL.Image.open(picture)
    pix = im.load()
    [width, height] = im.size
    mask = numpy.zeros((width,height))
    while 0 in mask:
        x = random.randint(0, width)
        y = random.randint(0, height)
        if mask[x, y] == 0:
            point = pix[x,y]
            to_average = {(x, y): pix[x, y]}
            start = time.clock()
            for row in range(0, width):
                for column in range(0, height):
                    if euclid_dist(point, pix[row,column]) <= threshold:
                        to_average[(row,column)] = pix[row, column]
            #to_average = in_sphere(pix, point)
            end = time.clock()
            print(end - start)
            to_average_sum = (0, 0, 0)
            for value in to_average.values():
                to_average_sum = tuple_sum(to_average_sum, value)
            average = tuple_divide(to_average_sum, len(to_average.values()))
            for coordinate in to_average.keys():
                pix[coordinate] = average
                mask[coordinate] = 1
                unique, counts = numpy.unique(mask, return_counts=True)
                progress = dict(zip(unique, counts))
                print((progress[1] / progress[0])*100, '%')
    im.save()
    return im
def euclid_dist(tuple1, tuple2):
    """
    Finds euclidian distance between two points in n dimensional sapce
    """
    tot_sq = 0
    for num1, num2 in zip(tuple1, tuple2):
        tot_sq += (num1 + num2)**2
    return math.sqrt(tot_sq)
def tuple_sum(tuple1, tuple2):
    """
    Returns tuple comprised of sums of input tuples 
    """
    sums = []
    for num1, num2 in zip(tuple1, tuple2):
        sums.append(num1 + num2)
    return tuple(sums)
def tuple_divide(tuple1, divisor):
    """
    Divides numerical values of tuples by divisisor, yielding integer results
    """
    quotients = []
    for value in tuple1:
        quotients.append(int(round(value/divisor)))
    return tuple(quotients)
Any information on how to incorporate the described boolean matrix, or any other ideas on how to speed this up, would be greatly appreciated.
