Re: improved __contourc__.cc + filled contours

[Shai Ayal]

After installing guplot 4.2.2 and octave 2.9.15+ it looks to be
working. gnuplot still complains from time to time more so with the
x11 terminal and very little with the wxt or ps  terminal

All the examples work, and I suggest this goes into 3.0 -- contourf is
quite a handy function, and even if it still has some undiscovered
bugs, it is not a critical function and could be fixed later.

Shai

p.s. The output with the wxt terminal is gorgeous !


On 10/12/07, Kai Habel <address@hidden> wrote:
> 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
>
>
>