Python中的曲面曲率Matlab等效
Python中的曲面曲率Matlab等效
提问于 2012-07-03 19:11:26
我试图计算由点数组(x，y，z)给出的曲面的曲率。首先，我试图拟合一个多项式方程z=a + bx + cx^2 + dy + exy + fy^2)，然后计算高斯曲率。
$K= \frac{F_{xx}\cdot F_{yy}-{F_{xy}^2}{(1+{F_x}^2+{F_y}^2)^2}$
然而，如果表面是复杂的，那么这个问题是合适的。我找到了这个Matlab程序来数值计算曲率。我想知道如何在Python中做同样的事情。
function [K,H,Pmax,Pmin] = surfature(X,Y,Z),
% SURFATURE - COMPUTE GAUSSIAN AND MEAN CURVATURES OF A SURFACE
% [K,H] = SURFATURE(X,Y,Z), WHERE X,Y,Z ARE 2D ARRAYS OF POINTS ON THE
% SURFACE. K AND H ARE THE GAUSSIAN AND MEAN CURVATURES, RESPECTIVELY.
% SURFATURE RETURNS 2 ADDITIONAL ARGUEMENTS,
% [K,H,Pmax,Pmin] = SURFATURE(...), WHERE Pmax AND Pmin ARE THE MINIMUM
% AND MAXIMUM CURVATURES AT EACH POINT, RESPECTIVELY.
% First Derivatives
[Xu,Xv] = gradient(X);
[Yu,Yv] = gradient(Y);
[Zu,Zv] = gradient(Z);
% Second Derivatives
[Xuu,Xuv] = gradient(Xu);
[Yuu,Yuv] = gradient(Yu);
[Zuu,Zuv] = gradient(Zu);
[Xuv,Xvv] = gradient(Xv);
[Yuv,Yvv] = gradient(Yv);
[Zuv,Zvv] = gradient(Zv);
% Reshape 2D Arrays into Vectors
Xu = Xu(:); Yu = Yu(:); Zu = Zu(:);
Xv = Xv(:); Yv = Yv(:); Zv = Zv(:);
Xuu = Xuu(:); Yuu = Yuu(:); Zuu = Zuu(:);
Xuv = Xuv(:); Yuv = Yuv(:); Zuv = Zuv(:);
Xvv = Xvv(:); Yvv = Yvv(:); Zvv = Zvv(:);
Xu = [Xu Yu Zu];
Xv = [Xv Yv Zv];
Xuu = [Xuu Yuu Zuu];
Xuv = [Xuv Yuv Zuv];
Xvv = [Xvv Yvv Zvv];
% First fundamental Coeffecients of the surface (E,F,G)
E = dot(Xu,Xu,2);
F = dot(Xu,Xv,2);
G = dot(Xv,Xv,2);
m = cross(Xu,Xv,2);
p = sqrt(dot(m,m,2));
n = m./[p p p];
% Second fundamental Coeffecients of the surface (L,M,N)
L = dot(Xuu,n,2);
M = dot(Xuv,n,2);
N = dot(Xvv,n,2);
[s,t] = size(Z);
% Gaussian Curvature
K = (L.*N - M.^2)./(E.*G - F.^2);
K = reshape(K,s,t);
% Mean Curvature
H = (E.*N + G.*L - 2.*F.*M)./(2*(E.*G - F.^2));
H = reshape(H,s,t);
% Principal Curvatures
Pmax = H + sqrt(H.^2 - K);
Pmin = H - sqrt(H.^2 - K);回答 6
Stack Overflow用户
发布于 2013-03-19 20:27:02
我希望我在这里不会太晚。我处理的问题完全相同(这是我工作的公司的产品)。
首先要考虑的是，点必须表示矩形网格。X是2D数组，Y是2D数组，Z是2D数组。如果你有一个非结构化的云点，一个矩阵形状的Nx3 (第一列是X，第二列是Y，第三列是Z)，那么你就不能应用这个matlab函数。
我开发了一个Python等价于这个Matlab脚本，其中我只计算了Z矩阵的平均曲率(我猜你可以从这个脚本中得到灵感，并调整它以得到所有想要的曲率)，而忽略X和Y，假设网格是正方形的。我认为你可以“掌握”我在做什么和如何做，并根据你的需要调整它：
注意:这假设你的数据点相隔一个单位。
def mean_curvature(Z):
Zy, Zx = numpy.gradient(Z)
Zxy, Zxx = numpy.gradient(Zx)
Zyy, _ = numpy.gradient(Zy)
H = (Zx**2 + 1)*Zyy - 2*Zx*Zy*Zxy + (Zy**2 + 1)*Zxx
H = -H/(2*(Zx**2 + Zy**2 + 1)**(1.5))
return HStack Overflow用户
发布于 2013-11-12 04:43:03
如果其他人偶然发现了这个问题，为了完整起见，我提供了以下代码，这是由heltonbiker启发的。
这里有一些python程序来计算高斯曲率，如方程(3)所描述的“利用几何内在权重计算距离图像的表面曲率”*，T. Kurita和P. Boulanger，1992。
import numpy as np
def gaussian_curvature(Z):
Zy, Zx = np.gradient(Z)
Zxy, Zxx = np.gradient(Zx)
Zyy, _ = np.gradient(Zy)
K = (Zxx * Zyy - (Zxy ** 2)) / (1 + (Zx ** 2) + (Zy **2)) ** 2
return K注意：
- heltonbiker方法本质上是本文中的方程(4)
- heltonbiker的方法在维基百科的“三维曲面，平均曲率”上也是一样的：曲率)
- 如果你同时需要K和H，那么在直升机自行车代码中包括计算"K“(高斯曲率)，返回K和H会节省一些处理时间。
- 假设曲面定义为两个坐标的函数，例如z= Z(x，y)。在我的例子中，Z是一个范围图像。
Stack Overflow用户
发布于 2015-06-09 22:46:02
虽然很晚了，但在贴上没有害处。我修改了在Python中使用的"surfature“函数。免责声明:我不是作者最初的"surfature.m“代码。无论他们在哪里都要学分。只是介绍Python实现。
def surfature(X,Y,Z):
# where X, Y, Z matrices have a shape (lr+1,lb+1)
#First Derivatives
Xv,Xu=np.gradient(X)
Yv,Yu=np.gradient(Y)
Zv,Zu=np.gradient(Z)
#Second Derivatives
Xuv,Xuu=np.gradient(Xu)
Yuv,Yuu=np.gradient(Yu)
Zuv,Zuu=np.gradient(Zu)
Xvv,Xuv=np.gradient(Xv)
Yvv,Yuv=np.gradient(Yv)
Zvv,Zuv=np.gradient(Zv)
#Reshape to 1D vectors
nrow=(lr+1)*(lb+1) #total number of rows after reshaping
Xu=Xu.reshape(nrow,1)
Yu=Yu.reshape(nrow,1)
Zu=Zu.reshape(nrow,1)
Xv=Xv.reshape(nrow,1)
Yv=Yv.reshape(nrow,1)
Zv=Zv.reshape(nrow,1)
Xuu=Xuu.reshape(nrow,1)
Yuu=Yuu.reshape(nrow,1)
Zuu=Zuu.reshape(nrow,1)
Xuv=Xuv.reshape(nrow,1)
Yuv=Yuv.reshape(nrow,1)
Zuv=Zuv.reshape(nrow,1)
Xvv=Xvv.reshape(nrow,1)
Yvv=Yvv.reshape(nrow,1)
Zvv=Zvv.reshape(nrow,1)
Xu=np.c_[Xu, Yu, Zu]
Xv=np.c_[Xv, Yv, Zv]
Xuu=np.c_[Xuu, Yuu, Zuu]
Xuv=np.c_[Xuv, Yuv, Zuv]
Xvv=np.c_[Xvv, Yvv, Zvv]
#% First fundamental Coeffecients of the surface (E,F,G)
E=np.einsum('ij,ij->i', Xu, Xu)
F=np.einsum('ij,ij->i', Xu, Xv)
G=np.einsum('ij,ij->i', Xv, Xv)
m=np.cross(Xu,Xv,axisa=1, axisb=1)
p=sqrt(np.einsum('ij,ij->i', m, m))
n=m/np.c_[p,p,p]
#% Second fundamental Coeffecients of the surface (L,M,N)
L= np.einsum('ij,ij->i', Xuu, n)
M= np.einsum('ij,ij->i', Xuv, n)
N= np.einsum('ij,ij->i', Xvv, n)
#% Gaussian Curvature
K=(L*N-M**2)/(E*G-L**2)
K=K.reshape(lr+1,lb+1)
#% Mean Curvature
H = (E*N + G*L - 2*F*M)/(2*(E*G - F**2))
H = H.reshape(lr+1,lb+1)
#% Principle Curvatures
Pmax = H + sqrt(H**2 - K)
Pmin = H - sqrt(H**2 - K)
return Pmax,Pmin页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/11317579
复制相关文章
相似问题
腾讯云开发者
