python skimage圖像處理(三)

python skimage圖像處理(三)

 

This blog is from: https://www.jianshu.com/p/7693222523c0 

 

霍夫線變換

在圖片處理中，霍夫變換主要是用來檢測圖片中的幾何形狀，包括直線、圓、橢圓等。
在skimage中，霍夫變換是放在tranform模塊內，本篇主要講解霍夫線變換。
對於平面中的一條直線，在笛卡爾座標系中，可用y=mx+b來表示，其中m爲斜率，b爲截距。可是若是直線是一條垂直線，則m爲無窮大，全部一般咱們在另外一座標系中表示直線，即極座標系下的r=xcos(theta)+ysin(theta)。便可用（r,theta）來表示一條直線。其中r爲該直線到原點的距離，theta爲該直線的垂線與x軸的夾角。以下圖所示。

 

 

 

對於一個給定的點（x0,y0), 咱們在極座標下繪出全部經過它的直線（r,theta)，將獲得一條正弦曲線。若是將圖片中的全部非0點的正弦曲線都繪製出來，則會存在一些交點。全部通過這個交點的正弦曲線，說明都擁有一樣的(r,theta), 意味着這些點在一條直線上。

 

 

發上圖所示，三個點(對應圖中的三條正弦曲線）在一條直線上，由於這三個曲線交於一點，具備相同的（r, theta)。霍夫線變換就是利用這種方法來尋找圖中的直線。
函數：

skimage.transform.hough_line(img) 

返回三個值：
h: 霍夫變換累積器
theta: 點與x軸的夾角集合，通常爲0-179度
distance: 點到原點的距離，即上面的所說的r.
例：

import skimage.transform as st import numpy as np import matplotlib.pyplot as plt #matplotlib inline # 構建測試圖片 image = np.zeros((100, 100)) #背景圖 idx = np.arange(25, 75) #25-74序列 image[idx[::-1], idx] = 255 # 線條\ image[idx, idx] = 255 # 線條/ # hough線變換 h, theta, d = st.hough_line(image) #生成一個一行兩列的窗口（可顯示兩張圖片）. fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 6)) plt.tight_layout() #顯示原始圖片 ax0.imshow(image, plt.cm.gray) ax0.set_title('Input image') ax0.set_axis_off() #顯示hough變換所得數據 ax1.imshow(np.log(1 + h)) ax1.set_title('Hough transform') ax1.set_xlabel('Angles (degrees)') ax1.set_ylabel('Distance (pixels)') ax1.axis('image') plt.show() 

 

 

從右邊那張圖能夠看出，有兩個交點，說明原圖像中有兩條直線。
若是咱們要把圖中的兩條直線繪製出來，則須要用到另一個函數：

