表面曲率Matlab等效于Python

Question

我试图计算由点阵列(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);

Answer 1

我希望我在这里不会太晚.我正在处理同样的问题(我工作的公司的产品).

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

我已经开发了一个类似于这个Matlab脚本的Python,我只计算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

Answer 2

如果其他人偶然发现这个问题,为了完整起见,我提供了以下代码,受到heltonbiker的启发.

这是一些用于计算高斯曲率的python代码,如"使用几何内在权重从范围图像计算表面曲率"*,T.Kurita和P.Boulanger,1992中的等式(3)所描述的.

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

注意:

1. heltonbiker的方法本质上是方程式(4)
2. heltonbiker的方法在维基百科上的"3D空间中的曲面,平均曲率"上也是相同的:http: //en.wikipedia.org/wiki/Mean_curvature)
3. 如果你需要K和H,那么在heltonbiker代码中包含"K"(高斯曲率)的计算并返回K和H.节省一点处理时间.
4. 我假设表面被定义为两个坐标的函数,例如z = Z(x,y).在我的情况下,Z是范围图像.

Answer 3

虽然很晚，但发帖无害。我修改了在 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

Answer 4

Python 中的点积

Python 中的派生类

用 Python 重塑

奇怪的是，所有这些都是SO问题。下次去看看，也许就能找到答案。另请注意，您需要使用 Python 的 NumPy 来执行此操作。使用起来相当直观。Matlibplot（或类似的东西）可能对你也有帮助！