本文轉(zhuǎn)自 python數(shù)字圖像處理

霍夫線變換

在圖片處理中雪标，霍夫變換主要是用來檢測圖片中的幾何形狀贡珊，包括直線隧魄、圓、橢圓等淹冰。
在skimage中，霍夫變換是放在tranform模塊內(nèi)巨柒，本篇主要講解霍夫線變換樱拴。
對于平面中的一條直線，在笛卡爾坐標(biāo)系中潘拱，可用y=mx+b來表示疹鳄，其中m為斜率，b為截距芦岂。但是如果直線是一條垂直線瘪弓，則m為無窮大，所有通常我們在另一坐標(biāo)系中表示直線禽最，即極坐標(biāo)系下的r=xcos(theta)+ysin(theta)腺怯。即可用（r,theta）來表示一條直線。其中r為該直線到原點(diǎn)的距離川无，theta為該直線的垂線與x軸的夾角呛占。如下圖所示。

對于一個給定的點(diǎn)（x0,y0), 我們在極坐標(biāo)下繪出所有通過它的直線（r,theta)懦趋，將得到一條正弦曲線晾虑。如果將圖片中的所有非0點(diǎn)的正弦曲線都繪制出來，則會存在一些交點(diǎn)仅叫。所有經(jīng)過這個交點(diǎn)的正弦曲線帜篇，說明都擁有同樣的(r,theta), 意味著這些點(diǎn)在一條直線上。

skimage.transform.hough_line(img)

返回三個值：
h: 霍夫變換累積器
theta: 點(diǎn)與x軸的夾角集合掏呼，一般為0-179度
distance: 點(diǎn)到原點(diǎn)的距離坏快，即上面的所說的r.
例：

import skimage.transform as st
import numpy as np
import matplotlib.pyplot as plt
#matplotlib inline

# 構(gòu)建測試圖片
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變換所得數(shù)據(jù)
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()

從右邊那張圖可以看出，有兩個交點(diǎn)憎夷，說明原圖像中有兩條直線假消。
如果我們要把圖中的兩條直線繪制出來，則需要用到另外一個函數(shù)：

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

用這個函數(shù)可以取出峰值點(diǎn)岭接，即交點(diǎn)富拗，也即原圖中的直線臼予。
返回的參數(shù)與輸入的參數(shù)一樣。我們修改一下上邊的程序啃沪，在原圖中將兩直線繪制出來粘拾。

import skimage.transform as st
import numpy as np
import matplotlib.pyplot as plt

# 構(gòu)建測試圖片
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變換所得數(shù)據(jù)
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()

注意，繪制線條的時候创千，要從極坐標(biāo)轉(zhuǎn)換為笛卡爾坐標(biāo)缰雇，公式為：

skimage還提供了另外一個檢測直線的霍夫變換函數(shù)，
概率霍夫線變換：

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

參數(shù)：

img: 待檢測的圖像追驴。
threshold： 閾值械哟，可先項，默認(rèn)為10
line_length: 檢測的最短線條長度殿雪，默認(rèn)為50
line_gap: 線條間的最大間隙暇咆。增大這個值可以合并破碎的線條。默認(rèn)為10

返回：

lines: 線條列表, 格式如((x0, y0), (x1, y0))丙曙，標(biāo)明開始點(diǎn)和結(jié)束點(diǎn)爸业。

下面，我們用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))
# 創(chuàng)建顯示窗口.
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()

在極坐標(biāo)中，圓的表示方式為：

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

圓心為(x0,y0),r為半徑索抓，θ為旋轉(zhuǎn)度數(shù)钧忽，值范圍為0-359
如果給定圓心點(diǎn)和半徑，則其它點(diǎn)是否在圓上逼肯，我們就能檢測出來了惰瓜。在圖像中，我們將每個非0像素點(diǎn)作為圓心點(diǎn)汉矿，以一定的半徑進(jìn)行檢測，如果有一個點(diǎn)在圓上备禀，我們就對這個圓心累加一次洲拇。如果檢測到一個圓，那么這個圓心點(diǎn)就累加到最大曲尸，成為峰值赋续。因此，在檢測結(jié)果中另患，一個峰值點(diǎn)纽乱，就對應(yīng)一個圓心點(diǎn)。
霍夫圓檢測的函數(shù)：

skimage.transform.hough_circle(image, radius)

radius是一個數(shù)組昆箕，表示半徑的集合鸦列，如[3租冠，4，5薯嗤，6]
返回一個3維的數(shù)組（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 = []  #保存所有圓心點(diǎn)坐標(biāo)
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')

結(jié)果圖如下：原圖中的圓用白色繪制玻褪，檢測出的圓用紅色繪制肉渴。

