手把手教你用傅立叶变换画可达鸭

手把手教你用傅立叶变换画可达鸭

傅立叶变换能干啥?

多图预警

上面的心就是通过傅里叶变换画出来的!

接下来介绍方法。

离散傅立叶变换

假设我有N个复数,表示为向量x。离散傅立叶变换可以把x转换成N个新数,表示为向量X,满足如下等式:

X_k = \sum_{n=0}^{N-1} x_n\cdot e^{-\frac {i 2\pi}{N}kn} = \sum_{n=0}^{N-1} x_n\cdot [\cos(2 \pi k n / N) - i\cdot \sin(2 \pi k n / N)]

X又可以反过来变回x

x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k\cdot e^{i 2 \pi k n / N}

公式都是维基抄的。数学就到此为止。

下面给具体的例子。

我们来看这8个点组成的粗糙的心形

如果我把这些点做傅里叶变换,可以得到这样一些点

看起来很随机没啥规律。把他们的坐标看成向量加起来看看

In [42]: X
Out[42]: 
array([24.        +22.j        , -1.41421356 -1.41421356j,
       -2.         +0.j        ,  0.24264069 -0.24264069j,
        0.         +2.j        ,  1.41421356 +1.41421356j,
        2.         +0.j        , -8.24264069 +8.24264069j])

In [43]: x
Out[43]: array([2.+4.j, 3.+3.j, 4.+4.j, 5.+3.j, 4.+2.j, 3.+1.j, 2.+2.j, 1.+3.j])

In [44]: fft(x)
Out[44]: 
array([24.        +22.j        , -1.41421356 -1.41421356j,
       -2.         +0.j        ,  0.24264069 -0.24264069j,
        0.         +2.j        ,  1.41421356 +1.41421356j,
        2.         +0.j        , -8.24264069 +8.24264069j])

In [45]: sum(fft(x))
Out[45]: (16+32j)

看起来蛮整的。是巧合么?

其实它正好是第一个点的8倍。

其实就对应了上面公式

x_0 = \frac{1}{8}(X_0\cdot e^0 + X_1\cdot e^0 + X_2\cdot e^0 + X_3\cdot e^0 + X_4\cdot e^0 + X_5\cdot e^0 + X_6\cdot e^0 + X_7\cdot e^0)

x_1 呢?

x_1 = \frac{1}{8}(X_0\cdot e^{0} + X_1\cdot e^{i2\pi/8} + X_2\cdot e^{{i4\pi/8}} + X_3\cdot e^{{i6\pi/8}} + X_4\cdot e^{{i8\pi/8}} + X_5\cdot e^{i10\pi/8} + X_6\cdot e^{i12\pi/8} + X_7\cdot e^{i14\pi/8})

验证一下

In [51]: sum(fft(x) * np.exp(np.array(range(8)) * (1j * 2 * np.pi / 8)) / 8)
Out[51]: (3+3j)

确实如此。

其实 x_1 的等式相当于让 X_0 不动, X1 逆时针转45度, X2 逆时针转90度,以此类推。X的下表越大,转的幅度越大。

相似的 x_2 ,相当于在 x_1 的基础上, X_1 再转45度, X_2 再转90度。。。

也就是说,x下标的增长,对应着X向量的旋转。

如果x的下标每秒增长1,X里面的向量就会像下面一样转起来旋转

第一个圆一动不动,决定了图像的中心。每个元的起始相位和半径不同,这都是由X决定的。下面四个圆转的太快了以至于看起来在反着转。

这些圆如果首尾相接,第二个元的圆心定在第一个元的向量末端,第三个圆的圆心定在第二个元的向量末端,那么最后一个圆的向量指向的位置就是所有向量的和,也就是心形上的那8个点。


理论上来说,只要找出图片的轮廓,我们就能通过N个首尾相接的圆来逼近这个轮廓,N的大小等于轮廓上采样的点数。

剩下的问题就是怎么提取轮廓了。

先准备好图


用Pillow打开图片

from PIL import Image

img = Image.open('psyduck.png')

