handle singular matrices

brick/processing/best_fit_ellipsoid.py

@@ -42,17 +42,22 @@ def ls_ellipsoid(xx: array_t, yy: array_t, zz: array_t) -> array_t:  # finds the
 
     JT = J.transpose()
     JTJ = np_.dot(JT, J)
-    InvJTJ = np_.linalg.inv(JTJ);
-    ABC = np_.dot(InvJTJ, np_.dot(JT, K))
+    try:
+        InvJTJ = np_.linalg.inv(JTJ);
+        ABC = np_.dot(InvJTJ, np_.dot(JT, K))
 
-    # Rearrange, move the 1 to the other side
-    #  Ax^2 + By^2 + Cz^2 +  Dxy +  Exz +  Fyz +  Gx +  Hy +  Iz - 1 = 0
-    #    or
-    #  Ax^2 + By^2 + Cz^2 +  Dxy +  Exz +  Fyz +  Gx +  Hy +  Iz + J = 0
-    #  where J = -1
-    ellipsoid_coef = np_.append(ABC, -1)
+        # Rearrange, move the 1 to the other side
+        #  Ax^2 + By^2 + Cz^2 +  Dxy +  Exz +  Fyz +  Gx +  Hy +  Iz - 1 = 0
+        #    or
+        #  Ax^2 + By^2 + Cz^2 +  Dxy +  Exz +  Fyz +  Gx +  Hy +  Iz + J = 0
+        #  where J = -1
+        ellipsoid_coef = np_.append(ABC, -1)
 
-    return (ellipsoid_coef)
+        return (ellipsoid_coef)
+
+    except:
+        print("Singular matrix, cannot find the best fitting ellipsoid.")
+        return (np_.zeros(9))
 
 
 def polyToParams3D(ellipsoid_coef: array_t, printMe: bool = False) -> tuple:
@@ -70,74 +75,78 @@ def polyToParams3D(ellipsoid_coef: array_t, printMe: bool = False) -> tuple:
     if printMe:
         print('\npolynomial\n', ellipsoid_coef)
 
