function T=TriCorr2(x)
%
% TriCorr2 2-D Triple Correlation function
%
% Example s(m1,n1,m2,n2)=TriCorr2(x);
%

[M,N]=size(x);
T=zeros(M,N,M,N);
for dx1=0:M-1
   for dy1=0:N-1
      for dx2=dx1:M-1
	 for dy2=dy1:N-1
            for i=0:M-1-dx2
	       for j=0:N-1-dy2
	          T(dx1+1,dy1+1,dx2+1,dy2+1)=T(dx1+1,dy1+1,dx2+1,dy2+1)+x(i+1,j+1)*x(i+1+dx1,j+1+dy1)*x(i+1+dx2,j+1+dy2);
               end;
            end;
            T(dx1+1,dy1+1,dx2+1,dy2+1)=T(dx1+1,dy1+1,dx2+1,dy2+1)/(M-dx2)/(N-dy2);
	 end;
      end;	  
   end;
end;
 
for dx1=0:M-1
   for dy2=0:N-1
      for dx2=dx1:M-1
	 for dy1=dy2+1:N-1
            for i=0:M-1-dx2
	       for j=0:N-1-dy1
	          T(dx1+1,dy1+1,dx2+1,dy2+1)=T(dx1+1,dy1+1,dx2+1,dy2+1)+x(i+1,j+1)*x(i+1+dx1,j+1+dy1)*x(i+1+dx2,j+1+dy2);
               end;
            end;
            T(dx1+1,dy1+1,dx2+1,dy2+1)=T(dx1+1,dy1+1,dx2+1,dy2+1)/(M-dx2)/(N-dy1);
	 end;
      end;	  
   end;
end; 
 
for dx2=0:M-1
   for dy1=0:N-1
      for dx1=dx2+1:M-1
	 for dy2=dy1:N-1
            for i=0:M-1-dx1
	       for j=0:N-1-dy2
	          T(dx1+1,dy1+1,dx2+1,dy2+1)=T(dx1+1,dy1+1,dx2+1,dy2+1)+x(i+1,j+1)*x(i+1+dx1,j+1+dy1)*x(i+1+dx2,j+1+dy2);
               end;
            end;
            T(dx1+1,dy1+1,dx2+1,dy2+1)=T(dx1+1,dy1+1,dx2+1,dy2+1)/(M-dx1)/(N-dy2);
	 end;
      end;	  
   end;
end;
  
for dx2=0:M-1
   for dy2=0:N-1
      for dx1=dx2+1:M-1
	 for dy1=dy2+1:N-1
            for i=0:M-1-dx1
	       for j=0:N-1-dy1
	          T(dx1+1,dy1+1,dx2+1,dy2+1)=T(dx1+1,dy1+1,dx2+1,dy2+1)+x(i+1,j+1)*x(i+1+dx1,j+1+dy1)*x(i+1+dx2,j+1+dy2);
               end;
            end;
            T(dx1+1,dy1+1,dx2+1,dy2+1)=T(dx1+1,dy1+1,dx2+1,dy2+1)/(M-dx1)/(N-dy1);
	 end;
      end;	  
   end;
end;