Python中等效的Surface Curvature Matlab [英] Surface Curvature Matlab equivalent in Python

问题描述

我试图计算由点(x，y，z)数组给出的曲面的曲率.最初我试图拟合多项式方程z = a + bx + cx ^ 2 + dy + exy + fy ^ 2) 然后计算高斯曲率

I was trying to calculate the curvature of a surface given by array of points (x,y,z). Initially I was trying to fit a polynomial equation z=a + bx + cx^2 + dy + exy + fy^2) and then calculate the gaussian curvature

$ K = \ frac {F_ {xx} \ cdot F_ {yy}-{F_ {xy}} ^ 2} {(1+ {F_x} ^ 2 + {F_y} ^ 2)^ 2} $

$ K = \frac{F_{xx}\cdot F_{yy}-{F_{xy}}^2}{(1+{F_x}^2+{F_y}^2)^2} $

但是，如果表面复杂，则问题很明显.我发现此Matlab代码可通过数值计算曲率.我想知道如何在Python中做同样的事情.

However the problem is fitting if the surface is complex. I found this Matlab code to numerically calculate curvature. I wonder how to do the same in 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);

推荐答案

我希望我在这里还不算太晚.我遇到的问题完全相同(我所工作的公司的产品).

I hope I'm not too late here. I work with exactely the same problem (a product for the company I work to).

您必须考虑的第一件事是这些点必须代表一个矩形网格. X是2D数组，Y是2D数组，Z是2D数组.如果您具有非结构化的浊点，并且具有单个矩阵形状的Nx3(第一列为X，第二列为Y，第三列为Z)，则无法应用此matlab函数.

The first thing you must consider is that the points must represent a rectangular mesh. X is a 2D array, Y is a 2D array, and Z is a 2D array. If you have an unstructured cloudpoint, with a single matrix shaped Nx3 (the first column being X, the second being Y and the third being Z) then you can't apply this matlab function.

我已经开发了与此Matlab脚本等效的Python，在这里我只计算Z矩阵的平均曲率(我想您可能会受到脚本的启发，并对其进行调整以获取所有所需的曲率)，而忽略了X和Y通过假设网格是正方形.我认为您可以掌握"我的工作方式和方式，并使其适应您的需求:

I have developed a Python equivalent of this Matlab script, where I only calculate Mean curvature (I guess you can get inspired by the script and adapt it to get all your desired curvatures) for the Z matrix, ignoring the X and Y by assuming the grid is square. I think you can "grasp" what and how I am doing, and adapt it for your needs:

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

这篇关于Python中等效的Surface Curvature Matlab的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！