★ FEATURED ARTICLE 简介这份资源是面向计算机、电子信息工程、数学等专业学生的锥形束CTCBCT图像重建Matlab示例包聚焦投影背投、FDK解析重建与迭代重建两类核心算法可用于课程设计、期末大作业与毕业设计也适合作为医学影像算法入门与原型验证的实践素材。压缩包共15个文件以12个m脚本为主体另含2个txt说明与1个mat案例数据整体约173KB脚本按参数设置、投影生成、滤波、加权、FDK重建及MLEM、SART、SQS等迭代重建分模块组织便于对照理解算法流程。资源采用参数化编程注释清晰配套案例数据可直接运行读者能借此掌握从投影数据到三维体积数据的完整重建链路并对比解析法与迭代法在噪声和数据不完整情形下的表现差异。目前已有148人学习下载适合需要快速上手CBCT重建实验的学习者。1. 从一张投影图到三维体数据这套 FDK 与迭代重建 Matlab 示例到底能跑出什么手里只有一圈二维投影图却要还原出物体内部的三维密度分布这是锥形束 CT 重建最核心的诉求。这套 Matlab 示例把两条主流路线都摆了出来一条是解析法里的 FDK 投影背投另一条是迭代重建。前者快、结构清晰适合先把几何关系和滤波背投链路跑通后者慢但在投影数少、噪声大、金属伪影明显的场景里往往能救回图像质量。它适合正在做 CBCT 算法验证、想拿现成脚本改参数看效果的人也适合刚接触三维重建、需要一份能直接运行、能逐行打断点看中间变量的工程参考。整套代码围绕锥束几何、滤波函数、体素遍历和迭代更新展开不是纯理论推导而是能落到矩阵运算和循环里的实现。2. FDK 解析重建滤波、锥束加权与背投的完整链路2.1 为什么先跑 FDK 而不是直接上迭代FDK 是锥形束圆轨迹重建里最经典的解析方法本质是把三维锥束问题拆成一系列二维扇束问题再沿探测器行方向做加权和滤波最后把滤波后的投影按锥束几何背投回三维体素。它的优势是计算量可控一次重建基本就是“加权—滤波—背投”三步不需要反复迭代适合先验证投影数据、几何参数和体素网格是否对齐。如果 FDK 都跑不出结构迭代重建大概率也救不回来因为几何或数据本身就有问题。常见做法是先用小尺寸体数据跑通 FDK确认重建切片里能看到物体轮廓再换迭代。示例里一般会提供投影生成、滤波和背投三个独立函数方便你单独替换某一环做对比。2.2 投影数据的组织与几何参数CBCT 投影数据通常组织成三维矩阵探测器行、探测器列、投影角度。假设投影数为N_proj探测器尺寸为N_u × N_v数据就是N_v × N_u × N_proj。几何参数包括源到旋转中心距离、源到探测器距离、探测器像素尺寸、旋转角度步长。这些参数一旦有一个对不上重建结果就会模糊、重影甚至完全错位。下面是一段常见的投影数据读取与几何参数初始化代码用来把数据整理成后续滤波和背投能直接用的形式% 读取投影数据假设为单精度浮点尺寸 [N_v, N_u, N_proj] proj load(projection_data.mat); proj proj.data; % 几何参数按实际扫描设置修改 geo.SDD 1000; % 源到探测器距离单位 mm geo.SOD 600; % 源到旋转中心距离单位 mm geo.du 0.8; % 探测器像素尺寸 u 方向mm geo.dv 0.8; % 探测器像素尺寸 v 方向mm geo.angles linspace(0, 2*pi, size(proj,3)1); geo.angles(end) []; % 去掉重复的 2*pi % 体素网格 vol.Nx 256; vol.Ny 256; vol.Nz 256; vol.dx 1.0; vol.dy 1.0; vol.dz 1.0; % 体素尺寸 mm逻辑说明SDD和SOD决定锥束张角du、dv决定探测器采样密度angles决定背投时每个投影对应的旋转位置。体素网格的Nx/Ny/Nz和dx/dy/dz要和重建目标匹配太大跑得慢太小看不清细节。参数改错最典型的表现是重建图像整体放大或缩小或者边缘出现明显错层。2.3 滤波函数的选择与实现FDK 里的滤波通常是斜坡滤波加窗常见的有 Ram-Lak、Shepp-Logan、Cosine、Hamming 等。Ram-Lak 分辨率最高但噪声也最明显Shepp-Logan 和 Cosine 会压高频、降噪但牺牲一点细节。示例里一般会提供一个filter_projection函数对每个角度的投影沿 u 方向做一维滤波。function proj_f filter_projection(proj, du, filter_type) % 对每个角度的投影沿探测器 u 方向做斜坡滤波 [N_v, N_u, N_proj] size(proj); proj_f zeros(size(proj), like, proj); % 构造频域滤波窗 n 2^nextpow2(2*N_u); f linspace(-1/(2*du), 1/(2*du), n); ramp abs(f); switch lower(filter_type) case ram-lak window ones(size(ramp)); case shepp-logan window sinc(f / (2*max(f))); case cosine window cos(pi * f / (2*max(f))); otherwise error(未知滤波类型); end H ramp .* window; H ifftshift(H); for k 1:N_proj for j 1:N_v row squeeze(proj(j,:,k)); row_f fft(row, n); row_f row_f .* H(:); row_b real(ifft(row_f, n)); proj_f(j,:,k) row_b(1:N_u); end end end逻辑说明先把投影沿 u 方向补零到n点转到频域乘上斜坡窗再反变换回空域。filter_type控制窗函数du影响频率轴刻度。参数改错时滤波太强会让图像发糊滤波太弱会留下明显条纹。注意补零长度要足够否则会出现环绕伪影。2.4 锥束加权与背投滤波之后要做锥束加权通常按SOD / sqrt(SDD^2 u^2 v^2)之类的因子对投影做逐像素修正然后再背投。背投就是遍历每个体素计算它在每个角度下投影到探测器的位置做双线性插值累加。% 锥束加权 [u_grid, v_grid] meshgrid( ... ((1:N_u) - (N_u1)/2) * du, ... ((1:N_v) - (N_v1)/2) * dv); weight geo.SOD ./ sqrt(geo.SDD^2 u_grid.^2 v_grid.^2); weight reshape(weight, [N_v, N_u, 1]); proj_w proj_f .* weight; % 背投简化版逐体素遍历 vol_rec zeros(vol.Nx, vol.Ny, vol.Nz, single); for iz 1:vol.Nz z (iz - (vol.Nz1)/2) * vol.dz; for iy 1:vol.Ny y (iy - (vol.Ny1)/2) * vol.dy; for ix 1:vol.Nx x (ix - (vol.Nx1)/2) * vol.dx; acc 0; for k 1:N_proj theta geo.angles(k); % 旋转到探测器坐标系 u x*cos(theta) y*sin(theta); v z; % 投影到探测器像素 u_idx u * geo.SDD / (geo.SOD x*sin(theta) - y*cos(theta)) / geo.du (N_u1)/2; v_idx v * geo.SDD / (geo.SOD x*sin(theta) - y*cos(theta)) / geo.dv (N_v1)/2; if u_idx 1 u_idx N_u v_idx 1 v_idx N_v acc acc interp2(proj_w(:,:,k), u_idx, v_idx, linear, 0); end end vol_rec(ix,iy,iz) acc; end end end逻辑说明加权是为了补偿锥束路径长度差异背投是把每个体素在所有角度下的贡献累加。u_idx、v_idx是体素投影到探测器的连续坐标interp2做双线性插值。循环顺序可以调整但要注意内存访问模式Matlab 里按列优先遍历会快一些。这段代码跑小体数据可以大体数据建议向量化或分块。3. 迭代重建从 ART 到 SIRT 的落地写法3.1 迭代重建适合什么场景当投影角度少、噪声大、或者有金属植入物时FDK 往往会出现明显条纹和伪影。迭代重建把重建问题当成一个线性方程组Ax b来解A是系统矩阵x是待重建体素b是投影数据。常见的有 ART、SIRT、SART 等。ART 每次只用一个投影更新收敛快但不稳定SIRT 用所有投影的残差做平均更新稳定但慢SART 折中按角度分块更新。示例里一般会提供至少一种迭代方法的实现方便和 FDK 做对比。迭代重建的代价是计算量大一次重建可能要几十到几百次迭代每次迭代都要做一次前投影和一次背投影。3.2 SIRT 的更新公式与实现SIRT 的更新可以写成x^(n1) x^n C A^T R (b - A x^n)其中C和R是归一化矩阵通常取每行或每列和倒数。实现时前投影和背投影可以用和 FDK 类似的射线驱动或体素驱动方式。% SIRT 迭代重建简化版 x zeros(vol.Nx, vol.Ny, vol.Nz, single); num_iter 50; relax 0.5; % 松弛因子 % 预计算归一化这里用简化常数实际应按射线长度累加 C 1 / N_proj; R 1 / (vol.Nx * vol.Ny * vol.Nz); for iter 1:num_iter % 前投影Ax Ax forward_project(x, geo, vol); % 残差 residual proj - Ax; % 背投影A^T residual back back_project(residual, geo, vol); % 更新 x x relax * C * R * back; % 可选非负约束 x(x 0) 0; fprintf(迭代 %d残差范数 %.4f\n, iter, norm(residual(:))); end逻辑说明forward_project和back_project需要和 FDK 里的几何一致否则迭代不收敛。relax控制步长太大震荡太小收敛慢。非负约束在 CT 里很常用因为密度不能为负。残差范数可以用来判断是否收敛但不要只看它还要看图像视觉质量。3.3 迭代次数与正则化的取舍迭代次数不是越多越好。次数太少图像欠重建次数太多噪声和伪影会被放大。常见做法是跑 20 到 100 次中途保存几个中间结果对比。如果噪声明显可以加 TV 正则化或早停。示例里如果没带正则化可以自己加一个简单的平滑项但要注意别把细节也抹掉。提示迭代重建跑之前先用 FDK 确认几何参数正确否则迭代只会把错误几何放大成更奇怪的伪影。4. 避坑与排查几何、滤波和内存这三处最容易翻车4.1 重建图像整体错位或重影现象重建切片里物体位置偏移或者出现双影。原因旋转中心、源到探测器距离、角度步长有一个对不上。解决先用一个已知位置的球体或点源做仿真检查投影上的位置和重建位置是否一致。角度步长如果是2*pi/N_proj注意不要多算一个重复角度。4.2 滤波后图像发糊或条纹明显现象FDK 重建结果模糊或者有明显条纹。原因滤波窗选择不当或者补零长度不够导致环绕。解决换 Ram-Lak 对比确认是不是窗函数问题把补零长度加到2*N_u以上检查du是否和探测器实际像素尺寸一致。4.3 迭代重建不收敛或震荡现象残差范数不降反升或者图像出现棋盘格。原因前投影和背投影的几何不一致或者松弛因子太大。解决先用一个简单仿真数据验证forward_project和back_project是否互为转置把relax降到 0.1 到 0.3 再试检查归一化矩阵是否合理。4.4 内存爆掉或跑得太慢现象Matlab 报内存不足或者循环跑几个小时。原因体素网格太大或者背投用了逐体素四重循环。解决先把Nx/Ny/Nz降到 128 或 64 跑通把背投向量化或者用parfor按切片并行投影数据用single而不是double。4.5 投影数据和体素网格单位不统一现象重建结果尺寸明显不对。原因du、dv、dx单位混用比如一个用 mm 一个用 cm。解决把所有几何参数统一到同一单位通常用 mm检查SOD和SDD的比例是否合理。5. 进阶技巧用仿真数据验证重建链路再换真实投影5.1 先造一个已知体模跑通全链路拿真实投影直接调重建出了问题很难判断是数据问题还是代码问题。我一般会先用 Shepp-Logan 体模或几个球体造仿真投影跑一遍 FDK 和 SIRT看重建结果和原始体模的差异。这样几何、滤波、背投、迭代每一环都能单独验证。% 造一个简单体模中心球 两个偏心球 vol_true zeros(vol.Nx, vol.Ny, vol.Nz, single); [xx, yy, zz] ndgrid( ... (1:vol.Nx) - (vol.Nx1)/2, ... (1:vol.Ny) - (vol.Ny1)/2, ... (1:vol.Nz) - (vol.Nz1)/2); vol_true(sqrt(xx.^2 yy.^2 zz.^2) 20) 1.0; vol_true(sqrt((xx-30).^2 yy.^2 zz.^2) 10) 0.8; vol_true(sqrt((xx30).^2 (yy-20).^2 zz.^2) 8) 0.6; % 前投影生成仿真数据 proj_sim forward_project(vol_true, geo, vol); % 用仿真数据跑 FDK 和 SIRT对比 vol_fdk fdk_reconstruct(proj_sim, geo, vol); vol_sirt sirt_reconstruct(proj_sim, geo, vol, 50, 0.3); % 计算均方误差 mse_fdk mean((vol_fdk(:) - vol_true(:)).^2); mse_sirt mean((vol_sirt(:) - vol_true(:)).^2); fprintf(FDK MSE: %.6f, SIRT MSE: %.6f\n, mse_fdk, mse_sirt);逻辑说明vol_true是已知真值forward_project生成仿真投影再分别用 FDK 和 SIRT 重建最后算 MSE。MSE 不是唯一指标但能快速判断链路是否跑通。如果 FDK 的 MSE 比 SIRT 还大很多说明几何或滤波有问题。5.2 换真实投影时的检查清单仿真跑通后换真实投影按这个顺序检查投影数据是否做了暗场和亮场校正角度范围是否覆盖 360 度探测器像素尺寸和几何参数是否来自同一份标定文件投影数据是否需要取负对数。常见做法是先重建一个中心切片看有没有明显伪影再决定要不要调滤波或加迭代。5.3 参数对比表参数FDK 常用值SIRT 常用值影响滤波类型Ram-Lak / Shepp-Logan不适用影响分辨率和噪声迭代次数120~100次数越多越慢噪声可能放大松弛因子不适用0.1~0.5太大震荡太小收敛慢体素尺寸0.5~2 mm0.5~2 mm越小越慢细节越多投影数360~72060~360越少迭代优势越明显5.4 我踩过的一个坑有一次换了一套投影数据FDK 重建出来全是条纹查了半天以为是滤波问题最后发现是角度步长用了2*pi/(N_proj-1)而数据实际是2*pi/N_proj。从那以后我每次换数据都强制先跑一遍仿真体模确认几何参数和角度定义一致再动真实投影。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
