[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: contour.m by Shai Ayal
|
From: |
Shai Ayal |
|
Subject: |
Re: contour.m by Shai Ayal |
|
Date: |
Sun, 24 Jul 2005 11:06:18 +0300 |
|
User-agent: |
Mozilla Thunderbird 1.0 (X11/20041206) |
Henry,
please try the attached files. They include a newer version of the contour
algorithm, and some renaming.
put them all in a fresh directory and than
mkoctfile __contourc__.cc
and then try again.
They work for me in ver 2.1.71 in fedora with gnuplot 4
Shai
Henry F. Mollet wrote:
I'm trying to use Shai's contour.m (see attached 3 files). I've recompiled
contourl.cc to get the oct file to be on the save side and then tried
contour (x,y,z):
error: can't perform indexing operations for <unknown type> type. When I
changed to the directory with the Octave contour.m it worked but this does
not use an oct file.
What am I doing incorrectly? It worked for Shai in March 2005 in GNU
Octave, version 2.1.57 (i686-pc-linux-gnu).
Henry
[~/Documents/aNewMethodsPaperReview] -bash-2.05b 506$ mkdir Shai
[~/Documents/aNewMethodsPaperReview] -bash-2.05b 507$ cd Shai
[~/Documents/aNewMethodsPaperReview/Shai] -bash-2.05b 508$ ls
[~/Documents/aNewMethodsPaperReview/Shai] -bash-2.05b 509$ ls
contour.m contourc.m contourl.cc
[~/Documents/aNewMethodsPaperReview/Shai] -bash-2.05b 510$ mkoctfile
contourl.cc
/usr/bin/ld: warning multiple definitions of symbol _xerbla_
/usr/local/lib/octave-2.1.71/liboctinterp.dylib(single module) definition of
_xerbla_
/usr/local/lib/octave-2.1.71/libcruft.dylib(single module) definition of
_xerbla_
/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib
.framework/Versions/A/libBLAS.dylib(single module) definition of _xerbla_
/usr/bin/ld: warning multiple definitions of symbol _UP
/usr/local/lib/libreadline.dylib(terminal.so) definition of _UP
/usr/lib/gcc/powerpc-apple-darwin8/4.0.0/../../../libncurses.dylib(lib_termc
ap.o) definition of _UP
/usr/bin/ld: warning multiple definitions of symbol _BC
/usr/local/lib/libreadline.dylib(terminal.so) definition of _BC
/usr/lib/gcc/powerpc-apple-darwin8/4.0.0/../../../libncurses.dylib(lib_termc
ap.o) definition of _BC
/usr/bin/ld: warning multiple definitions of symbol _PC
/usr/local/lib/libreadline.dylib(terminal.so) definition of _PC
/usr/lib/gcc/powerpc-apple-darwin8/4.0.0/../../../libncurses.dylib(lib_tputs
.o) definition of _PC
[~/Documents/aNewMethodsPaperReview/Shai] -bash-2.05b 511$ ls
contour.m contourc.m contourl.cc contourl.o contourl.oct
[~/Documents/aNewMethodsPaperReview/Shai] -bash-2.05b 512$ octave
GNU Octave, version 2.1.71 (powerpc-apple-darwin8.1.0).
octave:1> pwd
/Users/hfm/Documents/aNewMethodsPaperReview/Shai
octave:2> x=linspace(0,1.2,35);
octave:3> y=linspace(1,35,35);
octave:4> [xx,yy]=meshgrid(x,y);
octave:5> z=xx./(yy.*(1-xx)) -
(3*yy-yy+1).*xx.^(3*yy-yy+1)./yy./(1-xx.^(3*yy-yy+1)) + eps;
octave:6> ls
contour.m contourc.m contourl.cc contourl.o contourl.oct
octave:7> contour(x,y,z)
error: can't perform indexing operations for <unknown type> type
octave:7> cd /usr/local/share/octave/2.1.71/m/plot
octave:8> ls
contour.m top_title.m
octave:9> contour (x,y,z) % This now works.
octave:10>
--
Shai Ayal, Ph.D.
Head of Research
BioControl Medical BCM
3 Geron St.
Yehud 56100
ISRAEL
Tel: + 972 3 6322 126 ext 223
Fax: + 972 3 6322 125
email: address@hidden
/* Contour lines for function evaluated on a grid.
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"
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
cl_add_point (float x, float 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++;
}
/**********************************************************************
cl_end_contour();
Adds contents of current contour to contourc.
this_contour.cols () - 1;
**********************************************************************/
void
cl_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;
}
/**********************************************************************
cl_start_contour(flev,x,y);
Start a new contour, and adds contents of current one to contourc
**********************************************************************/
void
cl_start_contour (float flev, float x, float y)
{
cl_end_contour ();
this_contour.resize (2, 0);
cl_add_point (flev, flev);
cl_add_point (x, y);
}
void
cl_drawcn (RowVector & X, RowVector & Y, Matrix & Z, float flev,
int krow, int kcol, float lastx, float lasty, int startedge,
Matrix & ipts)
{
int kx = 0, lx = Z.cols () - 1, ky = 0, ly = Z.rows () - 1;
float f[4];
float px[4], py[4], locx[4], locy[4];
int iedge[4];
int i, j, k, num, first, inext, kcolnext, krownext;
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);
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;
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;
/* 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;
/* Calculate intersection points */
num = 0;
if (startedge < 0)
{
first = 1;
}
else
{
locx[num] = lastx;
locy[num] = lasty;
num++;
first = 0;
}
for (k = 0, i = (startedge < 0 ? 0 : startedge); k < 4;
k++, i = (i + 1) % 4)
{
if (i == startedge)
continue;
/* 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))
{
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);
}
}
if (first == 1)
{
/* Move back to first point */
cl_start_contour (flev, locx[num], locy[num]);
first = 0;
}
else
{
first = 1;
}
num++;
}
}
}
}
void
cl_cntr (RowVector & X, RowVector & Y, Matrix & Z, float flev)
{
Matrix ipts (Z.rows (), Z.cols (), 0);
for (int krow = 0; krow < Z.rows () - 1; krow++)
{
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);
}
}
}
}
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 ().transpose ();
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));
}
cl_end_contour ();
retval (0) = contourc;
return retval;
}
## 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} = contour
(@var{x},@var{y},@var{z},@var{vv})
## This function computes isolines (countour lines) of the matrix @var{z}.
## parameters @var{x}, @var{y} and @var{vn} are optional.
##
## The return value @var{c} is a 2xn matrix containing the contour lines
## as dexcribed in the help to the contourc function
##
## if @var{x} and @var{y} are ommited they are taken as the row/column
## index of @var{z}. @var{vn} is either a scalar denoting the number of lines
## to compute or a vector contianing the values of the lines. If only one
## value is wanted, set @var{vn}=[val,val];
## If @var{vn} is omitted it defualts to 10
##
## @example
## contour (@var{x}, @var{y}, @var{z}, linspace(0,2*pi,10))
## @end example
## @end deftypefn
## @seealso{contourc,line,plot}
## Author: shaia,pkienzle
function ret = contour(varargin)
[c,lev] = contourc(varargin{:});
## find minimal distance between levels for coloring
dl = min(abs(diff(lev)))/10;
held = ishold;
unwind_protect
## decode contourc output format
i1 = 1;
while(i1<length(c))
ii = i1+1:i1+c(2,i1);
## finds the apropriate color for this level
line_color = mod(find(abs(lev-c(1,i1))<dl),6)+1;
plot(c(1,ii),c(2,ii),sprintf("-%d;;",line_color)); hold on;
i1 += c(2,i1)+1;
endwhile
unwind_protect_cleanup
if ~held,
hold off;
endif
end_unwind_protect
if nargout>0,
ret = c;
endif
endfunction
## 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} = contourc
(@var{x},@var{y},@var{z},@var{vv})
## This function computes isolines (countour lines) of the matrix @var{z}.
## parameters @var{x}, @var{y} and @var{vn} are optional.
##
## The return value @var{c} is a 2xn matrix containing the contour lines
## in the following format
##
## @example
## @var{c} = [value1 , x1 , x2 , ... , valuen , x1 , x2 , ...
## len1 , y1 , y2 , ... , lenn , y1 , y2 , ... ]
## @end example
##
## @noindent
## where contour line n has a value (height) of valun and length of
## lenn.
##
## if @var{x} and @var{y} are ommited they are taken as the row/column
## index of @var{z}. @var{vn} is either a scalar denoting the number of lines
## to compute or a vector contianing the values of the lines. If only one
## value is wanted, set @var{vn}=[val,val];
## If @var{vn} is omitted it defualts to 10
##
## @example
## @var{c}=contourc (@var{x}, @var{y}, @var{z}, linspace(0,2*pi,10))
## @end example
## @end deftypefn
## @seealso{contour}
## Author: shaia
function [c,vv] = contourc(varargin)
if (nargin==1)
vn = 10;
z = varargin{1};
x = 1:size(z,1);
y = 1:size(z,2);
elseif (nargin==2)
vn = varargin{2};
z = varargin{1};
x = 1:size(z,1);
y = 1:size(z,2);
elseif (nargin==3)
vn = 10;
x = varargin{1};
y = varargin{2};
z = varargin{3};
elseif (nargin==4)
vn = varargin{4};
x = varargin{1};
y = varargin{2};
z = varargin{3};
else
error("Wrong number of arguments");
endif
if (isscalar(vn))
vv=linspace(min(min(z)),max(max(z)),vn+2)(2:end-1);
else
vv = unique(sort(vn));
end
## vectorize the x,y vectors, assuming they are output from meshgrid
if ~isvector(x),
x = x(1,:);
endif
if ~isvector(y),
y = y(:,1);
endif
## make everyone the right dimensions
if(size(x,2)==1)
x = x';
endif
if(size(y,2)==1)
y = y';
endif
## now call __contourc__ for the real work ...
c=__contourc__(x,y,z,vv);
endfunction