[OpenCV-Python] OpenCV 中的圖像處理 部分 IV (六)

部分 IV
OpenCV 中的圖像處理

 OpenCV-Python 中文教程（搬運）目錄

23 圖像變換

23.1 傅里葉變換
目標
本小節我們將要學習：
　　• 使用 OpenCV 對圖像進行傅里葉變換
　　• 使用 Numpy 中 FFT（快速傅里葉變換）函數
　　• 傅里葉變換的一些用處
　　• 我們將要學習的函數有：cv2.dft()，cv2.idft() 等
原理
　　傅里葉變換經常被用來分析不同濾波器的頻率特性。我們可以使用 2D 離散傅里葉變換 (DFT) 分析圖像的頻域特性。實現 DFT 的一個快速算法被稱為快速傅里葉變換（FFT）。關於傅里葉變換的細節知識可以在任意一本圖像處理或信號處理的書中找到。請查看本小節中更多資源部分。
對於一個正弦信號：x(t) = Asin(2πft), 它的頻率為 f，如果把這個信號轉到它的頻域表示，我們會在頻率 f 中看到一個峰值。如果我們的信號是由采樣產生的離散信號好組成，我們會得到類似的頻譜圖，只不過前面是連續的，現在是離散。你可以把圖像想象成沿着兩個方向采集的信號。所以對圖像同時進行 X 方向和 Y 方向的傅里葉變換，我們就會得到這幅圖像的頻域表示（頻譜圖）。
更直觀一點，對於一個正弦信號，如果它的幅度變化非常快，我們可以說他是高頻信號，如果變化非常慢，我們稱之為低頻信號。你可以把這種想法應用到圖像中，圖像那里的幅度變化非常大呢？邊界點或者噪聲。所以我們說邊界和噪聲是圖像中的高頻分量（注意這里的高頻是指變化非常快，而非出現的次數多）。如果沒有如此大的幅度變化我們稱之為低頻分量。
現在我們看看怎樣進行傅里葉變換。

23.1.1 Numpy 中的傅里葉變換
　　首先我們看看如何使用 Numpy 進行傅里葉變換。Numpy 中的 FFT 包可以幫助我們實現快速傅里葉變換。函數 np.fft.fft2() 可以對信號進行頻率轉換，輸出結果是一個復雜的數組。本函數的第一個參數是輸入圖像，要求是灰度格式。第二個參數是可選的, 決定輸出數組的大小。輸出數組的大小和輸入圖像大小一樣。如果輸出結果比輸入圖像大，輸入圖像就需要在進行 FFT 前補0。如果輸出結果比輸入圖像小的話，輸入圖像就會被切割。

現在我們得到了結果，頻率為 0 的部分（直流分量）在輸出圖像的左上角。如果想讓它（直流分量）在輸出圖像的中心，我們還需要將結果沿兩個方向平移
N
2
。函數 np.fft.fftshift() 可以幫助我們實現這一步。（這樣更容易分析）。
進行完頻率變換之后，我們就可以構建振幅譜了。

import cv2
import numpy as np
from matplotlib import pyplot as plt
img = cv2.imread('messi5.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# 這里構建振幅圖的公式沒學過
magnitude_spectrum = 20*np.log(np.abs(fshift))
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

結果如下：

　　　　
我們可以看到輸出結果的中心部分更白（亮），這說明低頻分量更多。
現在我們可以進行頻域變換了，我們就可以在頻域對圖像進行一些操作了，例如高通濾波和重建圖像（DFT 的逆變換）。比如我們可以使用一個60x60 的矩形窗口對圖像進行掩模操作從而去除低頻分量。然后再使用函數np.fft.ifftshift() 進行逆平移操作，所以現在直流分量又回到左上角了，左后使用函數 np.ifft2() 進行 FFT 逆變換。同樣又得到一堆復雜的數字，我們可以對他們取絕對值：

rows, cols = img.shape
crow,ccol = rows/2 , cols/2
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])

plt.show()

結果如下：

　　　　

上圖的結果顯示高通濾波其實是一種邊界檢測操作。這就是我們在前面圖像梯度那一章看到的。同時我們還發現圖像中的大部分數據集中在頻譜圖的低頻區域。我們現在已經知道如何使用 Numpy 進行 DFT 和 IDFT 了，接着我們來看看如何使用 OpenCV 進行這些操作。
如果你觀察仔細的話，尤其是最后一章 JET 顏色的圖像，你會看到一些不自然的東西（如我用紅色箭頭標出的區域）。看上圖那里有些條帶裝的結構，這被成為振鈴效應。這是由於我們使用矩形窗口做掩模造成的。這個掩模被轉換成正弦形狀時就會出現這個問題。所以一般我們不適用矩形窗口濾波。最好的選擇是高斯窗口。