skimage.transform.hough_line_peaks(hspace, angles, dists） 

用這個函數能夠取出峯值點，即交點，也即原圖中的直線。
返回的參數與輸入的參數同樣。咱們修改一下上邊的程序，在原圖中將兩直線繪製出來。

import skimage.transform as st import numpy as np import matplotlib.pyplot as plt # 構建測試圖片 image = np.zeros((100, 100)) #背景圖 idx = np.arange(25, 75) #25-74序列 image[idx[::-1], idx] = 255 # 線條\ image[idx, idx] = 255 # 線條/ # hough線變換 h, theta, d = st.hough_line(image) #生成一個一行三列的窗口（可顯示三張圖片）. fig, (ax0, ax1,ax2) = plt.subplots(1, 3, figsize=(8, 6)) plt.tight_layout() #顯示原始圖片 ax0.imshow(image, plt.cm.gray) ax0.set_title('Input image') ax0.set_axis_off() #顯示hough變換所得數據 ax1.imshow(np.log(1 + h)) ax1.set_title('Hough transform') ax1.set_xlabel('Angles (degrees)') ax1.set_ylabel('Distance (pixels)') ax1.axis('image') #顯示檢測出的線條 ax2.imshow(image, plt.cm.gray) row1, col1 = image.shape for _, angle, dist in zip(*st.hough_line_peaks(h, theta, d)): y0 = (dist - 0 * np.cos(angle)) / np.sin(angle) y1 = (dist - col1 * np.cos(angle)) / np.sin(angle) ax2.plot((0, col1), (y0, y1), '-r') ax2.axis((0, col1, row1, 0)) ax2.set_title('Detected lines') ax2.set_axis_off() plt.show() 

 

 

注意，繪製線條的時候，要從極座標轉換爲笛卡爾座標，公式爲：

 

 

skimage還提供了另一個檢測直線的霍夫變換函數，
機率霍夫線變換：

skimage.transform.probabilistic_hough_line(img, threshold=10, line_length=5,line_gap=3)

參數：

img: 待檢測的圖像。 threshold： 閾值，可先項，默認爲10 line_length: 檢測的最短線條長度，默認爲50 line_gap: 線條間的最大間隙。增大這個值能夠合併破碎的線條。默認爲10

返回：

lines: 線條列表, 格式如((x0, y0), (x1, y0))，標明開始點和結束點。

下面，咱們用canny算子提取邊緣，而後檢測哪些邊緣是直線？

import skimage.transform as st import matplotlib.pyplot as plt from skimage import data,feature #使用Probabilistic Hough Transform. image = data.camera() edges = feature.canny(image, sigma=2, low_threshold=1, high_threshold=25) lines = st.probabilistic_hough_line(edges, threshold=10, line_length=5,line_gap=3) print(len(lines)) # 建立顯示窗口. fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16, 6)) plt.tight_layout() #顯示原圖像 ax0.imshow(image, plt.cm.gray) ax0.set_title('Input image') ax0.set_axis_off() #顯示canny邊緣 ax1.imshow(edges, plt.cm.gray) ax1.set_title('Canny edges') ax1.set_axis_off() #用plot繪製出全部的直線 ax2.imshow(edges * 0) for line in lines: p0, p1 = line ax2.plot((p0[0], p1[0]), (p0[1], p1[1])) row2, col2 = image.shape ax2.axis((0, col2, row2, 0)) ax2.set_title('Probabilistic Hough') ax2.set_axis_off() plt.show() 

 

---

在極座標中，圓的表示方式爲：

x=x0+rcosθ y=y0+rsinθ

圓心爲(x0,y0),r爲半徑，θ爲旋轉度數，值範圍爲0-359
若是給定圓心點和半徑，則其它點是否在圓上，咱們就能檢測出來了。在圖像中，咱們將每一個非0像素點做爲圓心點，以必定的半徑進行檢測，若是有一個點在圓上，咱們就對這個圓心累加一次。若是檢測到一個圓，那麼這個圓心點就累加到最大，成爲峯值。所以，在檢測結果中，一個峯值點，就對應一個圓心點。
霍夫圓檢測的函數：

skimage.transform.hough_circle(image, radius) 

radius是一個數組，表示半徑的集合，如[3，4，5，6]
返回一個3維的數組（radius index, M, N), 第一維表示半徑的索引，後面兩維表示圖像的尺寸。
例1：繪製兩個圓形，用霍夫圓變換將它們檢測出來。

import numpy as np import matplotlib.pyplot as plt from skimage import draw,transform,feature img = np.zeros((250, 250,3), dtype=np.uint8) rr, cc = draw.circle_perimeter(60, 60, 50) #以半徑50畫一個圓 rr1, cc1 = draw.circle_perimeter(150, 150, 60) #以半徑60畫一個圓 img[cc, rr,:] =255 img[cc1, rr1,:] =255 fig, (ax0,ax1) = plt.subplots(1,2, figsize=(8, 5)) ax0.imshow(img) #顯示原圖 ax0.set_title('origin image') hough_radii = np.arange(50, 80, 5) #半徑範圍 hough_res =transform.hough_circle(img[:,:,0], hough_radii) #圓變換 centers = [] #保存全部圓心點座標 accums = [] #累積值 radii = [] #半徑 for radius, h in zip(hough_radii, hough_res): #每個半徑值，取出其中兩個圓 num_peaks = 2 peaks =feature.peak_local_max(h, num_peaks=num_peaks) #取出峯值 centers.extend(peaks) accums.extend(h[peaks[:, 0], peaks[:, 1]]) radii.extend([radius] * num_peaks) #畫出最接近的圓 image =np.copy(img) for idx in np.argsort(accums)[::-1][:2]: center_x, center_y = centers[idx] radius = radii[idx] cx, cy =draw.circle_perimeter(center_y, center_x, radius) image[cy, cx] =(255,0,0) ax1.imshow(image) ax1.set_title('detected image') 

