From 386207c5ef7d580c7e9121515dae9b62b469765c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20M=C3=BCller?=
 <38459088+jo-mueller@users.noreply.github.com>
Date: Sat, 22 Jul 2023 11:24:52 +0200
Subject: [PATCH] fixed ellipsoid vector sorting

---
 .../_approximation/fit_ellipsoid.py            |  8 +++++---
 src/napari_stress/_tests/test_approximation.py | 18 +++++++++++-------
 2 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/src/napari_stress/_approximation/fit_ellipsoid.py b/src/napari_stress/_approximation/fit_ellipsoid.py
index 7fe4d3ab..dd36f244 100644
--- a/src/napari_stress/_approximation/fit_ellipsoid.py
+++ b/src/napari_stress/_approximation/fit_ellipsoid.py
@@ -32,11 +32,13 @@ def least_squares_ellipsoid(points: PointsData
     center, axes, R, R_inverse = polynomial_to_parameters3D(
         coefficients=coefficients)
     direction = R * axes[:, None]
-    origin = np.stack(3 * [center])  # cheap repeat
-    vector = np.stack([origin, direction]).transpose((1, 0, 2))
 
     # sort from longest to shortes vector
-    vector = vector[np.argsort(np.linalg.norm(vector[:, 1], axis=1))]
+    direction = direction[np.argsort(np.linalg.norm(direction, axis=1))][::-1]
+
+    # assemble vector
+    origin = np.stack(3 * [center])  # cheap repeat
+    vector = np.stack([origin, direction]).transpose((1, 0, 2))
 
     return vector
 
diff --git a/src/napari_stress/_tests/test_approximation.py b/src/napari_stress/_tests/test_approximation.py
index ff54939b..5e43380f 100644
--- a/src/napari_stress/_tests/test_approximation.py
+++ b/src/napari_stress/_tests/test_approximation.py
@@ -9,6 +9,10 @@ def test_lsq_ellipsoid():
 
     fitted_points = approximation.expand_points_on_ellipse(ellipsoid, pointcloud)
 
+    assert np.linalg.norm(ellipsoid[:, 1, 1:], axis=1)[0] > np.linalg.norm(ellipsoid[:, 1, 1:], axis=1)[1]
+    assert np.linalg.norm(ellipsoid[:, 1, 1:], axis=1)[1] > np.linalg.norm(ellipsoid[:, 1, 1:], axis=1)[2]
+
+
 def test_lsq_ellipsoid2():
     import vedo
     from napari_stress import approximation
@@ -19,10 +23,10 @@ def test_lsq_ellipsoid2():
 
     from scipy.spatial.transform import Rotation as R
 
-    center = np.zeros((3,3))
-    l1 = 1
+    center = np.zeros((3, 3))
+    l1 = 3
     l2 = 2
-    l3 = 3
+    l3 = 1
     a1 = np.array([l1, 0, 0])
     a2 = np.array([0, l2, 0])
     a3 = np.array([0, 0, l3])
@@ -47,9 +51,9 @@ def test_lsq_ellipsoid2():
     lengths = _axes_lengths_from_ellipsoid(fitted_ellipsoid)
 
     assert np.allclose(center, x0, atol=0.01)
-    assert np.allclose(l2/2, lengths[0])
-    assert np.allclose(l1/2, lengths[1])
-    assert np.allclose(l3/2, lengths[2])
+    assert np.allclose(l1 / 2, lengths[0])
+    assert np.allclose(l2 / 2, lengths[1])
+    assert np.allclose(l3 / 2, lengths[2])
 
     rotation_matrix = R.from_euler('xz', [45, 45], degrees=True).as_matrix()
     rotated_points = np.dot(points, rotation_matrix)
@@ -60,6 +64,7 @@ def test_lsq_ellipsoid2():
 
     assert np.allclose(rotated_points, expanded_points)
 
+
 def test_ellipse_normals():
     from napari_stress import approximation, get_droplet_point_cloud
     pointcloud = get_droplet_point_cloud()[0][0][:, 1:]
@@ -84,4 +89,3 @@ def test_curvature_on_ellipsoid(make_napari_viewer):
     assert types._METADATAKEY_H_E123_ELLIPSOID in results[1]['metadata'].keys()
     assert types._METADATAKEY_H0_ELLIPSOID in results[1]['metadata'].keys()
     assert types._METADATAKEY_MEAN_CURVATURE in results[1]['features'].keys()
-