improved __contourc__.cc + filled contours

[Kai Habel]

Hello,

since octave has support for patch objects since some time now, I have 
ported contourf.m (filled contours) from octplot to octave itself.

As a first step a modification for __contourc__.cc is required. For some 
inputs contourc.m (__contourc__.cc) "loses" the contour, that means e.g. 
a closed contour consists of two or more segments. For the contour.m
function this is no problem, but if you draw patches with this input you
get disordered patches.

See the following examples:

##this test fails (two contour segements)
test1 = [0 0 1;0 1 1;0 0 0]
contourc(test1,[0.5 0.5])
input('press enter');

##this test fails (two contour segements)
test2 = [0 0 0 0 0;0 1 1 1 0;0 1 0.5 0 0;0 1 1 1 0;0 0 0 0 0]
contourc(test2,[0.5 0.5])
input('press enter');

##this test does not fail (only one contour segements)
contourc(test2,[0.6 0.6])
input('press enter');

##this test fails (two contour segements)
test3 = [0 0 0 0 0;0 0 1 1 0;0 1 0 1 0;0 1 1 1 0;0 0 0 0 0]
contourc(test3,[0.5 0.5])
input('press enter');

If you look at the matrices test[1-3] you would expect that a single (or 
two) closed contour should be returned from contourc. But as you can 
see, this is not the case. Therefore I propose the following heavily 
rewritten __contourc__.cc. (I hope it is o.k. to send the new file only, 
since a patch would be larger)

The new contourf function works only with this modified __contourc__.cc.
I would like to ask you to test both __contour__.cc and contourf.m. I 
did a lot of tests, but maybe your input matrix causes one of the 
functions to fail. In addtion I have attached small test script, which 
shows some problems with the contour function (to be tackled later)

Since these functions need some review and we are close to a 3.0 release 
(I think) we should target a later inclusion.


Kai

[x,y,z]=peaks(75);
figure(1),clf
subplot(2,1,1)
contour(z);
axis([0 80 0 80])
subplot(2,1,2)
contourf3(z);
axis([0 80 0 80])
disp('contour(z)');
input('press enter');

figure(1),clf
subplot(2,1,1)
contour(x,y,z);
axis([-3 3 -3 3])
subplot(2,1,2)
contourf3(x,y,z);
axis([-3 3 -3 3])
disp('contour(x,y,z)');
input('press enter');

figure(1),clf
subplot(2,1,1)
contour(x,y,z,16);
axis([-3 3 -3 3])
subplot(2,1,2)
contourf3(x,y,z,16);
axis([-3 3 -3 3])
disp('contour(x,y,z,16)');
input('press enter');

figure(1),clf
subplot(2,1,1)
#contour(x,y,z,[0 0]);
axis([-3 3 -3 3])
subplot(2,1,2)
contourf3(x,y,z,[0 0]);
axis([-3 3 -3 3])
disp('contour(x,y,z,[0 0])');
input('press enter');

figure(1),clf
subplot(2,1,1)
#contour(x,y,z,[0 1],'EdgeColor','c');
axis([-3 3 -3 3])
subplot(2,1,2)
contourf3(x,y,z,[0 1],'EdgeColor','c');
axis([-3 3 -3 3])
disp("contour(x,y,z,\'EdgeColor\',\'c\');")
input('press enter');

figure(1),clf
subplot(2,1,1)
#contour(z,[0 1],'EdgeColor','c');
axis([0 80 0 80])
subplot(2,1,2)
contourf3(z,[0 1],'EdgeColor','c');
axis([0 80 0 80])
disp("contour(z,[0 1],\'EdgeColor\',\'c\');")
input('press enter');

figure(1),clf
subplot(2,1,1)
#contour(z,'EdgeColor','m');
axis([0 80 0 80])
subplot(2,1,2)
contourf3(z,'EdgeColor','m');
axis([0 80 0 80])
disp("contour(z,\'EdgeColor\',\'m\');")
input('press enter');

figure(1),clf
subplot(2,1,1)
#contour(z,16,'EdgeColor','k');
axis([0 80 0 80])
subplot(2,1,2);
contourf3(z,16,'EdgeColor','k');
axis([0 80 0 80])
disp("contour(z,16,\'EdgeColor\',\'k\');")
input('press enter');

