137 void * coords_data,
size_t coords_size,
size_t coords_count,
141 double balance_point[3] = {0.0,0.0,0.0};
144 for (
size_t i = 0; i < coords_count; ++i) {
147 (
double*)(((
unsigned char*)coords_data) + i * coords_size);
148 balance_point[0] += coords[0];
149 balance_point[1] += coords[1];
150 balance_point[2] += coords[2];
153 if ((fabs(balance_point[0]) > 1e-9) ||
154 (fabs(balance_point[1]) > 1e-9) ||
155 (fabs(balance_point[2]) > 1e-9)) {
158 balance_point[0] = prev_gc_norm_vector[2];
159 balance_point[1] = prev_gc_norm_vector[0];
160 balance_point[2] = prev_gc_norm_vector[1];
186 for (
size_t j = begin; j < end; ++j)
200 size_t I_begin = I_FULL_size;
201 size_t U_begin = I_FULL_size + I_size;
202 size_t T_begin = I_FULL_size + I_size + U_size;
204 size_t I_end = I_begin + I_size;
205 size_t U_end = U_begin + U_size;
206 size_t T_end = T_begin + T_size;
208 while (I_end > I_begin && part_data[I_end-1].
node_type ==
I_NODE) I_end--;
209 while (U_end > U_begin && part_data[U_end-1].
node_type ==
U_NODE) U_end--;
210 while (T_end > T_begin && part_data[T_end-1].
node_type ==
T_NODE) T_end--;
212 for (
size_t i = 0; i < I_FULL_size; ++i)
220 I_begin = I_FULL_size;
221 for (
size_t i = I_begin; i < I_end; ++i)
233 U_begin = I_FULL_size + I_size;
234 T_begin = I_FULL_size + I_size + U_size;
235 for (
size_t i = U_begin; i < U_end; ++i)
246 size_t num_cell_ids,
size_t threshold,
struct sphere_part_node * parent_node,
247 double prev_gc_norm_vector[]) {
249 if (num_cell_ids == 0) {
255 part_data,
sizeof(*part_data), num_cell_ids,
261 size_t I_FULL_size = 0;
268 for (
size_t i = 0; i < num_cell_ids; ++i) {
306 max_inc_angle = inc_angle;
310 }
else if (angle.
cos < 0.0) {
326 I_size += I_FULL_size;
327 parent_node->
I_size = I_size;
328 parent_node->
U_size = U_size;
329 parent_node->
T_size = T_size;
335 parent_node->
I_angle = max_inc_angle;
349 for (
size_t i = 0; i < I_FULL_size; ++i) {
356 for (
size_t i = I_FULL_size; i < I_size; ++i) {
358 double GCp[3], bVp[3];
367 (fabs(bVp[0]) > 1e-9) || (fabs(bVp[1]) > 1e-9) ||
368 (fabs(bVp[2]) > 1e-9),
369 "ERROR(partition_data): "
370 "projected vector is nearly identical to gc_norm_vector")
381 double bnd_circle_lat_cos =
384 M_PI_2 - acos(curr_bnd_circle.
inc_angle.
sin/bnd_circle_lat_cos);
401 for (
size_t i = 0; i < I_size; ++i)
402 local_cell_ids[i] = part_data[i].local_id;
403 parent_node->
I.
list = (
void*)local_cell_ids;
406 parent_node->
I.
list = NULL;
409 local_cell_ids += I_size;
412 if (U_size <= threshold) {
414 for (
size_t i = 0; i < U_size; ++i)
415 local_cell_ids[i] = part_data[i].local_id;
416 parent_node->
U = (
void*)local_cell_ids;
425 local_cell_ids += U_size;
428 if (T_size <= threshold) {
430 for (
size_t i = 0; i < T_size; ++i)
431 local_cell_ids[i] = part_data[i].local_id;
432 parent_node->
T = (
void*)local_cell_ids;
434 local_cell_ids += T_size;
451 double prev_gc_norm_vector[],
size_t curr_tree_depth,
452 size_t * max_tree_depth) {
454 if (curr_tree_depth > *max_tree_depth) *max_tree_depth = curr_tree_depth;
474 while (left <= right) {
476 double dot = curr_coordinates_xyz[0] * gc_norm_vector[0] +
477 curr_coordinates_xyz[1] * gc_norm_vector[1] +
478 curr_coordinates_xyz[2] * gc_norm_vector[2];
481 if (dot > 0.0)
break;
486 while (left < right) {
487 double * curr_coordinates_xyz = &(right->coordinates_xyz[0]);
488 double dot = curr_coordinates_xyz[0] * gc_norm_vector[0] +
489 curr_coordinates_xyz[1] * gc_norm_vector[1] +
490 curr_coordinates_xyz[2] * gc_norm_vector[2];
493 if (dot <= 0.0)
break;
508 size_t U_size = left -
points;
516 if ((U_size <= threshold) || (U_size ==
num_points)) {
525 curr_tree_depth + 1, max_tree_depth);
528 if ((T_size <= threshold) || (T_size ==
num_points)) {
538 curr_tree_depth + 1, max_tree_depth);
549 double gc_norm_vector[3] = {0.0,0.0,1.0};
557 xmalloc(num_circles *
sizeof(*part_data));
558 for (
size_t i = 0; i < num_circles; ++i) {
572 const void * a,
const void * b) {
579 for (
int i = 0; i < 3; ++i)
589 yac_int const * ids,
int const * mask) {
594 double const scale = (double)(2 << 21);
596 size_t num_unmasked_points;
600 for (
size_t i = 0; i < num_unmasked_points; ++i) {
602 points_int32[i].
idx = i;
603 for (
size_t j = 0; j < 3; ++j)
605 (int32_t)round(coordinates_xyz[i][j] * scale);
608 num_unmasked_points = 0;
611 if (!mask[i])
continue;
612 points_int32[num_unmasked_points].
idx = i;
613 for (
size_t j = 0; j < 3; ++j)
615 (int32_t)round(coordinates_xyz[i][j] * scale);
616 num_unmasked_points++;
621 qsort(points_int32, num_unmasked_points,
625 dummy.
idx = SIZE_MAX;
631 size_t new_num_points = 0;
632 for (
size_t i = 0; i < num_unmasked_points; ++i, ++curr) {
634 size_t curr_idx = curr->
idx;
636 prev = points_int32 + new_num_points++;
637 if (prev != curr) *prev = *curr;
638 prev_id = ids[curr_idx];
640 yac_int curr_id = ids[curr_idx];
641 if (curr_id > prev_id) {
643 prev->
idx = curr_idx;
649 for (
size_t i = 0; i < new_num_points; ++i) {
650 size_t curr_idx = points_int32[i].
idx;
671 size_t max_tree_depth = 0;
688 yac_int const * ids,
int const * mask) {
704 (
double[3]){0.0,0.0,1.0}, 1, &max_tree_depth);
715 size_t ** restrict overlap_cells,
size_t * overlap_cells_array_size,
716 size_t * restrict num_overlap_cells,
717 struct overlaps * search_interval_tree_buffer,
double prev_gc_norm_vector[]) {
724 angle.
cos = fabs(angle.
cos);
731 (*overlap_cells)[(*num_overlap_cells)+i] =
737 double GCp[3], bVp[3];
744 (fabs(bVp[0]) > 1e-9) || (fabs(bVp[1]) > 1e-9) || (fabs(bVp[2]) > 1e-9),
745 "ERROR(search_bnd_circle_I_node): "
746 "projected vector is nearly identical to gc_norm_vector")
752 double bnd_circle_lat_cos =
755 M_PI_2 - acos(bnd_circle.
inc_angle.
sin/bnd_circle_lat_cos);
765 .left = base_angle - inc_angle,
766 .right = base_angle + inc_angle},
767 search_interval_tree_buffer);
773 for (
size_t i = 0; i < search_interval_tree_buffer->
num_overlaps; ++i)
774 (*overlap_cells)[(*num_overlap_cells)+i] =
778 *num_overlap_cells += search_interval_tree_buffer->
num_overlaps;
783 *num_overlap_cells + node->
I_size);
784 memcpy(*overlap_cells + *num_overlap_cells, node->
I.
list,
785 node->
I_size *
sizeof(**overlap_cells));
786 *num_overlap_cells += node->
I_size;
793 size_t ** restrict overlap_cells,
size_t * overlap_cells_array_size,
794 size_t * restrict num_overlap_cells,
795 struct overlaps * search_interval_tree_buffer,
double prev_gc_norm_vector[]) {
800 *num_overlap_cells + node->
T_size);
801 memcpy(*overlap_cells + *num_overlap_cells, node->
T,
802 node->
T_size *
sizeof(**overlap_cells));
803 *num_overlap_cells += node->
T_size;
807 node->
T, bnd_circle, overlap_cells, overlap_cells_array_size,
808 num_overlap_cells, search_interval_tree_buffer, node->
gc_norm_vector);
814 *num_overlap_cells + node->
U_size);
815 memcpy(*overlap_cells + *num_overlap_cells, node->
U,
816 node->
U_size *
sizeof(**overlap_cells));
817 *num_overlap_cells += node->
U_size;
821 node->
U, bnd_circle, overlap_cells, overlap_cells_array_size,
822 num_overlap_cells, search_interval_tree_buffer, node->
gc_norm_vector);
826 node, bnd_circle, overlap_cells, overlap_cells_array_size,
827 num_overlap_cells, search_interval_tree_buffer, prev_gc_norm_vector);
832 size_t ** restrict overlap_cells,
size_t * overlap_cells_array_size,
833 size_t * restrict num_overlap_cells,
834 struct overlaps * search_interval_tree_buffer,
double prev_gc_norm_vector[]) {
846 *num_overlap_cells + node->
T_size);
847 memcpy(*overlap_cells + *num_overlap_cells, node->
T,
848 node->
T_size *
sizeof(**overlap_cells));
849 *num_overlap_cells += node->
T_size;
853 node->
T, bnd_circle, overlap_cells, overlap_cells_array_size,
854 num_overlap_cells, search_interval_tree_buffer,
865 *num_overlap_cells + node->
U_size);
866 memcpy(*overlap_cells + *num_overlap_cells, node->
U,
867 node->
U_size *
sizeof(**overlap_cells));
868 *num_overlap_cells += node->
U_size;
872 node->
U, bnd_circle, overlap_cells, overlap_cells_array_size,
873 num_overlap_cells, search_interval_tree_buffer,
898 if (((angle_sum.
sin < 0.0) || (angle_sum.
cos <= 0.0)) ||
899 (fabs(dot) <= angle_sum.
sin)) {
901 node, bnd_circle, overlap_cells, overlap_cells_array_size,
902 num_overlap_cells, search_interval_tree_buffer, prev_gc_norm_vector);
908 size_t ** restrict overlap_cells,
909 size_t * overlap_cells_array_size,
910 size_t * restrict num_overlap_cells,
911 struct overlaps * search_interval_tree_buffer,
912 double prev_gc_norm_vector[]) {
917 node, bnd_circle, overlap_cells, overlap_cells_array_size,
918 num_overlap_cells, search_interval_tree_buffer, prev_gc_norm_vector);
921 node, bnd_circle, overlap_cells, overlap_cells_array_size,
922 num_overlap_cells, search_interval_tree_buffer, prev_gc_norm_vector);
927 double * point_coordinates_xyz,
struct sin_cos_angle * best_angle,
928 double (**result_coordinates_xyz)[3],
929 size_t * result_coordinates_xyz_array_size,
size_t ** local_point_ids,
930 size_t * local_point_ids_array_size,
size_t total_num_local_point_ids,
931 size_t * num_local_point_ids) {
933 size_t * local_point_ids_ = *local_point_ids;
934 size_t local_point_ids_array_size_ = *local_point_ids_array_size;
935 double (*result_coordinates_xyz_)[3];
936 size_t result_coordinates_xyz_array_size_;
937 size_t num_local_point_ids_ = *num_local_point_ids;
939 if (result_coordinates_xyz != NULL) {
940 result_coordinates_xyz_ = *result_coordinates_xyz;
941 result_coordinates_xyz_array_size_ = *result_coordinates_xyz_array_size;
943 result_coordinates_xyz_, result_coordinates_xyz_array_size_,
944 total_num_local_point_ids + num_local_point_ids_ +
num_points);
945 *result_coordinates_xyz = result_coordinates_xyz_;
946 *result_coordinates_xyz_array_size = result_coordinates_xyz_array_size_;
947 result_coordinates_xyz_ += total_num_local_point_ids;
950 local_point_ids_, local_point_ids_array_size_,
951 total_num_local_point_ids + num_local_point_ids_ +
num_points);
952 *local_point_ids = local_point_ids_;
953 *local_point_ids_array_size = local_point_ids_array_size_;
954 local_point_ids_ += total_num_local_point_ids;
962 points[i].coordinates_xyz, point_coordinates_xyz);
966 if (compare > 0)
continue;
971 *best_angle = curr_angle;
972 num_local_point_ids_ = 1;
973 if (result_coordinates_xyz != NULL) {
974 result_coordinates_xyz_[0][0] =
points[i].coordinates_xyz[0];
975 result_coordinates_xyz_[0][1] =
points[i].coordinates_xyz[1];
976 result_coordinates_xyz_[0][2] =
points[i].coordinates_xyz[2];
978 local_point_ids_[0] =
points[i].idx;
982 if (result_coordinates_xyz != NULL) {
983 result_coordinates_xyz_[num_local_point_ids_][0] =
984 points[i].coordinates_xyz[0];
985 result_coordinates_xyz_[num_local_point_ids_][1] =
986 points[i].coordinates_xyz[1];
987 result_coordinates_xyz_[num_local_point_ids_][2] =
988 points[i].coordinates_xyz[2];
990 local_point_ids_[num_local_point_ids_] =
points[i].idx;
991 num_local_point_ids_++;
995 *num_local_point_ids = num_local_point_ids_;
999 struct bounding_circle * bnd_circle,
double (**result_coordinates_xyz)[3],
1000 size_t * result_coordinates_xyz_array_size,
size_t ** local_point_ids,
1001 size_t * local_point_ids_array_size,
size_t total_num_local_point_ids,
1002 size_t * num_local_point_ids,
double * dot_stack,
1004 int * flags,
size_t curr_tree_depth) {
1006 double * point_coordinates_xyz = bnd_circle->
base_vector;
1009 double dot = dot_stack[curr_tree_depth];
1021 if ((dot < best_angle.
sin) | (best_angle.
cos <= 0.0)) {
1027 point_coordinates_xyz, &best_angle, result_coordinates_xyz,
1028 result_coordinates_xyz_array_size, local_point_ids,
1029 local_point_ids_array_size, total_num_local_point_ids,
1030 num_local_point_ids);
1040 dot_stack[curr_tree_depth] = dot;
1041 node_stack[curr_tree_depth] = node;
1042 flags[curr_tree_depth] = 0;
1055 if ((dot > - best_angle.
sin) || (best_angle.
cos <= 0.0)) {
1061 point_coordinates_xyz, &best_angle, result_coordinates_xyz,
1062 result_coordinates_xyz_array_size, local_point_ids,
1063 local_point_ids_array_size, total_num_local_point_ids,
1064 num_local_point_ids);
1074 dot_stack[curr_tree_depth] = dot;
1075 node_stack[curr_tree_depth] = node;
1076 flags[curr_tree_depth] = 0;
1084 if (curr_tree_depth == 0)
break;
1089 dot = dot_stack[curr_tree_depth];
1090 node = node_stack[curr_tree_depth];
1102 size_t ** local_point_ids,
size_t * local_point_ids_array_size,
1103 double (**result_coordinates_xyz)[3],
1104 size_t * result_coordinates_xyz_array_size,
1105 size_t total_num_local_point_ids,
size_t * num_local_point_ids) {
1113 total_num_local_point_ids + 1);
1114 (*local_point_ids)[total_num_local_point_ids] =
points[i].idx;
1115 if (result_coordinates_xyz != NULL) {
1117 *result_coordinates_xyz_array_size,
1118 total_num_local_point_ids + 1);
1119 memcpy((*result_coordinates_xyz) + total_num_local_point_ids,
1120 points[i].coordinates_xyz, 3 *
sizeof(
double));
1122 *num_local_point_ids = 1;
1132 double (*coordinates_xyz)[3],
1133 double * cos_angles,
1134 double (**result_coordinates_xyz)[3],
1135 size_t * result_coordinates_xyz_array_size,
1136 size_t ** local_point_ids,
1137 size_t * local_point_ids_array_size,
1138 size_t * num_local_point_ids) {
1140 memset(num_local_point_ids, 0,
num_points *
sizeof(*num_local_point_ids));
1142 if (search == NULL)
return;
1146 size_t total_num_local_point_ids = 0;
1157 double * curr_coordinates_xyz = coordinates_xyz[i];
1159 size_t curr_tree_depth = 0;
1161 size_t curr_num_points = 0;
1166 double dot = curr_node->
gc_norm_vector[0]*curr_coordinates_xyz[0] +
1170 dot_stack[curr_tree_depth] = dot;
1171 node_stack[curr_tree_depth] = curr_node;
1172 flags[curr_tree_depth] = 0;
1177 flags[curr_tree_depth] =
U_FLAG;
1180 if (curr_node->
U_size > 0) {
1182 curr_num_points = curr_node->
U_size;
1185 flags[curr_tree_depth] =
T_FLAG;
1188 "ERROR(yac_point_sphere_part_search_NN): "
1189 "if one branch is empty, the other has to be a leaf");
1191 curr_num_points = curr_node->
T_size;
1194 }
else curr_node = curr_node->
U;
1199 flags[curr_tree_depth] =
T_FLAG;
1202 if (curr_node->
T_size > 0) {
1204 curr_num_points = curr_node->
T_size;
1207 flags[curr_tree_depth] =
U_FLAG;
1210 "ERROR(yac_point_sphere_part_search_NN): "
1211 "if one branch is empty, the other has to be a leaf");
1213 curr_num_points = curr_node->
U_size;
1216 }
else curr_node = curr_node->
T;
1224 points, curr_num_points, curr_coordinates_xyz, local_point_ids,
1225 local_point_ids_array_size, result_coordinates_xyz,
1226 result_coordinates_xyz_array_size, total_num_local_point_ids,
1227 num_local_point_ids + i)) {
1229 if (cos_angles != NULL) cos_angles[i] = 1.0;
1234 bnd_circle.
base_vector[0] = curr_coordinates_xyz[0];
1235 bnd_circle.
base_vector[1] = curr_coordinates_xyz[1];
1236 bnd_circle.
base_vector[2] = curr_coordinates_xyz[2];
1238 bnd_circle.
sq_crd = DBL_MAX;
1242 result_coordinates_xyz, result_coordinates_xyz_array_size,
1243 local_point_ids, local_point_ids_array_size, total_num_local_point_ids,
1244 num_local_point_ids + i);
1248 &bnd_circle, result_coordinates_xyz, result_coordinates_xyz_array_size,
1249 local_point_ids, local_point_ids_array_size, total_num_local_point_ids,
1250 num_local_point_ids + i, dot_stack, node_stack, flags, curr_tree_depth);
1252 if (cos_angles != NULL) cos_angles[i] = bnd_circle.
inc_angle.
cos;
1255 total_num_local_point_ids += num_local_point_ids[i];
1271 if (ret != 0)
return ret;
1279 size_t * results_array_size) {
1294#pragma _NEC novector
1300 points[i].coordinates_xyz[0] * point_coordinates_xyz[0] +
1301 points[i].coordinates_xyz[1] * point_coordinates_xyz[1] +
1302 points[i].coordinates_xyz[2] * point_coordinates_xyz[2]);
1310 double min_cos_angle = results_[n - 1].
cos_angle;
1312 for (num_results = n;
1314 !(fabs(min_cos_angle - results_[num_results].
cos_angle) > 0.0);
1321 size_t n, double * point_coordinates_xyz,
1326 size_t num_results_ = *num_results;
1332 double min_cos_angle = results_[num_results_-1].
cos_angle;
1337 double curr_cos_angle =
1338 points[i].coordinates_xyz[0] * point_coordinates_xyz[0] +
1339 points[i].coordinates_xyz[1] * point_coordinates_xyz[1] +
1340 points[i].coordinates_xyz[2] * point_coordinates_xyz[2];
1343 if (curr_cos_angle < min_cos_angle)
continue;
1346 {.
point =
points[i], .cos_angle = curr_cos_angle};
1350 for (j = 0; j < num_results_; ++j) {
1353 &point, results_ + num_results_ - j - 1) < 0) {
1354 results_[num_results_ - j] = results_[num_results_ - j - 1];
1359 results_[num_results_ - j] = point;
1367 if (num_results_ > n) {
1369 size_t new_num_results;
1370 min_cos_angle = results_[n - 1].
cos_angle;
1372 for (new_num_results = n;
1373 (new_num_results < num_results_) &&
1374 !(fabs(min_cos_angle - results_[new_num_results].cos_angle) > 0.0);
1376 num_results_ = new_num_results;
1378 *num_results = num_results_;
1382 results_[num_results_-1].point.coordinates_xyz, point_coordinates_xyz);
1383 }
else return curr_angle;
1387 size_t n,
double * point_coordinates_xyz,
1389 size_t * num_results,
double * dot_stack,
1391 size_t curr_tree_depth) {
1395 (*results)[(*num_results)-1].point.coordinates_xyz,
1396 point_coordinates_xyz);
1401 double dot = dot_stack[curr_tree_depth];
1413 if ((dot < angle.
sin) | (angle.
cos <= 0.0)) {
1419 node->
U_size, results, results_array_size, num_results, angle);
1429 dot_stack[curr_tree_depth] = dot;
1430 node_stack[curr_tree_depth] = node;
1431 flags[curr_tree_depth] = 0;
1444 if ((dot > - angle.
sin) || (angle.
cos <= 0.0)) {
1450 node->
T_size, results, results_array_size, num_results, angle);
1460 dot_stack[curr_tree_depth] = dot;
1461 node_stack[curr_tree_depth] = node;
1462 flags[curr_tree_depth] = 0;
1470 if (curr_tree_depth == 0)
break;
1475 dot = dot_stack[curr_tree_depth];
1476 node = node_stack[curr_tree_depth];
1485 double (*coordinates_xyz)[3],
size_t n,
1486 double ** cos_angles,
1487 size_t * cos_angles_array_size,
1488 double (**result_coordinates_xyz)[3],
1489 size_t * result_coordinates_xyz_array_size,
1490 size_t ** local_point_ids,
1491 size_t * local_point_ids_array_size,
1492 size_t * num_local_point_ids) {
1496 if (cos_angles != NULL)
1501 search,
num_points, coordinates_xyz, (cos_angles!=NULL)?*cos_angles:NULL,
1502 result_coordinates_xyz, result_coordinates_xyz_array_size,
1503 local_point_ids, local_point_ids_array_size, num_local_point_ids);
1505 size_t total_num_local_points = 0;
1507 total_num_local_points += num_local_point_ids[i];
1509 if ((cos_angles != NULL) && (total_num_local_points >
num_points)) {
1512 total_num_local_points);
1514 for (
size_t i =
num_points - 1, offset = total_num_local_points - 1;
1515 i != (size_t)-1; i--) {
1517 for (
size_t j = 0; j < num_local_point_ids[i]; ++j, --offset)
1518 (*cos_angles)[offset] = (*cos_angles)[i];
1525 if (search == NULL) {
1526 memset(num_local_point_ids, 0,
num_points *
sizeof(*num_local_point_ids));
1532 size_t total_num_local_point_ids = 0;
1540 size_t results_array_size = 0;
1546 double * curr_coordinates_xyz = coordinates_xyz[i];
1548 size_t curr_tree_depth = 0;
1550 size_t curr_num_points = 0;
1555 double dot = curr_node->
gc_norm_vector[0]*curr_coordinates_xyz[0] +
1559 dot_stack[curr_tree_depth] = dot;
1560 node_stack[curr_tree_depth] = curr_node;
1561 flags[curr_tree_depth] = 0;
1566 if (curr_node->
U_size < n) {
1569 curr_num_points = curr_node->
U_size + curr_node->
T_size;
1573 flags[curr_tree_depth] =
U_FLAG;
1574 curr_num_points = curr_node->
U_size;
1578 flags[curr_tree_depth] =
U_FLAG;
1579 curr_node = curr_node->
U;
1584 if (curr_node->
T_size < n) {
1587 curr_num_points = curr_node->
U_size + curr_node->
T_size;
1592 flags[curr_tree_depth] =
T_FLAG;
1593 curr_num_points = curr_node->
T_size;
1598 flags[curr_tree_depth] =
T_FLAG;
1599 curr_node = curr_node->
T;
1606 assert(curr_num_points > 0);
1608 size_t num_results =
1610 n,
points, curr_num_points, curr_coordinates_xyz,
1611 &results, &results_array_size);
1615 n, curr_coordinates_xyz, &results, &results_array_size, &num_results,
1616 dot_stack, node_stack, flags, curr_tree_depth);
1620 total_num_local_point_ids + num_results);
1621 size_t * local_point_ids_ =
1622 (*local_point_ids) + total_num_local_point_ids;
1623 double * cos_angles_;
1624 if (cos_angles != NULL) {
1626 total_num_local_point_ids + num_results);
1627 cos_angles_ = (*cos_angles) + total_num_local_point_ids;
1631 double (*result_coordinates_xyz_)[3];
1632 if (result_coordinates_xyz != NULL) {
1634 *result_coordinates_xyz_array_size,
1635 total_num_local_point_ids + num_results);
1636 result_coordinates_xyz_ =
1637 (*result_coordinates_xyz) + total_num_local_point_ids;
1639 result_coordinates_xyz_ = NULL;
1642 for (
size_t j = 0; j < num_results; ++j) {
1644 local_point_ids_[j] = results[j].
point.
idx;
1645 if (cos_angles_ != NULL) cos_angles_[j] = results[j].
cos_angle;
1646 if (result_coordinates_xyz_ != NULL) {
1653 num_local_point_ids[i] = num_results;
1654 total_num_local_point_ids += num_results;
1666 size_t n,
size_t ** local_point_ids,
size_t * local_point_ids_array_size,
1667 size_t * num_local_point_ids) {
1669 if (num_bnd_circles == 0)
return;
1671 if (search == NULL) {
1673 num_local_point_ids, 0, num_bnd_circles *
sizeof(*num_local_point_ids));
1679 size_t total_num_local_point_ids = 0;
1687 size_t results_array_size = 0;
1689 for (
size_t i = 0; i < num_bnd_circles; ++i) {
1693 double * curr_coordinates_xyz = bnd_circles[i].
base_vector;
1695 size_t curr_tree_depth = 0;
1697 size_t curr_num_points = 0;
1702 double dot = curr_node->
gc_norm_vector[0]*curr_coordinates_xyz[0] +
1706 dot_stack[curr_tree_depth] = dot;
1707 node_stack[curr_tree_depth] = curr_node;
1708 flags[curr_tree_depth] = 0;
1713 if (curr_node->
U_size < n) {
1716 curr_num_points = curr_node->
U_size + curr_node->
T_size;
1720 flags[curr_tree_depth] =
U_FLAG;
1721 curr_num_points = curr_node->
U_size;
1725 flags[curr_tree_depth] =
U_FLAG;
1726 curr_node = curr_node->
U;
1731 if (curr_node->
T_size < n) {
1734 curr_num_points = curr_node->
U_size + curr_node->
T_size;
1739 flags[curr_tree_depth] =
T_FLAG;
1740 curr_num_points = curr_node->
T_size;
1745 flags[curr_tree_depth] =
T_FLAG;
1746 curr_node = curr_node->
T;
1754 curr_num_points > 0,
1755 "ERROR(yac_point_sphere_part_search_NNN_bnd_circle): "
1756 "insufficient number of points");
1758 size_t num_results =
1760 n,
points, curr_num_points, curr_coordinates_xyz,
1761 &results, &results_array_size);
1765 n, curr_coordinates_xyz, &results, &results_array_size, &num_results,
1766 dot_stack, node_stack, flags, curr_tree_depth);
1768 for (; num_results > 0; --num_results)
1769 if (results[num_results-1].cos_angle >= bnd_circles[i].inc_angle.cos)
1774 total_num_local_point_ids + num_results);
1775 size_t * local_point_ids_ =
1776 (*local_point_ids) + total_num_local_point_ids;
1778 for (
size_t j = 0; j < num_results; ++j)
1779 local_point_ids_[j] = results[j].point.idx;
1781 num_local_point_ids[i] = num_results;
1782 total_num_local_point_ids += num_results;
1804 "ERRROR(yac_point_sphere_part_search_NNN_ubound): "
1805 "invalid point sphere part search (has to be != NULL)");
1807 n > 0,
"ERROR(yac_point_sphere_part_search_NNN_ubound): "
1808 "invalid n (has to be > 0)")
1811 size_t temp_angles_array_size = 0;
1819 double * curr_coordinates_xyz = coordinates_xyz[i];
1822 size_t curr_num_points = 0;
1827 double dot = curr_node->
gc_norm_vector[0]*curr_coordinates_xyz[0] +
1834 if (curr_node->
U_size < n) {
1836 curr_num_points = curr_node->
U_size + curr_node->
T_size;
1840 curr_num_points = curr_node->
U_size;
1844 curr_node = curr_node->
U;
1849 if (curr_node->
T_size < n) {
1851 curr_num_points = curr_node->
U_size + curr_node->
T_size;
1856 curr_num_points = curr_node->
T_size;
1861 curr_node = curr_node->
T;
1867 curr_num_points >= n,
1868 "ERROR(yac_point_sphere_part_search_NNN_ubound): "
1869 "failed to find a sufficient number of points");
1878 curr_coordinates_xyz,
points[0].coordinates_xyz);
1879 for (
size_t j = 1; j < curr_num_points; ++j) {
1882 curr_coordinates_xyz,
points[j].coordinates_xyz);
1884 best_angle = curr_angle;
1886 angles[i] = best_angle;
1892 temp_angles, temp_angles_array_size, curr_num_points);
1893 for (
size_t j = 0; j < curr_num_points; ++j)
1896 curr_coordinates_xyz,
points[j].coordinates_xyz);
1899 temp_angles, curr_num_points,
sizeof(*temp_angles),
1902 angles[i] = temp_angles[n-1];
1911 size_t ** overlap_cells,
1912 size_t * overlap_cells_array_size,
1913 size_t * num_overlap_cells,
1914 struct overlaps * search_interval_tree_buffer,
1915 double prev_gc_norm_vector[]) {
1927 *num_overlap_cells + node->
T_size);
1928 memcpy(*overlap_cells + *num_overlap_cells, node->
T,
1929 node->
T_size *
sizeof(**overlap_cells));
1930 *num_overlap_cells += node->
T_size;
1934 overlap_cells_array_size, num_overlap_cells,
1945 *num_overlap_cells + node->
U_size);
1946 memcpy(*overlap_cells + *num_overlap_cells, node->
U,
1947 node->
U_size *
sizeof(**overlap_cells));
1948 *num_overlap_cells += node->
U_size;
1952 overlap_cells_array_size, num_overlap_cells,
1964 double GCp[3], bVp[3];
1971 if ((fabs(bVp[0]) > 1e-9) || (fabs(bVp[1]) > 1e-9) ||
1972 (fabs(bVp[2]) > 1e-9)) {
1975 search_interval.
left=base_angle;
1976 search_interval.
right=base_angle;
1978 search_interval.
left = -M_PI;
1979 search_interval.
right = M_PI;
1985 search_interval, search_interval_tree_buffer);
1988 *num_overlap_cells +
1991 for (
size_t i = 0; i < search_interval_tree_buffer->
num_overlaps;
1993 (*overlap_cells)[(*num_overlap_cells)+i] =
1998 *num_overlap_cells += search_interval_tree_buffer->
num_overlaps;
2003 *num_overlap_cells + node->
I_size);
2004 memcpy(*overlap_cells + *num_overlap_cells, node->
I.
list,
2005 node->
I_size *
sizeof(**overlap_cells));
2006 *num_overlap_cells += node->
I_size;
2013 size_t count,
size_t ** cells,
size_t * num_cells_per_coordinate) {
2017 size_t * temp_search_results = NULL;
2018 size_t temp_search_results_array_size = 0;
2019 size_t num_temp_search_results = 0;
2021 struct overlaps search_interval_tree_buffer = {0, 0, NULL};
2024 num_cells_per_coordinate, 0, count *
sizeof(*num_cells_per_coordinate));
2025 size_t * cells_ = NULL;
2026 size_t cells_array_size = 0;
2027 size_t num_cells = 0;
2030 for (
size_t i = 0; i < count; ++i) {
2032 double * curr_coordinates_xyz = &(coordinates_xyz[i][0]);
2034 num_temp_search_results = 0;
2036 double gc_norm_vector[3] = {0.0,0.0,1.0};
2038 search_point(base_node, curr_coordinates_xyz, &temp_search_results,
2039 &temp_search_results_array_size, &num_temp_search_results,
2040 &search_interval_tree_buffer, gc_norm_vector);
2043 cells_, cells_array_size, num_cells + num_temp_search_results);
2045 memcpy(cells_ + num_cells, temp_search_results,
2046 num_temp_search_results *
sizeof(*temp_search_results));
2047 num_cells_per_coordinate[i] = num_temp_search_results;
2048 num_cells += num_temp_search_results;
2051 free(temp_search_results);
2052 free(search_interval_tree_buffer.
overlap_iv);
2054 *cells =
xrealloc(cells_, num_cells *
sizeof(*cells_));
2059 size_t count,
size_t ** cells,
size_t * num_cells_per_bnd_circle) {
2063 size_t * temp_search_results = NULL;
2064 size_t temp_search_results_array_size = 0;
2065 size_t num_temp_search_results = 0;
2067 struct overlaps search_interval_tree_buffer = {0, 0, NULL};
2070 num_cells_per_bnd_circle, 0, count *
sizeof(*num_cells_per_bnd_circle));
2071 size_t * cells_ = NULL;
2072 size_t cells_array_size = 0;
2073 size_t num_cells = 0;
2076 for (
size_t i = 0; i < count; ++i) {
2078 num_temp_search_results = 0;
2080 double gc_norm_vector[3] = {0.0,0.0,1.0};
2083 &temp_search_results_array_size, &num_temp_search_results,
2084 &search_interval_tree_buffer, gc_norm_vector);
2087 cells_, cells_array_size, num_cells + num_temp_search_results);
2089 memcpy(cells_ + num_cells, temp_search_results,
2090 num_temp_search_results *
sizeof(*temp_search_results));
2091 num_cells_per_bnd_circle[i] = num_temp_search_results;
2092 num_cells += num_temp_search_results;
2095 free(temp_search_results);
2096 free(search_interval_tree_buffer.
overlap_iv);
2098 *cells =
xrealloc(cells_, num_cells *
sizeof(*cells_));
2134 if (search == NULL)
return;
2143 if (search == NULL)
return;
#define ENSURE_ARRAY_SIZE(arrayp, curr_array_size, req_size)
static struct sin_cos_angle get_vector_angle_2(double const a[3], double const b[3])
static double clamp_abs_one(double val)
static const struct sin_cos_angle SIN_COS_ZERO
static int points_are_identically(double const *a, double const *b)
static const struct sin_cos_angle SIN_COS_M_PI
static struct sin_cos_angle sum_angles_no_check(struct sin_cos_angle a, struct sin_cos_angle b)
static void crossproduct_kahan(double const a[], double const b[], double cross[])
static const struct sin_cos_angle SIN_COS_M_PI_2
static int compare_angles(struct sin_cos_angle a, struct sin_cos_angle b)
static void normalise_vector(double v[])
static double get_vector_angle(double const a[3], double const b[3])
static struct sin_cos_angle sin_cos_angle_new(double sin, double cos)
void yac_search_interval_tree(struct interval_node tree[], size_t num_nodes, struct interval query, struct overlaps *overlaps)
void yac_generate_interval_tree(struct interval_node intervals[], size_t num_nodes)
#define xrealloc(ptr, size)
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)
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 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[])
static int compare_point_idx_xyz(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 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 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 struct point_id_xyz * get_unique_points(size_t *num_points, yac_const_coordinate_pointer coordinates_xyz, yac_int const *ids, int const *mask)
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 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 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)
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)
struct bnd_sphere_part_search * yac_bnd_sphere_part_search_new(struct bounding_circle *circles, size_t num_circles)
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 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[])
void yac_delete_point_sphere_part_search(struct point_sphere_part_search *search)
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)
void yac_bnd_sphere_part_search_delete(struct bnd_sphere_part_search *search)
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)
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)
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 int compare_angles_(void const *a, void const *b)
static struct sphere_part_node * get_sphere_part_node()
static void free_sphere_part_tree(struct sphere_part_node tree)
static int compare_point_id_xyz_angle(const void *a, const void *b)
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)
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 swap_partition_data(struct temp_partition_data *a, struct temp_partition_data *b)
static void init_sphere_part_node(struct sphere_part_node *node)
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[])
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)
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)
static void free_point_sphere_part_tree(struct point_sphere_part_node *tree)
static int compare_points_int32_coord(const void *a, const void *b)
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)
algorithm for searching cells and points on a grid
struct sphere_part_node base_node
struct sin_cos_angle inc_angle
angle between the middle point and the boundary of the spherical cap
struct point_id_xyz point
int32_t coordinate_xyz[3]
double coordinates_xyz[3]
struct point_sphere_part_node base_node
struct point_id_xyz * points
struct sin_cos_angle I_angle
struct bounding_circle bnd_circle
struct interval_node * head_node
static struct user_input_data_points ** points
#define YAC_ASSERT(exp, msg)
double const (* yac_const_coordinate_pointer)[3]
double(* yac_coordinate_pointer)[3]