實作轉換矩陣

November 27, 2021

先前的文件透過 SciPy 實作了 hull、Voronoi 與 Delaunay 三角分割等,SciPy 是基於 NumPy,實際上,CadQuery 也有使用 NumPy,也就是,在安裝了 CadQuery 後,你就有 NumPy 可以使用了,想知道 NumPy 的基礎,可以參考〈笨學資料運算〉。

只不過,SciPy 的 spatial 模組跟 3D 建模有關係還可以理解,NumPy 跟 3D 建模又有什麼直接關係呢?當然有關係,3D 建模本身常涉及大量點、線、面資料的處理,CadQuery 底層會基於 NumPy,是可以想像的,不過就使用者層面,可能還是覺得不太會直接運用到 NumPy?

實作 Matrix3D

那就來實作個轉換矩陣吧!因為 CadQuery 雖然有個 Matrix 類別,不過功能有限,而且主要用於處理一些底層的 Open CASCADE 物件,我希望的是,有個通用的轉換矩陣,可以直接處理頂點或向量。

NumPy 陣列本身可以作為陣列來使用,而且實作了 @(Python 3.5 新增,對應的特殊方法是 __matmul__),可用於矩陣相乘:

>>> import numpy
>>> a = numpy.array([
...     [1, 2],
...     [3, 4]
... ])
>>> b = numpy.array([
...     [5, 6],
...     [7, 8]
... ])
>>> a @ b # 矩陣相乘
array([[19, 22],
       [43, 50]])
>>> a * b # 逐元素相乘
array([[ 5, 12],
       [21, 32]])
>>>

記得別錯用了 * 運算,那是用於逐元素相乘。

雖然直接使用 NumPy 的 ndarray 來操作也是可以,不過,將之封裝會是更好的作法:

import numpy

from math import cos, sin, radians

class Matrix3D:
    def __init__(self, m):
            # self.wrapped 是包裹的 numpy.ndarray 實例
        if isinstance(m, numpy.ndarray):
            self.wrapped = m
        else:
            self.wrapped = numpy.array(m)


    # Post-Multiplication (Right-Multiplication)
    def __matmul__(self, that):
        return Matrix3D(self.wrapped @ that.wrapped)


    def transform(self, v):
        vt = (v.x, v.y, v.z, 1) if isinstance(v, Vector) else v + (1,)
        return tuple((self.wrapped @ vt)[:-1])

這邊實作了 __matmul__,為了撰寫矩陣相乘時,看來像個數學公式,採用了後乘(Post-Multiplication)transform 方法會根據目前轉換矩陣的定義,將一個指定的點或向量進行線性轉換。

# 單位矩陣        
_identity = [
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, 1, 0],
    [0, 0, 0, 1]
]

def identity():
    return numpy.array(_identity)

# 平移矩陣
def translation(v):
    return Matrix3D(_translation(v))

# 繞 x 軸旋轉矩陣
def rotationX(angle):
    return Matrix3D(_rotationX(angle))

# 繞 y 軸旋轉矩陣
def rotationY(angle):
    return Matrix3D(_rotationY(angle))

# 繞 z 軸旋轉矩陣
def rotationZ(angle):
    return Matrix3D(_rotationZ(angle))

# 繞指定軸旋轉矩陣
def rotation(direction, angle):
    return Matrix3D(_rotation(direction, angle))


def _translation(v):
    vt = (v.x, v.y, v.z) if isinstance(v, Vector) else v
    return numpy.array([
        [1, 0, 0, vt[0]],
        [0, 1, 0, vt[1]],
        [0, 0, 1, vt[2]],
        [0, 0, 0, 1]
    ]) 

def _rotationX(angle):
    rad = radians(angle)
    c = cos(rad)
    s = sin(rad)
    return numpy.array([
        [1, 0, 0, 0],
        [0, c, -s, 0],
        [0, s, c, 0],
        [0, 0, 0, 1]
    ]) 

def _rotationY(angle):
    rad = radians(angle)
    c = cos(rad)
    s = sin(rad)
    return numpy.array([
        [c, 0, s, 0],
        [0, 1, 0, 0],
        [-s, 0, c, 0],
        [0, 0, 0, 1]
    ]) 

def _rotationZ(angle):
    rad = radians(angle)
    c = cos(rad)
    s = sin(rad)
    return numpy.array([
        [c, -s, 0, 0],
        [s, c, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1]
    ]) 

