Merge pull request #17 from avinashlalotra/masterHEADmaster

Abinash's Work

1 files changed, 69 insertions, 33 deletions

diff --git a/macros/xcorr2.sci b/macros/xcorr2.sci
index 688dc3a..40d3b57 100644
--- a/macros/xcorr2.sci
+++ b/macros/xcorr2.sci
@@ -1,34 +1,70 @@
-function c = xcorr2 (a, b, biasflag)
-//
-//Calling Sequence
-//c = xcorr2 (a)
-//c = xcorr2 (a, b)
-//c = xcorr2 (a, b, biasflag)
-//Parameters 
-//a:
-//b:
-//biasflag:
-//Description
-//This is an Octave function.
-
-//Examples
-//xcorr2(5,0.8,'coeff')
-//ans =  1
-
-funcprot(0);
-
-rhs = argn(2)
-if(rhs<1 | rhs>3)
-error("Wrong number of input arguments.");
-end
-
-	select(rhs)
-	case 1 then
-	c = callOctave("xcorr2",a);
-	case 2 then
-	c = callOctave("xcorr2",a,b);
-	case 3 then
-	c = callOctave("xcorr2",a,b,biasflag);
-	end;
+/*Calling Sequence
+        c = xcorr2 (a)
+        c = xcorr2 (a, b)
+        c = xcorr2 (a, b, scale)
+Description:
+        Compute the 2D cross-correlation of matrices a and b.
+        If b is not specified, computes autocorrelation of a, i.e., same as xcorr (a, a).
+        The optional argument scale, defines the type of scaling applied to the cross-correlation matrix. Possible values are:
+        "none" (default)
+            No scaling.
+        "biased"
+            Scales the raw cross-correlation by the maximum number of elements of a and b involved in the generation of any element of c.
+        "unbiased"
+            Scales the raw correlation by dividing each element in the cross-correlation matrix by the number of products a and b used to generate that element.
+        "coeff"
+            Scales the normalized cross-correlation on the range of [0 1] so that a value of 1 corresponds to a correlation coefficient of 1.
+  Examples
+        xcorr2(5,0.8,'coeff')
+            ans =  1 */
+function c = xcorr2 (a, b, scale)
+  funcprot(0);
+  nargin=argn(2);
+  if nargin < 3 then
+    scale = "none"
+  end
+  if (nargin < 1 || nargin > 3)
+    error("Wrong number of inputs")
+  elseif (nargin == 2 && type (b) == 10 )
+    scale = b;
+    b        = a;
+  elseif (nargin == 1)
+    // we have to set this case here instead of the function line, because if
+    // someone calls the function with zero argument, since a is never set, we
+    // will fail with "`a' undefined" error rather that print_usage
+    b = a;
+  end
+  if (ndims (a) ~= 2 || ndims (b) ~= 2)
+    error ("xcorr2: input matrices must must have only 2 dimensions");
+  end
+  // compute correlation
+  [ma,na] = size(a);
+  [mb,nb] = size(b);
+  c = conv2 (a, conj (b (mb:-1:1, nb:-1:1)));
+  // bias routines by Dave Cogdell (cogdelld@asme.org)
+  // optimized by Paul Kienzle (pkienzle@users.sf.net)
+  // coeff routine by Carnë Draug (carandraug+dev@gmail.com)
+  switch  (scale)
+    case {"none"}
+      // do nothing, it's all done
+    case {"biased"}
+      c = c / ( min ([ma, mb]) * min ([na, nb]) );
+    case {"unbiased"}
+      lo   = min ([na,nb]);
+      hi   = max ([na, nb]);
+      row  = [ 1:(lo-1), lo*ones(1,hi-lo+1), (lo-1):-1:1 ];
+      lo   = min ([ma,mb]);
+      hi   = max ([ma, mb]);
+      col  = [ 1:(lo-1), lo*ones(1,hi-lo+1), (lo-1):-1:1 ]';
+      bias = col*row;
+      c    = c./bias;
+    case {"coeff"}
+      a = double (a);
+      b = double (b);
+      a = conv2 (a.^2, ones (size (b,1) , size( b ,2)));
+      b = sum(b(:).* conj(b(:)));
+      c(:,:) = c(:,:) ./ sqrt (a(:,:) * b);
+    else
+      error ("xcorr2: invalid type of scale %s", scale);
+  end
 endfunction
-