% calculate Boundary error

warning off
n=50; N=128; ideal=avsread('huber/phtm.fld')';
%x1=ideal(1:N-1,1:N-1); x2=ideal(2:N,1:N-1); x3=ideal(1:N-1,2:N);
%x4=ideal(2:N,2:N); bb=max(max(x1,x2),max(x3,x4))-min(min(x1,x2),min(x3,x4));
%[a1x,a1y]=find(abs(bb(:,1:58)-10000)<0.01); a1x=a1x+0.5; a1y=a1y+0.5;
%[a2x,a2y]=find(abs(bb(:,80:N-1)-10000)<0.01); a2y=a2y+79+0.5; a2x=a2x+0.5;
%[a3x,a3y]=find(abs(bb-22000)<0.01|abs(bb-16400)<0.01);a3x=a3x+0.5; a3y=a3y+0.5;
%[a4x,a4y]=find(abs(bb(:,59:75)-10000)<0.01); a4y=a4y+58+0.5; a4x=a4x+0.5;
load idealbdry
n1=zeros(length(a1x),2); n2=zeros(length(a2x),2); 
n3=zeros(length(a3x),2); n4=zeros(length(a4x),2); 

%figure(1),image(ideal/200); colormap(gray(256));hold on;
%for j=1:length(a4x),plot(a4y(j),a4x(j),'r.');end
%for j=1:length(a3x),plot(a3y(j),a3x(j),'r.');end
%for j=1:length(a2x),plot(a2y(j),a2x(j),'r.');end
%for j=1:length(a1x),plot(a1y(j),a1x(j),'r.');end
%hold off; %axis([20 80 30 90]);
%axis('square');drawnow;

gennv
mag1=zeros(length(a1x),n);
mag2=zeros(length(a2x),n);
mag3=zeros(length(a3x),n);
mag4=zeros(length(a4x),n);
ee1x=zeros(length(a1x),n);
ee1y=zeros(length(a1y),n);
ee2x=zeros(length(a2x),n);
ee2y=zeros(length(a2y),n);
ee3x=zeros(length(a3x),n);
ee3y=zeros(length(a3y),n);
ee4x=zeros(length(a4x),n);
ee4y=zeros(length(a4y),n);
e1x=zeros(length(a1x),1);
e1y=zeros(length(a1y),1);
e2x=zeros(length(a2x),1);
e2y=zeros(length(a2y),1);
e3x=zeros(length(a3x),1);
e3y=zeros(length(a3y),1);
e4x=zeros(length(a4x),1);
e4y=zeros(length(a4y),1);

load datalc/c1;

for i=1:n,
	j=i;
	[ee1x(:,i),ee1y(:,i)]=behelp(a1x,a1y,l1(j,ch),c1(:,:,j,ch));
       	[ee2x(:,i),ee2y(:,i)]=behelp(a2x,a2y,l2(j,ch),c2(:,:,j,ch));
       	[ee3x(:,i),ee3y(:,i)]=behelp(a3x,a3y,l3(j,ch),c3(:,:,j,ch));
       	[ee4x(:,i),ee4y(:,i)]=behelp(a4x,a4y,l4(j,ch),c4(:,:,j,ch));
	mag1(:,i)=ee1x(:,i).*n1(:,1)+ee1y(:,i).*n1(:,2);
	ee1x(:,i)=mag1(:,i).*n1(:,1); ee1y(:,i)=mag1(:,i).*n1(:,2);
	mag2(:,i)=ee2x(:,i).*n2(:,1)+ee2y(:,i).*n2(:,2);
	ee2x(:,i)=mag2(:,i).*n2(:,1); ee2y(:,i)=mag2(:,i).*n2(:,2);
	mag3(:,i)=ee3x(:,i).*n3(:,1)+ee3y(:,i).*n3(:,2);
	ee3x(:,i)=mag3(:,i).*n3(:,1); ee3y(:,i)=mag3(:,i).*n3(:,2);
	mag4(:,i)=ee4x(:,i).*n4(:,1)+ee4y(:,i).*n4(:,2);
	ee4x(:,i)=mag4(:,i).*n4(:,1); ee4y(:,i)=mag4(:,i).*n4(:,2);
end

e1x=mean(ee1x,2); e1y=mean(ee1y,2);
e2x=mean(ee2x,2); e2y=mean(ee2y,2);
e3x=mean(ee3x,2); e3y=mean(ee3y,2);
e4x=mean(ee4x,2); e4y=mean(ee4y,2);
bimgx=zeros(N,N); bimgy=zeros(N,N);[xxm,yym]=meshgrid(1:128,1:128);
a1x=floor(a1x); a1y=floor(a1y);
a2x=floor(a2x); a2y=floor(a2y);
a3x=floor(a3x); a3y=floor(a3y);
a4x=floor(a4x); a4y=floor(a4y);
for i=1:length(e1x),bimgx(a1x(i),a1y(i))=e1x(i);bimgy(a1x(i),a1y(i))=e1y(i);end
for i=1:length(e2x),bimgx(a2x(i),a2y(i))=e2x(i);bimgy(a2x(i),a2y(i))=e2y(i);end
for i=1:length(e3x),bimgx(a3x(i),a3y(i))=e3x(i);bimgy(a3x(i),a3y(i))=e3y(i);end
for i=1:length(e4x),bimgx(a4x(i),a4y(i))=e4x(i);bimgy(a4x(i),a4y(i))=e4y(i);end
quiver(xxm,yym,bimgy,bimgx,5); zoom on; axis([1 128 1 128]); axis ij;

ttt=[[mean(mean(mag1)),mean(sqrt(e1x.^2+e1y.^2)),std(sqrt(e1x.^2+e1y.^2)),...
	mean(std((sqrt(ee1x.^2+ee1y.^2))'))]
	[mean(mean(mag2)),mean(sqrt(e2x.^2+e2y.^2)),std(sqrt(e2x.^2+e2y.^2)),...
        mean(std((sqrt(ee2x.^2+ee2y.^2))'))]
	[mean(mean(mag3)),mean(sqrt(e3x.^2+e3y.^2)),std(sqrt(e3x.^2+e3y.^2))...
        mean(std((sqrt(ee3x.^2+ee3y.^2))'))]
	[mean(mean(mag4)),mean(sqrt(e4x.^2+e4y.^2)),std(sqrt(e4x.^2+e4y.^2)),...
        mean(std((sqrt(ee4x.^2+ee4y.^2))'))]];
disp(ch);
disp(ttt);


%for i=1:n,clf;image(ideal/150); colormap(gray(256)); axis('square');hold on;i
%for j=1:l1(i,ch), plot(c1(j,2,i,ch),c1(j,1,i,ch),'r.'); end
%for j=1:l2(i,ch), plot(c2(j,2,i,ch),c2(j,1,i,ch),'r.'); end
%for j=1:l3(i,ch), plot(c3(j,2,i,ch),c3(j,1,i,ch),'r.'); end
%for j=1:l4(i,ch), plot(c4(j,2,i,ch),c4(j,1,i,ch),'r.'); end
%drawnow; hold off; pause; end