## 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} {} address@hidden, @var{H}] = contourf 
(@var{x},@var{y},@var{z},@var{lvl})
## @deftypefnx {Function File} {} address@hidden, @var{H}] = contourf 
(@var{x},@var{y},@var{z},@var{n})
## @deftypefnx {Function File} {} address@hidden, @var{H}] = contourf 
(@var{x},@var{y},@var{z})
## @deftypefnx {Function File} {} address@hidden, @var{H}] = contourf 
(@var{z},@var{n})
## @deftypefnx {Function File} {} address@hidden, @var{H}] = contourf 
(@var{z},@var{lvl})
## @deftypefnx {Function File} {} address@hidden, @var{H}] = contourf (@var{z})
## @deftypefnx {Function File} {} address@hidden, @var{H}] = contourf 
(@var{ax},...)
## @deftypefnx {Function File} {} address@hidden, @var{H}] = contourf 
(...,@var{"property"},@var{val})
## This function computes 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 @var{ax} is provided the filled contours are added can be used to 
## @example
## address@hidden, @var{y}, @var{z}] = peaks(50);
## contourf (@var{x}, @var{y}, @var{z}, -7:9)
## @end example
## @end deftypefn
## @seealso{contour,contourc,patch}

## Authors: Kai Habel, shaia

