简介:面向光学仿真学习者的MATLAB光线追迹程序包,聚焦太阳能热发电、光伏(PV)系统与通用光学系统的光路分析,适合需要快速开展几何光学仿真与二次开发的工程师和高校学生。压缩包共7个文件,以5个m脚本为主体,涵盖主程序main.m、反射/折射计算reflect.m与refract.m、交点求解intersection.m等核心模块,另附README说明与LICENSE许可文件;整套资源仅6KB,结构精简便于阅读。已有565人学习。程序利用Chebyshev多项式对光路进行近似,允许用户设置材料折射率、反射率、吸收率等参数,并通过MATLAB绘图展示光线路径与强度分布,便于直观理解聚光集热和光伏面板受光过程。对于学者而言,既可作为光线追迹入门的学习样板,也可在此基础上调整模型,适应不同光学设计场景,兼具教学与实用价值。 很多人一听到“光线追迹”四个字,第一反应就是“这玩意很硬核,非图形学大佬玩不转”。我第一次动手写的时候也差点被劝退,但当我真正把第一版程序跑通、看到屏幕上出现几个带有柔和渐变的球体时,才发现核心逻辑远比想象中简单。这篇文章是我计划里“光线追迹·相关程序”系列的第一篇,目标很明确:不用现成渲染引擎,从零手写一个最小可运行的光线追迹器,最终输出一张 PPM 格式的图片。适合想搞懂渲染原理、或者正在啃 Ray Tracing 系列教程但代码跑不起来的同学。
先说明白一件事:我这里说的“光线追迹”,指的是最基础的 Whitted 风格的光线追踪 / 路径追踪雏形。也就是说,我们从相机出发,对每个像素射出一条光线,光线打到球体后根据材质发生散射,散射后的光线继续和场景交互,最后把颜色写回像素。这一篇先不管性能,只追求逻辑闭环。
1. 先把核心逻辑拆开:为什么我们从相机“反向”射光线,而不是从光源射
如果你以前看过光线追迹的介绍,大概率听过类似的话:“渲染方程是对半球方向积分,理论上应该从光源出发追踪光子,但实际实现时为了节省计算量,我们反向追踪。” 这句话容易让人懵,我换个说法:光线从光源出发,经过无数次反射进入人眼,这个路径可以完全反过来走——从眼睛出发,一路采样,最后连到光源。物理上光路可逆,所以两种追踪方式结果一样,但反向追踪有一个巨大的好处:只计算那些真正进入相机的光线。正向追踪要模拟大量光子,绝大多数光子根本不会落进相机,纯属浪费算力。
这个程序的大致流程是固定的:
- 对输出图片的每一个像素,构造一条从相机位置出发、穿过该像素中心的光线;
- 让这条光线与场景中所有几何体求交,找出最近的交点;
- 在交点处根据材质决定继续散射的方向,生成新光线,重复求交;
- 当光线射到背景或达到最大弹射次数时停止,把累积的颜色写入像素。
你可以把相机想象成一只眼睛,像素想象成视网膜上的感光细胞。所谓“渲染”,其实就是逆推每一个感光细胞感知到的色彩是从哪个物体表面上反射过来的。这套逻辑的巧妙之处在于,它把全局光照问题拆解成了一条条可独立计算的光线路径,场景复杂了无非是增加“求交”的开销,但框架不变。
我建议新手先不要一上来就啃物理光学那一套,先把我上面的四步在纸上画出来,哪怕只是画一个球体、一条光线、一根法线,也比盲目抄代码强得多。
2. 从零写一个向量类和球体求交:问题全在判别式
有了流程,第一个要解决的就是数学工具。光线追迹的每一个计算几乎都离不开向量,所以先写一个最精简的vec3类。不需要 Numpy,因为我们要手动控制每一步,Python 原生类最直观。
import math class vec3: def __init__(self, x=0, y=0, z=0): self.x = x self.y = y self.z = z def __add__(self, o): return vec3(self.x + o.x, self.y + o.y, self.z + o.z) def __sub__(self, o): return vec3(self.x - o.x, self.y - o.y, self.z - o.z) def __mul__(self, s): return vec3(self.x * s, self.y * s, self.z * s) def dot(self, o): return self.x * o.x + self.y * o.y + self.z * o.z def cross(self, o): return vec3( self.y * o.z - self.z * o.y, self.z * o.x - self.x * o.z, self.x * o.y - self.y * o.x ) def norm(self): return math.sqrt(self.dot(self)) def normalize(self): n = self.norm() return vec3(self.x / n, self.y / n, self.z / n)然后是光线类。光线的数学定义很简单:经过点 (O)、方向为 (D) 的射线,可以写成 (P(t) = O + t \cdot D),其中 (t) 可以理解成光线走了多远。物体求交的本质就是求这个射线方程和物体表面方程的共同解。
class Ray: def __init__(self, origin, direction): self.origin = origin self.direction = direction.normalize() def at(self, t): return self.origin + self.direction * t这一篇的场景物体只做球体,因为球体求交有解析解,非常适合讲清原理。球体的隐式方程为 (|P - C|^2 = r^2),把 (P(t) = O + tD) 代进去:
[ |D|^2 t^2 + 2D \cdot (O - C) t + |O - C|^2 - r^2 = 0 ]
如果光线方向做了归一化,(|D|^2 = 1),这就是一个标准一元二次方程。判别式 (b^2 - 4ac) 决定光线是否穿过球体。
def hit_sphere(sphere, ray, t_min=0.001, t_max=1e9): oc = ray.origin - sphere.center a = ray.direction.dot(ray.direction) b = 2.0 * oc.dot(ray.direction) c = oc.dot(oc) - sphere.radius * sphere.radius disc = b * b - 4 * a * c if disc < 0: return None sqrt_d = math.sqrt(disc) # 先取较小的根,也就是光线进入球体的那一面 t = (-b - sqrt_d) / (2 * a) if t_min < t < t_max: return t # 如果较小的根不在有效区间,试较大的根 t = (-b + sqrt_d) / (2 * a) if t_min < t < t_max: return t return None这里要注意两个细节。第一,t_min不能设为 0,因为光线可能在某个交点处再次发射新光线,如果新光线的起点紧贴物体表面,浮点误差会算出微小的t,导致光线立刻和同一个物体相交,出现“自交黑斑”。第二,为什么先取较小根?因为在相机视角下,我们看到的球体永远是“进入点”那一侧,除非相机在球体内部才会看到“出点”。把所有交点收集起来取最近的那个,是后面场景管理和遮挡判断的基础。
别小看几十行代码,整个渲染器的“骨架”就在这里了。后面加场景管理、材质、相机,都是在为这个求交函数增加上下游环节。
3. 构建虚拟相机坐标系:让你的画面有“透视感”
求交解决了,但你得知道往哪个方向射光线。这个方向由虚拟相机决定。最简单的方式是把相机看成一个刚体,它有三个互相垂直的基向量:向右的 (u)、向上的 (v)、背对观察方向的 (w)。这套坐标系称为右手相机坐标系。
先给定lookfrom(相机位置)、lookat(观察目标)、up(世界坐标系中的上方向),构建方法如下:
w = (lookfrom - lookat).normalize() # 指向观察反方向 u = up.cross(w).normalize() # 右方向 v = w.cross(u) # 上方向然后用视场角(FOV)和输出图片宽高比确定视平面的大小。假设视平面距离相机为 1,那么视平面的垂直半高度为tan(fov / 2),水平半宽度为垂直半高度 * aspect。对于像素坐标(i, j),把它在宽高范围内归一化到[-1, 1],再乘以对应的半宽和半高,就得到光线的方向向量:
def ray_for_pixel(cam, i, j, width, height, fov): aspect = width / height px = (2 * ((i + 0.5) / width) - 1) * math.tan(fov / 2) * aspect py = (1 - 2 * ((j + 0.5) / height)) * math.tan(fov / 2) direction = (cam.u * px + cam.v * py - cam.w).normalize() return Ray(cam.position, direction)这里最容易被新手忽略的是py的符号。我们的像素坐标系是左上角为(0, 0),从上往下j递增;但世界坐标系的v方向是“向上”的。如果直接py = 2 * (j / height) - 1,出来的图像会是上下颠倒的。用1 - 2 * ((j + 0.5) / height)相当于把图像纵轴翻转,让第一行对应世界坐标系中最上面的扫描线。
构建好相机后,可以先做一个调试:不渲染任何物体,只输出背景颜色。背景色可以是基于光线方向y分量的蓝色渐变,比如天空到地平线的过渡梯度。这一步是用来确认每一行像素都对应正确的方向,避免整个画面的方向和预期不符。
相机模型带来的“透视感”也值得多说一句:为什么远处的物体会变小?因为同样大小的物体在视平面上占的角度更小,像素覆盖的面积更小。这个效果完全由透视投影的选择带来,和物体本身的物理尺寸无关。如果你希望做正交投影,只要把方向统一固定为相机朝向,不再根据像素坐标偏移即可,但那样就没有近大远小了。
4. 漫反射材质和伽马校正:虽然简单,但有无数人栽在颜色上
光线打中物体之后,下一步是决定它往哪里散射。物理真实的光照是 BRDF 积分,但这一篇先实现最简单的完全漫反射,也就是 Lambert 材质:入射光和出射方向无关,散射方向在交点法线所在的半球内均匀分布。这个模型虽然粗糙,但效果非常接近现实中未上釉的泥土、墙面、石膏表面。
实现原理很朴素:在交点处取法线 (N),把它看作球心,在单位球表面随机取一个点 (S),然后新的光线方向就是 (S) 到交点的方向。如果随机点落在法线下方,就把它翻到法线一侧。数学上等价于:随机单位方向向量d,如果d.dot(N) < 0,则取-d作为散射方向。这一步保证了散射只发生在表面外侧,光线不会穿进物体内部。
def random_in_hemisphere(normal): # 标准方法:先取单位球内的均匀点,再归一化 while True: p = vec3(random.uniform(-1, 1), random.uniform(-1, 1), random.uniform(-1, 1)) if p.dot(p) >= 1: continue return p.normalize() def scatter(ray, hit_point, normal): direction = random_in_hemisphere(normal) if direction.dot(normal) < 0: direction = direction * -1 return Ray(hit_point, direction)彩色输出我采用最简单的“光线能量衰减”策略:每次散射把颜色乘以一个小于 1 的衰减系数,比如albedo(物体表面反照率)。这样,光线弹射次数越多,能量越低,颜色越暗。一个典型的递归跟踪函数是这样的:
def trace(ray, depth=0): if depth > 50: return vec3(0, 0, 0) hit = scene_hit(ray) # 遍历所有球体,返回最近交点 if not hit: return background_color(ray.direction) scattered = scatter(ray, hit.point, hit.normal) # 0.5是漫反射反照率,可以换成球体自身的albedo return vec3(0.5, 0.5, 0.5) * trace(scattered, depth + 1)所有像素的颜色都需要多次随机采样然后取平均,这就是最简单的蒙特卡洛路径追踪雏形。每个像素采样几十次,噪点就会明显减少。采样数越高,画面越平滑,但耗时线性增加,所以调试时先开到最小分辨率、每像素采样 8 到 16 次就够了。
我现在要重点强调一个几乎所有第一次写渲染器的人都会踩的坑:伽马校正。屏幕和图片格式显示的颜色并不是线性光强,大部分显示器的输出亮度和像素值之间有一个约为 2.2 的幂次关系。我们计算出的颜色是“线性空间”的颜色,直接写进 PPM 图片再打开,整体会明显偏暗。解决办法是把颜色值做一次pow(color, 1/2.2),再乘 255 量化写入文件。不做这一步,你调出来的光照参数可能全都不准,这是最隐蔽的坑。
def write_color(f, color): r = max(0.0, min(1.0, color.x)) g = max(0.0, min(1.0, color.y)) b = max(0.0, min(1.0, color.z)) # 伽马校正 r = math.pow(r, 1.0 / 2.2) g = math.pow(g, 1.0 / 2.2) b = math.pow(b, 1.0 / 2.2) f.write(f"{int(r * 255)} {int(g * 255)} {int(b * 255)}\n")我用 PPM 格式做输出,因为它结构极简:文件头P3,然后是宽度、高度、最大像素值 255,后面每行三个整数表示一个像素的 RGB。几乎所有图像查看器都支持 PPM,Python 也能用 Pillow 轻松转 PNG。不要先急着学 HDR、EXR 那套,格式越简单越容易定位问题。
5. 第一版渲染器的实测效果与三个“隐形杀手”调试记录
如果真的照着上面的代码去跑,我保证你的第一张渲染图大概率不是完美的,而是会遇到这么几个问题。我说的不是网上理论上的“可能的坑”,而是我实际调试时踩过、并且花了不少时间才定位到原因的坑。
第一个坑:球体表面出现一圈黑色“焊点”或斑点。原因就是我前面提到的t_min设得不够大。光线在物体表面散射时,新的交点起点和原交点理论上是同一个点,但由于浮点误差,起点可能稍微陷入球体内部。如果t_min=0,第二条光线会在极短的t处再次与同一个球体相交,算出的颜色几乎全黑。解决方法是把t_min设为0.001甚至0.0001,给光线一个“安全距离”。这个小参数直接决定渲染画面干不干净。
第二个坑:画面并不是想象中的渐变,而是一团噪点,或者亮一块暗一块。这个问题八成出在随机单位向量的实现上。一种常见错误是直接用(x,y,z)在区间[-1,1]内均匀取点,然后只判断dot(p,p) < 1就归一化。这个写法的确可以生成单位球内的点,但我在初版代码里犯了个懒,没有做“方向在半球内”的判断,导致一半的散射光线往表面内部走,碰到球体内部又求交到奇怪的位置,最终结果就是颜色斑斓。用我上面给的random_in_hemisphere函数后,噪点会大幅度减少。
第三个坑:图像上下颠倒或者镜像。这个问题第 3 节提过,但实际发生的时候,很多人会先去怀疑相机坐标系建错了。我有一个快速验证方法:让相机稍稍向右偏离目标点,观察图像中球体的位置。如果球体应该在画面左侧,但渲染出来跑到右侧,那说明u方向反了;如果上下颠倒,则是py符号反了。这类方向问题用“单球体+背景渐变”调试最快,一次就能定位。
调试渲染器的时候,我强烈建议输出一个固定尺寸的小图,比如160x90,并且把递归深度限制从 50 改到 20。Python 纯循环做路径追踪非常慢,160 宽的小图,每像素 8 次采样,在我的电脑上大约要十几秒,这个量级足够我们频繁修改代码并观察结果。如果你一上来就跑 800px 的大图,每像素上百次采样,跑十分钟才出图,调试体验会非常糟糕。
还有一个和性能有关的细节:场景里只有一个球体时,trace函数每次弹射都要把所有球体遍历一遍。这个场景规模根本无所谓,但如果你按照教程加了 20 个球,每像素几十次采样,每帧几十次弹射,总的求交次数就会爆炸。不要急着上 BVH 或八叉树,先把输出分辨率降到160x90,把逻辑验证完,再考虑加速结构。
个人写渲染器时有个习惯:永远给程序加两个“自检开关”。第一个是输出每条光线的交点深度值,转成伪彩色,用来快速确认几何信息是否正确;第二个是输出随机种子,这样可以快速复现某个像素上的噪点。别觉得这是在浪费时间,调试光路程序时,能稳定复现一个 Bug 比什么都重要。
到目前为止,这个最小程序已经能画出带有柔和明暗变化的漫反射球体了。接下来可以继续加反射、折射、多物体、阴影光线、抗锯齿,这些内容不是往现有逻辑里硬塞,而是在“光线求交—散射”这条主链路上增加分支和递归策略。有了这一版作为地基,后面每一步都是水到渠成。
本文还有配套的精品资源,点击获取