前兩天參加了北師的數學建模校賽，B題是一道圖像處理的題，于是趁機練習了一下OpenCV，現在把做的東西移植過來。

（2020.5.31補充：此方法在競賽中取得二等獎。這次的參賽論文的確存在一些問題，例如沒有對結果進行量化評估、對處理方式的具體細節敘述得不夠明確、參考文獻不夠豐富（好吧，其實沒有引用參考文獻）等。）

題目大意

給出一張生產線上拍攝的手機鏡頭的圖像（如下），要求解決三個問題：

建立模型構造出一種分割方法，可以將左右兩個鏡頭的待測區域（白色環形內區域）準確地分離出來。

建立模型構造一種檢測方法，自動地在待測區域之內將所有缺陷點找出，缺陷點為人眼可識別的白點，最小可為一個像素點。要求給出缺陷點的數學描述，并根據該描述建立檢測模型，自動確定每個缺陷點的位置和像素大小。給出右側鏡頭中按像素大小排序的最大的前五個缺陷點的位置坐標。

由于在實際拍照中鏡頭可能會在模具中抖動，所以拍攝的圖片可能并不是正對鏡頭的，此為圖像的偏心現象。比如圖中左側圖像就是正對的情況，右側就是不正對（偏心）的情況。建立模型構造一種校正方法，校正右側圖像的偏心現象。呈現校正效果，并給出第2問所求五個缺陷點校正后的位置坐標。

問題求解

問題一

這個問題是目標檢測，并且需求十分明確：提取出白色圓環中的區域的圖像。觀察圖像可以發現圖中白色的部分幾乎只有需要檢測的白色圓環，其他的白色區域基本上都是不規則圖形以及一些噪點。一種比較簡單的處理方式是直接選取一個合適的閾值二值化，把除了需要的白色圓環之外的區域全部置位黑色。不過為了魯棒性我們并沒有使用這種簡單粗暴的方式。

我們的預處理方法是二值化去除多余細節開運算去除噪點高斯濾波減小像素間梯度，完成預處理后再進行輪廓搜索。二值化采取了全局二值化，主要是在最大類間方差法（OTSU法）與三角形法兩者之間進行選取，實驗發現后者會使黑白區域邊界模糊且曲折，并且很多白色噪點（第二問要檢測）受到了干擾，因此選擇了前者作為二值化方法。開運算的卷積核為20×20矩形卷積核，進行了一次效果就很好了。高斯濾波的直徑（嚴格來說并不能叫做直徑）經驗性地確定為了5。預處理后的效果如下圖所示。

預處理結束之后直接使用OpenCV內置的findContours()尋找邊界，這個函數非常方便的一點是它可以根據輪廓之間的嵌套關系，對各個輪廓構造層次關系，函數原型為：

cv2.findContours(image, mode, method[, contours[, hierarchy[, offset ]]])

其中的必要參數為：

image：輸入的圖像

mode：輪廓的檢索模式，共有四種取值。

cv2.RETR_EXTERNAL：只檢索外輪廓。

cv2.RETR_LIST：檢索輪廓不建立層次關系。

cv2.RETR_CCOMP：建立兩層的層次關系（即父子關系），父輪廓為外輪廓，子輪廓為相應內孔。若內孔中還有內孔，則下一層的內孔作為另一個父輪廓。

cv2.RETR_TREE：對輪廓建立等級樹的層次關系。

method：輪廓的逼近方法，共有四種取值。

cv2.CHAIN_APPROX_NONE：儲存所有輪廓點，相鄰兩個點的橫縱坐標之差均不超過1。

cv2.CHAIN_APPROX_SIMPLE：壓縮水平方向，垂直方向，對角線方向的元素，只保留該方向的終點坐標，例如一個矩形輪廓只需4個點來保存輪廓信息。

cv2.CHAIN_APPROX_TC89_L1和cv2.CHAIN_APPROX_TC89_KCOS都是使用teh-Chinl chain近似算法。

對于圓環來說直接選取CCOMP模式，圖中的所有輪廓中只有圓環的外輪廓有子輪廓，從而找出所有的目標邊界。根據此邊界創建遮罩，再將遮罩與原圖做按位與即可分割出目標圖像。最后做出來的結果相當不錯。（就不放在這里了，有興趣的可以自己拿原圖跑一下代碼）

