Generating boolean matrix from image

Question

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.

Answer 1

Just load the image as a numpy array, and then use array operations instead of looping over pixels:

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

def zone_map(picture, threshold, show=True):

    with PIL.Image.open(picture) as img:
        rgb = np.array(img, dtype=np.float)

    height, width, _ = rgb.shape

    mask = np.zeros_like(rgb)
    while not np.any(mask):

        # get random pixel
        position = np.random.randint(height), np.random.randint(width)
        color = rgb[position]

        # get euclidean distance of all pixels in colour space
        distance = np.sqrt(np.sum((rgb - color)**2, axis=-1))

        # threshold
        mask = distance < threshold
        if show: # show output
            fig, (ax1, ax2) = plt.subplots(1,2)
            ax1.imshow(rgb.astype(np.uint8))
            ax2.imshow(mask, cmap='gray')
            fig.suptitle('Random color: {}'.format(color))

    return mask

def test():
    zone_map("Lenna.jpg", threshold=20)
    plt.show()