matplotlib.pyplot.contour要求输入像素为单个字节,所以我们先转换成灰度像素

img = img.convert('L')

再用matplotlib提取轮廓

import pylab
fig, ax = pylab.subplots()
ax.contour(img, origin='image', levels=[100])
pylab.show()

轮廓只是一些像素点围成的多边形。接下来我们要把轮廓转换成画笔的路径。

不过这个问题比看起来要困难不少。

  • 孤岛

画笔从一个多边形跳转到另一个多边形,一定会留下多余的线条,影响效果。但有的跳转是必须的,比如可达鸭的眼睛和鼻子都是孤岛。

  • 样本点的顺序

我们虽然知道样本点要从轮廓上取,但是每个轮廓从哪里开始?是顺时针还是逆时针?是画完一个多边形再画第二个,还是尽可能“一笔画”?

  • 圆圈太多

傅里叶变换需要n个圆圈来拟合n个样本点。如果样本点取很多,虽然轮廓细节可以勾勒清晰,但是动画渲染会很慢。如何用尽可能少的圆圈勾勒出逼真的轮廓?

孤岛问题

画笔画出来的轨迹是连续的,也就是说画笔将所有样本点连成一块连通分量。要减少孤岛的违和感,需要用尽可能短的边将其连入联通分量中。

而最小生成树的Kruskal算法也用了类似的思想。

用相似的贪心法,我们不难把孤岛连入整个轮廓。

样本点的顺序

理想情况下,不带重复的一笔画是最佳的画笔轨迹。一笔画问题就是寻找欧拉路径。然而,只有2或0个点的度数为基数的连通图存在欧拉路径。连通图尚且不论,度数为3的点在可达鸭的轮廓中就有好多个,不可能存在欧拉路径。也就是说一笔画不适用于一般的图片。

退而求其次,如果我们允许画笔走回头路,同一条边画多次,那么是不是就有办法了呢?答案是肯定的。只要在连通图中做一次深度优先遍历(DFS),遍历的顺序就可以看作是画笔的轨迹了。

然而怎么把轮廓变成连通图呢?

刚刚的最小生成树给了我们一些新的启发。

如果我们把所有的轮廓上的每个像素点看作我们的样本点,那么它组成的最小生成树几乎肯定包含了大多数轮廓上的边,同时又把孤岛包含了进去,一箭双雕。

现在理论上我们已经能够用傅里叶变换画出可达鸭了。

最小生成树的效率

我们需要在一个完全图里找最小生成树,这恰巧是Kruskal算法最不擅长的情况,复杂度基本上达到 O(n^2logn) 。可达鸭的轮廓有上万个点,用python算一次最小生成树可能需要几分钟。

好在二维平面欧氏距离的最小生成树有更好的算法,只需要 O(n logn) 复杂度,我的电脑一秒以内就能算出来。大体的办法是先对点集做一个叫Delaunay的三角剖分,可以证明二位平面欧氏距离的最小生成树在这个三角剖分中。Delaunay需要 O(nlogn) ,三角剖分的边是 O(n) ,只对这些边跑Kruskal就只用 O(nlogn) 了。

Delaunay在scipy里有现成实现。

Delaunay 三角剖分后的可达鸭。最小生成树就藏在里面。
import numpy as np
from scipy.spatial import Delaunay
import heapq

def dis(samples, i, j):
        return np.linalg.norm(samples[i] - samples[j])

def mst(samples):
        n = len(samples)
        tri = Delaunay(samples)
        g = [[] for i in range(n)]   # 邻接表
        edges = {}
        nodes = set()
        # 收集Delaunay三角剖分中的边
        for simplex in tri.simplices:
            nodes |= set(simplex)
            for k in range(3):
                i, j = simplex[k - 1], simplex[k]
                edge = min(i, j), max(i, j)
                if edge not in edges:
                    edges[edge] = dis(samples, i, j)
        pq = [(d, i, j) for ((i, j), d) in edges.items()]
        heapq.heapify(pq)
        p = list(range(n))

        # 并查集
        def union(i, j):
            p[find(i)] = find(j)

        def find(i):
            if p[i] == i:
                return i
            p[i] = find(p[i])
            return p[i]

        # Kruskal
        cc = len(nodes)
        while cc > 1:
            d, i, j = heapq.heappop(pq)
            if find(i) != find(j):
                union(i, j)
                g[i].append((j, d))
                g[j].append((i, d))
                cc -= 1
        return g
左边是样本点,右边是最小生成树。由于样本很多,已经看不出树的模样。
样本减少后可以看出生成树的模样。失去灵魂的可达鸭。

遍历的起点和终点

最小生成树没有根,任何一个点都可以开始遍历。无论从哪里开始遍历,因为回溯的关系,走完的路径都是所有边的和的两倍。

有些回头路走的看起来就很傻:


有没有办法解决呢?

显然,当所有的边都画过至少一次之后,最后一次回溯是没有必要的。最后一次回溯的长短取决于遍历的起始点和子节点遍历的顺序。

如果我们能找出树中最长的路径,从其一端出发,遍历到另一端的时候正好所有边都画过一边,就可以省掉这个最长路了。

树的最长路径可以用DFS找到。只要调整一下子结点的顺序,就可以保证走到终点时所有子树都已遍历。

DFS是递归算法, 调用前需要调整python的recursion limit防止爆栈。

def find_farthest_leaf_pair(g):
    def dfs(i, parent):
        """
        Return
            - farthest leaf id in thissubtree and distance to root i
            - farthest leave pair in this subtree and distance between them
        """
        farthest_leaf = i
        farthest_leaf_dis = 0
        farthest_leaf_pair = None
        farthest_leaf_pair_dis = -1
        leave_dis = []
        for j, _ in g[i]:
            if j == parent:
                continue
            l, ld, pair, pair_dis = dfs(j, i)
            leave_dis.append((ld + 1, l))
            if ld + 1 > farthest_leaf_dis:
                farthest_leaf_dis = ld + 1
                farthest_leaf = l
            if farthest_leaf_pair_dis < pair_dis:
                farthest_leaf_pair = pair
                farthest_leaf_pair_dis = pair_dis
        if len(leave_dis) >= 2:
            (d1, l1), (d2, l2) = sorted(leave_dis)[-2:]
            if d1 + d2 > farthest_leaf_pair_dis:
                farthest_leaf_pair_dis = d1 + d2
                farthest_leaf_pair = l1, l2
        return farthest_leaf, farthest_leaf_dis, farthest_leaf_pair, farthest_leaf_pair_dis

    for i in range(len(g)):
        if len(g[i]):
            l, ld, pair, pair_dis = dfs(i, -1)
            if len(g[i]) == 1 and ld > pair_dis:
                # 根是叶子
                return i, l
            return pair


def rearange_children_order(g, st, ed):
    vis = set()

    def dfs(i):
        vis.add(i)
        if i == ed:
            return True
        for j in range(len(g[i])):
            if g[i][j][0] not in vis:
                if dfs(g[i][j][0]):
                    g[i][j], g[i][-1] = g[i][-1], g[i][j]
                    return True
        return False
    dfs(st)
    return st, ed

def generate_path(g, st, ed):
    res = []
    vis = set()

    def dfs(i):
        vis.add(i)
        res.append(samples[i])
        if i == ed:
            return True
        leaf = True
        for j, _ in g[i]:
            if j not in vis:
                leaf = False
                if dfs(j):
                    return True
        if not leaf:
            # 叶子只加一次
            res.append(samples[i])
        return False
    dfs(st)
    return res

圆圈太多

可达鸭的轮廓有上万个像素点组成,傅里叶变换后转换成上万个圆圈,效果并不理想。

要减少点的数量,最简单的方法就是取随机样本。

但是随机样本无法捕捉细节信息。当样本数量达到219的时候最小生成树画出的图已经有些难以辨认了。

随机样本遇到了瓶颈,我们需要另想办法。

注意到可达鸭的轮廓有的地方平滑,细节少;有的地方曲折,细节多。

几乎是直线的地方信息含量少,只需要少量样本就能描绘,而细节丰富的地方需要更多的样本。有什么办法能按照这样的规律来调整我们的路径呢?

Ramer–Douglas–Peucker算法就是为此就是为此设计。它能够简化几乎排一条直线上的样本点,只保留首尾;曲折较多的样本点则被更完整的保留。

【更新 07/16/2019】我原先使用了这个python的实现。

fhirschmann/rdpgithub.com图标

由于python的效率较低,我自己用cython重写了一个版本。

biran0079/crdpgithub.com图标

这样大多数图片不需要过度降取样,也能快速的计算路径。

def rdp_downsample(path, epsilon):
    before_n = len(path)
    path = rdp(path, epsilon=epsilon)
    after_n = len(path)
    print(f'samples on path path after rdp: {before_n} -> {after_n}')
    return path

Ramer–Douglas–Peucker可以将成千上万个点组成的路径减少1-2个数量级,同时不丢失太多轮廓细节。

下图中可达鸭的路径只用可182个样本点,效果远好于上面241个随机的样本点组成的路径。

有了样本路径之后,傅里叶变换就能帮我们得到圆圈的信息。

用scipy中快速傅里叶变换的实现,只需一行代码

from scipy.fftpack import fft
X = fft(path[..., 0] + path[..., 1] * 1j)

然后推算出每个时间点每个圆圈的位置,最后用matplotlib画出动画

frames = len(X)  # number of animation frames
n = frames       # number of ciecles
coef = np.arange(frames)
# pos[i][t] is the position of ith circle's center at t
# pos[0][t] = 0 for all t
# pos[i][t] = pos[i - 1][t] + factor[i][t]
# where factor[i][t] = x[i] * exp(2*pi*t*sqrt(-1)*i/n)
factor = []
for i in range(len(coef)):
  factor.append(x[i] * np.exp(2 * np.pi * 1j * coef[i] * np.arange(frames) / frames) / frames)
pos = [[0 for i in range(frames)]]
for i in range(len(coef)):
  pos.append(pos[-1] + factor[i])

fig, ax = pylab.subplots()
lines = []
circles = []
def _init_lines():
    lines += [ax.plot([], [], alpha=0.3)[0] for i in range(n)]

def _init_circles():
    for i in range(n):
        radius = np.abs(pos[i+1][0] - pos[i][0])
        circle = plt.Circle((-1, -1), radius, alpha=0.3, fill=False)
        circles.append(circle)
        ax.add_artist(circle)

path_x = []
path_y = []

def init_frame():
    ax.clear()
    ax.set(xlim=[0, 1], ylim=[0, 1])
    ax.axis('off')
    _init_lines(ax)
    _init_circles(ax)
    global path
    path = ax.plot([], [])[0]

def _render_lines(k):
    for i in range(n):
        p, q = pos[i][k], pos[i+1][k]
        lines[i].set_data([p.real, q.real], [p.imag, q.imag])

def _render_circles(k):
    for i in range(n):
        p = pos[i][k]
        circles[i].center = (p.real, p.imag)

def render_frame(k):
    _render_lines(k)
    _render_circles(k)
    p = pos[n][k]
    path_x.append(p.real)
    path_y.append(p.imag)
    path.set_data(path_x, path_y)

animation = FuncAnimation(fig,
                          render_frame,
                          frames=range(frames),
                          interval=100,
                          repeat_delay=1000,
                          init_func=init_frame,
                          repeat=True)
pylab.show()

代码和使用方法都在这里。里面有python的notebook,可以在binder上运行。

还有一些细节的东西本文里没有提到,比如用savgol_filter来平滑路径,大家可以到代码里去探索。


biran0079/circlegithub.com图标

其他的作品:

https://www.zhihu.com/video/1135153125211041792

https://www.bilibili.com/video/av59229141


附上网页版(JS较多,第一次打开较慢):circle.fancy-coder.com/

biran0079/circle-jsgithub.com图标

编辑于 2019-07-30

文章被以下专栏收录