-    Amat = np_.array(
-        [
-            [ellipsoid_coef[0], ellipsoid_coef[3] / 2.0, ellipsoid_coef[4] / 2.0, ellipsoid_coef[6] / 2.0],
-            [ellipsoid_coef[3] / 2.0, ellipsoid_coef[1], ellipsoid_coef[5] / 2.0, ellipsoid_coef[7] / 2.0],
-            [ellipsoid_coef[4] / 2.0, ellipsoid_coef[5] / 2.0, ellipsoid_coef[2], ellipsoid_coef[8] / 2.0],
-            [ellipsoid_coef[6] / 2.0, ellipsoid_coef[7] / 2.0, ellipsoid_coef[8] / 2.0, ellipsoid_coef[9]]
-        ])
-
-    if printMe:
-        print('\nAlgebraic form of polynomial\n', Amat)
-
-    # See B.Bartoni, Preprint SMU-HEP-10-14 Multi-dimensional Ellipsoidal Fitting
-    # equation 20 for the following method for finding the center
-    A3 = Amat[0:3, 0:3]
-    A3inv = inv(A3)
-    ofs = ellipsoid_coef[6:9] / 2.0
-    center = -np_.dot(A3inv, ofs)
-    if printMe:
-        print('\nCenter at:', center)
-
-    # Center the ellipsoid at the origin
-    Tofs = np_.eye(4)
-    Tofs[3, 0:3] = center
-    R = np_.dot(Tofs, np_.dot(Amat, Tofs.T))
-    if printMe:
-        print('\nAlgebraic form translated to center\n', R, '\n')
-
-    R3 = R[0:3, 0:3]
-    s1 = -R[3, 3]
-    R3S = R3 / s1
-
-    (el, ec) = eig(R3S)
-    el = np_.abs(el)
-
-    # Sorting in descending order of the eigenvalues and eigenvectors
-    sorting_vec = np_.argsort(-el)
-    el = el[sorting_vec]
-    ec = ec[:, sorting_vec]
-
-    # Calculating cartesian axes
-    recip = 1.0 / el
-    axes = np_.sqrt(recip)  # axes are in ascending order
-    if printMe:
-        print('\nAxes are\n', axes, '\n')
-
-    # Calculating the orientation
-    inverse = inv(ec)  # inverse is actually the transpose here
-    if printMe:
-        print('\nRotation matrix\n', inverse)
-
-    # Calculating spherical axes
-    # # eigenvector 1
-    r1 = np_.sqrt(sum(axe**2 for axe in ec[:, 0]))
-    phi1 = np_.arctan2(ec[1, 0], ec[0, 0])
-    theta1 = np_.arccos(ec[2, 0] / r1)
-    spherical_axes1 = (r1, theta1, phi1)
-    # # eigenvector 2
-    r2 = np_.sqrt(sum(axe**2 for axe in ec[:, 1]))
-    phi2 = np_.arctan2(ec[1, 1], ec[0, 1])
-    theta2 = np_.arccos(ec[2, 1] / r2)
-    spherical_axes2 = (r2, theta2, phi2)
-
-    spherical_axes = (spherical_axes1, spherical_axes2)
-
-    if printMe:
-        print('\nAxes in spherical coord are\n', spherical_axes, '\n')
-
-    return center, axes, inverse, spherical_axes
+    if ellipsoid_coef == (np_.zeros(9)):
+        return "NaN", "NaN", "NaN", "NaN"
+
+    else:
+        Amat = np_.array(
+            [
+                [ellipsoid_coef[0], ellipsoid_coef[3] / 2.0, ellipsoid_coef[4] / 2.0, ellipsoid_coef[6] / 2.0],
+                [ellipsoid_coef[3] / 2.0, ellipsoid_coef[1], ellipsoid_coef[5] / 2.0, ellipsoid_coef[7] / 2.0],
+                [ellipsoid_coef[4] / 2.0, ellipsoid_coef[5] / 2.0, ellipsoid_coef[2], ellipsoid_coef[8] / 2.0],
+                [ellipsoid_coef[6] / 2.0, ellipsoid_coef[7] / 2.0, ellipsoid_coef[8] / 2.0, ellipsoid_coef[9]]
+            ])
+
+        if printMe:
+            print('\nAlgebraic form of polynomial\n', Amat)
+
+        # See B.Bartoni, Preprint SMU-HEP-10-14 Multi-dimensional Ellipsoidal Fitting
+        # equation 20 for the following method for finding the center
+        A3 = Amat[0:3, 0:3]
+        A3inv = inv(A3)
+        ofs = ellipsoid_coef[6:9] / 2.0
+        center = -np_.dot(A3inv, ofs)
+        if printMe:
+            print('\nCenter at:', center)
+
+        # Center the ellipsoid at the origin
+        Tofs = np_.eye(4)
+        Tofs[3, 0:3] = center
+        R = np_.dot(Tofs, np_.dot(Amat, Tofs.T))
+        if printMe:
+            print('\nAlgebraic form translated to center\n', R, '\n')
+
+        R3 = R[0:3, 0:3]
+        s1 = -R[3, 3]
+        R3S = R3 / s1
+
+        (el, ec) = eig(R3S)
+        el = np_.abs(el)
+
+        # Sorting in descending order of the eigenvalues and eigenvectors
+        sorting_vec = np_.argsort(-el)
+        el = el[sorting_vec]
+        ec = ec[:, sorting_vec]
+
+        # Calculating cartesian axes
+        recip = 1.0 / el
+        axes = np_.sqrt(recip)  # axes are in ascending order
+        if printMe:
+            print('\nAxes are\n', axes, '\n')
+
+        # Calculating the orientation
+        inverse = inv(ec)  # inverse is actually the transpose here
+        if printMe:
+            print('\nRotation matrix\n', inverse)
+
+        # Calculating spherical axes
+        # # eigenvector 1
+        r1 = np_.sqrt(sum(axe**2 for axe in ec[:, 0]))
+        phi1 = np_.arctan2(ec[1, 0], ec[0, 0])
+        theta1 = np_.arccos(ec[2, 0] / r1)
+        spherical_axes1 = (r1, theta1, phi1)
+        # # eigenvector 2
+        r2 = np_.sqrt(sum(axe**2 for axe in ec[:, 1]))
+        phi2 = np_.arctan2(ec[1, 1], ec[0, 1])
+        theta2 = np_.arccos(ec[2, 1] / r2)
+        spherical_axes2 = (r2, theta2, phi2)
+
+        spherical_axes = (spherical_axes1, spherical_axes2)
+
+        if printMe:
+            print('\nAxes in spherical coord are\n', spherical_axes, '\n')
+
+        return center, axes, inverse, spherical_axes
 
 
 def GetConvexHull3D(soma_sites: site_h) -> tuple:
@@ -169,7 +178,12 @@ def FindBestFittingEllipsoid3D(soma: soma_t) -> tuple:
     # fit ellipsoid on the convex hull
     # # get ellipsoid polynomial coefficients
     ellipsoid_coef = ls_ellipsoid(convex_hull[0], convex_hull[1], convex_hull[2])
-    # # get ellipsoid 3D parameters
-    center, axes, orientation, spherical_coor = polyToParams3D(ellipsoid_coef, False)
 
-    return ellipsoid_coef, center, axes, orientation, spherical_coor, convex_hull, volume_hull
+    if ellipsoid_coef == (np_.zeros(9)):
+        return ellipsoid_coef, "NaN", "NaN", "NaN", "NaN", "NaN", "NaN"
+
+    else:
+        # # get ellipsoid 3D parameters
+        center, axes, orientation, spherical_coor = polyToParams3D(ellipsoid_coef, False)
+
+        return ellipsoid_coef, center, axes, orientation, spherical_coor, convex_hull, volume_hull