function varargout = contourf(varargin)

  [X, Y, Z, lvl, ax, patch_props] = parse_args(varargin);
  [nr, nc] = size(Z);
  [minx, maxx] = deal(min(min(X)), max(max(X)));
  [miny, maxy] = deal(min(min(Y)), max(max(Y)));

  if (diff(lvl) < 10*eps) 
    lvl_eps = 1e-6;
  else
    lvl_eps = min(diff(lvl)) / 1000.0;
  end

  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:length(lev)
    lvl_idx = find(abs(cont_lev - lev(k)) < lvl_eps);
    len = length(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 += 1;
      end
      if (length(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;
        end
      end
    end
  end

  # the algoritm 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 (length(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) = [];
  end

  # now we have everything together
  # and can start plotting the patches
  # beginning with largest area
  [tmp, svec] = sort(cont_area);
  len = ncont - length(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{:});
  end

  if (min(lev) == max(lev))
    set(gca(), "clim", [min(lev) - 1 max(lev) + 1]);
  else
    set(gca(), "clim", [min(lev) max(lev)]);
  end

  if (nargout > 0)
    varargout{1} = c;
  end
  if (nargout > 1)
    varargout{2} = h;
  end
endfunction

function [X, Y, Z, lvl, ax, patch_props] = parse_args(arg)

  patch_props = {};
  nolvl = false;

  if (isinteger(arg{1}) && ishandle(arg{1}) && 
strncmp(get(arg{1},"Type"),"axis",4))
    ax = arg{1};
    arg{1} = [];
  else
    ax = gca();
  end

  for n = 1 : length(arg)
    if (ischar(arg{n}))
      patch_props = arg(n:end);
      arg(n:end) = [];
      break;
    endif
  endfor

  if (mod(length(patch_props), 2) != 0)
    error("property value is missing");
  endif

  if (length(arg) < 3)
    Z = arg{1};
    [X, Y] = meshgrid(1 : columns(Z), 1 : rows(Z));
  end

  if (length(arg) == 1)
    nolvl = true;
    arg(1) = [];
  elseif (length(arg) == 2)
    lvl = arg{2};
    arg(1:2) = [];
  elseif (length(arg) == 3)
    arg{1:3};
    [X, Y, Z] = deal(arg{1:3});
    arg(1:3) = [];
    nolvl = true;
  elseif (length(arg) == 4)
    [X, Y, Z, lvl] = deal(arg{1:4});
    arg(1:4) = [];
  endif

  if (any(size(X) != size(Y)))
    usage("X and Y must be of same size")
  end

  if (isvector(X) || isvector(Y))
    [X, Y] = meshgrid(X, Y);
  end

  Z_no_inf = Z(!isinf(Z));
  [minz, maxz] = deal(min(min(Z_no_inf)), max(max(Z_no_inf)));
  Z(isnan(Z)) = -inf;

  if (nolvl)
    lvl = linspace(minz, maxz, 12);
  end

  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];
    end
  end
endfunction

/* Contour lines for function evaluated on a grid.

  Rewrite of some functions to avoid segmented contours
  Copyright (C) 2007 Kai Habel

  Copyright (C) 2004 Shai Ayal

  Adapted to an oct file from the stand alone contourl by Victro Munoz
  Copyright (C) 2004 Victor Munoz

  Based on contour plot routine (plcont.c) in PLPlot package
  http://plplot.org/

  Copyright (C) 1995, 2000, 2001 Maurice LeBrun
  Copyright (C) 2000, 2002 Joao Cardoso
  Copyright (C) 2000, 2001, 2002, 2004  Alan W. Irwin
  Copyright (C) 2004  Andrew Ross

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Library Public License as published
  by the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

  This program 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 Library General Public License for more details.

  You should have received a copy of the GNU Library General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

*/

#include "oct.h"
#include <cfloat>
static Matrix this_contour;
static Matrix contourc;
static int elem;

// 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

**********************************************************************/
void
add_point (double x, double y)
{
  if (elem % CONTOUR_QUANT == 0)
    this_contour = this_contour.append (Matrix (2, CONTOUR_QUANT, 0));

  this_contour (0, elem) = x;
  this_contour (1, elem) = y;
  elem++;
}

/**********************************************************************

 end_contour();

 Adds contents of current contour to contourc.
 this_contour.cols () - 1;
**********************************************************************/

void
end_contour ()
{
  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;
}

/**********************************************************************

  start_contour(lvl, x, y);

  Start a new contour, and adds contents of current one to contourc

**********************************************************************/

void
start_contour (double lvl, double x, double y)
{
  end_contour ();
  this_contour.resize (2, 0);
  add_point (lvl, 0);
  add_point (x, y);
}

void
drawcn (RowVector & X, RowVector & Y, Matrix & Z, double lvl,
           int r, int c, double ct_x, double ct_y, uint start_edge, bool first,
           charMatrix & mark)
{
  double px[4], py[4], pz[4], tmp;
  uint stop_edge, next_edge, pt[2];
  int next_r, next_c;

  //get x, y, and z - lvl for current facet
  px[0] = px[3] = X (c);
  px[1] = px[2] = X (c + 1);
  py[0] = py[1] = Y (r);
  py[2] = py[3] = Y (r + 1);
  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;

  // facet edge and point naming assignment
  //  0-----1   .-0-.
  //  |     |   |   |
  //  |     |   3   1
  //  |     |   |   |
  //  3-----2   .-2-.

  // get mark value of current facet
  char id = static_cast<char>(mark(r, c));

  //check startedge s
  if (start_edge == 255)
    {
      //find start edge
      for (uint k = 0; k < 4; k++) {
        if (static_cast<char>(pow(2, k)) & id)
          start_edge = k;
      }
    }

  if (start_edge == 255)
    return;

  //decrease mark value of current facet for start edge
  mark(r, c) -= static_cast<char>(pow(2, start_edge));

  // 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
      {
        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 (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()))
    if (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);
      }
}

void mark_facets(Matrix &Z, charMatrix &mark, double lvl)
{
  uint c, r, nr = mark.rows(), nc = mark.cols();
  double f[4];

  for (c = 0; c < nc; c++)
    for (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 (r = 0; r < nr; r++)
    for (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;
      }
}

void cl_cntr (RowVector & X, RowVector & Y, Matrix & Z, double lvl)
{
  uint r, c, nr = Z.rows(), nc = Z.cols();
  charMatrix mark (nr - 1, nc - 1, 0);
  mark_facets(Z, mark, lvl);
  // find contours that start at a domain edge
  for (c = 0; c < nc - 1; c++)
    {
      //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 (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 (r = 0; r < nr - 1; r++)
    for (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, nargout, "")
{
  octave_value_list retval;
  int nargin = args.length ();

  if (nargin != 4)
    {
      error ("must have 4 inputs (x,y,z,levels)");
      return retval;
    }

  RowVector X = args (0).row_vector_value ();
  RowVector Y = args (1).row_vector_value ();
  Matrix Z = args (2).matrix_value ();
  RowVector L = args (3).row_vector_value ();

  contourc.resize (2, 0);
  for (int i = 0; i < L.length (); i++)
    {
      cl_cntr (X, Y, Z, L (i));
    }
  end_contour ();

  retval (0) = contourc;
  return retval;
}