--- a/scripts/ChangeLog	Fri Oct 19 19:05:31 2007 +0000
+++ b/scripts/ChangeLog	Fri Oct 19 20:43:33 2007 +0000
@@ -1,3 +1,8 @@
+2007-10-19  Kai Habel  <kai.habel@gmx.de>
+
+	* plot/contourf.m: New function.
+	* plot/Makefile.in (SOURCES): Add it to the list.
+
 2007-10-19  John W. Eaton  <jwe@octave.org>
 
 	* plot/subplot.m: Doc fix.

--- a/scripts/plot/Makefile.in	Fri Oct 19 19:05:31 2007 +0000
+++ b/scripts/plot/Makefile.in	Fri Oct 19 20:43:33 2007 +0000
@@ -70,6 +70,7 @@
   closereq.m \
   contour.m \
   contourc.m \
+  contourf.m \
   drawnow.m \
   errorbar.m \
   figure.m \

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/contourf.m	Fri Oct 19 20:43:33 2007 +0000
@@ -0,0 +1,259 @@
+## Copyright (C) 2007 Kai Habel
+## Copyright (C) 2003 Shai Ayal
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## OctPlot is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with OctPlot; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+## 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{c}, @var{h}] =} contourf (@var{x}, @var{y}, @var{z}, @var{lvl})
+## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{x}, @var{y}, @var{z}, @var{n})
+## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{x}, @var{y}, @var{z})
+## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{z}, @var{n})
+## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{z}, @var{lvl})
+## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{z})
+## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{ax}, @dots{})
+## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@dots{}, @var{"property"}, @var{val})
+## Compute and plot filled contours of the matrix @var{z}.
+## Parameters @var{x}, @var{y} and @var{n} or @var{lvl} are optional.
+##
+## The return value @var{c} is a 2xn matrix containing the contour lines
+## as described in the help to the contourc function.
+##
+## The return value @var{h} is handle-vector to the patch objects creating
+## the filled contours.
+##
+## If @var{x} and @var{y} are ommited they are taken as the row/column
+## index of @var{z}.  @var{n} is a scalar denoting the number of lines
+## to compute.  Alternatively @var{lvl} is a vector containing the
+## contour levels. If only one value (e.g. lvl0) is wanted, set
+## @var{lvl} to [lvl0, lvl0].  If both @var{n} or @var{lvl} are omitted
+## a default value of 10 contour level is assumed.
+##
+## If provided, the filled contours are added to the axes object
+## @var{ax} instead of the current axis.
+##
+## The following example plots filled contours of the @code{peaks}
+## function.
+## @example
+## [x, y, z] = peaks (50);
+## contourf (x, y, z, -7:9)
+## @end example
+## @seealso{contour, contourc, patch}
+## @end deftypefn
+
+## Author: Kai Habel <kai.habel@gmx.de>
+## Author: shaia
+
+function varargout = contourf (varargin)
+
+  [X, Y, Z, lvl, ax, patch_props] = parse_args (varargin);
+  [nr, nc] = size (Z);
+  [minx, maxx] = deal (min (X(:)), max (X(:)));
+  [miny, maxy] = deal (min (Y(:)), max (Y(:)));
+
+  if (diff (lvl) < 10*eps) 
+    lvl_eps = 1e-6;
+  else
+    lvl_eps = min (diff (lvl)) / 1000.0;
+  endif
+
+  X0 = prepad(X, nc+1, 2 * X(1, 1) - X(1, 2), 2);
+  X0 = postpad(X0, nc+2, 2 * X(1, nc) - X(1, nc - 1), 2);
+  X0 = [X0(1, :); X0; X0(1, :)];
+  Y0 = prepad(Y, nr+1, 2 * Y(1, 1) - Y(2, 1), 1);
+  Y0 = postpad(Y0, nr+2, 2 * Y(nr, 1) - Y(nr - 1, 1));
+  Y0 = [Y0(:, 1), Y0, Y0(:, 1)];
+
+  Z0 = -Inf(nr+2, nc+2);
+  Z0(2:nr+1, 2:nc+1) = Z;
+  [c, lev] = contourc (X0, Y0, Z0, lvl);
+  cmap = colormap ();
+
+  levx = linspace (min (lev), max (lev), size (cmap, 1));
+
+  newplot ();
+
+  ## Decode contourc output format.
+  i1 = 1;
+  ncont = 0;
+  while (i1 < columns (c))
+    ncont++;
+    cont_lev(ncont) = c(1, i1);
+    cont_len(ncont) = c(2, i1);
+    cont_idx(ncont) = i1+1;
+
+    ii = i1+1:i1+cont_len(ncont);
+    cur_cont = c(:, ii);
+    c(:, ii);
+    startidx = ii(1);
+    stopidx = ii(end);
+    cont_area(ncont) = polyarea (c(1, ii), c(2, ii));
+    i1 += c(2, i1) + 1;
+  endwhile
+
+  ## Handle for each level the case where we have (a) hole(s) in a patch.
+  ## Those are to be filled with the color of level below or with the
+  ## background colour.
+  for k = 1:numel (lev)
+    lvl_idx = find (abs (cont_lev - lev(k)) < lvl_eps);
+    len = numel (lvl_idx);
+    if (len > 1)
+      ## mark = logical(zeros(size(lvl_idx)));
+      mark = false (size (lvl_idx));
+      a = 1;
+      while (a < len)
+        # take 1st patch
+        b = a + 1;
+        pa_idx = lvl_idx(a);
+        # get pointer to contour start, and contour length
+        curr_ct_idx = cont_idx(pa_idx);
+        curr_ct_len = cont_len(pa_idx);
+        # get contour
+        curr_ct = c(:, curr_ct_idx:curr_ct_idx+curr_ct_len-1);
+        b_vec = (a+1):len;
+        next_ct_pt_vec = c(:, cont_idx(lvl_idx(b_vec)));
+        in = inpolygon (next_ct_pt_vec(1,:), next_ct_pt_vec(2,:),
+			curr_ct(1, :), curr_ct(2, :));
+        mark(b_vec(in)) = !mark(b_vec(in));
+        a++;
+      endwhile
+      if (numel (mark) > 0)
+	## All marked contours describe a hole in a larger contour of
+	## the same level and must be filled with colour of level below.
+        ma_idx = lvl_idx(mark);
+        if (k > 1)
+	  ## Find color of level below.
+          tmp = find(abs(cont_lev - lev(k - 1)) < lvl_eps);
+          lvl_bel_idx = tmp(1);
+	  ## Set color of patches found.
+	  cont_lev(ma_idx) = cont_lev(lvl_bel_idx);
+        else
+	  ## Set lowest level contour to NaN.
+	  cont_lev(ma_idx) = NaN;
+        endif
+      endif
+    endif
+  endfor
+
+  ## The algorithm can create patches with the size of the plotting
+  ## area, we would like to draw only the patch with the highest level.
+  del_idx = [];
+  max_idx = find (cont_area == max (cont_area));
+  if (numel (max_idx) > 1)
+    # delete double entries
+    del_idx = max_idx(1:end-1);
+    cont_area(del_idx) = cont_lev(del_idx) = [];
+    cont_len(del_idx) = cont_idx(del_idx) = [];
+  endif
+
+  ## Now we have everything together and can start plotting the patches
+  ## beginning with largest area.
+  [tmp, svec] = sort (cont_area);
+  len = ncont - numel (del_idx);
+  h = zeros (1, len);
+  for n = len:-1:1
+    idx = svec(n);
+    ii = cont_idx(idx):cont_idx(idx) + cont_len(idx) - 2;
+    h(n) = patch (c(1, ii), c(2, ii), cont_lev(idx), patch_props{:});
+  endfor
+
+  if (min (lev) == max (lev))
+    set (gca (), "clim", [min(lev)-1, max(lev)+1]);
+  else
+    set (gca(), "clim", [min(lev), max(lev)]);
+  endif
+
+  if (nargout > 0)
+    varargout{2} = h;
+    varargout{1} = c;
+  endif
+
+endfunction
+
+function [X, Y, Z, lvl, ax, patch_props] = parse_args (arg)
+
+  patch_props = {};
+  nolvl = false;
+
+  if (isinteger (arg{1}) && ishandle (arg{1})
+      && strncmpi (get (arg{1}, "type"), "axis", 4))
+    ax = arg{1};
+    arg{1} = [];
+  else
+    ax = gca ();
+  endif
+
+  for n = 1:numel (arg)
+    if (ischar (arg{n}))
+      patch_props = arg(n:end);
+      arg(n:end) = [];
+      break;
+    endif
+  endfor
+
+  if (mod (numel (patch_props), 2) != 0)
+    error ("patch: property value is missing");
+  endif
+
+  if (numel (arg) < 3)
+    Z = arg{1};
+    [X, Y] = meshgrid (1:columns (Z), 1:rows (Z));
+  endif
+
+  if (numel (arg) == 1)
+    nolvl = true;
+    arg(1) = [];
+  elseif (numel (arg) == 2)
+    lvl = arg{2};
+    arg(1:2) = [];
+  elseif (numel (arg) == 3)
+    arg{1:3};
+    [X, Y, Z] = deal (arg{1:3});
+    arg(1:3) = [];
+    nolvl = true;
+  elseif (numel (arg) == 4)
+    [X, Y, Z, lvl] = deal (arg{1:4});
+    arg(1:4) = [];
+  endif
+
+  if (any (size (X) != size (Y)))
+    error ("patch: X and Y must be of same size")
+  endif
+
+  if (isvector (X) || isvector (Y))
+    [X, Y] = meshgrid (X, Y);
+  endif
+
+  Z_no_inf = Z(!isinf (Z));
+  [minz, maxz] = deal (min (Z_no_inf(:)), max (Z_no_inf(:)));
+  Z(isnan (Z)) = -Inf;
+
+  if (nolvl)
+    lvl = linspace (minz, maxz, 12);
+  endif
+
+  if (isscalar (lvl))
+    lvl = linspace (minz, maxz, lvl + 2)(1:end-1);
+  else
+    idx1 = find(lvl < minz);
+    idx2 = find(lvl > maxz);
+    lvl(idx1(1:end-1)) = [];
+    lvl(idx2) = [];
+    if (isempty (lvl))
+      lvl = [minz, minz];
+    endif
+  endif
+
+endfunction

