基于图像的裂缝分割与裂缝宽度计算（正交骨架线法）

分割

可以使用 opencv 的阈值分割，或者使用 CNN 等深度学习的方法进行裂缝分割，一般得到的分割结果如下，这里不再赘述。

寻找边缘

寻找边缘有两种方法，如下

• scikit-image skimage.measure.find_contours, 官方的示例如下：

• 另一种可以使用 Delaunay 三角网加上 alpha shape 进行边缘的提取，可以参考： Delaunay三角网 alpha shape 2D点集边缘线提取

中轴变换

为了计算裂缝的宽度，一般使用正交骨架线法，所以还需要计算裂缝的骨架线，这里一般使用中轴变换，同样可以使用 scikit-image 库中的算法，参见：Skeletonize — skimage v0.18.0 docs (scikit-image.org)，官方示例：

宽度计算

下面是最重要的一部分，主要包括法向量估计和宽度计算，主要参见下面的两个函数 estimate_normal_for_pos（） 和 get_crack_ctrlpts（）

• 法向量估计

  这里主要是查找给定点的 K 近邻点，可以使用 kd-tree 算法，然后使用 奇异值分解 （SVD）计算骨架线的法向量。然后，裂缝的宽度计算就是在这个方向量的方向上计算的。

• 宽度计算

  这里的主要思想是： 对弈给定的一个位置，根据上方，估计此位置骨架线的方向量，并将此作为局部坐标系的y轴方向，x轴方向与此垂直，然后把裂缝的边缘线变换到这个局部坐标系中，然后在局部坐标系中，查找到骨架线法向量最近的两个裂缝边缘线点，并计算这两个点所形成的线段与骨架线法向量的交点。对于左右侧的两个裂缝边缘线，分别计算。（这一部分具体可参见代码实现，后续会补充具体的图解）

• 代码解释

如上图所示，两条圆点线为边缘线，叉号线为骨架线，也即中心线，中心线上任意一点的法向量可以通过对邻域点集拟合直线得到。对于Condition1，中心线上一点 c c c, 根据 c c c 点的法向量建立局部坐标系，如上所示，根据 c c c 点的局部坐标系，可将裂缝边缘线上的点分为四个象限，为了找到每一个象限内离 y y y 轴最近的点，定义参数hband（见下面代码），然后只保留点到y轴的距离小于hband的点，也就是上图中三条虚线内的点，即保留中心点 c c c 两侧一定带宽内的点，然后vband是为了规定边缘线的粗细的参数，当发现 c c c 上方或者下方（两条边缘线）任意一侧点集的 y y y 坐标的极差（np.ptp）大于vband，说明这时候检测到的边缘线存在噪点，因为噪点为导致边缘线在 y y y 方向的极差变大，如图中Condition2中所标识处的噪点，噪点在vband带宽外，需要进一步过滤，过滤的方法是对存在噪点的地方，计算 y y y 坐标的均值，然后保留 y y y 坐标大于均值的点。通过hband,vband连个参数将中心点 c c c 处的边缘线分为4块区域，即图中的 b l t , b l b , b r t , b r b blt,blb,brt,brb blt,blb,brt,brb，最后，再在4块去区域中分别计算到 y y y 轴最近的点，即可以从来计算交点啦。

注：使用hband，vband主要是将c点附近的边缘线先初步框定出来。

代码

import numpy as np
from skimage import io
from skimage.morphology import medial_axis, skeletonize
from skimage import measure
from skimage import data
import matplotlib.pyplot as plt
from sklearn.neighbors import KDTree


def show_2dpoints(pointcluster,s=None,quivers=None,qscale=1):
    # pointcluster should be a list of numpy ndarray
    # This functions would show a list of pint cloud in different colors
    n = len(pointcluster)
    nmax = n
    if quivers is not None:
        nq = len(quivers)
        nmax = max(n,nq)
    
    colors = ['r','g','b','c','m','y','k','tomato','gold']
    if nmax < 10:
        colors = np.array(colors[0:nmax])
    else: 
        colors = np.random.rand(nmax,3)

    fig = plt.figure(num=1)
    ax = fig.add_subplot(1,1,1)

    if s is None:
        s = np.ones(n)*2

    for i in range(n):
        ax.scatter(pointcluster[i][:,0],pointcluster[i][:,1],s=s[i],c=[colors[i]],alpha=0.6)

    if quivers is not None:
        for i in range(nq):
            ax.quiver(quivers[i][:,0],quivers[i][:,1],quivers[i][:,2],quivers[i][:,3],color=[colors[i]],scale=qscale)

    plt.show()


