YAC 3.12.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_point_in_cell.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#include "tests.h"
6#include "grids/basic_grid.h"
7#include "geometry.h"
8#include "test_common.h"
9
15#define YAC_RAD (0.01745329251994329576923690768489) // M_PI / 180
16
17int main (void) {
18
19 enum yac_edge_type gc_edges[8] =
20 {YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
21 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
22 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
23 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE};
24 enum yac_edge_type latlon_edges[4] =
25 {YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE};
26
27 double point_3d[3];
28 struct yac_grid_cell cell;
29
30 { // simple inside test
31 LLtoXYZ(0, 0, point_3d);
32 cell = generate_cell_deg((double[]){-1,1,1,-1},
33 (double[]){-1,-1,1,1}, gc_edges, 4);
34
35 if (!yac_point_in_cell(point_3d, cell))
36 PUT_ERR("error in point_in_cell (simple inside test)\n")
37 yac_free_grid_cell(&cell);
38 }
39
40 { // simple outside test
41 LLtoXYZ_deg(0, 1.5, point_3d);
42 cell = generate_cell_deg((double[]){-1,1,1,-1},
43 (double[]){-1,-1,1,1}, gc_edges, 4);
44
45 if (yac_point_in_cell(point_3d, cell))
46 PUT_ERR("error in point_in_cell (simple outside test)\n")
47 yac_free_grid_cell(&cell);
48 }
49
50 { // test point on edge
51 LLtoXYZ_deg(0, 1, point_3d);
52 cell = generate_cell_deg((double[]){-1,1,1,-1},
53 (double[]){-1,-1,1,1}, gc_edges, 4);
54 if (!yac_point_in_cell(point_3d, cell))
55 PUT_ERR("error in point_in_cell (test point on edge)\n")
56 yac_free_grid_cell(&cell);
57 }
58
59 { // test near the pole
60 LLtoXYZ_deg(22.5, 85.1, point_3d);
61 cell = generate_cell_deg((double[]){0,45,45,0},
62 (double[]){85,85,88,88}, latlon_edges, 4);
63
64 if (!yac_point_in_cell(point_3d, cell))
65 PUT_ERR("error in point_in_cell (lon-lat cell near the pole)\n")
66 yac_free_grid_cell(&cell);
67 }
68
69 { // test near the pole
70 LLtoXYZ_deg(22.5, 85.1, point_3d);
71 cell = generate_cell_deg((double[]){0,45,45,0},
72 (double[]){85,85,88,88}, gc_edges, 4);
73
74 if (yac_point_in_cell(point_3d, cell))
75 PUT_ERR("error in point_in_cell (great circle cell near the pole)\n")
76 yac_free_grid_cell(&cell);
77 }
78
79 { // simple test at the pole
80 LLtoXYZ_deg(45, 88.3, point_3d);
81 cell = generate_cell_deg((double[]){0,90,180,270},
82 (double[]){88,88,88,88},
83 (enum yac_edge_type[]){
84 YAC_LAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE,
85 YAC_LAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE}, 4);
86
87 if (!yac_point_in_cell(point_3d, cell))
88 PUT_ERR("error in point_in_cell (lon-lat cell on the pole)\n")
89 yac_free_grid_cell(&cell);
90 }
91
92 { // simple test at the pole
93 LLtoXYZ_deg(45, 88.3, point_3d);
94 cell = generate_cell_deg((double[]){0,90,180,270},
95 (double[]){88,88,88,88}, gc_edges, 4);
96
97 if (yac_point_in_cell(point_3d, cell))
98 PUT_ERR("error in point_in_cell (great circle cell on the pole)\n")
99 yac_free_grid_cell(&cell);
100 }
101
102 { // point and cell at pole
103 for (int i = -1; i <= 1; i += 2) {
104 for (int j = -1; j <= 1; j += 2) {
105 LLtoXYZ_deg(45, 90.0 * (double)i, point_3d);
106 double lat_coords[5] = {88.0, 88.0, 88.0, 88.0, 88.0};
107 for (int k = 0; k < 5; ++k) lat_coords[k] *= (double)j;
108 cell = generate_cell_deg((double[]){0,60,120,180,240},
109 lat_coords, gc_edges, 5);
110
111 if (yac_point_in_cell(point_3d, cell) ^ (i == j))
112 PUT_ERR("error in point_in_cell (great circle cell on the pole)\n")
113 yac_free_grid_cell(&cell);
114 }
115 }
116 }
117
118 { // test point on cell corner and edge
119 LLtoXYZ_deg(1, 1, point_3d);
120 cell = generate_cell_deg((double[]){0,1,1,0},
121 (double[]){0,0,1,1}, latlon_edges, 4);
122
123 if (!yac_point_in_cell(point_3d, cell))
124 PUT_ERR("error in point_in_cell (point on corner of lon-lat cell)\n")
125
126 LLtoXYZ_deg(0.5, 0, point_3d);
127
128 if (!yac_point_in_cell(point_3d, cell))
129 PUT_ERR("error in point_in_cell (point on edge of lon-lat cell)\n")
130 yac_free_grid_cell(&cell);
131 }
132
133 { // test point on cell corner
134 LLtoXYZ_deg(1, 1, point_3d);
135 cell = generate_cell_deg((double[]){0,1,0,-1},
136 (double[]){0,1,2,1}, gc_edges, 4);
137
138 if (!yac_point_in_cell(point_3d, cell))
139 PUT_ERR("error in point_in_cell (point on corner of great circle cell)\n")
140
141 LLtoXYZ_deg(0, 2, point_3d);
142
143 if (!yac_point_in_cell(point_3d, cell))
144 PUT_ERR("error in point_in_cell (point on corner of great circle cell)\n")
145
146 LLtoXYZ_deg(0, 1, point_3d);
147
148 if (!yac_point_in_cell(point_3d, cell))
149 PUT_ERR("error in point_in_cell (point in cell)\n")
150
151 LLtoXYZ_deg(0, 3, point_3d);
152
153 if (yac_point_in_cell(point_3d, cell))
154 PUT_ERR("error in point_in_cell (point in cell)\n")
155 yac_free_grid_cell(&cell);
156 }
157
158 { // test point outside of cell
159 LLtoXYZ(1.5, 0.5, point_3d);
160 cell = generate_cell_deg((double[]){0,1,1,0},
161 (double[]){0,0,1,1}, latlon_edges, 4);
162
163 if (yac_point_in_cell(point_3d, cell))
164 PUT_ERR("error in point_in_cell\n")
165 yac_free_grid_cell(&cell);
166 }
167
168 { // test point on longitude edge of cell
169
170
171 for (unsigned i = 0; i < 2; ++i) {
172 for (int j = 0; j < 3; ++j) {
173
174 enum yac_edge_type edge_type[2][3] =
175 {{YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE},
176 {YAC_LON_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}};
177
178 double point_lat[3] = {40, 45, 50};
179
180 LLtoXYZ_deg(0.0, point_lat[j], point_3d);
181 cell = generate_cell_deg((double[]){0,0,5},
182 (double[]){40,50,45}, edge_type[i], 3);
183
184 if (!yac_point_in_cell(point_3d, cell))
185 PUT_ERR("error in point_in_cell\n")
186 yac_free_grid_cell(&cell);
187 }
188 }
189 }
190
191 { // test point on latitude edge of cell
192
193 for (unsigned i = 0; i < 2; ++i) {
194
195 LLtoXYZ_deg(5.0, 10, point_3d);
196 double y_coords[2][3] = {{10,10,15}, {10,10,5}};
197 cell = generate_cell_deg((double[]){0,10,0}, y_coords[i],
198 (enum yac_edge_type[])
199 {YAC_LAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
200 YAC_GREAT_CIRCLE_EDGE}, 3);
201
202 if (!yac_point_in_cell(point_3d, cell))
203 PUT_ERR("error in point_in_cell\n")
204 yac_free_grid_cell(&cell);
205 }
206 }
207
208 { // test point equator edge of cell
209
210 for (unsigned i = 0; i < 2; ++i) {
211
212 LLtoXYZ_deg(5.0, 0, point_3d);
213 double y_coords[2][3] = {{0,0,5}, {0,0,-5}};
214 cell = generate_cell_deg((double[]){0,10,0},
215 y_coords[i], gc_edges, 3);
216
217 if (!yac_point_in_cell(point_3d, cell))
218 PUT_ERR("error in point_in_cell\n")
219 yac_free_grid_cell(&cell);
220 }
221 }
222
223 { // test point outside of cell
224
225
226 for (unsigned i = 0; i < 2; ++i) {
227
228 LLtoXYZ_deg(0.0, 11, point_3d);
229 enum yac_edge_type edge_type[2][3] =
230 {{YAC_LON_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE},
231 {YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}};
232 cell = generate_cell_deg((double[]){0,0,5},
233 (double[]){5,10,20},
234 edge_type[i], 3);
235
236 if (yac_point_in_cell(point_3d, cell))
237 PUT_ERR("error in point_in_cell\n")
238 yac_free_grid_cell(&cell);
239 }
240 }
241
242 { // test point inside of concave cell
243
244 for (unsigned i = 0; i < 2; ++i) {
245
246 LLtoXYZ_deg(0.0, 7.0, point_3d);
247 enum yac_edge_type edge_type[2][4] =
248 {{YAC_LON_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
249 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE},
250 {YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
251 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}};
252 cell = generate_cell_deg((double[]){0, 0, -5, 5},
253 (double[]){0, 5, 10, 10},
254 edge_type[i], 4);
255
256 if (!yac_point_in_cell(point_3d, cell))
257 PUT_ERR("error in point_in_cell\n")
258 yac_free_grid_cell(&cell);
259 }
260 }
261
262 { // test point outside of concave cell
263
264 for (unsigned i = 0; i < 2; ++i) {
265
266 LLtoXYZ_deg(-2.5, 6.0, point_3d);
267 enum yac_edge_type edge_type[2][4] =
268 {{YAC_LON_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
269 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE},
270 {YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
271 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE}};
272 cell = generate_cell_deg((double[]){0, 0, -5, 5},
273 (double[]){0, 5, 10, 10},
274 edge_type[i], 4);
275
276 if (yac_point_in_cell(point_3d, cell))
277 PUT_ERR("error in point_in_cell\n")
278 yac_free_grid_cell(&cell);
279 }
280 }
281
282 { // from Renes clipping tests
283
284 double x_coords[] = { 0.0, 0.5, 0.0,-0.5};
285 double y_coords[] = {-0.5, 0.0, 0.5, 0.0};
286 cell = generate_cell_deg(x_coords, y_coords, gc_edges, 4);
287
288 for (unsigned i = 0; i < 4; ++i) {
289 LLtoXYZ_deg(x_coords[i], y_coords[i], point_3d);
290 if (!yac_point_in_cell(point_3d, cell))
291 PUT_ERR("error in point_in_cell\n")
292 }
293 yac_free_grid_cell(&cell);
294 }
295
296 { // example from toy
297
298 cell = generate_cell_rad((double[]){-2.4523357834544868,
299 -2.5132741228718340,
300 -2.5742124622891822},
301 (double[]){-0.8955349680874743,
302 -0.8292526951666540,
303 -0.8955349680874748},
304 gc_edges, 3);
305
306 LLtoXYZ(-2.5132741229357842, -0.8511807007280388, point_3d);
307 if (!yac_point_in_cell(point_3d, cell))
308 PUT_ERR("error in point_in_cell\n");
309 yac_free_grid_cell(&cell);
310 }
311
312 { // example from scrip toy
313
314 cell = generate_cell_3d((double[][3]){{-0.98109801397208807,
315 -1.4373158006630782e-09,
316 -0.19351146472502487},
317 {-0.98296571213463857,
318 -0.016290939681434788,
319 -0.18306560040580719},
320 {-0.98550150905720235,
321 -0.016447047624685553,
322 -0.16886761166786302},
323 {-0.98739170499602125,
324 -1.9584538881393463e-09,
325 -0.15829599143708675},
326 {-0.98550150905854683,
327 0.016447044383749162,
328 -0.16886761197567171},
329 {-0.98296571214686901,
330 0.016290936203869018,
331 -0.18306560064960375}},
332 gc_edges, 6);
333
334 if (!yac_point_in_cell((double[]){-0.98437791422272558,
335 -1.2055152618042087e-16,
336 -0.17606851504603632}, cell))
337 PUT_ERR("error in point_in_cell\n");
338 yac_free_grid_cell(&cell);
339 }
340
341 { // test degenerated triangles (with one very short edge)
342
343 double lon[3] = {0,0 + yac_angle_tol,0};
344 double lat[3] = {0,0,1};
345 struct {
346 double lon, lat, xyz[3];
347 int is_inside;
348 } test_points[] =
349 {{.lon = lon[0], .lat = lat[0], .is_inside = 1},
350 {.lon = lon[1], .lat = lat[1], .is_inside = 1},
351 {.lon = lon[2], .lat = lat[2], .is_inside = 1},
352 {.lon = yac_angle_tol * 0.5, .lat = 0, .is_inside = 1},
353 {.lon = yac_angle_tol * 0.5, .lat = 0.5, .is_inside = 1},
354 {.lon = -0.01, .lat = 0, .is_inside = 0},
355 {.lon = 0.01, .lat = 0, .is_inside = 0},
356 {.lon = 0, .lat = 1.001, .is_inside = 0},
357 {.lon = -0.01, .lat = 0.5, .is_inside = 0}};
358 enum{num_test_points = sizeof(test_points) / sizeof(test_points[0])};
359 for (size_t i = 0; i < num_test_points; ++i)
360 LLtoXYZ_deg(test_points[i].lon, test_points[i].lat, test_points[i].xyz);
361
362 // bounding circle that contains all points and the cell
363 struct bounding_circle bnd_circle;
364 LLtoXYZ_deg(0.0, 0.5, bnd_circle.base_vector);
365 bnd_circle.inc_angle =
366 sin_cos_angle_new(sin(1.0 * YAC_RAD), cos(1.0 * YAC_RAD));
367 bnd_circle.sq_crd = DBL_MAX;
368
369 for (int direction = -1; direction <= 1; direction += 2) {
370
371 for (int start = 0; start < 3; ++start) {
372
373 double temp_lon[3], temp_lat[3];
374 for (int edge_idx = start, i = 0; i < 3; ++i, edge_idx += direction) {
375 temp_lon[i] = lon[(edge_idx + 3) % 3];
376 temp_lat[i] = lat[(edge_idx + 3) % 3];
377 }
378 cell = generate_cell_deg(temp_lon, temp_lat, gc_edges, 3);
379
380 for (size_t i = 0; i < num_test_points; ++i) {
381 if (yac_point_in_cell(test_points[i].xyz, cell) !=
382 test_points[i].is_inside)
383 PUT_ERR("error in point_in_cell\n")
384 if (yac_point_in_cell2(test_points[i].xyz, cell, bnd_circle) !=
385 test_points[i].is_inside)
386 PUT_ERR("error in point_in_cell2\n")
387 }
388
389 yac_free_grid_cell(&cell);
390 }
391 }
392 }
393
394 { // test degenerated triangles (with three very short edges)
395
396 double lon[3] = {0,0 + yac_angle_tol,0};
397 double lat[3] = {0,0,0 + yac_angle_tol};
398 struct {
399 double lon, lat, xyz[3];
400 int is_inside;
401 } test_points[] =
402 {{.lon = lon[0], .lat = lat[0], .is_inside = 1},
403 {.lon = lon[1], .lat = lat[1], .is_inside = 1},
404 {.lon = lon[2], .lat = lat[2], .is_inside = 1},
405 {.lon = yac_angle_tol * 0.5, .lat = 0, .is_inside = 1},
406 {.lon = 0, .lat = yac_angle_tol * 0.5,
407 .is_inside = 1},
408 {.lon = -0.001, .lat = 0, .is_inside = 0},
409 {.lon = 0.001, .lat = 0, .is_inside = 0},
410 {.lon = 0, .lat = 0.001, .is_inside = 0},
411 {.lon = 0, .lat = -0.001, .is_inside = 0}};
412 size_t num_test_points = sizeof(test_points) / sizeof(test_points[0]);
413 for (size_t i = 0; i < num_test_points; ++i)
414 LLtoXYZ_deg(test_points[i].lon, test_points[i].lat, test_points[i].xyz);
415
416 // bounding circle that contains all points and the cell
417 struct bounding_circle bnd_circle;
418 LLtoXYZ_deg(0.0, 0.0, bnd_circle.base_vector);
419 bnd_circle.inc_angle =
420 sin_cos_angle_new(sin(0.5 * YAC_RAD), cos(0.5 * YAC_RAD));
421 bnd_circle.sq_crd = DBL_MAX;
422
423 for (int direction = -1; direction <= 1; direction += 2) {
424
425 for (int start = 0; start < 3; ++start) {
426
427 double temp_lon[3], temp_lat[3];
428 for (int edge_idx = start, i = 0; i < 3; ++i, edge_idx += direction) {
429 temp_lon[i] = lon[(edge_idx + 3) % 3];
430 temp_lat[i] = lat[(edge_idx + 3) % 3];
431 }
432 cell = generate_cell_deg(temp_lon, temp_lat, gc_edges, 3);
433
434 for (size_t i = 0; i < num_test_points; ++i) {
435 if (yac_point_in_cell(test_points[i].xyz, cell) !=
436 test_points[i].is_inside)
437 PUT_ERR("error in point_in_cell\n")
438 if (yac_point_in_cell2(test_points[i].xyz, cell, bnd_circle) !=
439 test_points[i].is_inside)
440 PUT_ERR("error in point_in_cell2\n")
441 }
442
443 yac_free_grid_cell(&cell);
444 }
445 }
446 }
447
448 { // example from HD
449
450 cell = generate_cell_3d((double[][3]){{0.39177807766083478,
451 0.48669359266172357,
452 0.78079400915120067},
453 {0.39106981008544894,
454 0.48726288480996138,
455 0.78079400915120067},
456 {0.39178019297633498,
457 0.48814800354785914,
458 0.77988448312789571},
459 {0.39248974712806789,
460 0.48757767727376555,
461 0.77988448312789571}},
462 latlon_edges, 4);
463
464 if (!yac_point_in_cell((double[]){0.39178019303094858,
465 0.48814800355996663,
466 0.77988448309288172}, cell))
467 PUT_ERR("error in point_in_cell\n");
468 yac_free_grid_cell(&cell);
469 }
470
471 return TEST_EXIT_CODE;
472}
473
basic_grid.h
yac_point_in_cell
int yac_point_in_cell(double point_coords[3], struct yac_grid_cell cell)
Definition check_overlap.c:465
yac_point_in_cell2
int yac_point_in_cell2(double point_coords[3], struct yac_grid_cell cell, struct bounding_circle bnd_circle)
Definition check_overlap.c:474
LLtoXYZ_deg
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:269
yac_angle_tol
#define yac_angle_tol
Definition geometry.h:26
sin_cos_angle_new
static struct sin_cos_angle sin_cos_angle_new(double sin, double cos)
Definition geometry.h:351
yac_free_grid_cell
void yac_free_grid_cell(struct yac_grid_cell *cell)
Definition grid_cell.c:44
yac_edge_type
yac_edge_type
Definition grid_cell.h:12
YAC_GREAT_CIRCLE_EDGE
@ YAC_GREAT_CIRCLE_EDGE
great circle
Definition grid_cell.h:13
YAC_LAT_CIRCLE_EDGE
@ YAC_LAT_CIRCLE_EDGE
latitude circle
Definition grid_cell.h:14
YAC_LON_CIRCLE_EDGE
@ YAC_LON_CIRCLE_EDGE
longitude circle
Definition grid_cell.h:15
bounding_circle
Definition geometry.h:47
bounding_circle::inc_angle
struct sin_cos_angle inc_angle
angle between the middle point and the boundary of the spherical cap
Definition geometry.h:53
bounding_circle::base_vector
double base_vector[3]
Definition geometry.h:51
bounding_circle::sq_crd
double sq_crd
Definition geometry.h:56
yac_grid_cell
Definition grid_cell.h:18
yac_grid_cell::edge_type
enum yac_edge_type * edge_type
Definition grid_cell.h:20
generate_cell_rad
struct yac_grid_cell generate_cell_rad(double *lon, double *lat, enum yac_edge_type *edge_type, size_t num_corners)
Definition test_common.c:42
generate_cell_3d
struct yac_grid_cell generate_cell_3d(yac_coordinate_pointer coords, enum yac_edge_type *edge_type, size_t num_corners)
Definition test_common.c:48
generate_cell_deg
struct yac_grid_cell generate_cell_deg(double *lon, double *lat, enum yac_edge_type *edge_type, size_t num_corners)
Definition test_common.c:36
test_common.h
gc_edges
static enum yac_edge_type gc_edges[]
Definition test_compute_overlap_area.c:16
latlon_edges
static enum yac_edge_type latlon_edges[4]
Definition test_compute_overlap_area.c:26
YAC_RAD
#define YAC_RAD
Definition test_point_in_cell.c:15
TEST_EXIT_CODE
#define TEST_EXIT_CODE
Definition tests.h:14
PUT_ERR
#define PUT_ERR(string)
Definition tests.h:10
LLtoXYZ
static void LLtoXYZ(double lon, double lat, double p_out[])
Definition toy_scrip.c:587