例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 = []  #保存中心點(diǎn)坐標(biāo)
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')

橢圓變換是類似的，使用函數(shù)為：

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

輸入?yún)?shù)：

img: 待檢測圖像庸诱。
accuracy: 使用在累加器上的短軸二進(jìn)制尺寸捻浦，是一個double型的值，默認(rèn)為1
thresh: 累加器閾值桥爽，默認(rèn)為4
min_size: 長軸最小長度朱灿，默認(rèn)為4
max_size: 短軸最大長度，默認(rèn)為None,表示圖片最短邊的一半钠四。

返回一個 [(accumulator, y0, x0, a, b, orientation)] 數(shù)組盗扒，accumulator表示累加器，（y0,x0)表示橢圓中心點(diǎn)缀去，（a,b)分別表示長短軸侣灶，orientation表示橢圓方向

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

import matplotlib.pyplot as plt
from skimage import data,draw,color,transform,feature

#加載圖片，轉(zhuǎn)換成灰度圖并檢測邊緣
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)

#執(zhí)行橢圓變換
result =transform.hough_ellipse(edges, accuracy=20, threshold=250,min_size=100, max_size=120)
result.sort(order='accumulator') #根據(jù)累加器排序

#估計橢圓參數(shù)
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) #在原圖中用藍(lán)色表示檢測出的橢圓

#分別用白色表示canny邊緣褥影，用紅色表示檢測出的橢圓，進(jìn)行對比
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()

霍夫橢圓變換速度非常慢咏雌，應(yīng)避免圖像太大凡怎。

在前面的python數(shù)字圖像處理（10）：圖像簡單濾波 中，我們已經(jīng)講解了很多算子用來檢測邊緣赊抖，其中用得最多的canny算子邊緣檢測统倒。
本篇我們講解一些其它方法來檢測輪廓。
1氛雪、查找輪廓（find_contours)
measure模塊中的find_contours()函數(shù)房匆，可用來檢測二值圖像的邊緣輪廓。
函數(shù)原型為：

skimage.measure.find_contours(array, level)

array: 一個二值數(shù)組圖像
level: 在圖像中查找輪廓的級別值
返回輪廓列表集合，可用for循環(huán)取出每一條輪廓浴鸿。
例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()

結(jié)果如下：不同的輪廓用不同的顏色顯示

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

2井氢、逼近多邊形曲線
逼近多邊形曲線有兩個函數(shù)：subdivide_polygon（)和 approximate_polygon（）
subdivide_polygon（)采用B樣條（B-Splines)來細(xì)分多邊形的曲線，該曲線通常在凸包線的內(nèi)部赚楚。
函數(shù)格式為：

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

coords: 坐標(biāo)點(diǎn)序列毙沾。
degree: B樣條的度數(shù)，默認(rèn)為2
preserve_ends: 如果曲線為非閉合曲線宠页，是否保存開始和結(jié)束點(diǎn)坐標(biāo)左胞，默認(rèn)為false
返回細(xì)分為的坐標(biāo)點(diǎn)序列。
approximate_polygon（）是基于Douglas-Peucker算法的一種近似曲線模擬举户。它根據(jù)指定的容忍值來近似一條多邊形曲線鏈烤宙，該曲線也在凸包線的內(nèi)部。
函數(shù)格式為:

skimage.measure.approximate_polygon(coords, tolerance)

coords: 坐標(biāo)點(diǎn)序列
tolerance: 容忍值
返回近似的多邊形曲線坐標(biāo)序列俭嘁。
例：

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

高級形態(tài)學(xué)處理

形態(tài)學(xué)處理躺枕，除了最基本的膨脹、腐蝕供填、開/閉運(yùn)算拐云、黑/白帽處理外，還有一些更高級的運(yùn)用近她，如凸包叉瘩，連通區(qū)域標(biāo)記，刪除小塊區(qū)域等粘捎。
1薇缅、凸包
凸包是指一個凸多邊形，這個凸多邊形將圖片中所有的白色像素點(diǎn)都包含在內(nèi)攒磨。
函數(shù)為：

skimage.morphology.convex_hull_image(image)

輸入為二值圖像泳桦，輸出一個邏輯二值圖像。在凸包內(nèi)的點(diǎn)為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()是將圖片中的所有目標(biāo)看作一個整體娩缰，因此計算出來只有一個最小凸多邊形灸撰。如果圖中有多個目標(biāo)物體，每一個物體需要計算一個最小凸多邊形拼坎，則需要使用convex_hull_object（）函數(shù)浮毯。
函數(shù)格式：

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

輸入?yún)?shù)image是一個二值圖像，neighbors表示是采用4連通還是8連通演痒，默認(rèn)為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()

