建立網格面


在〈曲線與曲面〉中談過的 parametricSurface,無法指定厚度,而 interpPlate 雖然可以指定厚度,不過其內插求曲面的方式並不適用多數場合。

這時可以試著自行構造曲面,當然,就算是指定面的點,也有不同的方式,最通用的一種應該是像 matplotlib 的 plot_surface,以網格的方式指定面上每個點(可參考〈Matplotlib 立體圖〉)。

在這邊規範,若有個函式 f(x, y),產生點的方式要是:

points = []
for y in range(min_value, max_value, step):
    row = []
    for x in range(min_value, max_value, step)
        row.append([x, y, f(x, y)])
    points.append(row)

先採用 f(x, y) 的形式,只是為了便於後續後說明,實際上只要能產生相對應格式的 points,你想用這一篇文章即將實現的函式,來建立一個環或球,基本上也不是問題(可以參考〈NumPy 與環面(一)〉、〈NumPy 與環面(二)〉)。

因為在〈實作 polyhedron〉中建立的 polyhedron,可以接受點與面索引來建立 3D 實體,只要能從方才的 points 自動產生這兩個資料就可以了,因為 points 採用網格的方式,切割三角形就很容易。例如,你有以下的點:

建立網格面

每個點往右與往右上各取一個點,就可以構成三角形:

建立網格面

每個點往右上與往上各取一個點,也可以構成三角形:

建立網格面

這樣就處理完一個方格了,接著就是使用迴圈處理完每個方格:

建立網格面

這邊規範看向曲面為正面,另一面為背面,頂點的順序為看向曲面時的逆時針方向。

那麼厚度呢?雖然對於 f(x, y) 之類的函式,單純地將面上每個點上升、下降一個高度,可以建立起厚度,這種作法計算上也簡單許多,不過,這會限縮曲面函式的用途,例如,如果你的面並不是上下有厚度,而是左右有厚度呢?

比較好的方式是,採用頂點法向量,來位移面上的點,這一方面是因為,points 指定的點,實際上就是用來作為 3D 實體的頂點。

在實作上,照理來說,是要收集某頂點周圍的面,將這些面各自的法向量加總後平均,作為某頂點的法向量,不過,CadQuery 的 Face 提供了 makeSplineApprox,你可以提供控制點,建立一個基於 B-spline 的曲面,它有個 normalAt,可以計算指定點的法向量,這可以省去不少自行實作求頂點法向量的功夫。

有了正面與背面的點,剩下的就是計算面索引的功夫而已,這部份就只有一個原則,別粗心算錯就好了,底下程式碼是包含了〈實作 polyhedron〉中建立的 polyhedron 函式,以及這邊新實作的 surface 函式:

from cadquery import Vector, Vertex, 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 surface(points, thickness):
    leng_row = len(points)
    leng_col = len(points[0])
    leng_pts = leng_col * leng_row

    def _all_pts():
        half_thickness = thickness / 2

        # 轉為向量   
        vectors = [[Vector(*p) for p in row] for row in points]

        # 用於查詢指定位置的法向量
        # makeSplineApprox 的控制點以每行(column)為一組
        # 因此需要將 points 轉置
        face = Face.makeSplineApprox([[
                     Vector(*points[ri][ci]) 
                for ri in range(leng_row)
            ] for ci in range(leng_col)]
        )

        front_thicken_pts = [] # 正面的點
        back_thicken_pts = []  # 背面的點
        for row in vectors:
            for vt in row:
                # 頂點處的單位法向量
                n = face.normalAt(vt).normalized()
                # 計算正面的點並收集
                v = vt + n.multiply(half_thickness)
                front_thicken_pts.append([v.x, v.y, v.z])
                # 計算背面的點並收集
                v = vt + n.multiply(-half_thickness)
                back_thicken_pts.append([v.x, v.y, v.z])
        return front_thicken_pts + back_thicken_pts

    def _all_faces():
        # 正面索引
        front_faces = []
        for ri in range(leng_row - 1):
            for ci in range(leng_col - 1):
                front_faces.append([ci + leng_col * ri, (ci + 1) + leng_col * ri, (ci + 1) + leng_col * (ri + 1)])
                front_faces.append([ci + leng_col * ri, (ci + 1) + leng_col * (ri + 1), ci + leng_col * (ri + 1)])

        # 背面索引        
        back_faces = [[f[2] + leng_pts, f[1] + leng_pts, f[0] + leng_pts] for f in front_faces]

        # 四個邊的索引
        side_faces1 = []
        for ci in range(leng_col - 1):
            side_faces1.append([ci, ci + leng_pts, ci + 1])
            side_faces1.append([ci + leng_pts, ci + leng_pts + 1, ci + 1])

        side_faces2 = []
        rx = leng_col - 1
        for ri in range(leng_row - 1):
            side_faces2.append([rx + (ri + 1) * leng_col + leng_pts, rx + (ri + 1) * leng_col, rx + ri * leng_col])
            side_faces2.append([rx + ri * leng_row + leng_pts, rx + (ri + 1) * leng_col + leng_pts, rx + ri * leng_col])

        side_faces3 = []
        for ci in range(leng_pts - leng_col, leng_pts - 1):
            side_faces3.append([ci + 1, ci + leng_pts, ci])
            side_faces3.append([ci + 1, ci + leng_pts + 1, ci + leng_pts])

        side_faces4 = []
        for ri in range(leng_row - 1):
            side_faces4.append([ri * leng_col, (ri + 1) * leng_col, (ri + 1) * leng_col + leng_pts])
            side_faces4.append([ri * leng_col, (ri + 1) * leng_col + leng_pts, ri * leng_row + leng_pts])

        return front_faces + back_faces + side_faces1 + side_faces2 + side_faces3 + side_faces4

    return polyhedron(_all_pts(), _all_faces())

先來個簡單的雙拋物面:

# 要結合以上的 polyhedron 及 surface

def paraboloid(x, y):
    return [x, y, ((y ** 2) - (x ** 2)) / 4]

min_value = -30
max_value = 30
step = 5
thickness = 0.5

points = [[
        paraboloid(x / 10, y / 10) 
    for x in range(min_value, max_value, step)
] for y in range(min_value, max_value, step)]

solid = surface(points, thickness)
show_object(solid)

這可以建立以下的模型,只要改變 step 的大小,就能改變細緻的程度:

建立網格面

來看看另一個漣漪模型:

from math import cos, sqrt, radians

# 要結合以上的 polyhedron 及 surface

def ripple(x, y):
    n = radians(sqrt(x ** 2 + y ** 2))
    return [x, y, 30 * (cos(n) + cos(3 * n))]

min_value = -200
max_value = 200
step = 10
thickness = 5

points = [[
        ripple(x, y) 
    for x in range(min_value, max_value, step)
] for y in range(min_value, max_value, step)]

solid = surface(points, thickness)
show_object(solid)

建立網格面

在厚度建立的部份,往頂點法向量兩個方向各位移半個厚度,只是其中一種方式,畢竟面本來就沒有厚度的概念,你要怎麼方向位移,其實都可以,或許可以有個 direction 參數,能自行指定位移的方向,例如,僅往正面或背面的法向量位移,或是可指定一個向量作為統一的位移方向之類的,這就留給你自己當練習了。