三角形 (Triangles)

三角形坐标 (Triangle Coordinates)

三角形内的区域可以用三个非共线的点$A$,$B$,$C$ 定义，并可以用代数形式写为三个点的凸组合:

$$\begin{gathered} P(w,u,v)=wA+uB+vC \\ s.t.0\le w,u,v \\ 1=w+u+v \end{gathered}$$

或者等价于

$$\begin{gathered} P(u,v)=(1-u-v)A+uB+vC \\ =A+u(B-A)+v(C-A) \\ s.t.0\le u,v \\ u+v\le 1 \end{gathered}$$

这里$u$,$v$ 被称为 "重心坐标". 如果我们删除$u$ 和$v$ 的范围，我们就得到了三角形所在平面的表达式。用下面的小部件来了解$u$ 和$v$ 的行为.

one_triangle = t.tensor([[0, 0, 0], [3, 0.5, 0], [2, 3, 0]])
A, B, C = one_triangle
x, y, z = one_triangle.T

fig = setup_widget_fig_triangle(x, y, z)

@interact(u=(-0.5, 1.5, 0.01), v=(-0.5, 1.5, 0.01))
def response(u=0.0, v=0.0):
    P = A + u * (B - A) + v * (C - A)
    fig.data[2].update({"x": [P[0]], "y": [P[1]]})

display(fig)

三角 - 射线相交

给定一条射线的原点$O$ 和方向$D$, 我们的相交算法将包含以下两个步骤:

• 通过联立方程$P(u,v)=P(s)$ 来计算交点坐标.
• 检查$u$ 和$v$ 是否满足范围.

展开等式$P(u,v)=P(s)$, 我们有:

$$\begin{gathered} A+u(B-A)+v(C-A)=O+sD \\ \Rightarrow \left(\begin{array}{ccc}-D& (B-A)& (C-A)\end{array}\right)\left(\begin{array}{c}s\\ u\\ v\end{array}\right)=(O-A) \\ \Rightarrow \left(\begin{array}{ccc}-D_x& (B-A{)}_x& (C-A{)}_x\\ -D_y& (B-A{)}_y& (C-A{)}_z\\ -D_z& (B-A{)}_z& (C-A{)}_z\end{array}\right)\left(\begin{array}{c}s\\ u\\ v\end{array}\right)=\left(\begin{array}{c}(O-A{)}_x\\ (O-A{)}_y\\ (O-A{)}_z\end{array}\right) \end{gathered}$$

因此我们可以通过求解这个线性方程解出交点的坐标`s`,`u`和`v`.

Exercise - 完成 triangle_ray_intersects

Difficulty: 🔴🔴🔴⚪⚪
Importance: 🔵🔵🔵⚪⚪
你应该花最多15-20分钟在这个练习上.

使用 torch.linalg.solve 和 torch.stack , 完成 triangle_ray_intersects(A, B, C, O, D)

一些提示:

• 如果你有一个零维的张量，形状为 () , 只储存了单个值，请使用 .item() 方法把他转换为普通的 Python 值.
• 如果你有一个形状为 tensor.shape = (3, ...) 的张量，那么你可以用类似 s, u, v = tensor 的方法沿着第一个维度把这个张量分解成三个独立的张量，就和你分解 python 中的列表一样.
  • 注意，如果你想要分解的维度不在第一个维度，有一个很好的替换方法是 s, u, v = tensor.unbind(dim) , 其中 dim 指定了你想要拆分的维度.
• 如果你函数没正常工作，尝试用漂亮的整数制作一个简单的射线和三角形，手动计算是否相交，然后从这开始慢慢调试.

Point = Float[Tensor, "points=3"]

@jaxtyped
@typeguard.typechecked
def triangle_ray_intersects(A: Point, B: Point, C: Point, O: Point, D: Point) -> bool:
    '''
    A: shape (3,), one vertex of the triangle
    B: shape (3,), second vertex of the triangle
    C: shape (3,), third vertex of the triangle
    O: shape (3,), origin point
    D: shape (3,), direction point

    Return True if the ray and the triangle intersect.
    '''
    pass


tests.test_triangle_ray_intersects(triangle_ray_intersects)

