CT重建(POCS-TVM)matlab实现

 主程序:

clc;clear;close all;
N = 180;  %图像大小
N2 = N^2;
I = phantom(N);  %图像

theta = linspace(0,180,61);  %%算头算尾61个点
theta = theta(1:60);  %投影角度
P_num = 260;
P = medfuncParallelBeamForwardProjection(theta, N, P_num);  %投影数据
%P = radon(I, theta);

delta = 1;
[W_ind, W_dat] = medfuncSystemMatrix(theta, N, P_num, delta);  %投影矩阵

irt_num = 5;
F0 = zeros(N2, 1);
num_TVM = 4;
lambda = 0.25;
alpha = 0;
F = medfuncPOCS_TVM(N, W_ind, W_dat, P, irt_num, F0, num_TVM, lambda, alpha);
F = reshape(F, N, N)';  % F排成矩阵时,一行一行的来

figure(1);
imshow(I);xlabel('(a)180*180头模型图像');
figure(2);
imshow(F);xlabel('(b)POCS-TVM算法重建图像');

投影矩阵计算:

function [W_ind, W_dat] = medfuncSystemMatrix(theta, N, P_num, delta)
% W_fun 计算投影矩阵
%---------------------------------
% 输入参数:
% theta:(平行)投影角度,适用于0=<theta<180
% N:矩阵(图像)大小
% P_num:每个投影角度下的射线条数(探测器通道个数)
% delta:网格大小
%---------------------------------
% 输出参数:
% 以W_ind和W_dat表示的投影矩阵
% W_ind:存储投影射线所穿过网格的编号,M*(2*N)
% W_dat:存储投影射线所穿过网格的长度,M*(2*N)
%=====================================================
N2 = N^2;
M = length(theta)*P_num;  %投影射线总条数
W_ind = zeros(M,2*N);  %存放射线穿过网格的编号
W_dat = zeros(M,2*N);  %存放射线穿过网格的长度
% t_max = sqrt(2)*N*delta;
% t = linspace(-t_max/2, t_max/2, P_num);
t = (-(P_num-1)/2:(P_num-1)/2)*delta;  %探测器坐标
%####当图像较小、射线条数较少时,画出扫描结构图(网格图)####%
if N<=10 && length(theta)<=5
    x = (-N/2:1:N/2)*delta;
    y = (-N/2:1:N/2)*delta;
    plot(x, meshgrid(y,x), 'k','LineStyle','--');
    hold on;
    plot(meshgrid(x,y), y, 'k','LineStyle','--');
    hold on;
    axis([-N/2-5, N/2+5, -N/2-5, N/2+5]);
    text(0, -0.4*delta, '0');
end  
%##########################################
% ===============投影矩阵的计算================ %
for jj = 1:length(theta)  
    for ii = 1:1:P_num      %每次完成一条射线权因子向量的计算
        u = zeros(1, 2*N);  v = zeros(1, 2*N);
        th = theta(jj);
        if th>=180 || th<0
             error('输入角度必须在0~180之间');
        %==============投影角度等于90时(平行于x轴)=============%
        elseif th == 90
            %##########画出对应的射线图############
            if N<=10 && length(theta)<=5
                xx =
评论 6
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值