[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
improved __contourc__.cc + filled contours
From: |
Kai Habel |
Subject: |
improved __contourc__.cc + filled contours |
Date: |
Thu, 11 Oct 2007 23:39:17 +0000 |
User-agent: |
IceDove 1.5.0.12 (X11/20070607) |
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;
}
- improved __contourc__.cc + filled contours,
Kai Habel <=
- Re: improved __contourc__.cc + filled contours, Shai Ayal, 2007/10/14
- Re: improved __contourc__.cc + filled contours, Shai Ayal, 2007/10/15
- wxt (was: Re: improved __contourc__.cc + filled contours), Søren Hauberg, 2007/10/15
- Re: wxt (was: Re: improved __contourc__.cc + filled contours), Shai Ayal, 2007/10/15
- Re: wxt (was: Re: improved __contourc__.cc + filled contours), Ben Abbott, 2007/10/15
- Re: wxt (was: Re: improved __contourc__.cc + filled contours), John W. Eaton, 2007/10/15
- Re: wxt (was: Re: improved __contourc__.cc + filled contours), Ben Abbott, 2007/10/15
- Re: wxt (was: Re: improved __contourc__.cc + filled contours), John W. Eaton, 2007/10/18
- Re: wxt (was: Re: improved __contourc__.cc + filled contours), John W. Eaton, 2007/10/15