YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
interp_method_avg.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 <string.h>
11
13#include "grid_cell.h"
14#include "interp_method_avg.h"
16
17#define SWAP(x, y, T) do { T SWAP = x; x = y; y = SWAP; } while (0)
18
19static size_t do_search_avg(struct interp_method * method,
20 struct yac_interp_grid * interp_grid,
21 size_t * tgt_points, size_t count,
22 struct yac_interp_weights * weights);
23static void delete_avg(struct interp_method * method);
24
25typedef int (*func_compute_weights)(
26 double[3], size_t, yac_const_coordinate_pointer, int*, double*,
27 enum yac_cell_type cell_type);
28
29static struct interp_method_vtable
33
39
41 double * point, size_t field_cell, size_t cell_size,
42 size_t const * vertex_to_cell, size_t const * vertex_to_cell_offsets,
43 yac_const_coordinate_pointer field_coordinates,
44 struct yac_grid_cell * cell_buffer) {
45
46 if (cell_size == 0) return 0;
47
48 cell_buffer->num_corners = cell_size;
49
50 yac_coordinate_pointer coordinates_xyz;
51
52 if (cell_size > cell_buffer->array_size) {
53 cell_buffer->coordinates_xyz =
54 ((coordinates_xyz =
56 cell_buffer->coordinates_xyz,
57 cell_size * sizeof(*coordinates_xyz))));
58 cell_buffer->edge_type =
60 cell_buffer->edge_type, cell_size * sizeof(cell_buffer->edge_type));
61 for (size_t i = 0; i < cell_size; ++i)
62 cell_buffer->edge_type[i] = YAC_GREAT_CIRCLE_EDGE;
63 cell_buffer->array_size = cell_size;
64 } else {
65 coordinates_xyz = cell_buffer->coordinates_xyz;
66 }
67
68 size_t const * field_cell_vertices =
69 vertex_to_cell + vertex_to_cell_offsets[field_cell];
70 for (size_t i = 0; i < cell_size; ++i)
71 memcpy(coordinates_xyz[i], field_coordinates[field_cell_vertices[i]],
72 sizeof(*coordinates_xyz));
73
74 return yac_point_in_cell(point, *cell_buffer);
75}
76
77#define IS_GC(x) (((x) == YAC_GREAT_CIRCLE_EDGE) || ((x) == YAC_LON_CIRCLE_EDGE))
78#define IS_LON_LAT(x) (((x) == YAC_LAT_CIRCLE_EDGE) || ((x) == YAC_LON_CIRCLE_EDGE))
79
81 enum yac_edge_type const * edge_types,
82 size_t const * edge_indices, size_t num_edges) {
83
84 if (edge_types == NULL) return YAC_GREAT_CIRCLE_CELL;
85
86 int edge_type_flag = 0;
87
88 if (num_edges == 4) {
89
90 enum yac_edge_type temp_edges[4] =
91 {edge_types[edge_indices[0]],
92 edge_types[edge_indices[1]],
93 edge_types[edge_indices[2]],
94 edge_types[edge_indices[3]]};
95 if (IS_LON_LAT(temp_edges[0]) &&
96 IS_LON_LAT(temp_edges[1]) &&
97 (temp_edges[0] == temp_edges[2]) &&
98 (temp_edges[1] == temp_edges[3]) &&
99 (temp_edges[0] != temp_edges[1]))
100 return YAC_LON_LAT_CELL;
101 }
102
103 for (size_t i = 0; i < num_edges; ++i)
104 edge_type_flag |= (1 << (edge_types[edge_indices[i]] == YAC_LAT_CIRCLE_EDGE));
105
106 enum yac_cell_type cell_types[4] =
111 return cell_types[edge_type_flag];
112}
113#undef IS_LON_LAT
114#undef IS_GC
115
116static size_t do_search_avg (struct interp_method * method,
117 struct yac_interp_grid * interp_grid,
118 size_t * tgt_points, size_t count,
119 struct yac_interp_weights * weights) {
120
121 struct interp_method_avg * method_avg = (struct interp_method_avg *)method;
123
125 yac_interp_grid_get_num_src_fields(interp_grid) == 1,
126 "ERROR(do_search_avg): invalid number of source fields")
127
128 enum yac_location src_field_location =
130
132 (src_field_location == YAC_LOC_CORNER) ||
133 (src_field_location == YAC_LOC_CELL),
134 "ERROR(do_search_avg): unsupported source field location type")
135
136 // get coordinates of target points
137 yac_coordinate_pointer tgt_coords = xmalloc(count * sizeof(*tgt_coords));
139 interp_grid, tgt_points, count, tgt_coords);
140
141 size_t * size_t_buffer = xmalloc(2 * count * sizeof(*size_t_buffer));
142 size_t * src_field_cells = size_t_buffer;
143 size_t * reorder_idx = size_t_buffer + count;
144
145 // get matching source cells for all target points
147 interp_grid, tgt_coords, count, src_field_cells);
148
149 size_t const * src_field_cell_to_vertex;
150 size_t const * src_field_cell_to_vertex_offsets;
151 size_t const * src_field_cell_to_edge;
152 size_t const * src_field_cell_to_edge_offsets;
153 int const * src_field_num_vertices_per_cell;
154 enum yac_edge_type const * src_edge_types;
155
156 // if the source field is location at cell points
157 if (src_field_location == YAC_LOC_CELL) {
158
159 // generate auxiliary grid for all search result cells
161 interp_grid, src_field_cells, count,
162 (size_t**)&src_field_cell_to_vertex,
163 (size_t**)&src_field_cell_to_vertex_offsets,
164 (int**)&src_field_num_vertices_per_cell);
165
166 struct yac_const_basic_grid_data * src_grid_data =
168 yac_const_coordinate_pointer src_field_coordinates =
170
171 struct yac_grid_cell cell_buffer;
172 yac_init_grid_cell(&cell_buffer);
173
174 // for all target point
175 for (size_t i = 0; i < count; ++i) {
176
177 size_t curr_src_cell = src_field_cells[i];
178 if (curr_src_cell == SIZE_MAX) continue;
179
180 double * curr_tgt_coord = tgt_coords[i];
181 size_t const * curr_vertices =
182 src_grid_data->cell_to_vertex +
183 src_grid_data->cell_to_vertex_offsets[curr_src_cell];
184 size_t curr_num_vertices =
185 src_grid_data->num_vertices_per_cell[curr_src_cell];
186
187 size_t result_src_field_cell = SIZE_MAX;
188
189 // for all auxiliary of the current source result cell
190 for (size_t j = 0; j < curr_num_vertices; ++j) {
191
192 size_t src_field_cell = curr_vertices[j];
193 size_t src_field_cell_size =
194 (size_t)src_field_num_vertices_per_cell[src_field_cell];
195
196 // check whether the target point is in the current source field cell
198 curr_tgt_coord, src_field_cell, src_field_cell_size,
199 src_field_cell_to_vertex, src_field_cell_to_vertex_offsets,
200 src_field_coordinates, &cell_buffer)) {
201 result_src_field_cell = src_field_cell;
202 break;
203 }
204 }
205
206 src_field_cells[i] = result_src_field_cell;
207 }
208 yac_free_grid_cell(&cell_buffer);
209
210 src_edge_types = NULL;
211 src_field_cell_to_edge = NULL;
212 src_field_cell_to_edge_offsets = NULL;
213 } else {
214
215 struct yac_const_basic_grid_data * grid_data =
217 src_field_cell_to_vertex = grid_data->cell_to_vertex;
218 src_field_cell_to_vertex_offsets = grid_data->cell_to_vertex_offsets;
219 src_field_num_vertices_per_cell = grid_data->num_vertices_per_cell;
220 src_edge_types = grid_data->edge_type;
221 src_field_cell_to_edge = grid_data->cell_to_edge;
222 src_field_cell_to_edge_offsets = grid_data->cell_to_edge_offsets;
223 }
224 yac_const_coordinate_pointer src_field_coordinates =
226 const_yac_int_pointer src_global_ids =
228
229 // sort target points, for which we found a source cell, to the beginning of
230 // the array
231 for (size_t i = 0; i < count; ++i) reorder_idx[i] = i;
232 yac_quicksort_index_size_t_size_t(src_field_cells, count, reorder_idx);
233
234 size_t result_count = 0;
235 for (result_count = 0; result_count < count; ++result_count)
236 if (src_field_cells[result_count] == SIZE_MAX) break;
237
238 size_t total_num_weights = 0;
239 size_t max_num_vertices_per_cell = 0;
240 size_t * num_weights_per_tgt =
241 xmalloc(result_count * sizeof(*num_weights_per_tgt));
242
243 // get the number of vertices per target cell and the total number of vertices
244 // required for the interpolation (may contain duplicated vertices)
245 for (size_t i = 0; i < result_count; ++i) {
246 size_t curr_num_vertices =
247 (size_t)(src_field_num_vertices_per_cell[src_field_cells[i]]);
248 num_weights_per_tgt[i] = curr_num_vertices;
249 total_num_weights += curr_num_vertices;
250 if (curr_num_vertices > max_num_vertices_per_cell)
251 max_num_vertices_per_cell = curr_num_vertices;
252 }
253 if (src_field_location == YAC_LOC_CELL)
254 free((void*)src_field_num_vertices_per_cell);
255
256 const_int_pointer src_field_mask =
258
259 yac_coordinate_pointer src_coord_buffer =
260 xmalloc(max_num_vertices_per_cell * sizeof(*src_coord_buffer));
261 int * mask_buffer =
262 (src_field_mask != NULL)?
263 xmalloc(max_num_vertices_per_cell * sizeof(*mask_buffer)):NULL;
264 double * w = xmalloc(total_num_weights * sizeof(*w));
265 size_t * src_points = xmalloc(total_num_weights * sizeof(*src_points));
266
267 // for each target point, extract relevant source data and compute the weights
268 // based on that
269 total_num_weights = 0;
270 for (size_t i = 0; i < result_count;) {
271
272 size_t curr_num_vertices = num_weights_per_tgt[i];
273
274 const_size_t_pointer curr_cell_to_vertex =
275 src_field_cell_to_vertex +
276 src_field_cell_to_vertex_offsets[src_field_cells[i]];
277 const_size_t_pointer curr_cell_to_edge =
278 (src_edge_types != NULL)?
279 (src_field_cell_to_edge +
280 src_field_cell_to_edge_offsets[src_field_cells[i]]):NULL;
281 double * curr_weights = w + total_num_weights;
282
283 // get the index of the source cell corner with the lowest global id
284 // (this is being used, that the order in which the source corners are
285 // processes is decomposition independent)
286 size_t lowest_global_id_idx = 0;
287 {
288 yac_int lowest_global_id = src_global_ids[curr_cell_to_vertex[0]];
289 for (size_t j = 1; j < curr_num_vertices; ++j) {
290 if (src_global_ids[curr_cell_to_vertex[j]] < lowest_global_id) {
291 lowest_global_id = src_global_ids[curr_cell_to_vertex[j]];
292 lowest_global_id_idx = j;
293 }
294 }
295 }
296
297 // get the mask and the source coordinates of the source vertices required
298 // for the current target point
299 for (size_t j = 0, l = lowest_global_id_idx; j < curr_num_vertices;
300 ++j, ++total_num_weights, ++l) {
301
302 if (l == curr_num_vertices) l = 0;
303
304 size_t curr_vertex_idx = curr_cell_to_vertex[l];
305 src_points[total_num_weights] = curr_vertex_idx;
306 for (size_t k = 0; k < 3; ++k)
307 src_coord_buffer[j][k] = src_field_coordinates[curr_vertex_idx][k];
308 if (src_field_mask != NULL)
309 mask_buffer[j] = src_field_mask[curr_vertex_idx];
310 }
311
312 enum yac_cell_type cell_type =
314 src_edge_types, curr_cell_to_edge, curr_num_vertices);
315
316 // compute the weights
317 if (compute_weights(tgt_coords[reorder_idx[i]], curr_num_vertices,
318 (yac_const_coordinate_pointer)src_coord_buffer,
319 mask_buffer, curr_weights, cell_type)) {
320
321 // if weight computation was successful
322 num_weights_per_tgt[i] = curr_num_vertices;
323 ++i;
324 } else {
325
326 // if no weights could be computed
327 total_num_weights -= curr_num_vertices;
328 result_count--;
329 src_field_cells[i] = src_field_cells[result_count];
330 size_t temp_reorder_idx = reorder_idx[i];
331 reorder_idx[i] = reorder_idx[result_count];
332 reorder_idx[result_count] = temp_reorder_idx;
333 }
334 }
335 if (src_field_location == YAC_LOC_CELL) {
336 free((void*)src_field_cell_to_vertex);
337 free((void*)src_field_cell_to_vertex_offsets);
338 }
339
340 // move the non-interpolated target points to the end
341 for (size_t i = 0; i < count; ++i) src_field_cells[reorder_idx[i]] = i;
342 yac_quicksort_index_size_t_size_t(src_field_cells, count, tgt_points);
343
344 free(tgt_coords);
345 free(mask_buffer);
346 free(src_coord_buffer);
347 free(size_t_buffer);
348
349 struct remote_points tgts = {
350 .data =
352 interp_grid, tgt_points, result_count),
353 .count = result_count};
354 struct remote_point * srcs =
356 interp_grid, 0, src_points, total_num_weights);
357
358 // store weights
360 weights, &tgts, num_weights_per_tgt, srcs, w);
361
362 free(tgts.data);
363 free(srcs);
364 free(w);
365 free(num_weights_per_tgt);
366 free(src_points);
367
368 return result_count;
369}
370
372 double tgt_coords[3], size_t num_vertices,
373 yac_const_coordinate_pointer src_coords, int * src_mask, double * weights,
374 enum yac_cell_type cell_type) {
375
376 UNUSED(tgt_coords);
377 UNUSED(cell_type);
378 UNUSED(src_coords);
379
380 size_t num_unmasked_points = 0;
381
382 // if we have a source mask, check whether any points is masked
383 if (src_mask != NULL) {
384 for (size_t i = 0; i < num_vertices; ++i)
385 if (src_mask[i]) num_unmasked_points++;
386 if (num_unmasked_points == 0) return 0;
387
388 double weight = 1.0 / (double)num_unmasked_points;
389 for (size_t i = 0; i < num_vertices; ++i)
390 weights[i] = (src_mask[i])?weight:0.0;
391
392 } else {
393 double weight = 1.0 / (double)num_vertices;
394 for (size_t i = 0; i < num_vertices; ++i) weights[i] = weight;
395 }
396 return 1;
397}
398
400 double tgt_coords[3], size_t num_vertices,
401 yac_const_coordinate_pointer src_coords, int * src_mask, double * weights,
402 enum yac_cell_type cell_type) {
403
404 UNUSED(tgt_coords);
405 UNUSED(cell_type);
406 UNUSED(src_coords);
407
408 // if we have a source mask, check whether any points is masked
409 if (src_mask != NULL)
410 for (size_t i = 0; i < num_vertices; ++i)
411 if (!src_mask[i]) return 0;
412
413 double weight = 1.0 / (double)num_vertices;
414 for (size_t i = 0; i < num_vertices; ++i) weights[i] = weight;
415
416 return 1;
417}
418
420 double tgt_coords[3], size_t num_vertices,
421 yac_const_coordinate_pointer src_coords, int * src_mask, double * weights,
422 enum yac_cell_type cell_type) {
423
424 UNUSED(cell_type);
425
426 // if there is a source mask, check if there are unmasked vertices
427 if (src_mask != NULL) {
428 int has_unmasked = 0;
429 for (size_t i = 0; i < num_vertices; ++i) has_unmasked |= src_mask[i];
430 if (!has_unmasked) return 0;
431 }
432
433 if (src_mask != NULL) {
434 for (size_t i = 0; i < num_vertices; ++i) {
435
436 if (src_mask[i]) {
437
438 double distance =
439 get_vector_angle(tgt_coords, (double*)(src_coords[i]));
440
441 // if the target and source point are nearly identical
442 if (distance < yac_angle_tol) {
443 for (size_t j = 0; j < num_vertices; ++j) weights[j] = 0.0;
444 weights[i] = 1.0;
445 return 1;
446 }
447
448 weights[i] = 1.0 / distance;
449 } else {
450 weights[i] = 0.0;
451 }
452 }
453 } else {
454 for (size_t i = 0; i < num_vertices; ++i) {
455
456 double distance =
457 get_vector_angle(tgt_coords, (double*)(src_coords[i]));
458
459 // if the target and source point are nearly identical
460 if (distance < yac_angle_tol) {
461 for (size_t j = 0; j < num_vertices; ++j) weights[j] = 0.0;
462 weights[i] = 1.0;
463 return 1;
464 }
465
466 weights[i] = 1.0 / distance;
467 }
468 }
469
470 // compute scaling factor for the weights
471 double inv_distance_sum = 0.0;
472 for (size_t i = 0; i < num_vertices; ++i)
473 inv_distance_sum += weights[i];
474 double scale = 1.0 / inv_distance_sum;
475
476 for (size_t i = 0; i < num_vertices; ++i) weights[i] *= scale;
477
478 return 1;
479}
480
482 double tgt_coords[3], size_t num_vertices,
483 yac_const_coordinate_pointer src_coords, int * src_mask, double * weights,
484 enum yac_cell_type cell_type) {
485
486 UNUSED(cell_type);
487
488 for (size_t i = 0; i < num_vertices; ++i) {
489
490 double distance =
491 get_vector_angle(tgt_coords, (double*)(src_coords[i]));
492
493 // if the target and source point are nearly identical
494 if (distance < yac_angle_tol) {
495 if ((src_mask != NULL) && !src_mask[i]) return 0;
496 for (size_t j = 0; j < num_vertices; ++j) weights[j] = 0.0;
497 weights[i] = 1.0;
498 return 1;
499 }
500
501 weights[i] = 1.0 / distance;
502 }
503
504 // if there is a source mask, check if there are masked vertices
505 // (we do this here because there may be a matching source point
506 // that is not masked)
507 if (src_mask != NULL)
508 for (size_t i = 0; i < num_vertices; ++i) if(!src_mask[i]) return 0;
509
510 // compute scaling factor for the weights
511 double inv_distance_sum = 0.0;
512 for (size_t i = 0; i < num_vertices; ++i)
513 inv_distance_sum += weights[i];
514 double scale = 1.0 / inv_distance_sum;
515
516 for (size_t i = 0; i < num_vertices; ++i) weights[i] *= scale;
517
518 return 1;
519}
520
521// compute the spherical barycentric coordinates of the given vertices with
522// respect to the given triangle
524 double * barycentric_coords, size_t triangle_indices[3],
525 yac_const_coordinate_pointer grid_coords) {
526
527 double A[3][3];
528 lapack_int n = 3, nrhs = 1, lda = n, ldx = n, ipiv[3];
529 for (int i = 0; i < 3; ++i)
530 memcpy(A[i], grid_coords[triangle_indices[i]], sizeof(*grid_coords));
531
532 // for a vertex v the spherical barycentric coordinates b are defined as
533 // follows
534 // A * b = v
535 // where: A is the matrix consisting of the vertex coordinates of the three
536 // corners of the triangle
537 // we compute b by solving this linear system using LAPACK
539 !LAPACKE_dgesv(
540 LAPACK_COL_MAJOR, n, nrhs, &A[0][0], lda,
541 ipiv, barycentric_coords, ldx),
542 "ERROR: internal error (could not solve linear 3x3 system)")
543}
544
545// returns 0 if edge_direction_vec points south; 1 otherwise
546static inline int get_lat_edge_ordering(double const * a, double const * b) {
547
548 return (a[0] * b[1] - a[1] * b[0]) < 0.0;
549}
550
553 double lon[2], double lat[2], int reorder[4]) {
554
555 // determine whether the lat neighbour to vertex 0 is at position 1 or 3
556 int lat_neigh_offset =
557 (fabs(coords[0][2] - coords[1][2]) <
558 fabs(coords[0][2] - coords[3][2]))?1:3;
559 // get the index of the coord which is closes to the pole
560 // (the one with the lowest absolut z coordinate)
561 int closest_to_equator_idx =
562 (fabs(coords[0][2]) < fabs(coords[1][2]))?0:2;
563 int upper_edge_idx = ((coords[0][2]) > coords[2][2])?0:2;
564 int lower_edge_idx = upper_edge_idx^2;
565 int upper_edge_is_closer_to_pole =
566 closest_to_equator_idx == upper_edge_idx;
567 int closest_to_equator_edge_ordering =
569 coords[closest_to_equator_idx],
570 coords[(closest_to_equator_idx+lat_neigh_offset)%4]);
571 int cell_ordering =
572 closest_to_equator_edge_ordering ^ upper_edge_is_closer_to_pole;
573
574 XYZtoLL(
575 coords[closest_to_equator_idx], &lon[closest_to_equator_edge_ordering],
576 &lat[upper_edge_is_closer_to_pole]);
577 XYZtoLL(
578 coords[(closest_to_equator_idx + lat_neigh_offset)%4],
579 &lon[closest_to_equator_edge_ordering^1],
580 &lat[upper_edge_is_closer_to_pole]);
581 lat[!upper_edge_is_closer_to_pole] =
582 M_PI_2 - acos(coords[closest_to_equator_idx ^ 2][2]);
583
584 reorder[cell_ordering] = lower_edge_idx;
585 reorder[cell_ordering^1] = (lower_edge_idx+lat_neigh_offset)%4;
586 reorder[2+cell_ordering] = upper_edge_idx;
587 reorder[2+(cell_ordering^1)] = (upper_edge_idx+lat_neigh_offset)%4;
588
589 if (lon[1] < lon[0]) lon[1] += 2.0 * M_PI;
590}
591
593 double point_coord[3], double cell_lon[2], double cell_lat[2],
594 double * point_lon, double * point_lat) {
595
596 UNUSED(cell_lat);
597
598 double lon, lat;
599 XYZtoLL(point_coord, &lon, &lat);
600
601 // adjust target point lon to source cell lon
602 if (lon < cell_lon[0]) {
603 while (fabs(cell_lon[0] - lon) > M_PI) lon += 2.0 * M_PI;
604 } else {
605 while (fabs(lon - cell_lon[1]) > M_PI) lon -= 2.0 * M_PI;
606 }
607
608 *point_lon = lon;
609 *point_lat = lat;
610}
611
613 double point_lon, double point_lat, double cell_lon[2], double cell_lat[2]) {
614
615 return ((cell_lat[1] - cell_lat[0]) * (point_lat - cell_lat[0]) -
616 (cell_lat[1] - cell_lat[0]) * (point_lon - cell_lon[0])) > 0.0;
617}
618
620 double tgt_coords[3],
621 yac_const_coordinate_pointer src_coords, double * weights) {
622
623 // get lon lat bounds of the cell and the target point
624 double src_lon[2], src_lat[2], tgt_lon, tgt_lat;
625 int src_reorder[4];
626 get_cell_lon_lat_bounds(src_coords, src_lon, src_lat, src_reorder);
627 get_point_lon_lat(tgt_coords, src_lon, src_lat, &tgt_lon, &tgt_lat);
628
629 int triangle_idx =
630 determine_triangle_idx(tgt_lon, tgt_lat, src_lon, src_lat);
631
632 if (triangle_idx) {
633
634 double w[2];
635 w[0] = (tgt_lat - src_lat[1]) / (src_lat[0] - src_lat[1]);
636 w[1] = ((src_lat[1] - src_lat[0]) * (tgt_lon - src_lon[0])) /
637 ((src_lon[0] - src_lon[1]) * (src_lat[0] - src_lat[1]));
638
639 weights[src_reorder[0]] = w[0];
640 weights[src_reorder[2]] = w[1];
641 weights[src_reorder[3]] = 1.0 - (w[0] + w[1]);
642 weights[src_reorder[1]] = 0.0;
643 } else {
644
645 double w[2];
646 w[0] = (tgt_lon - src_lon[1]) / (src_lon[0] - src_lon[1]);
647 w[1] = ((src_lat[1] - src_lat[0]) * (tgt_lon - src_lon[1]) +
648 (src_lon[0] - src_lon[1]) * (tgt_lat - src_lat[1])) /
649 ((src_lat[0] - src_lat[1]) * (src_lon[0] - src_lon[1]));
650
651 weights[src_reorder[0]] = w[0];
652 weights[src_reorder[1]] = w[1];
653 weights[src_reorder[2]] = 1.0 - (w[0] + w[1]);
654 weights[src_reorder[3]] = 0.0;
655 }
656}
657
659 double tgt_coords[3], size_t num_vertices,
660 yac_const_coordinate_pointer src_coords, int * src_mask, double * weights,
661 enum yac_cell_type cell_type) {
662
663 // if all source points are masked
664 if (src_mask != NULL) {
665 size_t i;
666 for (i = 0; i < num_vertices; ++i) if (src_mask[i] != 0) break;
667 if (i == num_vertices) return 0;
668 }
669
670 double const tol = 1e-9;
671
673 (cell_type == YAC_LON_LAT_CELL) ||
674 (cell_type == YAC_GREAT_CIRCLE_CELL),
675 "ERROR(compute_weights_bary_yes): barycentric coordinates "
676 "are only supported for great circle edge cells and lon lat cells")
677
678 if (cell_type == YAC_LON_LAT_CELL) {
679
680 // compute weights
681 compute_weights_bary_reg(tgt_coords, src_coords, weights);
682
683 // apply mask handling
684 if (src_mask != NULL) {
685 double weight_sum = 1.0;
686
687 // check if any of the source points is masked
688 for (int i = 0; i < 4; ++i) if (!src_mask[i]) {
689 weight_sum -= weights[i];
690 weights[i] = 0.0;
691 }
692
693 // if all relevant points are masked
694 if (weight_sum < tol) return 0;
695
696 // if only some of the relevant points are masked
697 if (weight_sum < 1.0) {
698 double scale = 1.0 / weight_sum;
699 for (int i = 0; i < 4; ++i) weights[i] *= scale;
700 }
701 }
702
703 return 1;
704
705 } else {
706
707 size_t corner_indices[num_vertices];
708 size_t triangle_indices[num_vertices-2][3];
709
710 for (size_t i = 0; i < num_vertices; ++i) corner_indices[i] = i;
711
712 // triangulate source polygon
713 if (num_vertices > 3) {
715 corner_indices, num_vertices, 0, triangle_indices);
716 } else {
717 for (size_t i = 0; i < 3; ++i) triangle_indices[0][i] = i;
718 }
719
720 // find the best matching triangle
721 double min_barycentric_coord = -DBL_MAX;
722 double barycentric_coords[3];
723 size_t match_index = SIZE_MAX;
724
725 for (size_t i = 0; i < num_vertices - 2; ++i) {
726
727 double temp_barycentric_coords[3];
728 memcpy(temp_barycentric_coords, tgt_coords, 3 * sizeof(double));
729
730 // compute barycentric coordinates for the current triangle
732 temp_barycentric_coords, triangle_indices[i], src_coords);
733
734 double curr_min_barycentric_coord =
735 MIN(temp_barycentric_coords[0],
736 MIN(temp_barycentric_coords[1],
737 temp_barycentric_coords[2]));
738
739 if (curr_min_barycentric_coord > min_barycentric_coord) {
740 min_barycentric_coord = curr_min_barycentric_coord;
741 match_index = i;
742 for (int j = 0; j < 3; ++j)
743 barycentric_coords[j] = MAX(0.0, temp_barycentric_coords[j]);
744 }
745 }
746
747 // initialise weights
748 if (num_vertices > 3)
749 for (size_t j = 0; j < num_vertices; ++j) weights[j] = 0.0;
750
751 // apply mask
752 if (src_mask != NULL)
753 for (int j = 0; j < 3; ++j)
754 if (src_mask[triangle_indices[match_index][j]] == 0)
755 barycentric_coords[j] = 0.0;
756
757 double weight_sum = 0.0;
758 for (int j = 0; j < 3; ++j) {
759 if (barycentric_coords[j] < 1e-4) barycentric_coords[j] = 0.0;
760 else weight_sum += barycentric_coords[j];
761 }
762
763 if (weight_sum < tol) return 0;
764
765 double scale = 1.0 / weight_sum;
766
767 for (int j = 0; j < 3; ++j)
768 weights[triangle_indices[match_index][j]] =
769 barycentric_coords[j] * scale;
770
771 return 1;
772 }
773}
774
776 double tgt_coords[3], size_t num_vertices,
777 yac_const_coordinate_pointer src_coords, int * src_mask, double * weights,
778 enum yac_cell_type cell_type) {
779
780 // if all source points are masked
781 if (src_mask != NULL) {
782 size_t i;
783 for (i = 0; i < num_vertices; ++i) if (src_mask[i] != 0) break;
784 if (i == num_vertices) return 0;
785 }
786
787 double const tol = 1e-9;
788
790 (cell_type == YAC_LON_LAT_CELL) || (cell_type == YAC_GREAT_CIRCLE_CELL),
791 "ERROR(compute_weights_bary_no): barycentric coordinates "
792 "are only supported for great circle edge cells and lon lat cells")
793
794 if (cell_type == YAC_LON_LAT_CELL) {
795
796 // compute weights
797 compute_weights_bary_reg(tgt_coords, src_coords, weights);
798
799 for (int i = 0; i < 4; ++i) if (weights[i] < tol) weights[i] = 0.0;
800
801 // check if any of the source points is masked
802 if (src_mask != NULL)
803 for (int i = 0; i < 4; ++i)
804 if (!src_mask[i] && (weights[i] != 0.0))
805 return 0;
806 return 1;
807
808 } else {
809
810 size_t corner_indices[num_vertices];
811 size_t triangle_indices[num_vertices-2][3];
812
813 for (size_t i = 0; i < num_vertices; ++i) corner_indices[i] = i;
814
815 // triangulate source polygon
816 if (num_vertices > 3) {
818 corner_indices, num_vertices, 0, triangle_indices);
819 } else {
820 for (size_t i = 0; i < 3; ++i) triangle_indices[0][i] = i;
821 }
822
823 // find the best matching triangle
824 double min_barycentric_coord = -DBL_MAX;
825 double barycentric_coords[3];
826 size_t match_index = SIZE_MAX;
827
828 for (size_t i = 0; i < num_vertices - 2; ++i) {
829
830 double temp_barycentric_coords[3];
831 memcpy(temp_barycentric_coords, tgt_coords, 3 * sizeof(double));
832
833 // compute barycentric coordinates for the current triangle
835 temp_barycentric_coords, triangle_indices[i], src_coords);
836
837 double curr_min_barycentric_coord =
838 MIN(temp_barycentric_coords[0],
839 MIN(temp_barycentric_coords[1],
840 temp_barycentric_coords[2]));
841
842 if (curr_min_barycentric_coord > min_barycentric_coord) {
843 min_barycentric_coord = curr_min_barycentric_coord;
844 match_index = i;
845 for (int j = 0; j < 3; ++j)
846 barycentric_coords[j] = MAX(0.0, temp_barycentric_coords[j]);
847 }
848 }
849
850 // initialise weights
851 if (num_vertices > 3)
852 for (size_t j = 0; j < num_vertices; ++j) weights[j] = 0.0;
853
854 // apply mask
855 if (src_mask != NULL)
856 for (int j = 0; j < 3; ++j)
857 if (!src_mask[triangle_indices[match_index][j]] &&
858 (barycentric_coords[j] >= 1e-4)) return 0;
859
860 double weight_sum = 0.0;
861 for (int j = 0; j < 3; ++j) {
862 if (barycentric_coords[j] < 1e-4) barycentric_coords[j] = 0.0;
863 else weight_sum += barycentric_coords[j];
864 }
865
866 if (weight_sum < tol) return 0;
867
868 double scale = 1.0 / weight_sum;
869
870 for (int j = 0; j < 3; ++j)
871 weights[triangle_indices[match_index][j]] =
872 barycentric_coords[j] * scale;
873
874 return 1;
875 }
876}
877
879 enum yac_interp_avg_weight_type weight_type, int partial_coverage) {
880
885 "ERROR(select_compute_weight_routine): invalid weight type")
886
887 switch(weight_type) {
889 return (partial_coverage)?compute_weights_avg_yes:compute_weights_avg_no;
891 return (partial_coverage)?compute_weights_dist_yes:compute_weights_dist_no;
893 default:
894 return (partial_coverage)?compute_weights_bary_yes:compute_weights_bary_no;
895 };
896}
897
900 int partial_coverage) {
901
902 struct interp_method_avg * method = xmalloc(1 * sizeof(*method));
903
905 method->compute_weights =
907
908 return (struct interp_method*)method;
909}
910
911static void delete_avg(struct interp_method * method) {
912 free(method);
913}
int yac_point_in_cell(double point_coords[3], struct yac_grid_cell cell)
static double const tol
enum yac_interp_ncc_weight_type weight_type
#define UNUSED(x)
Definition core.h:71
int const * const_int_pointer
size_t const *const const_size_t_pointer
yac_int const * const_yac_int_pointer
void yac_triangulate_cell_indices(size_t const *corner_indices, size_t num_corners, size_t start_corner, size_t(*triangle_indices)[3])
Definition grid_cell.c:144
#define yac_angle_tol
Definition geometry.h:34
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:364
static void XYZtoLL(double const p_in[], double *lon, double *lat)
Definition geometry.h:318
void yac_init_grid_cell(struct yac_grid_cell *cell)
Definition grid_cell.c:14
void yac_free_grid_cell(struct yac_grid_cell *cell)
Definition grid_cell.c:44
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
void yac_interp_grid_do_points_search(struct yac_interp_grid *interp_grid, yac_coordinate_pointer search_coords, size_t count, size_t *src_cells)
size_t yac_interp_grid_get_num_src_fields(struct yac_interp_grid *interp_grid)
const_int_pointer yac_interp_grid_get_src_field_mask(struct yac_interp_grid *interp_grid, size_t src_field_idx)
void yac_interp_grid_get_aux_grid_src(struct yac_interp_grid *interp_grid, size_t *cells, size_t count, size_t **vertex_to_cell, size_t **vertex_to_cell_offsets, int **num_cells_per_vertex)
const_yac_int_pointer yac_interp_grid_get_src_field_global_ids(struct yac_interp_grid *interp_grid, size_t src_field_idx)
struct remote_point * yac_interp_grid_get_tgt_remote_points(struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count)
enum yac_location yac_interp_grid_get_src_field_location(struct yac_interp_grid *interp_grid, size_t src_field_idx)
struct remote_point * yac_interp_grid_get_src_remote_points(struct yac_interp_grid *interp_grid, size_t src_field_idx, size_t *src_points, size_t count)
yac_const_coordinate_pointer yac_interp_grid_get_src_field_coords(struct yac_interp_grid *interp_grid, size_t src_field_idx)
void yac_interp_grid_get_tgt_coordinates(struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, yac_coordinate_pointer tgt_coordinates)
struct yac_const_basic_grid_data * yac_interp_grid_get_basic_grid_data_src(struct yac_interp_grid *interp_grid)
static int compute_weights_avg_no(double tgt_coords[3], size_t num_vertices, yac_const_coordinate_pointer src_coords, int *src_mask, double *weights, enum yac_cell_type cell_type)
static int compute_weights_avg_yes(double tgt_coords[3], size_t num_vertices, yac_const_coordinate_pointer src_coords, int *src_mask, double *weights, enum yac_cell_type cell_type)
static size_t do_search_avg(struct interp_method *method, struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, struct yac_interp_weights *weights)
static int check_src_field_cell(double *point, size_t field_cell, size_t cell_size, size_t const *vertex_to_cell, size_t const *vertex_to_cell_offsets, yac_const_coordinate_pointer field_coordinates, struct yac_grid_cell *cell_buffer)
static struct interp_method_vtable interp_method_avg_vtable
static void get_point_lon_lat(double point_coord[3], double cell_lon[2], double cell_lat[2], double *point_lon, double *point_lat)
int(* func_compute_weights)(double[3], size_t, yac_const_coordinate_pointer, int *, double *, enum yac_cell_type cell_type)
struct interp_method * yac_interp_method_avg_new(enum yac_interp_avg_weight_type weight_type, int partial_coverage)
static func_compute_weights select_compute_weight_routine(enum yac_interp_avg_weight_type weight_type, int partial_coverage)
static void compute_barycentric_coords(double *barycentric_coords, size_t triangle_indices[3], yac_const_coordinate_pointer grid_coords)
static int compute_weights_dist_no(double tgt_coords[3], size_t num_vertices, yac_const_coordinate_pointer src_coords, int *src_mask, double *weights, enum yac_cell_type cell_type)
static int compute_weights_bary_no(double tgt_coords[3], size_t num_vertices, yac_const_coordinate_pointer src_coords, int *src_mask, double *weights, enum yac_cell_type cell_type)
static int determine_triangle_idx(double point_lon, double point_lat, double cell_lon[2], double cell_lat[2])
#define IS_LON_LAT(x)
static int compute_weights_dist_yes(double tgt_coords[3], size_t num_vertices, yac_const_coordinate_pointer src_coords, int *src_mask, double *weights, enum yac_cell_type cell_type)
static void delete_avg(struct interp_method *method)
static int get_lat_edge_ordering(double const *a, double const *b)
static void get_cell_lon_lat_bounds(yac_const_coordinate_pointer coords, double lon[2], double lat[2], int reorder[4])
static enum yac_cell_type determine_cell_type(enum yac_edge_type const *edge_types, size_t const *edge_indices, size_t num_edges)
static int compute_weights_bary_yes(double tgt_coords[3], size_t num_vertices, yac_const_coordinate_pointer src_coords, int *src_mask, double *weights, enum yac_cell_type cell_type)
static void compute_weights_bary_reg(double tgt_coords[3], yac_const_coordinate_pointer src_coords, double *weights)
yac_interp_avg_weight_type
@ YAC_INTERP_AVG_DIST
@ YAC_INTERP_AVG_ARITHMETIC
@ YAC_INTERP_AVG_BARY
static void compute_weights(struct tgt_point_search_data *tgt_point_data, size_t num_tgt_points, struct edge_interp_data *edge_data, size_t num_edges, struct triangle_interp_data *triangle_data, size_t num_triangles, struct weight_vector_data **weights, size_t **num_weights_per_tgt, size_t *total_num_weights)
void yac_interp_weights_add_wsum(struct yac_interp_weights *weights, struct remote_points *tgts, size_t *num_src_per_tgt, struct remote_point *srcs, double *w)
yac_location
Definition location.h:12
@ YAC_LOC_CORNER
Definition location.h:15
@ YAC_LOC_CELL
Definition location.h:14
#define xrealloc(ptr, size)
Definition ppm_xfuncs.h:67
#define xmalloc(size)
Definition ppm_xfuncs.h:66
struct interp_method_vtable * vtable
func_compute_weights compute_weights
size_t(* do_search)(struct interp_method *method, struct yac_interp_grid *grid, size_t *tgt_points, size_t count, struct yac_interp_weights *weights)
struct remote_point * data
const_size_t_pointer cell_to_vertex_offsets
const_size_t_pointer cell_to_vertex
const const_int_pointer num_vertices_per_cell
const_size_t_pointer cell_to_edge
const_size_t_pointer cell_to_edge_offsets
const_yac_edge_type_pointer edge_type
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 MIN(a, b)
#define MAX(a, b)
void yac_quicksort_index_size_t_size_t(size_t *a, size_t n, size_t *idx)
#define YAC_ASSERT(exp, msg)
Definition yac_assert.h:16
Xt_int yac_int
Definition yac_types.h:15
double const (* yac_const_coordinate_pointer)[3]
Definition yac_types.h:20
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19