23.1.2 OpenCV 中的傅里葉變換
　　OpenCV 中相應的函數是 cv2.dft() 和 cv2.idft()。和前面輸出的結果一樣，但是是雙通道的。第一個通道是結果的實數部分，第二個通道是結果的虛數部分。輸入圖像要首先轉換成 np.float32 格式。我們來看看如何操作。

import numpy as np
import cv2
from matplotlib import pyplot as plt

img = cv2.imread('messi5.jpg',0)

dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

注意：你可以使用函數 cv2.cartToPolar()，它會同時返回幅度和相位。
現在我們來做逆 DFT。在前面的部分我們實現了一個 HPF（高通濾波），現在我們來做 LPF（低通濾波）將高頻部分去除。其實就是對圖像進行模糊操作。首先我們需要構建一個掩模，與低頻區域對應的地方設置為 1, 與高頻區域對應的地方設置為 0。

rows, cols = img.shape
crow,ccol = rows/2 , cols/2

# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1

# apply mask and inverse DFT
fshift = dft_shift*mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

結果如下：

　　　　

注意：OpenCV 中的函數 cv2.dft() 和 cv2.idft() 要比 Numpy 快。但是Numpy 函數更加用戶友好。關於性能的描述，請看下面的章節。

23.1.3 DFT 的性能優化
　　當數組的大小為某些值時 DFT 的性能會更好。當數組的大小是 2 的指數時 DFT 效率最高。當數組的大小是 2，3，5 的倍數時效率也會很高。所以如果你想提高代碼的運行效率時，你可以修改輸入圖像的大小（補 0）。對於OpenCV 你必須自己手動補 0。但是 Numpy，你只需要指定 FFT 運算的大小，它會自動補 0。
那我們怎樣確定最佳大小呢？OpenCV 提供了一個函數:cv2.getOptimalDFTSize()。
它可以同時被 cv2.dft() 和 np.fft.fft2() 使用。讓我們一起使用 IPython的魔法命令%timeit 來測試一下吧。

In [16]: img = cv2.imread('messi5.jpg',0)
In [17]: rows,cols = img.shape
In [18]: print rows,cols
342 548

In [19]: nrows = cv2.getOptimalDFTSize(rows)
In [20]: ncols = cv2.getOptimalDFTSize(cols)
In [21]: print nrows, ncols
360 576

看到了吧，數組的大小從（342，548）變成了（360，576）。現在我們為它補 0，然后看看性能有沒有提升。你可以創建一個大的 0 數組，然后把我們的數據拷貝過去，或者使用函數 cv2.copyMakeBoder()。

nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img

或者：

right = ncols - cols
bottom = nrows - rows
#just to avoid line breakup in PDF file
bordertype = cv2.BORDER_CONSTANT
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)

現在我們看看 Numpy 的表現：

In [22]: %timeit fft1 = np.fft.fft2(img)
10 loops, best of 3: 40.9 ms per loop
In [23]: %timeit fft2 = np.fft.fft2(img,[nrows,ncols])
100 loops, best of 3: 10.4 ms per loop

速度提高了 4 倍。我們再看看 OpenCV 的表現：

In [24]: %timeit dft1= cv2.dft(np.float32(img),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 13.5 ms per loop
In [27]: %timeit dft2= cv2.dft(np.float32(nimg),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 3.11 ms per loop

也提高了 4 倍，同時我們也會發現 OpenCV 的速度是 Numpy 的 3 倍。
你也可以測試一下逆 FFT 的表現。

23.1.4 為什么拉普拉斯算子是高通濾波器？
　　我在論壇中遇到了一個類似的問題。為什么拉普拉斯算子是高通濾波器？
　　為什么 Sobel 是 HPF？等等。對於第一個問題的答案我們以傅里葉變換的形式給出。我們一起來對不同的算子進行傅里葉變換並分析它們：

import cv2
import numpy as np
from matplotlib import pyplot as plt
# simple averaging filter without scaling parameter
mean_filter = np.ones((3,3))
# creating a guassian filter
x = cv2.getGaussianKernel(5,10)
#x.T 為矩陣轉置
gaussian = x*x.T
# different edge detecting filters
# scharr in x-direction
scharr = np.array([[-3, 0, 3],
[-10,0,10],
[-3, 0, 3]])
# sobel in x direction
sobel_x= np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]])
# sobel in y direction
sobel_y= np.array([[-1,-2,-1],
[0, 0, 0],
[1, 2, 1]])
# laplacian
laplacian=np.array([[0, 1, 0],
[1,-4, 1],
[0, 1, 0]])
filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
filter_name = ['mean_filter', 'gaussian','laplacian', 'sobel_x', \
'sobel_y', 'scharr_x']
fft_filters = [np.fft.fft2(x) for x in filters]
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
mag_spectrum = [np.log(np.abs(z)+1) for z in fft_shift]
for i in xrange(6):
plt.subplot(2,3,i+1),plt.imshow(mag_spectrum[i],cmap = 'gray')
plt.title(filter_name[i]), plt.xticks([]), plt.yticks([])
plt.show()

