Slow python image processing with PIL and numpy

Question

I'm trying to implement some image processing (finding regions of similar colour) in Python with PIL and Numpy. Can't figure out how to speed up this code. Could you help?

def findRegions(self, data):
#data is numpy.array
    ret = [[False for _ in range(self.width)] for _ in range(self.heigth)]

    for i in range(self.heigth):
        for j in range(self.width):
            k = 0
            acc = 0
            for x,y in [(-1,0),(0,-1),(0,1),(1,0)]:
                if (self.heigth>i+x>=0 and self.width>j+y>=0):
                    k = k+1
                    acc += math.sqrt(sum((data[i][j][c]-data[i+x][j+y][c])**2 for c in range(3)))
            if (acc/k<self.threshold):
                ret[i][j]= True
    return ret

PIL and other image libraries have got many filtering and processing functions which are really quick. But what is the best way to implement own image processing functions?

Get rid of the nested loops, and operate on the whole numpy array at once rather than per-pixel. — wim, Feb 02 '12 at 00:24
Use PyPy http://pypy.org/ and PIL http://stackoverflow.com/questions/4864674/how-to-use-pil-with-pypy — Mikko Ohtamaa, Feb 02 '12 at 01:13
Also look into scipy.ndimage It can be quite useful for this kind of thing. — Bi Rico, Feb 02 '12 at 06:10

score 4 · Accepted Answer · edited May 23 '17 at 12:16

Rather than looping over each row and column you can shift the array left, right, up, and down for the appropriate number of elements. On each shift you accumulate your values in a base array. After the shifting and accumulating you compute your average and apply your threshold to return a mask. See this post which has a general discussion on the topic. The idea is take advantage of numpy's broadcasting, which will apply a function or operator to all elements of an array in C rather than Python.

I've adapted the code from the linked post to fit what I believe you are trying to accomplish. In any case the general pattern should speed things up. You have to work out what to do with the edges in the return mask. Here I've simply set the return mask to False, but you could also eliminate the edges by expanding the input data by one pixel in each direction and filling with the nearest pixel, zeros, gray, etc.

def findRegions(self,data):
    #define the shifts for the kernel window
    shifts = [(-1,0),(0,-1),(0,1),(1,0)]

    #make the base array of zeros 
    #  array size by 2 in both dimensions
    acc = numpy.zeros(data.shape[:2])

    #compute the square root of the sum of squared color 
    # differences between a pixel and it's 
    # four cardinal neighbors
    for dx,dy in shifts:
        xstop = -1+dx or None
        ystop = -1+dy or None
        #per @Bago's comment, use the sum method to add up the color dimension
        #  instead of the list comprehension
        acc += ((data[1:-1,1:-1] - data[1+dx:xstop, 1+dy:ystop])**2).sum(-1)**.5

    #compute the average 
    acc /= (len(shifts) + 1)

    #build a mask array the same size as the original
    ret = numpy.zeros(data.shape[:2],dtype=numpy.bool)

    #apply the threshold
    #  note that the edges will be False
    ret[1:-1,1:-1] acc < self.threshold    

    return ret

You can use the arrays sum method to avoid looping over the color dimension, ie `((data[1:-1, 1:-1] - data[1+dx:xstop, 1+dy:ystop])**2).sum(-1)**.5` — Bi Rico, Feb 02 '12 at 03:04

score 2 · Answer 2 · answered Feb 07 '12 at 09:57

There are better segmentation algorithms included in http://scikits-image.org , but if you want to build your own you can look at this example, based on clustering, called ICM segmentation. Specify N=4 to identify four regions.

import numpy as np
from scipy.cluster.vq import kmeans2

def ICM(data, N, beta):
    print "Performing ICM segmentation..."

    # Initialise segmentation using kmeans
    print "K-means initialisation..."
    clusters, labels = kmeans2(np.ravel(data), N)

    print "Iterative segmentation..."
    f = data.copy()

    def _minimise_cluster_distance(data, labels, N, beta):
        data_flat = np.ravel(data)
        cluster_means = np.array(
            [np.mean(data_flat[labels == k]) for k in range(N)]
            )
        variance = np.sum((data_flat - cluster_means[labels])**2) \
                   / data_flat.size

        # How many of the 8-connected neighbouring pixels are in the
        # same cluster?
        count = np.zeros(data.shape + (N,), dtype=int)
        count_inside = count[1:-1, 1:-1, :]

        labels_img = labels.reshape(data.shape)
        for k in range(N):
            count_inside[..., k] += (k == labels_img[1:-1:, 2:])
            count_inside[..., k] += (k == labels_img[2:, 1:-1])
            count_inside[..., k] += (k == labels_img[:-2, 1:-1])
            count_inside[..., k] += (k == labels_img[1:-1, :-2])

            count_inside[..., k] += (k == labels_img[:-2, :-2])
            count_inside[..., k] += (k == labels_img[2:, 2:])
            count_inside[..., k] += (k == labels_img[:-2, 2:])
            count_inside[..., k] += (k == labels_img[2:, :-2])

        count = count.reshape((len(labels), N))
        cluster_measure = (data_flat[:, None] - cluster_means)**2 \
                          - beta * variance * count
        labels = np.argmin(cluster_measure, axis=1)

        return cluster_means, labels

    # Initialise segmentation
    cluster_means, labels = _minimise_cluster_distance(f, labels, N, 0)

    stable_counter = 0
    old_label_diff = 0
    i = 0
    while stable_counter < 3:
        i += 1

        cluster_means, labels_ = \
                       _minimise_cluster_distance(f, labels, N, beta)

        new_label_diff = np.sum(labels_ != labels)
        if  new_label_diff != old_label_diff:
            stable_counter = 0
        else:
            stable_counter += 1
        old_label_diff = new_label_diff

        labels = labels_

    print "Clustering converged after %d steps." % i

    return labels.reshape(data.shape)

Slow python image processing with PIL and numpy

2 Answers2