Re: improved __contourc__.cc + filled contours

[Shai Ayal]

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,

I tried your functions but all I get in the  lower subplot is a black
rectangle -- i.e. contourf does nto work for me. Presumably this is
beacuse I dont have the latest gnuplot?

I have gnuplot 4.0, octave 2.9.15

Shai

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