Re: contour plot

[Paul Kienzle]

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