Bilinear Remapping

[aaronzedwick]

trafficstars

Proposed new feature or change:

Implement bilinear remapping in UXarray. Good paper here by David Marsico and Paul Ullrich.

[aaronzedwick]

function bilinear_remap(source_grid, destination_grid, source_grid_data):
    values = []
    dual = source_grid.get_dual()
    if not source_grid.is_partial():
        # This tree will be used to find the polygon containing the point, as the first polygon that will be searched is one whose face center is nearest
        tree = source_grid.get_ball_tree(coordinates="nodes")

        for point in destination_grid:
            # Find a subset of polygons that contain the point
            polygons_subset = find_polygons_subset(tree, point)

            # Search the subset to find which one contains the point
            for polygon in polygons_subset:
                if point_in_polygon(polygon, point):
                    value = calculate_bilinear_weights(polygon, point)
                    values.append(value)
                    break

    else:  

        tree = source_grid.get_ball_tree(coordinates="nodes")

        for point in destination_grid:
            found = False
            dual_polygons_subset = find_polygons_subset(tree, point)
            
            for polygon in dual_polygons_subset:
                if point_in_polygon(polygon, point):
                    value = calculate_bilinear_weights(polygon, point)
                    values.append(value)
                    found = True
                    break

            # If it isn't found in the dual mesh, search the primal mesh for the point
            if not found:
                primal_tree = source_grid.get_ball_tree(coordinates="face centers")
                primal_polygons_subset = find_polygons_subset(primal_tree, point)
                
                for polygon in primal_polygons_subset:
                    # If found in the primal mesh, assign the nearest neighbor as the value
                    if point_in_polygon(polygon, point):
                        nearest_index = tree.query(point, k=1, return_distance=False)
                        values.append(source_grid_data[nearest_index])
                        found = True
                        break
                        
                # If it is not found in the primal either, assign 0 as the value
                if not found:
                    values.append(0.0)

    return values

function calculate_bilinear_weights(polygon, point):
    if polygon.nodes == 3:
        # Find the weights using barycentric coordinates
        return calculate_barycentric_coordinates(polygon, point)
    elif polygon.nodes == 4:
        # Splitting the quad into two triangles and using barycentric coordinates on the triangle containing the point might be better
        return calculate_weights_with_newton_quadrilaterals(polygon, point)
    else:
        # For polygons that aren't triangles or quads
        sub_triangle = find_sub_triangle(polygon, point)
        return calculate_barycentric_coordinates(sub_triangle, point)

@hongyuchen1030 This is my general idea for the remapping logic.

[aaronzedwick]

@philipc2 Depending on how we find the subset of polygons we might not use the tree structures at all, right? Except I suppose for the situation where the point is not in a polygon in the dual mesh, but one in the primal. Does our subsetting use tree structures at all?

[hongyuchen1030]

@philipc2 Depending on how we find the subset of polygons we might not use the tree structures at all, right? Except I suppose for the situation where the point is not in a polygon in the dual mesh, but one in the primal. Does our subsetting use tree structures at all?

I am curious for this PR:https://github.com/UXARRAY/uxarray/pull/938

If this is what I think it is doing, assigning each side of an edge a sign to see which side is within the given face, then we're basically doing a mesh-arrangement. And if we are doing mesh-arrangement, we don't need to do a brute force test about if a point is inside a given face

[aaronzedwick]

@hongyuchen1030 the algorithm will end up using ball tree, so it will be constructed before the call to point in polygon.