fix rotation / spin issue

ossdbs/electrodes/boston_scientific_vercise.py

@@ -4,6 +4,7 @@
 # SPDX-License-Identifier: GPL-3.0-or-later
 
 # Boston Scientific (Marlborough, Massachusetts, USA) vercise
+import logging
 from dataclasses import dataclass
 
 import netgen
@@ -12,7 +13,9 @@
 from scipy.spatial.transform import Rotation
 
 from .electrode_model_template import ElectrodeModel
-from .utilities import get_highest_edge, get_lowest_edge
+from .utilities import get_highest_edge, get_lowest_edge, get_signed_angle
+
+_logger = logging.getLogger(__name__)
 
 
 @dataclass
@@ -168,19 +171,22 @@ def _contacts(self) -> netgen.libngpy._NgOCC.TopoDS_Shape:
             desired_direction = (0, self._direction[2], -self._direction[1])
             rotate_vector = Rotation.from_rotvec(np.radians(angle) * np.array(rotation))
             current_direction = rotate_vector.apply((0, 1, 0))
-            print(desired_direction)
-            print(current_direction)
-            print(np.linalg.norm(desired_direction))
-            print(np.linalg.norm(current_direction))
-            rotation_angle = np.degrees(
-                np.arccos(np.dot(desired_direction, current_direction))
+            # get angle between current and desired direction
+            # current direction is normal
+            rotation_angle = get_signed_angle(
+                current_direction, desired_direction, current_direction
             )
-            print(rotation_angle)
+            if rotation_angle is None:
+                _logger.warning(
+                    "Could not determine rotation angle for "
+                    "correct spin as per Lead-DBS convention."
+                )
+                # to return unrotated geo
+                rotation_angle = 0.0
             if np.isclose(rotation_angle, 0):
                 return rotated_geo
-            print("Return rotated geo")
             return rotated_geo.Rotate(
-                occ.Axis(p=(0, 0, 0), d=self._direction), -rotation_angle
+                occ.Axis(p=(0, 0, 0), d=self._direction), rotation_angle
             )
 
     # ruff: noqa: C901

ossdbs/electrodes/utilities.py

@@ -1,7 +1,10 @@
 # Copyright 2023, 2024 Julius Zimmermann
 # SPDX-License-Identifier: GPL-3.0-or-later
 
+from typing import Union
+
 import netgen.occ as occ
+import numpy as np
 
 
 def get_lowest_edge(contact: occ.Face) -> occ.Edge:
@@ -22,3 +25,35 @@ def get_highest_edge(contact: occ.Face) -> occ.Edge:
             max_edge_val = edge.center.z
             max_edge = edge
     return max_edge
+
+
+def get_signed_angle(
+    v1: np.ndarray, v2: np.ndarray, vn: np.ndarray
+) -> Union[None, float]:
+    """Get signed angle between two vectors.
+
+    Parameters
+    ----------
+    v1: np.ndarray
+        First vector
+    v2: np.ndarray
+        Second vector
+    vn: np.ndarray
+        Normal vector in common plane
+
+    Notes
+    -----
+    See more details on StackOverflow:
+    https://stackoverflow.com/questions/5188561/signed-angle-between-two-3d-vectors-with-same-origin-within-the-same-plane
+    """
+    len_v1 = np.linalg.norm(v1)
+    len_v2 = np.linalg.norm(v2)
+    # catch zero-length
+    if np.isclose(len_v1, 0.0) or np.isclose(len_v2, 0.0):
+        return None
+
+    rotation_angle = np.degrees(np.arccos(np.dot(v1 / len_v1, v2 / len_v2)))
+    cross = np.cross(v1, v2)
+    if np.dot(vn, cross) < 0:
+        return -rotation_angle
+    return rotation_angle