6

我试图计算由点数组(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);
4

7 回答 7

9

我希望我来这里不会太晚。我处理完全相同的问题(我工作的公司的产品)。

您必须考虑的第一件事是这些点必须代表一个矩形网格。X 是二维数组,Y 是二维数组,Z 是二维数组。如果您有一个非结构化浊点,具有单个矩阵形状 Nx3(第一列是 X,第二列是 Y,第三列是 Z),那么您不能应用这个 matlab 函数。

我已经开发了这个 Matlab 脚本的 Python 等价物,我只计算 Z 矩阵的平均曲率(我想你可以从脚本中得到启发并调整它以获得所有你想要的曲率),通过假设忽略 X 和 Y网格是正方形的。我认为您可以“掌握”我在做什么和如何做,并根据您的需要进行调整:

注意:这假设您的数据点相隔 1 个单位。

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 H
于 2013-03-19T20:27:02.950 回答
7

万一其他人偶然发现了这个问题,为了完整起见,我提供了以下受 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 的方法在 Wikipedia 上的“Surfaces in 3D space, Mean Curvature”上也是一样的: http ://en.wikipedia.org/wiki/Mean_curvature )
  3. 如果您同时需要 K 和 H,则在 heltonbiker 代码中包含“K”(高斯曲率)的计算并返回 K 和 H。节省一点处理时间。
  4. 我假设曲面被定义为两个坐标的函数,例如 z = Z(x, y)。在我的例子中 Z 是一个范围图像。
于 2013-11-12T04:43:03.827 回答
6

虽然很晚,但发帖无妨。我修改了“surfature”函数以在 Python 中使用。免责声明:我不是作者原创的“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
于 2015-06-09T22:46:02.427 回答
1

heltonbiker 关于平均曲率的回答很好,但它假设二维数组中的数据点彼此相距 1 个单位。例如,如果您的数据点彼此相距 0.3 个单位,则您需要将答案中的每个数据点除以 0.3 平方 (0.09) 以解决此问题。

对于迈克尔回答中的高斯曲率,您需要将每个数据点乘以(1/.3**2)**2(123.45)

于 2020-07-29T21:01:46.107 回答
0

可在以下位置找到 BSD 许可的用于表面拟合的 Python 源代码

https://github.com/zunzun/pyeq2

(我是作者)。

于 2012-07-06T14:29:05.633 回答
0

我建议使用该igl软件包,请注意它仍处于测试阶段。以下是如何计算igl 上的主曲率

目前你必须通过 conda 安装,但以后会在 PyPi 上发布轮子。

于 2021-07-07T13:47:06.733 回答
-1

Python中的点积

Python中的派生词

在 Python 中重塑

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

于 2012-07-03T19:17:31.713 回答