几何建模与处理-实现平面点集CVT的Lloyd算法

老铁，点个在看呗 欢迎关注公众号:sumsmile /专注图形渲染、Android开发

内容主要参考中科大刘利刚老师的课程《几何建模与处理基础》

这篇文章的内容是独立，不需要多少数学知识，而且非常有趣。 最后实现的效果：

三角化

给出一些离散的点构建出几何形状，需要进行三角化处理，也称为"三角剖分"。

比如下图中离散的点，依次连接，形成三角形剖分

直观感觉左边的三角剖分更"合理"些，三角网格整体看起来更均匀。

网格质量是另一个话题，如何衡量三角剖分的优劣，有多个维度，暂不展开。

关于三角形剖分，必须提到两位大数学家

• Georgy F.Voronoi

• Boris N.Delaunay

基于Voronoi图和Delaunay算法的三角化，也叫DT(Delaunay triangle)

Voronoi图

Voronoi图的思路说明：

• 平面上某个区域里，给定一些离散的点
• 将该区域划分成多个网格，每个网格包含一个点
• 要求网格上的分割线刚好是每对点的中线，如下图所示，分割线v0v1上所有点到p0、p1的距离是相等

用严谨的数学公式描述： $V_i = \{ q;\forall_j \neq i \quad |qp_i| \leq |qp_j| \}$

每头狮子i占领一块领地用集合$V_i$表示，$V_i$表示这块领地的所有坐标。领地里任意一点q到狮子的距离$|qp_i|$ 小于该点到其他狮子的距离$|qp_j|$, $p_i$ 是狮子i当前的位置，$p_j$是其他狮子的位置

voronoi 图的生成有很多成熟的方法，比如 "A Voronoi tessellation emerges by radial growth from seeds outward" Voronoi_diagram

Delaunay 三角化

DT(Delaunay triangle) 实现原理比较简单，基于上一步生成的Voronoi图，将处于相邻网格中的点依次连起来即可。

Centroidal Voronoi Tessellation

Centroidal Voronoi Tessellation(中心化的Voronoi网格划分)

普通的Voroni图网格不一定很均匀，导致三角化也不均匀，需要做些优化。

如上图左1是原始的Voronoi图，网格大小不一。 处理的思路：

1. 生成Voronoi图
2. 计算每个网格的中心
3. 将每个点移到该网格的中心
4. 重复上面步骤，直至收敛

这个算法以Stuart P. Lloyd的名字命名。实现思路和k-means类似，都是求均值，逐步迭代。维基百科 Lloyd's algorithm

代码实现

我找了两个版本的实现 python C++

python实现代码看起来更简洁易懂

python 实现cvt

代码参考：

github.com/g1n0st/GAME…

代码不长，我加上了注释,代码中有一处有点问题： 求多边形的中心，不能简单的对所有点加和求均值，应该将多边形剖分成n个三角形，再对所有的三角形中心加和求均值，这里方便理解，暂不修改。参考 计算多边形中心 多边形centroid的计算方法

# coding=utf-8
# This is a sample Python script.

# Press ⌃R to execute it or replace it with your code.
# Press Double ⇧ to search everywhere for classes, files, tool windows, actions, and settings.


import numpy as np
import matplotlib.pyplot as plt
import platform
from scipy.spatial import Delaunay, Voronoi, voronoi_plot_2d
from scipy import argmin, inner


def cvt(samples, times, rnum):
    # rnum = 100 s_num = 10000
    s_num = rnum * rnum
    # 随机生成20个点，'2'表示是2纬的坐标
    p = np.random.random((samples, 2))
    # 10000个点均匀的分摊到0~1之间
    interval = np.linspace(0.0, 1.0, rnum)
    sx = np.zeros(s_num)
    sy = np.zeros(s_num)
    s = np.zeros([s_num, 2])

    # 将100 * 100个离散的点均匀的分布在 1 * 1的区域内
    k = 0
    for j in range(rnum):
        for i in range(rnum):
            sx[k], sy[k] = interval[i], interval[j]
            s[k, 0], s[k, 1] = interval[i], interval[j]
            k += 1

    # 迭代20次
    for it in range(times):
        # 源代码是每迭代3次，生成一张图，我这里为了生成完整的迭代过程，每一帧都会生成图片
        # if it % 3 == 0:
        #     vor = Voronoi(p)
        #     voronoi_plot_2d(vor)
        #     plt.title('Iteration Time: ' + str(it))
        #     plt.show()
        #

        # 生成voronoi图
        vor = Voronoi(p)
        # 绘制voronoi图
        voronoi_plot_2d(vor)
        plt.title('Iteration Time: ' + str(it))
        # 图片保存到本地
        plt.savefig('./Iteration Time: ' + str(it))
        # 展示图片
        plt.show()

        # 遍历100 * 100个点，计算每个点和之前随机成的采样点(20个)的距离
        # 生成一个100 * 100的数组，存储索引，索引表示离哪个采样点最近
        k = [argmin([inner(gg - ss, gg - ss) for gg in p]) for ss in s]

        # w是权重数组，长度为100 * 100，值均为1
        w = np.ones(s_num)
        # 计算每个点对应有多少个离它最近的点，m长度为20
        # 可以查下python bincount的用法，否则不好理解这个计算的意义.python 确实方便，一行代码放到别的语言可能得写十几行
        m = np.bincount(k, weights=w)

        #  用bincount的方式，计算Voronoi图中每个网格的中心，每个点等于周边像素的坐标和
        #  new_x、new_y为网格的中心，
        #  sx、sy 是像素点的坐标，这里作为权重，很巧妙，
        new_x = np.bincount(k, weights=sx)
        new_y = np.bincount(k, weights=sy)

        #  更新采样点的位置
        for i in range(samples):
            if m[i] > 0:
                new_x[i] = new_x[i] / float(m[i])
                new_y[i] = new_y[i] / float(m[i])

        p[:, 0], p[:, 1] = new_x[:], new_y[:]


# python的内置环境变量，判断是从当前文件执行，非import
if __name__ == '__main__':
    # 采样20个点，迭代20次，分辨率为100 * 100
    cvt(20, 20, 100)

c++ 实现cvt

篇幅有限，c++的实现和下篇文章《点画实线》一起吧

参考

[1] Voronoi_diagram-wikipedia: en.wikipedia.org/wiki/Vorono…

[2] 维基百科-Lloyd's algorithm: en.wikipedia.org/wiki/Lloyd%…

[3] 计算多边形中心: www.jianshu.com/p/39ef232ad…

[4] 计算多边形中心: jingsam.github.io/2016/10/05/…

老铁，点个在看呗 欢迎关注公众号:sumsmile /专注图形渲染、Android开发