Solution

def triangle_ray_intersects(A: Point, B: Point, C: Point, O: Point, D: Point) -> bool:
    '''
    A: shape (3,), one vertex of the triangle
    B: shape (3,), second vertex of the triangle
    C: shape (3,), third vertex of the triangle
    O: shape (3,), origin point
    D: shape (3,), direction point

    Return True if the ray and the triangle intersect.
    '''
    # SOLUTION
    s, u, v = t.linalg.solve(
        t.stack([-D, B - A, C - A], dim=1), 
        O - A
    )
    return ((u >= 0) & (v >= 0) & (u + v <= 1)).item()

渲染单个三角形 (Single-Triangle Rendering)

在仅调用 torch.linagl.solve 的前提下完成 raytrace_triangle .

重塑输出张量并使用 plt.imshow 进行可视化。边缘像素化和锯齿状是正常的 - 使用少量像素是快速调试的好方法.

如果你觉得你的代码已经能正常跑了，请增加像素的数量验证更高分辨率下边缘像素化程度是否降低.

视图和副本 (Views and Copies)

知道你什么时候创建了 Tensor 的副本而不是创建了与原始张量共享数据的视图是很重要的。尽可能使用视图是最好的，可以避免不必要的内存空间占用。但是另一方面，修改视图会修改原始张量，有时候可能会造成一些奇怪的结果。如果你不确定函数是否返回视图，请参阅文档。常用函数返回情况速查:

• torch.expand : 总是返回视图
• torch.view : 总是返回视图
• torch.detach : 总是返回视图
• torch.repeat : 总是复制
• torch.clone : 总是复制
• torch.flip : 总是复制 (和 numpy.flip 不同，后者总是返回视图)
• torch.tensor : 总是复制，但是 PyTorch 推荐使用 .clone().detach() 替换该函数
• torch.Tensor.contiguous : 如果可以就返回自身，否则就返回副本
• torch.transpose : 如果可以就返回视图，否则就返回副本
• torch.reshape : 如果可以就返回视图，否则就返回副本
• torch.flatten : 如果可以就返回视图，否则就返回副本 (和 numpy.flatten 不同，后者总是返回副本)
• einops.repeat : 如果可以就返回视图，否则就返回副本
• einops.rearrange : 如果可以就返回视图，否则就返回副本
• 基础索引会返回视图，高级索引会返回副本.

存储对象

在一个 Tensor 上调用 storage() 会返回一个包装了底层 C++ 数组的 Python 对象。无论 Tensor 的维数是多少，这都是个一维数组。这能让你看到 Tensor 抽象前的内容并了解到实际数据在内存中的排布方式.

请注意每次调用 storage() 都会生成一个新的 python 包装对象，并且 x.storage() == x.storage() 和 x.storage() is x.storage() 结果均为 False.

如果需要检查两个 Tensor 是否共享底层的 C++ 数组，可以比较他们的 storage().data_ptr() 字段，这是他们底层 C++ 数组在内存中的指针。这对于调试很有用.

Tensor._base

如果 x 是一个视图，你可以用 x._base 访问他的原始 Tensor . 这是一个没有写在文档里的内部功能，了解一下很有用。考虑下面这段代码:

x = t.zeros(1024*1024*1024)
y = x[0]
del x

这里， y 是通过基础索引创建的，所以 y 是一个视图且 y._base 指向 x . 这意味着执行 del x 后系统不会释放 4GB 内存，该空间仍然会被继续使用，这个结果可能会非常反直觉。你可以用 y = x[0].clone() 进行替换，使用这个方法后允许你回收 x 的内存.

Exercise - 完成 raytrace_triangle

Difficulty: 🔴🔴🔴🔴⚪
Importance: 🔵🔵🔵🔵⚪
你应该花最多15-20分钟在这个练习上.
这个练习和intersect_rays_1d难度差不多,但我还是希望你能完成的更熟练.

你需要完成函数 raytrace_triangle , 功能是能检测 ray 中的每条射线是否和一个给定的三角形相交.

