Bradley-Roth 自适应阈值算法 - 如何获得更好的性能？
Bradley-Roth Adaptive Thresholding Algorithm - How do I get better performance?
我有以下图像阈值处理代码，使用 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 的 window 时，我遇到了大约 5 或 6 秒的处理时间，但是如果我将我的 window 大小增加到 20，并且算法正在检查 20 个像素在平均值的每个方向上，我得到该图像的时间超过一分钟。
如果我用一张像2592x1936这样的图片，window尺寸只有5，需要将近10分钟才能完成。
那么，我该如何改进那些时间呢？ numpy 数组会更快吗？ im.getpixel 是否比将图像加载到像素访问模式更快？还有其他提高速度的技巧吗？提前致谢。
使用 %prun 对 IPython 中的代码进行性能分析显示：
ncalls tottime percall cumtime percall filename:lineno(function)
50246 2.009 0.000 2.009 0.000 <ipython-input-78-b628a43d294b>:15(<listcomp>)
50246 0.587 0.000 0.587 0.000 <ipython-input-78-b628a43d294b>:17(<listcomp>)
1 0.170 0.170 2.829 2.829 <ipython-input-78-b628a43d294b>:5(bradley_threshold)
50246 0.058 0.000 0.058 0.000 {built-in method sum}
50257 0.004 0.000 0.004 0.000 {built-in method len}
即几乎所有 运行 时间都是由于 Python 循环（慢）和非矢量化算术（慢）。因此，如果您使用 numpy 数组重写，我希望会有很大的改进；或者，如果您不知道如何对代码进行矢量化，则可以使用 cython。
参考我们的评论，我在这里编写了该算法的 MATLAB 实现：，它在大图像上速度相当快。
如果您想对算法有更好的解释，请在此处查看我的其他答案：。如果您想更好地理解我编写的代码，这可能是一个很好的起点。
因为 MATLAB 和 NumPy 相似，所以这是 Bradley-Roth 阈值算法的重新实现，但是在 NumPy 中。我将 PIL 图像转换为 NumPy 数组，对该图像进行处理，然后再转换回 PIL 图像。该函数接受三个参数：灰度图像 image、window 的大小 s 和阈值 t。这个阈值与你所拥有的不同，因为它完全遵循论文。阈值 t 是每个像素 window 总面积的 百分比 。如果总面积小于此阈值，则输出应为黑色像素 - 否则为白色像素。 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 测量了平均执行时间，这就是我在 post:
中上传的图片得到的结果
In [16]: %timeit bradley_roth_numpy(img)
100 loops, best of 3: 7.68 ms per loop
这意味着 运行在您上传的图片上重复此函数 100 次，最好的 3 次执行时间平均每 运行.
7.68 毫秒
当我对它进行阈值处理时，我也得到了这张图片：
好的，我来晚了。无论如何，让我分享一下我的想法：
您可以通过使用动态编程来计算均值来加快速度，但让 scipy 和 numpy 完成所有脏工作会更容易和更快。 （请注意，我的代码使用 Python3，因此 xrange 在您的代码中更改为 range）。
#!/usr/bin/env python3
import numpy as np
from scipy import ndimage
from PIL import Image
import copy
import time
def faster_bradley_threshold(image, threshold=75, window_r=5):
percentage = threshold / 100.
window_diam = 2*window_r + 1
# convert image to numpy array of grayscale values
img = np.array(image.convert('L')).astype(np.float) # float for mean precision
# matrix of local means with scipy
means = ndimage.uniform_filter(img, window_diam)
# result: 0 for entry less than percentage*mean, 255 otherwise
height, width = img.shape[:2]
result = np.zeros((height,width), np.uint8) # initially all 0
result[img >= percentage * means] = 255 # numpy magic :)
# convert back to PIL image
return Image.fromarray(result)
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 range(h):
for x in range(w):
#find neighboring pixels
neighbors =[(x+x2,y+y2) for x2 in range(-ws,ws) for y2 in range(-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
if __name__ == '__main__':
img = Image.open('test.jpg')
t0 = time.process_time()
threshed0 = bradley_threshold(img)
print('original approach:', round(time.process_time()-t0, 3), 's')
threshed0.show()
t0 = time.process_time()
threshed1 = faster_bradley_threshold(img)
print('w/ numpy & scipy :', round(time.process_time()-t0, 3), 's')
threshed1.show()
这使它在我的机器上更快：
$ python3 bradley.py
original approach: 3.736 s
w/ numpy & scipy : 0.003 s
PS：请注意，我在 scipy 中使用的均值在边界处的行为与您的代码中的均值略有不同（对于均值计算的 window 不是的位置完全包含在他的图像中）。但是，我认为这应该不是问题。
另一个小区别是来自 for 循环的 window 并不完全以像素为中心，因为 xrange(-ws,ws) 和 ws=5 的偏移量产生 -5,-4- ,...,3,4 的平均值为 -0.5。这可能不是故意的。
我有以下图像阈值处理代码，使用 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 的 window 时，我遇到了大约 5 或 6 秒的处理时间，但是如果我将我的 window 大小增加到 20，并且算法正在检查 20 个像素在平均值的每个方向上，我得到该图像的时间超过一分钟。
如果我用一张像2592x1936这样的图片，window尺寸只有5，需要将近10分钟才能完成。
那么，我该如何改进那些时间呢？ numpy 数组会更快吗？ im.getpixel 是否比将图像加载到像素访问模式更快？还有其他提高速度的技巧吗？提前致谢。
使用 %prun 对 IPython 中的代码进行性能分析显示：
ncalls tottime percall cumtime percall filename:lineno(function)
50246 2.009 0.000 2.009 0.000 <ipython-input-78-b628a43d294b>:15(<listcomp>)
50246 0.587 0.000 0.587 0.000 <ipython-input-78-b628a43d294b>:17(<listcomp>)
1 0.170 0.170 2.829 2.829 <ipython-input-78-b628a43d294b>:5(bradley_threshold)
50246 0.058 0.000 0.058 0.000 {built-in method sum}
50257 0.004 0.000 0.004 0.000 {built-in method len}
即几乎所有 运行 时间都是由于 Python 循环（慢）和非矢量化算术（慢）。因此，如果您使用 numpy 数组重写，我希望会有很大的改进；或者，如果您不知道如何对代码进行矢量化，则可以使用 cython。
参考我们的评论，我在这里编写了该算法的 MATLAB 实现：
如果您想对算法有更好的解释，请在此处查看我的其他答案：
因为 MATLAB 和 NumPy 相似，所以这是 Bradley-Roth 阈值算法的重新实现，但是在 NumPy 中。我将 PIL 图像转换为 NumPy 数组，对该图像进行处理，然后再转换回 PIL 图像。该函数接受三个参数：灰度图像 image、window 的大小 s 和阈值 t。这个阈值与你所拥有的不同，因为它完全遵循论文。阈值 t 是每个像素 window 总面积的 百分比 。如果总面积小于此阈值，则输出应为黑色像素 - 否则为白色像素。 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 测量了平均执行时间，这就是我在 post:
In [16]: %timeit bradley_roth_numpy(img)
100 loops, best of 3: 7.68 ms per loop
这意味着 运行在您上传的图片上重复此函数 100 次，最好的 3 次执行时间平均每 运行.
7.68 毫秒当我对它进行阈值处理时，我也得到了这张图片：
好的，我来晚了。无论如何，让我分享一下我的想法：
您可以通过使用动态编程来计算均值来加快速度，但让 scipy 和 numpy 完成所有脏工作会更容易和更快。 （请注意，我的代码使用 Python3，因此 xrange 在您的代码中更改为 range）。
#!/usr/bin/env python3
import numpy as np
from scipy import ndimage
from PIL import Image
import copy
import time
def faster_bradley_threshold(image, threshold=75, window_r=5):
percentage = threshold / 100.
window_diam = 2*window_r + 1
# convert image to numpy array of grayscale values
img = np.array(image.convert('L')).astype(np.float) # float for mean precision
# matrix of local means with scipy
means = ndimage.uniform_filter(img, window_diam)
# result: 0 for entry less than percentage*mean, 255 otherwise
height, width = img.shape[:2]
result = np.zeros((height,width), np.uint8) # initially all 0
result[img >= percentage * means] = 255 # numpy magic :)
# convert back to PIL image
return Image.fromarray(result)
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 range(h):
for x in range(w):
#find neighboring pixels
neighbors =[(x+x2,y+y2) for x2 in range(-ws,ws) for y2 in range(-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
if __name__ == '__main__':
img = Image.open('test.jpg')
t0 = time.process_time()
threshed0 = bradley_threshold(img)
print('original approach:', round(time.process_time()-t0, 3), 's')
threshed0.show()
t0 = time.process_time()
threshed1 = faster_bradley_threshold(img)
print('w/ numpy & scipy :', round(time.process_time()-t0, 3), 's')
threshed1.show()
这使它在我的机器上更快：
$ python3 bradley.py
original approach: 3.736 s
w/ numpy & scipy : 0.003 s
PS：请注意，我在 scipy 中使用的均值在边界处的行为与您的代码中的均值略有不同（对于均值计算的 window 不是的位置完全包含在他的图像中）。但是，我认为这应该不是问题。
另一个小区别是来自 for 循环的 window 并不完全以像素为中心，因为 xrange(-ws,ws) 和 ws=5 的偏移量产生 -5,-4- ,...,3,4 的平均值为 -0.5。这可能不是故意的。