結果圖以下：原圖中的圓用白色繪製，檢測出的圓用紅色繪製。

 

 

 

例2，檢測出下圖中存在的硬幣。

 

 

import numpy as np import matplotlib.pyplot as plt from skimage import data, color,draw,transform,feature,util image = util.img_as_ubyte(data.coins()[0:95, 70:370]) #裁剪原圖片 edges =feature.canny(image, sigma=3, low_threshold=10, high_threshold=50) #檢測canny邊緣 fig, (ax0,ax1) = plt.subplots(1,2, figsize=(8, 5)) ax0.imshow(edges, cmap=plt.cm.gray) #顯示canny邊緣 ax0.set_title('original iamge') hough_radii = np.arange(15, 30, 2) #半徑範圍 hough_res =transform.hough_circle(edges, hough_radii) #圓變換 centers = [] #保存中心點座標 accums = [] #累積值 radii = [] #半徑 for radius, h in zip(hough_radii, hough_res): #每個半徑值，取出其中兩個圓 num_peaks = 2 peaks =feature.peak_local_max(h, num_peaks=num_peaks) #取出峯值 centers.extend(peaks) accums.extend(h[peaks[:, 0], peaks[:, 1]]) radii.extend([radius] * num_peaks) #畫出最接近的5個圓 image = color.gray2rgb(image) for idx in np.argsort(accums)[::-1][:5]: center_x, center_y = centers[idx] radius = radii[idx] cx, cy =draw.circle_perimeter(center_y, center_x, radius) image[cy, cx] = (255,0,0) ax1.imshow(image) ax1.set_title('detected image') 

 

橢圓變換是相似的，使用函數爲：

skimage.transform.hough_ellipse(img,accuracy, threshold, min_size, max_size) 

輸入參數：

img: 待檢測圖像。 accuracy: 使用在累加器上的短軸二進制尺寸，是一個double型的值，默認爲1 thresh: 累加器閾值，默認爲4 min_size: 長軸最小長度，默認爲4 max_size: 短軸最大長度，默認爲None,表示圖片最短邊的一半。

返回一個 [(accumulator, y0, x0, a, b, orientation)] 數組，accumulator表示累加器，（y0,x0)表示橢圓中心點，（a,b)分別表示長短軸，orientation表示橢圓方向

例：檢測出咖啡圖片中的橢圓杯口

import matplotlib.pyplot as plt from skimage import data,draw,color,transform,feature #加載圖片，轉換成灰度圖並檢測邊緣 image_rgb = data.coffee()[0:220, 160:420] #裁剪原圖像，否則速度很是慢 image_gray = color.rgb2gray(image_rgb) edges = feature.canny(image_gray, sigma=2.0, low_threshold=0.55, high_threshold=0.8) #執行橢圓變換 result =transform.hough_ellipse(edges, accuracy=20, threshold=250,min_size=100, max_size=120) result.sort(order='accumulator') #根據累加器排序 #估計橢圓參數 best = list(result[-1]) #排完序後取最後一個 yc, xc, a, b = [int(round(x)) for x in best[1:5]] orientation = best[5] #在原圖上畫出橢圓 cy, cx =draw.ellipse_perimeter(yc, xc, a, b, orientation) image_rgb[cy, cx] = (0, 0, 255) #在原圖中用藍色表示檢測出的橢圓 #分別用白色表示canny邊緣，用紅色表示檢測出的橢圓，進行對比 edges = color.gray2rgb(edges) edges[cy, cx] = (250, 0, 0) fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(8, 4)) ax1.set_title('Original picture') ax1.imshow(image_rgb) ax2.set_title('Edge (white) and result (red)') ax2.imshow(edges) plt.show() 

 

 

霍夫橢圓變換速度很是慢，應避免圖像太大。

---

 

在前面的python數字圖像處理（10）：圖像簡單濾波 中，咱們已經講解了不少算子用來檢測邊緣，其中用得最多的canny算子邊緣檢測。
本篇咱們講解一些其它方法來檢測輪廓。
一、查找輪廓（find_contours)
measure模塊中的find_contours()函數，可用來檢測二值圖像的邊緣輪廓。
函數原型爲：

skimage.measure.find_contours(array, level) 