def SVD(points):
    # 二维，三维均适用
    # 二维直线，三维平面
    pts = points.copy()
    # 奇异值分解
    c = np.mean(pts, axis=0)
    A = pts - c # shift the points
    A = A.T #3*n
    u, s, vh = np.linalg.svd(A, full_matrices=False, compute_uv=True) # A=u*s*vh
    normal = u[:,-1]

    # 法向量归一化
    nlen = np.sqrt(np.dot(normal,normal))
    normal = normal / nlen
    # normal 是主方向的方向向量 与PCA最小特征值对应的特征向量是垂直关系
    # u 每一列是一个方向
    # s 是对应的特征值
    # c >>> 点的中心
    # normal >>> 拟合的方向向量
    return u,s,c,normal


def calcu_dis_from_ctrlpts(ctrlpts):
    if ctrlpts.shape[1]==4:
        return np.sqrt(np.sum((ctrlpts[:,0:2]-ctrlpts[:,2:4])**2,axis=1))
    else:
        return np.sqrt(np.sum((ctrlpts[:,[0,2]]-ctrlpts[:,[3,5]])**2,axis=1))


def estimate_normal_for_pos(pos,points,n):
    # estimate normal vectors at a given point
    pts = np.copy(points)
    tree = KDTree(pts, leaf_size=2)
    idx = tree.query(pos, k=n, return_distance=False, dualtree=False, breadth_first=False)
    #pts = np.concatenate((np.concatenate((pts[0].reshape(1,-1),pts),axis=0),pts[-1].reshape(1,-1)),axis=0)
    normals = []
    for i in range(0,pos.shape[0]):
        pts_for_normals = pts[idx[i,:],:]
        _,_,_,normal = SVD(pts_for_normals)
        normals.append(normal)
    normals = np.array(normals)
    return normals


def estimate_normals(points,n):
    pts = np.copy(points)
    tree = KDTree(pts, leaf_size=2)
    idx = tree.query(pts, k=n, return_distance=False, dualtree=False, breadth_first=False)
    #pts = np.concatenate((np.concatenate((pts[0].reshape(1,-1),pts),axis=0),pts[-1].reshape(1,-1)),axis=0)
    normals = []
    for i in range(0,pts.shape[0]):
        pts_for_normals = pts[idx[i,:],:]
        _,_,_,normal = SVD(pts_for_normals)
        normals.append(normal)
    normals = np.array(normals)
    return normals


def get_crack_ctrlpts(centers,normals,bpoints,hband=5,vband=2):
    # main algorithm to obtain crack width
    cpoints = np.copy(centers)
    cnormals = np.copy(normals)

    xmatrix = np.array([[0,1],[-1,0]])
    cnormalsx = np.dot(xmatrix,cnormals.T).T # the normal of x axis
    N = cpoints.shape[0]

    interp_segm = []
    widths = []
    for i in range(N):
        try:
            ny = cnormals[i]
            nx = cnormalsx[i]
            tform = np.array([nx,ny])
            bpoints_loc = np.dot(tform,bpoints.T).T
            cpoints_loc = np.dot(tform,cpoints.T).T
            ci = cpoints_loc[i]

            bl_ind = (bpoints_loc[:,0]-(ci[0]-hband))*(bpoints_loc[:,0]-ci[0])<0
            br_ind = (bpoints_loc[:,0]-ci[0])*(bpoints_loc[:,0]-(ci[0]+hband))<=0
            bl = bpoints_loc[bl_ind] # left points
            br = bpoints_loc[br_ind] # right points

            blt = bl[bl[:,1]>np.mean(bl[:,1])]
            if np.ptp(blt[:,1])>vband:
                blt = blt[blt[:,1]>np