裂缝分割与裂缝宽度计算（正交⾻架线法）
分割
可以使⽤ opencv 的阈值分割，或者使⽤ CNN 等深度学习的⽅法进⾏裂缝分割，⼀般得到的分割结果如下，这⾥不再赘述。
寻边缘有两种⽅法，如下
scikit-image , 官⽅的⽰例如下：
另⼀种可以使⽤ Delaunay 三⾓⽹加上 alpha shape 进⾏边缘的提取，可以参考：
中轴变换
为了计算裂缝的宽度，⼀般使⽤正交⾻架线法，所以还需要计算裂缝的⾻架线，这⾥⼀般使⽤中轴变换，同样可以使⽤ scikit-image 库中的算法，参见：，官⽅⽰例：
宽度计算
下⾯是最重要的⼀部分，主要包括法向量估计和宽度计算，主要参见下⾯的两个函数 estimate_normal_for_pos（） 和
get_crack_ctrlpts（）
法向量估计
这⾥主要是查给定点的 K 近邻点，可以使⽤ kd-tree 算法，然后使⽤ 奇异值分解 （SVD）计算⾻架线的法向量。然后，裂缝的宽度计算就是在这个⽅向量的⽅向上计算的。
宽度计算
这⾥的主要思想是： 对弈给定的⼀个位置，根据上⽅，估计此位置⾻架线的⽅向量，并将此作为局部坐标系的y轴⽅向，x轴⽅向与此垂直，然后把裂缝的边缘线变换到这个局部坐标系中，然后在局部坐标系中，查到⾻架线法向量最近的两个裂缝边缘线点，并计算这两个点所形成的线段与⾻架线法向量的交点。对于左右侧的两个裂缝边缘线，分别计算。（这⼀部分具体可参见代码实现，后续会补充具体的图解）
代码解释
如上图所⽰，两条圆点线为边缘线，叉号线为⾻架线，也即中⼼线，中⼼线上任意⼀点的法向量可以通过对邻域点集拟合直线得到。对于Condition1，中⼼线上⼀点 , 根据 点的法向量建⽴局部坐标系，如上所⽰，根据 点的局部坐标系，可将裂缝边缘线上的点分为四个象限，为了到每⼀个象限内离 轴最近的点，定义参数hband（见下⾯代码），然后只保留点到y轴的距离⼩于hband的点，也就是上图中三条虚线内的点，即保留中⼼点 两侧⼀定带宽内的点，然后vband是为了规定边缘线的粗细的参数，当发现 上⽅或者下⽅（两条边缘线）任意⼀侧点集的 坐标的极差（np.ptp）⼤于vband，说明这时候检测到的边缘线存在噪点，因为噪点为导致边缘线在 ⽅向的极差变⼤，如图中Condition2中所标识处的噪点，噪点在vband带宽外，需要进⼀步过滤，过滤的⽅法是对存在噪点的地⽅，计算 坐标的均值，然后保留 坐标⼤于均值的点。通过hband,vband连个参数将中⼼点 处的边缘线分为4块区域，即图中的 ，最后，再在4块去区域中分别计算到 轴最近的点，即可以从来计算交点啦。
注：使⽤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)
c c c y c c y y y y c blt ,blb ,brt ,brb y
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 = py()
# 奇异值分解
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,:],:]
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]&an(bl[:,1])]
if np.ptp(blt[:,1])>vband:
blt = blt[blt[:,1]&an(blt[:,1])]
blb = bl[bl[:,1]&an(bl[:,1])]
if np.ptp(blb[:,1])>vband:
blb = blb[blb[:,1]&an(blb[:,1])]
brt = br[br[:,1]&an(br[:,1])]
if np.ptp(brt[:,1])>vband:
brt = brt[brt[:,1]&an(brt[:,1])]
brb = br[br[:,1]&an(br[:,1])]
if np.ptp(brb[:,1])>vband:
brb = brb[brb[:,1]&an(brb[:,1])]
#bh = np.vstack((bl,br))
#bmax = np.max(bh[:,1])
#bmin = np.min(bh[:,1])
#blt = bl[bl[:,1]>bmax-vband] # left top points
#blb = bl[bl[:,1]<bmin+vband] # left bottom points
#brt = br[br[:,1]>bmax-vband] # right top points
#brb = br[br[:,1]<bmin+vband] # right bottom points
t1 = blt[np.argsort(blt[:,0])[-1]]
t2 = brt[np.argsort(brt[:,0])[0]]
b1 = blb[np.argsort(blb[:,0])[-1]]
b2 = brb[np.argsort(brb[:,0])[0]]
interp1 =(ci[0]-t1[0])*((t2[1]-t1[1])/(t2[0]-t1[0]))+t1[1]
interp2 =(ci[0]-b1[0])*((b2[1]-b1[1])/(b2[0]-b1[0]))+b1[1]
if interp1-ci[1]>0and interp2-ci[1]<0:
widths.append([i,interp1-ci[1],interp2-ci[1]])
interps = np.array([[ci[0],interp1],[ci[0],interp2]])
interps_rec = np.dot(np.linalg.inv(tform),interps.T).T
#show_2dpoints([bpointsxl_loc1,bpointsxl_loc2,bpointsxr_loc1,bpointsxr_loc2,np.array([ptsl_1,ptsl_2]),np.array([ptsr_1,ptsr_2]),shap e(1,-1)],s=[1,1,1,1,20,20,20,20])
interps_rec = shape(1,-1)[0,:]
interp_segm.append(interps_rec)
except:
print("the %d-th was wrong"% i)
continue
interp_segm = np.array(interp_segm)
widths = np.array(widths)
# check
# show_2dpoints([np.array([[ci[0],interp1],[ci[0],interp2]]),np.array([t1,t2,b1,b2]),cpoints_loc,bl,br],[10,20,15,2,2])
return interp_segm, widths
path ="E:/Users/SubChange/Desktop/"
image = io.imread(path+"7Q3A9060-18.png", as_gray=True)
iw,ih = image.shape
blobs = np.copy(image)
blobs[blobs<128]=0
blobs[blobs>128]=1
blobs = blobs.astype(np.uint8)
# Generate the data
#blobs = data.binary_blobs(200, blob_size_fraction=.2,
#volume_fraction=.35, seed=1)
# using scikit-image
## Compute the medial axis (skeleton) and the distance transform
#skel, distance = medial_axis(blobs, return_distance=True)
## Distance to the background for pixels of the skeleton
#dist_on_skel = distance * skel
# Compare with other skeletonization algorithms
skeleton = skeletonize(blobs)
#skeleton_lee = skeletonize(blobs, method='lee')
x, y = np.where(skeleton>0)
centers = np.hstack((x.reshape(-1,1),y.reshape(-1,1)))
normals = estimate_normals(centers,3)
# search contours of the crack
contours = measure.find_contours(blobs,0.8)
bl = contours[0]
br = contours[1]
bpoints = np.vstack((bl,br))
#interp_segm, widths = get_crack_ctrlpts(centers,normals,bpoints,hband=2,vband=2)
bpixel = np.zeros((iw,ih,3),dtype=np.uint8)
bpoints = bpoints.astype(np.int)
bpixel[bpoints[:,0],bpoints[:,1],0]=255
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议