home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / PLOTXYZ.DI$ / SURFL.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.4 KB  |  114 lines
  1. function hout=surfl(x,y,z,s,k)
  2. %SURFL    3-D shaded surface with lighting.
  3. %    SURFL(...) is the same as SURF(...) except that it draws the surface
  4. %    with highlights from a light source.
  5. %
  6. %    SURFL(Z), SURFL(X,Y,Z), SURFL(Z,S), and SURFL(X,Y,Z,S) are all 
  7. %    legal. S, if specified, is the three vector S = [Sx,Sy,Sz]
  8. %    that specifies the direction of the light source. S can also be
  9. %    specified in spherical coordinates, S = [AZ,EL].
  10. %
  11. %    The shading is based on a combination of diffuse, specular and 
  12. %    ambient lighting models.
  13. %
  14. %    The default value for S is 45 degrees counterclockwise from
  15. %    the current view direction.  Use CLA, HOLD ON, VIEW(AZ,EL),
  16. %    SURFL(...), HOLD OFF to plot the lighted surface with view
  17. %    direction (AZ,EL).
  18. %
  19. %    The relative contributions due to ambient light, diffuse
  20. %    reflection, specular reflection, and the specular spread
  21. %    coefficient can be set by using five arguments
  22. %    SURFL(X,Y,Z,S,K) where K=[ka,kd,ks,spread].
  23. %
  24. %    Relies on the ordering of points in the X,Y, and Z matrices
  25. %    to define the inside and outside of parametric surfaces.
  26. %    Try SURFL(X',Y',Z') if you don't like the results of
  27. %    this function.  Due to the way surface normal vectors are
  28. %    computed, SURFL requires matrices that are at least 3-by-3.
  29. %
  30. %    See also SURF.
  31.  
  32. %    Clay M. Thompson 4-24-91
  33. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  34.  
  35. error(nargchk(1,5,nargin));
  36.  
  37. if nargin==1,    % Define x,y
  38.   z = x;
  39.   [m,n] = size(z);
  40.   [x,y] = meshgrid(1:n,1:m);
  41.  
  42. elseif nargin==2,
  43.   z = x;
  44.   s = y;
  45.   [m,n] = size(z);
  46.   [x,y] = meshgrid(1:n,1:m);
  47. end
  48. if isstr(z) | isstr(x)
  49.     error('Must be 1, 3 or 4 numeric input arguments.')
  50. end
  51.  
  52. if nargin<5 % Define default weighting coefficients
  53.   k = [.55,.6,.4,10]; % Ambient,diffuse,specular,spread
  54. end 
  55. if length(k)~=4,
  56.   error('Weighting vector k must have four components.');
  57. end
  58.  
  59. [msg,x,y,z] = xyzchk(x,y,z); if ~isempty(msg), error(msg); end
  60. if any(size(z)<[3 3]), error('X, Y, and Z must be at least 3-by-3.'); end
  61.  
  62. stencil=ones(2,2)/4;
  63.  
  64. cax = newplot;
  65. next = lower(get(cax,'NextPlot'));
  66. if ~ishold, view(3); end    % Set graphics system for 3-D plot
  67. [vaz,vel] = view;
  68. vaz = vaz*pi/180; vel = vel*pi/180; % Convert to radians
  69.  
  70. if (nargin==1) | (nargin==3), % Use default S
  71.   phi = 45*pi/180;
  72.   s = zeros(1,3);
  73.   s(1) = cos(vaz)*sin(phi)+sin(vaz)*cos(vel)*cos(phi);
  74.   s(2) = sin(phi)*sin(vaz)-cos(vaz)*cos(vel)*cos(phi);
  75.   s(3) = sin(phi)*sin(vel);
  76. else
  77.   if (length(s)~=2) & (length(s)~=3),
  78.     error('S must be specified using [AZ,EL] or [Sx,Sy,Sz].');
  79.   end
  80. end
  81.  
  82. ms = length(s(:));
  83. if ms==2, % Compute source direction from [AZ,EL]
  84.   az = s(1)*pi/180; el = s(2)*pi/180; % Convert to radians
  85.   s = zeros(1,3);
  86.   s(1) =  sin(az)*cos(el);
  87.   s(2) = -cos(az)*cos(el);
  88.   s(3) =  sin(el);
  89. end
  90.  
  91. % Determine plot scaling factors for a cube-like plot domain.
  92. h = surf(x,y,z);
  93. a = [get(gca,'xlim') get(gca,'ylim') get(gca,'zlim')];
  94. Sx = a(2)-a(1);
  95. Sy = a(4)-a(3);
  96. Sz = a(6)-a(5);
  97. scale = max([Sx,Sy,Sz]);
  98. Sx = Sx/scale; Sy = Sy/scale; Sz = Sz/scale;
  99.  
  100. % Compute surface normals.  Rely on ordering to define inside or outside.
  101. xx = x/Sx; yy = y/Sy; zz = z/Sz;
  102. [nx,ny,nz] = surfnorm(xx,yy,zz);
  103.  
  104. % Compute Lambertian shading + specular + ambient light
  105. R = (k(1)+k(2)*diffuse(nx,ny,nz,s)+ ...
  106.     k(3)*specular(nx,ny,nz,s,[vaz,vel],k(4)))/ sum(k(1:3));
  107.  
  108. % Show surface
  109. set(h,'cdata',R);
  110. caxis([0,1]);     % Set color axis range.
  111. if nargout > 0
  112.     hout = h;
  113. end
  114.