def raytrace_triangle(
    rays: Float[Tensor, "nrays rayPoints=2 dims=3"],
    triangle: Float[Tensor, "trianglePoints=3 dims=3"]
) -> Bool[Tensor, "nrays"]:
    '''
    For each ray, return True if the triangle intersects that ray.
    '''
    pass


A = t.tensor([1, 0.0, -0.5])
B = t.tensor([1, -0.5, 0.0])
C = t.tensor([1, 0.5, 0.5])
num_pixels_y = num_pixels_z = 15
y_limit = z_limit = 0.5

# Plot triangle & rays
test_triangle = t.stack([A, B, C], dim=0)
rays2d = make_rays_2d(num_pixels_y, num_pixels_z, y_limit, z_limit)
triangle_lines = t.stack([A, B, C, A, B, C], dim=0).reshape(-1, 2, 3)
render_lines_with_plotly(rays2d, triangle_lines)

# Calculate and display intersections
intersects = raytrace_triangle(rays2d, test_triangle)
img = intersects.reshape(num_pixels_y, num_pixels_z).int()
imshow(img, origin="lower", width=600, title="Triangle (as intersected by rays)")

Solution

def raytrace_triangle(
    rays: Float[Tensor, "nrays rayPoints=2 dims=3"],
    triangle: Float[Tensor, "trianglePoints=3 dims=3"]
) -> Bool[Tensor, "nrays"]:
    '''
    For each ray, return True if the triangle intersects that ray.
    '''
    # SOLUTION
    NR = rays.size(0)

    # Triangle is [[Ax, Ay, Az], [Bx, By, Bz], [Cx, Cy, Cz]]
    A, B, C = einops.repeat(triangle, "pts dims -> pts NR dims", NR=NR)
    assert A.shape == (NR, 3)

    # Each element of `rays` is [[Ox, Oy, Oz], [Dx, Dy, Dz]]
    O, D = rays.unbind(dim=1)
    assert O.shape == (NR, 3)

    # Define matrix on left hand side of equation
    mat: Float[Tensor, "NR 3 3"] = t.stack([- D, B - A, C - A], dim=-1)

    # Get boolean of where matrix is singular, and replace it with the identity in these positions
    # Note - this works because mat[is_singular] has shape (NR_where_singular, 3, 3), so we
    # can broadcast the identity matrix to that shape.
    dets: Float[Tensor, "NR"] = t.linalg.det(mat)
    is_singular = dets.abs() < 1e-8
    mat[is_singular] = t.eye(3)

    # Define vector on the right hand side of equation
    vec = O - A

    # Solve eqns
    sol: Float[Tensor, "NR 3"] = t.linalg.solve(mat, vec)
    s, u, v = sol.unbind(dim=-1)

    # Return boolean of (matrix is nonsingular, and solution is in correct range implying intersection)
    return ((u >= 0) & (v >= 0) & (u + v <= 1) & ~is_singular)

调试 (Debugging)

调试代码是一件非常重要的事。就像用 GPT 辅助代码一样，他可以显著加快你的开发速度，让你不再 bug 上浪费过多的时间.

为了让你练习使用 VSCode 的内置调试器进行调试，下面我们提供了一个示例函数。这是 raytrace_triangle 的一个实现，其中有一些错误。你的任务是使用调试器找到错误并修复它 (当然我知道我们上面已经提供了 solution 你应该也已经看过了，但是这个部分是为了训练你使用 VSCode 内置调制器的能力，所以请忘了参考答案！).

有错误的函数

def raytrace_triangle_with_bug(
    rays: Float[Tensor, "nrays rayPoints=2 dims=3"],
    triangle: Float[Tensor, "trianglePoints=3 dims=3"]
) -> Bool[Tensor, "nrays"]:
    '''
    For each ray, return True if the triangle intersects that ray.
    '''
    NR = rays.size[0]

    A, B, C = einops.repeat(triangle, "pts dims -> pts NR dims", NR=NR)

    O, D = rays.unbind(-1)

    mat = t.stack([- D, B - A, C - A])

    dets = t.linalg.det(mat)
    is_singular = dets.abs() < 1e-8
    mat[is_singular] = t.eye(3)

    vec = O - A

    sol = t.linalg.solve(mat, vec)
    s, u, v = sol.unbind(dim=-1)

    return ((u >= 0) & (v >= 0) & (u + v <= 1) & ~is_singular)


