锥束CT（CBCT）三维重构算法：FDK算法详解与实现
锥束CT（CBCT）三维重构算法：FDK算法详解与实现
锥束CT（Cone Beam Computed Tomography）是一种广泛应用于医学影像和工业检测的三维成像技术。FDK（Feldkamp-Davis-Kress）算法是锥束CT重建中最经典的滤波反投影算法，由Feldkamp、Davis和Kress于1984年提出。
FDK算法原理
FDK算法是对传统二维CT中滤波反投影（FBP）算法的扩展，用于处理锥形束投影几何。其核心思想是：
- 加权投影数据：根据几何关系对投影数据进行加权
- 滤波处理：对加权后的投影数据进行斜坡滤波
- 三维反投影：将滤波后的数据反投影到三维空间
数学表达
重建图像\(f(x,y,z)\)可表示为：
\[f(x,y,z) = \frac{1}{2} \int_0^{2\pi} \frac{D^2}{(D + s)^2} \int_{-\infty}^{\infty} \frac{D}{\sqrt{D^2 + u'^2 + v'^2}} P_\beta(u',v') h(u - u') du' d\beta
\]
其中：
- \(D\)：源到旋转中心的距离
- \(s\)：体素点到源的距离在探测器平面上的投影
- \(P_\beta(u',v')\)：角度\(\beta\)处的投影数据
- \(h(u)\)：斜坡滤波器
MATLAB代码
function volume = fdk_reconstruction(projections, angles, DSD, DDR, vol_size, voxel_size)
% FDK锥束CT三维重建算法
% 输入：
% projections: 投影数据 [探测器行数, 探测器列数, 投影角度数]
% angles: 投影角度数组 (弧度)
% DSD: 源到探测器的距离 (mm)
% DDR: 源到旋转中心的距离 (mm)
% vol_size: 重建体积尺寸 [x, y, z]
% voxel_size: 体素尺寸 [dx, dy, dz] (mm)
% 输出：
% volume: 重建后的三维体积数据
% 参数设置
[num_rows, num_cols, num_angles] = size(projections);
vol_x = vol_size(1); vol_y = vol_size(2); vol_z = vol_size(3);
dx = voxel_size(1); dy = voxel_size(2); dz = voxel_size(3);
% 创建滤波器（斜坡滤波器加Hamming窗）
filter = design_filter(num_cols);
% 初始化重建体积
volume = zeros(vol_x, vol_y, vol_z, 'single');
% 计算探测器参数
du = (2 * DDR * tan(0.5 * atan(num_cols * 0.5 / DSD))) / num_cols;
dv = du; % 假设正方形像素
u_center = (num_cols + 1) / 2;
v_center = (num_rows + 1) / 2;
% 创建网格坐标
[x, y, z] = meshgrid(...
((0:vol_x-1) - (vol_x-1)/2) * dx, ...
((0:vol_y-1) - (vol_y-1)/2) * dy, ...
((0:vol_z-1) - (vol_z-1)/2) * dz);
% 主重建循环
for angle_idx = 1:num_angles
beta = angles(angle_idx);
fprintf('处理角度 %d/%d (%.1f°)\n', angle_idx, num_angles, rad2deg(beta));
% 获取当前角度投影
proj = projections(:, :, angle_idx);
% 1. 加权投影
[U, V] = meshgrid(((1:num_cols) - u_center) * du, ((1:num_rows) - v_center) * dv);
weight_factor = DSD ./ sqrt(DSD^2 + U.^2 + V.^2);
weighted_proj = proj .* weight_factor;
% 2. 滤波处理
filtered_proj = filter_projections(weighted_proj, filter);
% 3. 三维反投影
% 计算旋转矩阵
R = [cos(beta) sin(beta) 0;
-sin(beta) cos(beta) 0;
0 0 1];
% 坐标旋转
coords = [x(:)'; y(:)'; z(:)'];
rotated_coords = R * coords;
xr = reshape(rotated_coords(1,:), size(x));
yr = reshape(rotated_coords(2,:), size(y));
zr = reshape(rotated_coords(3,:), size(z));
% 计算投影坐标
u_proj = (DSD * yr) ./ (DDR - xr) / du + u_center;
v_proj = (DSD * zr) ./ (DDR - xr) / dv + v_center;
% 反投影权重
weight = DSD^2 ./ (DSD - xr).^2;
% 插值获取投影值
valid_idx = (u_proj >= 1) & (u_proj <= num_cols) & ...
(v_proj >= 1) & (v_proj <= num_rows);
proj_vals = zeros(size(x), 'single');
proj_vals(valid_idx) = interp2(filtered_proj, u_proj(valid_idx), v_proj(valid_idx), 'linear', 0);
% 累加到体积
volume = volume + proj_vals .* weight;
end
% 角度积分归一化
volume = volume * (2 * pi / num_angles);
end
%% 设计滤波器
function filter = design_filter(length)
% 设计斜坡滤波器（加Hamming窗）
n = 0:(length-1);
filter = zeros(1, length);
% 斜坡滤波器
filter(1) = 1 / (4); % 中心点
even = mod(length, 2) == 0;
if even
filter(2:end) = -1./(pi^2 * (n(2:end) - 0.5).^2);
else
filter(2:end) = -1./(pi^2 * n(2:end).^2);
end
% 加窗（Hamming窗）
window = hamming(length)';
filter = filter .* window;
% 归一化
filter = filter / sum(abs(filter));
end
%% 投影滤波
function filtered_proj = filter_projections(proj, filter)
% 对投影数据进行滤波
[num_rows, num_cols] = size(proj);
filtered_proj = zeros(num_rows, num_cols, 'single');
% 扩展投影数据（避免边界效应）
padded_proj = zeros(num_rows, 2 * num_cols);
padded_proj(:, 1:num_cols) = proj;
% 频域滤波
for i = 1:num_rows
row = padded_proj(i, :);
freq_row = fft(row);
filtered_freq = freq_row .* fftshift(filter);
filtered_row = real(ifft(filtered_freq));
filtered_proj(i, :) = filtered_row(1:num_cols);
end
end
%% 示例使用
function run_fdk_example()
% FDK重建示例
% 生成合成投影数据（Shepp-Logan 3D体模）
vol_size = [128, 128, 128]; % 体积尺寸
voxel_size = [0.5, 0.5, 0.5]; % 体素大小 (mm)
num_angles = 360; % 投影角度数
angles = linspace(0, 2*pi, num_angles); % 投影角度
DSD = 1000; % 源到探测器距离 (mm)
DDR = 500; % 源到旋转中心距离 (mm)
% 创建Shepp-Logan 3D体模
phantom_vol = phantom3d(vol_size(1));
% 生成投影数据
projections = simulate_projections(phantom_vol, angles, DSD, DDR, vol_size, voxel_size);
% 执行FDK重建
recon_vol = fdk_reconstruction(projections, angles, DSD, DDR, vol_size, voxel_size);
% 显示结果
figure;
subplot(1,3,1);
imshow(squeeze(phantom_vol(:,:,round(vol_size(3)/2))), []);
title('原始体模（中心切片）');
subplot(1,3,2);
imshow(squeeze(recon_vol(:,:,round(vol_size(3)/2))), []);
title('重建结果（中心切片）');
subplot(1,3,3);
imshow(projections(:,:,1), []);
title('投影数据示例');
% 显示三维重建结果
figure;
show_volume(recon_vol);
title('三维重建结果');
end
%% 辅助函数：三维体模生成
function vol = phantom3d(n)
% 生成3D Shepp-Logan体模
vol = zeros(n, n, n);
[x, y, z] = meshgrid(linspace(-1,1,n), linspace(-1,1,n), linspace(-1,1,n));
% 定义椭球参数 [A, a, b, c, x0, y0, z0, phi, theta, psi]
ellipsoids = [
2.0 0.69 0.92 0.90 0.0 0.0 0.0 0 0 0
-0.8 0.6624 0.874 0.88 0.0 0.0 0.0 0 0 0
-0.2 0.1100 0.310 0.22 0.22 0.0 0.0 -18 0 10
-0.2 0.1600 0.410 0.28 -0.22 0.0 0.0 18 0 10
0.1 0.2100 0.250 0.41 0.0 0.35 0.0 0 0 0
0.1 0.0460 0.046 0.05 0.0 0.1 0.0 0 0 0
0.1 0.0460 0.046 0.05 0.0 -0.1 0.0 0 0 0
0.1 0.0460 0.023 0.05 -0.08 -0.65 0.0 0 0 0
0.1 0.0230 0.023 0.02 0.0 -0.65 0.0 0 0 0
0.1 0.0230 0.046 0.02 0.06 -0.65 0.0 0 0 0
];
for i = 1:size(ellipsoids, 1)
A = ellipsoids(i,1);
a = ellipsoids(i,2);
b = ellipsoids(i,3);
c = ellipsoids(i,4);
x0 = ellipsoids(i,5);
y0 = ellipsoids(i,6);
z0 = ellipsoids(i,7);
phi = deg2rad(ellipsoids(i,8));
theta = deg2rad(ellipsoids(i,9));
psi = deg2rad(ellipsoids(i,10));
% 旋转矩阵
Rz = [cos(psi) -sin(psi) 0; sin(psi) cos(psi) 0; 0 0 1];
Ry = [cos(theta) 0 sin(theta); 0 1 0; -sin(theta) 0 cos(theta)];
Rx = [1 0 0; 0 cos(phi) -sin(phi); 0 sin(phi) cos(phi)];
R = Rz * Ry * Rx;
% 坐标变换
coords = [x(:)-x0, y(:)-y0, z(:)-z0] * R';
X = reshape(coords(:,1), size(x));
Y = reshape(coords(:,2), size(y));
Z = reshape(coords(:,3), size(z));
% 椭球方程
idx = (X/a).^2 + (Y/b).^2 + (Z/c).^2 <= 1;
vol(idx) = vol(idx) + A;
end
end
%% 辅助函数：投影模拟
function projections = simulate_projections(vol, angles, DSD, DDR, vol_size, voxel_size)
% 模拟锥束CT投影
[n, m, p] = size(vol);
num_angles = length(angles);
projections = zeros(m, p, num_angles, 'single');
% 计算探测器尺寸
max_dim = max(vol_size .* voxel_size);
detector_size = ceil(2 * max_dim * DSD / DDR);
% 创建探测器网格
[u, v] = meshgrid(linspace(-1,1,p), linspace(-1,1,m));
u = u * detector_size/2;
v = v * detector_size/2;
% 生成投影
for i = 1:num_angles
beta = angles(i);
R = [cos(beta) sin(beta) 0;
-sin(beta) cos(beta) 0;
0 0 1];
proj = zeros(m, p);
% 计算投影矩阵
for j = 1:m
for k = 1:p
% 计算射线方向
source = [0; -DDR; 0]; % 源位置
detector_point = [0; DDR; 0] + [u(j,k); 0; v(j,k)];
% 旋转坐标系
ray_dir = R * (detector_point - source);
ray_dir = ray_dir / norm(ray_dir);
% 计算射线与体模的交点
% (这里简化处理，实际应用中应使用Siddon算法等)
t_min = -max_dim;
t_max = max_dim;
% 步进累加
num_steps = 200;
step_size = (t_max - t_min) / num_steps;
value = 0;
for t = t_min:step_size:t_max
point = source + t * ray_dir;
% 转换为体素坐标
idx_x = round((point(1) + max_dim/2) / voxel_size(1));
idx_y = round((point(2) + max_dim/2) / voxel_size(2));
idx_z = round((point(3) + max_dim/2) / voxel_size(3));
if idx_x >= 1 && idx_x <= n && ...
idx_y >= 1 && idx_y <= m && ...
idx_z >= 1 && idx_z <= p
value = value + vol(idx_x, idx_y, idx_z) * step_size;
end
end
proj(j,k) = value;
end
end
projections(:,:,i) = proj;
end
end
%% 辅助函数：三维数据显示
function show_volume(vol)
% 显示三维体积
% 创建等值面
p = patch(isosurface(vol, 0.5));
isonormals(vol, p);
p.FaceColor = 'red';
p.EdgeColor = 'none';
% 设置视图
daspect([1 1 1]);
view(3);
axis tight;
camlight;
lighting gouraud;
end
FDK算法关键步骤详解
1. 投影数据加权
在锥束几何中，由于X射线路径长度不同，需要根据几何关系对投影数据进行加权：
weight_factor = DSD ./ sqrt(DSD^2 + U.^2 + V.^2);
weighted_proj = proj .* weight_factor;
其中：
DSD：源到探测器的距离U,V：探测器上的坐标
2. 滤波处理
FDK算法使用斜坡滤波器（Ram-Lak滤波器）进行滤波处理，通常加窗（如Hamming窗）以减少噪声：
function filter = design_filter(length)
% 斜坡滤波器核心
filter(1) = 1 / (4);
filter(2:end) = -1./(pi^2 * n(2:end).^2);
% 加Hamming窗
window = hamming(length)';
filter = filter .* window;
end
3. 三维反投影
反投影过程将滤波后的投影数据映射回三维空间：
% 计算投影坐标
u_proj = (DSD * yr) ./ (DDR - xr) / du + u_center;
v_proj = (DSD * zr) ./ (DDR - xr) / dv + v_center;
% 反投影权重
weight = DSD^2 ./ (DSD - xr).^2;
% 插值获取投影值
proj_vals(valid_idx) = interp2(filtered_proj, u_proj(valid_idx), v_proj(valid_idx), 'linear', 0);
% 累加到体积
volume = volume + proj_vals .* weight;
参考 锥束CT（CBCT）的三维重构算法（FDK）youwenfan.com/contentcna/65715.html
算法性能优化
-
并行计算：
parfor angle_idx = 1:num_angles % 处理每个角度的投影 end -
GPU加速：
% 将数据转移到GPU projections_gpu = gpuArray(projections); volume_gpu = gpuArray(zeros(vol_size, 'single')); % 在GPU上执行计算 % ... % 将结果传回CPU volume = gather(volume_gpu); -
插值优化：
- 使用最近邻插值加速（牺牲精度）
- 使用三次样条插值提高精度（增加计算量）
总结
FDK算法作为锥束CT三维重建的基石算法，具有实现简单、计算效率高的优点。尽管存在锥束伪影等局限性，通过合理的优化和改进，仍然是实际应用中广泛使用的重建方法。随着计算能力的提升和新型算法的出现，锥束CT重建技术将继续向更高精度、更低剂量、更快速度的方向发展。
本实现提供了FDK算法的完整MATLAB实现，包括投影模拟、滤波处理、三维反投影等关键步骤，以及结果可视化功能，可作为锥束CT研究和开发的基础框架。
浙公网安备 33010602011771号