結果：

　　　　
從圖像中我們就可以看出每一個算子允許通過那些信號。從這些信息中我
們就可以知道那些是 HPF 那是 LPF。
更多資源
1. An Intuitive Explanation of Fourier Theoryby Steven Lehar
2. Fourier Transformat HIPR
3. What does frequency domain denote in case of images?
練習

24 模板匹配

目標
在本節我們要學習：
　　1. 使用模板匹配在一幅圖像中查找目標
　　2. 函數：cv2.matchTemplate()，cv2.minMaxLoc()
原理
　　模板匹配是用來在一副大圖中搜尋查找模版圖像位置的方法。OpenCV 為我們提供了函數：cv2.matchTemplate()。和 2D 卷積一樣，它也是用模板圖像在輸入圖像（大圖）上滑動，並在每一個位置對模板圖像和與其對應的輸入圖像的子區域進行比較。OpenCV 提供了幾種不同的比較方法（細節請看文檔）。返回的結果是一個灰度圖像，每一個像素值表示了此區域與模板的匹配程度。
如果輸入圖像的大小是（WxH），模板的大小是（wxh），輸出的結果的大小就是（W-w+1，H-h+1）。當你得到這幅圖之后，就可以使用函數cv2.minMaxLoc() 來找到其中的最小值和最大值的位置了。第一個值為矩形左上角的點（位置），（w，h）為 moban 模板矩形的寬和高。這個矩形就是找到的模板區域了。
注意：如果你使用的比較方法是 cv2.TM_SQDIFF，最小值對應的位置才是匹配的區域。

24.1 OpenCV 中的模板匹配
　　我們這里有一個例子：我們在梅西的照片中搜索梅西的面部。所以我們要制作下面這樣一個模板：

　　　　
我們會嘗試使用不同的比較方法，這樣我們就可以比較一下它們的效果了。

import cv2
import numpy as np
from matplotlib import pyplot as plt
img = cv2.imread('messi5.jpg',0)
img2 = img.copy()
template = cv2.imread('messi_face.jpg',0)
w, h = template.shape[::-1]
# All the 6 methods for comparison in a list
methods = ['cv2.TM_CCOEFF', 'cv2.TM_CCOEFF_NORMED', 'cv2.TM_CCORR',
'cv2.TM_CCORR_NORMED', 'cv2.TM_SQDIFF', 'cv2.TM_SQDIFF_NORMED']
for meth in methods:
img = img2.copy()
#exec 語句用來執行儲存在字符串或文件中的 Python 語句。
# 例如，我們可以在運行時生成一個包含 Python 代碼的字符串，然后使用 exec 語句執行這些語句。
#eval 語句用來計算存儲在字符串中的有效 Python 表達式
method = eval(meth)
# Apply template Matching
res = cv2.matchTemplate(img,template,method)
min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(res)
# 使用不同的比較方法，對結果的解釋不同
# If the method is TM_SQDIFF or TM_SQDIFF_NORMED, take minimum
if method in [cv2.TM_SQDIFF, cv2.TM_SQDIFF_NORMED]:
top_left = min_loc
else:
top_left = max_loc
bottom_right = (top_left[0] + w, top_left[1] + h)
cv2.rectangle(img,top_left, bottom_right, 255, 2)
plt.subplot(121),plt.imshow(res,cmap = 'gray')
plt.title('Matching Result'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img,cmap = 'gray')
plt.title('Detected Point'), plt.xticks([]), plt.yticks([])
plt.suptitle(meth)
plt.show()

結果如下：
cv2.TM_CCOEFF
　　　　
cv2.TM_CCOEFF_NORMED
　　　　
cv2.TM_CCORR
　　　　
cv2.TM_CCORR_NORMED
　　　　
cv2.TM_SQDIFF
　　　　
cv2.TM_SQDIFF_NORMED
　　　　

我們看到 cv2.TM_CCORR 的效果不想我們想的那么好。

24.2 多對象的模板匹配
　　在前面的部分，我們在圖片中搜素梅西的臉，而且梅西只在圖片中出現了一次。假如你的目標對象只在圖像中出現了很多次怎么辦呢？函數cv.imMaxLoc() 只會給出最大值和最小值。此時，我們就要使用閾值了。
在下面的例子中我們要經典游戲 Mario 的一張截屏圖片中找到其中的硬幣。

import cv2
import numpy as np
from matplotlib import pyplot as plt

img_rgb = cv2.imread('mario.png')
img_gray = cv2.cvtColor(img_rgb, cv2.COLOR_BGR2GRAY)
template = cv2.imread('mario_coin.png',0)
w, h = template.shape[::-1]

