contourc patch

[Kai Habel]

Hello,

I propose to include the patch from Alexander and Peter [1]. I have
modified the original patch slightly and did some editorial changes.

With this patch we can plot something like:

[th, r] = meshgrid (linspace (0, 2*pi, 64), 0:.05:0.9);
[X, Y] = pol2cart (th, r);
f = ((X + i*Y).^4 - 1).^(1/4);
contour(X, Y, abs(f), 16)

I have tried this patch with a test script and have seen no regressions
for the normal case.

Kai

ChangeLog:

2007-12-19  Alexander Barth
  * scripts/contourc.m: Allow usage of irregular x,y data

[1]
http://www.cae.wisc.edu/pipermail/octave-maintainers/2007-October/004282.html

--- contourc.m.orig     2007-12-11 12:07:21.000000000 +0100
+++ contourc.m  2007-12-19 21:42:35.000000000 +0100
@@ -56,19 +56,20 @@
 ## @end deftypefn
 
 ## Author: shaia
+## Grid mapping contributed by Peter A. Gustafson (2007) and Alexander Barth 
 
 function [c, lev] = contourc (varargin)
 
   if (nargin == 1)
     vn = 10;
     z = varargin{1};
-    x = 1:size(z,2);
-    y = 1:size(z,1);
+    x = 1 : size(z, 2);
+    y = 1 : size(z, 1);
   elseif (nargin == 2)
     vn = varargin{2};
     z = varargin{1};
-    x = 1:size(z,2);
-    y = 1:size(z,1);
+    x = 1 : size(z, 2);
+    y = 1 : size(z, 1);
   elseif (nargin == 3)
     vn = 10;
     x = varargin{1};
@@ -89,25 +90,42 @@
     vv = unique (sort (vn));
   endif
 
-  ## Vectorize the x,y vectors, assuming they are output from meshgrid.
-  if (! isvector (x))
-    x = x(1,:);
-  endif
-
-  if (! isvector (y))
-    y = y(:,1);
-  endif
-
-  ## Make everyone the right dimensions.
-  if (size (x, 2) == 1)
-    x = x';
-  endif
-  if (size (y, 2) == 1)
-    y = y';
+  if (isvector (x) && isvector (y))
+    c = __contourc__ (x(:)', y(:)', z, vv);
+  else
+    ## Indexes x,y for the purpose of __contourc__.    
+    ii = 1:size (z,2);
+    jj = 1:size (z,1);
+  
+    ## Now call __contourc__ for the real work...
+    c = __contourc__ (ii, jj, z, vv);
+  
+    ## Map the contour lines from index space (i,j) back 
+    ## to the original grid (x,y)
+    i = 1;
+
+    while (i < size (c,2))
+      clen = c(2, i);      
+      ind = i + [1 : clen];
+
+      ci = c(1, ind);
+      cj = c(2,ind);
+
+      ## due to rounding errors some elements of ci and cj
+      ## can fall out of the range of ii and jj and interp2 would
+      ## return NA for those values. 
+      ## The permitted range is enforced here: 
+        
+      ci = max (ci, 1); ci = min (ci, size (z, 2));
+      cj = max (cj, 1); cj = min (cj, size (z, 1));
+        
+      c(1, ind) = interp2 (ii, jj, x, ci, cj);
+      c(2, ind) = interp2 (ii, jj, y, ci, cj);
+      
+      i = i + clen + 1;
+    endwhile
   endif
-
-  ## Now call __contourc__ for the real work...
-  c = __contourc__ (x, y, z, vv);
+    
   if (nargout == 2)
     lev = vv;
   endif