array: 一個二值數組圖像
level: 在圖像中查找輪廓的級別值
返回輪廓列表集合，可用for循環取出每一條輪廓。
例1：

import numpy as np import matplotlib.pyplot as plt from skimage import measure,draw #生成二值測試圖像 img=np.zeros([100,100]) img[20:40,60:80]=1 #矩形 rr,cc=draw.circle(60,60,10) #小圓 rr1,cc1=draw.circle(20,30,15) #大圓 img[rr,cc]=1 img[rr1,cc1]=1 #檢測全部圖形的輪廓 contours = measure.find_contours(img, 0.5) #繪製輪廓 fig, (ax0,ax1) = plt.subplots(1,2,figsize=(8,8)) ax0.imshow(img,plt.cm.gray) ax1.imshow(img,plt.cm.gray) for n, contour in enumerate(contours): ax1.plot(contour[:, 1], contour[:, 0], linewidth=2) ax1.axis('image') ax1.set_xticks([]) ax1.set_yticks([]) plt.show() 

結果以下：不一樣的輪廓用不一樣的顏色顯示

 

 

 

例2：

import matplotlib.pyplot as plt from skimage import measure,data,color #生成二值測試圖像 img=color.rgb2gray(data.horse()) #檢測全部圖形的輪廓 contours = measure.find_contours(img, 0.5) #繪製輪廓 fig, axes = plt.subplots(1,2,figsize=(8,8)) ax0, ax1= axes.ravel() ax0.imshow(img,plt.cm.gray) ax0.set_title('original image') rows,cols=img.shape ax1.axis([0,rows,cols,0]) for n, contour in enumerate(contours): ax1.plot(contour[:, 1], contour[:, 0], linewidth=2) ax1.axis('image') ax1.set_title('contours') plt.show() 

 

二、逼近多邊形曲線
逼近多邊形曲線有兩個函數：subdivide_polygon（)和 approximate_polygon（）
subdivide_polygon（)採用B樣條（B-Splines)來細分多邊形的曲線，該曲線一般在凸包線的內部。
函數格式爲：

skimage.measure.subdivide_polygon(coords, degree=2, preserve_ends=False) 

coords: 座標點序列。
degree: B樣條的度數，默認爲2
preserve_ends: 若是曲線爲非閉合曲線，是否保存開始和結束點座標，默認爲false
返回細分爲的座標點序列。
approximate_polygon（）是基於Douglas-Peucker算法的一種近似曲線模擬。它根據指定的容忍值來近似一條多邊形曲線鏈，該曲線也在凸包線的內部。
函數格式爲:

skimage.measure.approximate_polygon(coords, tolerance) 

coords: 座標點序列
tolerance: 容忍值
返回近似的多邊形曲線座標序列。
例：

import numpy as np import matplotlib.pyplot as plt from skimage import measure,data,color #生成二值測試圖像 hand = np.array([[1.64516129, 1.16145833], [1.64516129, 1.59375], [1.35080645, 1.921875], [1.375, 2.18229167], [1.68548387, 1.9375], [1.60887097, 2.55208333], [1.68548387, 2.69791667], [1.76209677, 2.56770833], [1.83064516, 1.97395833], [1.89516129, 2.75], [1.9516129, 2.84895833], [2.01209677, 2.76041667], [1.99193548, 1.99479167], [2.11290323, 2.63020833], [2.2016129, 2.734375], [2.25403226, 2.60416667], [2.14919355, 1.953125], [2.30645161, 2.36979167], [2.39112903, 2.36979167], [2.41532258, 2.1875], [2.1733871, 1.703125], [2.07782258, 1.16666667]]) #檢測全部圖形的輪廓 new_hand = hand.copy() for _ in range(5): new_hand =measure.subdivide_polygon(new_hand, degree=2) # approximate subdivided polygon with Douglas-Peucker algorithm appr_hand =measure.approximate_polygon(new_hand, tolerance=0.02) print("Number of coordinates:", len(hand), len(new_hand), len(appr_hand)) fig, axes= plt.subplots(2,2, figsize=(9, 8)) ax0,ax1,ax2,ax3=axes.ravel() ax0.plot(hand[:, 0], hand[:, 1],'r') ax0.set_title('original hand') ax1.plot(new_hand[:, 0], new_hand[:, 1],'g') ax1.set_title('subdivide_polygon') ax2.plot(appr_hand[:, 0], appr_hand[:, 1],'b') ax2.set_title('approximate_polygon') ax3.plot(hand[:, 0], hand[:, 1],'r') ax3.plot(new_hand[:, 0], new_hand[:, 1],'g') ax3.plot(appr_hand[:, 0], appr_hand[:, 1],'b') ax3.set_title('all') 

 