res = cv2.matchTemplate(img_gray,template,cv2.TM_CCOEFF_NORMED)
threshold = 0.8
loc = np.where( res >= threshold)
for pt in zip(*loc[::-1]):
    cv2.rectangle(img_rgb, pt, (pt[0] + w, pt[1] + h), (0,0,255), 2)

cv2.imwrite('res.png',img_rgb)

結果：
　　　　

25 Hough 直線變換

目標
• 理解霍夫變換的概念
• 學習如何在一張圖片中檢測直線
• 學習函數：cv2.HoughLines()，cv2.HoughLinesP()
原理
　　霍夫變換在檢測各種形狀的的技術中非常流行，如果你要檢測的形狀可以用數學表達式寫出，你就可以是使用霍夫變換檢測它。及時要檢測的形狀存在一點破壞或者扭曲也可以使用。我們下面就看看如何使用霍夫變換檢測直線。一條直線可以用數學表達式 y = mx + c 或者 ρ = xcosθ + y sinθ 表示。
ρ 是從原點到直線的垂直距離，θ 是直線的垂線與橫軸順時針方向的夾角（如果你使用的坐標系不同方向也可能不同，我是按 OpenCV 使用的坐標系描述的）。如下圖所示：

　　　　

所以如果一條線在原點下方經過，ρ 的值就應該大於 0，角度小於 180。但是如果從原點上方經過的話，角度不是大於 180，也是小於 180，但 ρ 的值小於 0。垂直的線角度為 0 度，水平線的角度為 90 度。讓我們來看看霍夫變換是如何工作的。每一條直線都可以用 (ρ,θ) 表示。
所以首先創建一個 2D 數組（累加器），初始化累加器，所有的值都為 0。行表示 ρ，列表示 θ。這個數組的大小決定了最后結果的准確性。如果你希望角度精確到 1 度，你就需要 180 列。對於 ρ，最大值為圖片對角線的距離。所以如果精確度要達到一個像素的級別，行數就應該與圖像對角線的距離相等。
想象一下我們有一個大小為 100x100 的直線位於圖像的中央。取直線上的第一個點，我們知道此處的（x，y）值。把 x 和 y 帶入上邊的方程組，然后遍歷 θ 的取值：0，1，2，3，...，180。分別求出與其對應的 ρ 的值，這樣我們就得到一系列（ρ,θ）的數值對，如果這個數值對在累加器中也存在相應的位置，就在這個位置上加 1。所以現在累加器中的（50，90）=1。（一個點可能存在與多條直線中，所以對於直線上的每一個點可能是累加器中的多個值同時加 1）。
現在取直線上的第二個點。重復上邊的過程。更新累加器中的值。現在累加器中（50,90）的值為 2。你每次做的就是更新累加器中的值。對直線上的每個點都執行上邊的操作，每次操作完成之后，累加器中的值就加 1，但其他地方有時會加 1, 有時不會。按照這種方式下去，到最后累加器中（50,90）的值肯定是最大的。如果你搜索累加器中的最大值，並找到其位置（50,90），這就說明圖像中有一條直線，這條直線到原點的距離為 50，它的垂線與橫軸的夾角為 90 度。下面的動畫很好的演示了這個過程（Image Courtesy: AmosStorkey ）。

　　　　
這就是霍夫直線變換工作的方式。很簡單，也許你自己就可以使用 Numpy搞定它。下圖顯示了一個累加器。其中最亮的兩個點代表了圖像中兩條直線的參數。（Image courtesy: Wikipedia）。

　　　　

25.1 OpenCV 中的霍夫變換
　　上面介紹的整個過程在 OpenCV 中都被封裝進了一個函數：cv2.HoughLines()。
返回值就是（ρ,θ）。ρ 的單位是像素，θ 的單位是弧度。這個函數的第一個參數是一個二值化圖像，所以在進行霍夫變換之前要首先進行二值化，或者進行Canny 邊緣檢測。第二和第三個值分別代表 ρ 和 θ 的精確度。第四個參數是閾值，只有累加其中的值高於閾值時才被認為是一條直線，也可以把它看成能檢測到的直線的最短長度（以像素點為單位）。

import cv2
import numpy as np

img = cv2.imread('dave.jpg')
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
edges = cv2.Canny(gray,50,150,apertureSize = 3)

lines = cv2.HoughLines(edges,1,np.pi/180,200)
for rho,theta in lines[0]:
    a = np.cos(theta)
    b = np.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))

    cv2.line(img,(x1,y1),(x2,y2),(0,0,255),2)

cv2.imwrite('houghlines3.jpg',img)

結果如下：

　　　　

