表面曲率Matlab等效Python

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

我试着计算由点(x,y,z)组成的曲面的曲率。最初我想适合一个多项式方程z = a + bx + cx ^ 2 + dy + exy + ^ 2)财政年度然后计算高斯曲率

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

$ K = \压裂{ f f { yy } - { xx } \ cdot f { xy } { } ^ 2 } {(1 + { F_x } ^ 2 + { F_y } ^ 2)^ 2 } $

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.

然而，如果表面是复杂的，问题就迎刃而解了。我找到了这个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);

5 个解决方案

#1

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

我希望我在这里不会太迟。我也遇到过同样的问题(我工作的公司的产品)。

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.

首先要考虑的是这些点必须代表一个矩形网格。X是2D数组，Y是2D数组，Z是2D数组。如果你有一个非结构化的cloudpoint，它有一个形状为Nx3的矩阵(第一列是X，第二列是Y，第三列是Z)那么你就不能应用这个matlab函数。

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:

我开发了一个Python版本的Matlab脚本，在这个脚本中，我只计算平均曲率(我猜你可以从脚本中得到灵感，调整它来得到你想要的所有曲率)。我认为你可以“掌握”我正在做什么和如何做，并根据你的需要进行调整:

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

#2

In case others stumble across this question, for completeness I offer the following code, inspired by heltonbiker.

如果其他人遇到这个问题，为了完整性，我提供以下代码，受heltonbiker的启发。

Here is some python code to calculate Gaussian curvature as described by equation (3) in "Computation of Surface Curvature from Range Images Using Geometrically Intrinsic Weights"*, T. Kurita and P. Boulanger, 1992.

下面是一些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

Note:

注意:

1. heltonbiker's method is essentially equation (4) from the paper
2. heltonbiker的方法本质上是来自于纸上的方程(4)。
3. heltonbiker's method is also the same on "Surfaces in 3D space, Mean Curvature" on Wikipedia: http://en.wikipedia.org/wiki/Mean_curvature)
4. heltonbiker的方法在*上的“三维空间表面，平均曲率”也是一样的:http://en.wikipedia.org/wiki/mean_曲率)
5. If you need both K and H then include the calculation of "K" (Gaussian curvature) in heltonbiker code and return K and H. Saves a little processing time.
6. 如果需要K和H，那么在heltonbiker代码中包含计算“K”(高斯曲率)，然后返回K和H，节省一点处理时间。
7. I assume the surface is defined as a function of two coordinates, e.g. z = Z(x, y). In my case Z is a range image.
8. 我假设曲面定义为两个坐标的函数，例如z = z (x, y)。

#3

Dot Product in Python

在Python中积

Derivates in Python

在Python中衍生

Reshaping in Python

在Python中重塑

Oddly enough all of these are SO questions. Take a gander around next time and you can likely find an answer. Also note that you'll want to be using NumPy for Python to do this. It's fairly intuitive to use. Matlibplot (or something like that) might be helpful for you too!

奇怪的是，所有这些都是问题。下次再找个借口，你可能会找到答案。还需要注意的是，您希望Python使用NumPy来实现这一点。使用起来很直观。Matlibplot(或类似的东西)可能对您也有帮助!

#4

Although very late, but no harm in posting. I modified the "surfature" function for use in Python. Disclaimer: I'm not the author original code. Credits wherever they are due.

虽然很晚了，但是发帖没有害处。我修改了在Python中使用的“surfature”函数。免责声明:我不是原作者。学分，无论什么时候到期。

    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

#5

BSD-licensed Python source code for surface fits can be found at

bsd许可的Python源代码的表面贴图可以在

https://github.com/zunzun/pyeq2

https://github.com/zunzun/pyeq2

(I'm the author).

作者(我)。