I have the following code for image thresholding, using the Bradley-Roth image thresholding method.

from PIL import Image
import copy
import time
def bradley_threshold(image, threshold=75, windowsize=5):ws = windowsizeimage2 = copy.copy(image).convert('L')w, h = image.sizel = image.convert('L').load()l2 = image2.load()threshold /= 100.0for y in xrange(h):for x in xrange(w):#find neighboring pixelsneighbors =[(x+x2,y+y2) for x2 in xrange(-ws,ws) for y2 in xrange(-ws, ws) if x+x2>0 and x+x2<w and y+y2>0 and y+y2<h]#mean of all neighboring pixelsmean = sum([l[a,b] for a,b in neighbors])/len(neighbors)if l[x, y] < threshold*mean:l2[x,y] = 0else:l2[x,y] = 255return image2i = Image.open('test.jpg')
windowsize = 5
bradley_threshold(i, 75, windowsize).show()

This works fine when windowsize is small and the image is small. I've been using this image for testing:

I'm experiencing processing times of about 5 or 6 seconds when using a window size of 5, but if I bump my window size up to 20 and the algorithm is checking 20 pixels in each direction for the mean value, i get times upwards of one minute for that image.

If I use an image with a size like 2592x1936 with a window size of only 5, it takes nearly 10 minutes to complete.

So, how can I improve those times? Would a numpy array be faster? Is im.getpixel faster than loading the image into pixel access mode? Are there any other tips for speed boosts? Thanks in advance.

Referencing our comments, I wrote a MATLAB implementation of this algorithm here: Extract a page from a uniform background in an image, and it was quite fast on large images.

If you'd like a better explanation of the algorithm, please see my other answer here: Bradley Adaptive Thresholding -- Confused (questions). This may be a good place to start if you want a better understanding of the code I wrote.

Because MATLAB and NumPy are similar, this is a re-implementation of the Bradley-Roth thresholding algorithm, but in NumPy. I convert the PIL image into a NumPy array, do the processing on this image, then convert back to a PIL image. The function takes in three parameters: the grayscale image image, the size of the window s and the threshold t. This threshold is different than what you have as this is following the paper exactly. The threshold t is a percentage of the total summed area of each pixel window. If the summed area is less than this threshold, then the output should be a black pixel - else it's a white pixel. The defaults for s and t are the number of columns divided by 8 and rounded, and 15% respectively:

import numpy as np
from PIL import Imagedef bradley_roth_numpy(image, s=None, t=None):# Convert image to numpy arrayimg = np.array(image).astype(np.float)# Default window size is round(cols/8)if s is None:s = np.round(img.shape[1]/8)# Default threshold is 15% of the total# area in the windowif t is None:t = 15.0# Compute integral imageintImage = np.cumsum(np.cumsum(img, axis=1), axis=0)# Define grid of points(rows,cols) = img.shape[:2](X,Y) = np.meshgrid(np.arange(cols), np.arange(rows))# Make into 1D grid of coordinates for easier accessX = X.ravel()Y = Y.ravel()# Ensure s is even so that we are able to index into the image# properlys = s + np.mod(s,2)# Access the four corners of each neighbourhoodx1 = X - s/2x2 = X + s/2y1 = Y - s/2y2 = Y + s/2# Ensure no coordinates are out of boundsx1[x1 < 0] = 0x2[x2 >= cols] = cols-1y1[y1 < 0] = 0y2[y2 >= rows] = rows-1# Ensures coordinates are integerx1 = x1.astype(np.int)x2 = x2.astype(np.int)y1 = y1.astype(np.int)y2 = y2.astype(np.int)# Count how many pixels are in each neighbourhoodcount = (x2 - x1) * (y2 - y1)# Compute the row and column coordinates to access# each corner of the neighbourhood for the integral imagef1_x = x2f1_y = y2f2_x = x2f2_y = y1 - 1f2_y[f2_y < 0] = 0f3_x = x1-1f3_x[f3_x < 0] = 0f3_y = y2f4_x = f3_xf4_y = f2_y# Compute areas of each windowsums = intImage[f1_y, f1_x] - intImage[f2_y, f2_x] - intImage[f3_y, f3_x] + intImage[f4_y, f4_x]# Compute thresholded image and reshape into a 2D gridout = np.ones(rows*cols, dtype=np.bool)out[img.ravel()*count <= sums*(100.0 - t)/100.0] = False# Also convert back to uint8out = 255*np.reshape(out, (rows, cols)).astype(np.uint8)# Return PIL image back to userreturn Image.fromarray(out)if __name__ == '__main__':img = Image.open('test.jpg').convert('L')out = bradley_roth_numpy(img)out.show()out.save('output.jpg')

The image is read in and converted to grayscale if required. The output image will be displayed, and it will be saved to the same directory where you ran the script to an image called output.jpg. If you want to override the settings, simply do:

out = bradley_roth_numpy(img, windowsize, threshold)

Play around with this to get good results. Using the default parameters and using IPython, I measured the average time of execution using timeit, and this is what I get for your image you uploaded in your post:

In [16]: %timeit bradley_roth_numpy(img)
100 loops, best of 3: 7.68 ms per loop

This means that running this function repeatedly 100 times on the image you uploaded, the best of 3 execution times gave on average 7.68 milliseconds per run.

I also get this image as a result when I threshold it: