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 =