2趋惨、連通區(qū)域標(biāo)記
在二值圖像中鸟顺，如果兩個像素點(diǎn)相鄰且值相同（同為0或同為1），那么就認(rèn)為這兩個像素點(diǎn)在一個相互連通的區(qū)域內(nèi)。而同一個連通區(qū)域的所有像素點(diǎn)讯嫂，都用同一個數(shù)值來進(jìn)行標(biāo)記蹦锋，這個過程就叫連通區(qū)域標(biāo)記。在判斷兩個像素是否相鄰時欧芽，我們通常采用4連通或8連通判斷莉掂。在圖像中，最小的單位是像素千扔，每個像素周圍有8個鄰接像素憎妙，常見的鄰接關(guān)系有2種：4鄰接與8鄰接。4鄰接一共4個點(diǎn)曲楚，即上下左右厘唾，如下左圖所示。8鄰接的點(diǎn)一共有8個龙誊，包括了對角線位置的點(diǎn)抚垃，如下右圖所示。

在skimage包中趟大，我們采用measure子模塊下的label（）函數(shù)來實現(xiàn)連通區(qū)域標(biāo)記鹤树。
函數(shù)格式：

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

參數(shù)中的image表示需要處理的二值圖像，connectivity表示連接的模式逊朽，1代表4鄰接罕伯，2代表8鄰接。
輸出一個標(biāo)記數(shù)組（labels), 從0開始標(biāo)記惋耙。

import numpy as np
import scipy.ndimage as ndi
from skimage import measure,color
import matplotlib.pyplot as plt

#編寫一個函數(shù)來生成原始二值圖像
def microstructure(l=256):
    n = 5
    x, y = np.ogrid[0:l, 0:l]  #生成網(wǎng)絡(luò)
    mask = np.zeros((l, l))
    generator = np.random.RandomState(1)  #隨機(jī)數(shù)種子
    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連通區(qū)域標(biāo)記
dst=color.label2rgb(labels)  #根據(jù)不同的標(biāo)記顯示不同的顏色
print('regions number:',labels.max()+1)  #顯示連通區(qū)域塊數(shù)(從0開始標(biāo)記)

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數(shù)組快速地轉(zhuǎn)換為int數(shù)組绽榛。
結(jié)果如圖：有10個連通的區(qū)域湿酸，標(biāo)記為0-9

如果想分別對每一個連通區(qū)域進(jìn)行操作，比如計算面積灭美、外接矩形推溃、凸包面積等，則需要調(diào)用measure子模塊的regionprops（）函數(shù)届腐。該函數(shù)格式為：

skimage.measure.regionprops(label_image)

返回所有連通區(qū)塊的屬性列表铁坎，常用的屬性列表如下表：

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

3、刪除小塊區(qū)域
有些時候犁苏，我們只需要一些大塊區(qū)域硬萍，那些零散的、小塊的區(qū)域围详，我們就需要刪除掉朴乖，則可以使用morphology子模塊的remove_small_objects（)函數(shù)祖屏。
函數(shù)格式：

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

參數(shù)：

ar: 待操作的bool型數(shù)組。
min_size: 最小連通區(qū)域尺寸买羞，小于該尺寸的都將被刪除袁勺。默認(rèn)為64.
connectivity: 鄰接模式，1表示4鄰接畜普，2表示8鄰接
in_place: bool型值期丰，如果為True,表示直接在輸入圖像中刪除小塊區(qū)域，否則進(jìn)行復(fù)制后再刪除吃挑。默認(rèn)為False.
返回刪除了小塊區(qū)域的二值圖像钝荡。

import numpy as np
import scipy.ndimage as ndi
from skimage import morphology
import matplotlib.pyplot as plt

#編寫一個函數(shù)來生成原始二值圖像
def microstructure(l=256):
    n = 5
    x, y = np.ogrid[0:l, 0:l]  #生成網(wǎng)絡(luò)
    mask = np.zeros((l, l))
    generator = np.random.RandomState(1)  #隨機(jī)數(shù)種子
    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的小塊區(qū)域刪除（由1變?yōu)?）儒鹿，結(jié)果如下圖：

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

4化撕、綜合示例：閾值分割+閉運(yùn)算+連通區(qū)域標(biāo)記+刪除小區(qū)塊+分色顯示

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)) #閉運(yùn)算

cleared = bw.copy()  #復(fù)制
cleared=segmentation.clear_border(cleared)  #清除與邊界相連的目標(biāo)物
cleared = sm.opening(cleared,sm.disk(2)) 
cleared = sm.closing(cleared,sm.disk(2))

label_image =measure.label(cleared)  #連通區(qū)域標(biāo)記
borders = np.logical_xor(bw, cleared) #異或
label_image[borders] = -1
image_label_overlay =color.label2rgb(label_image, image=image) #不同標(biāo)記用不同顏色顯示

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): #循環(huán)得到每一個連通區(qū)域?qū)傩约?    
    #忽略小區(qū)域
    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()

