Coverage for src/beamme/geometric_search/find_close_points.py: 95%
60 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-28 15:20 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2026 BeamMe 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 tolerance
23of each other will be considered as unique."""
25from enum import Enum as _Enum
26from enum import auto as _auto
28from beamme.geometric_search.scipy import (
29 find_close_points_scipy as _find_close_points_scipy,
30)
31from beamme.geometric_search.utils import arborx_is_available as _arborx_is_available
32from beamme.geometric_search.utils import cython_is_available as _cython_is_available
34if _cython_is_available():
35 from beamme.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 beamme.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 between the
55 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 """
64 unique_indices = []
65 inverse_indices = [-1 for i in range(len(point_partners))]
66 partner_id_to_unique_map = {}
67 i_partner = 0
68 for i_point, partner_index in enumerate(point_partners):
69 if partner_index == -1:
70 # This point does not have any partners, i.e., it is already a unique point.
71 unique_indices.append(i_point)
72 my_inverse_index = len(unique_indices) - 1
73 elif partner_index == i_partner:
74 # This point has partners and this is the first time that this partner index
75 # appears in the input list.
76 unique_indices.append(i_point)
77 my_inverse_index = len(unique_indices) - 1
78 partner_id_to_unique_map[partner_index] = my_inverse_index
79 i_partner += 1
80 elif partner_index < i_partner:
81 # This point has partners and the partner index has been previously found.
82 my_inverse_index = partner_id_to_unique_map[partner_index]
83 else:
84 raise ValueError(
85 "This should not happen, as the partners should be provided in order"
86 )
87 inverse_indices[i_point] = my_inverse_index
89 return unique_indices, inverse_indices
92def point_partners_to_partner_indices(point_partners, n_partners):
93 """Convert the partner indices for each point to a list of lists with the indices
94 for all partners."""
95 partner_indices = [[] for i in range(n_partners)]
96 for i, partner_index in enumerate(point_partners):
97 if partner_index != -1:
98 partner_indices[partner_index].append(i)
99 return partner_indices
102def partner_indices_to_point_partners(partner_indices, n_points):
103 """Convert the list of lists with the indices for all partners to the partner
104 indices for each point."""
105 point_partners = [-1 for _i in range(n_points)]
106 for i_partner, partners in enumerate(partner_indices):
107 for index in partners:
108 point_partners[index] = i_partner
109 return point_partners, len(partner_indices)
112def find_close_points(point_coordinates, *, algorithm=None, tol=1e-8, **kwargs):
113 """Find unique points in a point cloud, i.e., points that are within a certain
114 tolerance of each other will be considered as unique.
116 Args
117 ----
118 point_coordinates: _np.array(n_points x n_dim)
119 Point coordinates that are checked for partners. The number of spatial dimensions
120 does not have to be equal to 3.
121 algorithm: FindClosePointAlgorithm
122 Type of geometric search algorithm that should be used.
123 n_bins: list(int)
124 Number of bins in the first three dimensions.
125 tol: float
126 If the absolute distance between two points is smaller than tol, they
127 are considered to be equal, i.e., tol is the hyper sphere radius that
128 the point coordinates have to be within, to be identified as overlapping.
130 Return
131 ----
132 has_partner: array(int)
133 An array with integers, marking the partner index of each point. A partner
134 index of -1 means the node does not have a partner.
135 partner: int
136 Largest partner index.
137 """
138 n_points = len(point_coordinates)
140 if algorithm is None:
141 # Decide which algorithm to use
142 if n_points < 200 and _cython_is_available():
143 # For around 200 points the brute force cython algorithm is the fastest one
144 algorithm = FindClosePointAlgorithm.brute_force_cython
145 elif _arborx_is_available():
146 # For general problems with n_points > 200 the ArborX implementation is the fastest one
147 algorithm = FindClosePointAlgorithm.boundary_volume_hierarchy_arborx
148 else:
149 # The scipy implementation is slower than ArborX by a factor of about 2, but is scales
150 # the same
151 algorithm = FindClosePointAlgorithm.kd_tree_scipy
153 # Get list of closest pairs
154 if algorithm is FindClosePointAlgorithm.kd_tree_scipy:
155 has_partner, n_partner = _find_close_points_scipy(
156 point_coordinates, tol, **kwargs
157 )
158 elif algorithm is FindClosePointAlgorithm.brute_force_cython:
159 has_partner, n_partner = _find_close_points_brute_force_cython(
160 point_coordinates, tol, **kwargs
161 )
162 elif algorithm is FindClosePointAlgorithm.boundary_volume_hierarchy_arborx:
163 has_partner, n_partner = _find_close_points_arborx(
164 point_coordinates, tol, **kwargs
165 )
166 else:
167 raise TypeError(f"Got unexpected algorithm {algorithm}")
169 return has_partner, n_partner