Radon变换正向投影详解：从公式到C++/Matlab代码实现

Radon变换这个词，CT相关论文里几乎每篇都要出现，但真能把它从公式变成能跑的代码，并且代码结果还跟MATLAB自带radon函数对得上的人，其实不多。这篇是我CT断层成像系列的第三篇，专门讲Radon变换的正向投影是怎么实现的，也就是给定一张图像，怎么算出来它在不同角度下的投影数据，也就是俗称的sinogram，同时给出一套C++和一套Matlab可运行的代码。想做CT仿真数据生成、图像重建入门，或者准备笔试面试想快速过一遍Radon实现细节的同学，这篇应该能直接帮你省掉踩坑的时间。阅读前你只需要懂点积分、会写循环就行，剩下的我尽量用人话讲明白。

1. 正向投影到底在算什么：从Beer-Lambert定律到Radon公式

1.1 X射线投影的物理过程与数学表达

CT成像的第一步，不是重建，而是采集投影。X射线穿过人体时强度会衰减，衰减规律满足Beer-Lambert定律：I = I0 · e^(−∫ μ dl)。其中μ是物体对X射线的线性衰减系数，积分沿着射线路径走。探测器收到的不是这个积分本身，而是衰减后的强度I，所以实际使用中要对原始信号取对数：

g = ln(I0 / I) = ∫ μ dl

这个g就是我们常说的投影数据。也就是说，CT图像上的每一个像素值对应一个μ值，CT扫描得到的每个投影值，是沿着某条直线对μ的积分。从数学角度看，这个把二维函数f(x,y)变成一维积分集合的操作，就是Radon变换：

p(s, θ) = ∬ f(x,y) · δ(x·cosθ + y·sinθ − s) dxdy

你可以把它理解成：固定一个角度θ，用无数条平行线去扫这张图像，每条直线上的所有像素一起"叠加"成一个数值。不同的θ就得到不同角度的投影，把所有角度的投影按行堆积在一起，就是CT里最常见的二维数据图——sinogram。Sinogram的横轴是探测器位置s，纵轴是投影角度θ。

这里有一个初学者容易忽略的性质：平行束扫描时，角度覆盖180°就够了，不需要扫360°。因为p(s, θ+π)和p(−s, θ)描述的是同一条射线，只是方向相反，信息完全重复。所以后面所有代码里，角度范围都取0°~179°或者0°~180°，这是一条非常重要的约定。

1.2 离散化前必须决定的三个关键参数

连续公式看着简洁，但电脑里只有离散像素，所有变量都得换算成网格坐标。做正向投影之前，有三个参数必须拍板：

• 角度采样数量：决定sinogram有多少行。临床CT通常几千个角度，教学和仿真一般用180或360个，每度一个投影。角度间隔越小，投影数据越稠密，重建质量越高，但计算量和存储都会涨。
• 探测器单元数量：决定sinogram有多少列。这里有个非常容易被坑的细节：探测器数量不是随便取图像宽度N就行。因为当角度旋转到45°时，图像四角在投影方向上的最大距离会达到大约N·√2 / 2，而图像中心到边缘的距离只有N/2。如果你只准备N个探测器单元，大角度下靠近四个角的像素会直接飞出探测器覆盖范围，投影能量凭空丢失。所以稳妥做法是取 detector_count = ceil(N·√2)，再向上取一个偶数。
• 探测器间距与像素间距的比例：一般情况下我都设成1.0，也就是探测器单元尺寸和像素尺寸一致。这个比例不一定要是1，但保持1会让坐标换算最直观，也方便和MATLAB内置radon结果对照。

这三个参数定好，剩下就是坐标映射的事了。

1.3 像素驱动与射线驱动：两种实现思路的取舍

实现Radon正投影有两条路：像素驱动（pixel-driven）和射线驱动（ray-driven）。

像素驱动是数学公式最直接的翻译：遍历每一个像素，用当前像素中心坐标(x, y)算出它在探测器上的投影位置 s = x·cosθ + y·sinθ，然后把像素值累加到对应的探测器单元上。这个思路写起来简单，内层循环只有加减乘除，性能主要取决于图像尺寸。

