从 scipy.spatial.Delaunay 中过滤单纯形

Question

简短版本：

是否可以使用现有对象的三角形子集（2D 数据）创建新的 scipy.spatial.Delaunay 对象？

目标是在过滤掉单纯形的新对象上使用 find_simplex 方法。

相似但不完全相同

matplotlib轮廓/ **凹**非网格数据的轮廓

如何处理在 matplotlib 中使用三角测量时在几何图形边缘之间形成的（不需要的）三角形

长版本：

我正在查看我重新网格化的经纬度数据scipy.interpolate.griddata 就像下面的伪代码一样：

import numpy as np
from scipy.interpolate import griddata
from scipy.spatial import Delaunay
from scipy.interpolate.interpnd import _ndim_coords_from_arrays

#lat shape (a,b): 2D array of latitude values
#lon shape (a,b): 2D array of longitude values
#x shape (a,b): 2D array of variable of interest at lat and lon

# lat-lon data
nonan = ~np.isnan(lat)
flat_lat = lat[nonan]
flat_lon = lon[nonan]
flat_x = x[nonan]

# regular lat-lon grid for regridding
lon_ar = np.arange(loni,lonf,resolution)
lat_ar = np.arange(lati,latf,resolution)
lon_grid, lat_grid = np.meshgrid(lon_ar,lat_ar)

# regrid
x_grid = griddata((flat_lon,flat_lat),flat_x,(lon_grid,lat_grid), method='nearest')

# filter out extrapolated values
cloud_points = _ndim_coords_from_arrays((flat_lon,flat_lat))
regrid_points = _ndim_coords_from_arrays((lon_grid.ravel(),lat_grid.ravel()))
tri = Delaunay(cloud_points)
outside_hull = tri.find_simplex(regrid_points) < 0
x_grid[outside_hull.reshape(x_grid.shape)] = np.nan

# filter out large triangles ??
# it would be easy if I could "subset" tri into a new scipy.spatial.Delaunay object
# new_tri = ??
# outside_hull = new_tri.find_simplex(regrid_points) < 0

问题是凸包的质量低（非常大，在下面的示例中以蓝色显示）三角形，我想过滤掉它们，因为它们不代表数据出色地。我知道如何在输入点中过滤掉它们，但不知道如何在重新网格化的输出中过滤掉它们。这是过滤器函数：

def filter_large_triangles(
    points: np.ndarray, tri: Optional[Delaunay] = None, coeff: float = 2.0
):
    """
    Filter out triangles that have an edge > coeff * median(edge)
    Inputs:
        tri: scipy.spatial.Delaunay object
        coeff: triangles with an edge > coeff * median(edge) will be filtered out
    Outputs:
        valid_slice: boolean array that selects "normal" triangles
    """
    if tri is None:
        tri = Delaunay(points)

    edge_lengths = np.zeros(tri.vertices.shape)
    seen = {}
    # loop over triangles
    for i, vertex in enumerate(tri.vertices):
        # loop over edges
        for j in range(3):
            id0 = vertex[j]
            id1 = vertex[(j + 1) % 3]

            # avoid calculating twice for non-border edges
            if (id0,id1) in seen:
                edge_lengths[i, j] = seen[(id0,id1)]
            else:    
                edge_lengths[i, j] = np.linalg.norm(points[id1] - points[id0])

                seen[(id0,id1)] = edge_lengths[i, j]

    median_edge = np.median(edge_lengths.flatten())

    valid_slice = np.all(edge_lengths < coeff * median_edge, axis=1)

    return valid_slice

坏三角形如下图蓝色所示：

import matplotlib.pyplot as plt
no_large_triangles = filter_large_triangles(cloud_points,tri)
fig,ax = plt.subplot()
ax.triplot(points[:,0],points[:,1],tri.simplices,c='blue')
ax.triplot(points[:,0],points[:,1],tri.simplices[no_large_triangles],c='green')
plt.show()

是否可以创建一个新的 scipy.spatial。仅具有 no_large_triangles 单形的 Delaunay 对象？目标是在该新对象上使用 find_simplex 方法来轻松过滤掉点。

作为替代方案，我如何找到位于蓝色三角形内的 regrid_points 中的点的索引？ (tri.simplices[~no_large_triangles])

Answer 1

因此，可以修改 Delaunay 对象，以便在单纯形的子集上使用 find_simplex，但似乎只能使用暴力算法。

# filter out extrapolated values
cloud_points = _ndim_coords_from_arrays((flat_lon,flat_lat))
regrid_points = _ndim_coords_from_arrays((lon_grid.ravel(),lat_grid.ravel()))
tri = Delaunay(cloud_points)
outside_hull = tri.find_simplex(regrid_points) < 0

# filter out large triangles
large_triangles = ~filter_large_triangles(cloud_points,tri)
large_triangle_ids = np.where(large_triangles)[0]
subset_tri = tri # this doesn't preserve tri, effectively just a renaming
# the _find_simplex_bruteforce method only needs the simplices and neighbors
subset_tri.nsimplex = large_triangle_ids.size
subset_tri.simplices = tri.simplices[large_triangles]
subset_tri.neighbors = tri.neighbors[large_triangles]

# update neighbors
for i,triangle in enumerate(subset_tri.neighbors):
    for j,neighbor_id in enumerate(triangle):
        if neighbor_id in large_triangle_ids:
            # reindex the neighbors to match the size of the subset
            subset_tri.neighbors[i,j] = np.where(large_triangle_ids==neighbor_id)[0]
        elif neighbor_id>=0 and (neighbor_id not in large_triangle_ids):
            # that neighbor was a "normal" triangle that should not exist in the subset
            subset_tri.neighbors[i,j] = -1
inside_large_triangles = subset_tri.find_simplex(regrid_points,bruteforce=True) >= 0

invalid_slice = np.logical_or(outside_hull,inside_large_triangles)
x_grid[invalid_slice.reshape(x_grid.shape)] = np.nan

显示新的 Delaunay 对象仅具有大三角形的子集

import matplotlib.pyplot as plt
fig,ax = plt.subplot()
ax.triplot(cloud_points[:,0],cloud_points[:,1],subset_tri.simplices,color='red')
plt.show()

在过滤大三角形之前使用 pcolormesh 绘制 x_grid（放大上面的蓝色圆圈）：

过滤后：