octave_packages/plot-1.1.0/tricontour.m

   1 ## Copyright (C) 2008 Andreas Stahel <Andreas.Stahel@bfh.ch>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; if not, write to the Free Software
  15 ## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  16 ## 02110-1301  USA
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn {Function File} {} tricontour (@var{tri}, @var{x}, @var{y}, @var{z}, @var{levels})
  20 ## @deftypefnx {Function File} {} tricontour (@var{tri}, @var{x}, @var{y}, @var{z}, @var{levels}, @var{linespec})
  21 ## Plot level curves for the values of @code{@var{z}} on a triangular mesh in 2D.
  22 ##
  23 ## The variable @var{tri} is the triangular meshing of the points
  24 ## @code{(@var{x}, @var{y})} which is returned from @code{delaunay}. The
  25 ## variable @var{levels} is a vector with the values of the countour levels. If
  26 ## @var{levels} is a scalar, then it corresponds to the number of
  27 ## level curves to be drawn. If exactly one level curve is desired, list
  28 ## the level twice in the vector @var{levels}.
  29 ##
  30 ## If given, @var{linespec} determines the properties to use for the
  31 ## lines.
  32 ##
  33 ## The output argument @var{h} is the graphic handle to the plot.
  34 ## @seealso{plot, trimesh, delaunay}
  35 ## @end deftypefn
  36
  37 function h = tricontour (tri, x, y, z, levels, varargin)
  38   if (nargin < 5)
  39     print_usage ();
  40   endif
  41
  42   if isscalar(levels);
  43     dom=[min(z),max(z)];
  44     dom=mean(dom)+0.99*(dom-mean(dom));
  45     levels=linspace(dom(1),dom(2),levels);
  46   endif
  47
  48   levels=sort(levels);
  49   lmin=levels(1);
  50   lmax=levels(length(levels));
  51
  52   pData=[];                %% no preallocation
  53 %%  pData=zeros(12000,2);  %% preallocation
  54   pPoints=0;
  55
  56   for el=1:length(tri)
  57     values=[z(tri(el,1)),z(tri(el,2)),z(tri(el,3))];
  58     minval=min(values);
  59     maxval=max(values);
  60     locallevel=levels(minval<=levels+eps); # select the levels to be plotted
  61     if size(locallevel)>0
  62       locallevel=locallevel(locallevel<=maxval+eps);
  63     endif
  64     for level=locallevel
  65       points=zeros(1,2);
  66       npoints=1;
  67       dl=values-level;
  68
  69       if (abs(dl(1))<=10*eps)
  70         points(npoints,:)=[x(tri(el,1)),y(tri(el,1))];npoints++;endif
  71       if (abs(dl(2))<=10*eps)
  72         points(npoints,:)=[x(tri(el,2)),y(tri(el,2))];npoints++;endif
  73       if (abs(dl(3))<=10*eps)
  74         points(npoints,:)=[x(tri(el,3)),y(tri(el,3))];npoints++;endif
  75
  76       if (npoints<=2)
  77         if ((dl(1)*dl(2)) < 0)    # intersection between 1st and 2nd point
  78           points(npoints,:)= ( dl(2)*[x(tri(el,1)),y(tri(el,1))]...
  79                               -dl(1)*[x(tri(el,2)),y(tri(el,2))])/(dl(2)-dl(1));
  80           npoints++;
  81         endif
  82         if ((dl(1)*dl(3)) < 0)    # intersection between 1st and 3rd point
  83           points(npoints,:)= ( dl(3)*[x(tri(el,1)),y(tri(el,1))]...
  84                               -dl(1)*[x(tri(el,3)),y(tri(el,3))])/(dl(3)-dl(1));
  85           npoints++;
  86         endif
  87         if ((dl(3)*dl(2)) < 0)    # intersection between 2nd and 3rd point
  88           points(npoints,:)= ( dl(2)*[x(tri(el,3)),y(tri(el,3))]...
  89                               -dl(3)*[x(tri(el,2)),y(tri(el,2))])/(dl(2)-dl(3));
  90           npoints++;
  91         endif
  92       endif
  93       pData=[pData;points; NaN,NaN ];  %% no preallocation
  94   %%  pData(pPoints+1:pPoints+npoints,1:2)=[points; NaN,NaN ]; %% preallocation
  95       pPoints += npoints;
  96     endfor # level
  97   endfor # el
  98
  99   pData=pData(1:pPoints-1,:);
 100
 101   if (nargout>0)
 102     h= plot(pData(:,1),pData(:,2),varargin(:));
 103   else
 104     plot(pData(:,1),pData(:,2),varargin(:));
 105   endif
 106
 107 endfunction
 108
 109 %!demo
 110 %! rand ('state', 2)
 111 %! x = rand (100, 1)-0.5;
 112 %! y = rand (100, 1)-0.5;
 113 %! z= (x.*y);
 114 %! tri = delaunay (x, y);
 115 %! tricontour (tri, x, y, z, [-0.25:0.05:0.25]);
 116 %! axis equal
 117 %! grid on