1 ''' 2 * detection.py 3 Runtime environment: 4 python = 3.7.4 5 opencv-python = 4.1.1.26 6 numpy = 1.17.2 7 ''' 8 9 from cv2 import imread, IMREAD_GRAYSCALE, threshold, THRESH_BINARY, THRESH_OTSU,10 getStructuringElement, MORPH_RECT, erode, dilate, GaussianBlur, findContours,11 RETR_CCOMP, CHAIN_APPROX_SIMPLE, IMREAD_COLOR, drawContours, bitwise_and,12 imwrite13 from numpy import zeros, shape, uint814 15 def detection():16 original = imread('Original.png', IMREAD_GRAYSCALE)17 _, binary = threshold(original, 0, 255, THRESH_BINARY | THRESH_OTSU)18 imwrite('problem1_binary.png', binary)19 kernel = getStructuringElement(MORPH_RECT, (20, 20))20 eroded = erode(binary, kernel)21 dilated = dilate(eroded, kernel)22 blur = GaussianBlur(dilated, (5, 5), 0)23 imwrite('problem1_preprocess.png', blur)24 contours, hierarchies = findContours(blur, RETR_CCOMP, CHAIN_APPROX_SIMPLE)25 chromatic = imread('Original.png', IMREAD_COLOR)26 drawContours(chromatic, contours, -1, (0, 0, 255), 10)27 imwrite('problem1_contours.png', chromatic)28 chromatic = imread('Original.png', IMREAD_COLOR)29 for hierarchy in hierarchies[0, :]:30 if hierarchy[2] != -1:31 drawContours(chromatic, contours, hierarchy[2], (255, 0, 255), 15)32 imwrite('problem1_target_contours.png', chromatic)33 chromatic = imread('Original.png', IMREAD_COLOR)34 count = 035 for hierarchy in hierarchies[0, :]:36 if hierarchy[2] != -1:37 mask = zeros(shape(chromatic), dtype = uint8)38 drawContours(mask, contours, hierarchy[2], (255, 255, 255), -1)39 imwrite('mask' + str(count) + '.png', mask)40 imwrite('detection' + str(count) + '.png', bitwise_and(chromatic, mask))41 count += 142 43 if __name__ == '__main__':44 detection()

問題二

檢測缺陷點還要計算大小，這很明顯是一個圖搜索問題。把問題一預處理第一步，也就是二值化得到的圖像與遮罩疊加，所需要搜索的缺陷點就都顯現出來了。需要做的只是遍歷圖像中所有點，然后對每個點進行廣度優先搜索就可以了。這個問題也比較順利地解決了，唯一的缺點是遍歷廣搜運行起來有一點慢，要運行數十秒才能得到結果。

