Curvatura de superficie Matlab equivalente en Python

Estaba tratando de calcular la curvatura de una superficie dada por una matriz de puntos (x, y, z). Inicialmente intentaba ajustar una ecuación polinomial z = a + bx + cx ^ 2 + dy + exy + fy ^ 2) y luego calcular la curvatura gaussiana

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

Sin embargo, el problema es apropiado si la superficie es compleja. Encontré este código de Matlab para calcular numéricamente la curvatura. Me pregunto cómo hacer lo mismo en 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./[ppp]; % 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); 

Espero que no sea demasiado tarde aquí. Trabajo exactamente con el mismo problema (un producto para la empresa con la que trabajo).

Lo primero que debe considerar es que los puntos deben representar una malla rectangular. X es una matriz 2D, Y es una matriz 2D y Z es una matriz 2D. Si tiene un punto de nube no estructurado, con una matriz única con forma de Nx3 (la primera columna es X, la segunda es Y y la tercera es Z), entonces no puede aplicar esta función matlab.

He desarrollado un equivalente en Python de este script de Matlab, donde solo calculo la curvatura media (supongo que puede inspirarse en el script y adaptarlo para obtener todas las curvaturas deseadas) para la matriz Z, ignorando la X y la Y asumiendo que cuadrícula es cuadrada Creo que puedes “captar” qué estoy haciendo y cómo lo estoy haciendo, y adaptarlo a tus necesidades:

 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 

En caso de que otros tropiecen con esta pregunta, para completar, ofrezco el siguiente código, inspirado en heltonbiker.

Aquí hay un código de Python para calcular la curvatura gaussiana como se describe en la ecuación (3) en “Cálculo de la curvatura de la superficie a partir de imágenes de rango usando pesos geométricos intrínsecos” *, T. Kurita y 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 

Nota:

  1. El método de heltonbiker es esencialmente la ecuación (4) del artículo.
  2. El método de heltonbiker también es el mismo en “Superficies en espacio 3D, Curvatura media” en Wikipedia: http://en.wikipedia.org/wiki/Mean_curvature )
  3. Si necesita tanto K como H, incluya el cálculo de “K” (curvatura gaussiana) en el código de heltonbiker y devuelva K y H. Ahorra un poco de tiempo de procesamiento.
  4. Supongo que la superficie se define como una función de dos coordenadas, por ejemplo, z = Z (x, y). En mi caso Z es una imagen de rango.

Aunque es muy tarde, pero no hay daño en la publicación. Modifiqué la función “surfature” para usar en Python. Descargo de responsabilidad: No soy el código original del autor. Créditos a donde sean debidos.

  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*NM**2)/(E*GL**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 

Producto de punto en Python

Derivados en Python

Remodelación en Python

Por extraño que parezca, todas estas son tan preguntas. Eche un vistazo la próxima vez y probablemente pueda encontrar una respuesta. También tenga en cuenta que querrá usar NumPy for Python para hacer esto. Es bastante intuitivo de usar. ¡Matlibplot (o algo así) también podría ser útil para usted!

El código fuente Python con licencia BSD para ajustes de superficie se puede encontrar en

https://github.com/zunzun/pyeq2

(Soy el autor).