myncc fixes for zero variance

Improved handling of zero variance in regions in templatematch search
window for myncc.

templatematch.m

@@ -184,7 +184,6 @@
         mix=mix([2 1])./super;
         dxy(ii,:)=mix+dxyo(ii,1:2);
         Cout(ii,:)=[Cmax mean(abs(C(:)))];
-
     end
 
     if ~isempty(showprogress)
@@ -261,14 +260,17 @@
 sT=size(T); sB=size(B);
 sz=sB+sT-1;
 meanT=sum(T(:))/numel(T);
-sigmaT=sqrt(sum((T(:)-meanT).^2));
+sigmaT=realsqrt(sum((T(:)-meanT).^2));
 if sigmaT~=0
     fT=fft2(rot90(T,2),sz(1),sz(2));
     fB=fft2(B,sz(1),sz(2));
     C=real(ifft2(fB.*fT));
     lsumB=localsum(B,sT);
-    sigmaB=sqrt(max(localsum(B.*B,sT)-(lsumB.^2)/numel(T),0));
-    C=(C-lsumB*meanT)./(sigmaT*max(sigmaB,sigmaT/1e5)); %Uses 1e5 div0 trick from D.Kroon (thanks)
+    %sigmaB=sqrt(max(localsum(B.*B,sT)-(lsumB.^2)/numel(T),0));
+    %C=(C-lsumB*meanT)./(sigmaT*max(sigmaB,sigmaT/1e5)); %not 100% robust. Uses 1e5 div0 trick from D.Kroon (thanks)
+    sigmaB=realsqrt(max(localsum(B.*B,sT)-(lsumB.^2)/numel(T),0)); %is the max really necessary ?
+    C=(C-lsumB*meanT)./(sigmaT*sigmaB); 
+    C(abs(C)>1.1)=0; %this can happen if sigmaB almost 0, but we still allow C<1.1 to accomodate potential rounding issue for perfect correlation.
 else
     C=zeros(sz);
 end
@@ -279,6 +281,7 @@
 xx=-wkeep(2):wkeep(2);
 yy=-wkeep(1):wkeep(1);
 
+
 function lsum=localsum(A,sz)
 %Fast local sum of A. Local being within the a footprint of size sz
 A = cumsum(padarray(A,sz),2);