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} $
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);
Ёγ = κn n + κg n × ˙γ (1) Here κn is called the normal curvature and κg is the geodesic curvature of γ. n + κ2 g.
[K,H,P1,P2] = surfature(X,Y,Z) returns the gaussian curvature of a surface (K), mean curvature (H), and principal curvatures (P1,P2).
The Gaussian curvature of a surface at a point is defined as the product of the two principal normal curvatures; it is said to be positive if the principal normal curvatures curve in the same direction and negative if they curve in opposite directions.
The mean curvature H = (κ1 + κ2)/2, nevertheless, has its sign depending on the choice of n. The sign of K at a point p on a surface S has an important geometric meaning, which is detailed below.
I hope I'm not too late here. I work with exactly 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.
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:
note: This assumes that your data points are 1 unit apart.
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
In case others stumble across this question, for completeness I offer the following code, inspired by 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.
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:
Although very late, but no harm in posting. I modified the "surfature" function for use in Python. Disclaimer: I'm not the author original "surfature.m" code. Credits wherever they are due. Just presenting Python implementation.
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With