YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
clipping.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#ifdef HAVE_CONFIG_H
6// Get the definition of the 'restrict' keyword.
7#include "config.h"
8#endif
9
10#include <stdio.h>
11#include <stdlib.h>
12#include <string.h>
13#include <math.h>
14
15#include "geometry.h"
16#include "clipping.h"
17#include "area.h"
18#include "ensure_array_size.h"
19#include "utils_core.h"
20
21//#define YAC_VERBOSE_CLIPPING
22
23static double const tol = 1.0e-12;
24
32
38
39/* internal helper routines for working with linked lists of points */
40
41static void init_point_list(struct point_list * list);
42
43static void reset_point_list(struct point_list * list);
44
45static size_t generate_point_list(
46 struct point_list * list, struct yac_grid_cell cell, int cell_ordering,
47 struct yac_circle * circle_buffer);
48
49static struct point_list_element *
51
52static size_t remove_points(struct point_list * list);
53static size_t remove_zero_length_edges(struct point_list * list);
54
55static void free_point_list(struct point_list * list);
56
57static int get_cell_points_ordering(struct yac_grid_cell cell);
58
59static void generate_cell(struct point_list * list,
60 struct yac_grid_cell * cell);
61
62static enum yac_cell_type get_cell_type(struct yac_grid_cell target_cell);
63
64/* ------------------------- */
65
66static struct yac_grid_cell * overlap_cell_buffer = NULL;
67static size_t overlap_cell_buffer_size = 0;
68
69static inline struct yac_grid_cell * get_overlap_cell_buffer(size_t N) {
70
71 // ensure that there are enough buffer cells
72
74
75 size_t old_overlap_cell_buffer_size = overlap_cell_buffer_size;
76
78
79 for (; old_overlap_cell_buffer_size < overlap_cell_buffer_size;
80 ++old_overlap_cell_buffer_size)
81 yac_init_grid_cell(overlap_cell_buffer + old_overlap_cell_buffer_size);
82 }
83
85}
86
87/* ------------------------- */
88
89static double get_edge_direction(
90 double * ref_corner, double * corner_a, double * corner_b) {
91
92 double edge_norm[3];
93 crossproduct_kahan(corner_a, corner_b, edge_norm);
94 normalise_vector(edge_norm);
95
96 // sine of the angle between the edge and the reference corner
97 double angle =
98 edge_norm[0] * ref_corner[0] +
99 edge_norm[1] * ref_corner[1] +
100 edge_norm[2] * ref_corner[2];
101
102 // if the reference corner is directly on the edge
103 // (for small angles sin(x)==x)
104 if (fabs(angle) < yac_angle_tol) return 0.0;
105
106 return copysign(1.0, angle);
107}
108
110 struct yac_grid_cell * source_cell,
111 struct yac_grid_cell target_cell,
112 double * overlap_areas,
113 double (*overlap_barycenters)[3]) {
114
116 target_cell.num_corners > 2,
117 "ERROR(yac_compute_overlap_info): "
118 "target cell has too few corners")
119
120 struct yac_grid_cell * overlap_buffer = get_overlap_cell_buffer(N);
121 enum yac_cell_type target_cell_type;
122
123 // initialise barycenter coordinates if available
124 if (overlap_barycenters != NULL)
125 for (size_t i = 0; i < N; ++i)
126 for (int j = 0; j < 3; ++j)
127 overlap_barycenters[i][j] = 0.0;
128
129 // if the target is not a triangle
130 if ( target_cell.num_corners > 3 )
131 target_cell_type = get_cell_type (target_cell);
132
133 // special case: target is a triangle or a lon-lat cell -> no triangulation
134 // of the target is required
135 if ( target_cell.num_corners < 4 || target_cell_type == YAC_LON_LAT_CELL ) {
136 yac_cell_clipping ( N, source_cell, target_cell, overlap_buffer);
137 for (size_t i = 0; i < N; ++i) {
138 if (overlap_buffer[i].num_corners > 1) {
139 if (overlap_barycenters == NULL)
140 overlap_areas[i] = yac_huiliers_area (overlap_buffer[i]);
141 else {
142 overlap_areas[i] =
144 overlap_buffer[i], overlap_barycenters[i], 1.0);
146 (overlap_barycenters[i][0] != 0.0) ||
147 (overlap_barycenters[i][1] != 0.0) ||
148 (overlap_barycenters[i][2] != 0.0),
149 "ERROR(yac_compute_overlap_info): "
150 "overlap was computed, still barycenter is sphere origin");
151 normalise_vector(overlap_barycenters[i]);
152 if (overlap_areas[i] < 0.0) {
153 overlap_areas[i] = -overlap_areas[i];
154 overlap_barycenters[i][0] = -overlap_barycenters[i][0];
155 overlap_barycenters[i][1] = -overlap_barycenters[i][1];
156 overlap_barycenters[i][2] = -overlap_barycenters[i][2];
157 }
158 }
159 } else {
160 overlap_areas[i] = 0.0;
161 }
162 }
163 return;
164 }
165
166 // the triangulation algorithm only works for cells that only have great
167 // circle edges
168 // in order to also support other edge types, additional work would be
169 // required
171 target_cell_type == YAC_GREAT_CIRCLE_CELL,
172 "ERROR(yac_compute_overlap_info): invalid target cell type")
173
174 // data structure to hold the triangles of the target cell
175 struct yac_grid_cell target_partial_cell =
176 {.coordinates_xyz = (double[3][3]){{-1}},
177 .edge_type = (enum yac_edge_type[3]){YAC_GREAT_CIRCLE_EDGE,
180 .num_corners = 3};
181 // common node point to all partial target cells
182 double * base_corner = target_cell.coordinates_xyz[0];
183 target_partial_cell.coordinates_xyz[0][0] = base_corner[0];
184 target_partial_cell.coordinates_xyz[0][1] = base_corner[1];
185 target_partial_cell.coordinates_xyz[0][2] = base_corner[2];
186
187 // initialise overlap areas
188 for ( size_t n = 0; n < N; n++) overlap_areas[n] = 0.0;
189
190 // for all triangles of the target cell
191 // (triangles a formed by first corner of the target cells and each edge of
192 // the cell; the first and last edge of the cell already, contain the
193 // first corner, therefore we can skip them)
194 for ( size_t corner_idx = 1;
195 corner_idx < target_cell.num_corners - 1; ++corner_idx ) {
196
197 double * corner_a = target_cell.coordinates_xyz[corner_idx];
198 double * corner_b = target_cell.coordinates_xyz[corner_idx+1];
199
200 // if the current edge has a length of zero
201 if (points_are_identically(corner_a, corner_b)) continue;
202
203 double edge_direction =
204 get_edge_direction(base_corner, corner_a, corner_b);
205
206 target_partial_cell.coordinates_xyz[1][0] = corner_a[0];
207 target_partial_cell.coordinates_xyz[1][1] = corner_a[1];
208 target_partial_cell.coordinates_xyz[1][2] = corner_a[2];
209 target_partial_cell.coordinates_xyz[2][0] = corner_b[0];
210 target_partial_cell.coordinates_xyz[2][1] = corner_b[1];
211 target_partial_cell.coordinates_xyz[2][2] = corner_b[2];
212
213 // clip the current target cell triangle with all source cells
214 yac_cell_clipping(N, source_cell, target_partial_cell, overlap_buffer);
215
216 // for all source cells
217 for (size_t n = 0; n < N; n++) {
218
219 if (overlap_buffer[n].num_corners == 0) continue;
220
221 if (overlap_barycenters == NULL)
222 overlap_areas[n] +=
223 yac_huiliers_area(overlap_buffer[n]) * edge_direction;
224 else
225 overlap_areas[n] +=
227 overlap_buffer[n], overlap_barycenters[n], edge_direction);
228 }
229 }
230
231 for (size_t n = 0; n < N; n++) {
232
233 if (overlap_areas[n] < 0.0) {
234
235 overlap_areas[n] = -overlap_areas[n];
236
237 if (overlap_barycenters != NULL) {
238 overlap_barycenters[n][0] = -overlap_barycenters[n][0];
239 overlap_barycenters[n][1] = -overlap_barycenters[n][1];
240 overlap_barycenters[n][2] = -overlap_barycenters[n][2];
241 }
242 }
243 }
244
245 if (overlap_barycenters != NULL)
246 for (size_t n = 0; n < N; n++)
247 if ((overlap_areas[n] > 0.0) &&
248 ((overlap_barycenters[n][0] != 0.0) ||
249 (overlap_barycenters[n][1] != 0.0) ||
250 (overlap_barycenters[n][2] != 0.0)))
251 normalise_vector(overlap_barycenters[n]);
252
253#ifdef YAC_VERBOSE_CLIPPING
254 for (size_t n = 0; n < N; n++)
255 printf("overlap area %zu: %lf \n", n, overlap_areas[n]);
256#endif
257}
258
259static inline double dotproduct(double a[], double b[]) {
260
261 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
262}
263
264/* ------------------------- */
265
267 struct yac_grid_cell * source_cell,
268 struct yac_grid_cell target_cell,
269 double * partial_areas) {
270
272 N, source_cell, target_cell, partial_areas, NULL);
273}
274
275/* ------------------------- */
276
277static enum yac_cell_type get_cell_type(struct yac_grid_cell cell) {
278
279 if (cell.num_corners == 0) return YAC_GREAT_CIRCLE_CELL;
280
281 enum yac_cell_type cell_type = YAC_MIXED_CELL;
282
283 // if the cell is a typical lon-lat cell
284 if ((cell.num_corners == 4) &&
285 ((cell.edge_type[0] == YAC_LAT_CIRCLE_EDGE &&
286 cell.edge_type[1] == YAC_LON_CIRCLE_EDGE &&
287 cell.edge_type[2] == YAC_LAT_CIRCLE_EDGE &&
288 cell.edge_type[3] == YAC_LON_CIRCLE_EDGE) ||
289 (cell.edge_type[0] == YAC_LON_CIRCLE_EDGE &&
290 cell.edge_type[1] == YAC_LAT_CIRCLE_EDGE &&
291 cell.edge_type[2] == YAC_LON_CIRCLE_EDGE &&
292 cell.edge_type[3] == YAC_LAT_CIRCLE_EDGE))) {
293
294 cell_type = YAC_LON_LAT_CELL;
295
296 } else {
297
298 size_t count_lat_edges = 0, count_great_circle_edges = 0;
299
300 // count the number of edges for each type
301 // (lon edges are counted as gc ones)
302 for (size_t i = 0; i < cell.num_corners; ++i)
303 if (cell.edge_type[i] == YAC_LON_CIRCLE_EDGE ||
305 count_great_circle_edges++;
306 else
307 count_lat_edges++;
308
309 // if the cell is a lon lat cell with one lat edge having a length of zero
310 // due to being at a pole
311 if ((count_lat_edges == 1) && (count_great_circle_edges == 2)) {
312
313 size_t i;
314 for (i = 0; i < 3; ++i) if (cell.edge_type[i] == YAC_LAT_CIRCLE_EDGE) break;
315 size_t pol_index = (3 + i - 1) % 3;
316 if (fabs(fabs(cell.coordinates_xyz[pol_index][2])-1.0) <
318
319 // if the cell only consists of great circle edges
320 } else if (count_lat_edges == 0) cell_type = YAC_GREAT_CIRCLE_CELL;
321
322 // if the cell only consists of lat circle edges
323 else if (count_great_circle_edges == 0) cell_type = YAC_LAT_CELL;
324 }
325
327 cell_type != YAC_MIXED_CELL,
328 "invalid cell type (cell contains edges consisting "
329 "of great circles and circles of latitude)")
330
331 return cell_type;
332}
333
334int yac_circle_compare(void const * a, void const * b) {
335
336 struct yac_circle const * circle_a = *(struct yac_circle const **)a;
337 struct yac_circle const * circle_b = *(struct yac_circle const **)b;
338
340 ((circle_a->type == GREAT_CIRCLE) ||
341 (circle_a->type == LON_CIRCLE) ||
342 (circle_a->type == LAT_CIRCLE)/* ||
343 (circle_a->type == POINT)*/) &&
344 ((circle_b->type == GREAT_CIRCLE) ||
345 (circle_b->type == LON_CIRCLE) ||
346 (circle_b->type == LAT_CIRCLE)/* ||
347 (circle_b->type == POINT)*/),
348 "ERROR(yac_circle_compare): unsupported circle type")
349
350 // put lat circles to the end
351 int ret = (circle_a->type == LAT_CIRCLE) - (circle_b->type == LAT_CIRCLE);
352 if (!ret) ret = (int)(circle_a->type) - (int)(circle_b->type);
353 if (!ret) {
354 switch (circle_a->type) {
355 default:
356 case (GREAT_CIRCLE):
357 for (int i = 0; !ret && (i < 3); ++i) {
358 ret =
359 (circle_a->data.gc.norm_vector[i] >
360 circle_b->data.gc.norm_vector[i]) -
361 (circle_a->data.gc.norm_vector[i] <
362 circle_b->data.gc.norm_vector[i]);
363 }
364 break;
365 case (LON_CIRCLE):
366 for (int i = 0; !ret && (i < 3); ++i) {
367 ret =
368 (circle_a->data.lon.norm_vector[i] >
369 circle_b->data.lon.norm_vector[i]) -
370 (circle_a->data.lon.norm_vector[i] <
371 circle_b->data.lon.norm_vector[i]);
372 }
373 break;
374 case (LAT_CIRCLE):
375 ret = circle_a->data.lat.north_is_out - circle_b->data.lat.north_is_out;
376 if (!ret)
377 ret = (circle_a->data.lat.z > circle_b->data.lat.z) -
378 (circle_a->data.lat.z < circle_b->data.lat.z);
379 break;
380 // case (POINT):
381 // for (int i = 0; !ret && (i < 3); ++i)
382 // ret = (circle_a->data.p.vec[i] > circle_b->data.p.vec[i]) -
383 // (circle_a->data.p.vec[i] < circle_b->data.p.vec[i]);
384 // break;
385 }
386 }
387
388 return ret;
389}
390
391
393 double const a[3], double const b[3], double norm_vector[3]) {
394
396
397 double scale = 1.0 / sqrt(norm_vector[0] * norm_vector[0] +
398 norm_vector[1] * norm_vector[1] +
399 norm_vector[2] * norm_vector[2]);
400 norm_vector[0] *= scale;
401 norm_vector[1] *= scale;
402 norm_vector[2] *= scale;
403}
404
406 double const * a, double const * b, enum yac_edge_type type,
407 int edge_ordering, struct yac_circle * circle) {
408
409 if (points_are_identically(a, b)) {
410 circle->type = POINT;
411 circle->data.p.vec[0] = a[0];
412 circle->data.p.vec[1] = a[1];
413 circle->data.p.vec[2] = a[2];
414 } else {
415
416 // switch edges to ensure that "inside/outside" of the circle is computed
417 // correctly
418 if (edge_ordering < 0) {
419 double const * temp = a;
420 a = b;
421 b = temp;
422 }
423
428 "ERROR(yac_circle_generate): invalid edge type")
429
430 switch (type) {
431 default:
433 circle->type = GREAT_CIRCLE;
435 break;
437 circle->type = LON_CIRCLE;
439 break;
441 circle->type = LAT_CIRCLE;
442 circle->data.lat.north_is_out = (a[0] * b[1] - a[1] * b[0]) < 0.0;
443 if (circle->data.lat.north_is_out == 1337) exit(EXIT_FAILURE);
444 circle->data.lat.z = a[2];
445 break;
446 }
447 }
448}
449
450static inline struct yac_circle
452
453 struct yac_circle circle = {
454 .type = LAT_CIRCLE,
455 .data.lat.north_is_out = north_is_out,
456 .data.lat.z = z
457 };
458 return circle;
459}
460
462 double const point[3], struct yac_circle * circle) {
463
465 (circle->type == GREAT_CIRCLE) ||
466 (circle->type == LON_CIRCLE) ||
467 (circle->type == LAT_CIRCLE)/* ||
468 (circle->type == POINT)*/,
469 "ERROR(yac_circle_point_is_inside): unsupported circle type")
470
471 switch (circle->type) {
472 default:
473 case (GREAT_CIRCLE): {
474 double dot = point[0] * circle->data.gc.norm_vector[0] +
475 point[1] * circle->data.gc.norm_vector[1] +
476 point[2] * circle->data.gc.norm_vector[2];
477 if (dot < (- yac_angle_tol * 1e-3)) return 0;
478 else if (dot > (+ yac_angle_tol * 1e-3)) return 1;
479 else return 2;
480 }
481 case (LON_CIRCLE): {
482 double dot = point[0] * circle->data.lon.norm_vector[0] +
483 point[1] * circle->data.lon.norm_vector[1] +
484 point[2] * circle->data.lon.norm_vector[2];
485 if (dot < (- yac_angle_tol * 1e-3)) return 0;
486 else if (dot > (+ yac_angle_tol * 1e-3)) return 1;
487 else return 2;
488 }
489 case (LAT_CIRCLE): {
490 double diff_z = circle->data.lat.z - point[2];
491 if (fabs(diff_z) < yac_angle_tol * 1e-3) return 2;
492 else return (diff_z < 0.0) ^ circle->data.lat.north_is_out;
493 }
494 // case (POINT):
495 // return points_are_identically(point, circle->data.p.vec)?2:0;
496 }
497}
498
500 double const a[3], double const b[3], struct yac_circle * circle) {
501
503 (circle->type == GREAT_CIRCLE) ||
504 (circle->type == LON_CIRCLE) ||
505 (circle->type == LAT_CIRCLE) ||
506 (circle->type == POINT),
507 "ERROR(yac_circle_compare_distances): invalid circle type")
508
509 double dist_a, dist_b;
510 switch (circle->type) {
511 default:
512 case(GREAT_CIRCLE): {
513 double norm_vector[3] = {circle->data.gc.norm_vector[0],
514 circle->data.gc.norm_vector[1],
515 circle->data.gc.norm_vector[2]};
516 dist_a = fabs(a[0] * norm_vector[0] +
517 a[1] * norm_vector[1] +
518 a[2] * norm_vector[2]);
519 dist_b = fabs(b[0] * norm_vector[0] +
520 b[1] * norm_vector[1] +
521 b[2] * norm_vector[2]);
522 break;
523 }
524 case(LON_CIRCLE): {
525 double norm_vector[3] = {circle->data.lon.norm_vector[0],
526 circle->data.lon.norm_vector[1],
527 circle->data.lon.norm_vector[2]};
528 dist_a = fabs(a[0] * norm_vector[0] +
529 a[1] * norm_vector[1] +
530 a[2] * norm_vector[2]);
531 dist_b = fabs(b[0] * norm_vector[0] +
532 b[1] * norm_vector[1] +
533 b[2] * norm_vector[2]);
534 break;
535 }
536 case(LAT_CIRCLE): {
537 double circle_lat = acos(circle->data.lat.z);
538 dist_a = fabs(circle_lat - acos(a[2]));
539 dist_b = fabs(circle_lat - acos(b[2]));
540 break;
541 }
542 case(POINT): {
543 dist_a = get_vector_angle(circle->data.p.vec, a);
544 dist_b = get_vector_angle(circle->data.p.vec, b);
545 break;
546 }
547 }
548
549 return (dist_a > dist_b + yac_angle_tol) -
550 (dist_a + yac_angle_tol < dist_b);
551}
552
554
556 (circle->type == GREAT_CIRCLE) ||
557 (circle->type == LAT_CIRCLE) ||
558 (circle->type == LON_CIRCLE),
559 "ERROR(yac_circle_contains_north_pole): circle type has to be either "
560 "GREAT_CIRCLE or LAT_CIRCLE")
561
562 switch (circle->type) {
563 default:
564 case (GREAT_CIRCLE):
565 return circle->data.gc.norm_vector[2] > 0.0;
566 case (LAT_CIRCLE):
567 return !circle->data.lat.north_is_out;
568 case (LON_CIRCLE):
569 return 0;
570 }
571}
572
576static void circle_clipping(
577 struct point_list * cell, size_t num_cell_edges,
578 struct yac_circle ** clipping_circles, size_t num_clipping_circles) {
579
580 // to avoid some problems that can occur close the the pole, we process the
581 // target lat-circle edges at the end
582 qsort(clipping_circles, num_clipping_circles, sizeof(*clipping_circles),
584
585 // for all clipping circles
586 for (size_t i = 0; (i < num_clipping_circles) && (num_cell_edges > 1); ++i) {
587
588 struct point_list_element * cell_edge_start = cell->start;
589 struct point_list_element * cell_edge_end = cell_edge_start->next;
590
591 struct yac_circle * clipping_circle = clipping_circles[i];
592
593 int start_is_inside, first_start_is_inside;
594
595 start_is_inside =
597 cell_edge_start->vec_coords, clipping_circle);
598 first_start_is_inside = start_is_inside;
599
600 // for all edges of the cell
601 for (size_t cell_edge_idx = 0; cell_edge_idx < num_cell_edges;
602 ++cell_edge_idx) {
603
604 int end_is_inside =
605 (cell_edge_idx != num_cell_edges - 1)?
607 cell_edge_end->vec_coords, clipping_circle)):first_start_is_inside;
608
609 double * cell_edge_coords[2] =
610 {cell_edge_start->vec_coords, cell_edge_end->vec_coords};
611 struct yac_circle * cell_edge_circle = cell_edge_start->edge_circle;
612
613 double intersection[2][3];
614 int num_edge_intersections = -1;
615
616 // One intersection between the clipping circle and the cell edge is
617 // possible, if either the start or the end vertex of the current cell
618 // edge is inside and the respective other one is outside.
619 // Two intersections are possible, if either the clipping circle or the
620 // cell edge is a circle of latitude, while this other is not.
621 // Additionally, this is not possible if one of the two cell edge
622 // vertices is inside while the other is not or if both cell edge
623 // vertices are on the plane of the clipping edge.
624 int one_intersect_expected = (end_is_inside + start_is_inside == 1);
625 int two_intersect_possible =
626 ((cell_edge_circle->type == LAT_CIRCLE) ^
627 (clipping_circle->type == LAT_CIRCLE)) &&
628 (end_is_inside + start_is_inside != 1) &&
629 (end_is_inside + start_is_inside != 4);
630
631 // if we need to check for intersections
632 if (one_intersect_expected || two_intersect_possible) {
633
634 // get intersection points between the clipping circle and the circle
635 // of the cell edge
636 int num_circle_intersections =
638 *clipping_circle, *cell_edge_circle,
639 intersection[0], intersection[1]);
640
642 (num_circle_intersections >= -1) && (num_circle_intersections <= 2),
643 "ERROR(circle_clipping): invalid number of intersections")
644
645 // determine the intersections of the clipping circle with the cell
646 // edge based on the circle intersections
647 switch (num_circle_intersections) {
648
649 // special case:
650 // both circles are on the same plane
651 case (-1): {
652
653 // MoHa: this part of the code should not be reachable...but I
654 // leave it just in case...
655 start_is_inside = 2;
656 end_is_inside = 2;
657 num_edge_intersections = 0;
658
659 one_intersect_expected = 0;
660
661 // if this is the first iteration
662 if (cell_edge_idx == 0) first_start_is_inside = 2;
663 break;
664 }
665 // special case:
666 // no intersections between the two circles
667 // (can occure if one of the two circles is a latitude circle while
668 // the other is a great circle)
669 case (0): {
670 num_edge_intersections = 0;
671 break;
672 }
673 // standart case:
674 // two intersections between the two circles
675 // (one intersection between two circles is a special case, but is
676 // being handled here anyway)
677 default:
678 case (1):
679 case (2): {
680
681 // check the relation between the intersection points and the
682 // cell edge
683 int is_on_edge[2] = {
685 intersection[0], cell_edge_coords[0], cell_edge_coords[1],
686 cell_edge_circle->type),
687 (num_circle_intersections == 2)?
689 intersection[1], cell_edge_coords[0], cell_edge_coords[1],
690 cell_edge_circle->type):0};
691
692 // if both intersection points are on the edge
693 if (is_on_edge[0] && is_on_edge[1]) {
694
695 // if only one intersection was expected
697 !one_intersect_expected,
698 "ERROR: two intersections found, even "
699 "though no circle of latitude involed and "
700 "both edge vertices on different sides of "
701 "the cutting plane.\n"
702 "cell edge (%lf %lf %lf) (%lf %lf %lf) (edge type %d)\n"
703 "circle (gc: norm_vec %lf %lf %lf\n"
704 " lon: norm_vec %lf %lf %lf\n"
705 " lat: z %lf north_is_out %d\n"
706 " point: vec %lf %lf %lf) (circle type %d)\n"
707 "intersections points (%lf %lf %lf) (%lf %lf %lf)\n",
708 cell_edge_coords[0][0],
709 cell_edge_coords[0][1],
710 cell_edge_coords[0][2],
711 cell_edge_coords[1][0],
712 cell_edge_coords[1][1],
713 cell_edge_coords[1][2], (int)cell_edge_circle->type,
714 clipping_circle->data.gc.norm_vector[0],
715 clipping_circle->data.gc.norm_vector[1],
716 clipping_circle->data.gc.norm_vector[2],
717 clipping_circle->data.lon.norm_vector[0],
718 clipping_circle->data.lon.norm_vector[1],
719 clipping_circle->data.lon.norm_vector[2],
720 clipping_circle->data.lat.z,
721 clipping_circle->data.lat.north_is_out,
722 clipping_circle->data.p.vec[0],
723 clipping_circle->data.p.vec[1],
724 clipping_circle->data.p.vec[2], (int)(clipping_circle->type),
725 intersection[0][0], intersection[0][1], intersection[0][2],
726 intersection[1][0], intersection[1][1], intersection[1][2])
727
728 // if both cell edge vertices are on the same side of the
729 // clipping edge
730 if (end_is_inside == start_is_inside) {
731
732 // if the two intersection points are basically the same point
733 if (points_are_identically(intersection[0], intersection[1])) {
734
735 num_edge_intersections = 1;
736
737 } else {
738
739 // if the first intersection point is further away from the
740 // cell start edge vertex than the second one ->
741 // switch them
742 if (sq_len_diff_vec(cell_edge_coords[0], intersection[0]) >
743 sq_len_diff_vec(cell_edge_coords[0], intersection[1])) {
744
745 double temp_intersection[3];
746 temp_intersection[0] = intersection[1][0];
747 temp_intersection[1] = intersection[1][1];
748 temp_intersection[2] = intersection[1][2];
749 intersection[1][0] = intersection[0][0];
750 intersection[1][1] = intersection[0][1];
751 intersection[1][2] = intersection[0][2];
752 intersection[0][0] = temp_intersection[0];
753 intersection[0][1] = temp_intersection[1];
754 intersection[0][2] = temp_intersection[2];
755 }
756
757 num_edge_intersections = 2;
758 }
759
760 // one of the two cell edge vertices is on the clipping circle
761 // => one of the two intersection points must be more or less
762 // identical to a vertex of the cell edge
763 } else {
764
765 double * cell_edge_coord = cell_edge_coords[end_is_inside == 2];
766 double distances[2] = {
767 sq_len_diff_vec(cell_edge_coord, intersection[0]),
768 sq_len_diff_vec(cell_edge_coord, intersection[1])};
769
770 // if both intersection points are nearly identical to the cell
771 // edge vertex
772 if ((distances[0] < yac_sq_angle_tol) &&
773 (distances[1] < yac_sq_angle_tol)) {
774
775 num_edge_intersections = 0;
776
777 } else {
778
779 num_edge_intersections = 1;
780
781 // if the fist intersection points is closer than the second
782 if (distances[0] < distances[1]) {
783 intersection[0][0] = intersection[1][0];
784 intersection[0][1] = intersection[1][1];
785 intersection[0][2] = intersection[1][2];
786 }
787 }
788 }
789 // if only one intersection point is on the cell edge
790 } else if (is_on_edge[0] || is_on_edge[1]) {
791
792 // if one of the two cell edge vertices is on the clipping edge
793 if ((end_is_inside == 2) || (start_is_inside == 2)) {
794
795 // we assume that the intersection point is the respective
796 // cell edge vertex, which is on the clipping edge -> we do not
797 // need this intersection point
798 num_edge_intersections = 0;
799
800 } else {
801
802 // if the second intersection point is on the cell edge
803 if (is_on_edge[1]) {
804 intersection[0][0] = intersection[1][0];
805 intersection[0][1] = intersection[1][1];
806 intersection[0][2] = intersection[1][2];
807 }
808
809 num_edge_intersections = 1;
810 }
811
812 // if none of the two intersection points is on the cell edge
813 } else {
814 num_edge_intersections = 0;
815 }
816 break;
817 } // case one or two circle intersections
818 } // switch (num_circle_intersections)
819
820 // If an intersection was expected but none was found, the clipping
821 // circle was most propably to close to an edge vertex. We can now
822 // assume that the respective vertex is directly on the clipping circle.
823 if ((one_intersect_expected) && (num_edge_intersections == 0)) {
824
826 cell_edge_coords[0], cell_edge_coords[1],
827 clipping_circle) <= 0)
828 start_is_inside = 2;
829 else
830 end_is_inside = 2;
831
832 one_intersect_expected = 0;
833
834 // if this is the first iteration
835 if (cell_edge_idx == 0) first_start_is_inside = start_is_inside;
836 }
837 // else -> not edge intersection was excepted
838 } else {
839 num_edge_intersections = 0;
840 }
841
843 num_edge_intersections != -1,
844 "ERROR(circle_clipping): internal error");
845
846 // here we know the number of intersections and their location and we
847 // know the relation of the cell edge vertices to the clipping edge
848 // (start_is_inside and end_is_inside)
849
850 // if the start cell edge vertex is outside -> dump it after clipping
851 // is finished
852 cell_edge_start->to_be_removed = start_is_inside == 0;
853
854 // the easiest case is that we expected one intersection and got one
855 if (one_intersect_expected) {
856
857 // insert an intersection point in the cell point list in the
858 // current edge
859 struct point_list_element * intersect_point =
861 cell_edge_start->next = intersect_point;
862 intersect_point->next = cell_edge_end;
863
864 intersect_point->vec_coords[0] = intersection[0][0];
865 intersect_point->vec_coords[1] = intersection[0][1];
866 intersect_point->vec_coords[2] = intersection[0][2];
867
868 intersect_point->edge_circle =
869 (start_is_inside)?clipping_circle:cell_edge_circle;
870
871 // if the clipping edge goes through both of the two cell edge vertices
872 } else if ((start_is_inside == 2) && (end_is_inside == 2)) {
873
874 // if one of the two edges is a circle of latitude while the other is
875 // not, then we may have to replace the cell edge circle with the
876 // clipping circle
877 if ((cell_edge_circle->type == LAT_CIRCLE) ^
878 (clipping_circle->type == LAT_CIRCLE)) {
879
880 int clipping_circle_contains_north =
881 yac_circle_contains_north_pole(clipping_circle);
882 int same_inside_direction =
883 clipping_circle_contains_north ==
884 yac_circle_contains_north_pole(cell_edge_circle);
885 int cell_edge_is_on_south_hemisphere =
886 (cell_edge_end->vec_coords[2] < 0.0);
887 int clipping_circle_is_lat = clipping_circle->type == LAT_CIRCLE;
888
889 // after generating a truth table with these four logical values
890 // I came up with the following formula to determine whether
891 // the cell edge type should switch to the one of the clipping circle
892 // or not
893 if (same_inside_direction &&
894 (cell_edge_is_on_south_hemisphere ^
895 clipping_circle_contains_north ^
896 clipping_circle_is_lat))
897 cell_edge_start->edge_circle = clipping_circle;
898
899 }
900
901 // if we have no intersection, but the cell edge start vertex
902 // is on the clipping edge
903 } else if ((num_edge_intersections == 0) && (start_is_inside == 2)) {
904
905 // if the end cell edge vertex is outside
906 if (end_is_inside == 0)
907 cell_edge_start->edge_circle = clipping_circle;
908
909 // if we have two intersections (only happens if one of the two edges is
910 // a circle of latitude while the other is not)
911 } else if (num_edge_intersections == 2) {
912
913 struct point_list_element * intersect_points[2] =
916
917 // add two points between the current source edge vertices
918 cell_edge_start->next = intersect_points[0];
919 intersect_points[0]->next = intersect_points[1];
920 intersect_points[1]->next = cell_edge_end;
921
922 intersect_points[0]->vec_coords[0] = intersection[0][0];
923 intersect_points[0]->vec_coords[1] = intersection[0][1];
924 intersect_points[0]->vec_coords[2] = intersection[0][2];
925 intersect_points[1]->vec_coords[0] = intersection[1][0];
926 intersect_points[1]->vec_coords[1] = intersection[1][1];
927 intersect_points[1]->vec_coords[2] = intersection[1][2];
928
930 ((start_is_inside == 0) && (end_is_inside == 0)) ||
931 ((start_is_inside == 1) && (end_is_inside == 1)),
932 "ERROR: one cell edge vertex is on the clipping edge, therefore we "
933 "should not have two intersections.")
934
935 // if a and b are outside
936 if ((start_is_inside == 0) && (end_is_inside == 0)) {
937 intersect_points[0]->edge_circle = cell_edge_circle;
938 intersect_points[1]->edge_circle = clipping_circle;
939
940 // if a and b are inside
941 } else /*if ((start_is_inside == 1) && (end_is_inside == 1))*/ {
942 intersect_points[0]->edge_circle = clipping_circle;
943 intersect_points[1]->edge_circle = cell_edge_circle;
944 }
945
946 // if we have one intersection even though the two cell edge vertices
947 // are not on opposite sides of the clipping edge
948 } else if (two_intersect_possible && (num_edge_intersections == 1)) {
949
950 // ensure that both cell edge vertices are on the same side of the
951 // clipping edge or that at least one cell vertex is directly on
952 // the clipping edge
954 (start_is_inside == end_is_inside) ||
955 ((start_is_inside == 2) || (end_is_inside == 2)),
956 "ERROR: unhandled intersection case")
957
958 switch (MAX(start_is_inside, end_is_inside)) {
959
960 // if both cell edge vertices are outside -> circle of latitude and
961 // greate circle touch at a single vertex
962 default:
963 case (0): {
964 // insert an intersection point in the cell point list in the
965 // current edge
966 struct point_list_element * intersect_point =
968 cell_edge_start->next = intersect_point;
969 intersect_point->next = cell_edge_end;
970
971 intersect_point->vec_coords[0] = intersection[0][0];
972 intersect_point->vec_coords[1] = intersection[0][1];
973 intersect_point->vec_coords[2] = intersection[0][2];
974
975 intersect_point->edge_circle = clipping_circle;
976 break;
977 }
978
979 // if both cell edge vertices are inside -> nothing to be done
980 case (1): break;
981
982 // if one of the two cell edge vertices is on the clipping edge
983 // while the other is either inside or outside
984 case (2): {
985 // insert an intersection point in the cell point list in the
986 // current edge
987 struct point_list_element * intersect_point =
989 cell_edge_start->next = intersect_point;
990 intersect_point->next = cell_edge_end;
991
992 intersect_point->vec_coords[0] = intersection[0][0];
993 intersect_point->vec_coords[1] = intersection[0][1];
994 intersect_point->vec_coords[2] = intersection[0][2];
995
996 // if the start cell edge vertex is on the clipping edge
997 if (start_is_inside == 2) {
998 // if the end cell edge vertex in on the outside
999 if (end_is_inside == 0) {
1000 // cell_edge_start->edge_circle = cell_edge_circle;
1001 intersect_point->edge_circle = clipping_circle;
1002 } else {
1003 cell_edge_start->edge_circle = clipping_circle;
1004 intersect_point->edge_circle = cell_edge_circle;
1005 }
1006 // if the end cell edge vertex is on the clipping edge
1007 } else {
1008 // if the start cell edge vertex in on the outside
1009 if (start_is_inside == 0) {
1010 cell_edge_start->edge_circle = clipping_circle;
1011 intersect_point->edge_circle = cell_edge_circle;
1012 } else {
1013 // cell_edge_start->edge_circle = cell_edge_circle;
1014 intersect_point->edge_circle = clipping_circle;
1015 }
1016 }
1017 }
1018 }
1019 }
1020
1021 cell_edge_start = cell_edge_end;
1022 cell_edge_end = cell_edge_end->next;
1023 start_is_inside = end_is_inside;
1024
1025 } // for all cell edges
1026
1027 // remove all points that are to be deleted
1028 num_cell_edges = remove_points(cell);
1029 }
1030}
1031
1033 struct point_list * source_list, struct point_list target_list, size_t nct) {
1034
1035 struct yac_circle * clipping_circles[nct];
1036
1037 struct point_list_element * curr_tgt_point = target_list.start;
1038
1039 for (size_t i = 0; i < nct; ++i, curr_tgt_point = curr_tgt_point->next)
1040 clipping_circles[i] = curr_tgt_point->edge_circle;
1041
1042 YAC_ASSERT(
1043 source_list->start != NULL,
1044 "ERROR(point_list_clipping): source cell without corners")
1045
1046 // count the number of edges in the source cell
1047 size_t ncs = 1;
1048 for (struct point_list_element * curr = source_list->start->next,
1049 * end = source_list->start;
1050 curr != end; curr = curr->next, ++ncs);
1051
1052 circle_clipping(source_list, ncs, clipping_circles, nct);
1053}
1054
1055static void copy_point_list(struct point_list in, struct point_list * out) {
1056
1057 reset_point_list(out);
1058
1059 struct point_list_element * curr = in.start;
1060
1061 if (curr == NULL) return;
1062
1063 struct point_list_element * new_point_list = get_free_point_list_element(out);
1064 out->start = new_point_list;
1065 *new_point_list = *curr;
1066 curr = curr->next;
1067
1068 do {
1069
1070 new_point_list->next = get_free_point_list_element(out);
1071 new_point_list = new_point_list->next;
1072 *new_point_list = *curr;
1073 curr = curr->next;
1074
1075 } while (curr != in.start);
1076
1077 new_point_list->next = out->start;
1078}
1079
1080void yac_cell_clipping (size_t N,
1081 struct yac_grid_cell * source_cell,
1082 struct yac_grid_cell target_cell,
1083 struct yac_grid_cell * overlap_buffer) {
1084
1085 if (target_cell.num_corners < 2) {
1086 for (size_t n = 0; n < N; n++ ) overlap_buffer[n].num_corners = 0;
1087 return;
1088 }
1089
1090 enum yac_cell_type tgt_cell_type = get_cell_type(target_cell);
1091
1092 YAC_ASSERT(
1093 tgt_cell_type != YAC_MIXED_CELL,
1094 "invalid target cell type (cell contains edges consisting "
1095 "of great circles and circles of latitude)")
1096
1097 // determine ordering of target cell corners
1098 int target_ordering = get_cell_points_ordering(target_cell);
1099 // if all corners of the target cell are on the same great circle
1100 if (!target_ordering) {
1101 for (size_t n = 0; n < N; n++ ) overlap_buffer[n].num_corners = 0;
1102 return;
1103 }
1104
1105 size_t max_num_src_cell_corners = 0;
1106 for (size_t n = 0; n < N; ++n)
1107 if (source_cell[n].num_corners > max_num_src_cell_corners)
1108 max_num_src_cell_corners = source_cell[n].num_corners;
1109 struct yac_circle * circle_buffer =
1110 xmalloc(
1111 (target_cell.num_corners + max_num_src_cell_corners) *
1112 sizeof(*circle_buffer));
1113 struct yac_circle * src_circle_buffer =
1114 circle_buffer + target_cell.num_corners;
1115
1116 // generate point list for target cell (clip cell)
1117 struct point_list target_list;
1118 init_point_list(&target_list);
1119 size_t nct =
1121 &target_list, target_cell, target_ordering, circle_buffer);
1122
1123 struct point_list source_list, temp_list;
1124 init_point_list(&temp_list);
1125 init_point_list(&source_list);
1126
1127 // for all source cells
1128 for (size_t n = 0; n < N; n++ ) {
1129
1130 overlap_buffer[n].num_corners = 0;
1131
1132 enum yac_cell_type src_cell_type = get_cell_type(source_cell[n]);
1133
1134 YAC_ASSERT(
1135 src_cell_type != YAC_MIXED_CELL,
1136 "invalid source cell type (cell contains edges consisting "
1137 "of great circles and circles of latitude)")
1138
1139 if (source_cell[n].num_corners < 2) continue;
1140
1141 // determine ordering of source cell corners
1142 int source_ordering = get_cell_points_ordering(source_cell[n]);
1143
1144 // if all corners of the source cell are on the same great circle
1145 if (!source_ordering) continue;
1146
1147 // generate point list for current source list
1148 size_t ncs =
1150 &source_list, source_cell[n], source_ordering, src_circle_buffer);
1151
1152 struct point_list * overlap;
1153 double fabs_tgt_coordinate_z = fabs(target_cell.coordinates_xyz[0][2]);
1154 double fabs_src_coordinate_z = fabs(source_cell[n].coordinates_xyz[0][2]);
1155
1156 // in the case that source and target cell are both YAC_LAT_CELL's, than the
1157 // bigger one has to be the target cell
1158 // a similar problem occurs when the target cell is a YAC_LAT_CELL and the
1159 // source is a YAC_GREAT_CIRCLE_CELL which overlaps with the pole that is also
1160 // include in the target cell
1161 if (((tgt_cell_type == YAC_LAT_CELL) && (src_cell_type == YAC_GREAT_CIRCLE_CELL)) ||
1162 ((tgt_cell_type == YAC_LAT_CELL) && (src_cell_type == YAC_LAT_CELL) &&
1163 (fabs_tgt_coordinate_z > fabs_src_coordinate_z))) {
1164
1165 copy_point_list(target_list, &temp_list);
1166
1167 point_list_clipping(&temp_list, source_list, ncs);
1168
1169 overlap = &temp_list;
1170
1171 } else {
1172
1173 point_list_clipping(&source_list, target_list, nct);
1174
1175 overlap = &source_list;
1176 }
1177
1178 generate_cell(overlap, overlap_buffer + n);
1179
1180 }
1181 free(circle_buffer);
1182 free_point_list(&source_list);
1183 free_point_list(&target_list);
1184 free_point_list(&temp_list);
1185}
1186
1188 struct yac_grid_cell * cell, double z_upper_bound, double z_lower_bound,
1189 struct yac_grid_cell * overlap_buffer) {
1190
1191 double cell_upper_bound = cell->coordinates_xyz[0][2];
1192 double cell_lower_bound = cell->coordinates_xyz[2][2];
1193
1194 int upper_idx[2], lower_idx[2];
1195
1196 if (cell_upper_bound < cell_lower_bound) {
1197 double temp = cell_upper_bound;
1198 cell_upper_bound = cell_lower_bound;
1199 cell_lower_bound = temp;
1200 if (cell->edge_type[0] == YAC_LAT_CIRCLE_EDGE) {
1201 upper_idx[0] = 2;
1202 upper_idx[1] = 3;
1203 lower_idx[0] = 1;
1204 lower_idx[1] = 0;
1205 } else {
1206 upper_idx[0] = 2;
1207 upper_idx[1] = 1;
1208 lower_idx[0] = 3;
1209 lower_idx[1] = 0;
1210 }
1211 } else {
1212 if (cell->edge_type[0] == YAC_LAT_CIRCLE_EDGE) {
1213 upper_idx[0] = 0;
1214 upper_idx[1] = 1;
1215 lower_idx[0] = 3;
1216 lower_idx[1] = 2;
1217 } else {
1218 upper_idx[0] = 0;
1219 upper_idx[1] = 3;
1220 lower_idx[0] = 1;
1221 lower_idx[1] = 2;
1222 }
1223 }
1224
1225 // if z_upper_bound and z_lower_bound are identical or
1226 // if cell does not overlap with latitude band
1227 if ((z_upper_bound == z_lower_bound) ||
1228 (cell_upper_bound <= z_lower_bound) ||
1229 (cell_lower_bound >= z_upper_bound)) {
1230 overlap_buffer->num_corners = 0;
1231 return;
1232 }
1233
1234 if (overlap_buffer->array_size < 4) {
1235 overlap_buffer->coordinates_xyz =
1236 xmalloc(4 * sizeof(*(overlap_buffer->coordinates_xyz)));
1237 overlap_buffer->edge_type =
1238 xmalloc(4 * sizeof(*(overlap_buffer->edge_type)));
1239 overlap_buffer->array_size = 4;
1240 }
1241
1242 memcpy(overlap_buffer->coordinates_xyz, cell->coordinates_xyz,
1243 4 * sizeof(*(cell->coordinates_xyz)));
1244 memcpy(overlap_buffer->edge_type, cell->edge_type,
1245 4 * sizeof(*(cell->edge_type)));
1246 overlap_buffer->num_corners = 4;
1247
1248
1249 double tmp_scale;
1250 double * p[2];
1251 if (fabs(cell_lower_bound) < fabs(cell_upper_bound)) {
1252 p[0] = cell->coordinates_xyz[lower_idx[0]];
1253 p[1] = cell->coordinates_xyz[lower_idx[1]];
1254 tmp_scale = cell->coordinates_xyz[lower_idx[0]][0] *
1255 cell->coordinates_xyz[lower_idx[0]][0] +
1256 cell->coordinates_xyz[lower_idx[0]][1] *
1257 cell->coordinates_xyz[lower_idx[0]][1];
1258 } else {
1259 p[0] = cell->coordinates_xyz[upper_idx[0]];
1260 p[1] = cell->coordinates_xyz[upper_idx[1]];
1261 tmp_scale = cell->coordinates_xyz[upper_idx[0]][0] *
1262 cell->coordinates_xyz[upper_idx[0]][0] +
1263 cell->coordinates_xyz[upper_idx[0]][1] *
1264 cell->coordinates_xyz[upper_idx[0]][1];
1265 }
1266
1267 // if upper bound overlaps with cell
1268 if ((z_upper_bound < cell_upper_bound) &&
1269 (z_upper_bound > cell_lower_bound)) {
1270
1271 double scale = sqrt((1.0 - z_upper_bound * z_upper_bound) / tmp_scale);
1272
1273 overlap_buffer->coordinates_xyz[upper_idx[0]][0] = p[0][0] * scale;
1274 overlap_buffer->coordinates_xyz[upper_idx[0]][1] = p[0][1] * scale;
1275 overlap_buffer->coordinates_xyz[upper_idx[0]][2] = z_upper_bound;
1276 overlap_buffer->coordinates_xyz[upper_idx[1]][0] = p[1][0] * scale;
1277 overlap_buffer->coordinates_xyz[upper_idx[1]][1] = p[1][1] * scale;
1278 overlap_buffer->coordinates_xyz[upper_idx[1]][2] = z_upper_bound;
1279 }
1280
1281 // if lower bound overlaps with cell
1282 if ((z_lower_bound < cell_upper_bound) &&
1283 (z_lower_bound > cell_lower_bound)) {
1284
1285 double scale = sqrt((1.0 - z_lower_bound * z_lower_bound) / tmp_scale);
1286
1287 overlap_buffer->coordinates_xyz[lower_idx[0]][0] = p[0][0] * scale;
1288 overlap_buffer->coordinates_xyz[lower_idx[0]][1] = p[0][1] * scale;
1289 overlap_buffer->coordinates_xyz[lower_idx[0]][2] = z_lower_bound;
1290 overlap_buffer->coordinates_xyz[lower_idx[1]][0] = p[1][0] * scale;
1291 overlap_buffer->coordinates_xyz[lower_idx[1]][1] = p[1][1] * scale;
1292 overlap_buffer->coordinates_xyz[lower_idx[1]][2] = z_lower_bound;
1293 }
1294}
1295
1296static double get_closest_pole(struct point_list * cell_list) {
1297
1298 struct point_list_element * curr = cell_list->start;
1299 struct point_list_element * start = cell_list->start;
1300
1301 if (curr == NULL) return 1.0;
1302
1303 double max_z = 0.0;
1304
1305 do {
1306
1307 double curr_z = curr->vec_coords[2];
1308 if (fabs(curr_z) > fabs(max_z)) max_z = curr_z;
1309
1310 curr = curr->next;
1311 } while (curr != start);
1312
1313 return (max_z > 0.0)?1.0:-1.0;
1314}
1315
1316/*this routine is potentially being used in the CDO*/
1318 struct yac_grid_cell * cells,
1319 double lat_bounds[2], // lat in rad
1320 struct yac_grid_cell * overlap_buffer) {
1321
1322 double z_bounds[2] = {sin(lat_bounds[0]), sin(lat_bounds[1])};
1323 int upper_bound_idx = lat_bounds[0] < lat_bounds[1];
1324 int lower_bound_idx = upper_bound_idx ^ 1;
1325 int is_pole[2] = {fabs(fabs(lat_bounds[0]) - M_PI_2) < yac_angle_tol,
1326 fabs(fabs(lat_bounds[1]) - M_PI_2) < yac_angle_tol};
1327 int upper_is_north_pole =
1328 is_pole[upper_bound_idx] && (lat_bounds[upper_bound_idx] > 0.0);
1329 int lower_is_south_pole =
1330 is_pole[lower_bound_idx] && (lat_bounds[lower_bound_idx] < 0.0);
1331
1332 // if both bounds are nearly identical,
1333 // or if the upper bound is the south pole,
1334 // or if the lower bound is the north pole
1335 if ((fabs(lat_bounds[0] - lat_bounds[1]) < yac_angle_tol) ||
1336 (is_pole[upper_bound_idx] ^ upper_is_north_pole) ||
1337 (is_pole[lower_bound_idx] ^ lower_is_south_pole)) {
1338
1339 for (size_t n = 0; n < N; ++n)
1340 overlap_buffer[n].num_corners = 0;
1341 return;
1342 }
1343
1344 struct yac_circle lat_circle_buffer[2];
1345 struct yac_circle * lat_circles[2] =
1346 {&(lat_circle_buffer[0]), &(lat_circle_buffer[1])};
1347 size_t num_lat_circles = 0;
1348 if (!lower_is_south_pole)
1349 lat_circle_buffer[num_lat_circles++] =
1350 generate_lat_circle(z_bounds[lower_bound_idx], 0);
1351 if (!upper_is_north_pole)
1352 lat_circle_buffer[num_lat_circles++] =
1353 generate_lat_circle(z_bounds[upper_bound_idx], 1);
1354
1355 size_t max_num_cell_corners = 0;
1356 for (size_t n = 0; n < N; ++n)
1357 if (cells[n].num_corners > max_num_cell_corners)
1358 max_num_cell_corners = cells[n].num_corners;
1359 struct yac_circle * circle_buffer =
1360 xmalloc(max_num_cell_corners * sizeof(*circle_buffer));
1361
1362 struct point_list cell_list;
1363 init_point_list(&cell_list);
1364
1365 // for all source cells
1366 for (size_t n = 0; n < N; n++) {
1367
1368 if (cells[n].num_corners < 2) continue;
1369
1370 overlap_buffer[n].num_corners = 0;
1371
1372 enum yac_cell_type src_cell_type = get_cell_type(cells[n]);
1373
1374 YAC_ASSERT(
1375 src_cell_type != YAC_MIXED_CELL,
1376 "invalid source cell type (cell contains edges consisting "
1377 "of great circles and circles of latitude)\n")
1378
1379 if (src_cell_type == YAC_LON_LAT_CELL) {
1380
1382 cells + n, z_bounds[upper_bound_idx], z_bounds[lower_bound_idx],
1383 overlap_buffer + n);
1384
1385 } else {
1386
1387 int cell_ordering = get_cell_points_ordering(cells[n]);
1388
1389 // generate point list for current source list
1390 size_t num_corners =
1392 &cell_list, cells[n], cell_ordering, circle_buffer);
1393
1394 // if the cell contains a pole, we need to add this pole
1395 double pole = get_closest_pole(&cell_list);
1397 (double[3]){0.0, 0.0, pole}, cells[n])) {
1398
1399 int flag = 0;
1400
1401 // if the cell contains the north pole
1402 if (pole > 0.0) {
1403 // use lower bound if upper bound is north pole
1404 double z_bound = z_bounds[upper_bound_idx ^ upper_is_north_pole];
1405
1406 for (size_t i = 0; i < num_corners; ++i) {
1407 flag |= (z_bound < cells[n].coordinates_xyz[i][2]);
1408 }
1409 } else {
1410 // use upper bound if lower bound is south pole
1411 double z_bound = z_bounds[lower_bound_idx ^ lower_is_south_pole];
1412
1413 for (size_t i = 0; i < num_corners; ++i) {
1414 flag |= (z_bound > cells[n].coordinates_xyz[i][2]);
1415 }
1416 }
1417
1418 YAC_ASSERT(
1419 flag,
1420 "ERROR(yac_cell_lat_clipping): Latitude bounds are within a cell "
1421 "covering a pole, this is not supported. Increased grid resolution "
1422 "or widen lat bounds may help.")
1423 }
1424
1425 circle_clipping(&cell_list, num_corners, lat_circles, num_lat_circles);
1426
1427 generate_cell(&cell_list, overlap_buffer + n);
1428 }
1429 }
1430
1431 free(circle_buffer);
1432 free_point_list(&cell_list);
1433}
1434
1435/* ---------------------------------------------------- */
1436
1437void yac_correct_weights ( size_t nSourceCells, double * weight ) {
1438
1439 // number of iterations to get better accuracy of the weights
1440 enum {maxIter = 10};
1441
1442 for ( size_t iter = 1; iter < maxIter; iter++ ) {
1443
1444 double weight_sum = 0.0;
1445
1446 for ( size_t n = 0; n < nSourceCells; n++ ) weight_sum += weight[n];
1447
1448 if ( fabs(1.0 - weight_sum) < tol ) break;
1449
1450 double scale = 1.0 / weight_sum;
1451
1452 for ( size_t n = 0; n < nSourceCells; n++ ) weight[n] *= scale;
1453 }
1454}
1455
1456/* ---------------------------------------------------- */
1457
1459
1460 YAC_ASSERT(
1461 cell.num_corners > 2, "ERROR(get_cell_points_ordering): invalid cell")
1462
1463 double edge_norm_vectors[cell.num_corners][3];
1464 double vertices[cell.num_corners][3];
1465 size_t num_corners = 0;
1466
1467 for (size_t i = 0, j = cell.num_corners - 1; i < cell.num_corners; i++) {
1469 cell.coordinates_xyz[j], cell.coordinates_xyz[i])) {
1471 cell.coordinates_xyz[j], cell.coordinates_xyz[i],
1472 edge_norm_vectors[num_corners]);
1473 normalise_vector(edge_norm_vectors[num_corners]);
1474 vertices[num_corners][0] = cell.coordinates_xyz[i][0];
1475 vertices[num_corners][1] = cell.coordinates_xyz[i][1];
1476 vertices[num_corners][2] = cell.coordinates_xyz[i][2];
1477 ++num_corners;
1478 j = i;
1479 }
1480 }
1481
1482 struct sin_cos_angle angle_sum = SIN_COS_ZERO;
1483 int cell_direction = 0;
1484 int valid_angles_count = 0;
1485
1486 for (size_t i = 0, j = num_corners - 1; i < num_corners; j = i++) {
1487
1488 struct sin_cos_angle edge_angle =
1489 get_vector_angle_2(edge_norm_vectors[j], edge_norm_vectors[i]);
1490
1491 // if both edges are on the same plane
1492 if ((edge_angle.cos > 0.0) && (edge_angle.sin < yac_angle_tol))
1493 continue;
1494
1495 valid_angles_count++;
1496
1497 double edge_direction =
1498 edge_norm_vectors[j][0] * vertices[i][0] +
1499 edge_norm_vectors[j][1] * vertices[i][1] +
1500 edge_norm_vectors[j][2] * vertices[i][2];
1501
1502 struct sin_cos_angle temp = angle_sum;
1503 if (edge_direction < 0.0)
1504 cell_direction -= sub_angles(temp, edge_angle, &angle_sum);
1505 else
1506 cell_direction += sum_angles(temp, edge_angle, &angle_sum);
1507 }
1508
1509 if (valid_angles_count == 0) return 0;
1510 return (cell_direction >= 0)?1:-1;
1511}
1512
1513static void init_point_list(struct point_list * list) {
1514
1515 list->start = NULL;
1516 list->free_elements = NULL;
1517}
1518
1519static void reset_point_list(struct point_list * list) {
1520
1521 if (list->start != NULL) {
1522
1523 struct point_list_element * temp = list->start->next;
1524 list->start->next = list->free_elements;
1525 list->free_elements = temp;
1526
1527 list->start = NULL;
1528 }
1529}
1530
1531static size_t remove_points(struct point_list * list) {
1532
1533 struct point_list_element * curr = list->start;
1534 struct point_list_element * start = list->start;
1535 struct point_list_element * prev = NULL;
1536
1537 if (curr == NULL) return 0;
1538
1539 // find the first point that is not to be removed
1540 while(curr->to_be_removed) {
1541 prev = curr;
1542 curr = curr->next;
1543 if (curr == start) break;
1544 };
1545
1546 // there is no point to remain
1547 if (curr->to_be_removed) {
1548 reset_point_list(list);
1549 return 0;
1550 }
1551
1552 list->start = curr;
1553 start = curr;
1554 size_t num_remaining_points = 1;
1555
1556 prev = curr;
1557 curr = curr->next;
1558
1559 while (curr != start) {
1560
1561 if (curr->to_be_removed) {
1562 prev->next = curr->next;
1563 curr->next = list->free_elements;
1564 list->free_elements = curr;
1565 curr = prev;
1566 } else {
1567 num_remaining_points++;
1568 }
1569
1570 prev = curr;
1571 curr = curr->next;
1572
1573 };
1574
1575 if (list->start == list->start->next) {
1576
1577 list->start->next = list->free_elements;
1578 list->free_elements = list->start;
1579 list->start = NULL;
1580 num_remaining_points = 0;
1581 }
1582
1583 return num_remaining_points;
1584}
1585
1587static size_t remove_zero_length_edges(struct point_list * list) {
1588
1589 struct point_list_element * curr = list->start;
1590 struct point_list_element * start = list->start;
1591
1592 do {
1593
1594 if (!curr->to_be_removed) {
1595
1596 // if both points are nearly identical (angle between them is very small)
1597 if (points_are_identically(curr->vec_coords, curr->next->vec_coords)) {
1598 curr->to_be_removed = 1;
1599
1600 // check whether we can merge two lat circle edges
1601 } else if (curr->edge_circle->type == LAT_CIRCLE &&
1602 curr->next->edge_circle->type == LAT_CIRCLE) {
1603 double temp_a = atan2(curr->vec_coords[1], curr->vec_coords[0]);
1604 double temp_b = atan2(curr->next->next->vec_coords[1],
1605 curr->next->next->vec_coords[0]);
1606 if (fabs(get_angle(temp_a, temp_b)) < M_PI_2)
1607 curr->next->to_be_removed = 1;
1608 }
1609 }
1610
1611 curr = curr->next;
1612 } while (curr != start);
1613
1614 return remove_points(list);
1615}
1616
1618 struct point_list * list, struct yac_grid_cell cell, int cell_ordering,
1619 struct yac_circle * circle_buffer) {
1620
1621 reset_point_list(list);
1622
1623 YAC_ASSERT(
1624 cell.num_corners >= 2,
1625 "ERROR(generate_point_list): too few corners in cell");
1626
1627 struct point_list_element * curr = get_free_point_list_element(list);
1628
1629 list->start = curr;
1630
1631 for (size_t i = 0; i < cell.num_corners; ++i, curr = curr->next) {
1632
1633 curr->vec_coords[0] = cell.coordinates_xyz[i][0];
1634 curr->vec_coords[1] = cell.coordinates_xyz[i][1];
1635 curr->vec_coords[2] = cell.coordinates_xyz[i][2];
1637 cell.coordinates_xyz[i], cell.coordinates_xyz[(i+1)%(cell.num_corners)],
1638 cell.edge_type[i], cell_ordering,
1639 ((curr->edge_circle = circle_buffer + i)));
1640 curr->next =
1641 ((i + 1) == cell.num_corners)?
1642 (list->start):get_free_point_list_element(list);
1643 }
1644
1645 return remove_zero_length_edges(list);
1646}
1647
1648static struct point_list_element *
1650
1651 struct point_list_element * element;
1652
1653 if (list->free_elements == NULL) {
1654
1655 for (int i = 0; i < 7; ++i) {
1656
1657 element = (struct point_list_element *)xmalloc(1 * sizeof(*element));
1658
1659 element->next = list->free_elements;
1660 list->free_elements = element;
1661 }
1662
1663 element = (struct point_list_element *)xmalloc(1 * sizeof(*element));
1664
1665 } else {
1666
1667 element = list->free_elements;
1668 list->free_elements = list->free_elements->next;
1669 }
1670
1671 element->next = NULL;
1672 element->to_be_removed = 0;
1673
1674 return element;
1675}
1676
1677static void free_point_list(struct point_list * list) {
1678
1679 struct point_list_element * element;
1680
1681 if (list->start != NULL) {
1682
1683 struct point_list_element * temp = list->free_elements;
1684 list->free_elements = list->start->next;
1685 list->start->next = temp;
1686 }
1687
1688 while (list->free_elements != NULL) {
1689
1690 element = list->free_elements;
1691 list->free_elements = element->next;
1692 free(element);
1693 }
1694
1695 list->start = NULL;
1696 list->free_elements = NULL;
1697}
1698
1699static void compute_norm_vector(double a[], double b[], double norm[]) {
1700
1701 crossproduct_kahan(a, b, norm);
1702
1703 YAC_ASSERT(
1704 (fabs(norm[0]) >= tol) || (fabs(norm[1]) >= tol) || (fabs(norm[2]) >= tol),
1705 "ERROR: a and b are identical -> no norm vector")
1706
1707 normalise_vector(norm);
1708}
1709
1710static int is_empty_gc_cell(struct point_list * list, size_t num_edges) {
1711
1712 if (num_edges < 2) return 0;
1713
1714 // if the polygon only has two vertices
1715 if (num_edges == 2)
1716 // the cell is empty if both edges have the same type
1717 // (lon and gc are counted as the same type)
1718 return !((list->start->edge_circle->type == LAT_CIRCLE) ^
1719 (list->start->next->edge_circle->type == LAT_CIRCLE));
1720
1721 struct point_list_element * curr = list->start;
1722
1723 // if there is at least one lat circle edge, the cell is not empty
1724 for (size_t i = 0; i < num_edges; ++i, curr = curr->next)
1725 if (curr->edge_circle->type == LAT_CIRCLE) return 0;
1726
1727 // compute the norm vector for the first edge
1728 double norm_vec[3];
1729 compute_norm_vector(curr->vec_coords, curr->next->vec_coords, norm_vec);
1730 curr = curr->next->next;
1731
1732 // check whether at least one vertex of the polygon is on a diffent plane
1733 // than the one defined by the first edge of the polygon
1734 for (size_t i = 0; i < num_edges - 2; ++i, curr = curr->next)
1735 if (fabs(dotproduct(norm_vec, curr->vec_coords)) > yac_angle_tol)
1736 return 0;
1737
1738 return 1;
1739}
1740
1742 switch (type) {
1743 default:
1744 case (GREAT_CIRCLE): return YAC_GREAT_CIRCLE_EDGE;
1745 case (LON_CIRCLE): return YAC_LON_CIRCLE_EDGE;
1746 case (LAT_CIRCLE): return YAC_LAT_CIRCLE_EDGE;
1747 }
1748}
1749
1750static void generate_cell(struct point_list * list,
1751 struct yac_grid_cell * cell) {
1752
1753 if (list->start == NULL) {
1754
1755 reset_point_list(list);
1756 cell->num_corners = 0;
1757 return;
1758 }
1759
1760 size_t num_edges = remove_zero_length_edges(list);
1761
1762 if (is_empty_gc_cell(list, num_edges)) {
1763
1764 reset_point_list(list);
1765 cell->num_corners = 0;
1766 return;
1767 }
1768
1769 if (num_edges > cell->array_size) {
1770 free(cell->coordinates_xyz);
1771 free(cell->edge_type);
1772 cell->coordinates_xyz = (double(*)[3])xmalloc(num_edges * sizeof(*cell->coordinates_xyz));
1773 cell->edge_type = (enum yac_edge_type *)xmalloc(num_edges * sizeof(*cell->edge_type));
1774 cell->array_size = num_edges;
1775 }
1776 cell->num_corners = num_edges;
1777
1778 struct point_list_element * curr = list->start;
1779
1780 for (size_t i = 0; i < num_edges; ++i) {
1781
1782 cell->coordinates_xyz[i][0] = curr->vec_coords[0];
1783 cell->coordinates_xyz[i][1] = curr->vec_coords[1];
1784 cell->coordinates_xyz[i][2] = curr->vec_coords[2];
1785 cell->edge_type[i] = circle2edge_type(curr->edge_circle->type);
1786
1787 curr = curr->next;
1788 }
1789}
double yac_huiliers_area(struct yac_grid_cell cell)
Area calculation on a unit sphere taken from ESMF based on L'Huilier's Theorem.
Definition area.c:396
double yac_huiliers_area_info(struct yac_grid_cell cell, double *barycenter, double sign)
Definition area.c:482
Structs and interfaces for area calculations.
static int edge_direction(double *a, double *b)
int yac_point_in_cell(double point_coords[3], struct yac_grid_cell cell)
int yac_circle_point_is_inside(double const point[3], struct yac_circle *circle)
Definition clipping.c:461
static size_t overlap_cell_buffer_size
Definition clipping.c:67
static int is_empty_gc_cell(struct point_list *list, size_t num_edges)
Definition clipping.c:1710
void yac_correct_weights(size_t nSourceCells, double *weight)
correct interpolation weights
Definition clipping.c:1437
static void point_list_clipping(struct point_list *source_list, struct point_list target_list, size_t nct)
Definition clipping.c:1032
static void reset_point_list(struct point_list *list)
Definition clipping.c:1519
static struct yac_grid_cell * overlap_cell_buffer
Definition clipping.c:66
static void copy_point_list(struct point_list in, struct point_list *out)
Definition clipping.c:1055
static struct yac_grid_cell * get_overlap_cell_buffer(size_t N)
Definition clipping.c:69
static void circle_clipping(struct point_list *cell, size_t num_cell_edges, struct yac_circle **clipping_circles, size_t num_clipping_circles)
Definition clipping.c:576
void yac_cell_lat_clipping(size_t N, struct yac_grid_cell *cells, double lat_bounds[2], struct yac_grid_cell *overlap_buffer)
cell clipping to get the cells describing the intersections
Definition clipping.c:1317
void yac_cell_clipping(size_t N, struct yac_grid_cell *source_cell, struct yac_grid_cell target_cell, struct yac_grid_cell *overlap_buffer)
cell clipping to get the cells describing the intersections
Definition clipping.c:1080
static size_t remove_points(struct point_list *list)
Definition clipping.c:1531
static struct point_list_element * get_free_point_list_element(struct point_list *list)
Definition clipping.c:1649
void yac_compute_overlap_info(size_t N, struct yac_grid_cell *source_cell, struct yac_grid_cell target_cell, double *overlap_areas, double(*overlap_barycenters)[3])
calculates partial areas for all overlapping parts of the source cells with arbitrary target cells,...
Definition clipping.c:109
static size_t generate_point_list(struct point_list *list, struct yac_grid_cell cell, int cell_ordering, struct yac_circle *circle_buffer)
Definition clipping.c:1617
void yac_compute_overlap_areas(size_t N, struct yac_grid_cell *source_cell, struct yac_grid_cell target_cell, double *partial_areas)
calculates partial areas for all overlapping parts of the source cells with arbitrary target cells,...
Definition clipping.c:266
int yac_circle_compare(void const *a, void const *b)
Definition clipping.c:334
static double dotproduct(double a[], double b[])
Definition clipping.c:259
static void init_point_list(struct point_list *list)
Definition clipping.c:1513
static enum yac_edge_type circle2edge_type(enum yac_circle_type type)
Definition clipping.c:1741
static void free_point_list(struct point_list *list)
Definition clipping.c:1677
static int get_cell_points_ordering(struct yac_grid_cell cell)
Definition clipping.c:1458
static void yac_lon_lat_cell_lat_clipping(struct yac_grid_cell *cell, double z_upper_bound, double z_lower_bound, struct yac_grid_cell *overlap_buffer)
Definition clipping.c:1187
static enum yac_cell_type get_cell_type(struct yac_grid_cell target_cell)
Definition clipping.c:277
static struct yac_circle generate_lat_circle(double z, int north_is_out)
Definition clipping.c:451
static size_t remove_zero_length_edges(struct point_list *list)
returns number of edges/corners
Definition clipping.c:1587
static void compute_norm_vector(double a[], double b[], double norm[])
Definition clipping.c:1699
static void generate_cell(struct point_list *list, struct yac_grid_cell *cell)
Definition clipping.c:1750
int yac_circle_contains_north_pole(struct yac_circle *circle)
Definition clipping.c:553
void yac_circle_generate(double const *a, double const *b, enum yac_edge_type type, int edge_ordering, struct yac_circle *circle)
Definition clipping.c:405
static double get_closest_pole(struct point_list *cell_list)
Definition clipping.c:1296
static double const tol
Definition clipping.c:23
static double get_edge_direction(double *ref_corner, double *corner_a, double *corner_b)
Definition clipping.c:89
int yac_circle_compare_distances(double const a[3], double const b[3], struct yac_circle *circle)
Definition clipping.c:499
static void compute_norm_vector_kahan(double const a[3], double const b[3], double norm_vector[3])
Definition clipping.c:392
#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])
Definition geometry.h:415
static const struct sin_cos_angle SIN_COS_ZERO
Definition geometry.h:44
int yac_circle_intersect(struct yac_circle a, struct yac_circle b, double p[3], double q[3])
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_d(const double a[], const double b[], double cross[])
Definition geometry.h:347
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
static int sub_angles(struct sin_cos_angle a, struct sin_cos_angle b, struct sin_cos_angle *restrict sub)
Definition geometry.h:558
int yac_point_on_edge(double p[3], double const a[3], double const b[3], enum yac_circle_type circle_type)
#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
static double get_angle(double a_lon, double b_lon)
Definition geometry.h:127
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:364
static int sum_angles(struct sin_cos_angle a, struct sin_cos_angle b, struct sin_cos_angle *restrict sum)
Definition geometry.h:534
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
void yac_init_grid_cell(struct yac_grid_cell *cell)
Definition grid_cell.c:14
yac_cell_type
Definition grid_cell.h:26
@ YAC_LAT_CELL
Definition grid_cell.h:28
@ YAC_LON_LAT_CELL
Definition grid_cell.h:27
@ YAC_GREAT_CIRCLE_CELL
Definition grid_cell.h:29
@ YAC_MIXED_CELL
Definition grid_cell.h:30
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
enum callback_type type
#define xmalloc(size)
Definition ppm_xfuncs.h:66
struct yac_circle * edge_circle
Definition clipping.c:28
struct point_list_element * next
Definition clipping.c:30
double vec_coords[3]
Definition clipping.c:27
struct point_list_element * free_elements
Definition clipping.c:36
struct point_list_element * start
Definition clipping.c:35
double sin
Definition geometry.h:41
double cos
Definition geometry.h:41
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
int north_is_out
Definition geometry.h:89
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
size_t num_corners
Definition grid_cell.h:21
enum yac_edge_type * edge_type
Definition grid_cell.h:20
size_t array_size
Definition grid_cell.h:22
double(* coordinates_xyz)[3]
Definition grid_cell.h:19
#define MAX(a, b)
#define YAC_ASSERT_F(exp, format,...)
Definition yac_assert.h:19
#define YAC_ASSERT(exp, msg)
Definition yac_assert.h:16