25.2 Probabilistic Hough Transform
　　從上邊的過程我們可以發現：僅僅是一條直線都需要兩個參數，這需要大量的計算。Probabilistic_Hough_Transform 是對霍夫變換的一種優化。它不會對每一個點都進行計算，而是從一幅圖像中隨機選取（是不是也可以使用圖像金字塔呢？）一個點集進行計算，對於直線檢測來說這已經足夠了。但是使用這種變換我們必須要降低閾值（總的點數都少了，閾值肯定也要小呀！）。
下圖是對兩種方法的對比。（Image Courtesy : Franck Bettinger’s homepage）

　　　　
OpenCV 中使用的 Matas, J. ，Galambos, C. 和 Kittler, J.V. 提出的Progressive Probabilistic Hough Transform。這個函數是 cv2.HoughLinesP()。
它有兩個參數。
　　• minLineLength - 線的最短長度。比這個短的線都會被忽略。
　　• MaxLineGap - 兩條線段之間的最大間隔，如果小於此值，這兩條直線就被看成是一條直線。
更加給力的是，這個函數的返回值就是直線的起點和終點。而在前面的例子中，我們只得到了直線的參數，而且你必須要找到所有的直線。而在這里一切都很直接很簡單。

import cv2
import numpy as np

img = cv2.imread('dave.jpg')
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
edges = cv2.Canny(gray,50,150,apertureSize = 3)
minLineLength = 100
maxLineGap = 10
lines = cv2.HoughLinesP(edges,1,np.pi/180,100,minLineLength,maxLineGap)
for x1,y1,x2,y2 in lines[0]:
    cv2.line(img,(x1,y1),(x2,y2),(0,255,0),2)

cv2.imwrite('houghlines5.jpg',img)

結果如下：
　　　　

 

26 Hough 圓環變換

目標
　　• 學習使用霍夫變換在圖像中找圓形（環）。
　　• 學習函數：cv2.HoughCircles()。
原理
　　圓形的數學表達式為 (x − x center ) 2 +(y − y center ) 2 = r 2 ，其中（x center ,y center ）為圓心的坐標，r 為圓的直徑。從這個等式中我們可以看出：一個圓環需要 3個參數來確定。所以進行圓環霍夫變換的累加器必須是 3 維的，這樣的話效率就會很低。所以 OpenCV 用來一個比較巧妙的辦法，霍夫梯度法，它可以使用邊界的梯度信息。
我們要使用的函數為 cv2.HoughCircles()。文檔中對它的參數有詳細的解釋。這里我們就直接看代碼吧。

import cv2
import numpy as np
img = cv2.imread('opencv_logo.png',0)
img = cv2.medianBlur(img,5)
cimg = cv2.cvtColor(img,cv2.COLOR_GRAY2BGR)
circles = cv2.HoughCircles(img,cv2.HOUGH_GRADIENT,1,20,
param1=50,param2=30,minRadius=0,maxRadius=0)
circles = np.uint16(np.around(circles))
for i in circles[0,:]:
# draw the outer circle
cv2.circle(cimg,(i[0],i[1]),i[2],(0,255,0),2)
# draw the center of the circle
cv2.circle(cimg,(i[0],i[1]),2,(0,0,255),3)
cv2.imshow('detected circles',cimg)
cv2.waitKey(0)
cv2.destroyAllWindows()
#Python: cv2.HoughCircles(image, method, dp, minDist, circles, param1, param2, minRadius, maxRadius)
#Parameters:
#image – 8-bit, single-channel, grayscale input image.
# 返回結果為 Output vector of found circles. Each vector is encoded as a
#3-element floating-point vector (x, y, radius) .
#circle_storage – In C function this is a memory storage that will contain
#the output sequence of found circles.
#method – Detection method to use. Currently, the only implemented method is
#CV_HOUGH_GRADIENT , which is basically 21HT , described in [Yuen90].
#dp – Inverse ratio of the accumulator resolution to the image resolution.
#For example, if dp=1 , the accumulator has the same resolution as the input image.
#If dp=2 , the accumulator has half as big width and height.
#minDist – Minimum distance between the centers of the detected circles.
#If the parameter is too small, multiple neighbor circles may be falsely
#detected in addition to a true one. If it is too large, some circles may be missed.
#param1 – First method-specific parameter. In case of CV_HOUGH_GRADIENT ,
#it is the higher threshold of the two passed to the Canny() edge detector
# (the lower one is twice smaller).
#param2 – Second method-specific parameter. In case of CV_HOUGH_GRADIENT ,
# it is the accumulator threshold for the circle centers at the detection stage.
#The smaller it is, the more false circles may be detected. Circles,
# corresponding to the larger accumulator values, will be returned first.
#minRadius – Minimum circle radius.
#maxRadius – Maximum circle radius.

結果如下：

　　　　

更多資源
練習

27 分水嶺算法圖像分割