# 四元數旋轉公式實作
def _rotation(direction, angle):
    dir = direction if isinstance(direction, Vector) else Vector(*direction)
    axis = dir.normalized()
    half_a = radians(angle / 2)
    s = sin(half_a)
    x = s * axis.x
    y = s * axis.y
    z = s * axis.z
    w = cos(half_a)
    x2 = x + x
    y2 = y + y
    z2 = z + z
    xx = x * x2
    yx = y * x2
    yy = y * y2
    zx = z * x2
    zy = z * y2
    zz = z * z2
    wx = w * x2
    wy = w * y2
    wz = w * z2  
    return numpy.array([
        [1 - yy - zz, yx - wz, zx + wy, 0],
        [yx + wz, 1 - xx - zz, zy - wx, 0],
        [zx - wy, zy + wx, 1 - xx - yy, 0],
        [0, 0, 0, 1]
    ]) 

旋轉四面體

來個簡單的範例,將〈實作 polyhedron〉中實作的四面體,繞 z 軸轉一圈,為了便於檢視完整的程式,底下列出全部的程式碼:

import numpy
from math import cos, sin, radians
from cadquery import Vector, Edge, Wire, Solid, Shell, Face

class Matrix3D:
    def __init__(self, m):
        if isinstance(m, numpy.ndarray):
            self.wrapped = m
        else:
            self.wrapped = numpy.array(m)


    # Post-Multiplication (Right-Multiplication)
    def __matmul__(self, that):
        return Matrix3D(self.wrapped @ that.wrapped)


    def transform(self, v):
        vt = (v.x, v.y, v.z, 1) if isinstance(v, Vector) else v + (1,)
        return tuple((self.wrapped @ vt)[:-1])


# 單位矩陣        
_identity = [
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, 1, 0],
    [0, 0, 0, 1]
]

def identity():
    return numpy.array(_identity)

# 平移矩陣
def translation(v):
    return Matrix3D(_translation(v))

# 繞 x 軸旋轉矩陣
def rotationX(angle):
    return Matrix3D(_rotationX(angle))

# 繞 y 軸旋轉矩陣
def rotationY(angle):
    return Matrix3D(_rotationY(angle))

# 繞 z 軸旋轉矩陣
def rotationZ(angle):
    return Matrix3D(_rotationZ(angle))

# 繞指定軸旋轉矩陣
def rotation(direction, angle):
    return Matrix3D(_rotation(direction, angle))

def _translation(v):
    vt = (v.x, v.y, v.z) if isinstance(v, Vector) else v
    return numpy.array([
        [1, 0, 0, vt[0]],
        [0, 1, 0, vt[1]],
        [0, 0, 1, vt[2]],
        [0, 0, 0, 1]
    ]) 

def _rotationX(angle):
    rad = radians(angle)
    c = cos(rad)
    s = sin(rad)
    return numpy.array([
        [1, 0, 0, 0],
        [0, c, -s, 0],
        [0, s, c, 0],
        [0, 0, 0, 1]
    ]) 

def _rotationY(angle):
    rad = radians(angle)
    c = cos(rad)
    s = sin(rad)
    return numpy.array([
        [c, 0, s, 0],
        [0, 1, 0, 0],
        [-s, 0, c, 0],
        [0, 0, 0, 1]
    ]) 

def _rotationZ(angle):
    rad = radians(angle)
    c = cos(rad)
    s = sin(rad)
    return numpy.array([
        [c, -s, 0, 0],
        [s, c, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1]
    ]) 

# 四元數旋轉公式實作
def _rotation(direction, angle):
    dir = direction if isinstance(direction, Vector) else Vector(*direction)
    axis = dir.normalized()
    half_a = radians(angle / 2)
    s = sin(half_a)
    x = s * axis.x
    y = s * axis.y
    z = s * axis.z
    w = cos(half_a)
    x2 = x + x
    y2 = y + y
    z2 = z + z
    xx = x * x2
    yx = y * x2
    yy = y * y2
    zx = z * x2
    zy = z * y2
    zz = z * z2
    wx = w * x2
    wy = w * y2
    wz = w * z2  
    return numpy.array([
        [1 - yy - zz, yx - wz, zx + wy, 0],
        [yx + wz, 1 - xx - zz, zy - wx, 0],
        [zx - wy, zy + wx, 1 - xx - yy, 0],
        [0, 0, 0, 1]
    ]) 
    
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
        )
    )

points = (
    (5, -5, -5), (-5, 5, -5), (5, 5, 5), (-5, -5, 5)
)

faces = (
    (0, 1, 2), (0, 3, 1), (1, 3, 2), (0, 2, 3)
)

n = 6
a_step = 360 / n
for i in range(n):
    # 建立轉換矩陣
    m = rotationZ(i * a_step) @ translation((15, 0, 0)) @ rotationZ(45)
    # 將 points 逐一轉換並建立多面體
    solid = polyhedron([m.transform(p) for p in points], faces)
    show_object(solid)    

這會顯示以下的結果:

實作轉換矩陣