【大作业】太极图形课S1-太极光线追踪（扁平kd树实现）

January 3, 2022 · View on GitHub

成功效果展示

先贴效果图。刚学的 blender， 建模有点粗超，请见谅。因为没有显卡，所以是纯CPU运行（6核 x 3.00GHz）。

mushroom

crystal

简介

基于 taichi 编程语言，实现基于扁平kd树的光线追踪。

主要解决的问题：

3D模型:

  • 目前只支持 PLY ascii 格式。详见 plyreader2.py
  • 因为坐标轴问题，贴图要旋转 −π/2-\pi/2。

文件结构

- raytrace-c.py (晶体模型的渲染程序)
- raytrace-m.py (蘑菇模型的渲染程序)
- Scene.py (场景，摄像机，光线追踪，...)
- plyread2.py (读取PLY模型)
- view_uv.py (看贴图)
- mesh (模型与贴图的文件夹)

运行方式：

  • 晶体 python3 raytrace-c.py
  • 蘑菇 python3 raytrace-m.py
  • 看贴图 python3 view_uv.py

符号说明

  • 射线： p=o+tdp = o + t d. ee 为射线的端点， tt是时间， dd 是方向， 在代码中 ∥d∥F=1\|d\|_F = 1 长度固定成 $1$；
  • 平面： n(x−y)=0n(x - y) = 0. nn 是平面的法向量， xx 为平面上的任何一点， yy 是另一个点。

三角形面片

3D 模型最常用的就是三角形或四边形面片。那么一个核心的问题是如何判断射线与一个三角形面片相交。

一个方法是直接计算。 虎书（Steve Marschner and Peter Shirley. Fundamentals of Computer Graphics 4th ed.）的 4.4.2 有详细讲解。 但这种算法每一次计算相交都要算一个 $3 \times 3$ 的行列式。复杂度比较高。

这里提供另一个思路。 三角形面片是一个二位的物体，射线跟它相交必然先跟三角形所在的平面相交。 我们可以先计算射线与平面相交，然后在看交点是不是在三角形里面。

计算射线p=o+tdp = o + td 到平面 n(x−(o+td))=0n(x - (o + td)) = 0 相交的时间为 t=n(x−o)/(nd)t = n(x - o) / (nd)。 如果tt大于前面最短相交的时间，那么这个三角形不是需要的，可以直接略过。 如果tt小于前面最短相交的时间,那就判断交点是不是在三角形里面。

对于一个三角形：

   a
  / \
 /   \
x --- b

有交点坐标 pp， 参数α,β\alpha, \beta 满足

p=x+α(a−x)+β(b−x)p = x + \alpha (a-x) + \beta (b-x)

把xx移到左边，有:

p−x=α(a−x)+β(b−x)p - x = \alpha (a-x) + \beta (b-x)

为了表示方便，令p^=p−x\hat{p} = p - x, a^=a−x\hat{a} = a - x, b^=b−x\hat{b} = b - x， 展开来就是

p^[0]=αx^[0]+βb^[0]p^[1]=αx^[1]+βb^[1]p^[2]=αx^[2]+βb^[2]\begin{aligned} \hat{p}[0] &= \alpha \hat{x}[0] + \beta \hat{b}[0]\\ \hat{p}[1] &= \alpha \hat{x}[1] + \beta \hat{b}[1]\\ \hat{p}[2] &= \alpha \hat{x}[2] + \beta \hat{b}[2]\\ \end{aligned}

这里有三条等式，但只有两条的线性独立的（线性代数| 可以大概理解成3式等于1式+2式）。因此求解的时候只需用到两个方程。

那么选那两条方程呢？ 我们可以先把他们两两组合，如果能求解的话就求解，如果求不了（如选等式1和等式2，但发现等式1=等式2 乘 某个数），就换一个组合。因为三角形面积不为0， 一定是有一个组合能解的。

具体来说就是: 判断三个行列式是否为$0$

det⁡(A0)=∣a^[1]a^[2]b^[1]b^[2]∣=a^[1] b^[2]−a^[2] b^[1]det⁡(A1)=∣a^[0]a^[2]b^[0]b^[2]∣=a^[0] b^[2]−a^[2] b^[0]det⁡(A2)=∣a^[0]a^[1]b^[0]b^[1]∣=a^[0] b^[1]−a^[1] b^[0]\begin{aligned} \det(A_0) &= \left|\begin{matrix} \hat{a}[1] & \hat{a}[2] \\ \hat{b}[1] & \hat{b}[2] \\ \end{matrix}\right| &= \hat{a}[1] ~ \hat{b}[2] - \hat{a}[2] ~ \hat{b}[1] \\ \det(A_1) &= \left|\begin{matrix} \hat{a}[0] & \hat{a}[2] \\ \hat{b}[0] & \hat{b}[2] \\ \end{matrix}\right| &= \hat{a}[0] ~ \hat{b}[2] - \hat{a}[2] ~ \hat{b}[0] \\ \det(A_2) &= \left|\begin{matrix} \hat{a}[0] & \hat{a}[1] \\ \hat{b}[0] & \hat{b}[1] \\ \end{matrix}\right| &= \hat{a}[0] ~ \hat{b}[1] - \hat{a}[1] ~ \hat{b}[0] \\ \end{aligned}

发现三个行列式都与射线无关，所以可以在初始化阶段计算，节省时间。

