Add two methods: ConvexHullOffsets and CentroidOffsets class

point_utils/offsetter.py

@@ -1,14 +1,23 @@
 """
-Module for calculating offset vectors for points in a point cloud using different methods.
+Module for calculating offset points for the selected points in a point cloud
+using different methods.
 """
-
+import logging
 from abc import ABC, abstractmethod
 import numpy as np
 from scipy.spatial import KDTree
 from scipy.spatial import ConvexHull
 
 __all__ = ['offset_factory']
 
+logger = logging.getLogger(__name__)
+logger.setLevel(logging.INFO)
+# Create StreamHandler for console output
+handler = logging.StreamHandler()
+handler.setFormatter(logging.Formatter(fmt='[%(name)s] %(message)s'))
+# Add the handler to the logger
+logger.addHandler(handler)
+
 
 class OffsetsInterface(ABC):
     """
@@ -26,11 +35,8 @@ class OffsetsInterface(ABC):
     :ivar new_data_label: Label for the new offset points.
     """
 
-    def __init__(self,
-                 all_coordinates: np.ndarray,
-                 labels: list[str],
-                 data_label_to_offset: str,
-                 offset_magnitude: float,
+    def __init__(self, all_coordinates: np.ndarray, labels: list[str],
+                 data_label_to_offset: str, offset_magnitude: float,
                  new_data_label: str):
         self.all_coordinates = all_coordinates
         self.labels = np.array(labels)
@@ -53,17 +59,41 @@ def add_offset_points(self, **kwargs):
         Add new offset points to the original point cloud.
         """
         offset_vectors = self.get_offset_vecs(**kwargs)
-        new_points_coords = self.all_coordinates[self.point_indices] + offset_vectors
+        new_points_coords = self.all_coordinates[
+            self.point_indices] + offset_vectors
         # Add coordinates of new points to the original point cloud
-        self.all_coordinates = np.concatenate((self.all_coordinates, new_points_coords))
+        self.all_coordinates = np.concatenate(
+            (self.all_coordinates, new_points_coords))
 
         # Extend the labels list with the labels for new points
-        self.labels = np.concatenate((
-            self.labels, [self.new_data_label for _ in range(len(new_points_coords))]))
+        self.labels = np.concatenate(
+            (self.labels,
+             [self.new_data_label for _ in range(len(new_points_coords))]))
 
         assert self.all_coordinates.shape[0] == len(self.labels), \
             "The number of labels and coordinates mismatch."
 
+    @staticmethod
+    def normalizer(vec: np.ndarray) -> np.ndarray:
+        """
+        Normalize the input vector. If the norm of the input vector is zero,
+        a random unit vector is returned.
+
+        :param direction_vec: Numpy array of shape (1, 3) representing a direction vector.
+
+        :return: Normalized direction vector.
+        """
+        norm = np.linalg.norm(vec)
+
+        if np.isclose(norm, 0., atol=1e-18):
+            logger.warning(
+                "The norm of the input vector is zero. Use a random unit vector."
+            )
+            vec = np.random.randn(3)
+            norm = np.linalg.norm(vec)
+
+        return vec / norm
+
 
 class KDTreeOffsets(OffsetsInterface):
 
@@ -80,6 +110,11 @@ def get_offset_vecs(self, num_neighbors: int = 10) -> np.ndarray:
         use the mean of displacement vectors from this point to its nearest neighbors
         (computed by a K-D Tree) to determine the direction offset vector.
 
+        Note that this method is not suitable for highly symmetric point configurations
+        where the neighboring points are uniformly distributed around each point
+        since the mean displacement vector for each point is zero, e.g.
+        an infinite regular cubic lattice in 3D space.
+
         :param num_neighbors: The number of nearest neighbors to include when calculating
                               displacement vectors. Small values will have a localized offset
                               direction.
@@ -99,10 +134,16 @@ def get_offset_vecs(self, num_neighbors: int = 10) -> np.ndarray:
             avg_displacement = np.mean(
                 self.all_coordinates[nearest_neighbor_indices] - point, axis=0)
 
+            norm = np.linalg.norm(avg_displacement)
+            if np.isclose(norm, 0):
+                logger.warning(
+                    f"The mean displacement vector for {num_neighbors} nearest neighbors of point {idx} is zero. "
+                    "Consider increasing the number of neighbors or using a different method."
+                )
+
             # Flip the vector to point away from densely populated space,
             # normalize and scale by the input offset magnitude.
-            offset_vector = -avg_displacement / np.linalg.norm(
-                avg_displacement) * self.offset_magnitude
+            offset_vector = -avg_displacement / norm * self.offset_magnitude
 
             offset_vectors.append(offset_vector)
 
@@ -111,28 +152,80 @@ def get_offset_vecs(self, num_neighbors: int = 10) -> np.ndarray:
 
 class ConvexHullOffsets(OffsetsInterface):
 
-    def __init__(self, all_coordinates: np.ndarray, point_indices: list[int],
-                 offset_magnitude: float):
-        super().__init__(all_coordinates, point_indices, offset_magnitude)
+    def __init__(self, **settings):
+        super().__init__(**settings)
+
+    def name(self):
+        return self.__class__.__name__
 
     def get_offset_vecs(self, **kwargs):
         """
-        Calculates offset vectors for the input points using convex hulls
+        Calculates offset vector for each selected point in the direction of the
+        outward normal of the closest triangle/facet on the convex hull.
         """
-        pass
+        # Compute the convex hull of the point cloud
+        hull = ConvexHull(self.all_coordinates)
+        hull_points = self.all_coordinates[hull.vertices]
+        centroid = np.mean(self.all_coordinates[hull.vertices], axis=0)
+        tree = KDTree(hull_points)
+
+        offset_vectors = []
+        for idx in self.point_indices:
+            point = self.all_coordinates[idx]
+
+            # Find the nearest facet (triangle) of the convex hull
+            _, indices = tree.query(point.reshape(1, -1), k=3)
+            # Get the coordinates of the facet vertices
+            facet_coordinates = self.all_coordinates[indices[0]]
+            # Compute the outward normal vector of the facet
+            direction = np.cross(facet_coordinates[1] - facet_coordinates[0],
+                                 facet_coordinates[2] - facet_coordinates[0])
+            to_centroid = centroid - point
+            if np.dot(direction, to_centroid) > 0:
+                direction = -direction
+
+            # The following method doesn't work well, so commented out
+            # Approximate the outward normal vector of facet by the vector from
+            # the target point to the nearest hull vertex
+            # _, indices = tree.query(point.reshape(1, -1), k=1)
+            # nearest_vertex = hull_points[indices[0]]
+            # direction = point - nearest_vertex
+
+            # Normalize and scale the direction vector by the offset magnitude
+            offset_vectors.append(
+                self.normalizer(direction) * self.offset_magnitude)
+
+        return np.array(offset_vectors)
 
 
 class CentroidOffsets(OffsetsInterface):
 
-    def __init__(self, all_coordinates: np.ndarray, point_indices: list[int],
-                 offset_magnitude: float):
-        super().__init__(all_coordinates, point_indices, offset_magnitude)
+    def __init__(self, **settings):
+        super().__init__(**settings)
+
+    def name(self):
+        return self.__class__.__name__
 
     def get_offset_vecs(self, **kwargs):
         """
-        Calculates offset vectors for the input points using centroids
+        Calculates offset vectors for the specified points that point away from
+        the centroid of the entire point cloud.
         """
-        pass
+        # Compute the centroid of the point cloud
+        centroid = np.mean(self.all_coordinates, axis=0)
+
+        offset_vectors = []
+
+        for idx in self.point_indices:
+            point = self.all_coordinates[idx]
+            # Compute the direction vector from the centroid to the point
+            direction = point - centroid
+
+            # Normalize and scale the direction vector by the offset magnitude
+            offset_vectors.append(
+                self.normalizer(direction) * self.offset_magnitude)
+
+        return np.array(offset_vectors)
 
 
 OFFSET_METHOD_TO_CLASS = {

tests/offsetter_test.py

@@ -4,12 +4,15 @@
 from point_utils import offset_factory, get_data_from_txt
 from point_utils.offsetter import OFFSET_METHOD_TO_CLASS
 
+
 @pytest.fixture
 def data_dir():
     # Return the path to the data directory
     return Path(__file__).parent / "data"
 
-@pytest.mark.parametrize("offset_method", ["KDTreeOffsets"])
+
+@pytest.mark.parametrize(
+    "offset_method", ["KDTreeOffsets", "CentroidOffsets", "ConvexHullOffsets"])
 def test_offset_factory(data_dir, offset_method):
     inp_file = data_dir / "cdd.txt"
     all_coordinates, labels = get_data_from_txt(inp_file)
@@ -32,5 +35,7 @@ def test_offset_factory(data_dir, offset_method):
     offset_calculator.add_offset_points(**kwargs)
 
     # Check if the number of new points is as expected
-    assert offset_calculator.all_coordinates.shape[0] == len(labels) + expected_num_new_points
-    assert offset_calculator.labels.shape[0] == len(labels) + expected_num_new_points
+    assert offset_calculator.all_coordinates.shape[0] == len(
+        labels) + expected_num_new_points
+    assert offset_calculator.labels.shape[0] == len(
+        labels) + expected_num_new_points