python skimage图像处理(三)

作者: jiandanjinxin | 来源:发表于2017-09-26 16:57 被阅读1507次

    本文转自 python数字图像处理


    霍夫线变换

    在图片处理中,霍夫变换主要是用来检测图片中的几何形状,包括直线、圆、椭圆等。
    在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数字图像处理

    相关文章

      网友评论

        本文标题:python skimage图像处理(三)

        本文链接:https://www.haomeiwen.com/subject/qbmsextx.html