CT算法，radon变换基于MATLAB的CT算法，radon变换的三维建模仿真

Posted 2022-03-30 fpga&matlab

1.软件版本

MATLAB2021a

2.本算法理论知识

1.输入：T（x,y,z）

使用stl读取函数完成T的导入工作

2.做Radon变换，得投影图：P

正常Radon变换即可。

3.对P：应用斜坡滤波器，截断负数 p+(i)

 

该式为打印成型原理，核心是CT算法，I为阈值函数，α吸光率，Ω旋转速率，Dc为成型阈值，暂不考虑这些参数。f（r，z）即是模型的几何形状。

则对于一个横截面z，有：

有 f(r)= [g](r),  f(r)即截面图形状，是指数型反向radon变换，[g](r)是光线投影的线积分（即检测到的投影强度）。

现在我们的算法是逆过来的，即已知模型几何形状f（r）,去求所需的投影强度g的过程，由下式给出：

此式为傅里叶中心切片定理，F-1为傅里叶逆变换，其中

为f()函数的二维傅里叶变换。

，是一斜坡滤波器，同时去除频域中频率k在阈值α以下的部分。

由投影强度图g（r，z），通过反向radon变换（即反向投影）可以得到几何图形。由于这一理论值存在负数，实际的投影强度却只能为正，这里采用Otsu（）方法迭代优化：

算法的优化部分见补充材料S12，S13,

具体的流程图（见补充材料figure.s6）如下:

1.对输入的截面信息做radon变换

2.滤波+截断（傅里叶频域）

3.反radon变换

4.迭代优化Otsu（）部分

5.输出投影所用图像

3.核心代码

clc;
clear
close all;
warning off;
addpath 'func\\'

%DLP
theta = 1;%设置1~180之间






load R3.mat

%这里，根据各截面投影图；然后反投影得到物体几何结构，非常消耗内存，我这里进行间隔抽样，降低计算机资源消耗
%如果电脑配置高，则直接设置1即可
Skip = 2;

V=ones(512/Skip,512/Skip,520/Skip);

idx=0;
for k = 1:Skip:520
    idx=idx+1;
    tmps = imresize(Pnewk,1/Skip);
    [X,Y]= find(tmps==0);
    for j = 1:length(X)
        v(X(j),Y(j),idx) = NaN;
    end
end

[x,y,z]=meshgrid(1:Skip:512,1:Skip:512,520:-Skip:1);
figure;
slice(x,y,z,v,[1:Skip:512],[1:Skip:512],[520:-Skip:1]) 
shading interp                             
colorbar                                   
                         
set(gca,'zdir','reverse');                 
grid on 
box on       
colormap(bone);
axis equal
view([135,15]);
camlight('infinite');
 



for i = SETS
    E   = P0i;
    F   = Pi;
    for j = 1:length(E(:,theta))
        E2(521-i,j)= [E(end+1-j,theta)];
        F2(521-i,j)= [F(end+1-j,theta)]; 
    end
end

C=[2,1000];
figure;
subplot(121)
imagesc(E2,C);title('E');
colormap('hot');

subplot(122)
imagesc(F2,C);title('F');
colormap('hot');

function level=func_ostu(I); 

I        = im2uint8(I(:));
num_bins = 256;
counts   = imhist(I,num_bins);
num_bins = numel(counts);
counts   = double( counts(:) );
p        = counts / sum(counts);
omega    = cumsum(p);
mu       = cumsum(p .* (1:num_bins)');
mu_t     = mu(end);

sigma_b_squared = (mu_t * omega - mu).^2 ./ (omega .* (1 - omega));
 
maxval          = max(sigma_b_squared);
isfinite_maxval = isfinite(maxval);

if isfinite_maxval
    idx  = mean(find(sigma_b_squared == maxval));
    level = (idx - 1) / (num_bins - 1);
else
    level = 0;
end
 
end

4.操作步骤与仿真结论

 

 

5.参考文献

 

A19-26

6.完整源码获得方式

方式1：微信或者QQ联系博主

方式2：订阅MATLAB/FPGA教程，免费获得教程案例以及任意2份完整源码