對numpy的一系列實驗

05-09

本文主要總結了如何把numpy計算速度優化100倍的過程。

簡介

最近要對圖片做二值化，參照Adaptive Thresholding Using the Integral Image 這篇論文，可是Python的for循環非常的慢，單純的翻譯偽代碼實現起來速度要比C慢接近100倍，論文上的實現對640x480的圖片用時大概為15ms，接下來我們要一步步從最初的慢100倍優化到和C一樣快，以及最後再比它快一點點。

演算法

首先，介紹一下 Adaptive Thresholding Using the Integral Image 這篇論文，論文很簡單，對一個點周圍的一片區域取平均值，如果該點低於平均值的85%，就設為黑，否則為白色。論文主要是用累積和來改進取一片區域平均值的速度。

過程為水平方向做一次累積和，垂直方向做一次累積和，這樣每一個點就表示其左上角的區域的和，於是要計算一個矩形（a, b, c, d), 其區域上的和就是 d - b -c + a.

原始的掃描區域的每一個點求和，時間是O(n * n * boxsize * boxsize)，現在變成了常數時間。

直觀實現

首先我們直接照著原理實現，理論是最快的

（然而實際是最慢的）

def intuitive_method(mat, w, h , set_value, threshold, boxsize): mat = mat.T integral_map = np.zeros_like(mat, dtype=int) ret = np.zeros_like(mat, dtype=int) for i in range(w): csum=0 for j in range(h): csum += mat[i, j] if i > 0: integral_map[i,j] = integral_map[i-1, j] + csum else: integral_map[i,j] = csum half_s = (boxsize-1) // 2 for i in range(w): for j in range(h): x1 = max(i - half_s, 0) x2 = min(i+half_s, w - 1) y1 = max(j-half_s, 0) y2 = min(j + half_s, h - 1) field_sum = integral_map[x2, y2] - integral_map[x1, y2] - integral_map[x2, y1] + integral_map[x1, y1] count = (x2-x1)*(y2-y1) if (mat[i,j] * count)*100 < field_sum*(100-threshold): ret[i,j] = set_value else: ret[i,j] = set_value ^ 255 return ret.astype(np.uint8).T

mat是我們的圖片矩陣， boxsize是區域的長度，必須是奇數。

每個點都在周圍取一個boxsize長度的矩形，要特別處理邊界，讓他不要超出累積矩陣。

很直觀的代碼，純翻譯自論文裡面的偽代碼。

看看速度：

? python3 implement.pyImage shape 640x480intuitive_method take 1658.280885 ms

1600ms, （怎麼和說好的15ms不一樣）。

改進的直觀方法

def optimized_intuitive_method(mat, w, h , set_value, threshold, boxsize): mat = mat.T integral_map = np.zeros_like(mat, dtype=int) ret = np.zeros_like(mat, dtype=int) np.cumsum(mat, 0, out=integral_map) np.cumsum(integral_map, 1, out=integral_map) half_s = boxsize >> 1 inv_value = set_value ^ 255 for i in range(w): for j in range(h): x1 = max(i - half_s, 0) x2 = min(i + half_s, w - 1) y1 = max(j - half_s, 0) y2 = min(j + half_s, h - 1) field_sum = integral_map.item(x2, y2) - integral_map.item(x1, y2) - integral_map.item(x2, y1) + integral_map.item(x1, y1) count = (x2-x1)*(y2-y1) if (mat.item(i,j) * count)*100 < field_sum*(100-threshold): ret.itemset(i, j, set_value) else: ret.itemset(i, j, inv_value) return ret.astype(np.uint8).T

Numpy提供了item和itemset方法來加速對矩陣上單個點的操作。

同時兩次累積和實際可以用可以用np.cumsum來計算，但是這裡有一個小細節，

單純的使用兩次a=np.cumsum(...)會額外進行一次內部的內存拷貝，我們要指定out參數來阻止多餘的內存拷貝.

在計算half_s的時候，使用移位boxsize >> 1來加速。

? python3 implement.pyImage shape 640x480optimized_intuitive_method take 823.928799 ms

快了一倍，但是還是不夠

猥瑣法

不妨讓我們直接退化到累積和的優化之前，純粹的暴力計算，所有計算都交給opencv，看看會發生什麼事情。

def opencv_way(mat, w, h , set_value, threshold, boxsize): mat = mat.astype(int) inv_value = set_value ^ 255 kernel = np.ones((boxsize, boxsize)) * (100-threshold)/100 / (boxsize*boxsize) result = cv.filter2D(mat, -1, kernel) ret = np.where(mat < result, [set_value], [inv_value]) return ret.astype(np.uint8)

