Coverage for src/meshpy/geometric_search/find_close_points.py: 95%
60 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-28 04:21 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2025 MeshPy Authors
4#
5# Permission is hereby granted, free of charge, to any person obtaining a copy
6# of this software and associated documentation files (the "Software"), to deal
7# in the Software without restriction, including without limitation the rights
8# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9# copies of the Software, and to permit persons to whom the Software is
10# furnished to do so, subject to the following conditions:
11#
12# The above copyright notice and this permission notice shall be included in
13# all copies or substantial portions of the Software.
14#
15# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
21# THE SOFTWARE.
22"""Find unique points in a point cloud, i.e., points that are within a certain
23tolerance of each other will be considered as unique."""
25from enum import Enum as _Enum
26from enum import auto as _auto
28from meshpy.geometric_search.scipy import (
29 find_close_points_scipy as _find_close_points_scipy,
30)
31from meshpy.geometric_search.utils import arborx_is_available as _arborx_is_available
32from meshpy.geometric_search.utils import cython_is_available as _cython_is_available
34if _cython_is_available():
35 from meshpy.geometric_search.cython import (
36 find_close_points_brute_force_cython as _find_close_points_brute_force_cython,
37 )
39if _arborx_is_available():
40 from meshpy.geometric_search.arborx import (
41 find_close_points_arborx as _find_close_points_arborx,
42 )
45class FindClosePointAlgorithm(_Enum):
46 """Enum for different find_close_point algorithms."""
48 kd_tree_scipy = _auto()
49 brute_force_cython = _auto()
50 boundary_volume_hierarchy_arborx = _auto()
53def point_partners_to_unique_indices(point_partners):
54 """Convert the partner indices to lists that can be used for converting
55 between the full and unique coordinates.
57 Returns
58 ----
59 unique_indices: list(int)
60 Indices that result in the unique point coordinate array.
61 inverse_indices: list(int)
62 Indices of the unique array that can be used to reconstruct of the original points coordinates.
63 """
65 unique_indices = []
66 inverse_indices = [-1 for i in range(len(point_partners))]
67 partner_id_to_unique_map = {}
68 i_partner = 0
69 for i_point, partner_index in enumerate(point_partners):
70 if partner_index == -1:
71 # This point does not have any partners, i.e., it is already a unique point.
72 unique_indices.append(i_point)
73 my_inverse_index = len(unique_indices) - 1
74 elif partner_index == i_partner:
75 # This point has partners and this is the first time that this partner index
76 # appears in the input list.
77 unique_indices.append(i_point)
78 my_inverse_index = len(unique_indices) - 1
79 partner_id_to_unique_map[partner_index] = my_inverse_index
80 i_partner += 1
81 elif partner_index < i_partner:
82 # This point has partners and the partner index has been previously found.
83 my_inverse_index = partner_id_to_unique_map[partner_index]
84 else:
85 raise ValueError(
86 "This should not happen, as the partners should be provided in order"
87 )
88 inverse_indices[i_point] = my_inverse_index
90 return unique_indices, inverse_indices
93def point_partners_to_partner_indices(point_partners, n_partners):
94 """Convert the partner indices for each point to a list of lists with the
95 indices for all partners."""
96 partner_indices = [[] for i in range(n_partners)]
97 for i, partner_index in enumerate(point_partners):
98 if partner_index != -1:
99 partner_indices[partner_index].append(i)
100 return partner_indices
103def partner_indices_to_point_partners(partner_indices, n_points):
104 """Convert the list of lists with the indices for all partners to the
105 partner indices for each point."""
106 point_partners = [-1 for _i in range(n_points)]
107 for i_partner, partners in enumerate(partner_indices):
108 for index in partners:
109 point_partners[index] = i_partner
110 return point_partners, len(partner_indices)
113def find_close_points(point_coordinates, *, algorithm=None, tol=1e-8, **kwargs):
114 """Find unique points in a point cloud, i.e., points that are within a
115 certain tolerance of each other will be considered as unique.
117 Args
118 ----
119 point_coordinates: _np.array(n_points x n_dim)
120 Point coordinates that are checked for partners. The number of spatial dimensions
121 does not have to be equal to 3.
122 algorithm: FindClosePointAlgorithm
123 Type of geometric search algorithm that should be used.
124 n_bins: list(int)
125 Number of bins in the first three dimensions.
126 tol: float
127 If the absolute distance between two points is smaller than tol, they
128 are considered to be equal, i.e., tol is the hyper sphere radius that
129 the point coordinates have to be within, to be identified as overlapping.
131 Return
132 ----
133 has_partner: array(int)
134 An array with integers, marking the partner index of each point. A partner
135 index of -1 means the node does not have a partner.
136 partner: int
137 Largest partner index.
138 """
140 n_points = len(point_coordinates)
142 if algorithm is None:
143 # Decide which algorithm to use
144 if n_points < 200 and _cython_is_available():
145 # For around 200 points the brute force cython algorithm is the fastest one
146 algorithm = FindClosePointAlgorithm.brute_force_cython
147 elif _arborx_is_available():
148 # For general problems with n_points > 200 the ArborX implementation is the fastest one
149 algorithm = FindClosePointAlgorithm.boundary_volume_hierarchy_arborx
150 else:
151 # The scipy implementation is slower than ArborX by a factor of about 2, but is scales
152 # the same
153 algorithm = FindClosePointAlgorithm.kd_tree_scipy
155 # Get list of closest pairs
156 if algorithm is FindClosePointAlgorithm.kd_tree_scipy:
157 has_partner, n_partner = _find_close_points_scipy(
158 point_coordinates, tol, **kwargs
159 )
160 elif algorithm is FindClosePointAlgorithm.brute_force_cython:
161 has_partner, n_partner = _find_close_points_brute_force_cython(
162 point_coordinates, tol, **kwargs
163 )
164 elif algorithm is FindClosePointAlgorithm.boundary_volume_hierarchy_arborx:
165 has_partner, n_partner = _find_close_points_arborx(
166 point_coordinates, tol, **kwargs
167 )
168 else:
169 raise TypeError(f"Got unexpected algorithm {algorithm}")
171 return has_partner, n_partner