射线驱动则更贴近物理：遍历每一条射线，沿射线与图像的交线段去采样或积分，把射线上经过的所有像素按长度权重加总。这个方法需要用到求线段与矩形网格交点的算法，比如Liang-Barsky，实现复杂一些，但更精确，也更容易扩展成扇形束和锥形束模型。

两种方法结果在极限情况下等价。区别是：像素驱动速度快、实现容易，适合反向投影时的对称利用；射线驱动精度高、物理一致性更好，工业CT和蒙特卡洛仿真领域用得更多。我下面的代码全部用像素驱动加线性插值，理由很简单：代码短、逻辑清楚、便于对照公式。射线驱动看完这篇理解了坐标系统之后，自己也能写出来。

2. C++实现：像素驱动法带你一步步调通

2.1 坐标约定与探测器索引映射：所有偏移的根源

写代码前，先花一分钟把坐标系钉死。图像在电脑里是二维数组，行下标i从0到rows−1，列下标j从0到cols−1，第0行在显示上是顶部。我约定物理坐标x向右为正，y向上为正，原点在图像中心。于是：

x = j − cols/2 + 0.5
y = rows/2 − i − 0.5

为什么要加0.5？因为像素值其实是一个小方格内的均值，我们通常把这个方格的值浓缩到中心点。第0列的像素中心并不是x = −cols/2，而是x = −cols/2 + 0.5。这个半像素偏移虽然小，但在和内置radon做数值对比时会暴露出来，千万不能省。

有了物理坐标，任一像素在角度θ下的投影位置就是：

s = x·cosθ + y·sinθ

接下来把连续位置s换算成离散探测器索引。探测器阵列间距为1，物理覆盖范围是[−detector_count/2, detector_count/2]，所以：

s_index = s / spacing + detector_count / 2

在C++数组中，s_index是浮点数，直接取整就会把像素值一股脑塞进一个单元，产生锯齿。更稳的做法是用线性插值：找到s_index前后两个整数索引，按距离倒数分配权重，一个像素最多贡献给两个探测器单元。这样投影数据平滑，也更接近真实测量。

2.2 核心函数实现：线性插值像素驱动法

下面这段代码我完整跑过，拿64×64的仿真图像测，角度0°、45°、90°的结果都能对上手工计算。完整工程包含一个forwardProjectPixelDriven函数和一个demo用的main函数。

#include <iostream> #include <vector> #include <cmath> using Image = std::vector<std::vector<double>>; using Sinogram = std::vector<std::vector<double>>; constexpr double kPi = 3.14159265358979323846; struct Geometry { int detector_count; double detector_spacing; std::vector<double> angles_deg; }; Sinogram forwardProjectPixelDriven(const Image& img, const Geometry& geo) { int rows = img.size(); int cols = img[0].size(); int n_angles = geo.angles_deg.size(); int n_det = geo.detector_count; Sinogram sino(n_angles, std::vector<double>(n_det, 0.0)); double half_cols = cols / 2.0; double half_rows = rows / 2.0; double half_det = n_det / 2.0; for (int a = 0; a < n_angles; ++a) { double theta = geo.angles_deg[a] * kPi / 180.0; double ct = cos(theta); double st = sin(theta); for (int i = 0; i < rows; ++i) { double y = half_rows - i - 0.5; for (int j = 0; j < cols; ++j) { double x = j - half_cols + 0.5; double s = x * ct + y * st; double s_idx = s / geo.detector_spacing + half_det; int idx0 = (int)floor(s_idx); double frac = s_idx - idx0; if (idx0 >= 0 && idx0 < n_det) { sino[a][idx0] += img[i][j] * (1.0 - frac); } if (idx0 + 1 >= 0 && idx0 + 1 < n_det) { sino[a][idx0 + 1] += img[i][j] * frac; } } } } return sino; } int main() { const int N = 64; Image img(N, std::vector<double>(N, 0.0)); // 构造测试图：中心一个 20x20 的亮块 for (int i = 0; i < N; ++i) { for (int j = 0; j < N; ++j) { if (i >= 22 && i < 42 && j >= 22 && j < 42) { img[i][j] = 1.0; } } } Geometry geo; geo.detector_count = 2 * (int)ceil(N * sqrt(2.0) / 2.0); // 向上取偶数 geo.detector_spacing = 1.0; geo.angles_deg = {0.0, 45.0, 90.0, 135.0}; auto sino = forwardProjectPixelDriven(img, geo); for (int a = 0; a < (int)geo.angles_deg.size(); ++a) { std::cout << "angle " << geo.angles_deg[a] << ": "; for (int d = 0; d < geo.detector_count; ++d) { std::cout << sino[a][d] << " "; } std::cout << "\n"; } return 0; }

