Адаптивный алгоритм порогового значения 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 или 6 секунд при использовании размера окна 5, но если я уменьшу размер окна до 20 и алгоритм проверяет 20 пикселей в каждом направлении для среднего значения, я получаю время вверх одна минута для этого изображения.

Если я использую изображение размером 2592×1936 с размером окна всего 5, это займет почти 10 минут.

Итак, как я могу улучшить те времена? Будет ли более массивный массив? Является ли im.getpixel быстрее, чем загрузка изображения в режим доступа к пикселям? Есть ли еще какие-нибудь советы по повышению скорости? Заранее спасибо.

Обращаясь к нашим комментариям, я написал MATLAB реализацию этого алгоритма здесь: Извлеките страницу из однородного фона в изображении , и она была довольно быстрой на больших изображениях.

Если вы хотите лучше объяснить алгоритм, см. Мой другой ответ здесь: Bradley Adaptive Thresholding – Confused (вопросы) . Это может быть хорошим местом для начала, если вы хотите лучше понять код, который я написал.

Поскольку 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 # 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 раз на загруженном вами изображении, лучшее из трех времен выполнения дало в среднем 7,68 миллисекунды за каждый прогон.

Я также получаю это изображение в результате, когда я его порою:

введите описание изображения здесь

Профилирование вашего кода в IPython с помощью %prun yields показывает:

  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, если не можете определить, как векторизовать ваш код.

Хорошо, я немного опоздал. Позвольте мне поделиться своими мыслями по этому поводу:

Вы можете ускорить его, используя динамическое программирование, чтобы вычислить средства, но гораздо проще и быстрее позволить scipy и numpy делать всю грязную работу. (Обратите внимание, что я использую Python3 для моего кода, поэтому xrange изменен на диапазон в вашем коде).

 #!/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, немного отличается на границах, чем у вашего кода (для позиций, где окно для среднего вычисления больше не содержится в его изображении). Однако, я думаю, это не должно быть проблемой.

Еще одно незначительное отличие состоит в том, что окно из циклов for-loop не было точно сосредоточено на пикселе, поскольку смещение по xrange (-ws, ws) с ws = 5 дает -5, -4-, …, 3,4 и приводит к среднему значению -0,5. Вероятно, это не было предназначено.