NSDT工具推荐Three.js AI纹理开发包 - YOLO合成数据生成器 - GLTF/GLB在线编辑 - 3D模型格式在线转换 - 可编程3D场景编辑器 - REVIT导出3D模型插件 - 3D模型语义搜索引擎 - AI模型在线查看 - Three.js虚拟轴心开发包 - 3D模型在线减面 - STL模型在线切割 - 3D道路快速建模

Matplotlib 有一个非常漂亮的 3D 界面,具有许多功能(和一些限制),在用户中非常受欢迎。 然而,对于某些用户(或者可能对于大多数用户)来说,3D 仍然被认为是某种黑魔法。 因此,我想在这篇文章中解释一下,一旦你理解了一些概念,3D 渲染就会变得非常简单。 为了证明这一点,我们将使用 60 行 Python 代码和一个 Matplotlib 调用来渲染上面的兔子,而不使用 3D 轴。

如果你手头的模型不是.OBJ格式,可以用NSDT 3DConvert这个在线3D格式转换工具将其转换为.OBJ格式

https://3dconvert.nsdt.cloud

1、加载兔子

首先,我们需要加载模型。 我们将使用斯坦福兔子的简化版本。 该文件使用wavefront .ob格式,这是最简单的格式之一,所以让我们制作一个非常简单(但容易出错)的加载器,它将完成这篇文章(和这个模型)的工作:

V, F = [], []
with open("bunny.obj") as f:
   for line in f.readlines():
       if line.startswith('#'):
           continue
       values = line.split()
       if not values:
           continue
       if values[0] == 'v':
           V.append([float(x) for x in values[1:4]])
       elif values[0] == 'f':
           F.append([int(x) for x in values[1:4]])
V, F = np.array(V), np.array(F)-1

V 现在是一组顶点(如果你愿意,也可以是 3D 点), F 是一组面(= 三角形)。 每个三角形由相对于顶点数组的 3 个索引来描述。 现在,让我们标准化顶点,使整个兔子适合单位框:

V = (V-(V.max(0)+V.min(0))/2)/max(V.max(0)-V.min(0))

现在,我们可以通过仅获取顶点的 x,y 坐标并去掉 z 坐标来初步查看模型。 为此,我们可以使用强大的 PolyCollection 对象,它可以有效地渲染非规则多边形的集合。 因为我们想要渲染一堆三角形,所以这是一个完美的匹配。 因此,我们首先提取三角形并去掉 z 坐标:

T = V[F][...,:2]

我们现在可以渲染它:

fig = plt.figure(figsize=(6,6))
ax = fig.add_axes([0,0,1,1], xlim=[-1,+1], ylim=[-1,+1],
                  aspect=1, frameon=False)
collection = PolyCollection(T, closed=True, linewidth=0.1,
                            facecolor="None", edgecolor="black")
ax.add_collection(collection)
plt.show()

你应该得到这样的东西(bunny-1.py):

2、透视投影

我们刚刚所做的渲染实际上是正交投影,而顶部的兔子使用透视投影:

在这两种情况下,定义投影的正确方法是首先定义观看体积,即我们想要在屏幕上渲染的 3D 空间中的体积。 为此,我们需要考虑 6 个剪裁平面(左、右、上、下、远、近),它们相对于相机封闭观察体积(视锥体)。 如果我们定义相机位置和观察方向,则每个平面都可以用单个标量来描述。 一旦我们有了这个观看体积,我们就可以使用正交投影或透视投影投影到屏幕上。

对我们来说幸运的是,这些投影是众所周知的并且可以使用 4x4 矩阵来表示:

def frustum(left, right, bottom, top, znear, zfar):
    M = np.zeros((4, 4), dtype=np.float32)
    M[0, 0] = +2.0 * znear / (right - left)
    M[1, 1] = +2.0 * znear / (top - bottom)
    M[2, 2] = -(zfar + znear) / (zfar - znear)
    M[0, 2] = (right + left) / (right - left)
    M[2, 1] = (top + bottom) / (top - bottom)
    M[2, 3] = -2.0 * znear * zfar / (zfar - znear)
    M[3, 2] = -1.0
    return M

