updated method for finding barycentric weights

Author: aaronzedwick

benchmarks/mpas_ocean.py

@@ -138,6 +138,9 @@ def time_nearest_neighbor_remapping(self):
     def time_inverse_distance_weighted_remapping(self):
         self.uxds_480["bottomDepth"].remap.inverse_distance_weighted(self.uxds_120.uxgrid)
 
+    def time_bilinear_remapping(self):
+        self.uxds_480["bottomDepth"].remap.bilinear(self.uxds_120.uxgrid)
+
 
 class HoleEdgeIndices(DatasetBenchmark):
     def time_construct_hole_edge_indices(self, resolution):

uxarray/grid/geometry.py

@@ -24,6 +24,7 @@
 from uxarray.grid.intersections import (
     gca_gca_intersection,
 )
+from uxarray.grid.utils import _get_cartesian_face_edge_nodes
 from uxarray.utils.computing import allclose
 
 POLE_POINTS_XYZ = {
@@ -1251,7 +1252,7 @@ def barycentric_coordinates_cartesian(polygon, point):
     Parameters
     ----------
     polygon: np.array
-        Cartesian coordinates of the polygons' nodes
+        Cartesian coordinates of the polygon nodes
     point: np.array
         Cartesian coordinates of the point.
 
@@ -1262,7 +1263,7 @@ def barycentric_coordinates_cartesian(polygon, point):
     >>> point_xyz = np.array([0.0, 0.0, 1.0], dtype=np.float64)
 
     Define a cartesian polygon:
-    >>> point_xyz = np.array(
+    >>> polygon = np.array(
     ...     [[0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 0.0, -1.0]], dtype=np.float64
     ... )
 
@@ -1272,11 +1273,13 @@ def barycentric_coordinates_cartesian(polygon, point):
     Returns
     -------
     weights: np.ndarray
-        Barycentric coordinates corresponding to each vertex in the polygon.
+        Barycentric coordinates corresponding to each vertex in the triangle.
+    nodes: np.ndarray
+        Indices of which nodes of the polygon form the triangle containing the point
     """
+
     n = len(polygon)
-    weights = np.zeros(n)
-    total_area = 0.0
+    weights = np.zeros(3)
 
     # Use the first vertex as a reference point to form triangles
     p0 = polygon[0]
@@ -1285,41 +1288,64 @@ def barycentric_coordinates_cartesian(polygon, point):
         # Get the two other points of the triangle
         p1, p2 = polygon[i], polygon[i + 1]
 
-        # Find the inside portion of the equation: p - p0 =(p1 - p0)b + (p2 - p0)c
-        v1, v2 = p1 - p0, p2 - p0
+        # Get the faces in terms of their edges for the point_in_face check
+        face_edge = _get_cartesian_face_edge_nodes(
+            face_idx=0,
+            face_node_connectivity=np.array([[0, 1, 2]]),
+            n_edges_per_face=np.array([3]),
+            node_x=np.array([p0[0], p1[0], p2[0]], dtype=np.float64),
+            node_y=np.array([p0[1], p1[1], p2[1]], dtype=np.float64),
+            node_z=np.array([p0[2], p1[2], p2[2]], dtype=np.float64),
+        )
+
+        # Check to see if the point lies within the current triangle
+        contains_point = point_in_face(
+            face_edge,
+            point,
+            inclusive=True,
+        )
+
+        # If the point is in the current triangle, get the weights for that triangle
+        if contains_point:
+            # Find the inside portion of the equation: p - p0 =(p1 - p0)b + (p2 - p0)c
+            v1, v2 = p1 - p0, p2 - p0
+
+            # Solve the left side of the equation
+            vp = point - p0
 
-        # Solve the left side of the equation
-        vp = point - p0
+            # Find the largest coordinate, we will ignore this one and use the other two
+            abs_det = np.abs(np.cross(v1, v2))
+            max_coord = np.argmax(abs_det)
 
-        # Find the largest coordinate, we will ignore this one and use the other two
-        abs_det = np.abs(np.cross(v1, v2))
-        max_coord = np.argmax(abs_det)
+            # Select the two coordinates to find corresponding weights for
+            idx = [j for j in range(3) if j != max_coord]
+            array = np.array([v1[idx[0]], v1[idx[1]]])
+            array2 = np.array([v2[idx[0]], v2[idx[1]]])
 
-        # Select the two coordinates to find corresponding weights for
-        idx = [j for j in range(3) if j != max_coord]
-        array = np.array([v1[idx[0]], v1[idx[1]]])
-        array2 = np.array([v2[idx[0]], v2[idx[1]]])
+            A = np.column_stack((array, array2))  # 2x2 matrix
+            b = np.array([vp[idx[0]], vp[idx[1]]])  # Right-hand side vector
 
-        A = np.column_stack((array, array2))  # 2x2 matrix
-        b = np.array([vp[idx[0]], vp[idx[1]]])  # Right-hand side vector
+            # Solve for weights b, c
+            b_c = np.linalg.solve(A, b)
 
-        # Solve for weights b, c
-        b_c = np.linalg.solve(A, b)
+            # Solve for a
+            a = 1 - b_c[0] - b_c[1]
 
-        # Solve for a
-        a = 1 - b_c[0] - b_c[1]
+            # Compute triangle area
+            triangle_area = 0.5 * np.linalg.norm(np.cross(v1, v2))
 
-        # Compute triangle area
-        triangle_area = 0.5 * np.linalg.norm(np.cross(v1, v2))
-        total_area += triangle_area
+            # Calculate weights
+            weights[0] += a * triangle_area
+            weights[1] += b_c[0] * triangle_area
+            weights[2] += b_c[1] * triangle_area
 
-        # Accumulate weights
-        weights[0] += a * triangle_area
-        weights[i] += b_c[0] * triangle_area
-        weights[i + 1] += b_c[1] * triangle_area
+            # Normalize weights so they sum to 1
+            if triangle_area > 0:
+                weights /= triangle_area
 
-    # Normalize weights so they sum to 1
-    if total_area > 0:
-        weights /= total_area
+            # The nodes of the triangle, so it can be known which values to apply the weights to
+            nodes = np.array([0, i, i + 1])
+            return weights, nodes
 
-    return weights
+    # If the point doesn't reside in the polygon, raise an error
+    raise ValueError("Point does not reside in polygon")

uxarray/remap/bilinear.py

@@ -9,10 +9,11 @@
     from uxarray.core.dataset import UxDataset
 
 import numpy as np
+from numba import njit
 
 import uxarray.core.dataarray
 import uxarray.core.dataset
-from uxarray.constants import INT_FILL_VALUE
+from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
 from uxarray.grid import Grid
 
 
@@ -239,7 +240,33 @@ def _get_values(point_xyz, point_lonlat, dual, source_data, data_size, source_gr
                 # Create the data array on the polygon
                 data[ind] = source_data[node]
 
-            weights = barycentric_coordinates_cartesian(polygon, point_xyz[i])
-            values[i] = np.sum(weights * data, axis=-1)
+            # Check to see if the point lies on a node of the polygon
+            matching_indices = _find_matching_node_index(
+                polygon, point_xyz[i], tolerance=ERROR_TOLERANCE
+            )
+
+            # If the point lies on a node, assign it a weight of 1 and return the node index
+            if matching_indices[0] != -1:
+                weights, nodes = np.array([1]), matching_indices
+            else:
+                weights, nodes = barycentric_coordinates_cartesian(
+                    polygon, point_xyz[i]
+                )
+
+            values[i] = np.sum(weights * data[nodes], axis=-1)
 
     return values
+
+
+@njit(cache=True)
+def _find_matching_node_index(nodes, point, tolerance=ERROR_TOLERANCE):
+    for i in range(nodes.shape[0]):
+        match = True
+        for j in range(3):  # Compare each coordinate
+            diff = abs(nodes[i, j] - point[j])
+            if diff > tolerance + tolerance * abs(point[j]):
+                match = False
+                break
+        if match:
+            return np.array([i])  # Return first matching index
+    return np.array([-1])  # Not found