--- a/src/ChangeLog	Fri Oct 19 19:05:31 2007 +0000
+++ b/src/ChangeLog	Fri Oct 19 20:43:33 2007 +0000
@@ -1,3 +1,14 @@
+2007-10-19  Kai Habel  <kai.habel@gmx.de>
+
+	* DLD-FUNCTIONS/__contourc__.cc (add_point): Rename from
+	cl_add_point.  Change all uses.
+	(end_contour): Rename from cl_end_contour.  Change all uses.
+	(start_contour): Rename from cl_start_contour.  Change all uses.
+	(drawcn): Rename from cl_drawcn.  New algorithm for locating contours.
+	(mark_facets): New function.
+	(cntr): Rename from cl_cntr.  Change all uses.  New algorithm for
+	locating contours.
+
 2007-10-19  John W. Eaton  <jwe@octave.org>
 
 	* pt-idx.cc (tree_index_expression::lvalue): Correctly compute

--- a/src/DLD-FUNCTIONS/__contourc__.cc	Fri Oct 19 19:05:31 2007 +0000
+++ b/src/DLD-FUNCTIONS/__contourc__.cc	Fri Oct 19 20:43:33 2007 +0000
@@ -1,5 +1,6 @@
 /* Contour lines for function evaluated on a grid.
 
+Copyright (C) 2007 Kai Habel
 Copyright (C) 2004, 2007 Shai Ayal
 
 Adapted to an oct file from the stand alone contourl by Victro Munoz
@@ -35,6 +36,8 @@
 #include <config.h>
 #endif
 
+#include <cfloat>
+
 #include "quit.h"
 
 #include "defun-dld.h"
@@ -45,15 +48,13 @@
 static Matrix contourc;
 static int elem;
 
-// this is the quanta in which we increase this_contour
+// This is the quanta in which we increase this_contour.
 #define CONTOUR_QUANT 50
 
-// cl_add_point(x,y);
-//
-// Add a coordinate point (x,y) to this_contour 
+// Add a coordinate point (x,y) to this_contour.
 
 static void
-cl_add_point (double x, double y)
+add_point (double x, double y)
 {
   if (elem % CONTOUR_QUANT == 0)
     this_contour = this_contour.append (Matrix (2, CONTOUR_QUANT, 0));
@@ -63,237 +64,238 @@
   elem++;
 }
 
-// cl_end_contour();
-//
-// Adds contents of current contour to contourc.
+// Add contents of current contour to contourc.
 // this_contour.cols () - 1;
 
 static void
-cl_end_contour (void)
+end_contour (void)
 {
   if (elem > 2)
     {
       this_contour (1, 0) = elem - 1;
       contourc = contourc.append (this_contour.extract_n (0, 0, 2, elem));
     }
+
   this_contour = Matrix ();
   elem = 0;
 }
 
-// cl_start_contour(flev,x,y);
-//
-// Start a new contour, and adds contents of current one to contourc
+// Start a new contour, and add contents of current one to contourc.
 
 static void
-cl_start_contour (double flev, double x, double y)
+start_contour (double lvl, double x, double y)
 {
-  cl_end_contour ();
+  end_contour ();
   this_contour.resize (2, 0);
-  cl_add_point (flev, flev);
-  cl_add_point (x, y);
+  add_point (lvl, 0);
+  add_point (x, y);
 }
 
 static void
-cl_drawcn (RowVector & X, RowVector & Y, Matrix & Z, double flev,
-	   int krow, int kcol, double lastx, double lasty, int startedge,
-	   Matrix & ipts)
+drawcn (const RowVector& X, const RowVector& Y, const Matrix& Z,
+	double lvl, int r, int c, double ct_x, double ct_y,
+	uint start_edge, bool first, charMatrix& mark)
 {
-
-  int kx = 0, lx = Z.cols () - 1, ky = 0, ly = Z.rows () - 1;
-
-  double f[4];
-  double px[4], py[4], locx[4], locy[4];
-  int iedge[4];
-  int num, first, inext, kcolnext, krownext;
+  double px[4], py[4], pz[4], tmp;
+  uint stop_edge, next_edge, pt[2];
+  int next_r, next_c;
 
-  px[0] = X (krow + 1);
-  px[1] = X (krow);
-  px[2] = X (krow);
-  px[3] = X (krow + 1);
-  py[0] = Y (kcol);
-  py[1] = Y (kcol);
-  py[2] = Y (kcol + 1);
-  py[3] = Y (kcol + 1);
+  //get x, y, and z - lvl for current facet
+  px[0] = px[3] = X(c);
+  px[1] = px[2] = X(c+1);
 
-  f[0] = Z (krow + 1, kcol) - flev;
-  f[1] = Z (krow, kcol) - flev;
-  f[2] = Z (krow, kcol + 1) - flev;
-  f[3] = Z (krow + 1, kcol + 1) - flev;
+  py[0] = py[1] = Y(r);
+  py[2] = py[3] = Y(r+1);
 
-  for (int i = 0, j = 1; i < 4; i++, j = (j + 1) % 4)
-    {
-      iedge[i] = (f[i] * f[j] > 0.0) ? -1 : ((f[i] * f[j] < 0.0) ? 1 : 0);
-    }
-
-  // Mark this square as done
-  ipts(krow,kcol) = 1;
-
-  // Check if no contour has been crossed i.e. iedge[i] = -1
-  if (iedge[0] == -1 && iedge[1] == -1 && iedge[2] == -1 && iedge[3] == -1)
-    return;
+  pz[3] = Z(r+1, c) - lvl;
+  pz[2] = Z(r+1, c + 1) - lvl;
+  pz[1] = Z(r, c+1) - lvl;
+  pz[0] = Z(r, c) - lvl;
 
-  // Check if this is a completely flat square - in which case ignore it
-  if (f[0] == 0.0 && f[1] == 0.0 && f[2] == 0.0 && f[3] == 0.0)
-    return;
+  // Facet edge and point naming assignment.
+  //
+  //  0-----1   .-0-.
+  //  |     |   |   |
+  //  |     |   3   1
+  //  |     |   |   |
+  //  3-----2   .-2-.
 
-  // Calculate intersection points
-  num = 0;
-  if (startedge < 0)
+  // Get mark value of current facet.
+  char id = static_cast<char> (mark(r, c));
+
+  // Check startedge s.
+  if (start_edge == 255)
     {
-      first = 1;
-    }
-  else
-    {
-      locx[num] = lastx;
-      locy[num] = lasty;
-      num++;
-      first = 0;
+      // Find start edge.
+      for (uint k = 0; k < 4; k++)
+        if (static_cast<char> (pow(2, k)) & id)
+          start_edge = k;
     }
 
-  for (int k = 0, i = (startedge < 0 ? 0 : startedge); k < 4;
-       k++, i = (i + 1) % 4)
-    {
-      if (i == startedge)
-	continue;
+  if (start_edge == 255)
+    return;
+
+  // Decrease mark value of current facet for start edge.
+  mark(r, c) -= static_cast<char> (pow(2, start_edge));
 
-      // If the contour is an edge check it hasn't already been done
-      if (f[i] == 0.0 && f[(i + 1) % 4] == 0.0)
-	{
-	  kcolnext = kcol;
-	  krownext = krow;
-	  if (i == 0)
-	    kcolnext--;
-	  if (i == 1)
-	    krownext--;
-	  if (i == 2)
-	    kcolnext++;
-	  if (i == 3)
-	    krownext++;
-	  if (kcolnext < kx || kcolnext >= lx || krownext < ky
-	      || krownext >= ly || ipts(krownext,kcolnext) == 1)
-	    continue;
-	}
-      if (iedge[i] == 1 || f[i] == 0.0)
+  // Next point (clockwise).
+  pt[0] = start_edge;
+  pt[1] = (pt[0] + 1) % 4;
+
+  // Calculate contour segment start if first of contour.
+  if (first)
+    {
+      tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
+
+      if (xisnan (tmp))
+        ct_x = ct_y = 0.5;
+      else
 	{
-	  int j = (i + 1) % 4;
-	  if (f[i] != 0.0)
-	    {
-	      locx[num] =
-		(px[i] * fabs (f[j]) + px[j] * fabs (f[i])) / fabs (f[j] -
-								    f[i]);
-	      locy[num] =
-		(py[i] * fabs (f[j]) + py[j] * fabs (f[i])) / fabs (f[j] -
-								    f[i]);
-	    }
-	  else
-	    {
-	      locx[num] = px[i];
-	      locy[num] = py[i];
-	    }
-	  // If this is the start of the contour then move to the point
-	  if (first == 1)
-	    {
-	      cl_start_contour (flev, locx[num], locy[num]);
-	      first = 0;
-	    }
-	  else
-	    {
-	      // Link to the next point on the contour
-	      cl_add_point (locx[num], locy[num]);
-	      // Need to follow contour into next grid box
-	      // Easy case where contour does not pass through corner
-	      if (f[i] != 0.0)
-		{
-		  kcolnext = kcol;
-		  krownext = krow;
-		  inext = (i + 2) % 4;
-		  if (i == 0)
-		    kcolnext--;
-		  if (i == 1)
-		    krownext--;
-		  if (i == 2)
-		    kcolnext++;
-		  if (i == 3)
-		    krownext++;
-		  if (kcolnext >= kx && kcolnext < lx
-		      && krownext >= ky && krownext < ly
-		      && ipts(krownext,kcolnext) == 0)
-		    {
-		      cl_drawcn (X, Y, Z, flev, krownext, kcolnext,
-				 locx[num], locy[num], inext, ipts);
-		    }
-		}
-	      // Hard case where contour passes through corner.  This
-	      // is still not perfect - it may lose the contour  which
-	      // won't upset the contour itself (we can find it again
-	      // later) but might upset the labelling (which is only
-	      // relevant for the PLPlot implementation, since we
-	      // don't worry about labels---for now!)
-	      else
-		{
-		  kcolnext = kcol;
-		  krownext = krow;
-		  inext = (i + 2) % 4;
-		  if (i == 0)
-		    {
-		      kcolnext--;
-		      krownext++;
-		    }
-		  if (i == 1)
-		    {
-		      krownext--;
-		      kcolnext--;
-		    }
-		  if (i == 2)
-		    {
-		      kcolnext++;
-		      krownext--;
-		    }
-		  if (i == 3)
-		    {
-		      krownext++;
-		      kcolnext++;
-		    }
-		  if (kcolnext >= kx && kcolnext < lx
-		      && krownext >= ky && krownext < ly
-		      && ipts(krownext,kcolnext) == 0)
-		    {
-		      cl_drawcn (X, Y, Z, flev, krownext, kcolnext,
-				 locx[num], locy[num], inext, ipts);
-		    }
+	  ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
+	  ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
+	}
+
+      start_contour (lvl, ct_x, ct_y);
+    }
+
+  // Find stop edge FIXME: control flow --> while.
+  for (uint k = 1; k <= 4; k++)
+    {
+      if (start_edge == 0 || start_edge == 2)
+        stop_edge = (start_edge + k) % 4;
+      else
+        stop_edge = (start_edge - k) % 4;
+
+      if (static_cast<char> (pow(2, stop_edge)) & id)
+        break;
+    }
+
+  pt[0] = stop_edge;
+  pt[1] = (pt[0] + 1) % 4;
+  tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
 
-		}
-	      if (first == 1)
-		{
-		  // Move back to first point
-		  cl_start_contour (flev, locx[num], locy[num]);
-		  first = 0;
-		}
-	      else
-		{
-		  first = 1;
-		}
-	      num++;
-	    }
-	}
+  if (xisnan (tmp))
+    ct_x = ct_y = 0.5;
+  else
+    {
+      ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
+      ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
+    }
+
+  // Add point to contour.
+  add_point (ct_x, ct_y);
+
+  // Decrease id value of current facet for start edge.
+  mark(r, c) -= static_cast<char>(pow(2,stop_edge));
+
+  // Find next facet.
+  next_c = c;
+  next_r = r;
+
+  if (stop_edge == 0)
+    next_r--; 
+  else if (stop_edge == 1)
+    next_c++;
+  else if (stop_edge == 2)
+    next_r++;
+  else if (stop_edge == 3)
+    next_c--;
+
+  // Check if next facet is not done yet.
+  // Go to next facet.
+  if (next_r >= 0 && next_c >= 0 && next_r < mark.rows ()
+      && next_c < mark.cols () && mark(next_r, next_c) > 0)
+    {
+      next_edge = (stop_edge + 2) % 4;
+      drawcn (X, Y, Z, lvl, next_r, next_c, ct_x, ct_y, next_edge, false, mark);
     }
 }
 
 static void
-cl_cntr (RowVector & X, RowVector & Y, Matrix & Z, double flev)
+mark_facets (const Matrix& Z, charMatrix& mark, double lvl)
 {
-  Matrix ipts (Z.rows (), Z.cols (), 0);
+  uint nr = mark.rows ();
+  uint nc = mark.cols ();
+
+  double f[4];
+
+  for (uint c = 0; c < nc; c++)
+    for (uint r = 0; r < nr; r++)
+      {
+        f[0] = Z(r, c) - lvl;
+        f[1] = Z(r, c+1) - lvl;
+        f[3] = Z(r+1, c) - lvl;
+        f[2] = Z(r+1, c+1) - lvl;
+
+        for (uint i = 0; i < 4; i++)
+          if (fabs(f[i]) < DBL_EPSILON)
+            f[i] = DBL_EPSILON;
+
+        if (f[1] * f[2] < 0)
+	  mark(r, c) += 2;
+
+        if (f[0] * f[3] < 0)
+	  mark(r, c) += 8;
+      }
 
-  for (int krow = 0; krow < Z.rows () - 1; krow++)
+  for (uint r = 0; r < nr; r++)
+    for (uint c = 0; c < nc; c++)
+      {
+        f[0] = Z(r, c) - lvl;
+        f[1] = Z(r, c+1) - lvl;
+        f[3] = Z(r+1, c) - lvl;
+        f[2] = Z(r+1, c+1) - lvl;
+
+        for (uint i = 0; i < 4; i++)
+          if (fabs(f[i]) < DBL_EPSILON)
+            f[i] = DBL_EPSILON;
+
+        if (f[0] * f[1] < 0)
+	  mark(r, c) += 1;
+
+        if (f[2] * f[3] < 0)
+	  mark(r, c) += 4;
+      }
+}
+
+static void
+cntr (const RowVector& X, const RowVector& Y, const Matrix& Z, double lvl)
+{
+  uint nr = Z.rows ();
+  uint nc = Z.cols ();
+
+  charMatrix mark (nr - 1, nc - 1, 0);
+
+  mark_facets (Z, mark, lvl);
+
+  // Find contours that start at a domain edge.
+
+  for (uint c = 0; c < nc - 1; c++)
     {
-      for (int kcol = 0; kcol < Z.cols () - 1; kcol++)
-	{
-	  if (ipts(krow,kcol) == 0)
-	    {
-	      cl_drawcn (X, Y, Z, flev, krow, kcol, 0.0, 0.0, -2, ipts);
-	    }
-	}
+      // Top.
+      if (mark(0, c) & 1)
+        drawcn (X, Y, Z, lvl, 0, c, 0.0, 0.0, 0, true, mark);
+
+      // Bottom.
+      if (mark(nr - 2, c) & 4)
+	drawcn (X, Y, Z, lvl, nr - 2, c, 0.0, 0.0, 2, true, mark);
     }
+
+  for (uint r = 0; r < nr - 1; r++)
+    {
+      // Left.
+      if (mark(r, 0) & 8)
+        drawcn (X, Y, Z, lvl, r, 0, 0.0, 0.0, 3, true, mark);
+
+      // Right.
+      if (mark(r, nc - 2) & 2)
+        drawcn (X, Y, Z, lvl, r, nc - 2, 0.0, 0.0, 1, true, mark);
+    }
+
+  for (uint r = 0; r < nr - 1; r++)
+    for (uint c = 0; c < nc - 1; c++)
+      if (mark (r, c) > 0)
+        drawcn (X, Y, Z, lvl, r, c, 0.0, 0.0, 255, true, mark);
 }
 
 DEFUN_DLD (__contourc__, args, ,
@@ -308,7 +310,7 @@
     {
       RowVector X = args (0).row_vector_value ();
       RowVector Y = args (1).row_vector_value ();
-      Matrix Z = args (2).matrix_value ().transpose ();
+      Matrix Z = args (2).matrix_value ();
       RowVector L = args (3).row_vector_value ();
 
       if (! error_state)
@@ -316,9 +318,9 @@
 	  contourc.resize (2, 0);
 
 	  for (int i = 0; i < L.length (); i++)
-	    cl_cntr (X, Y, Z, L (i));
+	    cntr (X, Y, Z, L (i));
 
-	  cl_end_contour ();
+	  end_contour ();
 
 	  retval = contourc;
 	}