function s=dipole(x,level)
%
% dipole Image Dipole Calculation function
%
% Example s=dipole(x,level);
%           x: image
%           level: image gray level, ex: level=2 for binary image
%

[row,col]=size(x);
s=zeros(2*row-1,2*col-1,level,level);
for dx=0:(row-1)
   for dy=0:(col-1)
      for m=0:(row-1)-dx
	 for n=0:(col-1)-dy
	    s(dx+row,dy+col,x(m+1,n+1)+1,x(m+1+dx,n+1+dy)+1)=s(dx+row,dy+col,x(m+1,n+1)+1,x(m+1+dx,n+1+dy)+1)+1;
	 end;
      end;	  
   end;
end;

for dx=0:(row-1)
   for dy=-1:-1:-(col-1)
      for m=0:(row-1)-dx
	 for n=(col-1):-1:abs(dy)
	    s(dx+row,dy+col,x(m+1,n+1)+1,x(m+1+dx,n+1+dy)+1)=s(dx+row,dy+col,x(m+1,n+1)+1,x(m+1+dx,n+1+dy)+1)+1;
	 end;
      end;	  
   end;
end;

for dx=-1:-1:-(row-1)
   for dy=0:(col-1)
      for m=(row-1):-1:abs(dx)
	 for n=0:(col-1)-dy
	    s(dx+row,dy+col,x(m+1,n+1)+1,x(m+1+dx,n+1+dy)+1)=s(dx+row,dy+col,x(m+1,n+1)+1,x(m+1+dx,n+1+dy)+1)+1;
	 end;
      end;	  
   end;
end;

for dx=-1:-1:-(row-1)
   for dy=-1:-1:-(col-1)
      for m=(row-1):-1:abs(dx)
	 for n=(col-1):-1:abs(dy)
	    s(dx+row,dy+col,x(m+1,n+1)+1,x(m+1+dx,n+1+dy)+1)=s(dx+row,dy+col,x(m+1,n+1)+1,x(m+1+dx,n+1+dy)+1)+1;
	 end;
      end;	  
   end;
end;