Bradley-Roth自适应阈值算法 - 如何获得更好的性能？

Question

我有以下代码用于图像阈值处理,使用Bradley-Roth图像阈值处理方法.

from PIL import Image
import copy
import time
def bradley_threshold(image, threshold=75, windowsize=5):
    ws = windowsize
    image2 = copy.copy(image).convert('L')
    w, h = image.size
    l = image.convert('L').load()
    l2 = image2.load()
    threshold /= 100.0
    for y in xrange(h):
        for x in xrange(w):
            #find neighboring pixels
            neighbors =[(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 pixels
            mean = sum([l[a,b] for a,b in neighbors])/len(neighbors)
            if l[x, y] < threshold*mean:
                l2[x,y] = 0
            else:
                l2[x,y] = 255
    return image2

i = Image.open('test.jpg')
windowsize = 5
bradley_threshold(i, 75, windowsize).show()

当windowsize小而图像小时,这可以正常工作.我一直在使用这个图像进行测试:

当使用5的窗口大小时,我经历了大约5或6秒的处理时间,但是如果我将窗口大小提高到20并且算法检查每个方向上20个像素的平均值,我得到的时间超过一分钟的那个图像.

如果我使用尺寸为2592x1936且窗口大小仅为5的图像,则需要将近10分钟才能完成.

那么,我该如何改善这些时间呢？numpy数组会更快吗？im.getpixel比将图像加载到像素访问模式更快吗？还有其他提速提示吗？提前致谢.

Answer 1

参考我们的评论,我在这里编写了这个算法的MATLAB实现:从图像中的统一背景中提取页面,在大图像上它非常快.

如果您想要更好地解释算法,请在此处查看我的其他答案:Bradley Adaptive Thresholding - Confused(questions).如果您想要更好地理解我编写的代码,这可能是一个很好的起点.

因为MATLAB和NumPy是相似的,所以这是Bradley-Roth阈值算法的重新实现,但是在NumPy中.我将PIL图像转换为NumPy数组,对此图像进行处理,然后转换回PIL图像.该函数包含三个参数:灰度图像image,窗口大小s和阈值t.这个阈值与你所拥有的不同,因为这完全符合本文的要求.阈值t是每个像素窗口的总求和面积的百分比.如果求和面积小于此阈值,则输出应为黑色像素 - 否则为白色像素.的默认值s和t是除以8并四舍五入列的分别的数量,和15%:

import numpy as np
from PIL import Image

def bradley_roth_numpy(image, s=None, t=None):

    # Convert image to numpy array
    img = 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 window
    if t is None:
        t = 15.0

    # Compute integral image
    intImage = 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 access
    X = X.ravel()
    Y = Y.ravel()

    # Ensure s is even so that we are able to index into the image
    # properly
    s = s + np.mod(s,2)

    # Access the four corners of each neighbourhood
    x1 = X - s/2
    x2 = X + s/2
    y1 = Y - s/2
    y2 = Y + s/2

    # Ensure no coordinates are out of bounds
    x1[x1 < 0] = 0
    x2[x2 >= cols] = cols-1
    y1[y1 < 0] = 0
    y2[y2 >= rows] = rows-1

    # Ensures coordinates are integer
    x1 = 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 neighbourhood
    count = (x2 - x1) * (y2 - y1)

    # Compute the row and column coordinates to access
    # each corner of the neighbourhood for the integral image
    f1_x = x2
    f1_y = y2
    f2_x = x2
    f2_y = y1 - 1
    f2_y[f2_y < 0] = 0
    f3_x = x1-1
    f3_x[f3_x < 0] = 0
    f3_y = y2
    f4_x = f3_x
    f4_y = f2_y

    # Compute areas of each window
    sums = 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 grid
    out = np.ones(rows*cols, dtype=np.bool)
    out[img.ravel()*count <= sums*(100.0 - t)/100.0] = False

    # Also convert back to uint8
    out = 255*np.reshape(out, (rows, cols)).astype(np.uint8)

    # Return PIL image back to user
    return Image.fromarray(out)


if __name__ == '__main__':
    img = Image.open('test.jpg').convert('L')
    out = bradley_roth_numpy(img)
    out.show()
    out.save('output.jpg')

如果需要,读入图像并将其转换为灰度.将显示输出图像,它将保存到您将脚本运行到调用图像的同一目录中output.jpg.如果要覆盖设置,只需执行以下操作:

out = bradley_roth_numpy(img, windowsize, threshold)

玩这个以获得良好的结果.使用默认参数并使用IPython,我测量了使用的平均执行时间timeit,这是我在您的帖子中上传的图像所得到的:

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

这意味着在您上传的图像上重复运行此功能100次,最佳的3次执行时间平均每次运行7.68毫秒.

当我对其进行阈值处理时,我也会得到这个图像: