新建文件夹.zip_ART TV_ART 算法_MSE和SSIM_TV图像重建_art-tv

function [X_grad] = ART_TV_BIN_newYr_t(varargin) %sp0 源点坐标 %dp0 探测器横纵坐标 2*Ndetector % %grid 网格参数 %theta 投影角度 %pj 投影 %NART art次数 %NTV tv次数 %f 上一次的投影结果 %valid 有效探测器标记为1，坏探测器标记为0，针对坏探测器情形 pj=varargin{1}; SM=varargin{2}; imageSize=varargin{3}; X=varargin{4}; old_Yr=varargin{5}; % 新引入的 t = varargin{6}; s = varargin{7}; ND=256;%探测器个数 % X_grad = zeros(size(X)); % X_grad2 = X; Xsiz = [imageSize,imageSize,size(X,2)];%长 宽 能区 % Xcell = cell(1); % % X_grad = zeros(size(X)); % o = zeros(256,256); % % X_grad2 = zeros(size(X)); % for p = 1:Xsiz(3) % Xcell{p} = X(:, :, p); % % % Xcell_grad{p} = zeros(size(X(:, :, p))); % % % Xcell_grad2{p} = zeros(size(X(:, :, p))); % % % Xcell_grad{s} = X(:, :, s); % end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % for nart = 1:1 %ART迭代次数 % lastXcell = Xcell; X_grad = zeros(size(X)); % for t = 1:900 %每个角度 % % % % % for s = 1:ND %每个探测器，文章中的j=(t-1)*length（theta）+s if ~isempty(SM.isect{t}{s}) %若t，s有相交 % for p = 1:Xsiz(3) % Xcell{p} = X(:, :, p); % % end thsh = SM.table(t, s); %取出阈值 ll = length(SM.isect{t}{s}); %t，s相交的像素个数 f = zeros(ll,Xsiz(3)); %能区 相交数 f1 = zeros(ll,Xsiz(3)); % fnew = zeros(ll,Xsiz(3)); fnew2 = zeros(ll,Xsiz(3)); for k = 1:Xsiz(3) %每个能区 f(:,k) = X(SM.isect{t}{s},k); %取出每个能区下t，s有相交的投影图像值 f1(:,k) = (old_Yr(t,s,k)) * X(SM.isect{t}{s},k); % 新引入的，G(k,j)*f(k)，注意G(k,j)是数而f是图像向量（代码中只取出了f中不为0的值） end f(f<0) = 0; f1(f1<0) = 0; % if thsh == 1 %在动态低能范围内 f1对能区进行求和，即Gj(k)·f(k)。 % %由于ART算法要求图像是一维的（否则无法写出系统矩阵），这里先将G和f的 % %点积（一维向量）作为ART求解的图像，完成之后再进行高低能的反投影，最终解出f % f_L1 = f1(1, :);% 新修改的 % else f_L1 = sum(f1(:,1:thsh),2); % 新修改的 % end correctL = SM.lsect{t}{s}*(SM.lsect{t}{s}*f_L1 - pj(t, s, 1))/sum(SM.lsect{t}{s}.^2); % 新修改的 - ----> + %correctL为ART-TV算法中对fTV-DATA的更新项 的相反数 % f_L2 = f_L1 - correctL'; %即新的G·f f_L2t = correctL'; for i = 1:ll %在这里做低能的反投影，每个低能能区乘上了同一个常数，高能区乘上另一个常数。 if f_L1(i) == 0 if sum(old_Yr(t,s,1:thsh)) >0 % fnew(i,1:thsh) = f_L2(i) / sum(old_Yr(t,s,1:thsh)); % double(thsh) ; fnew2(i,1:thsh) = f_L2t(i) / sum(old_Yr(t,s,1:thsh)); % double(thsh) ; end else % fnew(i,1:thsh) = f_L2(i) / f_L1(i) * f(i,1:thsh); fnew2(i,1:thsh) = f(i,1:thsh)/ f_L1(i) * f_L2t(i); end end % if thsh+1 == Xsiz(3) % f_U1 = f1(Xsiz(3), :); % 新修改的 % else f_U1 = sum(f1(:,(thsh+1):Xsiz(3)),2); % 新修改的 % end correctU = SM.lsect{t}{s}*(SM.lsect{t}{s}*f_U1 - pj(t, s, 2))/sum(SM.lsect{t}{s}.^2); % 新修改的 - ----> + % f_U2 = f_U1 - correctU'; f_U2t = correctU'; for i = 1:ll if f_U1(i) == 0 if sum(old_Yr(t,s,(thsh+1):Xsiz(3))) >0 % fnew(i,(thsh+1):Xsiz(3)) = f_U2(i) / sum(old_Yr(t,s,(thsh+1):Xsiz(3))); % (Xsiz(3)-double(thsh)); fnew2(i,(thsh+1):Xsiz(3)) = f_U2t(i) / sum(old_Yr(t,s,(thsh+1):Xsiz(3))); % (Xsiz(3)-double(thsh)); end else % fnew(i,(thsh+1):Xsiz(3)) = f_U2(i) / f_U1(i) * f(i,(thsh+1):Xsiz(3)); fnew2(i,(thsh+1):Xsiz(3)) = f(i,(thsh+1):Xsiz(3)) / f_U1(i) * f_U2t(i); end end % X(SM.isect{t}{s},:) = f - fnew2; % X_grad(SM.isect{t}{s},:) = X_grad(SM.isect{t}{s},:) + fnew2; X_grad(SM.isect{t}{s},:) = fnew2; % end end % for p = 1:Xsiz(3) % X_grad(:, :, p)=Xcell{p}; % end % for p = 1:Xsiz(3) % X_grad2(:, :, p)=Xcell_grad{p}; % end % X = X - X_grad; % end % % end % for s = 1:Xsiz(3) % X(:, :, s)=Xcell{s}; % end % %以下为fTV-POS，与动态高低能无关 % for i = 1:Xsiz(3) % Xcell{i}(Xcell{i}<0)=0; % Xcell{i}(Xcell{i}>40)=40; % end % %以下部分为fTV-GRAD，与动态高低能无关。 % for i = 1:Xsiz(3) % f = Xcell{i}; % lastf = lastXcell{i}; % [s1,s2]=size(f); % eps=1e-8; % a=0.2; % d=norm(f(:)-lastf(:)); % % for ntv=1:NTV % %extend % fe=[zeros(1,s2+2);[zeros(s1,1) f zeros(s1,1)];zeros(1,s2+2)]; %把f框一圈0 % part1=((f-fe(1:s1,2:s2+1))+(f-fe(2:s1+1,1:s2)))... % ./sqrt((f-fe(1:s1,2:s2+1)).^2+(f-fe(2:s1+1,1:s2)).^2+eps); % part2=(-f+fe(3:s1+2,2:s2+1))... % ./sqrt((-f+fe(3:s1+2,2:s2+1)).^2+(fe(3:s1+2,2:s2+1)-fe(3:s1+2,1:s2)).^2+eps); % part3=(-f+fe(2:s1+1,3:s2+2))... % ./sqrt((-f+fe(2:s1+1,3:s2+2)).^2+(fe(2:s1+1,3:s2+2)-fe(1:s1,3:s2+2)).^2+eps); % v=part1-part2-part3; %每个像素与上下左右的TV % v=v/norm(v(:)); % f=f-a*d*v; % end; % % Xcell{i} = f; % end