Home > hiswa > contoursurf.m

contoursurf

PURPOSE ^

CONTOURSURF contouring over non-rectangular surface.

SYNOPSIS ^

function CS=contoursurf(arg1,arg2,arg3,arg4);

DESCRIPTION ^

  CONTOURSURF contouring over non-rectangular surface.
        This is an extension of contourc.
        CONTOURSURF calculates the contour matrix C for use by EXTCONTOUR
        to draw the actual contour plot.
        C = CONTOURSURF(Z) computes the contour matrix for a contour plot
        of matrix Z treating the values in Z as heights above a plane.
        C = CONTOURSURF(X,Y,Z), where X and Y are vectors, specifies the X- 
        and Y-axes for the contour computation. X and Y can also be matrices of 
        the same size as Z, in which case they specify a surface in an 
        identical manner as SURFACE.
        CONTOURSURF(Z,N) and CONTOURSURF(X,Y,Z,N) compute N contour lines, 
        overriding the default automatic value.
        CONTOURSURF(Z,V) and CONTOURSURF(X,Y,Z,V) compute LENGTH(V) contour 
        lines at the values specified in vector V.
  
        The contour matrix C is a two row matrix of contour lines. Each
        contiguous drawing segment contains the value of the contour, 
        the number of (x,y) drawing pairs, and the pairs themselves.  
        The segments are appended end-to-end as
  
            C = [level1 x1 x2 x3 ... level2 x2 x2 x3 ...;
                 pairs1 y1 y2 y3 ... pairs2 y2 y2 y3 ...]
  
        See also EXTCONTOUR.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function CS=contoursurf(arg1,arg2,arg3,arg4);
