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算子边缘检测。
本篇我们讲解一些其它方法来检测轮廓。
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) 

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

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

 

---

参考文献
python数字图像处理