intersects = raytrace_triangle_with_bug(rays2d, test_triangle)
img = intersects.reshape(num_pixels_y, num_pixels_z).int()
imshow(img, origin="lower", width=600, title="Triangle (as intersected by rays)")    

你可以通过点击单元格底部的 **Debug cell (调试单元格)** 来调试。你的单元格应该包含实际运行时导致错误的代码 (而非包含错误来源的函数). 在运行调试器之前，你可以通过单击行号左侧 (单机后将出现一个红点) 来设置断点。然后，你可以使用调试器运行时出现的按钮工具栏单步调试代码 (参阅这里以了解每个按钮的功能说明). 当程序运行到断点时，你可以使用以下工具:

• 在左侧边栏 VARIABLES (变量) 窗口中检查局部和全局变量.
• 在左侧边栏 WATCH (观察) 窗口中添加要观察的变量表达式.
  • 你可以在这里输入任何表达式，例如变量的类型或列表的长度，这将在代码单步运行之后实时更新.
• 通过将表达式输入到 DEBUG CONSOLE (调试控制台) (出现在屏幕底部) 来一次性计算表达式的值.

请注意，你的代码会在你打断点的那一行开始执行之前停止执行。因此如果你在某一行上出现报错，你需要做的就是直接在这一行打断点.

如果你在 VSCode 中使用 jupyter notebook, 那上面这些内容的基本工作原理都差不多，除了一些小更改，例如调试按钮在单元格左上角的下拉菜单中 (如果你找不到，那么你需要进入用户设置并添加一行 "notebook.consolidatedRunButton": true ).

我们还想讨论有关调试的更多细节，但是这些已经足够满足大多数需求。调试器通常是比 print 或者 assert 更有效的调试方式 (尽管这两个在某些情况下也很有帮助).

Answer - 这些bug是什么和如何修复这些bug.

NR = rays.size[0]

这里应该是rays.size(0)(或者等价的rays.shape[0])size是一个类方法,需要接受一个整数作为参数并返回这个维度的形状;shape是一个实例属性可以接受索引.
这个问题就算没有调试器也能很容易的解决,因为报错信息非常详细.

O, D = rays.unbind(-1)

我们在分解错误的维度.rays的形状是(nrays, points=2, dims=3),而我们实际上想沿着points维度分解.所以我们应该用rays.unbind(1).
我们可以通过在变量窗口观察rays实例发现这一点(你可以点击变量名的下拉箭头来观察变量的属性值,包括shape),或者你可以通过在调试控制台输入rays.shape来查看属性值.这个错误应该也是显而易见的(使用类型检查的好例子!).

mat = t.stack([- D, B - A, C - A])

这里需要在最后多一个参数dim=-1,因为torch.stack默认沿着第一个维度堆积张量.
这个错误可能比较难发现,因为报错的行在实际出现错误的行前一行.当然,我们也可以通过在变量窗口检查mat张量的形状来发现这个错误.

---

这些都是相对容易发现的错误(并非所有错误都会在代码中报错,有些可能只是意外的结果).但希望这个练习能让你了解如何使用调试器.这是一个非常强大的工具,可以节约你很多时间.

加载 Mesh

使用提供的代码来加载三角形组成的皮卡丘。按照惯例，使用 torch.save 写入的文件都以 .pt 结尾，但实际上这些只是 zip 文件罢了.

with open(section_dir / "pikachu.pt", "rb") as f:
    triangles = t.load(f)

渲染 Mesh

对于我们的目标，mesh 就是一组三角形，因此为了渲染它，我们将同时使所有光线和所有三角形相交。我们之前只是返回一个 bool 值判断给定的射线是否与三角形相交，但现在可能有多个三角形与给定的射线相交.

对于每一条射线 (像素), 如果可以，我们将返回一个表示到三角形最短距离的浮点值，否则返回表示无穷大的特殊值 float('inf') . 我们现在不会返回哪些三角形是相交的.

