手機(jī)鏡頭目標(biāo)提取、缺陷檢測(cè)與圖像畸變校正

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

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

題目大意

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

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

建立模型構(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問(wèn)所求五個(gè)缺陷點(diǎn)校正后的位置坐標(biāo)。

問(wèn)題求解

問(wèn)題一

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

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

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

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

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

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

對(duì)于圓環(huán)來(lái)說(shuō)直接選取CCOMP模式，圖中的所有輪廓中只有圓環(huán)的外輪廓有子輪廓，從而找出所有的目標(biāo)邊界。根據(jù)此邊界創(chuàng)建遮罩，再將遮罩與原圖做按位與即可分割出目標(biāo)圖像。最后做出來(lái)的結(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()

問(wèn)題二

檢測(cè)缺陷點(diǎn)還要計(jì)算大小，這很明顯是一個(gè)圖搜索問(wèn)題。把問(wèn)題一預(yù)處理第一步，也就是二值化得到的圖像與遮罩疊加，所需要搜索的缺陷點(diǎn)就都顯現(xiàn)出來(lái)了。需要做的只是遍歷圖像中所有點(diǎn)，然后對(duì)每個(gè)點(diǎn)進(jìn)行廣度優(yōu)先搜索就可以了。這個(gè)問(wèn)題也比較順利地解決了，唯一的缺點(diǎn)是遍歷廣搜運(yùn)行起來(lái)有一點(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í)不太行。

問(wèn)題三

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

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

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

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

問(wèn)題三的代碼如下，主要有兩部分，第一部分是求中間小的深色圓形陰影位置的代碼，第二部分是進(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é)

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