def perspective(fovy, aspect, znear, zfar):
    h = np.tan(0.5*radians(fovy)) * znear
    w = h * aspect
    return frustum(-w, w, -h, h, znear, zfar)

对于透视投影,我们还需要指定孔径角(或多或少)设置近平面相对于远平面的大小。 因此,对于高光圈,你会得到很多“变形”。

但是,如果查看上面的两个函数,你会发现它们返回 4x4 矩阵,而我们的坐标是 3D。 那么如何使用这些矩阵呢? 答案是齐次坐标。 长话短说,齐次坐标最适合处理 3D 中的变换和投影。 在我们的例子中,因为我们处理的是顶点(而不是向量),所以我们只需将 1 作为第四个坐标 (w) 添加到所有顶点。 然后我们可以使用点积应用透视变换。

V = np.c_[V, np.ones(len(V))] @ perspective(25,1,1,100).T

最后一步,我们需要重新标准化齐次坐标。 这意味着我们将每个变换后的顶点除以最后一个分量 (w),以便每个顶点始终具有 w=1。

V /= V[:,3].reshape(-1,1)

现在我们可以再次显示结果(bunny-2.py):

哦,奇怪的结果。 怎么回事? 问题是相机实际上在兔子体内。 为了获得正确的渲染效果,我们需要将兔子移离相机或将相机移离兔子。 我们再做后面的事吧。 相机当前位于 (0,0,0) 并沿 z 方向向上看(由于截锥体变换)。 因此,我们需要在透视变换之前将相机在 z 负方向上稍微移开一点:

V = V - (0,0,3.5)
V = np.c_[V, np.ones(len(V))] @ perspective(25,1,1,100).T
V /= V[:,3].reshape(-1,1)

现在你应该获得(bunny-3.py):

3、模型、视图、投影(MVP)

可能不太明显,但最后的渲染实际上是透视变换。 为了让它更明显,我们将旋转兔子。 为此,我们需要一些旋转矩阵(4x4),同时我们也可以定义平移矩阵:

def translate(x, y, z):
    return np.array([[1, 0, 0, x],
                     [0, 1, 0, y],
                     [0, 0, 1, z],
                     [0, 0, 0, 1]], dtype=float)

def xrotate(theta):
    t = np.pi * theta / 180
    c, s = np.cos(t), np.sin(t)
    return np.array([[1, 0,  0, 0],
                     [0, c, -s, 0],
                     [0, s,  c, 0],
                     [0, 0,  0, 1]], dtype=float)

def yrotate(theta):
    t = np.pi * theta / 180
    c, s = np.cos(t), np.sin(t)
    return  np.array([[ c, 0, s, 0],
                      [ 0, 1, 0, 0],
                      [-s, 0, c, 0],
                      [ 0, 0, 0, 1]], dtype=float)

现在,我们将根据模型(局部变换)、视图(全局变换)和投影来分解要应用的变换,以便我们可以计算一个可以同时完成所有操作的全局 MVP 矩阵:

model = xrotate(20) @ yrotate(45)
view  = translate(0,0,-3.5)
proj  = perspective(25, 1, 1, 100)
MVP   = proj  @ view  @ model

现在我们写:

V = np.c_[V, np.ones(len(V))] @ MVP.T
V /= V[:,3].reshape(-1,1)

你应该得到(bunny-4.py):

现在让我们稍微调整一下光圈,以便你可以看到差异。 请注意,我们还必须调整与相机的距离,以使兔子具有相同的外观尺寸(bunny-5.py):

4、深度排序

现在让我们尝试填充三角形(bunny-6.py)

