前兩天參加了北師的數(shù)學(xué)建模校賽，B題是一道圖像處理的題，于是趁機(jī)練習(xí)了一下OpenCV，現(xiàn)在把做的東西移植過來。

（2020.5.31補(bǔ)充：此方法在競賽中取得二等獎(jiǎng)。這次的參賽論文的確存在一些問題，例如沒有對(duì)結(jié)果進(jìn)行量化評(píng)估、對(duì)處理方式的具體細(xì)節(jié)敘述得不夠明確、參考文獻(xiàn)不夠豐富（好吧，其實(shí)沒有引用參考文獻(xiàn)）等。）

題目大意

給出一張生產(chǎn)線上拍攝的手機(jī)鏡頭的圖像（如下），要求解決三個(gè)問題：

建立模型構(gòu)造出一種分割方法，可以將左右兩個(gè)鏡頭的待測(cè)區(qū)域（白色環(huán)形內(nèi)區(qū)域）準(zhǔn)確地分離出來。

建立模型構(gòu)造一種檢測(cè)方法，自動(dòng)地在待測(cè)區(qū)域之內(nèi)將所有缺陷點(diǎn)找出，缺陷點(diǎn)為人眼可識(shí)別的白點(diǎn)，最小可為一個(gè)像素點(diǎn)。要求給出缺陷點(diǎn)的數(shù)學(xué)描述，并根據(jù)該描述建立檢測(cè)模型，自動(dòng)確定每個(gè)缺陷點(diǎn)的位置和像素大小。給出右側(cè)鏡頭中按像素大小排序的最大的前五個(gè)缺陷點(diǎn)的位置坐標(biāo)。

由于在實(shí)際拍照中鏡頭可能會(huì)在模具中抖動(dòng)，所以拍攝的圖片可能并不是正對(duì)鏡頭的，此為圖像的偏心現(xiàn)象。比如圖中左側(cè)圖像就是正對(duì)的情況，右側(cè)就是不正對(duì)（偏心）的情況。建立模型構(gòu)造一種校正方法，校正右側(cè)圖像的偏心現(xiàn)象。呈現(xiàn)校正效果，并給出第2問所求五個(gè)缺陷點(diǎn)校正后的位置坐標(biāo)。

問題求解

問題一

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

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

預(yù)處理結(jié)束之后直接使用OpenCV內(nèi)置的findContours()尋找邊界，這個(gè)函數(shù)非常方便的一點(diǎn)是它可以根據(jù)輪廓之間的嵌套關(guān)系，對(duì)各個(gè)輪廓構(gòu)造層次關(guān)系，函數(shù)原型為：

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

其中的必要參數(shù)為：

image：輸入的圖像

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

cv2.RETR_EXTERNAL：只檢索外輪廓。

cv2.RETR_LIST：檢索輪廓不建立層次關(guān)系。

cv2.RETR_CCOMP：建立兩層的層次關(guān)系（即父子關(guān)系），父輪廓為外輪廓，子輪廓為相應(yīng)內(nèi)孔。若內(nèi)孔中還有內(nèi)孔，則下一層的內(nèi)孔作為另一個(gè)父輪廓。

cv2.RETR_TREE：對(duì)輪廓建立等級(jí)樹的層次關(guān)系。

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

cv2.CHAIN_APPROX_NONE：儲(chǔ)存所有輪廓點(diǎn)，相鄰兩個(gè)點(diǎn)的橫縱坐標(biāo)之差均不超過1。

cv2.CHAIN_APPROX_SIMPLE：壓縮水平方向，垂直方向，對(duì)角線方向的元素，只保留該方向的終點(diǎn)坐標(biāo)，例如一個(gè)矩形輪廓只需4個(gè)點(diǎn)來保存輪廓信息。

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

對(duì)于圓環(huán)來說直接選取CCOMP模式，圖中的所有輪廓中只有圓環(huán)的外輪廓有子輪廓，從而找出所有的目標(biāo)邊界。根據(jù)此邊界創(chuàng)建遮罩，再將遮罩與原圖做按位與即可分割出目標(biāo)圖像。最后做出來的結(jié)果相當(dāng)不錯(cuò)。（就不放在這里了，有興趣的可以自己拿原圖跑一下代碼）

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()

問題二

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

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個(gè)缺陷點(diǎn)，雖然大多數(shù)都只有1~2個(gè)像素但不得不吐槽這個(gè)鏡頭的加工技術(shù)確實(shí)不太行。

問題三

這個(gè)問題是對(duì)鏡頭在模具內(nèi)抖動(dòng)造成的偏心畸變進(jìn)行修正，再重新計(jì)算缺陷點(diǎn)坐標(biāo)。修正畸變是本次各個(gè)問題中最為棘手的一個(gè)部分。查找了一下資料，偏心畸變是由于圖像中目標(biāo)的光軸與攝像機(jī)的光軸不重合造成的，這也是偏心畸變?cè)谟⑽闹斜环Q為decentering distortion的原因。在本問題中，大概是這樣：

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

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

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

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

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()

總結(jié)

雖然之前處理過一些圖像處理問題但從來沒有像這次一樣完整地做一次題目，也沒有深入地了解過各個(gè)運(yùn)算的內(nèi)在原理。這次的圖像處理問題前兩個(gè)比較基礎(chǔ)，最后一個(gè)比較有挑戰(zhàn)性，感覺對(duì)于學(xué)習(xí)OpenCV還是很有幫助的。