1 /**
2 * @file clipping.c
3 *
4 * @copyright Copyright (C) 2013 Moritz Hanke <hanke@dkrz.de>
5 * Rene Redler <rene.redler@mpimet.mpg.de>
6 *
7 * @version 1.0
8 * @author Moritz Hanke <hanke@dkrz.de>
9 * Rene Redler <rene.redler@mpimet.mpg.de>
10 */
11 /*
12 * Keywords:
13 * Maintainer: Moritz Hanke <hanke@dkrz.de>
14 * Rene Redler <rene.redler@mpimet.mpg.de>
15 * URL: https://dkrz-sw.gitlab-pages.dkrz.de/yac/
16 *
17 * This file is part of YAC.
18 *
19 * Redistribution and use in source and binary forms, with or without
20 * modification, are permitted provided that the following conditions are
21 * met:
22 *
23 * Redistributions of source code must retain the above copyright notice,
24 * this list of conditions and the following disclaimer.
25 *
26 * Redistributions in binary form must reproduce the above copyright
27 * notice, this list of conditions and the following disclaimer in the
28 * documentation and/or other materials provided with the distribution.
29 *
30 * Neither the name of the DKRZ GmbH nor the names of its contributors
31 * may be used to endorse or promote products derived from this software
32 * without specific prior written permission.
33 *
34 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
35 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
36 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
37 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
38 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
39 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
40 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
41 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
42 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
43 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
44 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
45 */
46
47 #ifdef HAVE_CONFIG_H
48 // Get the definition of the 'restrict' keyword.
49 #include "config.h"
50 #endif
51
52 #include <stdio.h>
53 #include <stdlib.h>
54 #include <string.h>
55 #include <math.h>
56
57 #include "geometry.h"
58 #include "clipping.h"
59 #include "area.h"
60 #include "ensure_array_size.h"
61 #include "utils.h"
62
63 //#define YAC_VERBOSE_CLIPPING
64
65 static double const tol = 1.0e-12;
66
67 enum cell_type {
68 LON_LAT_CELL,
69 LAT_CELL,
70 GREAT_CIRCLE_CELL,
71 MIXED_CELL
72 };
73
74 enum intersect_results {
75 P_ON_A = 1,
76 Q_ON_A = 2,
77 P_ON_B = 4,
78 Q_ON_B = 8,
79 SAME_PLANE = 16,
80 };
81
82 struct point_list_element {
83
84 double vec_coords[3];
85 enum yac_edge_type edge_type; // type of edge with next corner
86 int to_be_removed;
87 struct point_list_element * next;
88 };
89
90 struct point_list {
91
92 struct point_list_element * start;
93 struct point_list_element * free_elements;
94 };
95
96 /* internal helper routines for working with linked lists of points */
97
98 static void init_point_list(struct point_list * list);
99
100 static void reset_point_list(struct point_list * list);
101
102 static size_t generate_point_list(
103 struct point_list * list, struct grid_cell cell);
104
105 static struct point_list_element *
106 get_free_point_list_element(struct point_list * list);
107
108 static size_t remove_points(struct point_list * list);
109 static size_t remove_zero_length_edges(struct point_list * list);
110
111 static void free_point_list(struct point_list * list);
112
113 static int get_cell_points_ordering(struct point_list * cell);
114
115 static void generate_overlap_cell(struct point_list * list,
116 struct grid_cell * cell);
117
118 static enum cell_type get_cell_type(struct grid_cell target_cell);
119
120 /* ------------------------- */
121
122 static struct grid_cell * overlap_cell_buffer = NULL;
123 static size_t overlap_cell_buffer_size = 0;
124
get_overlap_cell_buffer(size_t N)125 static inline struct grid_cell * get_overlap_cell_buffer(size_t N) {
126
127 // ensure that there are enough buffer cells
128
129 if (overlap_cell_buffer_size < N) {
130
131 size_t old_overlap_cell_buffer_size = overlap_cell_buffer_size;
132
133 ENSURE_ARRAY_SIZE(overlap_cell_buffer, overlap_cell_buffer_size, N);
134
135 for (; old_overlap_cell_buffer_size < overlap_cell_buffer_size;
136 ++old_overlap_cell_buffer_size)
137 yac_init_grid_cell(overlap_cell_buffer + old_overlap_cell_buffer_size);
138 }
139
140 return overlap_cell_buffer;
141 }
142
143 /* ------------------------- */
144
yac_compute_overlap_areas(size_t N,struct grid_cell * source_cell,struct grid_cell target_cell,double * partial_areas)145 void yac_compute_overlap_areas (size_t N,
146 struct grid_cell * source_cell,
147 struct grid_cell target_cell,
148 double * partial_areas) {
149
150 struct grid_cell * overlap_buffer = get_overlap_cell_buffer(N);
151
152 /* Do the clipping and get the cell for the overlapping area */
153
154 yac_cell_clipping ( N, source_cell, target_cell, overlap_buffer);
155
156 /* Get the partial areas for the overlapping regions */
157
158 for (size_t n = 0; n < N; n++) {
159 partial_areas[n] = yac_huiliers_area (overlap_buffer[n]);
160 // we cannot use pole_area because it is rather inaccurate for great circle
161 // edges that are nearly circles of longitude
162 //partial_areas[n] = pole_area (overlap_buffer[n]);
163 }
164
165 #ifdef YAC_VERBOSE_CLIPPING
166 for (size_t n = 0; n < N; n++)
167 printf("overlap area : %lf\n", partial_areas[n]);
168 #endif
169 }
170
171 /* ------------------------- */
172
compute_cell_barycenter(struct grid_cell cell,double * barycenter)173 static void compute_cell_barycenter(
174 struct grid_cell cell, double * barycenter) {
175
176 barycenter[0] = 0.0;
177 barycenter[1] = 0.0;
178 barycenter[2] = 0.0;
179
180 for (size_t k = 0; k < cell.num_corners; ++k) {
181 barycenter[0] += cell.coordinates_xyz[k][0];
182 barycenter[1] += cell.coordinates_xyz[k][1];
183 barycenter[2] += cell.coordinates_xyz[k][2];
184 }
185 normalise_vector(barycenter);
186 }
187
yac_compute_concave_overlap_info(size_t N,struct grid_cell * source_cell,struct grid_cell target_cell,double target_node_xyz[3],double * overlap_areas,double (* overlap_barycenters)[3])188 void yac_compute_concave_overlap_info (size_t N,
189 struct grid_cell * source_cell,
190 struct grid_cell target_cell,
191 double target_node_xyz[3],
192 double * overlap_areas,
193 double (*overlap_barycenters)[3]) {
194
195 struct grid_cell * overlap_buffer = get_overlap_cell_buffer(N);
196 enum cell_type target_cell_type;
197
198 if ( target_cell.num_corners > 3 )
199 target_cell_type = get_cell_type (target_cell);
200
201 if ( target_cell.num_corners < 4 || target_cell_type == LON_LAT_CELL ) {
202 yac_cell_clipping ( N, source_cell, target_cell, overlap_buffer);
203 for (size_t i = 0; i < N; ++i) {
204 compute_cell_barycenter(overlap_buffer[i], overlap_barycenters[i]);
205 overlap_areas[i] = yac_huiliers_area (overlap_buffer[i]);
206 }
207 return;
208 }
209
210 if ( target_node_xyz == NULL )
211 yac_internal_abort_message(
212 "ERROR(yac_compute_concave_overlap_info): "
213 "missing target cell centers (required for "
214 "cells with more than 3 corners; lon-lat cells are an exception",
215 __FILE__ , __LINE__);
216
217 struct grid_cell target_partial_cell =
218 {.coordinates_xyz = (double[3][3]){{-1}},
219 .edge_type = (enum yac_edge_type[3]) {GREAT_CIRCLE},
220 .num_corners = 3};
221
222 /* Do the clipping and get the cell for the overlapping area */
223
224 for ( size_t n = 0; n < N; n++) {
225 overlap_areas[n] = 0.0;
226 overlap_barycenters[n][0] = 0.0;
227 overlap_barycenters[n][1] = 0.0;
228 overlap_barycenters[n][2] = 0.0;
229 }
230
231 // common node point to all partial target cells
232 target_partial_cell.coordinates_xyz[0][0] = target_node_xyz[0];
233 target_partial_cell.coordinates_xyz[0][1] = target_node_xyz[1];
234 target_partial_cell.coordinates_xyz[0][2] = target_node_xyz[2];
235
236 for ( size_t num_corners = 0;
237 num_corners < target_cell.num_corners; ++num_corners ) {
238
239 size_t corner_a = num_corners;
240 size_t corner_b = (num_corners+1)%target_cell.num_corners;
241
242 // skip clipping and area calculation for degenerated triangles
243 //
244 // If this is not sufficient, instead we can try something like:
245 //
246 // struct point_list target_list
247 // init_point_list(&target_list);
248 // generate_point_list(&target_list, target_cell);
249 // struct grid_cell temp_target_cell;
250 // generate_overlap_cell(target_list, temp_target_cell);
251 // free_point_list(&target_list);
252 //
253 // and use temp_target_cell for triangulation.
254 //
255 // Compared to the if statement below the alternative seems
256 // to be quite costly.
257
258 if ( ( ( fabs(target_cell.coordinates_xyz[corner_a][0] -
259 target_cell.coordinates_xyz[corner_b][0]) < tol ) &&
260 ( fabs(target_cell.coordinates_xyz[corner_a][1] -
261 target_cell.coordinates_xyz[corner_b][1]) < tol ) &&
262 ( fabs(target_cell.coordinates_xyz[corner_a][2] -
263 target_cell.coordinates_xyz[corner_b][2]) < tol ) ) ||
264 ( ( fabs(target_cell.coordinates_xyz[corner_a][0] -
265 target_partial_cell.coordinates_xyz[0][0]) < tol ) &&
266 ( fabs(target_cell.coordinates_xyz[corner_a][1] -
267 target_partial_cell.coordinates_xyz[0][1]) < tol ) &&
268 ( fabs(target_cell.coordinates_xyz[corner_a][2] -
269 target_partial_cell.coordinates_xyz[0][2]) < tol ) ) ||
270 ( ( fabs(target_cell.coordinates_xyz[corner_b][0] -
271 target_partial_cell.coordinates_xyz[0][0]) < tol ) &&
272 ( fabs(target_cell.coordinates_xyz[corner_b][1] -
273 target_partial_cell.coordinates_xyz[0][1]) < tol ) &&
274 ( fabs(target_cell.coordinates_xyz[corner_b][2] -
275 target_partial_cell.coordinates_xyz[0][2]) < tol ) ) )
276 continue;
277
278 target_partial_cell.coordinates_xyz[1][0] = target_cell.coordinates_xyz[corner_a][0];
279 target_partial_cell.coordinates_xyz[1][1] = target_cell.coordinates_xyz[corner_a][1];
280 target_partial_cell.coordinates_xyz[1][2] = target_cell.coordinates_xyz[corner_a][2];
281 target_partial_cell.coordinates_xyz[2][0] = target_cell.coordinates_xyz[corner_b][0];
282 target_partial_cell.coordinates_xyz[2][1] = target_cell.coordinates_xyz[corner_b][1];
283 target_partial_cell.coordinates_xyz[2][2] = target_cell.coordinates_xyz[corner_b][2];
284
285 yac_cell_clipping ( N, source_cell, target_partial_cell, overlap_buffer);
286
287 /* Get the partial areas for the overlapping regions as sum over the partial target cells. */
288
289 for (size_t n = 0; n < N; n++) {
290 if (overlap_buffer[n].num_corners == 0) continue;
291 double overlap_area = yac_huiliers_area (overlap_buffer[n]);
292 double temp_cell_barycenter[3];
293 overlap_areas[n] += overlap_area;
294 compute_cell_barycenter(overlap_buffer[n], temp_cell_barycenter);
295 overlap_barycenters[n][0] += overlap_area * temp_cell_barycenter[0];
296 overlap_barycenters[n][1] += overlap_area * temp_cell_barycenter[1];
297 overlap_barycenters[n][2] += overlap_area * temp_cell_barycenter[2];
298 }
299 }
300 for (size_t n = 0; n < N; n++) normalise_vector(overlap_barycenters[n]);
301
302 #ifdef YAC_VERBOSE_CLIPPING
303 for (size_t n = 0; n < N; n++)
304 printf("overlap area %i: %lf \n", n, overlap_areas[n]);
305 #endif
306 }
307
dotproduct(double a[],double b[])308 static inline double dotproduct(double a[], double b[]) {
309
310 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
311 }
312
gc_cell_is_convex(struct grid_cell cell)313 static int gc_cell_is_convex(struct grid_cell cell) {
314
315 double const sq_angle_tol = yac_angle_tol * yac_angle_tol;
316 double (* restrict coordinates_xyz)[3] = cell.coordinates_xyz;
317
318 // for all edges (we skip the edge 0->num_corners-1; because we do not need it
319 // to determine whether a polygon is convex or not
320 for (size_t i = 0; i < cell.num_corners-1; ++i) {
321
322 // compute norm vector for current edge
323 double edge_norm_vec[3];
324 crossproduct_d(
325 coordinates_xyz[i], coordinates_xyz[i+1], edge_norm_vec);
326 // (sin(edge_length))^2
327 double sq_sin_edge_length = edge_norm_vec[0] * edge_norm_vec[0] +
328 edge_norm_vec[1] * edge_norm_vec[1] +
329 edge_norm_vec[2] * edge_norm_vec[2];
330 // if the edge is too short -> skip edge
331 if (sq_sin_edge_length < sq_angle_tol) continue;
332
333 normalise_vector(edge_norm_vec);
334
335 size_t same_side_count = 0;
336
337 // for the corners of the triangle
338 // (from corner index 0 to start of the edge)
339 for (size_t j = 0; j < i; ++j)
340 same_side_count += dotproduct(edge_norm_vec, coordinates_xyz[j]) > 0.0;
341
342 // for the remaining corners
343 for (size_t j = i + 2; j < cell.num_corners; ++j)
344 same_side_count += dotproduct(edge_norm_vec, coordinates_xyz[j]) > 0.0;
345
346 // if not all corners are on the same side of the current edge -> concave
347 if ((same_side_count != 0) &&
348 (same_side_count != cell.num_corners - 2)) return 0;
349 }
350
351 // the cell is convex
352 return 1;
353 }
354
355 /* ------------------------- */
356
yac_compute_concave_overlap_areas(size_t N,struct grid_cell * source_cell,struct grid_cell target_cell,double target_node_xyz[3],double * partial_areas)357 void yac_compute_concave_overlap_areas (size_t N,
358 struct grid_cell * source_cell,
359 struct grid_cell target_cell,
360 double target_node_xyz[3],
361 double * partial_areas) {
362 enum cell_type target_cell_type;
363
364 if ( target_cell.num_corners > 3 )
365 target_cell_type = get_cell_type (target_cell);
366
367 if (( target_cell.num_corners < 4 || target_cell_type == LON_LAT_CELL ) ||
368 ((target_cell_type == GREAT_CIRCLE_CELL) &&
369 gc_cell_is_convex(target_cell))) {
370 yac_compute_overlap_areas (N, source_cell, target_cell, partial_areas);
371 return;
372 }
373
374 if ( target_node_xyz == NULL )
375 yac_internal_abort_message("ERROR: missing target point coordinates "
376 "(x_coordinates == NULL || y_coordinates == NULL)",
377 __FILE__ , __LINE__);
378
379 struct grid_cell target_partial_cell =
380 {.coordinates_xyz = (double[3][3]){{-1}},
381 .edge_type = (enum yac_edge_type[3]) {GREAT_CIRCLE},
382 .num_corners = 3};
383
384 struct grid_cell * overlap_buffer = get_overlap_cell_buffer(N);
385
386 /* Do the clipping and get the cell for the overlapping area */
387
388 for ( size_t n = 0; n < N; n++) partial_areas[n] = 0.0;
389
390 // common node point to all partial target cells
391 target_partial_cell.coordinates_xyz[0][0] = target_node_xyz[0];
392 target_partial_cell.coordinates_xyz[0][1] = target_node_xyz[1];
393 target_partial_cell.coordinates_xyz[0][2] = target_node_xyz[2];
394
395 for ( size_t num_corners = 0;
396 num_corners < target_cell.num_corners; ++num_corners ) {
397
398 size_t corner_a = num_corners;
399 size_t corner_b = (num_corners+1)%target_cell.num_corners;
400
401 // skip clipping and area calculation for degenerated triangles
402 //
403 // If this is not sufficient, instead we can try something like:
404 //
405 // struct point_list target_list
406 // init_point_list(&target_list);
407 // generate_point_list(&target_list, target_cell);
408 // struct grid_cell temp_target_cell;
409 // generate_overlap_cell(target_list, temp_target_cell);
410 // free_point_list(&target_list);
411 //
412 // and use temp_target_cell for triangulation.
413 //
414 // Compared to the if statement below the alternative seems
415 // to be quite costly.
416
417 if ( ( ( fabs(target_cell.coordinates_xyz[corner_a][0] -
418 target_cell.coordinates_xyz[corner_b][0]) < tol ) &&
419 ( fabs(target_cell.coordinates_xyz[corner_a][1] -
420 target_cell.coordinates_xyz[corner_b][1]) < tol ) &&
421 ( fabs(target_cell.coordinates_xyz[corner_a][2] -
422 target_cell.coordinates_xyz[corner_b][2]) < tol ) ) ||
423 ( ( fabs(target_cell.coordinates_xyz[corner_a][0] -
424 target_partial_cell.coordinates_xyz[0][0]) < tol ) &&
425 ( fabs(target_cell.coordinates_xyz[corner_a][1] -
426 target_partial_cell.coordinates_xyz[0][1]) < tol ) &&
427 ( fabs(target_cell.coordinates_xyz[corner_a][2] -
428 target_partial_cell.coordinates_xyz[0][2]) < tol ) ) ||
429 ( ( fabs(target_cell.coordinates_xyz[corner_b][0] -
430 target_partial_cell.coordinates_xyz[0][0]) < tol ) &&
431 ( fabs(target_cell.coordinates_xyz[corner_b][1] -
432 target_partial_cell.coordinates_xyz[0][1]) < tol ) &&
433 ( fabs(target_cell.coordinates_xyz[corner_b][2] -
434 target_partial_cell.coordinates_xyz[0][2]) < tol ) ) )
435 continue;
436
437 target_partial_cell.coordinates_xyz[1][0] = target_cell.coordinates_xyz[corner_a][0];
438 target_partial_cell.coordinates_xyz[1][1] = target_cell.coordinates_xyz[corner_a][1];
439 target_partial_cell.coordinates_xyz[1][2] = target_cell.coordinates_xyz[corner_a][2];
440 target_partial_cell.coordinates_xyz[2][0] = target_cell.coordinates_xyz[corner_b][0];
441 target_partial_cell.coordinates_xyz[2][1] = target_cell.coordinates_xyz[corner_b][1];
442 target_partial_cell.coordinates_xyz[2][2] = target_cell.coordinates_xyz[corner_b][2];
443
444 yac_cell_clipping ( N, source_cell, target_partial_cell, overlap_buffer);
445
446 /* Get the partial areas for the overlapping regions as sum over the partial target cells. */
447
448 for (size_t n = 0; n < N; n++) {
449 partial_areas[n] += yac_huiliers_area (overlap_buffer[n]);
450 // we cannot use pole_area because it is rather inaccurate for great circle
451 // edges that are nearly circles of longitude
452 //partial_areas[n] = pole_area (overlap_buffer[n]);
453 }
454 }
455
456 #ifdef YAC_VERBOSE_CLIPPING
457 for (size_t n = 0; n < N; n++)
458 printf("overlap area %i: %lf \n", n, partial_areas[n]);
459 #endif
460 }
461
462 /* ------------------------- */
463
compute_norm_vector(double a[],double b[],double norm[])464 static void compute_norm_vector(double a[], double b[], double norm[]) {
465
466 crossproduct_ld(a, b, norm);
467
468 if ((fabs(norm[0]) < tol) &&
469 (fabs(norm[1]) < tol) &&
470 (fabs(norm[2]) < tol))
471 yac_internal_abort_message("ERROR: a and b are identical -> no norm vector\n",
472 __FILE__, __LINE__);
473
474 double scale = 1.0 / sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
475
476 norm[0] *= scale;
477 norm[1] *= scale;
478 norm[2] *= scale;
479 }
480
compute_lat_circle_z_value(double a[],double b[],double z[])481 static void compute_lat_circle_z_value(double a[], double b[], double z[]) {
482
483 double temp[3];
484
485 crossproduct_ld(a, b, temp);
486
487 z[0] = 0;
488 z[1] = 0;
489
490 if (temp[2] > 0)
491 z[2] = 1.0 - a[2];
492 else
493 z[2] = -1.0 - a[2];
494 }
495
496 /**
497 * Determines whether a given point is on a hemisphere that is defined by a plane
498 * through the middle of the sphere.\n
499 * The plane is defined by its norm vector.
500 * @param[in] point point to be checked
501 * @param[in] norm_vec norm vector of the plane dividing the sphere
502 * @returns 0 if the point is not inside the hemisphere
503 * 1 if the point is inside the hemisphere
504 * 2 if the point is in the plane
505 */
is_inside_gc(double point[],double norm_vec[])506 static int is_inside_gc(double point[], double norm_vec[]) {
507
508 // the product is defined as follows
509 // a * b = |a| * |b| * cos(alpha)
510 // where alpha is the angle between a and b
511
512 double dot = dotproduct(point, norm_vec);
513
514 // if the point is on the line
515 if (fabs(dot) <= yac_angle_tol * 1e-3)
516 return 2;
517
518 return dot < 0;
519 }
520
is_inside_latc(double point[],double z)521 static int is_inside_latc(double point[], double z) {
522
523 double temp = fabs(point[2] + z);
524
525 if (fabs(1.0 - temp) <= yac_angle_tol * 1e-3) return 2;
526 else return temp < 1.0;
527 }
528
is_inside_lat_circle(double z,double ref_z,int north_is_out)529 static int is_inside_lat_circle(double z, double ref_z, int north_is_out) {
530
531 double diff_z = z - ref_z;
532
533 if (fabs(diff_z) <= yac_angle_tol * 1e-3) return 2;
534 else return (diff_z > 0.0) ^ north_is_out;
535 }
536
is_inside(double point[],double help_vec[],enum yac_edge_type edge_type,int cell_points_ordering)537 static int is_inside(double point[], double help_vec[],
538 enum yac_edge_type edge_type,
539 int cell_points_ordering) {
540
541 int ret_val = 0;
542
543 switch (edge_type) {
544
545 case (LON_CIRCLE) :
546 case (GREAT_CIRCLE) :
547 ret_val = is_inside_gc(point, help_vec);
548 break;
549 case (LAT_CIRCLE) :
550 ret_val = is_inside_latc(point, help_vec[2]);
551 break;
552 default:
553 yac_internal_abort_message("invalid edge type\n", __FILE__, __LINE__);
554 };
555
556 if (ret_val == 2) return 2;
557 else return ret_val ^ cell_points_ordering;
558 }
559
get_cell_type(struct grid_cell target_cell)560 static enum cell_type get_cell_type(struct grid_cell target_cell) {
561
562 size_t count_lat_edges = 0, count_great_circle_edges = 0;
563
564 if ((target_cell.num_corners == 4) &&
565 ((target_cell.edge_type[0] == LAT_CIRCLE &&
566 target_cell.edge_type[1] == LON_CIRCLE &&
567 target_cell.edge_type[2] == LAT_CIRCLE &&
568 target_cell.edge_type[3] == LON_CIRCLE) ||
569 (target_cell.edge_type[0] == LON_CIRCLE &&
570 target_cell.edge_type[1] == LAT_CIRCLE &&
571 target_cell.edge_type[2] == LON_CIRCLE &&
572 target_cell.edge_type[3] == LAT_CIRCLE)))
573 return LON_LAT_CELL;
574 else
575 for (size_t i = 0; i < target_cell.num_corners; ++i)
576 if (target_cell.edge_type[i] == LON_CIRCLE ||
577 target_cell.edge_type[i] == GREAT_CIRCLE)
578 count_great_circle_edges++;
579 else
580 count_lat_edges++;
581
582 if (count_lat_edges && count_great_circle_edges)
583 return MIXED_CELL;
584 else if (count_lat_edges)
585 return LAT_CELL;
586 else
587 return GREAT_CIRCLE_CELL;
588 }
589
590 /**
591 * cell clipping using Sutherland-Hodgman algorithm;
592 */
point_list_clipping(struct point_list * source_list,int source_ordering,struct point_list target_list,int target_ordering,size_t nct,double * tgt_edge_norm_vec)593 static void point_list_clipping (struct point_list * source_list, int source_ordering,
594 struct point_list target_list, int target_ordering,
595 size_t nct, double * tgt_edge_norm_vec) {
596
597 struct {
598 double * edge_norm_vec;
599 struct point_list_element * point;
600 } tgt_points[nct];
601
602 // to avoid some problems that can occur close the the pole, we process the
603 // target lat-circle edges at the end
604 {
605
606 size_t count = 0;
607
608 struct point_list_element * curr_tgt_point = target_list.start;
609
610 for (size_t i = 0; i < nct; ++i, curr_tgt_point = curr_tgt_point->next) {
611
612 if (curr_tgt_point->edge_type == LAT_CIRCLE) continue;
613
614 tgt_points[count].edge_norm_vec = tgt_edge_norm_vec + i * 3;
615 tgt_points[count++].point = curr_tgt_point;
616 }
617
618 if (count != nct) {
619 for (size_t i = 0; i < nct; ++i, curr_tgt_point = curr_tgt_point->next) {
620
621 if (curr_tgt_point->edge_type != LAT_CIRCLE) continue;
622
623 tgt_points[count].edge_norm_vec = tgt_edge_norm_vec + i * 3;
624 tgt_points[count++].point = curr_tgt_point;
625 }
626 }
627 }
628
629 // count the number of edges in the source cell
630 size_t ncs;
631 if (source_list->start != NULL) {
632 ncs = 1;
633 for (struct point_list_element * curr = source_list->start->next,
634 * end = source_list->start;
635 curr != end; curr = curr->next, ++ncs);
636 } else {
637 ncs = 0;
638 }
639
640 // for all edges of the target cell
641 for (size_t i = 0; (i < nct) && (ncs > 1); ++i) {
642
643 struct point_list_element * prev_src_point = source_list->start;
644 struct point_list_element * curr_src_point = prev_src_point->next;
645
646 double * tgt_edge_coords[2] = {tgt_points[i].point->vec_coords,
647 tgt_points[i].point->next->vec_coords};
648 enum yac_edge_type tgt_edge_type = tgt_points[i].point->edge_type;
649 double * curr_tgt_edge_norm_vec = tgt_points[i].edge_norm_vec;
650
651 int prev_is_inside, curr_is_inside, last_curr_is_inside;
652
653 prev_is_inside = is_inside(prev_src_point->vec_coords,
654 tgt_points[i].edge_norm_vec,
655 tgt_points[i].point->edge_type, target_ordering);
656 last_curr_is_inside = prev_is_inside;
657
658 // for all edges of the source cell
659 for (size_t j = 0; j < ncs; ++j) {
660
661 if (j != ncs - 1)
662 curr_is_inside = is_inside(curr_src_point->vec_coords,
663 curr_tgt_edge_norm_vec, tgt_edge_type,
664 target_ordering);
665 else
666 curr_is_inside = last_curr_is_inside;
667
668 double * src_edge_coords[2] = {prev_src_point->vec_coords,
669 curr_src_point->vec_coords};
670 enum yac_edge_type src_edge_type = prev_src_point->edge_type;
671
672 double intersection[2][3];
673 int num_intersections = 0;
674
675 // One intersection is possible, if either the previous or the current
676 // source point is inside and the respective other one is outside.
677 // Two intersections are possible, if one edge is a circle of latitude,
678 // while this other is not. Additionally, this is not possible if one of
679 // the two source points is inside while the other is not or if both
680 // source points are on the plane of the target edge.
681 int one_intersect_expected = (curr_is_inside + prev_is_inside == 1);
682 int two_intersect_possible =
683 ((src_edge_type == LAT_CIRCLE) ^ (tgt_edge_type == LAT_CIRCLE)) &&
684 (curr_is_inside + prev_is_inside != 1) &&
685 (curr_is_inside + prev_is_inside != 4);
686
687 // if we need to check for intersections
688 if (one_intersect_expected || two_intersect_possible) {
689
690 // get intersection points
691 int intersect =
692 yac_intersect_vec(
693 src_edge_type, src_edge_coords[0], src_edge_coords[1],
694 tgt_edge_type, tgt_edge_coords[0], tgt_edge_coords[1],
695 intersection[0], intersection[1]);
696
697 // if yac_intersect_vec found something -> determine the number of
698 // intersections and to some processing on them if necessary.
699 if (intersect != -1) {
700
701 // if both edges are in the same plane
702 if (intersect & SAME_PLANE) {
703
704 prev_is_inside = 2;
705 curr_is_inside = 2;
706 num_intersections = 0;
707
708 one_intersect_expected = 0;
709
710 // if this is the first iteration
711 if (j == 0) last_curr_is_inside = 2;
712
713 // if there are two intersections on the edge
714 } else if ((intersect & P_ON_A) && (intersect & Q_ON_A)) {
715
716 // if only one was expected
717 if (one_intersect_expected) {
718 char error_string[1024];
719 sprintf(error_string, "ERROR: two intersections found, even "
720 "though no circle of latitude involed and "
721 "both edge vertices on different sides of "
722 "the cutting plane.\n"
723 "edge A %lf %lf %lf (edge type %d)\n"
724 "edge B %lf %lf %lf (edge type %d)\n"
725 "yac_intersect_vec return value %d\n"
726 "p %lf %lf %lf\n"
727 "q %lf %lf %lf\n",
728 src_edge_coords[0][0], src_edge_coords[0][1],
729 src_edge_coords[0][2], (int)src_edge_type,
730 tgt_edge_coords[0][0], tgt_edge_coords[0][1],
731 tgt_edge_coords[0][2], (int)tgt_edge_type, intersect,
732 intersection[0][0], intersection[0][1],
733 intersection[0][2], intersection[1][0],
734 intersection[1][1], intersection[1][2]);
735 yac_internal_abort_message(error_string, __FILE__, __LINE__);
736 }
737
738 // if both source edge vertices are on the same side of the
739 // target edge
740 if (curr_is_inside == prev_is_inside) {
741
742 // if the two intersection points are basically the same point
743 if (compare_angles(
744 get_vector_angle_2(intersection[0],
745 intersection[1]), SIN_COS_TOL) <= 0) {
746
747 num_intersections = 1;
748
749 } else {
750
751 // compute the angle between the previous source edge vertex and
752 // the two intersection points
753 struct sin_cos_angle angle_a =
754 get_vector_angle_2(src_edge_coords[0], intersection[0]);
755 struct sin_cos_angle angle_b =
756 get_vector_angle_2(src_edge_coords[0], intersection[1]);
757 int compare_angles_ret_value = compare_angles(angle_a, angle_b);
758
759 // if the two angles are identical, we assume that both points
760 // are more or less identical to each other
761 if (compare_angles_ret_value == 0) {
762
763 num_intersections = 1;
764
765 } else {
766
767 // if the first intersection point is further away from the
768 // previous source edge vertex than the second one ->
769 // switch them
770 if (compare_angles_ret_value > 0) {
771
772 double temp_intersection[3];
773 temp_intersection[0] = intersection[1][0];
774 temp_intersection[1] = intersection[1][1];
775 temp_intersection[2] = intersection[1][2];
776 intersection[1][0] = intersection[0][0];
777 intersection[1][1] = intersection[0][1];
778 intersection[1][2] = intersection[0][2];
779 intersection[0][0] = temp_intersection[0];
780 intersection[0][1] = temp_intersection[1];
781 intersection[0][2] = temp_intersection[2];
782 }
783
784 num_intersections = 2;
785 }
786 }
787
788 // one of the two source edge vertices is on the target edge
789 // => one of the two intersection points must be more or less
790 // identical to a vertex of the source edge
791 } else {
792
793 int src_on_edge_idx = (curr_is_inside == 2);
794
795 struct sin_cos_angle angles[2] =
796 {get_vector_angle_2(src_edge_coords[src_on_edge_idx],
797 intersection[0]),
798 get_vector_angle_2(src_edge_coords[src_on_edge_idx],
799 intersection[1])};
800 int compare_angles_ret_value[2] =
801 {compare_angles(angles[0], SIN_COS_TOL),
802 compare_angles(angles[1], SIN_COS_TOL)};
803
804 // if both intersection points are nearly identical to the src
805 // edge vertex
806 if ((compare_angles_ret_value[0] <= 0) &&
807 (compare_angles_ret_value[1] <= 0)) {
808
809 num_intersections = 0;
810
811 } else {
812
813 num_intersections = 1;
814
815 // if the second intersection points is closer than the first
816 if (compare_angles(angles[0], angles[1]) < 0) {
817 intersection[0][0] = intersection[1][0];
818 intersection[0][1] = intersection[1][1];
819 intersection[0][2] = intersection[1][2];
820 }
821 }
822 }
823
824 // if one of the intersection points is on the source edge
825 } else if (intersect & (P_ON_A | Q_ON_A)) {
826
827 // if one of the two source edge vertices is on the target edge
828 if ((curr_is_inside == 2) || (prev_is_inside == 2)) {
829
830 // we assume that the intersection point is the respective
831 // source edge vertex, which is on the target edge -> we do not
832 // need this intersection point
833 num_intersections = 0;
834
835 } else {
836
837 // if q is on the source edge
838 if (intersect & Q_ON_A) {
839 intersection[0][0] = intersection[1][0];
840 intersection[0][1] = intersection[1][1];
841 intersection[0][2] = intersection[1][2];
842 }
843
844 num_intersections = 1;
845 }
846 }
847 }
848
849 // If an intersection was expected but none was found, the cutting edge
850 // was most propably to close to an edge vertex. We can now assume that
851 // the respective vertex is directly on the cutting edge.
852 if (one_intersect_expected && (num_intersections == 0)) {
853
854 double temp_a, temp_b;
855
856 // which of the two source edge vertices closer to the cutting edge
857
858 switch (tgt_points[i].point->edge_type) {
859
860 case (LON_CIRCLE) :
861 case (GREAT_CIRCLE) :
862 temp_a = fabs(dotproduct(src_edge_coords[0],
863 curr_tgt_edge_norm_vec));
864 temp_b = fabs(dotproduct(src_edge_coords[1],
865 curr_tgt_edge_norm_vec));
866 break;
867 case (LAT_CIRCLE) :
868 temp_a = fabs(1.0 - fabs(src_edge_coords[0][2] +
869 curr_tgt_edge_norm_vec[2]));
870 temp_b = fabs(1.0 - fabs(src_edge_coords[1][2] +
871 curr_tgt_edge_norm_vec[2]));
872 break;
873 default:
874 yac_internal_abort_message("invalid edge type\n", __FILE__, __LINE__);
875 return;
876 };
877
878 if (temp_a < temp_b) prev_is_inside = 2;
879 else curr_is_inside = 2;
880
881 one_intersect_expected = 0;
882
883 // if this is the first iteration
884 if (j == 0) last_curr_is_inside = prev_is_inside;
885 }
886 }
887
888 // here we know the number of intersections and their location and we
889 // know the relation of the source edge vertices to the cutting edge
890 // (prev_is_inside and curr_is_inside)
891
892 // if the previous source edge vertex is outside -> dump it after cutting
893 // is finished
894 prev_src_point->to_be_removed = prev_is_inside == 0;
895
896 // the easiest case is that we expected one intersection and got one
897 if (one_intersect_expected) {
898
899 // insert an intersection point in the source point list in the
900 // current edge
901 struct point_list_element * intersect_point =
902 get_free_point_list_element(source_list);
903 prev_src_point->next = intersect_point;
904 intersect_point->next = curr_src_point;
905
906 intersect_point->vec_coords[0] = intersection[0][0];
907 intersect_point->vec_coords[1] = intersection[0][1];
908 intersect_point->vec_coords[2] = intersection[0][2];
909
910 intersect_point->edge_type =
911 (prev_is_inside)?tgt_edge_type:src_edge_type;
912
913 // if the cutting edge goes through both of the two source edge vertices
914 } else if ((prev_is_inside == 2) && (curr_is_inside == 2)) {
915
916 // if one of the two edges is a circle of latitude while the other is
917 // not
918 if ((src_edge_type == LAT_CIRCLE) ^ (tgt_edge_type == LAT_CIRCLE)) {
919
920 double cross_src_z = (double)((long double)src_edge_coords[0][0] *
921 (long double)src_edge_coords[1][1] -
922 (long double)src_edge_coords[0][1] *
923 (long double)src_edge_coords[1][0]);
924 double cross_tgt_z = (double)((long double)tgt_edge_coords[0][0] *
925 (long double)tgt_edge_coords[1][1] -
926 (long double)tgt_edge_coords[0][1] *
927 (long double)tgt_edge_coords[1][0]);
928
929 int same_ordering = source_ordering == target_ordering;
930 int same_direction = (cross_src_z > 0) == (cross_tgt_z > 0);
931
932 // if source and target cell have the same ordering and both
933 // edges have the same direction or if both cells have different
934 // ordering and the edges have different directions, then we might
935 // have to change the edge type of the source edge
936 if (same_ordering == same_direction) {
937
938 // well...it works...do not ask ;-)
939 // ((edge is on south hemisphere) XOR (direction of source edge) XOR
940 // (ordering of source cell))
941 if ((curr_src_point->vec_coords[2] > 0) ^
942 (cross_src_z < 0) ^ source_ordering)
943 prev_src_point->edge_type = LAT_CIRCLE;
944 else
945 prev_src_point->edge_type = GREAT_CIRCLE;
946 }
947 } else {
948 // nothing to be done here
949 }
950
951 // if we have no intersection, but the previous source point edge vector
952 // is on the target edge
953 } else if ((num_intersections == 0) && (prev_is_inside == 2)) {
954
955 // if the current source edge vertex is outside
956 if (curr_is_inside == 0) prev_src_point->edge_type = tgt_edge_type;
957
958 // if we have two intersections (only happens if one of the two edges is
959 // a circle of latitude while the other is not
960 } else if (num_intersections == 2) {
961
962 struct point_list_element * intersect_points[2] =
963 {get_free_point_list_element(source_list),
964 get_free_point_list_element(source_list)};
965
966 // add two points between the current source edge vertices
967 prev_src_point->next = intersect_points[0];
968 intersect_points[0]->next = intersect_points[1];
969 intersect_points[1]->next = curr_src_point;
970
971 intersect_points[0]->vec_coords[0] = intersection[0][0];
972 intersect_points[0]->vec_coords[1] = intersection[0][1];
973 intersect_points[0]->vec_coords[2] = intersection[0][2];
974 intersect_points[1]->vec_coords[0] = intersection[1][0];
975 intersect_points[1]->vec_coords[1] = intersection[1][1];
976 intersect_points[1]->vec_coords[2] = intersection[1][2];
977
978 // if a and b are outside
979 if ((prev_is_inside == 0) && (curr_is_inside == 0)) {
980 intersect_points[0]->edge_type = src_edge_type;
981 intersect_points[1]->edge_type = tgt_edge_type;
982
983 // if a and b are inside
984 } else if ((prev_is_inside == 1) && (curr_is_inside == 1)) {
985 intersect_points[0]->edge_type = tgt_edge_type;
986 intersect_points[1]->edge_type = src_edge_type;
987 } else {
988 yac_internal_abort_message(
989 "ERROR: one source edge vertex is on the target edge, therefore we "
990 "should not have two intersections.", __FILE__, __LINE__);
991 }
992
993 // if we have one intersection even though the two source edge vertices
994 // are not on opposite sides of the target edge
995 } else if (two_intersect_possible && (num_intersections == 1)) {
996
997 // if both points are outside -> circle of latitude and greate circle
998 // touch at a single point
999 if ((prev_is_inside == 0) && (curr_is_inside == 0)) {
1000
1001 // insert an intersection point in the source point list in the
1002 // current edge
1003 struct point_list_element * intersect_point =
1004 get_free_point_list_element(source_list);
1005 prev_src_point->next = intersect_point;
1006 intersect_point->next = curr_src_point;
1007
1008 intersect_point->vec_coords[0] = intersection[0][0];
1009 intersect_point->vec_coords[1] = intersection[0][1];
1010 intersect_point->vec_coords[2] = intersection[0][2];
1011
1012 intersect_point->edge_type = tgt_edge_type;
1013
1014 // if one of the two source edge vertices is on the edge while the other
1015 // is either inside or outside
1016 } else if ((prev_is_inside == 2) || (curr_is_inside == 2)) {
1017
1018
1019 // insert an intersection point in the source point list in the
1020 // current edge
1021 struct point_list_element * intersect_point =
1022 get_free_point_list_element(source_list);
1023 prev_src_point->next = intersect_point;
1024 intersect_point->next = curr_src_point;
1025
1026 intersect_point->vec_coords[0] = intersection[0][0];
1027 intersect_point->vec_coords[1] = intersection[0][1];
1028 intersect_point->vec_coords[2] = intersection[0][2];
1029
1030 // if the previous source edge vertex is on the target edge
1031 if (prev_is_inside == 2) {
1032 // if the current source edge vertex in on the outside
1033 if (curr_is_inside == 0) {
1034 prev_src_point->edge_type = src_edge_type;
1035 intersect_point->edge_type = tgt_edge_type;
1036 } else {
1037 prev_src_point->edge_type = tgt_edge_type;
1038 intersect_point->edge_type = src_edge_type;
1039 }
1040 // if the current source edge vertex is on the target edge
1041 } else {
1042 // if the previous source edge vertex in on the outside
1043 if (prev_is_inside == 0) {
1044 prev_src_point->edge_type = tgt_edge_type;
1045 intersect_point->edge_type = src_edge_type;
1046 } else {
1047 prev_src_point->edge_type = src_edge_type;
1048 intersect_point->edge_type = tgt_edge_type;
1049 }
1050 }
1051
1052 // if both source edge vertices are inside
1053 } else if ((curr_is_inside == 1) && (prev_is_inside == 1)) {
1054 // nothing to be done here
1055 } else {
1056 yac_internal_abort_message(
1057 "ERROR: unhandled intersection case", __FILE__, __LINE__);
1058 }
1059 }
1060
1061 prev_src_point = curr_src_point;
1062 curr_src_point = curr_src_point->next;
1063 prev_is_inside = curr_is_inside;
1064
1065 } // for all source edges
1066
1067 // remove all points that are to be deleted
1068 ncs = remove_points(source_list);
1069 }
1070 }
1071
1072 /**
1073 * cell clipping along a single latitude circle using
1074 * Sutherland-Hodgman algorithm;
1075 */
point_list_lat_clipping(struct point_list * cell,size_t num_corners,double lat_bound,int north_is_out)1076 static size_t point_list_lat_clipping (
1077 struct point_list * cell, size_t num_corners,
1078 double lat_bound, int north_is_out) {
1079
1080 if (cell->start == NULL) return 0;
1081
1082 // if the latitude circle actually is a pole
1083 if (fabs(lat_bound - 1.0) < yac_angle_tol) {
1084
1085 // if the outside pole and the latitude bound not are the same
1086 if ((lat_bound > 0.0) ^ north_is_out) {
1087 reset_point_list(cell);
1088 return 0;
1089 } else {
1090 return num_corners;
1091 }
1092 }
1093
1094 struct point_list_element * prev_point = cell->start;
1095 struct point_list_element * curr_point = prev_point->next;
1096
1097 int prev_is_inside, curr_is_inside, last_curr_is_inside;
1098
1099 prev_is_inside =
1100 is_inside_lat_circle(
1101 prev_point->vec_coords[2], lat_bound, north_is_out);
1102 last_curr_is_inside = prev_is_inside;
1103
1104 double lon_bound = sqrt(1.0 - lat_bound * lat_bound);
1105 double dummy_lat_edge_coords[2][3] =
1106 {{lon_bound, 0.0, lat_bound}, {0.0, lon_bound, lat_bound}};
1107
1108 // for all edges of the source cell
1109 for (size_t j = 0; j < num_corners; ++j) {
1110
1111 if (j != num_corners - 1)
1112 curr_is_inside =
1113 is_inside_lat_circle(
1114 curr_point->vec_coords[2], lat_bound, north_is_out);
1115 else
1116 curr_is_inside = last_curr_is_inside;
1117
1118 double * cell_edge_coords[2] =
1119 {prev_point->vec_coords, curr_point->vec_coords};
1120
1121 double intersection[2][3];
1122 int num_intersections = 0;
1123
1124 // One intersection is possible, if either the previous or the current
1125 // point is inside and the respective other one is outside.
1126 // Otherwise, two intersections are possible unless both cell points
1127 // are on the plane of the latitude circle.
1128 int one_intersect_expected = (curr_is_inside + prev_is_inside == 1);
1129 int two_intersect_possible =
1130 !one_intersect_expected && (curr_is_inside + prev_is_inside != 4);
1131
1132 // if we need to check for intersections
1133 if (one_intersect_expected || two_intersect_possible) {
1134
1135 // get intersection points
1136 int intersect =
1137 yac_gcxlatc_vec(
1138 cell_edge_coords[0], cell_edge_coords[1],
1139 dummy_lat_edge_coords[0], dummy_lat_edge_coords[1],
1140 intersection[0], intersection[1]);
1141
1142 // if yac_intersect_lat_circle_vec found something -> determine
1143 // the number of intersections and to some processing on them
1144 // if necessary.
1145 if (intersect != -1) {
1146
1147 // if both edges are in the same plane
1148 if (intersect & SAME_PLANE) {
1149
1150 prev_is_inside = 2;
1151 curr_is_inside = 2;
1152 num_intersections = 0;
1153
1154 one_intersect_expected = 0;
1155
1156 // if this is the first iteration
1157 if (j == 0) last_curr_is_inside = 2;
1158
1159 // if there are two intersections on the edge
1160 } else if ((intersect & P_ON_A) && (intersect & Q_ON_A)) {
1161
1162 // if only one was expected
1163 if (one_intersect_expected) {
1164 char error_string[1024];
1165 sprintf(error_string, "ERROR: two intersections found, even "
1166 "though both edge vertices on different "
1167 "sides of the cutting plane.\n"
1168 "edge (%lf %lf %lf; %lf %lf %lf)\n"
1169 "lat %lf\n"
1170 "yac_gcxlatc_vec return value %d\n"
1171 "p %lf %lf %lf\n"
1172 "q %lf %lf %lf\n",
1173 cell_edge_coords[0][0], cell_edge_coords[0][1],
1174 cell_edge_coords[0][2], cell_edge_coords[1][0],
1175 cell_edge_coords[1][1], cell_edge_coords[1][2],
1176 lat_bound, intersect,
1177 intersection[0][0], intersection[0][1], intersection[0][2],
1178 intersection[1][0], intersection[1][1], intersection[1][2]);
1179 yac_internal_abort_message(error_string, __FILE__, __LINE__);
1180 }
1181
1182 // if both cell edge vertices are on the same side of the
1183 // target edge
1184 if (curr_is_inside == prev_is_inside) {
1185
1186 // if the two intersection points are basically the same point
1187 if (compare_angles(
1188 get_vector_angle_2(intersection[0],
1189 intersection[1]), SIN_COS_TOL) <= 0) {
1190
1191 num_intersections = 1;
1192
1193 } else {
1194
1195 // compute the angle between the previous cell edge vertex and
1196 // the two intersection points
1197 struct sin_cos_angle angle_a =
1198 get_vector_angle_2(cell_edge_coords[0], intersection[0]);
1199 struct sin_cos_angle angle_b =
1200 get_vector_angle_2(cell_edge_coords[0], intersection[1]);
1201 int compare_angles_ret_value = compare_angles(angle_a, angle_b);
1202
1203 // if the two angles are identical, we assume that both points
1204 // are more or less identical to each other
1205 if (compare_angles_ret_value == 0) {
1206
1207 num_intersections = 1;
1208
1209 } else {
1210
1211 // if the first intersection point is further away from the
1212 // previous cell edge vertex than the second one ->
1213 // switch them
1214 if (compare_angles_ret_value > 0) {
1215
1216 double temp_intersection[3];
1217 temp_intersection[0] = intersection[1][0];
1218 temp_intersection[1] = intersection[1][1];
1219 temp_intersection[2] = intersection[1][2];
1220 intersection[1][0] = intersection[0][0];
1221 intersection[1][1] = intersection[0][1];
1222 intersection[1][2] = intersection[0][2];
1223 intersection[0][0] = temp_intersection[0];
1224 intersection[0][1] = temp_intersection[1];
1225 intersection[0][2] = temp_intersection[2];
1226 }
1227
1228 num_intersections = 2;
1229 }
1230 }
1231
1232 // one of the two cell edge vertices is on the target edge
1233 // => one of the two intersection points must be more or less
1234 // identical to a vertex of the cell edge
1235 } else {
1236
1237 int on_edge_idx = (curr_is_inside == 2);
1238
1239 struct sin_cos_angle angles[2] =
1240 {get_vector_angle_2(cell_edge_coords[on_edge_idx],
1241 intersection[0]),
1242 get_vector_angle_2(cell_edge_coords[on_edge_idx],
1243 intersection[1])};
1244 int compare_angles_ret_value[2] =
1245 {compare_angles(angles[0], SIN_COS_TOL),
1246 compare_angles(angles[1], SIN_COS_TOL)};
1247
1248 // if both intersection points are nearly identical to the cell
1249 // edge vertex
1250 if ((compare_angles_ret_value[0] <= 0) &&
1251 (compare_angles_ret_value[1] <= 0)) {
1252
1253 num_intersections = 0;
1254
1255 } else {
1256
1257 num_intersections = 1;
1258
1259 // if the second intersection points is closer than the first
1260 if (compare_angles(angles[0], angles[1]) < 0) {
1261 intersection[0][0] = intersection[1][0];
1262 intersection[0][1] = intersection[1][1];
1263 intersection[0][2] = intersection[1][2];
1264 }
1265 }
1266 }
1267
1268 // if one of the intersection points is on the cell edge
1269 } else if (intersect & (P_ON_A | Q_ON_A)) {
1270
1271 // if one of the two cell edge vertices is on the target edge
1272 if ((curr_is_inside == 2) || (prev_is_inside == 2)) {
1273
1274 // we assume that the intersection point is the respective
1275 // cell edge vertex, which is on the target edge -> we do not
1276 // need this intersection point
1277 num_intersections = 0;
1278
1279 } else {
1280
1281 // if q is on the cell edge
1282 if (intersect & Q_ON_A) {
1283 intersection[0][0] = intersection[1][0];
1284 intersection[0][1] = intersection[1][1];
1285 intersection[0][2] = intersection[1][2];
1286 }
1287
1288 num_intersections = 1;
1289 }
1290 }
1291 }
1292
1293 // If an intersection was expected but none was found, the cutting edge
1294 // was most propably to close to an edge vertex. We can now assume that
1295 // the respective vertex is directly on the cutting edge.
1296 if (one_intersect_expected && (num_intersections == 0)) {
1297
1298 double temp_a, temp_b;
1299
1300 // which of the two cell edge vertices closer to the cutting edge
1301 temp_a = fabs(1.0 - fabs(cell_edge_coords[0][2] + lat_bound));
1302 temp_b = fabs(1.0 - fabs(cell_edge_coords[1][2] + lat_bound));
1303
1304 if (temp_a < temp_b) prev_is_inside = 2;
1305 else curr_is_inside = 2;
1306
1307 one_intersect_expected = 0;
1308
1309 // if this is the first iteration
1310 if (j == 0) last_curr_is_inside = prev_is_inside;
1311 }
1312 }
1313
1314 // here we know the number of intersections and their location and we
1315 // know the relation of the cell edge vertices to the cutting edge
1316 // (prev_is_inside and curr_is_inside)
1317
1318 // if the previous cell edge vertex is outside -> dump it after cutting
1319 // is finished
1320 prev_point->to_be_removed = prev_is_inside == 0;
1321
1322 // the easiest case is that we expected one intersection and got one
1323 if (one_intersect_expected) {
1324
1325 // insert an intersection point in the cell point list in the
1326 // current edge
1327 struct point_list_element * intersect_point =
1328 get_free_point_list_element(cell);
1329 prev_point->next = intersect_point;
1330 intersect_point->next = curr_point;
1331
1332 intersect_point->vec_coords[0] = intersection[0][0];
1333 intersect_point->vec_coords[1] = intersection[0][1];
1334 intersect_point->vec_coords[2] = intersection[0][2];
1335
1336 intersect_point->edge_type =
1337 (prev_is_inside)?LAT_CIRCLE:GREAT_CIRCLE;
1338
1339 // if the cutting edge goes through both of the two cell edge vertices
1340 } else if ((prev_is_inside == 2) && (curr_is_inside == 2)) {
1341
1342 prev_point->edge_type = LAT_CIRCLE;
1343
1344 // if we have no intersection, but the previous cell point edge vector
1345 // is on the target edge
1346 } else if ((num_intersections == 0) && (prev_is_inside == 2)) {
1347
1348 // if the current cell edge vertex is outside
1349 if (curr_is_inside == 0) prev_point->edge_type = LAT_CIRCLE;
1350
1351 // if we have two intersections
1352 } else if (num_intersections == 2) {
1353
1354 struct point_list_element * intersect_points[2] =
1355 {get_free_point_list_element(cell),
1356 get_free_point_list_element(cell)};
1357
1358 // add two points between the current cell edge vertices
1359 prev_point->next = intersect_points[0];
1360 intersect_points[0]->next = intersect_points[1];
1361 intersect_points[1]->next = curr_point;
1362
1363 intersect_points[0]->vec_coords[0] = intersection[0][0];
1364 intersect_points[0]->vec_coords[1] = intersection[0][1];
1365 intersect_points[0]->vec_coords[2] = intersection[0][2];
1366 intersect_points[1]->vec_coords[0] = intersection[1][0];
1367 intersect_points[1]->vec_coords[1] = intersection[1][1];
1368 intersect_points[1]->vec_coords[2] = intersection[1][2];
1369
1370 // if a and b are outside
1371 if ((prev_is_inside == 0) && (curr_is_inside == 0)) {
1372 intersect_points[0]->edge_type = GREAT_CIRCLE;
1373 intersect_points[1]->edge_type = LAT_CIRCLE;
1374
1375 // if a and b are inside
1376 } else if ((prev_is_inside == 1) && (curr_is_inside == 1)) {
1377 intersect_points[0]->edge_type = LAT_CIRCLE;
1378 intersect_points[1]->edge_type = GREAT_CIRCLE;
1379 } else {
1380 yac_internal_abort_message(
1381 "ERROR: one cell edge vertex is on the latitude circle, therefore "
1382 "we should not have two intersections.", __FILE__, __LINE__);
1383 }
1384
1385 // if we have one intersection even though the two cell edge vertices
1386 // are not on opposite sides of the latitude circle
1387 } else if (two_intersect_possible && (num_intersections == 1)) {
1388
1389 // if both points are outside -> circle of latitude and greate circle
1390 // touch at a single point
1391 if ((prev_is_inside == 0) && (curr_is_inside == 0)) {
1392
1393 // insert an intersection point in the cell point list in the
1394 // current edge
1395 struct point_list_element * intersect_point =
1396 get_free_point_list_element(cell);
1397 prev_point->next = intersect_point;
1398 intersect_point->next = curr_point;
1399
1400 intersect_point->vec_coords[0] = intersection[0][0];
1401 intersect_point->vec_coords[1] = intersection[0][1];
1402 intersect_point->vec_coords[2] = intersection[0][2];
1403
1404 intersect_point->edge_type = LAT_CIRCLE;
1405
1406 // if one of the two cell edge vertices is on the edge while the other
1407 // is either inside or outside
1408 } else if ((prev_is_inside == 2) || (curr_is_inside == 2)) {
1409
1410
1411 // insert an intersection point in the cell point list in the
1412 // current edge
1413 struct point_list_element * intersect_point =
1414 get_free_point_list_element(cell);
1415 prev_point->next = intersect_point;
1416 intersect_point->next = curr_point;
1417
1418 intersect_point->vec_coords[0] = intersection[0][0];
1419 intersect_point->vec_coords[1] = intersection[0][1];
1420 intersect_point->vec_coords[2] = intersection[0][2];
1421
1422 // if the previous cell edge vertex is on the latitude circle
1423 if (prev_is_inside == 2) {
1424 // if the current cell edge vertex in on the outside
1425 if (curr_is_inside == 0) {
1426 prev_point->edge_type = GREAT_CIRCLE;
1427 intersect_point->edge_type = LAT_CIRCLE;
1428 } else {
1429 prev_point->edge_type = LAT_CIRCLE;
1430 intersect_point->edge_type = GREAT_CIRCLE;
1431 }
1432 // if the current cell edge vertex is on the latitude circle
1433 } else {
1434 // if the previous cell edge vertex in on the outside
1435 if (prev_is_inside == 0) {
1436 prev_point->edge_type = LAT_CIRCLE;
1437 intersect_point->edge_type = GREAT_CIRCLE;
1438 } else {
1439 prev_point->edge_type = GREAT_CIRCLE;
1440 intersect_point->edge_type = LAT_CIRCLE;
1441 }
1442 }
1443
1444 // if both source edge vertices are inside
1445 } else if ((curr_is_inside == 1) && (prev_is_inside == 1)) {
1446 // nothing to be done here
1447 } else {
1448 yac_internal_abort_message(
1449 "ERROR: unhandled intersection case", __FILE__, __LINE__);
1450 }
1451 }
1452
1453 prev_point = curr_point;
1454 curr_point = curr_point->next;
1455 prev_is_inside = curr_is_inside;
1456
1457 } // for all edges of the cell
1458
1459 // remove all points that are to be deleted
1460 return remove_points(cell);
1461 }
1462
copy_point_list(struct point_list in,struct point_list * out)1463 static void copy_point_list(struct point_list in, struct point_list * out) {
1464
1465 reset_point_list(out);
1466
1467 struct point_list_element * curr = in.start;
1468
1469 if (curr == NULL) return;
1470
1471 struct point_list_element * new_point_list = get_free_point_list_element(out);
1472 out->start = new_point_list;
1473 *new_point_list = *curr;
1474 curr = curr->next;
1475
1476 do {
1477
1478 new_point_list->next = get_free_point_list_element(out);
1479 new_point_list = new_point_list->next;
1480 *new_point_list = *curr;
1481 curr = curr->next;
1482
1483 } while (curr != in.start);
1484
1485 new_point_list->next = out->start;
1486 }
1487
yac_cell_clipping(size_t N,struct grid_cell * source_cell,struct grid_cell target_cell,struct grid_cell * overlap_buffer)1488 void yac_cell_clipping (size_t N,
1489 struct grid_cell * source_cell,
1490 struct grid_cell target_cell,
1491 struct grid_cell * overlap_buffer) {
1492
1493 size_t ncs; /* number of vertices of source cell */
1494 size_t nct; /* number of vertices of target cell */
1495
1496 struct point_list target_list, source_list, temp_list;
1497
1498 int target_ordering; /* ordering of target cell corners */
1499 int source_ordering; /* ordering of source cell corners */
1500
1501
1502 enum cell_type tgt_cell_type = get_cell_type(target_cell);
1503
1504 if (tgt_cell_type == MIXED_CELL)
1505 yac_internal_abort_message("invalid target cell type (cell contains edges consisting "
1506 "of great circles and circles of latitude)\n",
1507 __FILE__, __LINE__);
1508
1509 init_point_list(&temp_list);
1510
1511 // generate point list for target cell (clip cell)
1512 init_point_list(&target_list);
1513 nct = generate_point_list(&target_list, target_cell);
1514
1515 // if there is no target cell (e.g. if all edges of target cell have a length
1516 // of zero)
1517 if (nct == 0) {
1518 free_point_list(&target_list);
1519 for (size_t i = 0; i < N; ++i) overlap_buffer[i].num_corners = 0;
1520 return;
1521 }
1522
1523 // compute target direction
1524 target_ordering = get_cell_points_ordering(&target_list);
1525
1526 // if all corners of the target cell are on the same great circle
1527 if (target_ordering == -1) {
1528 free_point_list(&target_list);
1529 for (size_t n = 0; n < N; n++ ) overlap_buffer[n].num_corners = 0;
1530 return;
1531 }
1532
1533 struct point_list_element * prev_tgt_point = target_list.start;
1534 struct point_list_element * curr_tgt_point = target_list.start->next;
1535
1536 /* norm vector for temporary target edge plane */
1537 double * norm_vec = (double *)xmalloc(3 * nct * sizeof(*norm_vec));
1538
1539 // compute norm vectors for all edges
1540 // or for lat circle edges a special z value
1541 for (size_t i = 0; i < nct; ++i) {
1542
1543 switch (prev_tgt_point->edge_type) {
1544
1545 case (LON_CIRCLE) :
1546 case (GREAT_CIRCLE) :
1547 compute_norm_vector(prev_tgt_point->vec_coords, curr_tgt_point->vec_coords,
1548 norm_vec + 3 * i);
1549 break;
1550 case (LAT_CIRCLE):
1551 compute_lat_circle_z_value(prev_tgt_point->vec_coords, curr_tgt_point->vec_coords,
1552 norm_vec + 3 * i);
1553 break;
1554 default:
1555 yac_internal_abort_message("invalid edge type\n", __FILE__, __LINE__);
1556 };
1557 prev_tgt_point = curr_tgt_point;
1558 curr_tgt_point = curr_tgt_point->next;
1559 }
1560
1561 init_point_list(&source_list);
1562
1563 // for all source cells
1564 for (size_t n = 0; n < N; n++ ) {
1565
1566 overlap_buffer[n].num_corners = 0;
1567
1568 enum cell_type src_cell_type = get_cell_type(source_cell[n]);
1569
1570 if (src_cell_type == MIXED_CELL)
1571 yac_internal_abort_message("invalid source cell type (cell contains edges consisting "
1572 "of great circles and circles of latitude)\n",
1573 __FILE__, __LINE__);
1574
1575 if (source_cell[n].num_corners < 2)
1576 continue;
1577
1578 // generate point list for current source list
1579 ncs = generate_point_list(&source_list, source_cell[n]);
1580
1581 // compute source direction
1582 source_ordering = get_cell_points_ordering(&source_list);
1583
1584 // if all corners of the source cell are on the same great circle
1585 if (source_ordering == -1) continue;
1586
1587 struct point_list * overlap;
1588 double fabs_tgt_coordinate_z = fabs(target_cell.coordinates_xyz[0][2]);
1589 double fabs_src_coordinate_z = fabs(source_cell[n].coordinates_xyz[0][2]);
1590
1591 // in this case the source an target cell are both LAT_CELL's than the
1592 // bigger one has to be the target cell
1593 // a similar problem occurs when the target cell is a LAT_CELL and the
1594 // source is a GREAT_CIRCLE_CELL which overlaps with the pole that is also
1595 // include in the target cell
1596 if (((tgt_cell_type == LAT_CELL) && (src_cell_type == GREAT_CIRCLE_CELL)) ||
1597 ((tgt_cell_type == LAT_CELL) && (src_cell_type == LAT_CELL) &&
1598 (fabs_tgt_coordinate_z > fabs_src_coordinate_z))) {
1599
1600 copy_point_list(target_list, &temp_list);
1601
1602 double temp_norm_vec[3*ncs];
1603 struct point_list_element * src_point = source_list.start;
1604
1605
1606 for (size_t i = 0; i < ncs; ++i) {
1607 switch (src_point->edge_type) {
1608
1609 case (LON_CIRCLE) :
1610 case (GREAT_CIRCLE) :
1611 compute_norm_vector(src_point->vec_coords,
1612 src_point->next->vec_coords,
1613 temp_norm_vec + 3 * i);
1614 break;
1615 case (LAT_CIRCLE):
1616 compute_lat_circle_z_value(src_point->vec_coords,
1617 src_point->next->vec_coords,
1618 temp_norm_vec + 3 * i);
1619 break;
1620 default:
1621 yac_internal_abort_message("invalid edge type\n", __FILE__, __LINE__);
1622 };
1623 src_point = src_point->next;
1624 }
1625
1626 point_list_clipping(&temp_list, target_ordering,
1627 source_list, source_ordering, ncs, temp_norm_vec);
1628
1629 overlap = &temp_list;
1630
1631 } else {
1632
1633 point_list_clipping(&source_list, source_ordering,
1634 target_list, target_ordering, nct, norm_vec);
1635
1636 overlap = &source_list;
1637 }
1638
1639 if (overlap->start != NULL)
1640 generate_overlap_cell(overlap, overlap_buffer + n);
1641 }
1642
1643 free(norm_vec);
1644 free_point_list(&source_list);
1645 free_point_list(&target_list);
1646 free_point_list(&temp_list);
1647 }
1648
yac_lon_lat_cell_lat_clipping(struct grid_cell * cell,double z_upper_bound,double z_lower_bound,struct grid_cell * overlap_buffer)1649 static void yac_lon_lat_cell_lat_clipping(
1650 struct grid_cell * cell, double z_upper_bound, double z_lower_bound,
1651 struct grid_cell * overlap_buffer) {
1652
1653 double cell_upper_bound = cell->coordinates_xyz[0][2];
1654 double cell_lower_bound = cell->coordinates_xyz[2][2];
1655
1656 int upper_idx[2], lower_idx[2];
1657
1658 if (cell_upper_bound < cell_lower_bound) {
1659 double temp = cell_upper_bound;
1660 cell_upper_bound = cell_lower_bound;
1661 cell_lower_bound = temp;
1662 if (cell->edge_type[0] == LAT_CIRCLE) {
1663 upper_idx[0] = 2;
1664 upper_idx[1] = 3;
1665 lower_idx[0] = 1;
1666 lower_idx[1] = 0;
1667 } else {
1668 upper_idx[0] = 2;
1669 upper_idx[1] = 1;
1670 lower_idx[0] = 3;
1671 lower_idx[1] = 0;
1672 }
1673 } else {
1674 if (cell->edge_type[0] == LAT_CIRCLE) {
1675 upper_idx[0] = 0;
1676 upper_idx[1] = 1;
1677 lower_idx[0] = 3;
1678 lower_idx[1] = 2;
1679 } else {
1680 upper_idx[0] = 0;
1681 upper_idx[1] = 3;
1682 lower_idx[0] = 1;
1683 lower_idx[1] = 2;
1684 }
1685 }
1686
1687 // if z_upper_bound and z_lower_bound are identical or
1688 // if cell does not overlap with latitude band
1689 if ((z_upper_bound == z_lower_bound) ||
1690 (cell_upper_bound <= z_lower_bound) ||
1691 (cell_lower_bound >= z_upper_bound)) {
1692 overlap_buffer->num_corners = 0;
1693 return;
1694 }
1695
1696 if (overlap_buffer->array_size < 4) {
1697 overlap_buffer->coordinates_xyz =
1698 xmalloc(4 * sizeof(*(overlap_buffer->coordinates_xyz)));
1699 overlap_buffer->edge_type =
1700 xmalloc(4 * sizeof(*(overlap_buffer->edge_type)));
1701 overlap_buffer->array_size = 4;
1702 }
1703
1704 memcpy(overlap_buffer->coordinates_xyz, cell->coordinates_xyz,
1705 4 * sizeof(*(cell->coordinates_xyz)));
1706 memcpy(overlap_buffer->edge_type, cell->edge_type,
1707 4 * sizeof(*(cell->edge_type)));
1708 overlap_buffer->num_corners = 4;
1709
1710
1711 double tmp_scale;
1712 double * p[2];
1713 if (fabs(cell_lower_bound) < fabs(cell_upper_bound)) {
1714 p[0] = cell->coordinates_xyz[lower_idx[0]];
1715 p[1] = cell->coordinates_xyz[lower_idx[1]];
1716 tmp_scale = cell->coordinates_xyz[lower_idx[0]][0] *
1717 cell->coordinates_xyz[lower_idx[0]][0] +
1718 cell->coordinates_xyz[lower_idx[0]][1] *
1719 cell->coordinates_xyz[lower_idx[0]][1];
1720 } else {
1721 p[0] = cell->coordinates_xyz[upper_idx[0]];
1722 p[1] = cell->coordinates_xyz[upper_idx[1]];
1723 tmp_scale = cell->coordinates_xyz[upper_idx[0]][0] *
1724 cell->coordinates_xyz[upper_idx[0]][0] +
1725 cell->coordinates_xyz[upper_idx[0]][1] *
1726 cell->coordinates_xyz[upper_idx[0]][1];
1727 }
1728
1729 // if upper bound overlaps with cell
1730 if ((z_upper_bound < cell_upper_bound) &&
1731 (z_upper_bound > cell_lower_bound)) {
1732
1733 double scale = sqrt((1.0 - z_upper_bound * z_upper_bound) / tmp_scale);
1734
1735 overlap_buffer->coordinates_xyz[upper_idx[0]][0] = p[0][0] * scale;
1736 overlap_buffer->coordinates_xyz[upper_idx[0]][1] = p[0][1] * scale;
1737 overlap_buffer->coordinates_xyz[upper_idx[0]][2] = z_upper_bound;
1738 overlap_buffer->coordinates_xyz[upper_idx[1]][0] = p[1][0] * scale;
1739 overlap_buffer->coordinates_xyz[upper_idx[1]][1] = p[1][1] * scale;
1740 overlap_buffer->coordinates_xyz[upper_idx[1]][2] = z_upper_bound;
1741 }
1742
1743 // if lower bound overlaps with cell
1744 if ((z_lower_bound < cell_upper_bound) &&
1745 (z_lower_bound > cell_lower_bound)) {
1746
1747 double scale = sqrt((1.0 - z_lower_bound * z_lower_bound) / tmp_scale);
1748
1749 overlap_buffer->coordinates_xyz[lower_idx[0]][0] = p[0][0] * scale;
1750 overlap_buffer->coordinates_xyz[lower_idx[0]][1] = p[0][1] * scale;
1751 overlap_buffer->coordinates_xyz[lower_idx[0]][2] = z_lower_bound;
1752 overlap_buffer->coordinates_xyz[lower_idx[1]][0] = p[1][0] * scale;
1753 overlap_buffer->coordinates_xyz[lower_idx[1]][1] = p[1][1] * scale;
1754 overlap_buffer->coordinates_xyz[lower_idx[1]][2] = z_lower_bound;
1755 }
1756 }
1757
get_closest_pole(struct point_list * cell_list)1758 static double get_closest_pole(struct point_list * cell_list) {
1759
1760 struct point_list_element * curr = cell_list->start;
1761 struct point_list_element * start = cell_list->start;
1762
1763 if (curr == NULL) return 1.0;
1764
1765 double max_z = 0.0;
1766
1767 do {
1768
1769 double curr_z = curr->vec_coords[2];
1770 if (fabs(curr_z) > fabs(max_z)) max_z = curr_z;
1771
1772 curr = curr->next;
1773 } while (curr != start);
1774
1775 return (max_z > 0.0)?1.0:-1.0;
1776 }
1777
1778 /*this routine is potentially being used in the CDO*/
yac_cell_lat_clipping(size_t N,struct grid_cell * cells,double lat_bounds[2],struct grid_cell * overlap_buffer)1779 void yac_cell_lat_clipping (size_t N,
1780 struct grid_cell * cells,
1781 double lat_bounds[2], // lat in rad
1782 struct grid_cell * overlap_buffer) {
1783
1784 struct point_list cell_list;
1785
1786 init_point_list(&cell_list);
1787
1788 double upper_bound, lower_bound;
1789 if (lat_bounds[0] > lat_bounds[1]) {
1790 upper_bound = sin(lat_bounds[0]);
1791 lower_bound = sin(lat_bounds[1]);
1792 } else {
1793 upper_bound = sin(lat_bounds[1]);
1794 lower_bound = sin(lat_bounds[0]);
1795 }
1796
1797 // if both bounds are at the same pole
1798 if ((fabs(-1.0 - upper_bound) < yac_angle_tol) ||
1799 (fabs(1.0 - lower_bound) < yac_angle_tol)) {
1800
1801 for (size_t n = 0; n < N; ++n)
1802 overlap_buffer[n].num_corners = 0;
1803 return;
1804 }
1805
1806 // for all source cells
1807 for (size_t n = 0; n < N; n++ ) {
1808
1809 if (cells[n].num_corners < 2) continue;
1810
1811 overlap_buffer[n].num_corners = 0;
1812
1813 enum cell_type src_cell_type = get_cell_type(cells[n]);
1814
1815 if (src_cell_type == MIXED_CELL)
1816 yac_internal_abort_message("invalid source cell type (cell contains edges consisting "
1817 "of great circles and circles of latitude)\n",
1818 __FILE__, __LINE__);
1819
1820 if (src_cell_type == LON_LAT_CELL) {
1821
1822 yac_lon_lat_cell_lat_clipping(
1823 cells + n, upper_bound, lower_bound, overlap_buffer + n);
1824
1825 } else {
1826
1827 // generate point list for current source list
1828 size_t num_corners = generate_point_list(&cell_list, cells[n]);
1829
1830 // if the cell contains a pole, we need to add this pole
1831 double pole = get_closest_pole(&cell_list);
1832 if (yac_point_in_cell(
1833 (double[3]){0.0, 0.0, pole}, cells[n])) {
1834
1835 int flag = 0;
1836
1837 // if the cell contains the north pole
1838 if (pole > 0.0) {
1839 double bound =
1840 (fabs(upper_bound - 1.0) > yac_angle_tol)?upper_bound:lower_bound;
1841
1842 for (size_t i = 0; i < num_corners; ++i) {
1843 flag |= (bound < cells[n].coordinates_xyz[i][2]);
1844 }
1845 } else {
1846 double bound =
1847 (fabs(lower_bound + 1.0) > yac_angle_tol)?lower_bound:upper_bound;
1848
1849 for (size_t i = 0; i < num_corners; ++i) {
1850 flag |= (bound > cells[n].coordinates_xyz[i][2]);
1851 }
1852 }
1853
1854 if (!flag) {
1855
1856 yac_internal_abort_message(
1857 "ERROR(yac_cell_lat_clipping): Latitude bounds are within a cell "
1858 "covering a pole, this is not supported. Increased grid resolution "
1859 "or widen lat bounds may help.", __FILE__, __LINE__);
1860 }
1861
1862 }
1863
1864 num_corners =
1865 point_list_lat_clipping(
1866 &cell_list, num_corners, upper_bound, 1);
1867 point_list_lat_clipping(&cell_list, num_corners, lower_bound, 0);
1868
1869 remove_zero_length_edges(&cell_list);
1870
1871 if (cell_list.start != NULL)
1872 generate_overlap_cell(&cell_list, overlap_buffer + n);
1873 }
1874 }
1875
1876 free_point_list(&cell_list);
1877 }
1878
1879 /* ---------------------------------------------------- */
1880
yac_correct_weights(size_t nSourceCells,double * weight)1881 void yac_correct_weights ( size_t nSourceCells, double * weight ) {
1882
1883 static const size_t maxIter = 10; // number of iterations to get better accuracy of the weights
1884
1885 double weight_diff;
1886 #ifdef YAC_VERBOSE_CLIPPING
1887 double weight_sum;
1888 #endif
1889
1890 for ( size_t iter = 1; iter < maxIter; iter++ ) {
1891
1892 weight_diff = 1.0;
1893
1894 for ( size_t n = 0; n < nSourceCells; n++ ) weight_diff -= weight[n];
1895
1896 #ifdef YAC_VERBOSE_CLIPPING
1897 for ( size_t n = 0; n < nSourceCells; n++ ) weight_sum += weight[n];
1898
1899 printf ("weight sum is %.15f \n", weight_sum);
1900 printf ("weights are ");
1901 for (size_t i = 0; i < nSourceCells; ++i)
1902 printf (" %.15f", weight[i]);
1903 printf("\n");
1904 #endif
1905
1906 if ( fabs(weight_diff) < tol ) break;
1907
1908 for ( size_t n = 0; n < nSourceCells; n++ )
1909 weight[n] += weight[n] * weight_diff;
1910 }
1911 #ifdef YAC_VERBOSE_CLIPPING
1912 if ( fabs(weight_diff) > tol )
1913 printf ("weight sum is %.15f \n", weight_sum);
1914 #endif
1915 }
1916
1917 /* ---------------------------------------------------- */
1918
get_cell_points_ordering(struct point_list * cell)1919 static int get_cell_points_ordering(struct point_list * cell) {
1920
1921 if ((cell->start == NULL) || (cell->start == cell->start->next))
1922 yac_internal_abort_message("ERROR: invalid cell\n", __FILE__, __LINE__);
1923
1924 double norm_vec[3];
1925 struct point_list_element * curr = cell->start;
1926
1927 compute_norm_vector(curr->vec_coords, curr->next->vec_coords, norm_vec);
1928
1929 curr = curr->next;
1930
1931 if (curr->next == cell->start)
1932 yac_internal_abort_message("ERROR: invalid cell\n", __FILE__, __LINE__);
1933
1934 do {
1935 curr = curr->next;
1936
1937 double dot = dotproduct(curr->vec_coords, norm_vec);
1938
1939 if (fabs(dot) > tol)
1940 return dot > 0;
1941
1942 } while (curr != cell->start);
1943
1944 return -1;
1945 }
1946
init_point_list(struct point_list * list)1947 static void init_point_list(struct point_list * list) {
1948
1949 list->start = NULL;
1950 list->free_elements = NULL;
1951 }
1952
reset_point_list(struct point_list * list)1953 static void reset_point_list(struct point_list * list) {
1954
1955 if (list->start != NULL) {
1956
1957 struct point_list_element * temp = list->start->next;
1958 list->start->next = list->free_elements;
1959 list->free_elements = temp;
1960
1961 list->start = NULL;
1962 }
1963 }
1964
remove_points(struct point_list * list)1965 static size_t remove_points(struct point_list * list) {
1966
1967 struct point_list_element * curr = list->start;
1968 struct point_list_element * start = list->start;
1969 struct point_list_element * prev = NULL;
1970
1971 if (curr == NULL) return 0;
1972
1973 // find the first point that is not to be removed
1974 while(curr->to_be_removed) {
1975 prev = curr;
1976 curr = curr->next;
1977 if (curr == start) break;
1978 };
1979
1980 // there is no point to remain
1981 if (curr->to_be_removed) {
1982 reset_point_list(list);
1983 return 0;
1984 }
1985
1986 list->start = curr;
1987 start = curr;
1988 size_t num_remaining_points = 1;
1989
1990 prev = curr;
1991 curr = curr->next;
1992
1993 while (curr != start) {
1994
1995 if (curr->to_be_removed) {
1996 prev->next = curr->next;
1997 curr->next = list->free_elements;
1998 list->free_elements = curr;
1999 curr = prev;
2000 } else {
2001 num_remaining_points++;
2002 }
2003
2004 prev = curr;
2005 curr = curr->next;
2006
2007 } while (curr != start);
2008
2009 if (list->start == list->start->next) {
2010
2011 list->start->next = list->free_elements;
2012 list->free_elements = list->start;
2013 list->start = NULL;
2014 num_remaining_points = 0;
2015 }
2016
2017 return num_remaining_points;
2018 }
2019
2020 //! returns number of edges/corners
remove_zero_length_edges(struct point_list * list)2021 static size_t remove_zero_length_edges(struct point_list * list) {
2022
2023 struct point_list_element * curr = list->start;
2024 struct point_list_element * start = list->start;
2025
2026 if (curr == NULL) return 0;
2027
2028 if (curr->next == curr) {
2029 reset_point_list(list);
2030 return 0;
2031 }
2032
2033 do {
2034
2035 // if both points are nearly identical (angle between them is very small)
2036 if (!curr->to_be_removed &&
2037 (compare_angles(
2038 get_vector_angle_2(
2039 curr->vec_coords, curr->next->vec_coords),
2040 SIN_COS_LOW_TOL) <= 0)) {
2041 curr->to_be_removed = 1;
2042 } else if (curr->edge_type == LAT_CIRCLE && curr->next->edge_type == LAT_CIRCLE) {
2043 double temp_a = atan2(curr->vec_coords[1], curr->vec_coords[0]);
2044 double temp_b = atan2(curr->next->next->vec_coords[1], curr->next->next->vec_coords[0]);
2045 if (fabs(get_angle(temp_a, temp_b)) < M_PI_2)
2046 curr->next->to_be_removed = 1;
2047 }
2048
2049 curr = curr->next;
2050 } while (curr != start);
2051
2052 return remove_points(list);
2053 }
2054
generate_point_list(struct point_list * list,struct grid_cell cell)2055 static size_t generate_point_list(struct point_list * list,
2056 struct grid_cell cell) {
2057
2058 reset_point_list(list);
2059
2060 if (cell.num_corners == 0) return 0;
2061
2062 struct point_list_element * curr = get_free_point_list_element(list);
2063
2064 list->start = curr;
2065 curr->vec_coords[0] = cell.coordinates_xyz[0][0];
2066 curr->vec_coords[1] = cell.coordinates_xyz[0][1];
2067 curr->vec_coords[2] = cell.coordinates_xyz[0][2];
2068
2069 for (size_t i = 1; i < cell.num_corners; ++i) {
2070
2071 curr->next = get_free_point_list_element(list);
2072 curr->edge_type = cell.edge_type[i - 1];
2073 curr = curr->next;
2074
2075 curr->vec_coords[0] = cell.coordinates_xyz[i][0];
2076 curr->vec_coords[1] = cell.coordinates_xyz[i][1];
2077 curr->vec_coords[2] = cell.coordinates_xyz[i][2];
2078 curr->edge_type = cell.edge_type[i];
2079 }
2080
2081 curr->next = list->start;
2082
2083 return remove_zero_length_edges(list);
2084 }
2085
2086 static struct point_list_element *
get_free_point_list_element(struct point_list * list)2087 get_free_point_list_element(struct point_list * list) {
2088
2089 struct point_list_element * element;
2090
2091 if (list->free_elements == NULL) {
2092
2093 for (int i = 0; i < 7; ++i) {
2094
2095 element = (struct point_list_element *)xmalloc(1 * sizeof(*element));
2096
2097 element->next = list->free_elements;
2098 list->free_elements = element;
2099 }
2100
2101 element = (struct point_list_element *)xmalloc(1 * sizeof(*element));
2102
2103 } else {
2104
2105 element = list->free_elements;
2106 list->free_elements = list->free_elements->next;
2107 }
2108
2109 element->next = NULL;
2110 element->to_be_removed = 0;
2111
2112 return element;
2113 }
2114
free_point_list(struct point_list * list)2115 static void free_point_list(struct point_list * list) {
2116
2117 struct point_list_element * element;
2118
2119 if (list->start != NULL) {
2120
2121 struct point_list_element * temp = list->free_elements;
2122 list->free_elements = list->start->next;
2123 list->start->next = temp;
2124 }
2125
2126 while (list->free_elements != NULL) {
2127
2128 element = list->free_elements;
2129 list->free_elements = element->next;
2130 free(element);
2131 }
2132
2133 list->start = NULL;
2134 list->free_elements = NULL;
2135 }
2136
is_empty_gc_cell(struct point_list * list,size_t num_edges)2137 static int is_empty_gc_cell(struct point_list * list, size_t num_edges) {
2138
2139 #define NORM_TOL (1e-6)
2140
2141 if ((num_edges == 2) &&
2142 !((list->start->edge_type == LAT_CIRCLE) ^
2143 (list->start->next->edge_type == LAT_CIRCLE)))
2144 return 1;
2145
2146 struct point_list_element * curr = list->start;
2147
2148 for (size_t i = 0; i < num_edges; ++i) {
2149
2150 if (curr->edge_type == LAT_CIRCLE) return 0;
2151 curr = curr->next;
2152 }
2153
2154 double ref_norm[3];
2155
2156 compute_norm_vector(curr->vec_coords, curr->next->vec_coords, ref_norm);
2157
2158 for (size_t i = 0; i < num_edges-1; ++i) {
2159
2160 curr = curr->next;
2161
2162 double norm[3];
2163
2164 compute_norm_vector(curr->vec_coords, curr->next->vec_coords, norm);
2165
2166 if (((fabs(ref_norm[0] - norm[0]) > NORM_TOL) ||
2167 (fabs(ref_norm[1] - norm[1]) > NORM_TOL) ||
2168 (fabs(ref_norm[2] - norm[2]) > NORM_TOL)) &&
2169 ((fabs(ref_norm[0] + norm[0]) > NORM_TOL) ||
2170 (fabs(ref_norm[1] + norm[1]) > NORM_TOL) ||
2171 (fabs(ref_norm[2] + norm[2]) > NORM_TOL)))
2172 return 0;
2173 }
2174
2175 return 1;
2176 #undef NORM_TOL
2177 }
2178
generate_overlap_cell(struct point_list * list,struct grid_cell * cell)2179 static void generate_overlap_cell(struct point_list * list,
2180 struct grid_cell * cell) {
2181
2182 //! \todo test whether all points of the cell are on a single
2183 //! great circle --> empty cell
2184
2185 size_t num_edges = remove_zero_length_edges(list);
2186
2187 if ((num_edges < 2) ||
2188 is_empty_gc_cell(list, num_edges)){
2189
2190 reset_point_list(list);
2191 cell->num_corners = 0;
2192 return;
2193 }
2194
2195 if (num_edges > cell->array_size) {
2196 free(cell->coordinates_xyz);
2197 free(cell->edge_type);
2198 cell->coordinates_xyz = (double(*)[3])xmalloc(num_edges * sizeof(*cell->coordinates_xyz));
2199 cell->edge_type = (enum yac_edge_type *)xmalloc(num_edges * sizeof(*cell->edge_type));
2200 cell->array_size = num_edges;
2201 }
2202 cell->num_corners = num_edges;
2203
2204 struct point_list_element * curr = list->start;
2205
2206 for (size_t i = 0; i < num_edges; ++i) {
2207
2208 cell->coordinates_xyz[i][0] = curr->vec_coords[0];
2209 cell->coordinates_xyz[i][1] = curr->vec_coords[1];
2210 cell->coordinates_xyz[i][2] = curr->vec_coords[2];
2211 cell->edge_type[i] = curr->edge_type;
2212
2213 curr = curr->next;
2214 }
2215 }
2216
2217