format short e
x=avsread('../huber/phtm.fld'); xx=(1:128); [X Y]=meshgrid(xx,xx);

r1=and((x==12000),Y<60); r2=and((x==12000),Y>60);
r3=zeros(128,128); r3(85:89,50:54)=1;
%r3=(x==32000); r4=or((x==22000),(x==26000));;
r4=(X>70 & X<72 & Y<72 & Y>70);

xx=(-15:15); [X Y]=meshgrid(xx,xx); R2=X.^2+Y.^2;
sigma=4; Gf=1/sqrt(2*pi)/sigma*exp(-R2/2/sigma/sigma); Gf=Gf/sum(Gf(:));

y=avsread('../huber/phtmtbl.fld');
tt=y.*r1; t1=mean(tt(find(tt~=0)));
tt=y.*r2; t2=mean(tt(find(tt~=0)));
tt=y.*r3; t3=mean(tt(find(tt~=0)));
tt=y.*r4; t4=mean(tt(find(tt~=0)));

N=50;
%load ../data4/myout3; z=myout3;
z=avsread('../huber/lc/ouv13.40.fld');
s1=zeros(1,N); s2=zeros(1,N); s3=zeros(1,N); s4=zeros(1,N); 
for i=1:50,
	%z(:,:,i)=filter2(Gf,z(:,:,i));
	ss=z(:,:,i).*r1; s1(i)=mean(ss(find(ss~=0)));
	ss=z(:,:,i).*r2; s2(i)=mean(ss(find(ss~=0)));
	ss=z(:,:,i).*r3; s3(i)=mean(ss(find(ss~=0)));
	ss=z(:,:,i).*r4; s4(i)=mean(ss(find(ss~=0)));
end
disp([mean(s1),mean(s2),mean(s3),std(s1),std(s2),std(s3)]);
%disp([mean(s1),mean(s2)])