代码逻辑概括起来就是三重循环：最外层遍历角度，中间层遍历像素行，最内层遍历像素列。每算出一个s，就用线性插值把img[i][j]按权重分给两个探测器单元。注意第二个if的边界判断一定要写成idx0 + 1 < n_det，如果你只写idx0 < n_det就去掉一半权重，总能量会悄悄减少。

2.3 编译运行与性能优化：O2、OpenMP和小图调试

编译命令很简单：

g++ -O2 forward_radon.cpp -o forward_radon

如果设备支持OpenMP，可以在角度循环前加一行#pragma omp parallel for，因为每个角度写的是sino里不同的一行，互不干扰，可以直接并行。编译命令改成：

g++ -O2 -fopenmp forward_radon.cpp -o forward_radon

实测下来，128×128图像180个角度，单线程大概一秒钟内跑完；512×512图像180个角度，单线程需要几十秒量级，开O2加OpenMP能压进几秒。如果你想加速又不引入OpenMP依赖，还有一个笨办法：把图像缩小到64×64调试逻辑，确认输出正确后再放大，不要在512×512上第一次跑就急着去验证数学推导，那会浪费大量时间。

调试阶段建议打印sinogram每一行的总能量，也就是所有探测器值求和。Radon变换本质是线积分，理论上每个角度下的投影总能量应该等于图像全部像素值之和。如果45°时总能量明显小于0°的总能量，多半是探测器数量不够导致边缘截断，而不是代码逻辑错误。这个自检方法我一直觉得比盯着数字看要高效率很多。

3. Matlab实现：从教学循环到向量化

3.1 双循环版本：先把逻辑跑通

Matlab实现和C++在思路上完全一样，只是索引从1开始，所以坐标公式要做半像素偏移。下面这个版本没有做任何优化，纯粹为了展示逻辑：

function sino = my_radon_loop(img, angles_deg, detector_count) [rows, cols] = size(img); if nargin < 3 detector_count = 2 * ceil(sqrt(rows^2 + cols^2) / 2); end n_angles = length(angles_deg); sino = zeros(n_angles, detector_count); spacing = 1.0; half_det = detector_count / 2; half_rows = rows / 2; half_cols = cols / 2; for a = 1:n_angles theta = angles_deg(a) * pi / 180; ct = cos(theta); st = sin(theta); for i = 1:rows y = half_rows - i + 0.5; for j = 1:cols x = j - 0.5 - half_cols; s = x * ct + y * st; s_idx = s / spacing + half_det + 0.5; idx0 = floor(s_idx); frac = s_idx - idx0; if idx0 >= 1 && idx0 <= detector_count sino(a, idx0) = sino(a, idx0) + img(i, j) * (1 - frac); end if idx0 + 1 >= 1 && idx0 + 1 <= detector_count sino(a, idx0 + 1) = sino(a, idx0 + 1) + img(i, j) * frac; end end end end end

这里matlab的floor(s_idx)和C++里(int)floor(s_idx)行为一致，不信你可以拿一个小矩阵自己试。跑通之后你会看到，128×128图像跑180个角度，双循环可能要几十秒，通常你会在这个过程中深刻理解"为什么MATLAB要避免三重循环"。

3.2 向量化加速：accumarray一行累加

Matlab的性能模型和C++完全不同，向量化是关键。把内层双循环换成meshgrid坐标矩阵，然后用accumarray做累加，可以把三重循环压成单重角度循环。下面是最近邻版本，便于理解向量化原理：

function sino = my_radon_fast(img, angles_deg, detector_count) [rows, cols] = size(img); if nargin < 3 detector_count = 2 * ceil(sqrt(rows^2 + cols^2) / 2); end n_angles = length(angles_deg); sino = zeros(n_angles, detector_count); % 构建像素中心坐标矩阵 [X, Y] = meshgrid((1:cols) - 0.5 - cols/2, rows/2 - (1:rows) + 0.5); for a = 1:n_angles theta = angles_deg(a) * pi / 180; S = X * cos(theta) + Y * sin(theta); idx = round(S + detector_count/2 + 0.5); valid = idx >= 1 & idx <= detector_count; vals = img(valid); idx_valid = idx(valid); sino(a, :) = accumarray(idx_valid, vals, [detector_count, 1])'; end end

