- android - 多次调用 OnPrimaryClipChangedListener
- android - 无法更新 RecyclerView 中的 TextView 字段
- android.database.CursorIndexOutOfBoundsException : Index 0 requested, 光标大小为 0
- android - 使用 AppCompat 时,我们是否需要明确指定其 UI 组件(Spinner、EditText)颜色
给定:
我在 csv 文件中有一组 3D 点。
X, Y, Z
x0, y0, z0
x1, y1, z1
...
xn, yn, zn
问题陈述:目标是根据最小二乘误差拟合平面。获取点与平面之间的误差距离。我可以自由选择Python或C++。我更喜欢 c++,但 python 也不错。
平面方程为:
Ax + By + Cz + D = 0 -- Equation 1
选项 1:C++ 方式
我在网上找到了这个链接,描述了一种解决平面拟合的c++方法 http://www.janssenprecisionengineering.com/downloads/Fit-plane-through-data-points.pdf但这里他们没有使用方程 1,而是使用方程 2
Z = A'x + B'y + C' -- Equation 2
where, A' = -A/C
B' = -B/C
C' = -D/C
一旦我知道了 A' B' C',我就可以根据方程 3 计算每个点到平面的距离(引用 -> Math Insight ):
Error = abs(A * x + B * y - z + C) / sqrt(pow(A, 2) + pow(B, 2) + 1); -- Equation 3
我开始用C++实现它
std::ofstream abc;
abc.open("Logs\\abc.csv");
cv::Mat Plane;
double Xi = 0;
double Yi = 0;
double Zi = 0;
double X2i = 0;
double Y2i = 0;
double XiYi = 0;
double XiZi = 0;
double YiZi = 0;
for (int o = 0; o < GL.PointX.size(); ++o) {
std::cout << "POINT X : " << GL.PointX[o] << std::endl;
std::cout << "POINT Y : " << GL.PointY[o] << std::endl;
std::cout << "POINT Z : " << GL.PointZ[o] << std::endl;
Xi = Xi + GL.PointX[o];
Yi = Yi + GL.PointY[o];
Zi = Zi + GL.PointZ[o];
X2i = X2i + (GL.PointX[o] * GL.PointX[o]);
Y2i = Y2i + (GL.PointY[o] * GL.PointY[o]);
XiYi = XiYi + (GL.PointX[o] * GL.PointY[o]);
XiZi = XiZi + (GL.PointX[o] * GL.PointZ[o]);
YiZi = YiZi + (GL.PointY[o] * GL.PointZ[o]);
}
cv::Mat PlaneA_Mat = (cv::Mat_<double>(3, 3) << X2i, XiYi, Xi, XiYi, Y2i, Yi, Xi, Yi, 1);
cv::Mat PlaneB_Mat(3, 1, CV_64FC1);
double R = 0, R2 = 0, FR2=0;
int Observation = 70;
PlaneB_Mat.at<double>(0, 0) = XiZi;
PlaneB_Mat.at<double>(1, 0) = YiZi;
PlaneB_Mat.at<double>(2, 0) = Zi;
Plane = PlaneA_Mat.inv() * PlaneB_Mat;
double A = Plane.at<double>(0, 0); //-A/C
double B = Plane.at<double>(1, 0); //-B/C
double C = Plane.at<double>(2, 0); //-D/C
abc << A << "," << B << "," << "-1" << "," << C;
abc <<"\n";
double Dsum = 0;
for (int o = 0; o < GL.PointX.size(); ++o) {
double Error = abs(A * GL.PointX[o] + B * GL.PointY[o] - GL.PointZ[o] + C) / sqrt(pow(A, 2) + pow(B, 2) + 1);
std::cout << "Error : " << Error << std::endl;
Error_projection << Error ;
Error_projection << "\n";
Dsum = Dsum + Error;
GL.D.push_back(Error);
}
R = (Observation * XiYi - Xi * Yi) / sqrt((Observation * X2i - Xi * Xi) * (Observation * Y2i - Yi * Yi));
R2 = pow(R, 2);
FR2 = pow(Zi - Dsum, 2) / pow(Zi - (1 / Observation)*Zi, 2);
std::cout << "PlaneA_Mat : " << PlaneA_Mat << std::endl;
std::cout << "PlaneB_Mat: " << PlaneB_Mat << std::endl;
std::cout << "FINAL PLANE: " << Plane << std::endl;
std::cout << "R: " << R << std::endl;
std::cout << " Correlation coefficient (R^2) : " << R2 << std::endl;
std::cout << " Final Correlation coefficient (FR^2) : " << FR2 << std::endl;
我得到的结果是:
A : 12.36346708893272
B : 0.07867292340114898
C : -3.714791111490779
并且还得到了每个点相对于平面的误差。 错误值非常大。
*选项 2:Python 方式*
然后我开始在网络上查找一些代码并发现了这个 --> 3D Plane Fit代码。我根据自己的需要进行了修改,效果非常好。
import numpy as np
import scipy.optimize
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import pandas
fig = plt.figure()
ax = fig.gca(projection='3d')
def fitPlaneLTSQ(XYZ):
(rows, cols) = XYZ.shape
print(rows, cols)
G = np.ones((rows, 3))
print(XYZ)
print(XYZ[0])
G[:, 0] = XYZ[0] # X
G[:, 1] = XYZ[1] # Y
print(G)
Z = XYZ[2]
(a, b, c), resid, rank, s = np.linalg.lstsq(G, Z)
print("a : ", a )
print("b : ", b )
print("c : ", c)
print("Residual : ", resid)
print("Rank : ", rank)
print("Singular Value", s)
normal = (a, b, -1) # I Don't Know WHY ?
print("normal : ", normal)
'''
normala = abs(pow(a,2)+pow(b,2)+pow(-1,2))
np.sqrt(normala)
abb=normal/normala
print("Normala : ",normala)
print("Normala : ", abb)
'''
nn = np.linalg.norm(normal) # I Don't Know WHY ?
print("nn : ", nn)
normal = normal/nn # I Don't Know WHY ?
print("Normal : ", normal)
return (c, normal)
#Import Data from CSV
result = pandas.read_csv("C:/Users/Logs/points_L.csv", header=None) # , names=['X', 'Y','Z']
#result =result.head(5)
print(result)
normal1 = pandas.read_csv("C:/Users/Logs/pose_left.csv", header=None) # , names=['X', 'Y','Z']
print(normal1)
abc = pandas.read_csv("C:/Users/Logs/abc.csv", header=None) # , names=['X', 'Y','Z']
print(abc)
#standard normal distribution / Bell.
#np.random.seed(seed=1)
data = result
#print(data)
print("NEW : ")
print(data)
c, normal = fitPlaneLTSQ(data)
print(c, normal)
# plot fitted plane
maxx = np.max(data[0])
maxy = np.max(data[1])
minx = np.min(data[0])
miny = np.min(data[1])
print(maxx,maxy, minx, miny)
point = np.array([0.0, 0.0, c]) # I Don't Know WHY ?
print("Point : ", point)
d = -point.dot(normal) # I Don't Know WHY ?
print("D : ", d)
# plot original points
ax.scatter(data[0], data[1], data[2])
ax.quiver(data[0], data[1], data[2], normal1[0], normal1[1], normal1[2], length=0.2)
# compute needed points for plane plotting
xx, yy = np.meshgrid([minx, maxx], [miny, maxy])
print(xx)
print(yy)
print("minx : ", minx)
print("maxx : ", maxx)
print("miny : ", miny)
print("maxy : ", maxy)
print("xx : ", xx)
print("yy : ", yy)
z = (-normal[0]*xx - normal[1]*yy - d)*1. / normal[2] # I Don't Know WHY ?
unit1 = np.sqrt(pow(normal[0], 2) + pow(normal[1],2) + pow(normal[2],2))
print("Unit 1 : ", unit1)
Error = abs(normal[0]*data[0] + normal[1]*data[1] + normal[2]*data[2] + d)/unit1
print("Error", Error)
Error_F = pandas.DataFrame(Error)
print("Print : ", Error_F)
Error_F.to_csv("C:/Users/Logs/Py_Error.csv")
print("Z : ", z)
# plot plane
ax.plot_surface(xx, yy, z, alpha=0.2)
ax.set_xlim(-1, 1)
ax.set_ylim(-1,1)
ax.set_zlim(1,2)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
如果您能回答以下问题,将会非常有帮助:
# I Dont Know WHY ?
您能否提供他们使用它背后的数学推理。请小心地解释它,就像你向菜鸟解释一样。 为什么我的最终相关系数 R^2 提供奇怪的结果?
Thank You Very Much !
最佳答案
嗯,这不是一个答案,但对于评论来说太长了。 “点与平面之间的距离”的另一种解释是,对于点 p
n.p + d
其中平面有方程
n.p + d = 0
我们假设
|n| = l
因此,如果我们有 N 个点 p[] 的集合,我们会寻找实数 d 和单位 vector n 来最小化
S = Sum{ i | (n.P[i] + d)*(n.P[i] + d) }
一些代数派生出以下计算方法:
a/通过计算 Pbar(均值点)和 Q(P[] 的“协方差”)
Pbar = Sum{ i | P[i]}/N
Q = Sum{ i | (P[i]-Pbar)*(P[i]-Pbar)' } / N
b/计算 n 和对应于最小特征值的 Q 特征向量。
c/计算
d = -n'*Pbar
关于python - 3D 平面拟合方法(C++ 与 Python),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/59881691/
我想了解 Ruby 方法 methods() 是如何工作的。 我尝试使用“ruby 方法”在 Google 上搜索,但这不是我需要的。 我也看过 ruby-doc.org,但我没有找到这种方法。
Test 方法 对指定的字符串执行一个正则表达式搜索,并返回一个 Boolean 值指示是否找到匹配的模式。 object.Test(string) 参数 object 必选项。总是一个
Replace 方法 替换在正则表达式查找中找到的文本。 object.Replace(string1, string2) 参数 object 必选项。总是一个 RegExp 对象的名称。
Raise 方法 生成运行时错误 object.Raise(number, source, description, helpfile, helpcontext) 参数 object 应为
Execute 方法 对指定的字符串执行正则表达式搜索。 object.Execute(string) 参数 object 必选项。总是一个 RegExp 对象的名称。 string
Clear 方法 清除 Err 对象的所有属性设置。 object.Clear object 应为 Err 对象的名称。 说明 在错误处理后,使用 Clear 显式地清除 Err 对象。此
CopyFile 方法 将一个或多个文件从某位置复制到另一位置。 object.CopyFile source, destination[, overwrite] 参数 object 必选
Copy 方法 将指定的文件或文件夹从某位置复制到另一位置。 object.Copy destination[, overwrite] 参数 object 必选项。应为 File 或 F
Close 方法 关闭打开的 TextStream 文件。 object.Close object 应为 TextStream 对象的名称。 说明 下面例子举例说明如何使用 Close 方
BuildPath 方法 向现有路径后添加名称。 object.BuildPath(path, name) 参数 object 必选项。应为 FileSystemObject 对象的名称
GetFolder 方法 返回与指定的路径中某文件夹相应的 Folder 对象。 object.GetFolder(folderspec) 参数 object 必选项。应为 FileSy
GetFileName 方法 返回指定路径(不是指定驱动器路径部分)的最后一个文件或文件夹。 object.GetFileName(pathspec) 参数 object 必选项。应为
GetFile 方法 返回与指定路径中某文件相应的 File 对象。 object.GetFile(filespec) 参数 object 必选项。应为 FileSystemObject
GetExtensionName 方法 返回字符串,该字符串包含路径最后一个组成部分的扩展名。 object.GetExtensionName(path) 参数 object 必选项。应
GetDriveName 方法 返回包含指定路径中驱动器名的字符串。 object.GetDriveName(path) 参数 object 必选项。应为 FileSystemObjec
GetDrive 方法 返回与指定的路径中驱动器相对应的 Drive 对象。 object.GetDrive drivespec 参数 object 必选项。应为 FileSystemO
GetBaseName 方法 返回字符串,其中包含文件的基本名 (不带扩展名), 或者提供的路径说明中的文件夹。 object.GetBaseName(path) 参数 object 必
GetAbsolutePathName 方法 从提供的指定路径中返回完整且含义明确的路径。 object.GetAbsolutePathName(pathspec) 参数 object
FolderExists 方法 如果指定的文件夹存在,则返回 True;否则返回 False。 object.FolderExists(folderspec) 参数 object 必选项
FileExists 方法 如果指定的文件存在返回 True;否则返回 False。 object.FileExists(filespec) 参数 object 必选项。应为 FileS
我是一名优秀的程序员,十分优秀!