python數(shù)字圖像處理之骨架提取與分水嶺算法

更新時(shí)間：2018年04月27日 09:21:01 作者：denny402

這篇文章主要介紹了python數(shù)字圖像處理之骨架提取與分水嶺算法，小編覺得挺不錯(cuò)的，現(xiàn)在分享給大家，也給大家做個(gè)參考。一起跟隨小編過來看看吧

骨架提取與分水嶺算法也屬于形態(tài)學(xué)處理范疇，都放在morphology子模塊內(nèi)。

1、骨架提取

骨架提取，也叫二值圖像細(xì)化。這種算法能將一個(gè)連通區(qū)域細(xì)化成一個(gè)像素的寬度，用于特征提取和目標(biāo)拓?fù)浔硎尽?/p>

morphology子模塊提供了兩個(gè)函數(shù)用于骨架提取，分別是Skeletonize（）函數(shù)和medial_axis（）函數(shù)。我們先來看Skeletonize（）函數(shù)。

格式為：skimage.morphology.skeletonize(image)

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

例1：

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

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

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

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

#實(shí)施骨架算法
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()

生成一幅測(cè)試圖像，上面有三個(gè)目標(biāo)對(duì)象，分別進(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 #反相
#實(shí)施骨架算法
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就是中軸的意思，利用中軸變換方法計(jì)算前景（1值）目標(biāo)對(duì)象的寬度，格式為：

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

mask: 掩模。默認(rèn)為None, 如果給定一個(gè)掩模，則在掩模內(nèi)的像素值才執(zhí)行骨架算法。

return_distance: bool型值，默認(rèn)為False. 如果為True, 則除了返回骨架，還將距離變換值也同時(shí)返回。這里的距離指的是中軸線上的所有點(diǎn)與背景點(diǎn)的距離。

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

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

#計(jì)算中軸和距離變換值
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é)上就是指一個(gè)山脊，水通常會(huì)沿著山脊的兩邊流向不同的“匯水盆”。分水嶺算法是一種用于圖像分割的經(jīng)典算法，是基于拓?fù)淅碚摰臄?shù)學(xué)形態(tài)學(xué)的分割方法。如果圖像中的目標(biāo)物體是連在一起的，則分割起來會(huì)更困難，分水嶺算法經(jīng)常用于處理這類問題，通常會(huì)取得比較好的效果。

分水嶺算法可以和距離變換結(jié)合，尋找“匯水盆地”和“分水嶺界限”，從而對(duì)圖像進(jìn)行分割。二值圖像的距離變換就是每一個(gè)像素點(diǎn)到最近非零值像素點(diǎn)的距離，我們可以使用scipy包來計(jì)算距離變換。

在下面的例子中，需要將兩個(gè)重疊的圓分開。我們先計(jì)算圓上的這些白色像素點(diǎn)到黑色背景像素點(diǎn)的距離變換，選出距離變換中的最大值作為初始標(biāo)記點(diǎn)（如果是反色的話，則是取最小值），從這些標(biāo)記點(diǎn)開始的兩個(gè)匯水盆越集越大，最后相交于分山嶺。從分山嶺處斷開，我們就得到了兩個(gè)分離的圓。

例1：基于距離變換的分山嶺圖像分割

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

#創(chuàng)建兩個(gè)帶有重疊圓的圖像
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)在我們用分水嶺算法分離兩個(gè)圓
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é)合，來實(shí)現(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)) #計(jì)算梯度
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()