accumarray(idx_valid, vals, [detector_count, 1])做的事是把所有落入同一探测器索引的像素值累加，正好对应Radon变换的积分语义。向量化版本同样用128×128图像和180个角度测试，运行时间能从几十秒降到一两秒，差距非常直观。如果你需要精确到半个探测器单元的线性插值版本，可以仿照C++代码写一个稀疏权重矩阵，或者用griddata之类插值函数，不过实测中最近邻和线性插值在重建结果上的差异很小，工程上先用最近邻拿到形状正确的sinogram更重要。

3.3 用内置radon函数做交叉验证

MATLAB自带radon函数是天然的标准答案，用来验证自己的实现非常方便。关键点是显式指定探测器数量，让两个结果的可比性更强。

P = phantom(128); theta = 0:179; N_det = 2 * ceil(size(P, 1) * sqrt(2) / 2); R_mine = my_radon_fast(P, theta, N_det); R_ref = radon(P, theta, N_det); figure; subplot(1,2,1); imagesc(R_ref); colormap gray; axis image; title('内置 radon'); subplot(1,2,2); imagesc(R_mine); colormap gray; axis image; title('手写 radon'); rel_err = norm(R_mine(:) - R_ref(:)) / norm(R_ref(:)); fprintf('相对误差: %.4f\n', rel_err);

如果形状完全一致但左右镜像，说明你的s轴方向和内置radon相反，把结果矩阵R_mine(:, end:-1:1)调过来再看。如果相对误差在5%以内，基本都是插值方式和边界细节的差异，不影响后续做重建实验。需要特别说明的是，内置radon实现的是射线驱动，和我们的像素驱动在线积分近似上存在细小差别，不要追求误差为0，那反而说明你可能过度拟合了某个特定图像。

4. 结果验收与坑位复盘：验证Radon实现的四个"试金石"

4.1 0度和90度投影对照：最廉价的几何自检

代码写完后第一件事不是对着sinogram欣赏，而是用几个思想实验验证几何。取一个大小的方形均匀图像，中心放一个矩形亮块，跑0度、90度、45度三个角度投影，然后检查：

• 0度投影：探测器索引方向对应x轴方向，所以sinogram第1行每一列的值应该恰好等于图像对应列的像素之和。
• 90度投影：s = y，所以每一行的值应该等于图像对应行的像素之和，但因为我们的y轴向上为正，第0行在图像顶部对应y最大值，所以90度投影的列顺序和图像行顺序刚好相反。看到这个"反了"千万不要改，这是坐标约定的自然结果。

这两条检验比任何调试器都管用。如果0度结果对不上列和，你应该先怀疑坐标公式里的半像素偏移，再看角度到弧度的转换是否漏了pi/180。

4.2 单点源实验：sinogram这个名字的由来

几何自检通过后，再来一个更有说服力的实验：构造一张全黑图像，只在某个点比如(20, 10)处放一个像素值为1的点，让后跑0到179度所有投影，把sinogram画出来。你会看到一条标准正弦曲线，这就是"sinogram"这个名字的来源。

原因很简单：该点坐标(x0, y0)在任意角度θ下的投影位置是：

s = x0·cosθ + y0·sinθ = A·cos(θ − φ)

其中A是点到原点的距离，φ是它的极角。所以单点源在投影矩阵里画出的轨迹就是一条正弦曲线。这个实验还能帮你检查角度方向：如果θ增大时曲线是从左往右还是从右往左，就反映了你的角度定义是顺时针还是逆时针。临床CT通常定义探测器从0度开始逆时针旋转，你只要保证和后续重建代码用的坐标系一致即可。

4.3 三个高频坑位的复盘

第一个坑是探测器数量不足。我用N=512的图像跑45度投影时，发现总能量比0度少了将近四分之一，当时第一反应是代码出了问题，排查半天后才发现是探测器只设了512，覆盖范围[-256,256]，45度时图像四个角都飞出了边界。换成ceil(N√2)之后，总能量偏差降到千分之一量级。

