从画出一束光开始：用 C 实现二维光线追踪与反射
这个项目的起点，是看到 Daniel Hirsch 的 《Ray Tracing in C》。一个光源、一个障碍物，再加上从光源向四周散开的光线，简单的画面就能把传播和遮挡表现得很直观。于是我也想自己写一个，在 AI 时代练习一下古法编程（笑）。
于是就有了这个小项目：
窗口和像素绘制交给 TIGR，光线如何前进、怎么判断碰撞、碰撞后往哪里走，由我自己实现。
在视频带来的基础思路上，我又加上了光线反射：光碰到窗口边缘会折返，碰到障碍物也会沿着新的方向继续传播。这里最值得展开的，是如何只靠碰撞位置附近的像素，估计出表面法向量。
demo效果：
当前程序里，按住鼠标左键可以移动光源，右键可以移动障碍物。默认发射 360 条光线，并绘制一次反射后的路径。这是一个二维光路演示，重点是让传播、碰撞和反射看得见；三维场景、材质着色和全局光照并没有涉及。
1. 程序流程
项目的主体在 ray_tracing.c 中。数据结构也很直接：Circle 保存光源和圆形障碍物的位置、半径及颜色，Vec2 表示二维向量，Mouse_Event 保存鼠标状态。
TIGR 轻量图形库提供创建窗口、读写像素、绘制圆形和读取鼠标输入的接口。每一帧，我先清空画面，读取鼠标并更新障碍物与光源的位置，再画障碍物和光源，最后计算并绘制所有光线并刷新窗口。
这里的绘制顺序有实际意义：当前版本直接读取画布颜色来识别障碍物。 is_obstacle() 用 tigrGet() 取得像素，再检查 RGB 是否等于 SILVER。因此，必须先把障碍物画出来，光线才有东西可碰。
几个函数基本就能串起这条流程：
| 函数 | 做什么 |
|---|---|
Process_Mouse_Evevt() |
根据鼠标位置移动光源或障碍物 |
Generate_Rays() |
发射光线，步进绘制，处理碰撞和反射 |
is_obstacle() |
根据像素颜色判断是否碰到障碍物 |
is_boundary() |
判断某个障碍物像素是否位于边界 |
Normal() |
采样边界像素，用 PCA 估计法向量 |
这种做法让场景数据非常简单，也方便从零搭起来。不过，显示和碰撞共用一张画布，会让算法依赖绘制顺序和颜色。后面如果增加材质、叠加效果或更多物体，应该可以通过障碍物掩码的方式来实现（虽然大概率不会修改了，毕竟是练习项目hhh）。
2. 绘制一条光线：其实就是不断往前走的点
先不考虑反射。给定起点和方向，只要沿方向不断移动，再把经过的位置画出来，就得到一条光线。
2.1 均匀获取光线发射角
设光线数量为 N，对点光源来说，相邻两条光线的角度差为：
\Delta\alpha=\frac{2\pi}{N}
第 i 条光线的角度是 \alpha_i=i\Delta\alpha。代码中使用：
double rad = i * rad_step;
double x_step = sin(rad);
double y_step = cos(rad);
这里把角度零点放在了向上的方向，所以是 sin 对应 x、cos 对应 y。它们组成的方向向量长度为 1。
还要记住一个贯穿全文的约定：方向向量的 y 向上为正，屏幕坐标的 y 向下为正。 因此，把方向应用到屏幕位置时，y 分量要用减法。
2.2 从光源表面出发，按 step 步进
光源在画面中是一个小圆。光线的起点取在圆周上：
photon_x = bulb.x + bulb.r * x_step;
photon_y = bulb.y - bulb.r * y_step;
随后重复下面的操作。为突出顺序，这里只保留绘制、推进和取像素坐标的部分：
tigrPlot(screen, (int)photon_x, (int)photon_y, ray_color);
photon_x += x_step;
photon_y -= y_step;
pixel_x = (int)photon_x;
pixel_y = (int)photon_y;
// 接下来检查新位置是否碰到边框或障碍物。
如果显式写出步长 h，更新公式就是：
x_{k+1}=x_k+h d_x,\qquad y_{k+1}=y_k-h d_y
当前实现相当于 h=1，每次在连续坐标中前进一个像素单位。photon_x、photon_y 使用 double 保存，所以小数不会随着每次绘制被丢掉；只有画点和查询像素时才转成 int。在屏幕内的非负坐标上，这相当于取所在的整数网格位置。
“前进一个单位”和“每次恰好画一个新的相邻像素”并不完全相同：斜线可能连续落入同一个像素，画出来也会有阶梯感。这里用固定步长来换取容易理解的实现。步长增大可能漏过细小障碍物，减小则会增加查询次数。
3. 加上反射：先找到表面朝向
边框的反射比较容易。左右边框是竖直的，撞到后让 x_step 变号；上下边框是水平的，让 y_step 变号。另一个分量保持不变。
圆形障碍物也有一个直接解法：从圆心指向碰撞点，再归一化，就是法向量。代码中保留了 Normal_circle()，但真正参与反射计算的是 Normal() 返回的 PCA 法向量。
我想尝试的是：如果只知道哪些像素属于障碍物，能不能从它的轮廓推断表面朝向？这样，计算过程就不必依赖圆心和半径。不过这仍然是局部估计，在尖角或多个轮廓混进同一个窗口时，不一定得到理想结果。
3.1 先在碰撞点附近找边界像素
并不是把窗口里的所有障碍物像素都拿来计算。我要找的是轮廓：一个像素本身属于障碍物，而且上下左右四个邻居中至少有一个不属于障碍物，它就是边界像素。
is_boundary() 的逻辑：
如果当前像素不是障碍物：返回否
如果左、右、上、下任一邻居不是障碍物：返回是
否则：返回否
以碰撞像素为中心，在 x、y 两个方向各取 SAMPLE_RADIUS 个像素。当前半径是 10，循环包含两端，因此采样范围是一个 21 × 21 的正方形，并不是圆形窗口。
Normal() 把这些边界点放入 b_points 数组。接下来要做的事情，可以先理解成：给这小段像素轮廓找一条最贴近它的直线。这条直线的方向近似于切线，垂直于它的方向就是法线。
3.2 PCA 为什么能找到切线
PCA 是 Principal Component Analysis，中文叫主成分分析。这里用到的是一个很具体的性质：一组点沿哪个方向最分散，第一主成分就指向哪个方向。对局部近似为直线的边界点，沿轮廓方向的分布更长，垂直方向的分布更窄，因此第一主成分可用来估计切线。关于“最大方差方向”的一般解释，可以参考 CMU 的 PCA 讲义。
下面把这个想法写成代码对应的计算。
假设采样得到 m 个边界像素 (x_i,y_i)。先求均值，也就是点云的中心：
\bar{x}=\frac{1}{m}\sum_i x_i,\qquad
\bar{y}=\frac{1}{m}\sum_i y_i
再把每个点减去中心，记作 u_i=x_i-\bar{x}、v_i=y_i-\bar{y}，累加三个量：
A=\sum_i u_i^2,\qquad B=\sum_i u_iv_i,\qquad D=\sum_i v_i^2
于是得到对称矩阵：
M=\begin{bmatrix}A&B\\B&D\end{bmatrix}
在源码里，它们分别叫 xx、xy 和 yy。严格说，代码算的是未归一化的散布矩阵；除以样本数后才得到一种常用定义下的协方差矩阵。整体乘一个正数不会改变特征向量，所以这里只求方向，可以省去这一步。
取一个候选方向 t=(\cos\theta,\sin\theta)，把中心化后的点投影到它上面，投影平方和是：
S(\theta)=t^TMt
=A\cos^2\theta+2B\sin\theta\cos\theta+D\sin^2\theta
使用二倍角公式整理：
S(\theta)=\frac{A+D}{2}
+\frac{A-D}{2}\cos(2\theta)+B\sin(2\theta)
当后面两项取得最大值时，就找到了点云最分散的方向。在方向不退化时，可以写成：
\theta=\frac{1}{2}\operatorname{atan2}(2B,A-D)
这就是源码里这行公式的来历：
theta = 0.5 * atan2(2.0 * xy, xx - yy);
它给出的是切线方向角，还不是法向量。由于输入的是屏幕像素坐标，屏幕坐标下的切向量是 (\cos\theta,\sin\theta)，将其旋转 90°，可得屏幕法向量 (-\sin\theta,\cos\theta)。
再把屏幕 y 转换成向上为正的方向坐标，就得到源码使用的结果：
normal.x = -sin(theta);
normal.y = -cos(theta);
这个向量天然是单位向量。这里如果漏掉 y 的符号转换，后续反射就会使用不一致的坐标系。
3.3 有了法向量，计算反射光线
设入射单位方向为 d，单位法向量为 n。把 d 拆成沿表面和垂直表面的两部分：
d_n=(d\cdot n)n,\qquad d_t=d-d_n
镜面反射时，沿表面的分量保留，垂直表面的分量反向，所以：
r=d_t-d_n=d-2(d\cdot n)n
对应的 C 代码只有几行：
double dot = x_step * normal.x + y_step * normal.y;
x_step = x_step - 2.0 * dot * normal.x;
y_step = y_step - 2.0 * dot * normal.y;
例如光线以 d=(0.6,-0.8) 撞向水平表面，取法向量 n=(0,1)，点积是 -0.8，反射方向就是 (0.6,0.8)：水平分量没变，竖直分量翻转，和前面的边框处理一致。
PCA 算出的法线可能朝向两侧中的任意一侧。代码在 dot > 0 时同时翻转法向量和点积，让法线迎着入射方向。这个判断来自点积的几何意义，与屏幕 y 轴本身无关；而对上述反射公式，n 和 -n 实际上会得到相同结果。
3.4 获取反射方向后，更改光线步进 step
改变方向后，代码把光线放到碰撞像素沿反射方向前进一个单位的位置：
photon_x = pixel_x + x_step;
photon_y = pixel_y - y_step;
这是为了让光线离开碰撞位置，减少马上再次命中同一个像素的情况。它是简单的偏移处理，并不能保证在所有曲率和入射角下都避开自相交。
每次反射还会把颜色的 alpha 乘以 0.75，让反射段看起来更暗。这里控制的是绘制不透明度，用来表现衰减，没有进行严格的光能计算。
REF_NUM 默认是 1。源码循环条件为 ref_cont <= ref_num，因此会绘制初始段和第一次反射后的光路；在下一次碰撞使计数超过限制后退出，不再画下一段。调整参数时，要把“处理了碰撞”和“绘制了反射段”区分开。
3.5 采样窗口不是越大越好
窗口小，得到的法向量更局部，但容易受像素阶梯影响；窗口大，点更多，方向通常更平滑，却也可能把弯曲轮廓甚至另一段边界混在一起。
这套方法隐含了一个前提：窗口中的轮廓能够近似成一条直线。尖角附近可能同时出现两个方向，PCA 会给出一个折中的方向，并不能自动识别应该反射在哪一条边上。
当前实现也还有一些值得继续补的地方：没有处理零采样点、点云退化或两个主方向同样显著的情况。默认 21 × 21 窗口最多包含 441 个像素，能放进 b_points[512]；但增大采样半径时，就需要同时检查数组容量，不能只改宏定义。
4. 为什么我觉得它适合拿来练 C
我喜欢这个项目的一点，是代码里的数值很快就能变成画面：改光线数量，会看到疏密变化；改法向量，反射方向就跟着变；改采样范围，则能观察到局部轮廓估计的差别。
它也把几个平时容易分开练习的内容放在了一起。结构体用来描述场景，指针用来更新对象状态，数组保存采样点，循环负责逐条、逐步推进光线，sin、cos、atan2 则把几何关系落成具体计算。每一部分都有明确的用途。
图形库方面，我用 TIGR 创建窗口，通过 tigrPlot() 画点、tigrGet() 读像素、tigrMouse() 读取输入，再用 tigrUpdate() 更新画面。这样能练到图形程序的基本循环，又不用一开始就处理底层窗口系统，让练习的范围比较合适。
同时了解一下 CMake
项目的 CMakeLists.txt 将 ray_tracing.c 和 tigr/tigr.c 一起编译为 raytracer，要求 C99，并设置头文件搜索路径。平台相关依赖则放在条件分支里：Windows 链接 opengl32、gdi32，macOS 链接 OpenGL、Cocoa 框架，Linux 分支链接 GLU、GL、X11 和数学库 m。
从这个例子能看出，构建一个 C 程序除了编译源文件，还要解决头文件搜索和链接依赖。CMake 把这些关系记录下来，再交给具体构建工具执行。
准备好 C 编译器、CMake 3.10 或以上版本以及对应平台依赖后，具体按照仓库的 中文 README编译即可。
还有个实际会遇到的小细节：仓库内置 TIGR 的 Windows 和 Linux/X11 后端，对鼠标右键使用的位值不同。当前源码按 Windows 的 4 判断，Linux/X11 对应 2，移植运行时需要相应调整。
源码放在 GitHub：ray-tracing-c。如果你也刚学完 C 的基础语法，可以从只画一条光线开始，接着做遮挡，再让它在边框上弹回来。等光线真的沿着自己算出的方向反射时，那几行向量运算就不再只是公式了。