1 ''' 2 Runtime environment: 3 python = 3.7.4 4 opencv-python = 4.1.1.26 5 numpy = 1.17.2 6 ''' 7 8 from cv2 import imread, IMREAD_GRAYSCALE, threshold, THRESH_BINARY, THRESH_OTSU, 9 imwrite, bitwise_and, IMREAD_COLOR, circle10 from numpy import shape, zeros, uint811 12 def findDefect():13 original = imread('Original.png', IMREAD_GRAYSCALE)14 _, binary = threshold(original, 0, 255, THRESH_BINARY | THRESH_OTSU)15 mask = imread('mask0.png', IMREAD_GRAYSCALE)16 target = bitwise_and(binary, mask)17 imwrite('problem2_target.png', target)18 flag = zeros(shape(target), dtype = uint8)19 defects = []20 for i in range(shape(target)[0]):21 for j in range(shape(target)[1]):22 if target[i][j] == 255 and flag[i][j] == 0:23 queue = []24 head, tail= 0, 025 x, y = i, j26 queue.append(None)27 queue[head] = (x, y)28 flag[x][y] = 129 head += 130 while head > tail:31 if x > 0 and target[x - 1][y] == 255 and flag[x - 1][y] == 0:32 queue.append(None)33 queue[head] = (x - 1, y)34 flag[x - 1][y] = 135 head += 136 if y > 0 and target[x][y - 1] == 255 and flag[x][y - 1] == 0:37 queue.append(None)38 queue[head] = (x, y - 1)39 flag[x][y - 1] = 140 head += 141 if x < shape(target)[0] - 1 and target[x + 1][y] == 255 and flag[x + 1][y] == 0:42 queue.append(None)43 queue[head] = (x + 1, y)44 flag[x + 1][y] = 145 head += 146 if y < shape(target)[1] - 1 and target[x][y + 1] == 255 and flag[x][y + 1] == 0:47 queue.append(None)48 queue[head] = (x, y + 1)49 flag[x][y + 1] = 150 head += 151 (x, y) = queue[tail]52 tail = tail + 153 size = len(queue)54 xsum, ysum = 0, 055 for (x, y) in queue:56 xsum += x57 ysum += y58 defects.append((size, xsum // size, ysum // size))59 defects.sort()60 print(defects[::-1], len(defects))61 print(defects[-5:])62 return defects[-5:]63 64 def visualize(defects):65 original = imread('Original.png', IMREAD_COLOR)66 for defect in defects:67 circle(original, (defect[2], defect[1]), 10, (0, 0, 255), -1)68 imwrite('defects.png', original)69 70 if __name__ == '__main__':71 defects = findDefect()72 visualize(defects)

最后得到了116個缺陷點，雖然大多數都只有1~2個像素但不得不吐槽這個鏡頭的加工技術確實不太行。

問題三

這個問題是對鏡頭在模具內抖動造成的偏心畸變進行修正，再重新計算缺陷點坐標。修正畸變是本次各個問題中最為棘手的一個部分。查找了一下資料，偏心畸變是由于圖像中目標的光軸與攝像機的光軸不重合造成的，這也是偏心畸變在英文中被稱為decentering distortion的原因。在本問題中，大概是這樣：

如果把鏡頭內表面看做圓錐面的話，偏心畸變的產生原因就是這個圓錐稍微“倒下”了一點。想要從幾何上對齊進行修正，就要將這個圓錐“扶正”，具體方式是將圓錐面上的每個點都映射到另一個正立的圓錐上，使得其在母線上的位置比例關系不變。

如上圖，這是從圓錐的底面看向頂點的視圖。目標是將紅色的圓錐母線映射到藍色的圓錐母線上，在左圖看來，就是對于圓O內任意一點P，連接O'P并延長交圓O于Q，連接OQ，在OQ上找到一點P'使得O'P/O'Q=OP'/OQ，P'即為所求。具體的公式推導就不推了，主要過程就是先將原坐標系中的坐標映射到這個圓O的坐標系中，得到目標點的坐標后再映射回去（因為線性代數很多都忘記了所以數學推導十分受苦QAQ）。最后的修正效果如下：

左圖是修正前的原圖，右圖是修正后的圖像。雖然直觀上看并沒有太大變化，但仔細觀察中間的深色原點以及深灰色圓形陰影的位置，就可以看出整幅圖像得到了一個從左下到右上的校正。最后的總體效果還是比較令人滿意的，在新的圖像上重復問題一、二的算法，問題即得解。

問題三的代碼如下，主要有兩部分，第一部分是求中間小的深色圓形陰影位置的代碼，第二部分是進行畸變校正的代碼（實現比較暴力，相應地運行效率也比較低）。

1 ''' 2 * locating.py 3 Runtime environment: 4 python = 3.7.4 5 opencv-python = 4.1.1.26 6 numpy = 1.17.2 7 ''' 8 9 from cv2 import imread, IMREAD_GRAYSCALE, threshold, THRESH_OTSU, THRESH_BINARY,10 imshow, waitKey, imwrite, THRESH_TRIANGLE, adaptiveThreshold, ADAPTIVE_THRESH_MEAN_C,11 ADAPTIVE_THRESH_GAUSSIAN_C, HoughCircles, HOUGH_GRADIENT, circle,12 getStructuringElement, MORPH_RECT, erode, dilate, medianBlur, GaussianBlur,13 Canny, findContours, RETR_CCOMP, CHAIN_APPROX_SIMPLE, drawContours,14 IMREAD_COLOR, RETR_TREE, minEnclosingCircle15 from numpy import uint1616 17 if __name__ == '__main__':18 original = imread('detection0.png', IMREAD_GRAYSCALE) # read original image as grayscale image19 kernel = getStructuringElement(MORPH_RECT, (20, 20))20 eroded = erode(original, kernel)21 dilated = dilate(eroded, kernel)22 dilated = dilate(dilated, kernel)23 eroded = erode(dilated, kernel)24 blur = GaussianBlur(eroded, (5, 5), 0)25 original = blur26 binary = adaptiveThreshold(original, 255, ADAPTIVE_THRESH_MEAN_C, THRESH_BINARY, 99, 2)27 imwrite('adaptive_mean.png', binary) # save image adaptive mean method(loc.)28 origin = imread('adaptive_mean.png', IMREAD_GRAYSCALE)29 kernel = getStructuringElement(MORPH_RECT, (40, 40))30 eroded = erode(origin, kernel)31 dilated = dilate(eroded, kernel)32 blur = GaussianBlur(dilated, (5, 5), 0)33 origin = blur34 contours, hierarchies = findContours(origin, RETR_TREE, CHAIN_APPROX_SIMPLE)35 print(hierarchies)36 chromatic = imread('Original.png', IMREAD_COLOR)37 for i in range(len(hierarchies[0])):38 if hierarchies[0][i][2] == -1:39 break40 length = len(contours[i])41 (x0, y0), r = minEnclosingCircle(contours[i])42 sum = [0, 0]43 for k in contours[i]:44 sum = sum + k45 print(sum // length)46 x, y = tuple(sum[0] // length)47 circle(chromatic, (int(x0), int(y0)), 5, (0, 255, 0), -1)48 circle(chromatic, (int(x0), int(y0)), int(r), (0, 255, 0), 10)49 X, Y, R = (2585, 1270, 433)50 circle(chromatic, (X, Y), 5, (0, 0, 255), -1)51 circle(chromatic, (X, Y), R, (0, 0, 255), 10)52 print(int(x0), int(y0), int(r))53 print(X, Y, R)54 imwrite('contours.png', chromatic)

1 """ 2 * calibrate.py 3 Runtime environment: 4 python = 3.7.4 5 opencv-python = 4.1.1.26 6 numpy = 1.17.2 7 """ 8 9 from math import sqrt10 from cv2 import imread, IMREAD_GRAYSCALE, imwrite, medianBlur11 from numpy import shape12 13 14 def dist(p1, p2):15 r = (float(p1[0] - p2[0]) ** 2 + float(p1[1] - p2[1]) ** 2) ** 0.516 return r17 18 19 def calibrate():20 x, y, r = 2567.0, 1289.0, 63.021 x0, y0, r0 = 2585.0, 1270.0, 433.022 dist0 = dist((x, y), (x0, y0))23 input_img = imread('Original.png', IMREAD_GRAYSCALE)24 output = imread('Original.png', IMREAD_GRAYSCALE)25 tan_theta = float(y - y0) / float(x0 - x)26 sin_theta = tan_theta / sqrt(1 + tan_theta * tan_theta)27 cos_theta = 1 / sqrt(1 + tan_theta * tan_theta)28 sin_theta, cos_theta = sin_theta.real, cos_theta.real29 for i in range(shape(input_img)[1]):30 for j in range(shape(input_img)[0]):31 original = (i, j)32 if dist(original, (x0, y0)) < r0:33 neo = (cos_theta * float(i - x0) - sin_theta * float(j - y0),34 -sin_theta * float(i - x0) - cos_theta * float(j - y0))35 a = float(neo[1]) ** 2 + (float(neo[0]) + dist0) ** 236 b = -2.0 * float(neo[1]) * dist0 * (float(neo[0]) + dist0)37 c = float(neo[1]) ** 2 * (dist0 ** 2 - r0 ** 2)38 delta = b ** 2 - 4 * a * c39 if delta < 0 or a == 0 or float(neo[1]) == 0:40 continue41 yr = (sqrt(delta) - b) / (2 * a)42 if (yr * float(neo[1])) < 0:43 yr = (0 - b - sqrt(delta)) / (2 * a)44 xr = ((float(neo[0]) + dist0) * yr / float(neo[1])) - dist045 x2, y2 = xr / yr * float(neo[1]), float(neo[1])46 real = (cos_theta * x2 - sin_theta * y2 + x0, -sin_theta * x2 - cos_theta * y2 + y0)47 output[int(real[1])][int(real[0])] = input_img[int(original[1])][int(original[0])]48 imwrite('problem3_after_mapping.png', output)49 medianed = medianBlur(output, 3)50 imwrite('Result3.png', medianed)51 52 53 if __name__ == '__main__':54 calibrate()

總結

雖然之前處理過一些圖像處理問題但從來沒有像這次一樣完整地做一次題目，也沒有深入地了解過各個運算的內在原理。這次的圖像處理問題前兩個比較基礎，最后一個比較有挑戰性，感覺對于學習OpenCV還是很有幫助的。