---

高級形態學處理

形態學處理，除了最基本的膨脹、腐蝕、開/閉運算、黑/白帽處理外，還有一些更高級的運用，如凸包，連通區域標記，刪除小塊區域等。
一、凸包
凸包是指一個凸多邊形，這個凸多邊形將圖片中全部的白色像素點都包含在內。
函數爲：

skimage.morphology.convex_hull_image(image) 

輸入爲二值圖像，輸出一個邏輯二值圖像。在凸包內的點爲True, 不然爲False
例：

import matplotlib.pyplot as plt from skimage import data,color,morphology #生成二值測試圖像 img=color.rgb2gray(data.horse()) img=(img<0.5)*1 chull = morphology.convex_hull_image(img) #繪製輪廓 fig, axes = plt.subplots(1,2,figsize=(8,8)) ax0, ax1= axes.ravel() ax0.imshow(img,plt.cm.gray) ax0.set_title('original image') ax1.imshow(chull,plt.cm.gray) ax1.set_title('convex_hull image') 

 

 

convex_hull_image()是將圖片中的全部目標看做一個總體，所以計算出來只有一個最小凸多邊形。若是圖中有多個目標物體，每個物體須要計算一個最小凸多邊形，則須要使用convex_hull_object（）函數。
函數格式：

skimage.morphology.convex_hull_object(image, neighbors=8)

輸入參數image是一個二值圖像，neighbors表示是採用4連通仍是8連通，默認爲8連通。
例：

import matplotlib.pyplot as plt from skimage import data,color,morphology,feature #生成二值測試圖像 img=color.rgb2gray(data.coins()) #檢測canny邊緣,獲得二值圖片 edgs=feature.canny(img, sigma=3, low_threshold=10, high_threshold=50) chull = morphology.convex_hull_object(edgs) #繪製輪廓 fig, axes = plt.subplots(1,2,figsize=(8,8)) ax0, ax1= axes.ravel() ax0.imshow(edgs,plt.cm.gray) ax0.set_title('many objects') ax1.imshow(chull,plt.cm.gray) ax1.set_title('convex_hull image') plt.show() 

 

 

二、連通區域標記
在二值圖像中，若是兩個像素點相鄰且值相同（同爲0或同爲1），那麼就認爲這兩個像素點在一個相互連通的區域內。而同一個連通區域的全部像素點，都用同一個數值來進行標記，這個過程就叫連通區域標記。在判斷兩個像素是否相鄰時，咱們一般採用4連通或8連通判斷。在圖像中，最小的單位是像素，每一個像素周圍有8個鄰接像素，常見的鄰接關係有2種：4鄰接與8鄰接。4鄰接一共4個點，即上下左右，以下左圖所示。8鄰接的點一共有8個，包括了對角線位置的點，以下右圖所示。

 

 

 

在skimage包中，咱們採用measure子模塊下的label（）函數來實現連通區域標記。
函數格式：

skimage.measure.label（image,connectivity=None) 

參數中的image表示須要處理的二值圖像，connectivity表示鏈接的模式，1表明4鄰接，2表明8鄰接。
輸出一個標記數組（labels), 從0開始標記。

import numpy as np import scipy.ndimage as ndi from skimage import measure,color import matplotlib.pyplot as plt #編寫一個函數來生成原始二值圖像 def microstructure(l=256): n = 5 x, y = np.ogrid[0:l, 0:l] #生成網絡 mask = np.zeros((l, l)) generator = np.random.RandomState(1) #隨機數種子 points = l * generator.rand(2, n**2) mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 mask = ndi.gaussian_filter(mask, sigma=l/(4.*n)) #高斯濾波 return mask > mask.mean() data = microstructure(l=128)*1 #生成測試圖片 labels=measure.label(data,connectivity=2) #8連通區域標記 dst=color.label2rgb(labels) #根據不一樣的標記顯示不一樣的顏色 print('regions number:',labels.max()+1) #顯示連通區域塊數(從0開始標記) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) ax1.imshow(data, plt.cm.gray, interpolation='nearest') ax1.axis('off') ax2.imshow(dst,interpolation='nearest') ax2.axis('off') fig.tight_layout() plt.show() 

