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