From 923fb350c5ae7a62831ce1b959164e3cca7cee89 Mon Sep 17 00:00:00 2001 From: Renke Huang Date: Mon, 18 Nov 2024 14:41:21 -0500 Subject: [PATCH 1/2] Add two methods: "ConvexHullOffsets" and "CentroidOffsets" class --- point_utils/offsetter.py | 130 ++++++++++++++++++++++++++++++++------- tests/offsetter_test.py | 11 +++- 2 files changed, 115 insertions(+), 26 deletions(-) diff --git a/point_utils/offsetter.py b/point_utils/offsetter.py index 1dda110..4d82888 100644 --- a/point_utils/offsetter.py +++ b/point_utils/offsetter.py @@ -1,7 +1,8 @@ """ -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 @@ -9,6 +10,14 @@ __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,71 @@ 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 convex hull. + """ - pass + # Compute the convex hull of the point cloud + hull = ConvexHull(self.all_coordinates) + hull_points = self.all_coordinates[hull.vertices] + tree = KDTree(hull_points) + + offset_vectors = [] + for idx in self.point_indices: + point = self.all_coordinates[idx] + + # Find the closest facet (triangle) of the convex hull + _, indices = tree.query(point.reshape(1, -1), k=1) + print(idx, indices) + + # Approximate the outward normal vector of facet by the vector from + # the target point to the nearest hull vertex + 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 = { diff --git a/tests/offsetter_test.py b/tests/offsetter_test.py index dd5461f..838dfb9 100644 --- a/tests/offsetter_test.py +++ b/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 \ No newline at end of file + 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 From d62530e8af6648d5e1b229401578ea551b166036 Mon Sep 17 00:00:00 2001 From: Renke Huang Date: Mon, 18 Nov 2024 15:02:39 -0500 Subject: [PATCH 2/2] Use exact method to calculate the outward normal vector of facet --- point_utils/offsetter.py | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/point_utils/offsetter.py b/point_utils/offsetter.py index 4d82888..5a0555b 100644 --- a/point_utils/offsetter.py +++ b/point_utils/offsetter.py @@ -161,26 +161,35 @@ def name(self): def get_offset_vecs(self, **kwargs): """ Calculates offset vector for each selected point in the direction of the - outward normal of the convex hull. - + outward normal of the closest triangle/facet on the convex hull. """ # 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 closest facet (triangle) of the convex hull - _, indices = tree.query(point.reshape(1, -1), k=1) - print(idx, indices) - + # 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 - nearest_vertex = hull_points[indices[0]] - direction = point - nearest_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(