✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,代码获取、论文复现及科研仿真合作可私信。
🍎个人主页:Matlab科研工作室
🍊个人信条:格物致知。
更多Matlab完整代码及仿真定制内容点击👇
🔥 内容介绍
Siddon 投影是一种用于计算机断层扫描 (CT) 成像的投影算法。它是一种快速且内存高效的算法,广泛用于医学成像和其他领域。本文将深入探讨基于 Siddon 投影的图像重建过程。
Siddon 投影算法
Siddon 投影算法基于以下原理:
-
射线穿过对象时,其强度会随着对象内部的衰减而减弱。
-
衰减量与射线路径上的对象密度成正比。
算法通过将对象离散化为体素网格来工作。对于每个体素,算法计算射线穿过体素的路径长度。然后,它根据路径长度和体素密度计算射线的衰减。最后,算法将所有射线的衰减值累加起来,得到投影图像。
图像重建
基于 Siddon 投影的图像重建涉及以下步骤:
-
**投影采集:**从不同角度采集对象的投影图像。
-
**反投影:**将每个投影图像反投影到图像空间中。
-
**滤波回投影:**对反投影图像进行滤波,以去除噪声和伪影。
-
**重建:**将滤波后的反投影图像相加,得到最终的重建图像。
优点
基于 Siddon 投影的图像重建具有以下优点:
-
**快速:**算法计算效率高,可以快速重建图像。
-
**内存高效:**算法只需要存储投影图像,不需要存储整个体素网格。
-
**易于实现:**算法的实现相对简单。
缺点
基于 Siddon 投影的图像重建也有一些缺点:
-
**伪影:**算法容易产生伪影,例如条纹和环状伪影。
-
**分辨率有限:**算法的分辨率受到投影图像分辨率的限制。
-
**噪声敏感:**算法对噪声敏感,这可能会导致图像质量下降。
应用
基于 Siddon 投影的图像重建广泛用于以下领域:
-
**医学成像:**CT 扫描、正电子发射断层扫描 (PET) 和单光子发射断层扫描 (SPECT)
-
**工业成像:**无损检测、材料表征
-
**科学成像:**显微成像、天文成像
结论
基于 Siddon 投影的图像重建是一种快速、内存高效且易于实现的算法。它广泛用于医学成像和其他领域。虽然算法有一些缺点,但其优点使其成为图像重建任务的宝贵工具。
📣 部分代码
function [p,ellipse]=phantom3dAniso(varargin)
%PHANTOM3D Three-dimensional analogue of MATLAB Shepp-Logan phantom
% P = PHANTOM3D(DEF,N) generates a 3D head phantom that can
% be used to test 3-D reconstruction algorithms.
%
% DEF is a string that specifies the type of head phantom to generate.
% Valid values are:
%
% 'Shepp-Logan' A test image used widely by researchers in
% tomography
% 'Modified Shepp-Logan' (default) A variant of the Shepp-Logan phantom
% in which the contrast is improved for better
% visual perception.
% 'yu-ye-wang' Another version of the modified Shepp-Logan
% phantom from "Katsevich-Type Algorithms for
% Variable Radius Spiral Cone-BeamCT"
%
% N specifies the 3D grid size of P
% If N is a scalar, P will have isotropic size [N, N, N]
% If N is a 3-vector, P will have size [N(1) N(2) N(3)]
% If you omit the argument, N defaults to [64 64 64].
%
% P = PHANTOM3D(E,N) generates a user-defined phantom, where each row
% of the matrix E specifies an ellipsoid in the image. E has ten columns,
% with each column containing a different parameter for the ellipsoids:
%
% Column 1: A the additive intensity value of the ellipsoid
% Column 2: a the length of the x semi-axis of the ellipsoid
% Column 3: b the length of the y semi-axis of the ellipsoid
% Column 4: c the length of the z semi-axis of the ellipsoid
% Column 5: x0 the x-coordinate of the center of the ellipsoid
% Column 6: y0 the y-coordinate of the center of the ellipsoid
% Column 7: z0 the z-coordinate of the center of the ellipsoid
% Column 8: phi phi Euler angle (in degrees) (rotation about z-axis)
% Column 9: theta theta Euler angle (in degrees) (rotation about x-axis)
% Column 10: psi psi Euler angle (in degrees) (rotation about z-axis)
%
% For purposes of generating the phantom, the domains for the x-, y-, and
% z-axes span [-1,1]. Columns 2 through 7 must be specified in terms
% of this range.
%
% [P,E] = PHANTOM3D(...) returns the matrix E used to generate the phantom.
%
% Class Support
% -------------
% All inputs must be of class double. All outputs are of class double.
%
% Remarks
% -------
% For any given voxel in the output image, the voxel's value is equal to the
% sum of the additive intensity values of all ellipsoids that the voxel is a
% part of. If a voxel is not part of any ellipsoid, its value is 0.
%
% The additive intensity value A for an ellipsoid can be positive or negative;
% if it is negative, the ellipsoid will be darker than the surrounding pixels.
% Note that, depending on the values of A, some voxels may have values outside
% the range [0,1].
%
% Example
% -------
% ph = phantom3d(128);
% figure, imshow(squeeze(ph(64,:,:)))
%
[ellipse,n] = parse_inputs(varargin{:});
nx = n(1); ny = n(2); nz = n(3);
p = zeros([nx ny nz]);
rngx = ( (0:nx-1)-(nx-1)/2 ) / ((nx-1)/2);
rngy = ( (0:ny-1)-(ny-1)/2 ) / ((ny-1)/2);
rngz = ( (0:nz-1)-(nz-1)/2 ) / ((nz-1)/2);
% PJB: Note the swap of the x and y with meshgrid parameters.
%[x,y,z] = meshgrid(rngx,rngy,rngz);
[x,y,z] = meshgrid(rngy,rngx,rngz);
x=single(x);y=single(y);z=single(z);
coord = [flatten(single(x)); flatten(single(y)); flatten(single(z))];
clear x y z;
p = flatten(p);
for k = 1:size(ellipse,1)
A = ellipse(k,1); % Amplitude change for this ellipsoid
function e = yu_ye_wang
%
% Yu H, Ye Y, Wang G, Katsevich-Type Algorithms for Variable Radius Spiral Cone-Beam CT
%
% A a b c x0 y0 z0 phi theta psi
% -----------------------------------------------------------------
e = [ 1 .6900 .920 .900 0 0 0 0 0 0
-.8 .6624 .874 .880 0 0 0 0 0 0
-.2 .4100 .160 .210 -.22 0 -.25 108 0 0
-.2 .3100 .110 .220 .22 0 -.25 72 0 0
.2 .2100 .250 .500 0 .35 -.25 0 0 0
.2 .0460 .046 .046 0 .1 -.25 0 0 0
.1 .0460 .023 .020 -.08 -.65 -.25 0 0 0
.1 .0460 .023 .020 .06 -.65 -.25 90 0 0
.2 .0560 .040 .100 .06 -.105 .625 90 0 0
-.2 .0560 .056 .100 0 .100 .625 0 0 0 ];
return;
⛳️ 运行结果
🔗 参考文献
[1] 陈涵天,郑建奎,崔索超,等.一种基于三维角锥的光束多程传输系统:202311366503[P][2024-03-08].
[2] 李志勇,王学涛,王笑宇,等.基于扩散模型的牙科锥束CT金属伪影校正方法及装置:202311274241[P][2024-03-08].
[3] 包宜骏,陶思玮,匡翠方,等.X射线锥束CT位姿矫正方法,系统及装置:202311684036[P][2024-03-08].
🎈 部分理论引用网络文献,若有侵权联系博主删除
🎁 关注我领取海量matlab电子书和数学建模资料
👇 私信完整代码和数据获取及论文数模仿真定制
1 各类智能优化算法改进及应用
生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱船配载优化、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化
2 机器学习和深度学习方面
2.1 bp时序、回归预测和分类
2.2 ENS声神经网络时序、回归预测和分类
2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类
2.4 CNN/TCN卷积神经网络系列时序、回归预测和分类
2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类
2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类
2.7 ELMAN递归神经网络时序、回归\预测和分类
2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类
2.9 RBF径向基神经网络时序、回归预测和分类