图像二值化之后,可以用cv2.CHAIN_APPROX_NONE或者cv2.CHAIN_APPROX_SIMPLE两种函数,NONE得到的是所有图像轮廓点,而SIMPLE得到的是“重要的点”,这个可以搜索该函数有详细的算法介绍。我用的NONE得到所有的像素点,并将其以csv格式保存。
再对上面保存的csv文件读取,进行B样条曲线拟合。可以修改拟合点的数量。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import math
import csv
from matplotlib.patches import Circle
from scipy.spatial.distance import pdist
import pandas as pd
from time import process_time
from time import time
#读取图像轮廓曲线(None/Simple),并进行B样条曲线拟合,可以设置确定数量的点来拟合(100个或者合适的点数量来拟合图像轮廓)
# 以准确率最高为准,并在此基础上求图像轮廓曲线的曲率。根据曲率的大小设置合适的参考值,也可称阈值来求得分叉点,输出分叉位置坐标
#以下代码是通过读取contours函数保存的csv文件进行B样条曲线拟合,并保存拟合后的csv取像,再读取相应的csv文件来求曲率,并得到曲率数值
#保存为csv文件,最终的程序可以选择通过读取保存csv文件作为中间步骤,但在时间上可以减少该步骤,从而减少程序运行的时间。
class BS_curve(object):
def __init__(self, n, p, cp=None, knots=None):
self.n = n # n+1 control points >>> p0,p1,,,pn
self.p = p
if cp:
self.cp = cp
self.u = knots
self.m = knots.shape[0] - 1 # m+1 knots >>> u0,u1,,,nm
else:
self.cp = None
self.u = None
self.m = None
self.paras = None
def check(self):
if self.m == self.n + self.p + 1:
return 1
else:
return 0
def coeffs(self, uq):
# n+1 control points >>> p0,p1,,,pn
# m+1 knots >>> u0,u1,,,nm
# algorithm is from https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve-coef.html
# N[] holds all intermediate and the final results
# in fact N is longer than control points,this is just to hold the intermediate value
# at last, we juest extract a part of N,that is N[0:n+1]
N = np.zeros(self.m + 1, dtype=np.float64)
# rule out special cases Important Properties of clamped B-spline curve
if uq == self.u[0]:
N[0] = 1.0
return N[0:self.n + 1]
elif uq == self.u[self.m]:
N[self.n] = 1.0
return N[0:self.n + 1]
# now u is between u0 and um
# first find k uq in span [uk,uk+1)
check = uq - self.u
ind = check >= 0
k = np.max(np.nonzero(ind))
# sk >>> multiplicity of u[k]
sk = np.sum(self.u == self.u[k])
N[k] = 1.0 # degree 0
# degree d goes from 1 to p
for d in range(1, self.p + 1):
r_max = self.m - d - 1 # the maximum subscript value of N in degree d,the minimum is 0
if k - d >= 0:
if self.u[k + 1] - self.u[k - d + 1]:
N[k - d] = (self.u[k + 1] - uq) / (self.u[k + 1] - self.u[k - d + 1]) * N[
k - d + 1] # right (south-west corner) term only
else:
N[k - d] = (self.u[k + 1] - uq) / 1 * N[k - d + 1] # right (south-west corner) term only
for i in range(k - d + 1, (k - 1) + 1):
if i >= 0 and i <= r_max:
Denominator1 = self.u[i + d] - self.u[i]
Denominator2 = self.u[i + d + 1] - self.u[i + 1]
# 0/0=0
if Denominator1 == 0:
Denominator1 = 1
if Denominator2 == 0:
Denominator2 = 1
N[i] = (uq - self.u[i]) / (Denominator1) * N[i] + (self.u[i + d + 1] - uq) / (Denominator2) * N[
i + 1]
if k <= r_max:
if self.u[k + d] - self.u[k]:
N[k] = (uq - self.u[k]) / (self.u[k + d] - self.u[k]) * N[k]
else:
N[k] = (uq - self.u[k]) / 1 * N[k]
return N[0:self.n + 1]
def De_Boor(self, uq):
# Input: a value u
# Output: the point on the curve, C(u)
# first find k uq in span [uk,uk+1)
check = uq - self.u
ind = check >= 0
k = np.max(np.nonzero(ind))
# inserting uq h times
if uq in self.u:
# sk >>> multiplicity of u[k]
sk = np.sum(self.u == self.u[k])
h = self.p - sk
else:
sk = 0
h = self.p
# rule out special cases
if h == -1:
if k == self.p:
return np.array(self.cp[0])
elif k == self.m:
return np.array(self.cp[-1])
# initial values of P(affected control points) >>> Pk-s,0 Pk-s-1,0 ... Pk-p+1,0
P = self.cp[k - self.p:k - sk + 1]
P = P.copy()
dis = k - self.p # the index distance between storage loaction and varibale i
# 1-h
for r in range(1, h + 1):
# k-p >> k-sk
temp = [] # uesd for Storing variables of the current stage
for i in range(k - self.p + r, k - sk + 1):
a_ir = (uq - self.u[i]) / (self.u[i + self.p - r + 1] - self.u[i])
temp.append((1 - a_ir) * P[i - dis - 1] + a_ir * P[i - dis])
P[k - self.p + r - dis:k - sk + 1 - dis] = np.array(temp)
# the last value is what we want
return P[-1]
def bs(self, us):
y = []
for x in us:
y.append(self.De_Boor(x))
y = np.array(y)
return y
def estimate_parameters(self, data_points, method="centripetal"):
pts = data_points.copy()
N = pts.shape[0]
w = pts.shape[1]
Li = []
for i in range(1, N):
Li.append(np.sum([pts[i, j] ** 2 for j in range(w)]) ** 0.5)
L = np.sum(Li)
t = [0]
for i in range(len(Li)):
Lki = 0
for j in range(i + 1):
Lki += Li[j]
t.append(Lki / L)
t = np.array(t)
self.paras = t
ind = t > 1.0
t[ind] = 1.0
return t
def get_knots(self, method="average"):
knots = np.zeros(self.p + 1).tolist()
paras_temp = self.paras.copy()
# m = n+p+1
self.m = self.n + self.p + 1
# we only need m+1 knots
# so we just select m+1-(p+1)-(p+1)+(p-1)+1+1 paras to average
num = self.m - self.p # select n+1 paras
ind = np.linspace(0, paras_temp.shape[0] - 1, num)
ind = ind.astype(int)
paras_knots = paras_temp[ind]
for j in range(1, self.n - self.p + 1):
k_temp = 0
# the maximun of variable i is n-1
for i in range(j, j + self.p - 1 + 1):
k_temp += paras_knots[i]
k_temp /= self.p
knots.append(k_temp)
add = np.ones(self.p + 1).tolist()
knots = knots + add
knots = np.array(knots)
self.u = knots
self.m = knots.shape[0] - 1
return knots
def set_paras(self, parameters):
self.paras = parameters
def set_knots(self, knots):
self.u = knots
def approximation(self, pts):
## Obtain a set of parameters t0, ..., tn
# pts_paras = self.estimate_parameters(pts)
## knot vector U;
# knots = self.get_knots()
num = pts.shape[0] - 1 # (num+1) is the number of data points
P = np.zeros((self.n + 1, pts.shape[1]), dtype=np.float64) # n+1 control points
P[0] = pts[0]
P[-1] = pts[-1]
# compute N
N = []
for uq in self.paras:
N_temp = self.coeffs(uq)
N.append(N_temp)
N = np.array(N)
Q = [0] # hold the location
for k in range(1, num - 1 + 1):
Q_temp = pts[k] - N[k, 0] * pts[0] - N[k, self.n] * pts[-1]
Q.append(Q_temp)
b = [0]
for i in range(1, self.n - 1 + 1):
b_temp = 0
for k in range(1, num - 1 + 1):
b_temp += N[k, i] * Q[k]
b.append(b_temp)
b = b[1::]
b = np.array(b)
N = N[:, 1:(self.n - 1) + 1]
A = np.dot(N.T, N)
cpm = np.linalg.solve(A, b)
P[1:self.n] = cpm
self.cp = P
return P
if __name__ == "__main__":
process_time()
bs = BS_curve(50, 3)
# xx = np.linspace(0, 4 * np.pi, 101)
# yy = np.sin(xx) + 0.6 * np.random.random(101)
fig = plt.figure(figsize=(5, 10))
ax = fig.add_subplot(111)
# ax.scatter(xx, yy)
circle = Circle(xy=(1415,748),radius=5,alpha=1,color='b')
ax.add_patch(circle)
import pandas as pd
# data = csv.reader(open('C:/Users/86185/Desktop/cef.csv'), 'r')
data = pd.read_csv('C:/Users/86185/Desktop/4bmp.csv')
data = np.array(data)
# data = np.array([xx, yy]).T
paras = bs.estimate_parameters(data)
knots = bs.get_knots()
if bs.check():
cp = bs.approximation(data)
uq = np.linspace(0, 1, 401)
y = bs.bs(uq)
其中uq = np.linspace(0,1,101)表示最终拟合剩下的点。
参考的这位大佬的文章
https://blog.csdn.net/Subtlechange/article/details/107170966?ops_request_misc=%25257B%252522request%25255Fid%252522%25253A%252522160856178416780277847055%252522%25252C%252522scm%252522%25253A%25252220140713.130102334.pc%25255Fall.%252522%25257D&request_id=160856178416780277847055&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2allfirst_rank_v2~rank_v29-1-107170966.pc_search_result_cache&utm_term=B%E6%A0%B7%E6%9D%A1%E6%9B%B2%E7%BA%BF%E6%8B%9F%E5%90%88%E8%BD%AE%E5%BB%93
得到数据点以csv格式保存
下面利用上面的csv文件进行曲率计算。
data1 = pd.read_csv('C:/Users/86185/Desktop/4bpine.csv')
data1 = np.array(data1)
start = time()
a, b= data1.shape
K = []
for j in range(a-2):
x1 = data1[j][0]
y1 = data1[j][1]
x2 = data1[j+1][0]
y2 = data1[j+1][1]
x3 = data1[j+2][0]
y3 = data1[j+2][1]
# print(x1, x2, x3)
# print(y1, y2, y3)
x = (x2 -x1, y2-y1)
y = (x3-x2,y3-y2)
d = 1 - pdist([x,y],'cosine')
sin = np.sqrt(1 - d**2)
dis = np.sqrt((x1-x3)**2 + (y1-y3)**2)
k = 2*sin/dis
# print(k[0])
K.append(k[0])
当然你也可以不用以csv格式,直接用数据进行计算,得到你想要的结果。
曲率计算参考的这位博主的文章。
https://blog.csdn.net/luxinfeng666/article/details/101475986?ops_request_misc=%25257B%252522request%25255Fid%252522%25253A%252522160938408016780302953073%252522%25252C%252522scm%252522%25253A%25252220140713.130102334…%252522%25257D&request_id=160938408016780302953073&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2allsobaiduend~default-4-101475986.first_rank_v2_pc_rank_v29&utm_term=python%E6%B1%82%E6%9B%B2%E7%BA%BF%E6%9B%B2%E7%8E%87