function c=correlation2D(x,y)
% c=correlation(x,y)

[n1,n2,nt]=size(x);

xm=mean(x,3);
ym=mean(y,3);

xa=x-xm(:,:,ones(1,nt));
ya=y-ym(:,:,ones(1,nt));

yv=mean(ya.^2,3);
xv=mean(xa.^2,3);

yok=abs(yv)>0;
xok=abs(xv)>0;
ok=xok&yok;

top=mean(xa.*ya,3);
bottom=sqrt(xv.*yv);
c=zeros(n1,n2);
c(ok)=top(ok)./bottom(ok);
c(~ok)=NaN;