在代碼中，有些地方乘以1，則能夠將bool數組快速地轉換爲int數組。
結果如圖：有10個連通的區域，標記爲0-9

 

 

 

若是想分別對每個連通區域進行操做，好比計算面積、外接矩形、凸包面積等，則須要調用measure子模塊的regionprops（）函數。該函數格式爲：

skimage.measure.regionprops(label_image) 

返回全部連通區塊的屬性列表，經常使用的屬性列表以下表：

屬性名稱 類型 描述 area int 區域內像素點總數 bbox tuple 邊界外接框(min_row, min_col, max_row, max_col) centroid array　　 質心座標 convex_area int 凸包內像素點總數 convex_image ndarray 和邊界外接框同大小的凸包　　 coords ndarray 區域內像素點座標 Eccentricity float 離心率 equivalent_diameter float 和區域面積相同的圓的直徑 euler_number int　　 區域歐拉數 extent float 區域面積和邊界外接框面積的比率 filled_area int 區域和外接框之間填充的像素點總數 perimeter float 區域周長 label int 區域標記

三、刪除小塊區域
有些時候，咱們只須要一些大塊區域，那些零散的、小塊的區域，咱們就須要刪除掉，則可使用morphology子模塊的remove_small_objects（)函數。
函數格式：

skimage.morphology.remove_small_objects(ar, min_size=64, connectivity=1, in_place=False) 

參數：

ar: 待操做的bool型數組。 min_size: 最小連通區域尺寸，小於該尺寸的都將被刪除。默認爲64. connectivity: 鄰接模式，1表示4鄰接，2表示8鄰接 in_place: bool型值，若是爲True,表示直接在輸入圖像中刪除小塊區域，不然進行復制後再刪除。默認爲False. 返回刪除了小塊區域的二值圖像。

import numpy as np import scipy.ndimage as ndi from skimage import morphology import matplotlib.pyplot as plt #編寫一個函數來生成原始二值圖像 def microstructure(l=256): n = 5 x, y = np.ogrid[0:l, 0:l] #生成網絡 mask = np.zeros((l, l)) generator = np.random.RandomState(1) #隨機數種子 points = l * generator.rand(2, n**2) mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 mask = ndi.gaussian_filter(mask, sigma=l/(4.*n)) #高斯濾波 return mask > mask.mean() data = microstructure(l=128) #生成測試圖片 dst=morphology.remove_small_objects(data,min_size=300,connectivity=1) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) ax1.imshow(data, plt.cm.gray, interpolation='nearest') ax2.imshow(dst,plt.cm.gray,interpolation='nearest') fig.tight_layout() plt.show() 

在此例中，咱們將面積小於300的小塊區域刪除（由1變爲0），結果以下圖：

 

 

Keep the labels with 2 largest areas.

label_image =measure.label(newpr) areas = [r.area for r in regionprops(label_image)] areas.sort() if len(areas) > 2: for region in regionprops(label_image): if region.area < areas[-2]: for coordinates in region.coords: label_image[coordinates[0], coordinates[1]] = 0 binary = label_image > 0 label_image = morphology.remove_small_holes(binary, areas[-2]) 

areas = [r.area for r in regionprops(label_image)] areas.sort() if len(areas) > 2: for region in regionprops(label_image): if region.area < areas[-2]: for coordinates in region.coords: label_image[coordinates[0], coordinates[1]] = 0 binary = label_image > 0 if plot == True: plots[3].axis('off') plots[3].imshow(binary, cmap=plt.cm.bone) 

https://github.com/vivek14632/LungCancerProject/blob/master/preprocessing2.py

四、綜合示例：閾值分割+閉運算+連通區域標記+刪除小區塊+分色顯示

