|
| static void | init_sphere_part_node (struct sphere_part_node *node) |
| |
| static struct sphere_part_node * | get_sphere_part_node () |
| |
| static void | compute_gc_norm_vector (void *coords_data, size_t coords_size, size_t coords_count, double prev_gc_norm_vector[], double gc_norm_vector[]) |
| |
| static void | swap_partition_data (struct temp_partition_data *a, struct temp_partition_data *b) |
| |
| static size_t | swap_node_type (struct temp_partition_data *part_data, size_t i, int node_type, size_t begin, size_t end) |
| |
| static void | sort_partition_data (struct temp_partition_data *part_data, size_t I_FULL_size, size_t I_size, size_t U_size, size_t T_size) |
| |
| static void | partition_data (size_t *local_cell_ids, struct temp_partition_data *part_data, size_t num_cell_ids, size_t threshold, struct sphere_part_node *parent_node, double prev_gc_norm_vector[]) |
| |
| static int | compare_point_idx_xyz (void const *a, void const *b) |
| |
| static struct point_sphere_part_node * | partition_point_data (struct point_id_xyz *points, size_t num_points, size_t threshold, double prev_gc_norm_vector[], size_t curr_tree_depth, size_t *max_tree_depth, int *list_flag) |
| |
| struct bnd_sphere_part_search * | yac_bnd_sphere_part_search_new (struct bounding_circle *circles, size_t num_circles) |
| |
| static int | compare_points_int32_coord (const void *a, const void *b) |
| |
| static struct point_id_xyz * | get_unique_points (size_t *num_points, yac_const_coordinate_pointer coordinates_xyz, yac_int const *ids, int const *mask) |
| |
| struct point_sphere_part_search * | yac_point_sphere_part_search_new (size_t num_points, yac_const_coordinate_pointer coordinates_xyz, yac_int const *ids) |
| |
| struct point_sphere_part_search * | yac_point_sphere_part_search_mask_new (size_t num_points, yac_const_coordinate_pointer coordinates_xyz, yac_int const *ids, int const *mask) |
| |
| static void | search_bnd_circle_I_node (struct sphere_part_node *node, struct bounding_circle bnd_circle, size_t **restrict overlap_cells, size_t *overlap_cells_array_size, size_t *restrict num_overlap_cells, struct overlaps *search_interval_tree_buffer, double prev_gc_norm_vector[]) |
| |
| static void | search_big_bnd_circle (struct sphere_part_node *node, struct bounding_circle bnd_circle, size_t **restrict overlap_cells, size_t *overlap_cells_array_size, size_t *restrict num_overlap_cells, struct overlaps *search_interval_tree_buffer, double prev_gc_norm_vector[]) |
| | TODO change to iterative implementation and allocate overlap_cells first.
|
| |
| static void | search_small_bnd_circle (struct sphere_part_node *node, struct bounding_circle bnd_circle, size_t **restrict overlap_cells, size_t *overlap_cells_array_size, size_t *restrict num_overlap_cells, struct overlaps *search_interval_tree_buffer, double prev_gc_norm_vector[]) |
| |
| static void | search_bnd_circle (struct sphere_part_node *node, struct bounding_circle bnd_circle, size_t **restrict overlap_cells, size_t *overlap_cells_array_size, size_t *restrict num_overlap_cells, struct overlaps *search_interval_tree_buffer, double prev_gc_norm_vector[]) |
| |
| static void | check_leaf_NN (struct point_id_xyz *points, size_t num_points, double *point_coordinates_xyz, struct sin_cos_angle *best_angle, double(**result_coordinates_xyz)[3], size_t *result_coordinates_xyz_array_size, size_t **local_point_ids, size_t *local_point_ids_array_size, size_t total_num_local_point_ids, size_t *num_local_point_ids) |
| |
| static void | point_search_NN (struct bounding_circle *bnd_circle, double(**result_coordinates_xyz)[3], size_t *result_coordinates_xyz_array_size, size_t **local_point_ids, size_t *local_point_ids_array_size, size_t total_num_local_point_ids, size_t *num_local_point_ids, double *dot_stack, struct point_sphere_part_node **node_stack, int *flags, size_t curr_tree_depth) |
| |
| static int | leaf_contains_matching_point (struct point_id_xyz *points, size_t num_points, double coordinate_xyz[3], size_t **local_point_ids, size_t *local_point_ids_array_size, double(**result_coordinates_xyz)[3], size_t *result_coordinates_xyz_array_size, size_t total_num_local_point_ids, size_t *num_local_point_ids) |
| |
| void | yac_point_sphere_part_search_NN (struct point_sphere_part_search *search, size_t num_points, double(*coordinates_xyz)[3], double *cos_angles, double(**result_coordinates_xyz)[3], size_t *result_coordinates_xyz_array_size, size_t **local_point_ids, size_t *local_point_ids_array_size, size_t *num_local_point_ids) |
| |
| static int | compare_point_id_xyz_angle (const void *a, const void *b) |
| |
| static size_t | initial_point_bnd_search_NNN (size_t n, struct point_id_xyz *points, size_t num_points, double *point_coordinates_xyz, struct point_id_xyz_angle **results, size_t *results_array_size) |
| |
| static struct sin_cos_angle | check_leaf_NNN (size_t n, double *point_coordinates_xyz, struct point_id_xyz *points, size_t num_points, struct point_id_xyz_angle **results, size_t *results_array_size, size_t *num_results, struct sin_cos_angle curr_angle) |
| |
| static void | point_search_NNN (size_t n, double *point_coordinates_xyz, struct point_id_xyz_angle **results, size_t *results_array_size, size_t *num_results, double *dot_stack, struct point_sphere_part_node **node_stack, int *flags, size_t curr_tree_depth) |
| |
| void | yac_point_sphere_part_search_NNN (struct point_sphere_part_search *search, size_t num_points, double(*coordinates_xyz)[3], size_t n, double **cos_angles, size_t *cos_angles_array_size, double(**result_coordinates_xyz)[3], size_t *result_coordinates_xyz_array_size, size_t **local_point_ids, size_t *local_point_ids_array_size, size_t *num_local_point_ids) |
| |
| void | yac_point_sphere_part_search_NNN_bnd_circle (struct point_sphere_part_search *search, size_t num_bnd_circles, struct bounding_circle *bnd_circles, size_t n, size_t **local_point_ids, size_t *local_point_ids_array_size, size_t *num_local_point_ids) |
| |
| static int | compare_angles_ (void const *a, void const *b) |
| |
| void | yac_point_sphere_part_search_NNN_ubound (struct point_sphere_part_search *search, size_t num_points, yac_coordinate_pointer coordinates_xyz, size_t n, struct sin_cos_angle *angles) |
| |
| static void | search_bnd_circle_points (struct point_sphere_part_node const *node, struct point_id_xyz const *points, struct bounding_circle const *bnd_circle, size_t **local_point_ids, size_t *local_point_ids_array_size, size_t *num_local_point_ids) |
| |
| void | yac_point_sphere_part_search_bnd_circle (struct point_sphere_part_search *search, size_t num_bnd_circles, const_bounding_circle_pointer bnd_circles, size_t **local_point_ids, size_t *local_point_ids_array_size, size_t *num_local_point_ids) |
| |
| static void | search_point (struct sphere_part_node *node, double point[], size_t **overlap_cells, size_t *overlap_cells_array_size, size_t *num_overlap_cells, struct overlaps *search_interval_tree_buffer, double prev_gc_norm_vector[]) |
| |
| void | yac_bnd_sphere_part_search_do_point_search (struct bnd_sphere_part_search *search, yac_coordinate_pointer coordinates_xyz, size_t count, size_t **cells, size_t *num_cells_per_coordinate) |
| |
| void | yac_bnd_sphere_part_search_do_bnd_circle_search (struct bnd_sphere_part_search *search, struct bounding_circle *bnd_circles, size_t count, size_t **cells, size_t *num_cells_per_bnd_circle) |
| |
| static void | free_sphere_part_tree (struct sphere_part_node tree) |
| |
| static void | free_point_sphere_part_tree (struct point_sphere_part_node *tree) |
| |
| void | yac_delete_point_sphere_part_search (struct point_sphere_part_search *search) |
| |
| void | yac_bnd_sphere_part_search_delete (struct bnd_sphere_part_search *search) |
| |
Implementation of the sphere partitioning algorithm for spatial search.
This file implements the Sphere Partitioning Algorithm; a hierarchical space partitioning algorithm for efficient spatial searches on the unit sphere. The algorithm subdivides the sphere using planes defined through great circles. These subdivisions are stored in a binary tree structure for efficient search.
Key features of the algorithm:
- Recursive binary space partitioning using great circles orthogonal to the great circles used in the previous stage; first great circle is orthogonal to the equator
- Each node partitions data (points on the sphere or bounding circles) into three sets:
- I: Data intersecting the partitioning great circle (stored as interval tree or list depending on size)
- U: Data on the negative side of the great circle (recursively partitioned)
- T: Data on the positive side of the great circle (recursively partitioned)
- Supports both point-based and bounding circle-based searches
- Interval trees used for efficient intersection queries along great circles
The partitioning plane at each level is chosen to:
- Be orthogonal to the previous partitioning plane
- Pass through the centroid of the data This ensures good spatial subdivision and balanced tree depth.
Two Types of Search Structures
This implementation provides two distinct search structures, each optimized for different use cases:
Point Search
The Point Search structure is designed for searching exact point locations. It uses simpler 2-way partitioning:
- U subtree: Points on the negative side of the great circle plane (recursively partitioned if above threshold)
- T subtree: Points on the positive side of the great circle plane (recursively partitioned if above threshold)
Each node stores U_size, T_size, and gc_norm_vector. This simpler structure results in faster tree construction and more efficient memory usage compared to Bounding Circle Search.
Use cases: Nearest neighbor searches (NN, NNN), points within bounding circles, and distance-bounded queries. The structure automatically removes duplicate points during construction.
Bounding Circle Search
The Bounding Circle Search structure is designed for searching bounding circles. It uses 3-way partitioning:
- I list: Bounding circles that intersect the partitioning great circle. Since bounding circles have spatial extent (inclusion angle), they can overlap with the partitioning plane. The I list is implemented as an interval tree for efficient angular range queries when the number of bounding circles exceeds a threshold.
- U subtree: Bounding circles on the negative side of the great circle plane (recursively partitioned if above threshold)
- T subtree: Bounding circles on the positive side of the great circle plane (recursively partitioned if above threshold)
Each node stores I_size, U_size, T_size, an I_angle defining the angular range covering all intersecting bounding circles, and the gc_norm_vector defining the partitioning plane.
Use cases: Finding which bounding circles contain or overlap with query points or bounding circles. Search operations include yac_bnd_sphere_part_search_do_point_search and yac_bnd_sphere_part_search_do_bnd_circle_search.
Why Two Different Structures?
The fundamental distinction arises from the geometric properties of the data:
- Points are dimensionless: They are always strictly on one side of any plane, allowing simpler 2-way partitioning.
- Bounding circles are 2D-shapes on the sphere: The inclusion angle means they can overlap with the partitioning plane, requiring the I list and interval tree machinery.
Using specialized structures for each case provides optimal performance for the respective search operations rather than forcing one approach to handle both cases suboptimally.
Search Query Implementation
Both search structures use tree traversal algorithms that exploit the geometric properties of great circle partitioning to prune the search space efficiently.
Geometric Foundation: The Dot Product Test
The core spatial test used throughout the search algorithms is the dot product between a query point and a node's gc_norm_vector (great circle normal vector):
dot = gc_norm_vector[0] * point[0] +
gc_norm_vector[1] * point[1] +
gc_norm_vector[2] * point[2]
Geometric interpretation: The great circle partitioning plane passes through the sphere's center (origin) with normal vector gc_norm_vector. For any point p on the unit sphere:
dot > 0: Point p lies on the positive side of the plane (same hemisphere as gc_norm_vector points toward). The angle between p and gc_norm_vector is less than PI/2.dot < 0: Point p lies on the negative side of the plane (opposite hemisphere from gc_norm_vector). The angle between p and gc_norm_vector is greater than PI/2.dot = 0: Point p lies nearly on the great circle itself. The angle between p and gc_norm_vector is approximately PI/2, meaning p is perpendicular to the normal vector.
Mathematical basis: For a plane through the origin with unit normal n, (gc_norm_vector) the signed distance of point p from the plane is n * p. On the unit sphere, this equals cos(alpha) where alpha is the angle between n and p. Thus:
cos(alpha) > 0 <==> alpha < PI/2 <==> same hemispherecos(alpha) < 0 <==> alpha > PI/2 <==> opposite hemisphere
Why this matters: The dot product provides an O(1) test to determine which child node ('UorT`) contains a query point, or whether both subtrees must be searched (for range queries). This enables logarithmic-time pruning of the search space.
Bounding Circle Search Queries
The Bounding Circle Search supports two types of queries:
Point-in-cells search (yac_bnd_sphere_part_search_do_point_search):
Finds all cells whose bounding circles contain a given point. The algorithm:
- Computes
dot = gc_norm_vector * query_point at each node (see Geometric Foundation: The Dot Product Test for geometric interpretation)
- Uses the dot product to determine which subtrees to search:
dot > 0 (point on positive side): search T subtreedot < 0 (point on negative side): search U subtree- Always searches I list, as objects intersecting the great circle may contain the point regardless of which side the point is on
- For I list queries with interval trees:
- Projects query point onto the great circle plane
- Computes angular position along the circle relative to reference point (
prev_gc_norm_vector)
- Performs interval tree query to efficiently find bounding circles whose angular ranges overlap the query point's position
- Recursively descends relevant subtrees until all potential cells are found
Bounding-circle-overlap search ( yac_bnd_sphere_part_search_do_bnd_circle_search):
Finds all cells whose bounding circles overlap with a query bounding circle. Distinguishes between two cases:
- Small circles (inc_angle < PI/2): Uses selective traversal. At each node, computes dot = gc_norm_vector * query_circle.base_vector. Then checks:
- Search T if:
dot > -inc_angle.sin (circle could extend to positive side)
- Search U if:
dot < inc_angle.sin (circle could extend to negative side)
- Search I if:
|dot| ≤ sin(I_angle + inc_angle) (circle overlaps great circle) This exploits the inclusion angle to determine spatial overlap possibilities.
- Large circles (
inc_angle ≥ PI/2): Searches entire tree since such circles span more than a hemisphere and can overlap with cells anywhere.
Point Search Queries
The point_sphere_part_search implements nearest neighbor searches with different strategies:
Nearest neighbor search (yac_point_sphere_part_search_NN):
Finds the closest point(s) to each query location.
- Initial descent: Traverses tree to find the leaf containing the query location:
- Computes
dot = gc_norm_vector * query_point at each node
- If
dot ≤ 0: descend to U subtree (point on negative side, angle ≥ PI/2)
- If
dot > 0: descend to T subtree (point on positive side, angle < PI/2)
- This follows the same-side principle until reaching a leaf
- Check for exact match: If query point is identical to any point in the leaf (within yac_angle_tol), return immediately with
cos_angle = 1.0.
- Initial best candidate: Searches the initial leaf for the closest point(s), establishing a search radius (best_angle = distance to best point found).
- Iterative expansion: Uses an explicit stack-based algorithm to traverse back up the tree and explore other branches that might contain closer points:
- At each node level, tests whether U or T subtrees could contain points closer than the current best_angle
- Pruning condition for U: Search U only if
dot < sin(best_angle) (meaning: angle between query and gc_norm_vector is such that points on the negative side could be within best_angle of the query)
- Pruning condition for T: Search T only if
dot > -sin(best_angle) (symmetric condition for positive side)
- Geometric insight: If
dot = cos(alpha) where alpha is angle from gc_norm_vector, then points on one side can only be closer than best_angle if the query point is within (alpha + best_angle) or (alpha - best_angle) of straddling the great circle. The sine term comes from converting angular distances.
- Updates best_angle as closer points are found, progressively tightening the pruning threshold
N-nearest neighbors search (yac_point_sphere_part_search_NNN):
Extends NN to find the n closest points. Key differences:
- Initial descent attempts to find a leaf with at least n points, or the best available leaf if no single leaf has n points
- Maintains a sorted array of current n best candidates (sorted by cos_angle in descending order, i.e., smallest angles first)
- Uses the n-th best distance as the pruning threshold instead of just the single best distance
- Handles ties: includes all points at exactly the n-th distance, so result size may exceed n
- Iterative expansion checks:
dot < sin(angle_to_nth_best) for U, dot > -sin(angle_to_nth_best) for T
Bounded N-nearest neighbors ( yac_point_sphere_part_search_NNN_bnd_circle):
Constrains NNN search to points within a bounding circle, useful for limiting search radius. After finding n nearest neighbors, removes those outside the bounding circle constraint.
For a detailed description of the algorithm, data structures, and search procedures, see Sphere Partitioning Algorithm in sphere_part.h.
Definition in file sphere_part.c.
1.11.0