help-octave
[Top][All Lists]

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

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

```