working on - curvature and tension- - phi theta ellipsoid -

brick/processing/best_fit_ellipsoid.py

@@ -55,7 +55,7 @@ def ls_ellipsoid(xx: array_t, yy: array_t, zz: array_t) -> array_t:  # finds the
     return (ellipsoid_coef)
 
 
-def polyToParams3D(ellipsoid_coef: array_t, printMe: bool) -> tuple:
+def polyToParams3D(ellipsoid_coef: array_t, spherical_coordinates: bool = False, printMe: bool = False) -> tuple:
     """
     gets 3D parameters of an ellipsoid.
     convert the polynomial form of the 3D-ellipsoid to parameters
@@ -94,24 +94,41 @@ def polyToParams3D(ellipsoid_coef: array_t, printMe: bool) -> tuple:
     if printMe: print('\nAlgebraic form translated to center\n', R, '\n')
 
     R3 = R[0:3, 0:3]
-    R3test = R3 / R3[0, 0]
-    # print('normed \n',R3test)
     s1 = -R[3, 3]
     R3S = R3 / s1
+
     (el, ec) = eig(R3S)
     print(el, ec)
+
     # Sorting in croissant order of the eigenvalues and eigenvectors
     sorting_vec = np_.argsort(el)
-    el = el[sorting_vec]
-    ec = ec[:, sorting_vec]
+    el2 = el[sorting_vec]
+    ec2 = ec[:, sorting_vec]
 
-    recip = 1.0 / np_.abs(el)
-    axes = np_.sqrt(recip)
-    if printMe: print('\nAxes are\n', axes, '\n')
+    # Calculating cartesian axes
+    recip2 = 1.0 / np_.abs(el2)
+    axes2 = np_.sqrt(recip2)
+    if printMe: print('\nAxes are\n', axes2, '\n')
 
-    inve = inv(ec)  # inverse is actually the transpose here
+    # Calculating the orientation
+    inve = inv(ec2)  # inverse is actually the transpose here
     if printMe: print('\nRotation matrix\n', inve)
-    return (center, axes, inve)
+
+    if spherical_coordinates:
+        # Calculating spherical axes
+        recip = 1.0 / np_.abs(el)
+        axes = np_.sqrt(recip)
+        r = np_.sqrt(sum(axe**2 for axe in axes))
+        phi = np_.arctan(axes[1]/axes[0])
+        theta = np_.arctan(axes[2]/r)
+        spherical_axes = (r, theta, phi)
+
+        if printMe: print('\nAxes in spherical coord are\n', spherical_axes, '\n')
+
+        return (center, axes2, inve, spherical_axes)
+
+    else:
+        return (center, axes2, inve)
 
 
 def GetConvexHull3D(soma_sites: site_h) -> tuple:

brick/processing/graph_feat_extraction.py

@@ -159,9 +159,10 @@ def ExtractFeaturesInDF(somas, size_voxel_in_micron: list, number_of_bins: int,
     somas_features_dict = {}  # Dict{soma 1: [features], soma 2: [features], ...}
     columns = [
         "Coef_V_soma__V_convex_hull",
-        "x_ellipsoid",
-        "y_ellipsoid",
-        "z_ellipsoid",
+        # "theta_a",
+        # "phi_a",
+        # "theta_b",
+        # "phi_b",
         "Coef_axes_ellips_y__x",
         "Coef_axes_ellips_z__x",
         #
@@ -292,6 +293,12 @@ def ExtractFeaturesInDF(somas, size_voxel_in_micron: list, number_of_bins: int,
             max_length = in_.ToMicron(soma.skl_graph.max_length, size_voxel_in_micron, decimals=decimals)
             std_lengths = np_.std(ext_lengths)
             entropy_lengths = st_.entropy(ext_lengths)
+            #
+            # Curvature
+            for _, _, edge in soma.skl_graph.edges.data("as_edge_t"):
+                if edge is not None:
+                    edge.SetEndPointDirections(size_voxel_in_micron)
+                    for point in
 
             # Find the thickness of the extensions
             for ___, ___, edge in soma.skl_graph.edges.data("as_edge_t"):
@@ -524,9 +531,6 @@ def ExtractFeaturesInDF(somas, size_voxel_in_micron: list, number_of_bins: int,
 
         somas_features_dict[f"soma {soma.uid}"] = [
             Coef_V_soma__V_convex_hull,
-            soma.axes_ellipsoid[0],
-            soma.axes_ellipsoid[1],
-            soma.axes_ellipsoid[2],
             Coef_axes_ellips_y__x,
             Coef_axes_ellips_z__x,
             N_nodes,

sklgraph/brick/edge.py

@@ -162,6 +162,41 @@ class edge_t:
                 f_dir, dtype=np_.float64
             ) / np_.linalg.norm(f_dir)
 
+    def FindCurvatureAndTorsion(self, size_voxel: list):
+        #
+        if self.as_curve is None:
+            self.SetCurveRepresentation(size_voxel=size_voxel)
+
+        if self.as_curve is not None:
+            derivates = {}
+            for d_idx in range(self.dim):
+                list = []
+                for point in self.as_curve[0].x:
+                    d1 = self.as_curve[d_idx](point, 1)
+                    d2 = self.as_curve[d_idx](point, 2)
+                    d3 = self.as_curve[d_idx](point, 3)
+                    list.append((d1, d2, d3))
+                derivates[f"x{d_idx+1}"] = list
+
+            curvatures = []
+            torsions = []
+            for idx in len(self.as_curve[0].x):
+
+                a = derivates["x3"][1] * derivates["x2"][0] - derivates["x2"][1] * derivates["x3"][0]
+                b = derivates["x1"][1] * derivates["x3"][0] - derivates["x3"][1] * derivates["x1"][0]
+                c = derivates["x2"][1] * derivates["x1"][0] - derivates["x1"][1] * derivates["x2"][0]
+
+                k = (np_.sqrt(a**2 + b**2 + c**2)) / (derivates["x1"][0]**2 + derivates["x2"][0]**2 + derivates["x3"][0]**2)**(3/2)
+
+                t = (derivates["x1"][2] * a +
+                     derivates["x2"][2] * b +
+                     derivates["x3"][2] * c) / (a**2 + b**2 + c**2)
+
+                curvatures.append(k)
+                torsions.append(t)
+
+        return curvatures, torsions
+
     @property
     def uid(self) -> str:
         """