目標
本節我們將要學習
　　• 使用分水嶺算法基於掩模的圖像分割
　　• 函數：cv2.watershed()
原理
　　任何一副灰度圖像都可以被看成拓撲平面，灰度值高的區域可以被看成是山峰，灰度值低的區域可以被看成是山谷。我們向每一個山谷中灌不同顏色的水。隨着水的位的升高，不同山谷的水就會相遇匯合，為了防止不同山谷的水匯合，我們需要在水匯合的地方構建起堤壩。不停的灌水，不停的構建堤壩知道所有的山峰都被水淹沒。我們構建好的堤壩就是對圖像的分割。這就是分水嶺算法的背后哲理。你可以通過訪問網站CMM webpage on watershed來加深自己的理解。
但是這種方法通常都會得到過度分割的結果，這是由噪聲或者圖像中其他不規律的因素造成的。為了減少這種影響，OpenCV 采用了基於掩模的分水嶺算法，在這種算法中我們要設置那些山谷點會匯合，那些不會。這是一種交互式的圖像分割。我們要做的就是給我們已知的對象打上不同的標簽。如果某個區域肯定是前景或對象，就使用某個顏色（或灰度值）標簽標記它。如果某個區域肯定不是對象而是背景就使用另外一個顏色標簽標記。而剩下的不能確定是前景還是背景的區域就用 0 標記。這就是我們的標簽。然后實施分水嶺算法。每一次灌水，我們的標簽就會被更新，當兩個不同顏色的標簽相遇時就構建堤壩，直到將所有山峰淹沒，最后我們得到的邊界對象（堤壩）的值為 -1。

27.1 代碼
下面的例子中我們將就和距離變換和分水嶺算法對緊挨在一起的對象進行分割。
如下圖所示，這些硬幣緊挨在一起。就算你使用閾值操作，它們任然是緊挨着的。

　　　　

我們從找到硬幣的近似估計開始。我們可以使用 Otsu's 二值化。

import numpy as np
import cv2
from matplotlib import pyplot as plt

img = cv2.imread('coins.png')
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
ret, thresh = cv2.threshold(gray,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)

結果

　　　　
現在我們要去除圖像中的所有的白噪聲。這就需要使用形態學中的開運算。
為了去除對象上小的空洞我們需要使用形態學閉運算。所以我們現在知道靠近對象中心的區域肯定是前景，而遠離對象中心的區域肯定是背景。而不能確定的區域就是硬幣之間的邊界。
所以我們要提取肯定是硬幣的區域。腐蝕操作可以去除邊緣像素。剩下就可以肯定是硬幣了。當硬幣之間沒有接觸時，這種操作是有效的。但是由於硬幣之間是相互接觸的，我們就有了另外一個更好的選擇：距離變換再加上合適的閾值。接下來我們要找到肯定不是硬幣的區域。這是就需要進行膨脹操作了。膨脹可以將對象的邊界延伸到背景中去。這樣由於邊界區域被去處理，我們就可以知道那些區域肯定是前景，那些肯定是背景。如下圖所示。

　　　　
剩下的區域就是我們不知道該如何區分的了。這就是分水嶺算法要做的。
這些區域通常是前景與背景的交界處（或者兩個前景的交界）。我們稱之為邊界。從肯定是不是背景的區域中減去肯定是前景的區域就得到了邊界區域。

# noise removal
kernel = np.ones((3,3),np.uint8)
opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN,kernel, iterations = 2)
# sure background area
sure_bg = cv2.dilate(opening,kernel,iterations=3)
# Finding sure foreground area
# 距離變換的基本含義是計算一個圖像中非零像素點到最近的零像素點的距離，也就是到零像素點的最短距離
# 個最常見的距離變換算法就是通過連續的腐蝕操作來實現，腐蝕操作的停止條件是所有前景像素都被完全
# 腐蝕。這樣根據腐蝕的先后順序，我們就得到各個前景像素點到前景中心??像素點的
# 距離。根據各個像素點的距離值，設置為不同的灰度值。這樣就完成了二值圖像的距離變換
#cv2.distanceTransform(src, distanceType, maskSize)
# 第二個參數 0,1,2 分別表示 CV_DIST_L1, CV_DIST_L2 , CV_DIST_C
dist_transform = cv2.distanceTransform(opening,1,5)
ret, sure_fg = cv2.threshold(dist_transform,0.7*dist_transform.max(),255,0)
# Finding unknown region
sure_fg = np.uint8(sure_fg)
unknown = cv2.subtract(sure_bg,sure_fg)

