当前位置：首页 > news >正文

三维CT图像重建算法

news 2025/9/23 14:55:11

三维CT图像重建是一个复杂的过程，涉及从多个二维投影重建三维体积。下面我将提供一个完整的MATLAB实现，包括滤波反投影（FBP）算法和迭代重建算法（如SART），以及相关的预处理和后处理步骤。

三维CT重建算法实现

% 三维CT图像重建算法
clear; clc; close all;%% 1. 参数设置
fprintf('设置重建参数...\n');% 扫描几何参数
num_angles = 360;              % 投影角度数
angles = linspace(0, 360-360/num_angles, num_angles); % 角度范围 (度)
source_to_detector = 1000;     % 源到探测器的距离 (mm)
source_to_object = 500;        % 源到物体的距离 (mm)
detector_size = 256;           % 探测器尺寸 (像素)
detector_pixel_size = 0.5;     % 探测器像素尺寸 (mm)% 重建参数
volume_size = 256;             % 重建体积尺寸 (体素)
voxel_size = 0.25;             % 体素尺寸 (mm)
num_iterations = 10;           % 迭代算法迭代次数% 滤波器参数
filter_type = 'ram-lak';       % 滤波器类型: 'ram-lak', 'shepp-logan', 'cosine', 'hamming'%% 2. 生成模拟三维Shepp-Logan模体
fprintf('生成模拟三维Shepp-Logan模体...\n');
phantom_3d = generate_3d_shepp_logan(volume_size);% 显示模体切片
figure;
subplot(2,2,1);
imagesc(squeeze(phantom_3d(volume_size/2, :, :)));
title('XY平面中心切片');
axis image; colormap gray; colorbar;subplot(2,2,2);
imagesc(squeeze(phantom_3d(:, volume_size/2, :)));
title('XZ平面中心切片');
axis image; colormap gray; colorbar;subplot(2,2,3);
imagesc(squeeze(phantom_3d(:, :, volume_size/2)));
title('YZ平面中心切片');
axis image; colormap gray; colorbar;subplot(2,2,4);
isosurface(phantom_3d, 0.5);
title('三维等值面');
axis equal; axis vis3d; camlight; lighting gouraud;%% 3. 模拟投影数据采集 (正向投影)
fprintf('模拟投影数据采集...\n');
sinogram = forward_project_3d(phantom_3d, angles, detector_size, ...source_to_detector, source_to_object, ...detector_pixel_size, voxel_size);% 显示正弦图
figure;
imagesc(sinogram);
title('正弦图 (所有角度)');
xlabel('探测器位置');
ylabel('投影角度');
colormap gray; colorbar;%% 4. 添加噪声 (模拟真实CT数据)
fprintf('添加噪声模拟真实数据...\n');
noise_level = 0.05; % 噪声水平
sinogram_noisy = add_ct_noise(sinogram, noise_level);% 显示带噪声的正弦图
figure;
imagesc(sinogram_noisy);
title('带噪声的正弦图');
xlabel('探测器位置');
ylabel('投影角度');
colormap gray; colorbar;%% 5. 滤波反投影 (FBP) 重建
fprintf('开始滤波反投影重建...\n');
recon_fbp = fbp_reconstruct_3d(sinogram_noisy, angles, volume_size, ...source_to_detector, source_to_object, ...detector_pixel_size, voxel_size, filter_type);% 显示FBP重建结果
figure;
subplot(2,2,1);
imagesc(squeeze(recon_fbp(volume_size/2, :, :)));
title('FBP: XY平面中心切片');
axis image; colormap gray; colorbar;subplot(2,2,2);
imagesc(squeeze(recon_fbp(:, volume_size/2, :)));
title('FBP: XZ平面中心切片');
axis image; colormap gray; colorbar;subplot(2,2,3);
imagesc(squeeze(recon_fbp(:, :, volume_size/2)));
title('FBP: YZ平面中心切片');
axis image; colormap gray; colorbar;subplot(2,2,4);
isosurface(recon_fbp, 0.3);
title('FBP: 三维等值面');
axis equal; axis vis3d; camlight; lighting gouraud;%% 6. 迭代重建算法 (SART)
fprintf('开始SART迭代重建...\n');
recon_sart = sart_reconstruct_3d(sinogram_noisy, angles, volume_size, ...source_to_detector, source_to_object, ...detector_pixel_size, voxel_size, num_iterations);% 显示SART重建结果
figure;
subplot(2,2,1);
imagesc(squeeze(recon_sart(volume_size/2, :, :)));
title('SART: XY平面中心切片');
axis image; colormap gray; colorbar;subplot(2,2,2);
imagesc(squeeze(recon_sart(:, volume_size/2, :)));
title('SART: XZ平面中心切片');
axis image; colormap gray; colorbar;subplot(2,2,3);
imagesc(squeeze(recon_sart(:, :, volume_size/2)));
title('SART: YZ平面中心切片');
axis image; colormap gray; colorbar;subplot(2,2,4);
isosurface(recon_sart, 0.3);
title('SART: 三维等值面');
axis equal; axis vis3d; camlight; lighting gouraud;%% 7. 重建质量评估
fprintf('评估重建质量...\n');% 计算均方根误差 (RMSE)
rmse_fbp = compute_rmse(phantom_3d, recon_fbp);
rmse_sart = compute_rmse(phantom_3d, recon_sart);% 计算结构相似性指数 (SSIM)
ssim_fbp = compute_ssim_3d(phantom_3d, recon_fbp);
ssim_sart = compute_ssim_3d(phantom_3d, recon_sart);fprintf('重建质量指标:\n');
fprintf('FBP - RMSE: %.4f, SSIM: %.4f\n', rmse_fbp, ssim_fbp);
fprintf('SART - RMSE: %.4f, SSIM: %.4f\n', rmse_sart, ssim_sart);% 绘制中心行剖面比较
center_slice = volume_size/2;
phantom_profile = squeeze(phantom_3d(center_slice, center_slice, :));
fbp_profile = squeeze(recon_fbp(center_slice, center_slice, :));
sart_profile = squeeze(recon_sart(center_slice, center_slice, :));figure;
plot(phantom_profile, 'k-', 'LineWidth', 2);
hold on;
plot(fbp_profile, 'r--', 'LineWidth', 1.5);
plot(sart_profile, 'b-.', 'LineWidth', 1.5);
legend('原始模体', 'FBP重建', 'SART重建');
title('中心行剖面比较');
xlabel('体素位置');
ylabel('强度值');
grid on;%% 8. 辅助函数定义function phantom_3d = generate_3d_shepp_logan(size)% 生成三维Shepp-Logan模体% 基于2D Shepp-Logan模体扩展为3D% 生成2D Shepp-Logan模体phantom_2d = phantom('Modified Shepp-Logan', size);% 扩展为3Dphantom_3d = zeros(size, size, size);for z = 1:size% 在Z方向上添加一些变化以创建3D结构scale = 1 - 0.5 * abs(z - size/2) / (size/2);phantom_3d(:, :, z) = phantom_2d * scale;% 添加一些球体以创建真正的3D结构if z > size/4 && z < 3*size/4% 在中心区域添加球体[x, y] = meshgrid(1:size, 1:size);center = size/2;radius = size/8;% 创建球体sphere = ((x - center).^2 + (y - center).^2) <= radius^2;phantom_3d(:, :, z) = phantom_3d(:, :, z) + 0.3 * sphere;endend% 归一化到[0, 1]范围phantom_3d = (phantom_3d - min(phantom_3d(:))) / (max(phantom_3d(:)) - min(phantom_3d(:)));
endfunction sinogram = forward_project_3d(volume, angles, detector_size, ...source_to_detector, source_to_object, ...detector_pixel_size, voxel_size)% 三维正向投影 - 模拟CT扫描过程num_angles = length(angles);sinogram = zeros(num_angles, detector_size);% 计算几何参数magnification = source_to_detector / source_to_object;detector_length = detector_size * detector_pixel_size;% 创建探测器坐标detector_coords = linspace(-detector_length/2, detector_length/2, detector_size);% 对每个角度进行投影for i = 1:num_anglesangle = angles(i);% 旋转体积rotated_volume = imrotate3(volume, angle, [0 0 1], 'crop', 'linear');% 计算投影 (沿X轴积分)projection = squeeze(sum(rotated_volume, 1));% 考虑几何放大projection = imresize(projection, [detector_size, detector_size] * magnification);% 提取中心行作为投影center_row = round(size(projection, 1)/2);sinogram(i, :) = projection(center_row, :);end
endfunction sinogram_noisy = add_ct_noise(sinogram, noise_level)% 添加CT噪声模拟真实数据% 确保正弦图为正值sinogram = max(sinogram, 0);% 模拟光子计数噪声 (泊松噪声)max_val = max(sinogram(:));sinogram_noisy = imnoise(sinogram / max_val, 'poisson') * max_val;% 添加高斯噪声sinogram_noisy = imnoise(sinogram_noisy, 'gaussian', 0, noise_level^2);
endfunction recon_volume = fbp_reconstruct_3d(sinogram, angles, volume_size, ...source_to_detector, source_to_object, ...detector_pixel_size, voxel_size, filter_type)% 三维滤波反投影重建num_angles = length(angles);recon_volume = zeros(volume_size, volume_size, volume_size);% 计算几何参数magnification = source_to_detector / source_to_object;detector_length = detector_size * detector_pixel_size;% 创建滤波器filter = create_filter(filter_type, detector_size, detector_pixel_size);% 对每个角度进行反投影for i = 1:num_anglesangle = angles(i);% 获取当前角度的投影projection = sinogram(i, :);% 滤波投影filtered_proj = filter_projections(projection, filter);% 反投影到体积recon_volume = backproject_3d(recon_volume, filtered_proj, angle, ...volume_size, voxel_size, ...source_to_detector, source_to_object, ...detector_pixel_size);end% 归一化重建体积recon_volume = recon_volume / num_angles;
endfunction filter = create_filter(filter_type, detector_size, detector_pixel_size)% 创建重建滤波器% 创建频率轴freq = linspace(-1, 1, detector_size) * (1/(2*detector_pixel_size));% 创建Ram-Lak滤波器 (斜坡滤波器)ramp_filter = abs(freq);% 根据滤波器类型应用窗函数switch lower(filter_type)case 'ram-lak'filter = ramp_filter;case 'shepp-logan'% Shepp-Logan滤波器: sinc窗斜坡滤波器sinc_window = sinc(freq / (2*max(freq)));filter = ramp_filter .* sinc_window;case 'cosine'% Cosine滤波器cosine_window = cos(pi * freq / (2*max(freq)));filter = ramp_filter .* cosine_window;case 'hamming'% Hamming窗滤波器hamming_window = 0.54 + 0.46 * cos(pi * freq / max(freq));filter = ramp_filter .* hamming_window;otherwiseerror('未知滤波器类型: %s', filter_type);end% 确保滤波器为行向量filter = filter(:)';
endfunction filtered_proj = filter_projections(projection, filter)% 滤波投影数据% 傅里叶变换proj_fft = fft(projection);% 应用滤波器filtered_fft = proj_fft .* filter;% 逆傅里叶变换filtered_proj = real(ifft(filtered_fft));
endfunction volume = backproject_3d(volume, projection, angle, volume_size, voxel_size, ...source_to_detector, source_to_object, detector_pixel_size)% 三维反投影% 计算几何参数magnification = source_to_detector / source_to_object;detector_length = detector_size * detector_pixel_size;% 创建体积坐标网格[x, y, z] = meshgrid(1:volume_size, 1:volume_size, 1:volume_size);% 转换为物理坐标 (以体积中心为原点)x = (x - volume_size/2) * voxel_size;y = (y - volume_size/2) * voxel_size;z = (z - volume_size/2) * voxel_size;% 旋转坐标系统以匹配当前角度theta = deg2rad(angle);rot_x = x * cos(theta) - y * sin(theta);rot_y = x * sin(theta) + y * cos(theta);rot_z = z; % Z轴不变% 计算探测器上的位置u = (source_to_detector * rot_x) ./ (source_to_object + rot_y);% 转换为探测器像素索引detector_coords = linspace(-detector_length/2, detector_length/2, length(projection));pixel_indices = interp1(detector_coords, 1:length(projection), u, 'linear', 'extrap');% 确保索引在有效范围内pixel_indices = max(1, min(length(projection), pixel_indices));% 插值投影值proj_values = interp1(1:length(projection), projection, pixel_indices, 'linear', 0);% 添加到体积中 (考虑几何权重)weight = source_to_detector^2 ./ (source_to_object + rot_y).^2;volume = volume + proj_values .* weight;
endfunction recon_volume = sart_reconstruct_3d(sinogram, angles, volume_size, ...source_to_detector, source_to_object, ...detector_pixel_size, voxel_size, num_iterations)% 三维同步代数重建技术 (SART)num_angles = length(angles);recon_volume = zeros(volume_size, volume_size, volume_size);% 预处理正弦图sinogram = max(sinogram, 0); % 确保非负% 迭代重建for iter = 1:num_iterationsfprintf('SART迭代 %d/%d\n', iter, num_iterations);% 对每个角度进行更新for i = 1:num_anglesangle = angles(i);% 计算当前角度的正向投影proj_estimate = forward_project_3d_single(recon_volume, angle, ...detector_size, source_to_detector, ...source_to_object, detector_pixel_size, ...voxel_size);% 计算投影误差proj_error = sinogram(i, :) - proj_estimate;% 反投影误差到体积空间error_volume = backproject_3d_single(zeros(size(recon_volume)), proj_error, angle, ...volume_size, voxel_size, ...source_to_detector, source_to_object, ...detector_pixel_size);% 计算权重体积 (用于归一化)weight_volume = backproject_3d_single(zeros(size(recon_volume)), ones(size(proj_error)), angle, ...volume_size, voxel_size, ...source_to_detector, source_to_object, ...detector_pixel_size);% 避免除以零weight_volume(weight_volume == 0) = 1;% 更新重建体积recon_volume = recon_volume + error_volume ./ weight_volume;% 施加非负约束recon_volume = max(recon_volume, 0);endend
endfunction projection = forward_project_3d_single(volume, angle, detector_size, ...source_to_detector, source_to_object, ...detector_pixel_size, voxel_size)% 单角度三维正向投影% 计算几何参数magnification = source_to_detector / source_to_object;detector_length = detector_size * detector_pixel_size;% 旋转体积rotated_volume = imrotate3(volume, angle, [0 0 1], 'crop', 'linear');% 计算投影 (沿X轴积分)projection = squeeze(sum(rotated_volume, 1));% 考虑几何放大projection = imresize(projection, [detector_size, detector_size] * magnification);% 提取中心行作为投影center_row = round(size(projection, 1)/2);projection = projection(center_row, :);
endfunction volume = backproject_3d_single(volume, projection, angle, volume_size, voxel_size, ...source_to_detector, source_to_object, detector_pixel_size)% 单角度三维反投影% 计算几何参数magnification = source_to_detector / source_to_object;detector_length = detector_size * detector_pixel_size;% 创建体积坐标网格[x, y, z] = meshgrid(1:volume_size, 1:volume_size, 1:volume_size);% 转换为物理坐标 (以体积中心为原点)x = (x - volume_size/2) * voxel_size;y = (y - volume_size/2) * voxel_size;z = (z - volume_size/2) * voxel_size;% 旋转坐标系统以匹配当前角度theta = deg2rad(angle);rot_x = x * cos(theta) - y * sin(theta);rot_y = x * sin(theta) + y * cos(theta);rot_z = z; % Z轴不变% 计算探测器上的位置u = (source_to_detector * rot_x) ./ (source_to_object + rot_y);% 转换为探测器像素索引detector_coords = linspace(-detector_length/2, detector_length/2, length(projection));pixel_indices = interp1(detector_coords, 1:length(projection), u, 'linear', 'extrap');% 确保索引在有效范围内pixel_indices = max(1, min(length(projection), pixel_indices));% 插值投影值proj_values = interp1(1:length(projection), projection, pixel_indices, 'linear', 0);% 添加到体积中 (考虑几何权重)weight = source_to_detector^2 ./ (source_to_object + rot_y).^2;volume = volume + proj_values .* weight;
endfunction rmse = compute_rmse(original, reconstructed)% 计算均方根误差diff = original - reconstructed;rmse = sqrt(mean(diff(:).^2));
endfunction ssim_val = compute_ssim_3d(original, reconstructed)% 计算三维结构相似性指数% 简化版本：计算所有切片的平均SSIMnum_slices = size(original, 3);ssim_vals = zeros(1, num_slices);for i = 1:num_slicesssim_vals(i) = ssim(squeeze(original(:, :, i)), squeeze(reconstructed(:, :, i)));endssim_val = mean(ssim_vals);
end