请注意，到三角形的距离特指沿 x 轴的距离，而不是欧几里得距离.

Exercise - 完成 raytrace_mesh

Difficulty: 🔴🔴🔴⚪⚪
Importance: 🔵🔵🔵🔵⚪
你应该花最多20-25分钟在这个练习上.
这是我们一直在构建的主要函数,并且完成他标志着这一节核心内容的完成.他涉及了大量上一个练习中应该重复利用的代码.

和先前一样完成 raytrace_mesh , 重塑并可视化输出。你的皮卡丘以 (0,0,0) 为中心点，所以你应该把你的光线出发点移到旁边，至少 x=-2 来完整的观察他.

请记住， t.linalg.solve (和大部分批操作) 可以接受以批处理的方式接受多个维度。先前你仅仅用 NR (the number of rays, 光线的数量) 表示批维度，但是现在你也可以使用 (NR, NT) (the number of rays and triangles, 光线和三角形的数量) 作为你的批维度，所以你可以一口气 solve 所有的光线和三角形.

def raytrace_mesh(
    rays: Float[Tensor, "nrays rayPoints=2 dims=3"],
    triangles: Float[Tensor, "ntriangles trianglePoints=3 dims=3"]
) -> Float[Tensor, "nrays"]:
    '''
    For each ray, return the distance to the closest intersecting triangle, or infinity.
    '''
    pass


num_pixels_y = 120
num_pixels_z = 120
y_limit = z_limit = 1

rays = make_rays_2d(num_pixels_y, num_pixels_z, y_limit, z_limit)
rays[:, 0] = t.tensor([-2, 0.0, 0.0])
dists = raytrace_mesh(rays, triangles)
intersects = t.isfinite(dists).view(num_pixels_y, num_pixels_z)
dists_square = dists.view(num_pixels_y, num_pixels_z)
img = t.stack([intersects, dists_square], dim=0)

fig = px.imshow(img, facet_col=0, origin="lower", color_continuous_scale="magma", width=1000)
fig.update_layout(coloraxis_showscale=False)
for i, text in enumerate(["Intersects", "Distance"]): 
    fig.layout.annotations[i]['text'] = text
fig.show()

Solution

def raytrace_mesh(
    rays: Float[Tensor, "nrays rayPoints=2 dims=3"],
    triangles: Float[Tensor, "ntriangles trianglePoints=3 dims=3"]
) -> Float[Tensor, "nrays"]:
    '''
    For each ray, return the distance to the closest intersecting triangle, or infinity.
    '''
    # SOLUTION
    NR = rays.size(0)
    NT = triangles.size(0)

    # Each triangle is [[Ax, Ay, Az], [Bx, By, Bz], [Cx, Cy, Cz]]
    triangles = einops.repeat(triangles, "NT pts dims -> pts NR NT dims", NR=NR)
    A, B, C = triangles
    assert A.shape == (NR, NT, 3)

    # Each ray is [[Ox, Oy, Oz], [Dx, Dy, Dz]]
    rays = einops.repeat(rays, "NR pts dims -> pts NR NT dims", NT=NT)
    O, D = rays
    assert O.shape == (NR, NT, 3)

    # Define matrix on left hand side of equation
    mat: Float[Tensor, "NR NT 3 3"] = t.stack([- D, B - A, C - A], dim=-1)
    # Get boolean of where matrix is singular, and replace it with the identity in these positions
    dets: Float[Tensor, "NR NT"] = t.linalg.det(mat)
    is_singular = dets.abs() < 1e-8
    mat[is_singular] = t.eye(3)

    # Define vector on the right hand side of equation
    vec: Float[Tensor, "NR NT 3"] = O - A

    # Solve eqns (note, s is the distance along ray)
    sol: Float[Tensor, "NR NT 3"] = t.linalg.solve(mat, vec)
    s, u, v = sol.unbind(-1)

    # Get boolean of intersects, and use it to set distance to infinity wherever there is no intersection
    intersects = ((u >= 0) & (v >= 0) & (u + v <= 1) & ~is_singular)
    s[~intersects] = float("inf") # t.inf

    # Get the minimum distance (over all triangles) for each ray
    return s.min(dim=-1).values