如結果所示，在閾值化之后的圖像中，我們得到了肯定是硬幣的區域，而且硬幣之間也被分割開了。（有些情況下你可能只需要對前景進行分割，而不需要將緊挨在一起的對象分開，此時就沒有必要使用距離變換了，腐蝕就足夠了。當然腐蝕也可以用來提取肯定是前景的區域。）

　　　　
現在知道了那些是背景那些是硬幣了。那我們就可以創建標簽（一個與原圖像大小相同，數據類型為 in32 的數組），並標記其中的區域了。對我們已經確定分類的區域（無論是前景還是背景）使用不同的正整數標記，對我們不確定的區域使用 0 標記。我們可以使用函數 cv2.connectedComponents()來做這件事。它會把將背景標記為 0，其他的對象使用從 1 開始的正整數標記。
但是，我們知道如果背景標記為 0，那分水嶺算法就會把它當成未知區域了。所以我們想使用不同的整數標記它們。而對不確定的區域（函數cv2.connectedComponents 輸出的結果中使用 unknown 定義未知區域）標記為 0。

# Marker labelling
ret, markers1 = cv2.connectedComponents(sure_fg)
# Add one to all labels so that sure background is not 0, but 1
markers = markers1+1
# Now, mark the region of unknown with zero
markers[unknown==255] = 0
結果使用 JET 顏色地圖表示。深藍色區域為未知區域。肯定是硬幣的區域使用不同的顏色標記。其余區域就是用淺藍色標記的背景了。
現在標簽准備好了。到最后一步：實施分水嶺算法了。標簽圖像將會被修改，邊界區域的標記將變為 -1.
markers3 = cv2.watershed(img,markers)
img[markers3 == -1] = [255,0,0]

結果如下。有些硬幣的邊界被分割的很好，也有一些硬幣之間的邊界分割的不好。
　　　　

Now our marker is ready. It is time for final step, apply watershed. Then marker image will be modified. The boundary region will be marked with -1.

markers = cv2.watershed(img,markers)
img[markers == -1] = [255,0,0]

 

See the result below. For some coins, the region where they touch are segmented properly and for some, they are not.

練習
1. OpenCV 自帶的示例中有一個交互式分水嶺分割程序：watershed.py。
自己玩玩吧。

28 使用 GrabCut 算法進行交互式前景提取

目標
在本節中我們將要學習：
　　• GrabCut 算法原理，使用 GrabCut 算法提取圖像的前景
　　• 創建一個交互是程序完成前景提取

原理
　　GrabCut 算法是由微軟劍橋研究院的 Carsten_Rother，Vladimir_Kolmogorov和 Andrew_Blake 在文章《GrabCut”: interactive foreground extraction using iterated graph cuts》中共同提出的。此算法在提取前景的操作過程中需要很少的人機交互，結果非常好。
　　從用戶的角度來看它到底是如何工作的呢？開始時用戶需要用一個矩形將前景區域框住（前景區域應該完全被包括在矩形框內部）。然后算法進行迭代式分割直達達到最好結果。但是有時分割的結果不夠理想，比如把前景當成了背景，或者把背景當成了前景。在這種情況下，就需要用戶來進行修改了。用戶只需要在不理想的部位畫一筆（點一下鼠標）就可以了。畫一筆就等於在告訴計算機：“嗨，老兄，你把這里弄反了，下次迭代的時候記得改過來呀！”。然后，在下一輪迭代的時候你就會得到一個更好的結果了。
如下圖所示。運動員和足球被藍色矩形包圍在一起。其中有我做的幾個修正，白色畫筆表明這里是前景，黑色畫筆表明這里是背景。最后我得到了一個很好的結果。

　　　　
在整個過程中到底發生了什么呢？
　　• 用戶輸入一個矩形。矩形外的所有區域肯定都是背景（我們在前面已經提到，所有的對象都要包含在矩形框內）。矩形框內的東西是未知的。同樣用戶確定前景和背景的任何操作都不會被程序改變。
　　• 計算機會對我們的輸入圖像做一個初始化標記。它會標記前景和背景像素。
　　• 使用一個高斯混合模型（GMM）對前景和背景建模。
　　• 根據我們的輸入，GMM 會學習並創建新的像素分布。對那些分類未知的像素（可能是前景也可能是背景），可以根據它們與已知分類（如背景）的像素的關系來進行分類（就像是在做聚類操作）。
　　• 這樣就會根據像素的分布創建一副圖。圖中的節點就是像素點。除了像素點做節點之外還有兩個節點：Source_node 和 Sink_node。所有的前景像素都和 Source_node 相連。所有的背景像素都和 Sink_node 相連。
　　• 將像素連接到 Source_node/end_node 的（邊）的權重由它們屬於同一類（同是前景或同是背景）的概率來決定。兩個像素之間的權重由邊的信息或者兩個像素的相似性來決定。如果兩個像素的顏色有很大的不同，那么它們之間的邊的權重就會很小。
　　• 使用 mincut 算法對上面得到的圖進行分割。它會根據最低成本方程將圖分為 Source_node 和 Sink_node。成本方程就是被剪掉的所有邊的權重之和。在裁剪之后，所有連接到 Source_node 的像素被認為是前景，所有連接到 Sink_node 的像素被認為是背景。
　　• 繼續這個過程直到分類收斂。
下圖演示了這個過程（Image Courtesy: http://www.cs.ru.ac.za/research/g02m1682/）：

　　　　