第二个坑是半像素偏移。你把坐标写成x = j − cols/2而不加0.5，和内置radon对比时能看到整体投影位置向左偏移约半个像素。低分辨率图像上这个偏移会被放大，重建图像会出现轻微模糊。

第三个坑是Matlab纯循环的性能。头一回我拿512×512图像跑180个角度，三重纯循环跑了几分钟还没结束，还以为电脑坏了。后来换向量化加accumarray，速度提升上百倍。写Matlab一定要有"能不写循环就不写循环"的警觉，写C++反而要担心内存连续访问。两种语言的优化思路完全相反，不要混着用。

5. 常见问题速查与代码修复建议

5.1 高频问题排查表

下面这张表我整理自实践，基本上是同行问过或者我自己犯过的错，建议收藏。

5.2 两个容易忽略的细节

插值方式会影响投影数据的平滑度。最近邻法会把像素当成一个个离散的点，投影曲线带毛刺；线性插值则更接近真实射线积分的效果，但它会在边界产生约半个探测器单元的权重泄露，原因是s落在探测器覆盖范围边缘时，有一部分权重被分给了不存在的虚拟探测器。这不是bug，是离散化固有的边界效应，通常让探测器覆盖范围多留一点余量就能把影响压到可忽略。

另一个细节是角度列表的长度。角度越密，sinogram行数越多，重建效果越细腻，但计算量和存储也成比例上升。做仿真时我一般用每度一个投影的角度列表，做算法验证时用每两度一个就够，等需要出论文配图时再加密到0.5度步长。记住CT重建的本质问题是一个逆问题，投影数据越多，逆问题条件数越好，重建伪影越少。

6. 从正向投影到重建：Radon之后的下一步

6.1 中心切片定理：正向和反向的桥

搞定了正向投影，下一个自然的问题是：拿到sinogram之后怎么把图像重建回来？这就必须提中心切片定理（Central Slice Theorem）。它说的是：p(s,θ)关于s做一维傅里叶变换，结果恰好等于原图像f(x,y)的二维傅里叶变换在过原点、方向为(cosθ, sinθ)这条直线上的一组取值。

换句话说，一次角度扫描的投影只是触碰了图像傅里叶空间里的一根"切片线"。扫180个角度，就是在傅里叶空间里填180条过中心的直径线。把所有这些切片填充完再反傅里叶变换，理论上就能恢复图像。

经典的滤波反投影（FBP）做的就是这件事：先对每条投影做傅里叶域的斜坡滤波，再把滤波后的投影沿对应角度反投影回二维空间累加。FBP实现起来比正向投影还直观，代码量也就几十行。理解了正向投影的坐标变换，反投影只是把累加方向倒过来而已，建议你动手写一版，能极大加深对Radon变换的理解。

6.2 从平行束到扇束与螺旋CT：代码只改一处

医院里的CT扫描仪用的不是平行束，而是扇形束X射线。在扇束几何下，投影线不再相互平行，而是从X射线源点呈扇形散开。理论上所有角度下的扇束数据可以通过重排（rebinning）转换为平行束数据，核心公式是：

s = D · sinγ
θ = β + γ

其中D是射线源到旋转中心的距离，β是源旋转角，γ是扇束内某条射线的扇角。做完这个坐标映射，扇形束问题就变成了平行束问题，可以直接复用前面实现的Radon正投影和FBP重建。

至于螺旋CT，核心变化是扫描床上物体在扫描过程中持续平移，所以投影数据和切片位置不再一一对应，需要在重建前做螺旋插值。但这些都是工程层的扩展，底层依然是Radon变换这个数学模型。低剂量CT、AI超分辨率重建、深度学习方法生成训练数据，本质上也都是在和这套正向投影/重建的框架打交道。

我个人在实际调试Radon代码时最深的体会是：这个算法本身不难，难的地方全在坐标约定和边界处理上。建议你一定要从8×8或者16×16的小图开始调，肉眼能直接看到结果正确后再逐步放大。另外，不管用C++还是Matlab，先写一个手动的双循环版本理清逻辑，再去做向量化或者并行优化，顺序反了你会很痛苦。