算法►渐变算法

不规则平面区域的渐变过渡算法

[toc]

一、问题描述

这是个两年前的工作了，想着如果不整理一下，哪天电脑坏了，估计也就跟着没了吧，这样的事情也不少发生。

当时在设计一个新算法，其中的一小步就是需要对不规则平面区域进行插值。不规则平面区域由封闭曲线构成，下图给出了一个例子：

当然，上面这张图片中的云朵形状还是过于复杂了，但是暂时没找到合适的示意图，就凑合吧。实际需求没有这么复杂，形状整体上更偏向于一个凸集。

这里要解决的问题是，对于由封闭曲线构成的不规则区域，如何设计一个类似双线性插值的算法，使得区域内部的点能够根据区域轮廓曲线的属性值，从上至下或者从左至右进行平滑过渡呢？

二、问题分析与算法实现

2.1 四边形区域的插值

如果不规则区域是一个四边形，只需要进行双线性插值即可满足需求。现在回去看当时写的代码，里面有一段magic code，完全不知道是在干嘛的：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
def solve(vertices, coordinates, axis=0):
    x0, y0 = vertices[0]
    x1, y1 = vertices[1]
    x2, y2 = vertices[2]
    x3, y3 = vertices[3]
    y, x = coordinates

    if axis == 0:
        a = y - y3
        b = y2 - y3
        c = x - x3
        d = x2 - x3
        e = y2 - y3 - y1 + y0
        f = y3 - y0
        g = x2 - x3 - x1 + x0
        h = x3 - x0
    else:
        a = y - y1
        b = y2 - y1
        c = x - x1
        d = x2 - x1
        e = y3 - y0 - y2 + y1
        f = y0 - y1
        g = x3 - x0 - x2 + x1
        h = x0 - x1

    A = d * e - b * g
    B = a * g - b * h - c * e + d * f
    C = a * h - c * f

    r = (-B + np.sqrt(B ** 2 - 4 * A * C)) / (2 * A)
    r[r < 0] = 0

    return r

所以说及时做笔记是相当重要的。后来细想一下，应该是这么推导得来的。设四边形的的四个端点分别为，，和，从左下角的端点A开始按照逆时针顺序进行排列：

将端点坐标记为，，，。

线段DC的参数方程为：

同理，线段AB的参数方程为：

设是线段DC上的一点，是线段AB上的点，对应的参数均为。则线段HI的参数方程为：

设区域内一点为在线段上，其坐标为，对应的参数为，则：

展开得方程组：

两个方程，两个未知数，有手就能解：

将代入第二等式，得：

等式两边同乘：

合并同类项：

可以看出这是个一元二次方程，根据求根公式容易求得的解。

最终的效果为：

具体的实现代码可参考：https://github.com/mo-vic/xy-interpolation/blob/main/quadrangle.py

2.2 不规则区域的插值

如果不规则区域是由封闭曲线构成的，可以仿照上面的思路，将区域的轮廓进行划分，构造一个曲边四边形。具体的划分方式是，将轮廓中的左上角，左下角，右下角，右上角四个点作为分界点，将封闭曲线划分为四段，对于每一段轮廓曲线，使用贝塞尔曲线进行拟合：

以从左至右的过渡方式为例，当点靠近蓝色曲线时，值应该接近于0，当点靠近红色曲线时，值应该接近于1。红蓝曲线的形状均由它们的控制点决定。因此，对于区域内点的取值，可以按比例混合红蓝曲线的控制点，并从橙色和绿色曲线中以同样的参数各取一个点分别作为起始点和结束点，构造一条新的贝塞尔曲线，令新构造的贝塞尔曲线上的点的取值均为即可。

可以看到，这样的插值方式有时候并不完美，新构造的贝塞尔曲线可能会越过区域边界，但是暂时没想到更好的方法了。

最终的效果为：

从上图可以看到，在插值的时候，区域内部有些点没有被采样到，可以通过提高采样分辨率来提高覆盖率。但是有另外一种方式，就是使用已有数据构造一棵KDTree，在查询的时候，通过KDTree取最近邻值即可。

算法实际应用时的效果如下：

代码可以参考：https://github.com/mo-vic/xy-interpolation/blob/main/curve.py

现在回想起来，其实可以试一试蒙特卡洛方法的，但是今年是2023年了，不是2020年了。