主程序:
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 =