這是時間：

? python3 implement.pyImage shape 640x480opencv_way take 35.920926 ms

好了，可以全劇終了。

簡單的6行代碼，純粹的暴力美學。速度快了40倍。（你學到了么？）

這裡我們藉助opencv 的 filter2D方法，自己定義了一個kernel，用來求一個點周圍區域的平均值，接下來的計算就全部交給opencv了。

關於opencv，有一點需要補充， opencv表示， numpy是辣雞，請多使用opencv自帶的函數。

向量化

還不能全劇終，接下來是重點， numpy是為矩陣計算而生的，我們前面一直按照C的思維來對單個點做操作，而且for循環的存在導致速度下降太多，接下來我們要尋求把for循環替換掉，改成向量的計算。

對於numpy來說，比如我有兩個2維矩陣A和B，要計算兩個矩陣上所有點的差，從速度上來說，

A - B > A[i] - B[i] > A[i,j] - B[i,j]

我們現在希望把所有的for循環，替換為矩陣的計算 A - B

一開始的時候我想到利用一次for循環，從上往下計算兩行之間的值，接著從左往右計算每一列的值，但是這樣只能去掉一個for循環，後來我想了2個晚上，想通了實際上所有等距離點操作都可以轉化為等距離的向量操作，進而可以通過構造矩陣來進行矩陣間的操作。

對於之前的代碼，計算累積和部分已經很好了，關鍵之後的計算，代碼冗餘，效率低下。

其中造成整個for循環存在的關鍵因素是下面這句代碼

integral_map[x2, y2] - integral_map[x1, y2] - integral_map[x2, y1] + integral_map[x1, y1]

而實際上，這些點之間是等距離的，（boxsize的距離），也就是說，每一個點對應的x1， x2, y1, y2 都是固定的，我們可以構造出4個矩陣，分別向左上，右上，左下，右下位移一段距離，

這樣就可以吧整個for循環計算化簡成

D + A - C - B

那麼怎麼構造出這些矩陣呢？

想像一下，把integral_map[x1, y1]當成一個集合，記為A， integral_map[x1, y2]的集合記為B， integral_map[x2, y1]的集合記為C，integral_map[x2, y2]記為D。

以構造A為例， A實際上是對原先的矩陣向上和向左平移了一半的boxsize的距離，比如對於矩陣右下角的點，A對應於 integral_map[h-half_boxsize, w-half_boxsize]

但是對於左上角的點，超出了累積矩陣的範圍，為了避免額外的對邊界的判斷，我們要額外構造出邊框。

想像一下，比如我們有一個5x5的矩陣， boxsize取3,

0 1 2 3 40 * * * * *1 * * * * *2 * * * * *3 * * * * *4 * * * * *

計算 (2, 2) 的時候，取 (1, 1), (1, 3)，（3，1），（3，3）四個點

計算(0,0)上的和的時候，右下角取（1，1），左上角超出了邊界，所以只能取（0，0）

計算（0，1）的時候，右下角是（1，2）但是左上角還是（0，0）

所以我們可以構造出一個(boxsize-1)/2 長度的邊框，邊框上的值和邊界的值一樣。這樣我們可以通過平移得到ABCD四個矩陣，而不用額外的邊界判斷。

def vector_method(mat, w, h, set_value, threshold, boxsize): inv_value = set_value ^ 255 integral_map = np.zeros_like(mat, dtype=int) np.cumsum(mat,0, out=integral_map) np.cumsum(integral_map, 1, out=integral_map) half_s = boxsize >> 1 big_mat = cv.copyMakeBorder(integral_map, half_s, half_s, half_s, half_s, borderType=cv.BORDER_REPLICATE) big_h, big_w = big_mat.shape # 構造4個矩陣 # top left A = big_mat[:h, :w] # bottom left B = big_mat[big_h - h:, :w] # top right C = big_mat[:h, big_w - w:] # bottom right D = big_mat[big_h - h:, big_w - w:] COUNT = get_count_matrix(w, h, boxsize) result = (D + A - B - C)*((100-threshold)/100) / COUNT ret = np.where(mat < result, [set_value], [inv_value]) return ret.astype(np.uint8)

但是，現在還有一道陰影籠罩著我們，我們沒有辦法得到窗口的實際大小。邊界附近的點上的窗口是變化的，我們得額外構造出一個矩陣COUNT來表示每個點上的窗口大小。