0002 %  CONTOURSURF contouring over non-rectangular surface.
0003 %        This is an extension of contourc.
0004 %        CONTOURSURF calculates the contour matrix C for use by EXTCONTOUR
0005 %        to draw the actual contour plot.
0006 %        C = CONTOURSURF(Z) computes the contour matrix for a contour plot
0007 %        of matrix Z treating the values in Z as heights above a plane.
0008 %        C = CONTOURSURF(X,Y,Z), where X and Y are vectors, specifies the X-
0009 %        and Y-axes for the contour computation. X and Y can also be matrices of
0010 %        the same size as Z, in which case they specify a surface in an
0011 %        identical manner as SURFACE.
0012 %        CONTOURSURF(Z,N) and CONTOURSURF(X,Y,Z,N) compute N contour lines,
0013 %        overriding the default automatic value.
0014 %        CONTOURSURF(Z,V) and CONTOURSURF(X,Y,Z,V) compute LENGTH(V) contour
0015 %        lines at the values specified in vector V.
0016 %
0017 %        The contour matrix C is a two row matrix of contour lines. Each
0018 %        contiguous drawing segment contains the value of the contour,
0019 %        the number of (x,y) drawing pairs, and the pairs themselves.
0020 %        The segments are appended end-to-end as
0021 %
0022 %            C = [level1 x1 x2 x3 ... level2 x2 x2 x3 ...;
0023 %                 pairs1 y1 y2 y3 ... pairs2 y2 y2 y3 ...]
0024 %
0025 %        See also EXTCONTOUR.
0026 
0027 % Author: R. Pawlowicz (IOS) rich@ios.bc.ca
0028 %         12/12/94
0029 
0030 
0031  
0032 if (nargin <=2 ),
0033  numarg_for_call='arg1';
0034  for ii=2:nargin,
0035   numarg_for_call=[numarg_for_call ',arg' int2str(ii)];
0036  end;
0037  zz=arg1;
0038 else
0039  numarg_for_call='arg3';
0040  for ii=4:nargin,
0041   numarg_for_call=[numarg_for_call ',arg' int2str(ii)];
0042  end;
0043  zz=arg3;
0044 end; 
0045 
0046 eval(['CS=contourc(' numarg_for_call ');']);
0047 [Ny,Nx]=size(zz);
0048  
0049 % Find data values and check curve orientation.
0050  
0051 ii=ones(1,size(CS,2));
0052 k=1;
0053 while (k < size(CS,2)),
0054   nl=CS(2,k);
0055   
0056   % Now this is a little bit of magic needed to make the filled contours
0057   % work. Essentially I draw the *closed* contours so that the "high" side is
0058   % always on the right. To test this, I take the cross product of the
0059   % first vector with a vector to a corner point and test the sign
0060   % against the elevation change. There are several special cases:
0061   % (1) If the contour line goes through a point (which happen when -Infs
0062   % are around), and (2) when the contour level equals the level on the high
0063   % side (this always seems to happen in 'simple test' cases!). We take
0064   % care of (1) by choosing other points, and we take care of (2) by adding
0065   % eps to the data before comparing with the contour data.
0066   
0067   if ( CS(:,k+1)==CS(:,k+nl) & nl>1 ),
0068     lev=CS(1,k);
0069     x1=CS(1,k+1); y1=CS(2,k+1);
0070     x2=CS(1,k+2); y2=CS(2,k+2);
0071     vx1=x2-x1; vy1=y2-y1;
0072     cpx=round(x1); cpy=round(y1);
0073     if ( [cpx cpy]==[x1 y1] ),
0074       cpx=round(x2); cpy=round(y2);
0075       if ( [cpx cpy]==[x2 y2]),
0076          if ( ~([cpx cpy]==round([x1 y1])) ),
0077            cpx=round(x1);
0078          else
0079            cpx=round(x1)+y2-y1;
0080            cpy=round(y1)-x2+x1;
0081          end;
0082        end;
0083     end;
0084     vx2=cpx-x1; vy2=cpy-y1;
0085 %    if (sign(zz(cpy,cpx)-lev+epslev)==0) disp('lev=0'); end;
0086 %    if (sign(vx1*vy2-vx2*vy1)==0) disp('cross=0'); end;
0087     if ( sign(zz(cpy,cpx)-lev+eps) == sign(vx1*vy2-vx2*vy1)  ),
0088       CS(:,k+[1:nl])=fliplr(CS(:,k+[1:nl]));
0089     end; 
0090   end;
0091   ii(k)=0;
0092   k=k+1+nl;
0093 end;
0094 
0095 % Data from integer coords to data coords. There are 3 cases
0096 % (1) Matrix X/Y
0097 % (2) Vector X/Y
0098 % (3) no X/Y. (do nothing);
0099 
0100 if (nargin>2 & min(size(arg1))>1 ),
0101  
0102  X=CS(1,ii)';   Y=CS(2,ii)';
0103  cX=ceil(X);    fX=floor(X);
0104  cY=ceil(Y);    fY=floor(Y);
0105  
0106  Ibl=cY+(fX-1)*Ny;    Itl=fY+(fX-1)*Ny;
0107  Itr=fY+(cX-1)*Ny;    Ibr=cY+(cX-1)*Ny;
0108  
0109  dy=cY-Y; dx=X-fX;
0110  
0111  % Correct for possible conflicts in matlabs [1 1 1 ] indexing. This
0112  % probably will *never* happen in real life, but turns up annoyingly
0113  % often in "simple" test cases.
0114  if (Nx*Ny == length(Ibl) ),
0115    Ibl=[1;Ibl];   Itl=[1;Itl];
0116    Itr=[1;Itr];   Ibr=[1;Ibr];
0117    dx=[0;dx];     dy=[0;dy];
0118    Csave=CS(:,1); 
0119    ii(1)=1;  
0120  end;
0121  
0122  CS(1,(ii)) = [ arg1(Ibl).*(1-dx).*(1-dy) + arg1(Itl).*(1-dx).*dy ...
0123              + arg1(Itr).*dx.*dy         + arg1(Ibr).*dx.*(1-dy)]';
0124  CS(2,(ii)) = [ arg2(Ibl).*(1-dx).*(1-dy) + arg2(Itl).*(1-dx).*dy ...
0125              + arg2(Itr).*dx.*dy         + arg2(Ibr).*dx.*(1-dy) ]';
0126 
0127  if (exist('Csave')),
0128   CS(:,1)=Csave;
0129  end;
0130  
0131 elseif (nargin>2 & min(size(arg1))==1 ),
0132  X=CS(1,ii);  Y=CS(2,ii);
0133  cX=ceil(X); fX=floor(X);
0134  cY=ceil(Y); fY=floor(Y);
0135 
0136  dy=cY-Y;    dx=X-fX;
0137 
0138  if (size(arg1,2)==1), 
0139    CS(1,ii)=[arg1(fX)'.*(1-dx)+arg1(cX)'.*dx];
0140  else
0141    CS(1,ii)=[arg1(fX).*(1-dx)+arg1(cX).*dx];
0142  end;
0143  
0144  if (size(arg2,2)==1), 
0145    CS(2,ii)=[arg2(fY)'.*dy+arg2(cY)'.*(1-dy)]; 
0146  else
0147    CS(2,ii)=[arg2(fY).*dy+arg2(cY).*(1-dy)]; 
0148  end;
0149  
0150 end;
0151              
0152               
0153

Generated on Wed 30-Nov-2005 15:36:39 by m2html © 2003