骨架提取與分水嶺算法

骨架提取與分水嶺算法也屬于形態(tài)學(xué)處理范疇，都放在morphology子模塊內(nèi)约炎。
1植阴、骨架提取
骨架提取，也叫二值圖像細(xì)化圾浅。這種算法能將一個連通區(qū)域細(xì)化成一個像素的寬度掠手，用于特征提取和目標(biāo)拓?fù)浔硎尽?br> morphology子模塊提供了兩個函數(shù)用于骨架提取，分別是Skeletonize（）函數(shù)和medial_axis（）函數(shù)狸捕。我們先來看Skeletonize（）函數(shù)喷鸽。
格式為：s

kimage.morphology.skeletonize(image)

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

from skimage import morphology,draw
import numpy as np
import matplotlib.pyplot as plt

#創(chuàng)建一個二值圖像用于測試
image = np.zeros((400, 400))

#生成目標(biāo)對象1(白色U型)
image[10:-10, 10:100] = 1
image[-100:-10, 10:-10] = 1
image[10:-10, -100:-10] = 1

#生成目標(biāo)對象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

#生成目標(biāo)對象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)

#顯示結(jié)果
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()

生成一幅測試圖像灸拍，上面有三個目標(biāo)對象做祝，分別進(jìn)行骨架提取，結(jié)果如下：

例2：利用系統(tǒng)自帶的馬圖片進(jìn)行骨架提取

from skimage import morphology,data,color
import matplotlib.pyplot as plt

image=color.rgb2gray(data.horse())
image=1-image #反相
#實施骨架算法
skeleton =morphology.skeletonize(image)

#顯示結(jié)果
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值）目標(biāo)對象的寬度混槐，格式為：

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

mask: 掩模。默認(rèn)為None, 如果給定一個掩模轩性，則在掩模內(nèi)的像素值才執(zhí)行骨架算法声登。
return_distance: bool型值，默認(rèn)為False. 如果為True, 則除了返回骨架揣苏，還將距離變換值也同時返回悯嗓。這里的距離指的是中軸線上的所有點(diǎn)與背景點(diǎn)的距離。

import numpy as np
import scipy.ndimage as ndi
from skimage import morphology
import matplotlib.pyplot as plt

#編寫一個函數(shù)卸察，生成測試圖像
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)

#中軸上的點(diǎn)到背景像素點(diǎn)的距離
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()

2脯厨、分水嶺算法
分水嶺在地理學(xué)上就是指一個山脊，水通常會沿著山脊的兩邊流向不同的“匯水盆”坑质。分水嶺算法是一種用于圖像分割的經(jīng)典算法合武，是基于拓?fù)淅碚摰臄?shù)學(xué)形態(tài)學(xué)的分割方法个少。如果圖像中的目標(biāo)物體是連在一起的，則分割起來會更困難眯杏，分水嶺算法經(jīng)常用于處理這類問題，通常會取得比較好的效果壳澳。
分水嶺算法可以和距離變換結(jié)合岂贩，尋找“匯水盆地”和“分水嶺界限”，從而對圖像進(jìn)行分割巷波。二值圖像的距離變換就是每一個像素點(diǎn)到最近非零值像素點(diǎn)的距離萎津，我們可以使用scipy包來計算距離變換。
在下面的例子中抹镊，需要將兩個重疊的圓分開锉屈。我們先計算圓上的這些白色像素點(diǎn)到黑色背景像素點(diǎn)的距離變換，選出距離變換中的最大值作為初始標(biāo)記點(diǎn)（如果是反色的話垮耳，則是取最小值）颈渊，從這些標(biāo)記點(diǎn)開始的兩個匯水盆越集越大，最后相交于分山嶺终佛。從分山嶺處斷開俊嗽，我們就得到了兩個分離的圓。
例1：基于距離變換的分山嶺圖像分割

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage as ndi
from skimage import morphology,feature

#創(chuàng)建兩個帶有重疊圓的圖像
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)

#現(xiàn)在我們用分水嶺算法分離兩個圓
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] #初始標(biāo)記點(diǎn)
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()

分水嶺算法也可以和梯度相結(jié)合铃彰，來實現(xiàn)圖像分割绍豁。一般梯度圖像在邊緣處有較高的像素值，而在其它地方則有較低的像素值牙捉，理想情況 下竹揍，分山嶺恰好在邊緣。因此邪铲，我們可以根據(jù)梯度來尋找分山嶺芬位。
例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的作為開始標(biāo)記點(diǎn)
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()

參考文獻(xiàn)
python數(shù)字圖像處理