[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: contour plot
From: |
Paul Kienzle |
Subject: |
Re: contour plot |
Date: |
Fri, 30 Aug 2002 12:16:10 -0400 |
Attached is an implementation of contourf based on this technique. It is
also available from octave-forge CVS.
## contourf(z,n,w)
## Plots a filled contour plot of n levels using the current
## colormap. The width w is the width of the convolution window
## which smooths the contours.
##
## E.g.,
## [x,y] = meshgrid(linspace(-2,2,200));
## z = sinc(sqrt(x.^2 + y.^2)) + 0.5*randn(size(x));
## filled_contour(z);
## This program is in the public domain
function contourf(z,n,w)
if nargin < 3, w = 16; end
if nargin < 2, n = 10; end
if nargin < 1 || nargin > 3
usage("contourf(z [, n [, w]])");
endif
## generate the gradient image from the original
M = imagesc(z);
## convolute the original with a gaussian if desired
if w > 0
[x,y] = meshgrid(2.5*linspace(-1,1,w));
B = exp(-.5*(x.^2+y.^2));
z = filter2(B,z);
endif
## find the contours
C = colormap;
colormap(rand(n+1,3));
z = filter2(ones(2)/4,imagesc(z));
M(z!=fix(z)) = 0;
## draw the gradient image with contours
colormap([0,0,0; C]);
image(flipud(M)+1);
## restore the colormap
colormap(C);
endfunction
Paul Kienzle
address@hidden
On Fri, Aug 30, 2002 at 10:50:59AM -0400, Paul Kienzle wrote:
> You can get a lot of the way there using imagesc.
>
> First, construct some data to contour:
>
> [x,y] = meshgrid(linspace(-2,2,200));
> z = sinc(sqrt(x.^2 + y.^2));
> colormap(rand(9,3));
> imagesc(z);
>
> That's unrealistic since real data is noisy. With
> the added noise, the contours aren't very pretty:
>
> Z = z + 0.05*randn(size(z));
> imagesc(Z);
>
> We can clean up the boundaries by convoluting with a gaussian:
>
> [x,y] = meshgrid(2.5*linspace(-1,1,16));
> B = exp(-.5*(x.^2+y.^2));
> imagesc(filter2(B,Z,"same"))
>
> It would be nice to be able to outline the individual
> contours, so let's make every point black if it is
> not the same as its neighbours. A simple heuristic
> is if the sum of the four points in a square divided by
> four is an integer, then you are in the middle of a
> contour, otherwise you are on the boundary. We set the
> boundary points to colour 0 which is before the beginning
> of the colormap, then we add black to the start of the
> colormap and redraw it:
>
> Q = filter2(ones(2)/4,imagesc(filter2(B,Z,"same")));
> Q(Q != fix(Q)) = 0;
> colormap([0,0,0; colormap]);
> image(Q+1);
>
> Now that we have found the contour lines, let's do a gradient
> fill on the original noisy data and plot the contours on top
> of it:
>
> colormap(hsv(64));
> M = imagesc(Z);
> M(Q==0) = 0;
> colormap([0,0,0; colormap]);
> image(M+1);
>
> You can use epstk to make a pretty plot with axis markers, etc.
>
> Hope this helps.
>
> Paul Kienzle
> address@hidden
>
> On Thu, Aug 29, 2002 at 05:47:12PM -0300, Afonso de Moraes Paiva wrote:
> > Hi
> > I am used to matlab but I am new to octave. Is there a way of making
> > decent filled-contour plots with octave? I've got a fill3.m script from
> > "somewhere" but can't make it work. Thanks for any help.
> > Afonso
> > --
> >
> > Afonso de Moraes Paiva
> > Programa de Engenharia Oceanica (PEnO/COPPE)
> > Universidade Federal do Rio de Janeiro (UFRJ)
> > Cidade Universitaria, CT, sala C203
> > Rio de Janeiro, RJ, 21945-970, Brasil
> > phone: (021) 2562-8729
> > email: address@hidden
> >
> >
> >
> > -------------------------------------------------------------
> > Octave is freely available under the terms of the GNU GPL.
> >
> > Octave's home on the web: http://www.octave.org
> > How to fund new projects: http://www.octave.org/funding.html
> > Subscription information: http://www.octave.org/archive.html
> > -------------------------------------------------------------
> >
>
>
>
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
>
> Octave's home on the web: http://www.octave.org
> How to fund new projects: http://www.octave.org/funding.html
> Subscription information: http://www.octave.org/archive.html
> -------------------------------------------------------------
>
-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------