實作轉換矩陣
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 方法會根據目前轉換矩陣的定義，將一個指定的點或向量進行線性轉換。
矩陣的實作會採用列為主（row-major），例如以下的矩陣：
實現時會逐列編寫，也就是如下看待陣列中的矩陣元素：
[
row1
row2
row3
]
每一列使用一個 list，構成二維的 list：
[
[1, 0, tx],
[0, 1, ty],
[0, 0, 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]
])
旋轉四面體
來個簡單的範例，將〈實作 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)
這會顯示以下的結果：