计算三个行列式后，就可以用下面三种情况的任意一种直接解算 α,β\alpha,\beta。计算量从原来的 $3 \times 3 行列式变成 \2 \times 2$ 行列式。

如果 det⁡(A0)≠0\det(A_0) \ne 0,

α=(p^[0]a^[1]+p^[1]a^[2])/det⁡(A0)β=(p^[0]b^[1]+p^[1]b^[2])/det⁡(A0)\begin{aligned} \alpha &= (\hat{p}[0] \hat{a}[1] + \hat{p}[1] \hat{a}[2]) / \det(A_0) \\ \beta &= (\hat{p}[0] \hat{b}[1] + \hat{p}[1] \hat{b}[2]) / \det(A_0) \\ \end{aligned}

如果 det⁡(A1)≠0\det(A_1) \ne 0,

α=(p^[0]a^[0]+p^[1]a^[2])/det⁡(A1)β=(p^[0]b^[0]+p^[1]b^[2])/det⁡(A1)\begin{aligned} \alpha &= (\hat{p}[0] \hat{a}[0] + \hat{p}[1] \hat{a}[2]) / \det(A_1) \\ \beta &= (\hat{p}[0] \hat{b}[0] + \hat{p}[1] \hat{b}[2]) / \det(A_1) \\ \end{aligned}

如果 det⁡(A2)≠0\det(A_2) \ne 0,

α=(p^[0]a^[0]+p^[1]a^[1])/det⁡(A2)β=(p^[0]b^[0]+p^[1]b^[1])/det⁡(A2)\begin{aligned} \alpha &= (\hat{p}[0] \hat{a}[0] + \hat{p}[1] \hat{a}[1]) / \det(A_2) \\ \beta &= (\hat{p}[0] \hat{b}[0] + \hat{p}[1] \hat{b}[1]) / \det(A_2) \\ \end{aligned}

解算出α,β\alpha, \beta 就可以判断点pp是不是在三角形里面

α>0; β>0; α+β<1.\alpha > 0;~ \beta > 0;~ \alpha + \beta < 1.

如果是四边形面片，前面的过程相同，判断α,β\alpha, \beta时，使用

α>0; β>0; α<1; β<1.\alpha > 0;~ \beta > 0;~ \alpha < 1; ~ \beta < 1.

(原理是向量加法)

    a
   /
  /---p
 /   /
x---------b

这里的 α,β\alpha, \beta 也是定位贴图的依据。

kd树生成

因为 taichi 既不支持递归又不支持在 ti.func里创建动态内存(超级大坑，很多树都用不了)，所以kd树是在Python环境里创建，然后再转化到 taichi。

先举一个例子： 比如要在全国范围内找一个人，可以先找他所在的省，然后是市，一级一级往下找。这就比地毯式搜索快很多。

这点在判断相交时也适用。我们可以把物体组成一系列的框，先找一个大框，然后再从大框里找小框，一直找下去。

怎么画框？常见的方法是画分割线，分割线左边的是一个框，右边的是另一个。在中间的就两边都放。

但在找框的时候是要递归或队列的。。。

既没有动态内存又没有递归... 很多常用的方法都用不上 ``/(ㄒ^ㄒ)/```。 在多次试错后，终于找到了一个方法。

根据taichi的特性，计算必须在循环里面完成，那么判断框框的位置关系必须在一次判断就出结果。对于好分离的框框，不用考虑太多。对于重叠的情况

|------
|     |
|   --|----
|   | |   |
|---|--   |
    |-----|

|---------|
|         |
|     ----|
|     |   |
|     ----|
|---------|

就不好办了。

那么我们将画3个框：一个是绝对的左边，一个是绝对的右边，还有一个是中间的。然后Python递归就可以创建树了。

分割线的确定

分割线要求左边的跟右边的物体数量接近，并且处在中间的物体要少。 根据这个，对于每一个维度，网格搜索（偷懒的做法），找一条分割线使得 abs(cl−cr)+0.5cm\text{abs}(c_l - c_r) + 0.5 c_m 最小， cl,cm,crc_l, c_m, c_r 分别为左中右物体数量。

扁平化kd树

因为 taichi 的特性。。。这里将前面建好的树压成一维。

数据格式为（前面的数字表示相对偏移量）：

# [kdflat]
# 0: corner index; 左下和右上角坐标的位置
# 1: mark; 标志符号
# if mark >= 0: # 偏移量区
#   1: left index; 左树的偏移量
#   2: right index; 右树的偏移量
#   3: mid index; 中间物体偏移量
#   4： lr data index; 左右树物体偏移量
# else: # 数据区
#   2: data len (n); 数据长度
#   3+0: triangle index 0; 三角形的索引。。。
#   3+1: triangle index 1
#   ...: .....
#   3+n-1: triangle index n-1
#
# next block

扁平的遍历

以 0 偏移节点为根

  • 先遍历中间
  • 如果只与左边相交：
    • 以左树偏移节点作为根
  • 如果只与右边相交
    • 以右树偏移节点作为根
  • 都相交
    • 遍历 lr data index;
  • 循环

只与左右相交很好理解。同时与左右相交就是下图的情况。虽然打到左框，但是射线没有与左框内任何物体相交。因为要一次计算内完成，所以只能把左右物体全遍历。

    |------|   |------|
    |   A  |   |      |
-----------------> B  |
    |   A  |   |      |
    |------|   |------|

致谢

超级感谢太极图形， 天添老师，所有助教们和所有为课程付出的人员！！！