Voronoi 2D/3D

我用 OpenSCAD 創作的作品中，有一些是基於 Voronoi，主要是基於 dotSCAD 中的 Voronoi 函式或模組來完成，想在 OpenSCAD 中實作 Voronoi 蠻辛苦的，用其他語言是簡單一些，不過也需要一點功夫，有興趣認識原理的話，可以參考〈玩轉 p5.js〉中 Voronoi 圖形、Delaunay 三角分割的部份。

2D Voronoi

SciPy 的 sptial 提供了 Voronoi 函式，在〈玩玩 spatial 模組〉有談到一些，現在我們在 CadQuery 中玩玩，因為不借助 spatial.voronoi_plot_2d 來畫 Voronoi 了，需要知道更多 spatial.Voronoi 的細節就是了。

先從簡單的開始：

import numpy as np
from scipy import spatial
import cadquery as cq

def polygon(points):
    return cq.Workplane().polyline(points).close()

n = 50
points = np.random.rand(n, 2)
voronoi = spatial.Voronoi(points)

r = cq.Workplane()
for region in voronoi.regions:
    if len(region) > 0 and not (-1 in region):
        r.add(polygon([voronoi.vertices[i] for i in region]))

spatial.Voronoi 接受一組點，通常是隨機點，建立的 Voronoi 圖，每一塊勢力稱為一個 region，spatial.Voronoi 傳回的物件上，vertices 是這些 region 用到的全部頂點，regions 包含了一組 region，每個 region 由一組頂點索引構成。

要注意的是，在 spatial.Voronoi 的預設值下，regions 一定會有個 []，而且 region 中會包含 -1 索引，這是因為其底層使用了 Qhull 程式庫，預設會有個無限遠的點，用來增加 Voronoi 計算時的精度，無限遠點的 regin 就是那個空 region，頂點索引是 -1。

這並不是什麼大問題，如上撰寫程式，只在 region 不為空且不包含 -1 時，使用其索引來取得頂點以建立多邊形就可以了：

Voronoi 2D/3D

如果你覺得這麼寫不漂亮，也可以使用 point_region，它包含了指定的隨機點所在的 region，是在 regions 的哪個索引處，至於那個 -1 索引，就看你怎麼處理了，方才是遇到 -1 索引，直接拋棄該 region 不畫，你也可以選擇忽略無限遠的點：

import numpy as np
from scipy import spatial
import cadquery as cq

def polygon(points):
    return cq.Workplane().polyline(points).close()

n = 50
points = np.random.rand(n, 2)
voronoi = spatial.Voronoi(points)

r = cq.Workplane()
for region_i in voronoi.point_region:
    region = voronoi.regions[region_i]
    r.add(polygon([voronoi.vertices[i] for i in region if i != -1]))

這畫出來會像是：

Voronoi 2D/3D

3D Voronoi

spatial.Voronoi 也可以指定 3D 的點，這會建立 3D 版本的 Voronoi，此時的 regions 中每個 region，代表著一個凸多面體的頂點索引，那麼該怎麼在 CadQuery 中建模呢？在〈建立 Convex hull〉中不是有實作了 hull3D 嗎？拿來用就好了！因為 3D Voronoi 的每個 region 就是凸多面體，用 hull3D 來畫還算是有效率：

import numpy as np
from scipy import spatial
import cadquery as cq

from random import randint
from scipy.spatial import ConvexHull
from cadquery import Vector, Edge, Wire, Solid, Shell, Face

def polyhedron(points, faces):
    def _edges(vectors, face_indices):
        leng_vertices = len(face_indices)   
        return (
            Edge.makeLine(
                vectors[face_indices[i]], 
                vectors[face_indices[(i + 1) % leng_vertices]]
            ) 
            for i in range(leng_vertices)
        )

    vectors = [Vector(*p) for p in points]

    return Solid.makeSolid(
        Shell.makeShell(
            Face.makeFromWires(
                Wire.assembleEdges(
                    _edges(vectors, face_indices)
                )
            )
            for face_indices in faces
        )
    )
    
def hull3D(points):
    hull = ConvexHull(points)
    vertices = [points[i] for i in hull.vertices]
    v_i_lookup = {v: i for i, v in enumerate(vertices)}
    faces = [
        [v_i_lookup[points[i]] for i in face]
        for face in hull.simplices
    ]

    return polyhedron(vertices, faces)

def random_color():
    return (randint(0, 255), randint(0, 255), randint(0, 255))

    # 3D Voronoi

n = 40
points = np.random.rand(n, 3)
voronoi = spatial.Voronoi(points)

box = cq.Workplane().box(1, 1, 1).translate((.5, .5, .5))

for region_i in voronoi.point_region:
    region = voronoi.regions[region_i]
    if not (-1 in region):
        p = hull3D([tuple(voronoi.vertices[i]) for i in region])
        show_object(box.intersect(p), options = { 'color' : random_color() })

為了比較容易呈現 3D Voronoi 的效果，選擇忽略無限遠點，並且用了個方塊來交集每個 region，然後用隨機顏色顯示出來：

Voronoi 2D/3D