global img t
img=avsread('phtmebl.fld');

ee=avsread('proj.fld');
pideal=avsread('../huber/proj.fld');

x=avsread('../huber/hc/out113.fld');
%load ../data4/myout3; x=myout3;
%x=avsread('../huber/lc/out113.fld');
%load ../myout3; x=myout3;
%x=avsread('../huber/lc/fbp.fld');
xmean=mean(x,3); avswrite('xmean.fld',xmean);
!i project pmean.fld xmean.fld ../huber/my.wtf
!rm xmean.fld
pmean=avsread('pmean.fld');
!rm pmean.fld

% For Huber hc 11.3 ***** Use this
gf1=exp(-((1:21)'-11).^2/10); gf2=exp(-((1:21)'-11).^2/5.5);
% For myout3 hc ***** Use this
%gf1=exp(-((1:21)'-11).^2/10); gf2=exp(-((1:21)'-11).^2/5.5);
% For Huber lc 11.3 ***** Use this
%gf1=exp(-((1:21)'-11).^2/18); gf2=exp(-((1:21)'-11).^2/11.4);
% For myout3 lc ***** Use this
%gf1=exp(-((1:21)'-11).^2/18); gf2=exp(-((1:21)'-11).^2/11);
% For FBP hc
%gf1=exp(-((1:21)'-11).^2/10); gf2=exp(-((1:21)'-11).^2/3);
% For FBP lc
%gf1=exp(-((1:21)'-11).^2/10); gf2=exp(-((1:21)'-11).^2/5);

gf1=gf1/sum(gf1); gf2=gf2/sum(gf2);

pout=filter2(gf1,ee.*exp(-pideal))./filter2(gf2,exp(-pmean));
avswrite('projn.fld',pout);
!rm fbp.fld
!i fbp dsc fbp.fld projn.fld ~/old/real1/921,thorax,attn,2.dsc boxcar,1
!rm projn.fld
ex=avsread('fbp.fld'); image(256-ex'/4e4*256);

findreso(ex)
