主程序：
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 = (-N/2-2:0.01:N/2+2)*delta;
yy = t(ii);
plot(xx, yy, 'b');
hold on;
end
%######################################
%如果超出网格范围，直接计算下一条
if t(ii)>= N/2*delta || t(ii)<= -N/2*delta
continue;
end
kout = N*ceil(N/2-t(ii)/delta); %出射点网格编号
kk = (kout-(N-1)):1:kout;
u(1:N) = kk;
v(1:N) = ones(1,N)*delta;
%===========投影角度等于0时==========%
elseif th == 0
%##############################
if N<=10 && length(theta)<=5
yy = (-N/2-2:0.01:N/2+2)*delta;
xx = t(ii);
plot(xx, yy, 'b');
hold on;
end
%##############################
%如果超出网格范围，直接计算下一条
if t(ii)>= N/2*delta || t(ii)<= -N/2*delta
continue;
end
kin = ceil(N/2+t(ii)/delta); %入射点网格编号
kk = kin:N:(kin+N*(N-1));
u(1:N) = kk;
v(1:N) = ones(1,N)*delta;
%====当90<th<180时，投影射线的角度（以x正半轴为起点）为th-90===
%====当 0<th<90 时，与投影射线关于y轴对称的射线的角度为90-th===
else
if th>90
th_temp = th-90;
elseif th<90
th_temp = 90-th;
end
th_temp = th_temp*pi/180; %角度转化为弧度
b = t/cos(th_temp);
m = tan(th_temp); %确定射线y=mx+b的m和b
y1d = -N/2*delta*m + b(ii); %入射点y坐标
y2d = N/2*delta*m + b(ii); %出射点y坐标
%################################
if N<=10 && length(theta)<=5
xx = (-N/2-2:0.01:N/2+2)*delta;
if th<90
yy = -m*xx+b(ii);
elseif th>90
yy = m*xx+b(ii);
end
plot(xx, yy, 'b');
hold on;
end
%###################################
%如果超出网格范围，直接计算下一条
if (y1d<-N/2*delta && y2d<-N/2*delta) || (y1d>N/2*delta && y2d>N/2*delta) %感觉书里误加了一个负号
continue;
end
%确定入射点（xin,yin）、出射点（xout,yout）及参数d1
if y1d<=N/2*delta && y1d>=-N/2*delta && y2d>N/2*delta
yin = y1d;
d1 = yin-floor(yin/delta)*delta; %floor向下取整
kin = N*floor(N/2-yin/delta)+1;
yout = N/2*delta;
xout = (yout-b(ii))/m;
kout = ceil(xout/delta)+N/2; %ceil向上取整
elseif y1d<=N/2*delta && y1d>=-N/2*delta && y2d>=-N/2*delta && y2d<=N/2*delta
yin = y1d;
d1 = yin-floor(yin/delta)*delta;
kin = N*floor(N/2-yin/delta)+1;
yout = y2d;
kout = N*floor(N/2-yout/delta)+N;
elseif y1d<-N/2*delta && y2d>N/2*delta
yin = -N/2*delta;
xin = (yin-b(ii))/m;
d1 = N/2*delta+(floor(xin/delta)*delta*m+b(ii));
kin = N*(N-1)+N/2+ceil(xin/delta);
yout = N/2*delta;
xout = (yout-b(ii))/m;
kout = ceil(xout/delta)+N/2;
elseif y1d<-N/2*delta && y2d>=-N/2*delta && y2d<=N/2*delta
yin = -N/2*delta;
xin = (yin-b(ii))/m;
d1 = N/2*delta+(floor(xin/delta)*delta*m+b(ii));
kin = N*(N-1)+N/2+ceil(xin/delta);
yout = y2d;
kout = N*floor(N/2-yout/delta)+N;
else
continue;
end
%计算射线i穿过像素的编号和长度
k = kin;
c = 0;
d2 = d1+m*delta;
while k>=1 && k<=N2
c = c+1;
if d1>=0 && d2>delta
u(c) = k;
v(c) = (delta-d1)*sqrt(m^2+1)/m;
if k>N && k~=kout
k = k-N;
d1 = d1-delta;
d2 = d1 +m*delta;
else
break
end
elseif d1>=0 && d2==delta
u(c) = k;
v(c) = delta*sqrt(m^2+1);
if k>N && k~=kout
k = k-N+1;
d1 = 0;
d2 = d1 +m*delta;
else
break
end
elseif d1>=0 && d2<delta
u(c) = k;
v(c) = delta*sqrt(m^2+1);
if k~=kout
k = k+1;
d1 = d2;
d2 = d1 +m*delta;
else
break
end
elseif d1<=0 && d2>=0 && d2<=delta
u(c) = k;
v(c) = d2*sqrt(m^2+1)/m;
if k~=kout
k = k+1;
d1 = d2;
d2 = d1 +m*delta;
else
break
end
elseif d1<=0 && d2>delta
u(c) = k;
v(c) = delta*sqrt(m^2+1)/m;
if k>N && k~=kout
k = k-N;
d1 = d1-delta;
d2 = d1 +m*delta;
else
break
end
end
end
%%如果投影角度小于90，利用投影射线关于y轴的对称性重新计算编号
if th<90
u_temp = zeros(1, 2*N);
if any(u)==0
continue %如果不经过任何网格，直接计算下一条
end
ind = u>0;
for k=1:length(u(ind))
r = rem(u(k), N);
if r==0
u_temp(k) = u(k)-N+1;
else
u_temp(k) = u(k)-2*r+N+1;
end
end
u = u_temp;
end
end
W_ind((jj-1)*P_num+ii, :) = u;
W_dat((jj-1)*P_num+ii, :) = v;
end
end平行束投影：
function P = medfuncParallelBeamForwardProjection(theta, N, N_d)
% Parallel beam forward projection(只针对头模型的解析算法)
%--------------------------
% 输入参数：
% theta：投影角度
% N：图像大小
% N_d：探测器通道个数
%--------------------------
% 输出参数：
% P：投影矩阵（N_d*theta_num）
%====================================
%
shep = [1 .69 .92 0 0 0
-.8 .6624 .8740 0 -.0184 0
-.2 .1100 .3100 .22 0 -18
-.2 .1600 .4100 -.22 0 18
.1 .2100 .2500 0 .35 0
.1 .0460 .0460 0 .1 0
.1 .0460 .0460 0 -.1 0
.1 .0460 .0230 -.08 -.605 0
.1 .0230 .0230 0 -.606 0
.1 .0230 .0460 .06 -.605 0];
theta_num = length(theta);
P = zeros(N_d, theta_num);
rho = shep(:,1).'; %椭圆对应密度
ae = 0.5*N*shep(:,2).'; %椭圆短半轴
be = 0.5*N*shep(:,3).'; %椭圆长半轴
xe = 0.5*N*shep(:,4).'; %椭圆中心x坐标
ye = 0.5*N*shep(:,5).'; %椭圆中心y坐标
alpha = shep(:,6).'; %椭圆旋转角度
alpha = alpha*pi/180;
theta = theta*pi/180; %角度换算为弧度
TT = -(N_d-1)/2:(N_d-1)/2;
for k1 = 1:theta_num
P_theta = zeros(1, N_d);
for k2 = 1:max(size(xe))
a = (ae(k2)*cos(theta(k1)-alpha(k2)))^2+(be(k2)*sin(theta(k1)-alpha(k2)))^2;
temp = a-(TT-xe(k2)*cos(theta(k1))-ye(k2)*sin(theta(k1))).^2;
ind = temp>0; %根号内值需为非负
P_theta(ind) = P_theta(ind)+rho(k2)*(2*ae(k2)*be(k2)*sqrt(temp(ind)))./a;
end
P(:, k1) = P_theta;
end
endPOCS-TVM算法：
function F = medfuncPOCS_TVM(N, W_ind, W_dat, P, irt_num, F0, num_TVM, lambda, alpha)
% POCS-TVM算法
%-----------------------------
%输入参数：
% N：图像大小
% W_ind：射线所穿过网格的编号，M*(2*N)
% w_dat：射线所穿过网格的长度，M*(2*N)
% P：投影数据
% irt_num：算法的总迭代次数
% F0：初始图像向量
% num_TVM：全变分最小化（TVM）过程的迭代次数，默认值为6
% lambda：松弛因子，默认值为0.25
% alpha：调节因子，默认值为0.2
%------------------------------
%输出参数：
% F：重建图像
%==============================
if nargin<6, F0 = zeros(N^2, 1); end
if nargin<7, num_TVM = 6; end
if nargin<8, lambda = 0.25; end
if nargin<9, alpha = 0.2; end
F = F0;
N2 = N^2;
[P_num, theta_num] = size(P); % P_num:每个投影的射线条数；theta_num:投影个数
e = 0.00000001;
k1 = 0; %循环控制变量
while(k1<irt_num)
TEMP1 = F;
%----------------------------
for j = 1:theta_num
for i = 1:1:P_num
%取得一条射线所穿过的网格编号和长度 W行向量
u = W_ind((j-1)*P_num+i,:); %编号
v = W_dat((j-1)*P_num+i,:); %长度
%如果射线不经过任何像素，不做计算
if any(u) == 0;
continue;
end
%恢复投影矩阵中与这一条射线对应的行向量w
w = zeros(1,N2);
ind = u>0;
w(u(ind)) = v(ind);
%图像进行一次ART迭代
PP = w*F; %前向投影
C = (P(i,j)-PP)/sum(w.^2)*w'; %修正项
F = F+lambda*C;
end
end
F(F<0) = 0; %非负约束
%------------------------------
d = sqrt((TEMP1-F)'*(TEMP1-F)); %增量因子
k2 = 0;
while(k2<num_TVM)
G = zeros(N2, 1);
for i=2:N-1
for j=2:N-1
G((j-1)*N+i) = (2*F((j-1)*N+i)-F((j-1)*N+i-1)-F((j-2)*N+i))/sqrt(e+(F((j-1)*N+i)-F((j-1)*N+i-1))^2+(F((j-1)*N+i)-F((j-2)*N+i))^2)-...
F(j*N+i)-F((j-1)*N+i)/sqrt(e+(F(j*N+i)-F((j-1)*N+i))^2+(F(j*N+i)-F(j*N+i-1))^2)-...
F((j-1)*N+i+1)-F((j-1)*N+i)/sqrt(e+(F((j-1)*N+i+1)-F((j-1)*N+i))^2+(F((j-1)*N+i+1)-F((j-2)*N+i+1))^2);
end
end
G = G/sqrt(G'*G);
F = F-alpha*d*G;
k2 = k2+1;
end
k1 = k1+1;
end
end理论略解以后会传（大概不会了，我好懒），大家也可以去读黄力宇老师出的《医学断层图像重建仿真实验》，写的真的很好！！文章代码全来自此书，根据我的理解，对一两处印刷错误进行了改进。
欢迎大家纠错与讨论(〃'▽'〃)
1268
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