def get_count_matrix(w, h, boxsize): """ 矩陣上每個點的值表示該點上的窗口覆蓋的總點數。 """ # faster than (boxsize-1) // 2 half_boxsize = (boxsize) >> 1 # 窗口的真正大小並不是boxsize！而是（boxsize - 1）！ A = np.ones((h, w)) * (boxsize-1) B = np.ones((h, w)) * (boxsize-1) for i in range(half_boxsize): value = i + half_boxsize A[i, :] = value A[-1-i, :] = value B[:, i] = value B[:, -1-i] = value return A * B

看一下所需的時間

? python3 implement.pyImage shape 640x480vector_method take 22.166661 ms

現在我們比最初快了80倍，已經接近論文上的理論速度了，但是get_count_matrix函數顯得非常的力不從心，還可以尋找優化的辦法。

優化後的向量化方法

回想之前的整個過程，我們對邊界的處理實際上過於保守了，導致計算窗口覆蓋點數的時候過於小心翼翼了(不然boxsize-1是怎麼得出來的，都是淚啊)。

我們可以在計算累積和之前先對邊界翻折來進行填充，之後再對擴展後的矩陣進行累積和，這樣可以省略掉get_count_matrix函數，統一採用(boxsize-1)*(boxsize-1)。

def optimized_vector_method(mat, w, h, set_value, threshold, boxsize): inv_value = set_value ^ 255 half_s = boxsize >> 1 integral_map = cv.copyMakeBorder(mat, half_s, half_s, half_s, half_s, borderType=cv.BORDER_REFLECT101) integral_map = integral_map.astype(np.uint32) integral_map.cumsum(0, out=integral_map) integral_map.cumsum(1, out=integral_map) big_h, big_w = integral_map.shape # 構造4個矩陣 # top left A = integral_map[:h, :w] # bottom left B = integral_map[big_h - h:, :w] # top right C = integral_map[:h, big_w - w:] # bottom right D = integral_map[big_h - h:, big_w - w:] # COUNT = get_count_matrix(w, h, boxsize) COUNT = (boxsize-1)**2 result = (D + A - B - C) * ((100-threshold) / 100) / COUNT ret = np.where(mat < result, [set_value], [inv_value]) return ret.astype(np.uint8)

看看計算時間，

? python3 implement.pyImage shape 640x480optimized_vector_method take 10.539570999999999 ms

比最初快了150倍，而且超過了論文上的C語言實現的15ms的時間。

總結分析

對於最後的這個方法，需要額外解釋一下為什麼會比論文快。

其中一個是CPU的原因，論文上用的是4核奔騰3.4Ghz的CPU，而我自己的是i5 2.7GHz。

i5應該是快一些的，（我猜）

之後是在對窗口大小的計算時間上，因為我們通過填充擴展矩陣，使得每一個點的窗口大小都是固定的值，比論文中的方法省去了多餘的窗口大小計算和多餘的邊界判斷。

但是擴展矩陣在計算累積和的時候有額外的時間消耗， boxsize默認採用邊長的1/8，

構造後的矩陣比原先大了 81/64 倍的大小。

不過在之後的計算中並沒有額外的消耗。

再次，雖然我們構造出了4個矩陣ABCD，但是實際上沒有對內存做拷貝， ABCD都是累積和矩陣上的映射，沒有對內存中的實際數據做拷貝。

最後，

integral_map = integral_map.astype(np.uint32)

我指定了矩陣中的元素的大小，統一為無符號32位，速度又得到了一定的提升。

但是uint32有其局限，只能存下2^32次方的數據，灰度圖的每一個點的最大取值是2^8，於是uint32的累積矩陣只能存下2^24個像素點，然而由於填充了邊框，矩陣比實際大了81/64，約1.266倍， sqrt(2^24/1.266) = 3640, 所以只能處理3640x3640的圖片，但是從平均情況來說，每一個點平均取2^7，於是我們可以處理最大5148x5148的圖片。

到此，完成了從最初的1600ms到10ms的優化。

彩蛋

混沌邪惡

import cv2 as cvdef wtf(mat, w, h , set_value, threshold, boxsize): # 這才叫Python嘛！ return cv.adaptiveThreshold(mat, 255, cv.ADAPTIVE_THRESH_GAUSSIAN_C, cv.THRESH_BINARY_INV, boxsize, 20)

咕咕咕？？？

我選擇opencv自帶的局部自適應高斯二值化

? python3 implement.pyImage shape 640x480wtf take 8.884578000000001 ms