YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
intersection.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 <stdlib.h>
6#include <math.h>
7#include <stdio.h>
8
9#include "utils_core.h"
10#include "geometry.h"
11
12// angle tolerance
13static double const tol = 1.0e-12;
14
15enum {
16 p_on_a = 1,
17 q_on_a = 2,
18 p_on_b = 4,
19 q_on_b = 8,
21};
22
23enum {
29};
30
39static int gcxgc(
40 double const norm_vector_a[3], double const norm_vector_b[3],
41 double p[3], double q[3]) {
42
43 // compute p and q
44 // p' = (a X b) X (c X d)
45 // p = p' / length(p')
46 // q = -p
47 double temp_p[3];
48 crossproduct_kahan(norm_vector_a, norm_vector_b, temp_p);
49
50 double sq_sin_plane_angle =
51 temp_p[0] * temp_p[0] + temp_p[1] * temp_p[1] + temp_p[2] * temp_p[2];
52
53 // if both circles are on the same plane
54 if (sq_sin_plane_angle <= yac_sq_angle_tol) return -1;
55
56 double scale = 1.0 / sqrt(sq_sin_plane_angle);
57 p[0] = (temp_p[0] * scale);
58 p[1] = (temp_p[1] * scale);
59 p[2] = (temp_p[2] * scale);
60 q[0] = -p[0], q[1] = -p[1], q[2] = -p[2];
61
62 return 2;
63}
64
65#if defined __INTEL_COMPILER
66#pragma intel optimization_level 0
67#elif defined _CRAYC
68#pragma _CRI noopt
69#endif
70
71static int loncxlatc(
72 double const gc_norm_vector[3], double const z,
73 double p[3], double q[3]) {
74
75 // based on gcxlatc with some optimisations for circles of longitude
76
77 // square of z of lat circle
78 double sq_z = z * z;
79
80 double b = sqrt(MAX(1.0 - sq_z,0.0));
81 double b_n0 = b * gc_norm_vector[0];
82 double b_n1 = b * gc_norm_vector[1];
83
84 p[0] = b_n1;
85 p[1] = -b_n0;
86 p[2] = z;
87 q[0] = -p[0];
88 q[1] = -p[1];
89 q[2] = z;
90
91 return 2;
92}
93
94static int gcxlatc(
95 double const gc_norm_vector[3], double const z,
96 double p[3], double q[3]) {
97
98/* How to compute intersection point(s) between great circle (given by its
99 * norm vector) and a circle of latitude (give by z coordinate):
100 *
101 * n = norm vector of great circle
102 * x = vector to intersection point(s)
103 * z = cos(lat) // where lat is the latitude of the circle of latitude
104 *
105 * I. n * x = 0 // intersection is orthogonal to the norm vector
106 * II. x * x = 1 // intersection vector has a length of 1
107 * III. x_2 = z // z-coordinate of intersection point is z
108 *
109 * from I we can derive:
110 * x_1=-(n_0*x_0+n_2*z)/n_1
111 *
112 * x_1 and x_2 in II:
113 * x_0^2+(n_0*x_0+n_2*z)^2/n_1^2+z^2=1
114 * x_0^2+(2*n_0*n_2*z)/(n_0^2+n_1^2)*x_0+
115 * (n_2^2*z^2+n_1^2*z^2-n_1^2)/(n_0^2+n_1^2)=0
116 *
117 * with M = n_0^2+n_1^2:
118 * x_0^2 + (2*n_0*n_2*z)/M*x_0 + (n_2^2*z^2+n_1^2*z^2-n_1^2)/M = 0
119 *
120 * x_0 = -(n_0*n_2*z)/M
121 * +- sqrt((n_0*n_2*z)^2/M^2 - (n_2^2*z^2+n_1^2*z^2-n_1^2)/M)
122 * = -(*n_2*z/M)*n_0 +- sqrt(M-z^2)/M*n_1
123 *
124 * with a = -n_2*z/M and b = sqrt(M-z^2)/M:
125 * x_0 = a * n_0 +- b * n_1
126 * x_1 = a * n_1 -+ b * n_0
127 */
128
129 // if the great circle is a circle of longitude
130 if (fabsl(gc_norm_vector[2]) < yac_angle_low_tol)
131 return loncxlatc(gc_norm_vector, z, p, q);
132
133 double sq_n[3] = {gc_norm_vector[0] * gc_norm_vector[0],
134 gc_norm_vector[1] * gc_norm_vector[1],
135 gc_norm_vector[2] * gc_norm_vector[2]};
136
137 // square of highest z-value of great circle
138 double max_sq_gc_z = sq_n[0] + sq_n[1];
139
140 // square of z of lat circle
141 double sq_z = z * z;
142
143 // if both circles are the equator -> identical circles
144 if ((max_sq_gc_z < yac_sq_angle_tol) && (sq_z < yac_sq_angle_tol)) return -1;
145
146 // determine number of intersections
147 int num_intersections =
148 ((max_sq_gc_z > (sq_z - tol)) - (max_sq_gc_z < (sq_z + tol))) + 1;
149
150 switch (num_intersections) {
151
152 // if no intersection is possible
153 case (0): break;
154
155 // the great circles touches the circle of latitude in a single point
156 // ==> (b = 0)
157 case (1): {
158
159 double a = -(gc_norm_vector[2] * z) / max_sq_gc_z;
160
161 p[0] = ((q[0] = gc_norm_vector[0] * a));
162 p[1] = ((q[1] = gc_norm_vector[1] * a));
163 p[2] = z, q[2] = z;
164 break;
165 }
166
167 default:
168 case (2): {
169
170 double a = -(gc_norm_vector[2] * z) / max_sq_gc_z;
171 double b = sqrt(max_sq_gc_z - sq_z)/max_sq_gc_z;
172 double a_n0 = a * gc_norm_vector[0];
173 double a_n1 = a * gc_norm_vector[1];
174 double b_n0 = b * gc_norm_vector[0];
175 double b_n1 = b * gc_norm_vector[1];
176
177 p[0] = a_n0 + b_n1;
178 p[1] = a_n1 - b_n0;
179 p[2] = z;
180 q[0] = a_n0 - b_n1;
181 q[1] = a_n1 + b_n0;
182 q[2] = z;
183 break;
184 }
185 }
186
187 return num_intersections;
188}
189
190#if defined _CRAYC
191#pragma _CRI opt
192#endif
193
194static int loncxlonc(
195 double const norm_vector_a[3], double const norm_vector_b[3],
196 double p[3], double q[3]) {
197
198 double sq_plane_diff =
199 MIN((norm_vector_a[0] - norm_vector_b[0]) *
200 (norm_vector_a[0] - norm_vector_b[0]) +
201 (norm_vector_a[1] - norm_vector_b[1]) *
202 (norm_vector_a[1] - norm_vector_b[1]),
203 (norm_vector_a[0] + norm_vector_b[0]) *
204 (norm_vector_a[0] + norm_vector_b[0]) +
205 (norm_vector_a[1] + norm_vector_b[1]) *
206 (norm_vector_a[1] + norm_vector_b[1]));
207
208 // if both circles are on the same circle of longitude
209 if (sq_plane_diff < yac_sq_angle_tol) return -1;
210
211 // the intersection points of the two circles of longitude are the poles
212 p[0] = 0, p[1] = 0; p[2] = 1;
213 q[0] = 0, q[1] = 0; q[2] = -1;
214
215 return 2;
216}
217
222static int pxgc(
223 double const vec[3], double const norm_vector[3], double p[3], double q[3]) {
224
225 // compute the angle between the point vector and the norm vector
226 // we use kahan summation to increase precision
227 double dot = vec[0] * norm_vector[0];
228 double err = fma(vec[0], norm_vector[0], -dot);
229 dot = fma(vec[1], norm_vector[1], dot);
230 err += fma(vec[1], norm_vector[1], -vec[1]*norm_vector[1]);
231 dot = fma(vec[2], norm_vector[2], dot);
232 err += fma(vec[2], norm_vector[2], -vec[2]*norm_vector[2]);
233 dot += err;
234
235 // if the point is not on the plane
236 if (fabsl(dot) > yac_angle_tol) return 0;
237
238 q[0] = -((p[0] = vec[0]));
239 q[1] = -((p[1] = vec[1]));
240 q[2] = -((p[2] = vec[2]));
241 return 1;
242}
243
244static int pxlatc(
245 double const vec[3], double z, double p[3], double q[3]) {
246
247 // if the point is not on the plane of the circle of latitude
248 if (fabs(vec[2] - z) > tol) return 0;
249
250 q[0] = -((p[0] = vec[0]));
251 q[1] = -((p[1] = vec[1]));
252 q[2] = ((p[2] = vec[2]));
253 return 1;
254}
255
256static int pxp(
257 double const vec_a[3], double const vec_b[3], double p[3], double q[3]) {
258
259 // if the points are identical
260 if (points_are_identically(vec_a, vec_b)) {
261 q[0] = -((p[0] = vec_a[0]));
262 q[1] = -((p[1] = vec_a[1]));
263 q[2] = -((p[2] = vec_a[2]));
264 return 1;
265 } else {
266 return 0;
267 }
268}
269
271 struct yac_circle a, struct yac_circle b, double p[3], double q[3]) {
272
274 ((a.type == GREAT_CIRCLE) ||
275 (a.type == LON_CIRCLE) ||
276 (a.type == LAT_CIRCLE) ||
277 (a.type == POINT)) &&
278 ((b.type == LON_CIRCLE) ||
279 (b.type == GREAT_CIRCLE) ||
280 (b.type == LAT_CIRCLE) ||
281 (b.type == POINT)),
282 "ERROR(yac_circle_intersect): unknown circle type.")
283
284 if (((a.type == GREAT_CIRCLE) && (b.type == GREAT_CIRCLE)) ||
285 ((a.type == GREAT_CIRCLE) && (b.type == LON_CIRCLE)) ||
286 ((a.type == LON_CIRCLE) && (b.type == GREAT_CIRCLE))) {
287 return gcxgc(a.data.gc.norm_vector, b.data.gc.norm_vector, p, q);
288 } else if (((a.type == GREAT_CIRCLE) && (b.type == LAT_CIRCLE))) {
289 return gcxlatc(a.data.gc.norm_vector, b.data.lat.z, p, q);
290 } else if (((a.type == LAT_CIRCLE) && (b.type == GREAT_CIRCLE))) {
291 return gcxlatc(b.data.gc.norm_vector, a.data.lat.z, p, q);
292 } else if ((a.type == LON_CIRCLE) && (b.type == LON_CIRCLE)) {
293 return loncxlonc(a.data.lon.norm_vector, b.data.lon.norm_vector, p, q);
294 } else if ((a.type == LON_CIRCLE) && (b.type == LAT_CIRCLE)) {
295 return loncxlatc(a.data.lon.norm_vector, b.data.lat.z, p, q);
296 } else if ((a.type == LAT_CIRCLE) && (b.type == LON_CIRCLE)) {
297 return loncxlatc(b.data.lon.norm_vector, a.data.lat.z, p, q);
298 } else if ((a.type == POINT) && (b.type == GREAT_CIRCLE)) {
299 return pxgc(a.data.p.vec, b.data.gc.norm_vector, p, q);
300 } else if ((a.type == POINT) && (b.type == LON_CIRCLE)) {
301 return pxgc(a.data.p.vec, b.data.gc.norm_vector, p, q);
302 } else if ((a.type == POINT) && (b.type == LAT_CIRCLE)) {
303 return pxlatc(a.data.p.vec, b.data.lat.z, p, q);
304 } else if ((a.type == POINT) && (b.type == POINT)) {
305 return pxp(a.data.p.vec, b.data.p.vec, p, q);
306 } else if ((a.type == GREAT_CIRCLE) && (b.type == POINT)) {
307 return pxgc(b.data.p.vec, a.data.gc.norm_vector, p, q);
308 } else if ((a.type == LON_CIRCLE) && (b.type == POINT)) {
309 return pxgc(b.data.p.vec, a.data.gc.norm_vector, p, q);
310 } else /*if ((a.type == LAT_CIRCLE) && (b.type == POINT))*/ {
311 return pxlatc(b.data.p.vec, a.data.lat.z, p, q);
312 }
313}
314
319static inline int vector_is_between (
320 double const a[], double const b[], double const p[],
321 double sq_len_diff_ab) {
322
323 // In case the angle between the vectors a and b is 180 degree, the angle
324 // between the vectors pa and pb is exactly 90 degree (Pythagorean theorem).
325 // In all other cases the angle between pa and pb must be bigger than 90
326 // degree for p to be between a and b
327 // from this we can deduce:
328 // ||ab||^2 >= ||ap||^2 + ||bp||^2 => (Pythagorean theorem)
329
330 // the tol is required in case p is very close to a or b
331 double sq_len_diff_vec_ap = sq_len_diff_vec(a,p);
332 double sq_len_diff_vec_bp = sq_len_diff_vec(b,p);
333 return sq_len_diff_ab + tol >= sq_len_diff_vec_ap + sq_len_diff_vec_bp;
334}
335
337 double const a[], double const b[], double const p[]) {
338
339/* determines whether p is between a and b
340 (a, b, p have the same latitude)*/
341
342/* I. a_0 * alpha + b_0 * beta = p_0
343 II. a_1 * alpha + b_1 * beta = p_1
344
345 if alpha > 0 and beta > 0 -> p is between a and b */
346
347 if (fabs(fabs(a[2]) - 1.0) < tol) return 1;
348
349 double fabs_a[2] = {fabs(a[0]), fabs(a[1])};
350 int flag = fabs_a[0] > fabs_a[1];
351 double a_0 = a[flag], a_1 = a[flag^1];
352 double b_0 = b[flag], b_1 = b[flag^1];
353 double p_0 = p[flag], p_1 = p[flag^1];
354
355 double temp = b_0 - (b_1 * a_0) / a_1;
356
358 (fabs(temp) > tol),
359 "ERROR(vector_is_between_lat): "
360 "routine does not support zero length edges")
361 /*// if a and b are nearly identical
362 if (fabs(temp) < tol)
363 return (fabs(a_0 - p_0) < tol) && (fabs(a_1 - p_1) < tol);*/
364
365 double beta = (p_0 - (p_1 * a_0) / a_1) / temp;
366
367 if (beta < -tol) return 0;
368
369 double alpha = (p_1 - b_1 * beta) / a_1;
370
371 return alpha > -tol;
372}
373
375 enum yac_edge_type edge_type, double const a[3], double const b[3],
376 double middle[3]) {
377
379 (edge_type == YAC_GREAT_CIRCLE_EDGE) ||
380 (edge_type == YAC_LAT_CIRCLE_EDGE) ||
381 (edge_type == YAC_LON_CIRCLE_EDGE),
382 "ERROR(compute_edge_middle_point_vec): invalide edge_type")
383
384 if (edge_type == YAC_LAT_CIRCLE_EDGE) {
385 middle[0] = a[0]+b[0], middle[1] = a[1]+b[1], middle[2] = a[2];
386 double temp = middle[0]*middle[0] + middle[1]*middle[1];
387 double scale = sqrt((1.0 - middle[2] * middle[2])/temp);
388 middle[0] *= scale;
389 middle[1] *= scale;
390 } else {
391 middle[0] = a[0]+b[0], middle[1] = a[1]+b[1], middle[2] = a[2]+b[2];
392 normalise_vector(middle);
393 }
394}
395
397 enum yac_edge_type edge_type, double const a[3], double const b[3],
398 double const c[3], double const d[3], double p[3], double q[3],
399 int intersection_type) {
400
402 (intersection_type == no_intersection) ||
403 (intersection_type == a_between_cd + b_between_cd) ||
404 (intersection_type == a_between_cd + b_between_cd + c_between_ab) ||
405 (intersection_type == a_between_cd + b_between_cd + d_between_ab) ||
406 (intersection_type == a_between_cd + c_between_ab) ||
407 (intersection_type == a_between_cd + c_between_ab + d_between_ab) ||
408 (intersection_type == a_between_cd + d_between_ab) ||
409 (intersection_type == b_between_cd + c_between_ab) ||
410 (intersection_type == b_between_cd + c_between_ab + d_between_ab) ||
411 (intersection_type == b_between_cd + d_between_ab) ||
412 (intersection_type == c_between_ab + d_between_ab) ||
413 (intersection_type == a_between_cd + b_between_cd +
415 "internal error")
416
417 switch (intersection_type) {
418 default:
419
420 // edges do not intersect
421 case (no_intersection):
422 // use middle points of both edges as intersection points
423 compute_edge_middle_point_vec(edge_type, a, b, p);
424 compute_edge_middle_point_vec(edge_type, c, d, q);
426
427 // both vertices of first edge overlap with the second edge or
428 // edges are identical
430 case (a_between_cd + b_between_cd):
433 p[0] = a[0], p[1] = a[1], p[2] = a[2];
434 q[0] = b[0], q[1] = b[1], q[2] = b[2];
435 break;
436
437 // both edges of the second edge overlap with the first edge
438 case (c_between_ab + d_between_ab):
441 p[0] = c[0], p[1] = c[1], p[2] = c[2];
442 q[0] = d[0], q[1] = d[1], q[2] = d[2];
443 break;
444
445 // only one vertex of each edge overlaps with the other edge
446 case (a_between_cd + c_between_ab):
447 p[0] = a[0], p[1] = a[1], p[2] = a[2];
448 q[0] = c[0], q[1] = c[1], q[2] = c[2];
449 break;
450 case (a_between_cd + d_between_ab):
451 p[0] = a[0], p[1] = a[1], p[2] = a[2];
452 q[0] = d[0], q[1] = d[1], q[2] = d[2];
453 break;
454 case (b_between_cd + c_between_ab):
455 p[0] = b[0], p[1] = b[1], p[2] = b[2];
456 q[0] = c[0], q[1] = c[1], q[2] = c[2];
457 break;
458 case (b_between_cd + d_between_ab):
459 p[0] = b[0], p[1] = b[1], p[2] = b[2];
460 q[0] = d[0], q[1] = d[1], q[2] = d[2];
461 break;
462 }
464}
465
466// computes intersection of great circle edges that are on
467// the same great circle
469 double const a[3], double const b[3], double const c[3], double const d[3],
470 double p[3], double q[3]) {
471
472 double sq_len_diff_ab = sq_len_diff_vec(a, b);
473 double sq_len_diff_cd = sq_len_diff_vec(c, d);
474
475 return
477 YAC_GREAT_CIRCLE_EDGE, a, b, c, d, p, q,
478 (vector_is_between(c, d, a, sq_len_diff_cd) << 0) +
479 (vector_is_between(c, d, b, sq_len_diff_cd) << 1) +
480 (vector_is_between(a, b, c, sq_len_diff_ab) << 2) +
481 (vector_is_between(a, b, d, sq_len_diff_ab) << 3));
482}
483
484// determines the return value for the provided intersection points p and q
486 double const a[3], double const b[3], double const c[3], double const d[3],
487 double const p[3], double const q[3]) {
488
489 // square of the Euclidean distance between the vertices of the two edges
490 double sq_len_diff_ab = sq_len_diff_vec(a, b);
491 double sq_len_diff_cd = sq_len_diff_vec(c, d);
492
493 return (vector_is_between(a, b, p, sq_len_diff_ab) << 0) +
494 (vector_is_between(a, b, q, sq_len_diff_ab) << 1) +
495 (vector_is_between(c, d, p, sq_len_diff_cd) << 2) +
496 (vector_is_between(c, d, q, sq_len_diff_cd) << 3);
497}
498
500 double const a[3], double const b[3], double norm_vector[3]) {
501
502 norm_vector[0] = a[1] * b[2] - a[2] * b[1];
503 norm_vector[1] = a[2] * b[0] - a[0] * b[2];
504 norm_vector[2] = a[0] * b[1] - a[1] * b[0];
505
506 double scale = 1.0 / sqrt(norm_vector[0] * norm_vector[0] +
507 norm_vector[1] * norm_vector[1] +
508 norm_vector[2] * norm_vector[2]);
509 norm_vector[0] *= scale;
510 norm_vector[1] *= scale;
511 norm_vector[2] *= scale;
512}
513
535static int yac_gcxgc_vec(
536 double const a[3], double const b[3], double const c[3], double const d[3],
537 double p[3], double q[3]) {
538
539 double norm_vector_ab[3], norm_vector_cd[3];
540
541 compute_norm_vector(a, b, norm_vector_ab);
542 compute_norm_vector(c, d, norm_vector_cd);
543
544 switch (gcxgc(norm_vector_ab, norm_vector_cd, p, q)) {
545 default:
546 case(2):
547 return yac_check_pq_gcxgc(a, b, c, d, p, q);
548 case(-1):
549 return yac_identical_gcxgc_vec(a, b, c, d, p, q);
550 };
551}
552
571 double const a[3], double const b[3], double const c[3], double const d[3],
572 double p[3], double q[3]) {
573
574 // two circles of latitude can only intersect if they are on the same latitude
575 if (fabs(a[2] - c[2]) > tol) return -1;
576 else
577 return
579 YAC_LAT_CIRCLE_EDGE, a, b, c, d, p, q,
580 (vector_is_between_lat(c, d, a) << 0) +
581 (vector_is_between_lat(c, d, b) << 1) +
582 (vector_is_between_lat(a, b, c) << 2) +
583 (vector_is_between_lat(a, b, d) << 3));
584}
585
587 double const a[3], double const b[3], double norm_vector[3]) {
588
589 norm_vector[0] = a[1] * b[2] - a[2] * b[1];
590 norm_vector[1] = a[2] * b[0] - a[0] * b[2];
591 norm_vector[2] = 0.0;
592
593 double scale = 1.0 / sqrt(norm_vector[0] * norm_vector[0] +
594 norm_vector[1] * norm_vector[1]);
595 norm_vector[0] *= scale;
596 norm_vector[1] *= scale;
597}
598
617 double const a[3], double const b[3], double const c[3], double const d[3],
618 double p[3], double q[3]) {
619
620 double norm_vector_ab[3], norm_vector_cd[3];
621 compute_lon_norm_vector(a, b, norm_vector_ab);
622 compute_lon_norm_vector(c, d, norm_vector_cd);
623
624 switch(loncxlonc(norm_vector_ab, norm_vector_cd, p, q)) {
625
626 default:
627 case(2):
628 return yac_check_pq_gcxgc(a, b, c, d, p, q);
629 case(-1):
630 return yac_identical_gcxgc_vec(a, b, c, d, p, q);
631 };
632}
633
634// determines the return value for the provided intersection points p and q
636 double const a[3], double const b[3], double const c[3], double const d[3],
637 double const p[3], double const q[3]) {
638
639 // square of the Euclidean distance between the vertices of the two edges
640 double sq_len_diff_ab = sq_len_diff_vec(a, b);
641
642 return
643 (vector_is_between(a, b, p, sq_len_diff_ab) << 0) +
644 (vector_is_between(a, b, q, sq_len_diff_ab) << 1) +
645 (vector_is_between_lat(c, d, p) << 2) +
646 (vector_is_between_lat(c, d, q) << 3);
647}
648
667 double const a[3], double const b[3], double const c[3], double const d[3],
668 double p[3], double q[3]) {
669
670 double norm_vector_ab[3];
671 compute_lon_norm_vector(a, b, norm_vector_ab);
672
673 int num_intersections = loncxlatc(norm_vector_ab, c[2], p, q);
674 YAC_ASSERT(num_intersections == 2, "ERROR(yac_loncxlatc_vec): internal error")
675
676 return yac_check_pq_gcxlatc(a, b, c, d, p, q);
677}
678
702#if defined __INTEL_COMPILER
703#pragma intel optimization_level 0
704#elif defined _CRAYC
705#pragma _CRI noopt
706#endif
707
709 double const a[3], double const b[3], double const c[3], double const d[3],
710 double p[3], double q[3]) {
711
712 double norm_vector_ab[3];
713 compute_norm_vector(a, b, norm_vector_ab);
714
715 int num_intersections = gcxlatc(norm_vector_ab, c[2], p, q);
716
717 switch (num_intersections) {
718 default:
719 case(2):
720 case(1):
721 return yac_check_pq_gcxlatc(a, b, c, d, p, q);
722 case(0):
723 return -1;
724 case(-1):
725 return yac_latcxlatc_vec(a, b, c, d, p, q);
726 }
727}
728
729#if defined _CRAYC
730#pragma _CRI opt
731#endif
732
733static int yac_pxp_vec(
734 double const a[3], double const b[3], double p[3], double q[3]) {
735
736 return (pxp(a, b, p, q))?(p_on_a + p_on_b):-1;
737}
738
739static int yac_pxgc_vec(
740 double const point[3], double const a[3], double const b[3],
741 double p[3], double q[3]) {
742
743 // compute norm vector for the plane of ab
744 double norm_ab[3];
745 compute_norm_vector(a, b, norm_ab);
746
747 // if the point is on the plane of ab
748 if (pxgc(point, norm_ab, p, q)) {
749
750 int ret_value = p_on_a;
751 double sq_len_diff_ab = sq_len_diff_vec(a, b);
752
753 if (vector_is_between(a, b, p, sq_len_diff_ab)) ret_value |= p_on_b;
754 else if (vector_is_between(a, b, q, sq_len_diff_ab)) ret_value |= q_on_b;
755
756 return ret_value;
757
758 } else {
759 return -1;
760 }
761}
762
763static int yac_pxlat_vec(
764 double const point[3], double const a[3], double const b[3],
765 double p[3], double q[3]) {
766
767 if (pxlatc(point, a[2], p, q)) {
768
769 int ret_value = p_on_a;
770
771 if (vector_is_between_lat(a, b, p)) ret_value |= p_on_b;
772 else if (vector_is_between_lat(a, b, q)) ret_value |= q_on_b;
773
774 return ret_value;
775 } else {
776 return -1;
777 }
778}
779
780static inline int adjust_ret_value(int ret_value) {
781 return (ret_value & (~(1 + 2 + 4 + 8))) +
782 ((ret_value & (1 + 2)) << 2) +
783 ((ret_value & (4 + 8)) >> 2);
784}
785
787 enum yac_edge_type edge_type_a, double const a[3], double const b[3],
788 enum yac_edge_type edge_type_b, double const c[3], double const d[3],
789 double p[3], double q[3]) {
790
792 ((edge_type_a == YAC_LON_CIRCLE_EDGE) ||
793 (edge_type_a == YAC_LAT_CIRCLE_EDGE) ||
794 (edge_type_a == YAC_GREAT_CIRCLE_EDGE)) &&
795 ((edge_type_b == YAC_LON_CIRCLE_EDGE) ||
796 (edge_type_b == YAC_LAT_CIRCLE_EDGE) ||
797 (edge_type_b == YAC_GREAT_CIRCLE_EDGE)), "ERROR: unknown edge type.")
798
799 int edge_a_is_point = points_are_identically(a, b);
800 int edge_b_is_point = points_are_identically(c, d);
801
802 // if both edges are points
803 if (edge_a_is_point && edge_b_is_point) {
804
805 return yac_pxp_vec(a, c, p, q);
806
807 // if one edge is a point
808 } else if (edge_a_is_point || edge_b_is_point) {
809
810 enum yac_edge_type edge_type[2] = {edge_type_a, edge_type_b};
811 double const * edge[2][2] = {{a,b}, {c,d}};
812 double const * point[2] = {c,a};
813
814 int ret_value;
815 switch (edge_type[edge_a_is_point]) {
816 case (YAC_LAT_CIRCLE_EDGE):
817 ret_value =
819 point[edge_a_is_point],
820 edge[edge_a_is_point][0], edge[edge_a_is_point][1], p, q);
821 break;
822 default:
823 case (YAC_LON_CIRCLE_EDGE):
825 ret_value =
827 point[edge_a_is_point],
828 edge[edge_a_is_point][0], edge[edge_a_is_point][1], p, q);
829 break;
830 }
831
832 if (edge_a_is_point) return ret_value;
833 else return adjust_ret_value(ret_value);
834
835 // if both edges are on circles of latitude
836 } else if (edge_type_a == YAC_LAT_CIRCLE_EDGE &&
837 edge_type_b == YAC_LAT_CIRCLE_EDGE) {
838
839 return yac_latcxlatc_vec(a, b, c, d, p, q);
840
841 // if both edges are on circle of longitude
842 } else if (edge_type_a == YAC_LON_CIRCLE_EDGE &&
843 edge_type_b == YAC_LON_CIRCLE_EDGE) {
844
845 return yac_loncxlonc_vec(a, b, c, d, p, q);
846
847 // if both edges are on great circles
848 } else if ((edge_type_a == YAC_GREAT_CIRCLE_EDGE &&
849 edge_type_b == YAC_GREAT_CIRCLE_EDGE) ||
850 (edge_type_a == YAC_LON_CIRCLE_EDGE &&
851 edge_type_b == YAC_GREAT_CIRCLE_EDGE) ||
852 (edge_type_a == YAC_GREAT_CIRCLE_EDGE &&
853 edge_type_b == YAC_LON_CIRCLE_EDGE)) {
854
855 return yac_gcxgc_vec(a, b, c, d, p, q);
856
857 // if one edge a is on a great circle and edge b on a circle of latitude
858 } else if (edge_type_a == YAC_GREAT_CIRCLE_EDGE &&
859 edge_type_b == YAC_LAT_CIRCLE_EDGE) {
860
861 return yac_gcxlatc_vec(a, b, c, d, p, q);
862
863 // if one edge a is on a circle of latitude and edge b on a great circle
864 } else if (edge_type_a == YAC_LAT_CIRCLE_EDGE &&
865 edge_type_b == YAC_GREAT_CIRCLE_EDGE ) {
866
867 return adjust_ret_value(yac_gcxlatc_vec(c, d, a, b, p, q));
868
869 // if one edge a is on a circle of longitude and edge b on a circle of latitude
870 } else if (edge_type_a == YAC_LON_CIRCLE_EDGE &&
871 edge_type_b == YAC_LAT_CIRCLE_EDGE) {
872
873 return yac_loncxlatc_vec(a, b, c, d, p, q);
874
875 // if one edge a is on a circle of latitude and edge b on a circle of longitude
876 } else /* if (edge_type_a == YAC_LAT_CIRCLE_EDGE &&
877 edge_type_b == YAC_LON_CIRCLE_EDGE )*/ {
878
879 return adjust_ret_value(yac_loncxlatc_vec(c, d, a, b, p, q));
880 }
881}
882
884 double p[3], double const a[3], double const b[3],
885 enum yac_circle_type circle_type) {
886
888 (circle_type == LON_CIRCLE) ||
889 (circle_type == LAT_CIRCLE) ||
890 (circle_type == GREAT_CIRCLE) ||
891 (circle_type == POINT),
892 "ERROR(yac_point_on_edge): unknown edge type.")
893
895 (circle_type != LAT_CIRCLE) || (fabs(a[2] - b[2]) < tol),
896 "ERROR(yac_point_on_edge): "
897 "LAT_CIRCLE but edge vertices do not have the same latitude");
898
899 switch (circle_type) {
900 default:
901 case (GREAT_CIRCLE):
902 case (LON_CIRCLE):
903 return vector_is_between(a, b, p, sq_len_diff_vec(a, b));
904 case (LAT_CIRCLE):
905 if (fabs(p[2] - a[2]) < tol) return vector_is_between_lat(a, b, p);
906 return 0;
907 case(POINT):
908 return points_are_identically(p, a);
909 }
910}
static int points_are_identically(double const *a, double const *b)
Definition geometry.h:648
static double sq_len_diff_vec(double const a[3], double const b[3])
computes square of the lenght of the vector ab
Definition geometry.h:249
static void crossproduct_kahan(double const a[], double const b[], double cross[])
Definition geometry.h:332
#define yac_sq_angle_tol
Definition geometry.h:35
#define yac_angle_low_tol
Definition geometry.h:37
#define yac_angle_tol
Definition geometry.h:34
static void normalise_vector(double v[])
Definition geometry.h:689
yac_circle_type
Definition geometry.h:72
@ POINT
Definition geometry.h:76
@ LAT_CIRCLE
Definition geometry.h:74
@ GREAT_CIRCLE
Definition geometry.h:73
@ LON_CIRCLE
Definition geometry.h:75
yac_edge_type
Definition grid_cell.h:12
@ YAC_GREAT_CIRCLE_EDGE
great circle
Definition grid_cell.h:13
@ YAC_LAT_CIRCLE_EDGE
latitude circle
Definition grid_cell.h:14
@ YAC_LON_CIRCLE_EDGE
longitude circle
Definition grid_cell.h:15
static int yac_pxlat_vec(double const point[3], double const a[3], double const b[3], double p[3], double q[3])
static int yac_check_pq_gcxlatc(double const a[3], double const b[3], double const c[3], double const d[3], double const p[3], double const q[3])
@ c_between_ab
@ a_between_cd
@ no_intersection
@ d_between_ab
@ b_between_cd
static int yac_latcxlatc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
compute the intersection point two circles of latitude
static int gcxgc(double const norm_vector_a[3], double const norm_vector_b[3], double p[3], double q[3])
int yac_circle_intersect(struct yac_circle a, struct yac_circle b, double p[3], double q[3])
static int yac_identical_gcxgc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
static int vector_is_between(double const a[], double const b[], double const p[], double sq_len_diff_ab)
static int yac_pxgc_vec(double const point[3], double const a[3], double const b[3], double p[3], double q[3])
static void compute_edge_middle_point_vec(enum yac_edge_type edge_type, double const a[3], double const b[3], double middle[3])
static int gcxlatc(double const gc_norm_vector[3], double const z, double p[3], double q[3])
static int vector_is_between_lat(double const a[], double const b[], double const p[])
static int yac_check_pq_gcxgc(double const a[3], double const b[3], double const c[3], double const d[3], double const p[3], double const q[3])
static int pxgc(double const vec[3], double const norm_vector[3], double p[3], double q[3])
static int adjust_ret_value(int ret_value)
static void compute_norm_vector(double const a[3], double const b[3], double norm_vector[3])
static int yac_gcxgc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
static int loncxlonc(double const norm_vector_a[3], double const norm_vector_b[3], double p[3], double q[3])
static int yac_loncxlatc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
compute the intersection point of a meridian and a parallel
int yac_intersect_vec(enum yac_edge_type edge_type_a, double const a[3], double const b[3], enum yac_edge_type edge_type_b, double const c[3], double const d[3], double p[3], double q[3])
int yac_point_on_edge(double p[3], double const a[3], double const b[3], enum yac_circle_type circle_type)
static int pxlatc(double const vec[3], double z, double p[3], double q[3])
@ p_on_a
@ q_on_b
@ q_on_a
@ circles_are_identical
@ p_on_b
static int yac_gcxlatc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
compute the intersection of a great circle with the parallel
static int yac_pxp_vec(double const a[3], double const b[3], double p[3], double q[3])
static int yac_loncxlonc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
compute the intersection point two circles of longitude
static double const tol
static int yac_identical_circles_vec(enum yac_edge_type edge_type, double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3], int intersection_type)
static void compute_lon_norm_vector(double const a[3], double const b[3], double norm_vector[3])
static int loncxlatc(double const gc_norm_vector[3], double const z, double p[3], double q[3])
static int pxp(double const vec_a[3], double const vec_b[3], double p[3], double q[3])
struct yac_circle::@3::@4 gc
union yac_circle::@3 data
double z
Definition geometry.h:90
enum yac_circle_type type
Definition geometry.h:82
struct yac_circle::@3::@6 p
struct yac_circle::@3::@5 lat
struct yac_circle::@3::@4 lon
double norm_vector[3]
Definition geometry.h:86
double vec[3]
Definition geometry.h:93
#define MIN(a, b)
#define MAX(a, b)
#define YAC_ASSERT(exp, msg)
Definition yac_assert.h:16