CT重建算法之统计迭代类算法;ML-EM算法（matlab的实现版本）

function [ W_ind, W_dat ] = medfuncSystemMatrix( theta, N, P_num, delta ) %UNTITLED Summary of this function goes here % Detailed explanation goes here N2 = N^2; M = length(theta)*P_num; W_ind = zeros(M, 2*N); W_dat = zeros(M, 2*N); 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'); hold on; plot(x,meshgrid(x,y),y,'k'); 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('输入角度有问题'); 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; 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; 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); y1d = -N/2*delta*m + b(ii); y2d = N/2*delta*m + b(ii); 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 continue; end if y1d>N/2*delta & y2d>-N/2*delta continue; end if y1d<=N/2*delta & y1d>=-N/2*delta & y2d>N/2*delta yin = y1d; d1 = yin - floor(yin/delta)*delta; kin = N*floor(N/2 - yin/delta) + 1; yout = N/2*delta; xout = (yout-b(ii))/m; kout = ceil(xout/delta) + N/2; 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 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 = -delta + d1; d2 = d1 + m*delta; else break; end end end 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