numpy implementation of barycentric interpolation from delaunay triangles

BarycentricInterpolation.py

import numpy as np
from scipy.spatial import Delaunay

def calc_dis(x1, y1, x2, y2, x3, y3):
    return np.sqrt(np.power(x1 - x2, 2) + np.power(y1 - y2, 2)) + \
           np.sqrt(np.power(x1 - x3, 2) + np.power(y1 - y3, 2)) + \
           np.sqrt(np.power(x2 - x3, 2) + np.power(y2 - y3, 2))

def calc_triangle_area(x1, y1, x2, y2, x3, y3):
    return 1 / 2 * ((x1 * y2 - x2 * y1) + (x2 * y3 - x3 * y2) + (x3 * y1 - x1 * y3))

def interp(triangles, values, width=1280, height=720, th=100):
    """
    重心插值
    :param triangles: scipy中的Delaunay object，由原始数据点生成
    :param values: 原始数据的value
    :param width: grid 宽度(默认grid粒度为1)
    :param height: grid 高度(默认grid粒度为1)
    :param th: 过滤Delaunay三角的周长阈值 防止边缘出现太大或者太扁的三角
    :return:
    """
    out = np.zeros((width, height))

    # 生成grid
    xi = np.array([[i, j] for i in range(width) for j in range(height)], dtype=np.int)

    # 找到点所属的三角形
    simplex_ids = triangles.find_simplex(xi)
    indices = np.where(simplex_ids != -1)[0]
    xi = xi[indices, :]
    simplex_ids = simplex_ids[indices]

    # 过滤周长
    # 找到simplex对应的三角形三个定点的坐标 以及顶点坐标对应的值
    V = triangles.points[triangles.simplices[simplex_ids, :].reshape(-1), :].reshape(-1, 6)
    values = values[triangles.simplices[simplex_ids, :].reshape(-1)].reshape(-1, 3)
    x1, y1, x2, y2, x3, y3 = V[:, 0], V[:, 1], V[:, 2], V[:, 3], V[:, 4], V[:, 5]
    # 求面积和周长
    areas = calc_triangle_area(x1, y1, x2, y2, x3, y3)
    C = calc_dis(x1, y1, x2, y2, x3, y3)
    # 过滤掉面积和周长太大的simplex 及其对应的value 和原始数据点坐标
    indices = np.where(C < th)[0]
    simplex_ids = simplex_ids[indices]
    xi = xi[indices, :]
    V = V[indices, :]
    values = values[indices, :]
    areas = areas[indices]

    # 进行重心插值算法
    x1, y1, x2, y2, x3, y3 = V[:, 0], V[:, 1], V[:, 2], V[:, 3], V[:, 4], V[:, 5]
    A1 = calc_triangle_area(xi[:, 0], xi[:, 1], x2, y2, x3, y3)
    A2 = calc_triangle_area(x1, y1, xi[:, 0], xi[:, 1], x3, y3)
    A3 = calc_triangle_area(x1, y1, x2, y2, xi[:, 0], xi[:, 1])
    # 插值算法结果
    interp_result = (A1 * values[:, 0] + A2 * values[:, 1] + A3 * values[:, 2]) / areas

    # 生成图像
    out[xi[:, 0], xi[:, 1]] = interp_result
    print(simplex_ids.shape, xi.shape, values.shape, V.shape)
    print(A1.shape, A2.shape, A3.shape, areas.shape)
    print(interp_result.shape)
    return out
    
if __name__ == '__main__':
    # generate data
    points = np.random.rand(20110, 2) * 512
    points[:5000, 1] = np.sqrt(np.maximum(200 ** 2 - (points[:5000, 0] - 256) ** 2, 0)) + 200
    points[5000:10000, 1] = np.sqrt(np.maximum(185 ** 2 - (points[5000:10000, 0] - 256) ** 2, 0)) + 200
    points[10000:15000, 1] = np.sqrt(np.maximum(175 ** 2 - (points[10000:15000, 0] - 256) ** 2, 0)) + 200
    points[15000:20000, 1] = np.sqrt(np.maximum(100 ** 2 - (points[15000:20000, 0] - 256) ** 2, 0)) + 200
    values = np.random.rand(20110) * 255
    
    # generate triangles
    triangles = Delaunay(points)
    
    # interpolation
    out = interp(triangles, values, width=512, height=512, th=40)
    print(out.shape)
    
    import matplotlib.pyplot as plt
    origin = np.zeros_like(out)
    origin[np.int_(points[:, 0]), np.int_(points[:, 1])] = values
    plt.subplot(1, 2, 1)
    plt.imshow(origin.T, cmap=plt.get_cmap('hot'))
    plt.title('origin')
    plt.subplot(1, 2, 2)
    plt.imshow(out.T, cmap=plt.get_cmap('hot'))
    plt.title('interpolation')
    plt.show()