python3维图形处理_python skimage图像处理(三)

霍夫线变换

在图片处理中，霍夫变换主要是用来检测图片中的几何形状，包括直线、圆、椭圆等。

在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算子边缘检测。

本篇我们讲解一些其它方法来检测轮廓。

1、查找轮廓(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()

2、逼近多边形曲线

逼近多边形曲线有两个函数：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')

高级形态学处理

形态学处理，除了最基本的膨胀、腐蚀、开/闭运算、黑/白帽处理外，还有一些更高级的运用，如凸包，连通区域标记，删除小块区域等。

1、凸包

凸包是指一个凸多边形，这个凸多边形将图片中所有的白色像素点都包含在内。

函数为：

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

2、连通区域标记

在二值图像中，如果两个像素点相邻且值相同(同为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 区域标记

3、删除小块区域

有些时候，我们只需要一些大块区域，那些零散的、小块的区域，我们就需要删除掉，则可以使用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)

4、综合示例：阈值分割+闭运算+连通区域标记+删除小区块+分色显示

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子模块内。

1、骨架提取

骨架提取，也叫二值图像细化。这种算法能将一个连通区域细化成一个像素的宽度，用于特征提取和目标拓扑表示。

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

2、分水岭算法

分水岭在地理学上就是指一个山脊，水通常会沿着山脊的两边流向不同的“汇水盆”。分水岭算法是一种用于图像分割的经典算法，是基于拓扑理论的数学形态学的分割方法。如果图像中的目标物体是连在一起的，则分割起来会更困难，分水岭算法经常用于处理这类问题，通常会取得比较好的效果。

分水岭算法可以和距离变换结合，寻找“汇水盆地”和“分水岭界限”，从而对图像进行分割。二值图像的距离变换就是每一个像素点到最近非零值像素点的距离，我们可以使用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()