正如你所看到的,结果很“有趣”并且完全错误。 问题是 PolyCollection 会按照给定的顺序绘制三角形,而我们希望从后到前绘制三角形。 这意味着我们需要根据它们的深度对它们进行排序。 好消息是,当我们应用 MVP 转换时,我们已经计算了这些信息。 它存储在新的 z 坐标中。 然而,这些 z 值是基于顶点的,而我们需要对三角形进行排序。 因此,我们将平均 z 值作为三角形深度的代表。 如果三角形相对较小且不相交,则效果很好:

T =  V[:,:,:2]
Z = -V[:,:,2].mean(axis=1)
I = np.argsort(Z)
T = T[I,:]

现在一切都渲染正确了(bunny-7.py):

让我们使用深度缓冲区添加一些颜色。 我们将根据每个三角形的深度为其着色。 PolyCollection 对象的美妙之处在于你可以使用 NumPy 数组指定每个三角形的颜色,所以让我们这样做:

zmin, zmax = Z.min(), Z.max()
Z = (Z-zmin)/(zmax-zmin)
C = plt.get_cmap("magma")(Z)
I = np.argsort(Z)
T, C = T[I,:], C[I,:]

现在一切都渲染正确了(bunny-8.py):

最终脚本有 57 行(但很难读):

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection

def frustum(left, right, bottom, top, znear, zfar):
    M = np.zeros((4, 4), dtype=np.float32)
    M[0, 0] = +2.0 * znear / (right - left)
    M[1, 1] = +2.0 * znear / (top - bottom)
    M[2, 2] = -(zfar + znear) / (zfar - znear)
    M[0, 2] = (right + left) / (right - left)
    M[2, 1] = (top + bottom) / (top - bottom)
    M[2, 3] = -2.0 * znear * zfar / (zfar - znear)
    M[3, 2] = -1.0
    return M
def perspective(fovy, aspect, znear, zfar):
    h = np.tan(0.5*np.radians(fovy)) * znear
    w = h * aspect
    return frustum(-w, w, -h, h, znear, zfar)
def translate(x, y, z):
    return np.array([[1, 0, 0, x], [0, 1, 0, y],
                     [0, 0, 1, z], [0, 0, 0, 1]], dtype=float)
def xrotate(theta):
    t = np.pi * theta / 180
    c, s = np.cos(t), np.sin(t)
    return np.array([[1, 0,  0, 0], [0, c, -s, 0],
                     [0, s,  c, 0], [0, 0,  0, 1]], dtype=float)
def yrotate(theta):
    t = np.pi * theta / 180
    c, s = np.cos(t), np.sin(t)
    return  np.array([[ c, 0, s, 0], [ 0, 1, 0, 0],
                      [-s, 0, c, 0], [ 0, 0, 0, 1]], dtype=float)
V, F = [], []
with open("bunny.obj") as f:
    for line in f.readlines():
        if line.startswith('#'):  continue
        values = line.split()
        if not values:            continue
        if values[0] == 'v':      V.append([float(x) for x in values[1:4]])
        elif values[0] == 'f' :   F.append([int(x) for x in values[1:4]])
V, F = np.array(V), np.array(F)-1
V = (V-(V.max(0)+V.min(0))/2) / max(V.max(0)-V.min(0))
MVP = perspective(25,1,1,100) @ translate(0,0,-3.5) @ xrotate(20) @ yrotate(45)
V = np.c_[V, np.ones(len(V))]  @ MVP.T
V /= V[:,3].reshape(-1,1)
V = V[F]
T =  V[:,:,:2]
Z = -V[:,:,2].mean(axis=1)
zmin, zmax = Z.min(), Z.max()
Z = (Z-zmin)/(zmax-zmin)
C = plt.get_cmap("magma")(Z)
I = np.argsort(Z)
T, C = T[I,:], C[I,:]
fig = plt.figure(figsize=(6,6))
ax = fig.add_axes([0,0,1,1], xlim=[-1,+1], ylim=[-1,+1], aspect=1, frameon=False)
collection = PolyCollection(T, closed=True, linewidth=0.1, facecolor=C, edgecolor="black")
ax.add_collection(collection)
plt.show()

原文链接:Custom 3D engine in Matplotlib

BimAnt翻译整理,转载请标明出处