import numpy as np import matplotlib.pyplot as plt import matplotlib.patches as mpatches from skimage import data,filters,segmentation,measure,morphology,color #加載並裁剪硬幣圖片 image = data.coins()[50:-50, 50:-50] thresh =filters.threshold_otsu(image) #閾值分割 bw =morphology.closing(image > thresh, morphology.square(3)) #閉運算 cleared = bw.copy() #複製 cleared=segmentation.clear_border(cleared) #清除與邊界相連的目標物 cleared = sm.opening(cleared,sm.disk(2)) cleared = sm.closing(cleared,sm.disk(2)) label_image =measure.label(cleared) #連通區域標記 borders = np.logical_xor(bw, cleared) #異或 label_image[borders] = -1 image_label_overlay =color.label2rgb(label_image, image=image) #不一樣標記用不一樣顏色顯示 fig,(ax0,ax1)= plt.subplots(1,2, figsize=(8, 6)) ax0.imshow(cleared,plt.cm.gray) ax1.imshow(image_label_overlay) for region in measure.regionprops(label_image): #循環獲得每個連通區域屬性集 #忽略小區域 if region.area < 100: continue #繪製外包矩形 minr, minc, maxr, maxc = region.bbox rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr, fill=False, edgecolor='red', linewidth=2) ax1.add_patch(rect) fig.tight_layout() plt.show() 

 

---

骨架提取與分水嶺算法

骨架提取與分水嶺算法也屬於形態學處理範疇，都放在morphology子模塊內。
一、骨架提取
骨架提取，也叫二值圖像細化。這種算法能將一個連通區域細化成一個像素的寬度，用於特徵提取和目標拓撲表示。
morphology子模塊提供了兩個函數用於骨架提取，分別是Skeletonize（）函數和medial_axis（）函數。咱們先來看Skeletonize（）函數。
格式爲：s

kimage.morphology.skeletonize(image) 

輸入和輸出都是一幅二值圖像。
例1：

from skimage import morphology,draw import numpy as np import matplotlib.pyplot as plt #建立一個二值圖像用於測試 image = np.zeros((400, 400)) #生成目標對象1(白色U型) image[10:-10, 10:100] = 1 image[-100:-10, 10:-10] = 1 image[10:-10, -100:-10] = 1 #生成目標對象2（X型） rs, cs = draw.line(250, 150, 10, 280) for i in range(10): image[rs + i, cs] = 1 rs, cs = draw.line(10, 150, 250, 280) for i in range(20): image[rs + i, cs] = 1 #生成目標對象3（O型） ir, ic = np.indices(image.shape) circle1 = (ic - 135)**2 + (ir - 150)**2 < 30**2 circle2 = (ic - 135)**2 + (ir - 150)**2 < 20**2 image[circle1] = 1 image[circle2] = 0 #實施骨架算法 skeleton =morphology.skeletonize(image) #顯示結果 fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(8, 4)) ax1.imshow(image, cmap=plt.cm.gray) ax1.axis('off') ax1.set_title('original', fontsize=20) ax2.imshow(skeleton, cmap=plt.cm.gray) ax2.axis('off') ax2.set_title('skeleton', fontsize=20) fig.tight_layout() plt.show() 

生成一幅測試圖像，上面有三個目標對象，分別進行骨架提取，結果以下：

 

 

 

例2：利用系統自帶的馬圖片進行骨架提取

from skimage import morphology,data,color import matplotlib.pyplot as plt image=color.rgb2gray(data.horse()) image=1-image #反相 #實施骨架算法 skeleton =morphology.skeletonize(image) #顯示結果 fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(8, 4)) ax1.imshow(image, cmap=plt.cm.gray) ax1.axis('off') ax1.set_title('original', fontsize=20) ax2.imshow(skeleton, cmap=plt.cm.gray) ax2.axis('off') ax2.set_title('skeleton', fontsize=20) fig.tight_layout() plt.show() 

 

 

medial_axis就是中軸的意思，利用中軸變換方法計算前景（1值）目標對象的寬度，格式爲：

skimage.morphology.medial_axis(image, mask=None, return_distance=False) 

mask: 掩模。默認爲None, 若是給定一個掩模，則在掩模內的像素值才執行骨架算法。
return_distance: bool型值，默認爲False. 若是爲True, 則除了返回骨架，還將距離變換值也同時返回。這裏的距離指的是中軸線上的全部點與背景點的距離。