28.1 演示
　　現在我們進入 OpenCV 中的 grabcut 算法。OpenCV 提供了函數：

cv2.grabCut()。我們來先看看它的參數：
　　• img - 輸入圖像
　　• mask-掩模圖像，用來確定那些區域是背景，前景，可能是前景/背景等。可以設置為：cv2.GC_BGD,cv2.GC_FGD,cv2.GC_PR_BGD,cv2.GC_PR_FGD，或者直接輸入 0,1,2,3 也行。
　　• rect - 包含前景的矩形，格式為 (x,y,w,h)
　　• bdgModel, fgdModel - 算法內部使用的數組. 你只需要創建兩個大小為 (1,65)，數據類型為 np.float64 的數組。
　　• iterCount - 算法的迭代次數
　　• mode 可以設置為 cv2.GC_INIT_WITH_RECT 或 cv2.GC_INIT_WITH_MASK，

也可以聯合使用。這是用來確定我們進行修改的方式，矩形模式或者掩模模式。
首先，我們來看使用矩形模式。加載圖片，創建掩模圖像，構建 bdgModel和 fgdModel。傳入矩形參數。都是這么直接。讓算法迭代 5 次。由於我們在使用矩形模式所以修改模式設置為 cv2.GC_INIT_WITH_RECT。運行grabcut。算法會修改掩模圖像，在新的掩模圖像中，所有的像素被分為四類：
背景，前景，可能是背景/前景使用 4 個不同的標簽標記（前面參數中提到過）。
然后我們來修改掩模圖像，所有的 0 像素和 1 像素都被歸為 0（例如背景），所有的 1 像素和 3 像素都被歸為 1（前景）。我們最終的掩模圖像就這樣准備好了。用它和輸入圖像相乘就得到了分割好的圖像。

import numpy as np
import cv2
from matplotlib import pyplot as plt
img = cv2.imread('messi5.jpg')
mask = np.zeros(img.shape[:2],np.uint8)
bgdModel = np.zeros((1,65),np.float64)
fgdModel = np.zeros((1,65),np.float64)
rect = (50,50,450,290)
# 函數的返回值是更新的 mask, bgdModel, fgdModel
cv2.grabCut(img,mask,rect,bgdModel,fgdModel,5,cv2.GC_INIT_WITH_RECT)
mask2 = np.where((mask==2)|(mask==0),0,1).astype('uint8')
img = img*mask2[:,:,np.newaxis]
plt.imshow(img),plt.colorbar(),plt.show()

結果如下：

　　　　
　　哎呀，梅西的頭發被我們弄沒了！讓我們來幫他找回頭發。所以我們要在那里畫一筆（設置像素為 1，肯定是前景）。同時還有一些我們並不需要的草地。我們需要去除它們，我們再在這個區域畫一筆（設置像素為 0，肯定是背景）。現在可以象前面提到的那樣來修改掩模圖像了。
實際上我是怎么做的呢？我們使用圖像編輯軟件打開輸入圖像，添加一個圖層，使用筆刷工具在需要的地方使用白色繪制（比如頭發，鞋子，球等）；使用黑色筆刷在不需要的地方繪制（比如，logo，草地等）。然后將其他地方用灰色填充，保存成新的掩碼圖像。在 OpenCV 中導入這個掩模圖像，根據新的掩碼圖像對原來的掩模圖像進行編輯。代碼如下：

# newmask is the mask image I manually labelled
newmask = cv2.imread('newmask.png',0)

# whereever it is marked white (sure foreground), change mask=1
# whereever it is marked black (sure background), change mask=0
mask[newmask == 0] = 0
mask[newmask == 255] = 1

mask, bgdModel, fgdModel = cv2.grabCut(img,mask,None,bgdModel,fgdModel,5,cv2.GC_INIT_WITH_MASK)

mask = np.where((mask==2)|(mask==0),0,1).astype('uint8')
img = img*mask[:,:,np.newaxis]
plt.imshow(img),plt.colorbar(),plt.show()

結果如下：

　　　　
就是這樣。你也可以不使用矩形初始化，直接進入掩碼圖像模式。使用 2像素和 3 像素（可能是背景/前景）對矩形區域的像素進行標記。然后象我們在第二個例子中那樣對肯定是前景的像素標記為 1 像素。然后直接在掩模圖像模式下使用 grabCut 函數。

練習
1. OpenCV 自帶的示例中有一個使用 grabcut 算法的交互式工具 grabcut.py，自己玩一下吧。
2. 你可以創建一個交互是取樣程序，可以繪制矩形，帶有滑動條（調節筆刷的粗細）等。