home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / COLOR.DI$ / SURFNORM.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.3 KB  |  111 lines
  1. function [nxout,nyout,nzout] = surfnorm(x,y,z)
  2. %SURFNORM 3-D surface normals.
  3. %    [Nx,Ny,Nz] = SURFNORM(X,Y,Z) returns the components of the 3-D 
  4. %    surface normal for the surface with components (X,Y,Z).  The 
  5. %    normal is unnormalized.  
  6. %
  7. %    [Nx,Ny,Nz] = SURFNORM(Z) returns the surface normal components
  8. %    for the surface Z.
  9. %
  10. %    Without lefthand arguments, SURFNORM(X,Y,Z) or SURFNORM(Z) 
  11. %    plots the surface with the normals emanating from it.  Normals 
  12. %    are not shown for surface elements that face away from the viewer.
  13. %
  14. %    The surface normals returned are based on a bicubic fit of
  15. %    the data.
  16.  
  17. %    Clay M. Thompson  1-15-91
  18. %    Revised 8-5-91, 9-17-91 by cmt.
  19. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  20.  
  21. if nargin==1,
  22.   z = x;
  23.   [m,n] = size(z);
  24.   [x,y] = meshgrid(1:n,1:m);
  25. elseif nargin==2,
  26.   error('Wrong number of input arguments.');
  27. end
  28.  
  29.  
  30. [m,n] = size(x);
  31. if size(y)~=[m,n], error('Y and Z must be the same size as X.'); end
  32. if size(z)~=[m,n], error('Y and Z must be the same size as X.'); end
  33.  
  34. stencil1 = [1 0 -1]/2;
  35. stencil2 =  [-1;0;1]/2;
  36.  
  37. if nargout==0, % If plotting, then scale to match plot aspect ratio.
  38.   % Determine plot scaling factors for a cube-like plot domain.
  39.   cax = newplot;
  40.   next = lower(get(cax,'NextPlot'));
  41.   hold_state =  ishold;
  42.   h = surf(x,y,z,'b');    % Blue surface
  43.   a = [get(gca,'xlim') get(gca,'ylim') get(gca,'zlim')];
  44.   Sx = a(2)-a(1);
  45.   Sy = a(4)-a(3);
  46.   Sz = a(6)-a(5);
  47.   scale = max([Sx,Sy,Sz]);
  48.   Sx = Sx/scale; Sy = Sy/scale; Sz = Sz/scale;
  49.  
  50.   % Scale surface
  51.   xx = x/Sx; yy = y/Sy; zz = z/Sz;
  52. else
  53.   xx = x; yy = y; zz = z;
  54. end
  55.  
  56. % Expand x,y,z so interpolation is valid at the boundaries.
  57. xx = [3*xx(1,:)-3*xx(2,:)+xx(3,:);xx;3*xx(m,:)-3*xx(m-1,:)+xx(m-2,:)];
  58. xx = [3*xx(:,1)-3*xx(:,2)+xx(:,3),xx,3*xx(:,n)-3*xx(:,n-1)+xx(:,n-2)];
  59. yy = [3*yy(1,:)-3*yy(2,:)+yy(3,:);yy;3*yy(m,:)-3*yy(m-1,:)+yy(m-2,:)];
  60. yy = [3*yy(:,1)-3*yy(:,2)+yy(:,3),yy,3*yy(:,n)-3*yy(:,n-1)+yy(:,n-2)];
  61. zz = [3*zz(1,:)-3*zz(2,:)+zz(3,:);zz;3*zz(m,:)-3*zz(m-1,:)+zz(m-2,:)];
  62. zz = [3*zz(:,1)-3*zz(:,2)+zz(:,3),zz,3*zz(:,n)-3*zz(:,n-1)+zz(:,n-2)];
  63.  
  64. rows = 2:m+1; cols = 2:n+1;
  65. ax = filter2(stencil1,xx); ax = ax(rows,cols);
  66. ay = filter2(stencil1,yy); ay = ay(rows,cols);
  67. az = filter2(stencil1,zz); az = az(rows,cols);
  68.  
  69. bx = filter2(stencil2,xx); bx = bx(rows,cols);
  70. by = filter2(stencil2,yy); by = by(rows,cols);
  71. bz = filter2(stencil2,zz); bz = bz(rows,cols);
  72.  
  73. nx = -(ay.*bz - az.*by);
  74. ny = -(az.*bx - ax.*bz);
  75. nz = -(ax.*by - ay.*bx);
  76.  
  77. if nargout==0,
  78.  
  79.   if ~hold_state, view(3); end    % Set graphics system for 3-D plot
  80.  
  81.   % Set the length of the surface normals
  82.   mag = sqrt(nx.*nx+ny.*ny+nz.*nz)*(10/scale);
  83.   d = (mag==0); mag(d) = eps*ones(size(mag(d)));
  84.   nx = nx ./mag;
  85.   ny = ny ./mag;
  86.   nz = nz ./mag;
  87.  
  88.   % Normal vector points
  89.   xc = x; yc = y; zc = z;
  90.  
  91.   hold on;
  92.   % use nan trick here
  93.   xp = [xc(:) Sx*nx(:)+xc(:) nan*ones(size(xc(:)))]';
  94.   yp = [yc(:) Sy*ny(:)+yc(:) nan*ones(size(yc(:)))]';
  95.   zp = [zc(:) Sz*nz(:)+zc(:) nan*ones(size(zc(:)))]';
  96.  
  97.   plot3(xp(:),yp(:),zp(:),'r-')
  98.  
  99.   if ~hold_state, set(cax,'NextPlot',next); end
  100.   return
  101.  
  102. end
  103.  
  104. % Normalize the length of the surface normals to 1.
  105. mag = sqrt(nx.*nx+ny.*ny+nz.*nz);
  106. d = (mag==0); mag(d) = eps*ones(size(mag(d)));
  107. nxout = nx ./mag;
  108. nyout = ny ./mag;
  109. nzout = nz ./mag;
  110.  
  111.