From be5ef5afbb82e689078baea6141b4995b5dad444 Mon Sep 17 00:00:00 2001
From: javicarron <javicarron@gmail.com>
Date: Wed, 12 Apr 2023 20:56:36 +0200
Subject: [PATCH] Created general_curvature for 3D

This allows the computation of V2 and V3 for 3D DataFields
---
 pynkowski/stats/minkowski.py | 22 +++++++++++++++++++---
 1 file changed, 19 insertions(+), 3 deletions(-)

diff --git a/pynkowski/stats/minkowski.py b/pynkowski/stats/minkowski.py
index 0490f49..8864d84 100644
--- a/pynkowski/stats/minkowski.py
+++ b/pynkowski/stats/minkowski.py
@@ -169,10 +169,26 @@ def general_curvature(field, order):
             den = field.first_der[0]**2. + field.first_der[1]**2.
             return num / den
             
-#     if field.dim == 3:
-#         if order == 1: # Mean curvature (times 2), 2H = κ_1 + κ_2   (times |∇f|)
+    if field.dim == 3:
+        if order == 1: # Mean curvature (times 2), 2H = κ_1 + κ_2   (times |∇f|)
+            mod_grad = np.sqrt(np.sum(field.first_der**2., axis=0))
+            hessian = np.array([[field.second_der[0], field.second_der[3], field.second_der[4]],
+                                [field.second_der[3], field.second_der[1], field.second_der[5]],
+                                [field.second_der[4], field.second_der[5], field.second_der[2]]])
+            norm_grad = (field.first_der / mod_grad)
+            mean_curvature = np.einsum('j...,jk...,k...->...', norm_grad, hessian, norm_grad) - np.trace(hessian, axis1=0, axis2=1)
+            return mean_curvature
             
-#         if order == 2: # Gaussian curvature, K = κ_1 · κ_2   (times |∇f|)
+        if order == 2: # Gaussian curvature, K = κ_1 · κ_2   (times |∇f|)
+            mod_grad = np.sqrt(np.sum(field.first_der**2., axis=0))
+            norm_grad = (field.first_der / mod_grad)
+            norm_second_der = (field.second_der / mod_grad)
+            extended_hessian_det = np.linalg.det(np.array([[norm_second_der[0], norm_second_der[3], norm_second_der[4], norm_grad[0]],
+                                                           [norm_second_der[3], norm_second_der[1], norm_second_der[5], norm_grad[1]],
+                                                           [norm_second_der[4], norm_second_der[5], norm_second_der[2], norm_grad[2]],
+                                                           [norm_grad[0], norm_grad[1], norm_grad[2], np.zeros_like(norm_grad[0])]]).T).T
+            return -extended_hessian_det*mod_grad
+
         
     if field.dim > 3:
         raise NotImplementedError("The computation of principal curvatures in not implemented for spaces with 4 or more dimensions. If you need this, please get in touch.")