scilab/modules/graphics/macros/contourf.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) INRIA
   3 // This file must be used under the terms of the CeCILL.
   4 // This source file is licensed as described in the file COPYING, which
   5 // you should have received as part of this distribution.  The terms
   6 // are also available at
   7 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
   8
   9 function contourf(x,y,z,nv,style,strf,leg,rect,nax)
  10
  11         [nout,nin]=argn(0);
  12
  13         if nin == 0 then   // demo
  14                 t = -%pi:0.1:%pi;
  15                 m = sin(t)'*cos(t);
  16                 contourf(t,t,m);
  17                 return;
  18   end
  19
  20 if nin <= 0 then x=1:10;end
  21 if nin <= 1 then y=1:10;end
  22 if nin <= 2 then z=rand(size(x,'*'),size(y,'*'));end
  23 if nin <= 3 then zmin=mini(z);zmax=maxi(z);nv = zmin + (1:10)*(zmax-zmin)/(11);end
  24 if nin <= 5 then strf="121";end
  25 if nin <= 6 then leg=" ";end
  26 if nin <= 7 then rect=[0,0,1,1];end
  27 if nin <= 8 then nax=[1,10,1,10];end
  28 if x==[] then x=1:size(z,'r');end
  29 if y==[] then y=1:size(z,'c');end
  30
  31 nvs=size(nv,'*') ;
  32 if nvs==1 then nvs=nv;zmin=mini(z);zmax=maxi(z);nv = zmin + (1:nvs)*(zmax-zmin)/(nvs+1);end;
  33 if nin <= 4 then style = -1*ones(1,nvs);end
  34 if nin <= 7 then rect=[mini(x),mini(y),maxi(x),maxi(y)]; end
  35 nv1=nv
  36 [mz,nz] = size(z);
  37 minz = min(z);
  38 maxz = max(z);
  39 // Surround the matrix by a very low region to get closed contours, and
  40 // replace any NaN with low numbers as well.
  41 zz=[ %nan*ones(1,nz+2); %nan*ones(mz,1),z,%nan*ones(mz,1);%nan*ones(1,nz+2)];
  42 kk=find(isnan(zz(:)));
  43 zz(kk)=minz-1e4*(maxz-minz)+zeros(kk);
  44
  45 xx = [2*x(1)-x(2); x(:); 2*x(mz)-x(mz-1)];
  46 yy = [2*y(1)-y(2); y(:); 2*y(nz)-y(nz-1)];
  47
  48 // Internal call to get the contours
  49 [x1,y1]=contour2di(xx,yy,zz,nv);
  50 CS=[x1;y1];
  51 // Find the indices of the curves in the c matrix, and get the
  52 // area of closed curves in order to draw patches correctly.
  53 ii = 1;
  54 ncurves = 0;
  55 I = [];
  56 Area=[];
  57
  58 while (ii < size(CS,2)),
  59   nl=CS(2,ii);
  60   ncurves = ncurves + 1;
  61   I(ncurves) = ii;
  62   xp=CS(1,ii+(1:nl));  // First patch
  63   yp=CS(2,ii+(1:nl));
  64   Area(ncurves)=sum( mtlb_diff(xp).*(yp(1:nl-1)+yp(2:nl))/2 );
  65   ii = ii + nl + 1;
  66 end
  67
  68 lp=xget('lastpattern');
  69
  70 if size(nv,'*') > 1 // case where nv is a vector defining the level curve values
  71   if  size(nv,'*') > lp ; write(%io(2),'Colormap too small');return ;end
  72 else
  73   if nv > lp ; write(%io(2),'Colormap too small');return ;end
  74 end
  75
  76 min_nv=mini(nv);
  77 max_nv=maxi(nv);
  78
  79 plot2d([mini(xx);maxi(xx)],[mini(yy);maxi(yy)],0,strf,leg,rect,nax);
  80
  81 // Plot patches in order of decreasing size. This makes sure that
  82 // all the lev1es get drawn, not matter if we are going up a hill or
  83 // down into a hole. When going down we shift levels though, you can
  84 // tell whether we are going up or down by checking the sign of the
  85 // area (since curves are oriented so that the high side is always
  86 // the same side). Lowest curve is largest and encloses higher data
  87 // always.
  88
  89 draw_min=1;
  90 H=[];
  91 [FA,IA]=sort(abs(Area));
  92
  93   drawlater(); // postpon the drawing here
  94   a=gca();
  95   old_foreground = a.foreground;
  96   pat=xget('pattern');
  97   for jj=IA',
  98     nl=CS(2,I(jj));
  99     lev1=CS(1,I(jj));
 100     if (lev1 ~= minz | draw_min),
 101       xp=CS(1,I(jj)+(1:nl));
 102       yp=CS(2,I(jj)+(1:nl));
 103       pat=size(find( nv <= lev1),'*');
 104       xset("pattern",pat);
 105       xfpoly(xp,yp)
 106     end;
 107   end
 108
 109   if style(1)<>-1 then
 110     contour2d(xx,yy,zz,nv,style,"000",leg,rect,nax);
 111   end
 112   a.foreground = old_foreground;
 113   drawnow(); // draw all now!
 114
 115 endfunction