import numpy as np import scipy.ndimage as ndi from skimage import morphology import matplotlib.pyplot as plt #編寫一個函數，生成測試圖像 def microstructure(l=256): n = 5 x, y = np.ogrid[0:l, 0:l] mask = np.zeros((l, l)) generator = np.random.RandomState(1) points = l * generator.rand(2, n**2) mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 mask = ndi.gaussian_filter(mask, sigma=l/(4.*n)) return mask > mask.mean() data = microstructure(l=64) #生成測試圖像 #計算中軸和距離變換值 skel, distance =morphology.medial_axis(data, return_distance=True) #中軸上的點到背景像素點的距離 dist_on_skel = distance * skel fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) ax1.imshow(data, cmap=plt.cm.gray, interpolation='nearest') #用光譜色顯示中軸 ax2.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest') ax2.contour(data, [0.5], colors='w') #顯示輪廓線 fig.tight_layout() plt.show() 

 

 

二、分水嶺算法
分水嶺在地理學上就是指一個山脊，水一般會沿着山脊的兩邊流向不一樣的「匯水盆」。分水嶺算法是一種用於圖像分割的經典算法，是基於拓撲理論的數學形態學的分割方法。若是圖像中的目標物體是連在一塊兒的，則分割起來會更困難，分水嶺算法常常用於處理這類問題，一般會取得比較好的效果。
分水嶺算法能夠和距離變換結合，尋找「匯水盆地」和「分水嶺界限」，從而對圖像進行分割。二值圖像的距離變換就是每個像素點到最近非零值像素點的距離，咱們可使用scipy包來計算距離變換。
在下面的例子中，須要將兩個重疊的圓分開。咱們先計算圓上的這些白色像素點到黑色背景像素點的距離變換，選出距離變換中的最大值做爲初始標記點（若是是反色的話，則是取最小值），從這些標記點開始的兩個匯水盆越集越大，最後相交於分山嶺。從分山嶺處斷開，咱們就獲得了兩個分離的圓。
例1：基於距離變換的分山嶺圖像分割

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage as ndi from skimage import morphology,feature #建立兩個帶有重疊圓的圖像 x, y = np.indices((80, 80)) x1, y1, x2, y2 = 28, 28, 44, 52 r1, r2 = 16, 20 mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2 mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2 image = np.logical_or(mask_circle1, mask_circle2) #如今咱們用分水嶺算法分離兩個圓 distance = ndi.distance_transform_edt(image) #距離變換 local_maxi =feature.peak_local_max(distance, indices=False, footprint=np.ones((3, 3)), labels=image) #尋找峯值 markers = ndi.label(local_maxi)[0] #初始標記點 labels =morphology.watershed(-distance, markers, mask=image) #基於距離變換的分水嶺算法 fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(8, 8)) axes = axes.ravel() ax0, ax1, ax2, ax3 = axes ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest') ax0.set_title("Original") ax1.imshow(-distance, cmap=plt.cm.jet, interpolation='nearest') ax1.set_title("Distance") ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest') ax2.set_title("Markers") ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest') ax3.set_title("Segmented") for ax in axes: ax.axis('off') fig.tight_layout() plt.show() 

 

 

分水嶺算法也能夠和梯度相結合，來實現圖像分割。通常梯度圖像在邊緣處有較高的像素值，而在其它地方則有較低的像素值，理想狀況 下，分山嶺剛好在邊緣。所以，咱們能夠根據梯度來尋找分山嶺。
例2：基於梯度的分水嶺圖像分割

import matplotlib.pyplot as plt from scipy import ndimage as ndi from skimage import morphology,color,data,filter image =color.rgb2gray(data.camera()) denoised = filter.rank.median(image, morphology.disk(2)) #過濾噪聲 #將梯度值低於10的做爲開始標記點 markers = filter.rank.gradient(denoised, morphology.disk(5)) <10 markers = ndi.label(markers)[0] gradient = filter.rank.gradient(denoised, morphology.disk(2)) #計算梯度 labels =morphology.watershed(gradient, markers, mask=image) #基於梯度的分水嶺算法 fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(6, 6)) axes = axes.ravel() ax0, ax1, ax2, ax3 = axes ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest') ax0.set_title("Original") ax1.imshow(gradient, cmap=plt.cm.spectral, interpolation='nearest') ax1.set_title("Gradient") ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest') ax2.set_title("Markers") ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest') ax3.set_title("Segmented") for ax in axes: ax.axis('off') fig.tight_layout() plt.show() 

 

---

參考文獻
python數字圖像處理