octave_packages/m/plot/surfl.m

   1 ## Copyright (C) 2009-2012 Kai Habel
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {} surfl (@var{x}, @var{y}, @var{z})
  21 ## @deftypefnx {Function File} {} surfl (@var{z})
  22 ## @deftypefnx {Function File} {} surfl (@var{x}, @var{y}, @var{z}, @var{L})
  23 ## @deftypefnx {Function File} {} surfl (@var{x}, @var{y}, @var{z}, @var{L}, @var{P})
  24 ## @deftypefnx {Function File} {} surfl (@dots{}, "light")
  25 ## Plot a lighted surface given matrices @var{x}, and @var{y} from
  26 ## @code{meshgrid} and
  27 ## a matrix @var{z} corresponding to the @var{x} and @var{y} coordinates of
  28 ## the mesh.  If @var{x} and @var{y} are vectors, then a typical vertex
  29 ## is (@var{x}(j), @var{y}(i), @var{z}(i,j)).  Thus, columns of @var{z}
  30 ## correspond to different @var{x} values and rows of @var{z} correspond
  31 ## to different @var{y} values.
  32 ##
  33 ## The light direction can be specified using @var{L}.  It can be
  34 ## given as 2-element vector [azimuth, elevation] in degrees or as 3-element
  35 ## vector [lx, ly, lz].
  36 ## The default value is rotated 45° counter-clockwise from the current view.
  37 ##
  38 ## The material properties of the surface can specified using a 4-element vector
  39 ## @var{P} = [@var{AM} @var{D} @var{SP} @var{exp}] which defaults to
  40 ## @var{p} = [0.55 0.6 0.4 10].
  41 ## @table @asis
  42 ## @item "AM" strength of ambient light
  43 ##
  44 ## @item "D" strength of diffuse reflection
  45 ##
  46 ## @item "SP" strength of specular reflection
  47 ##
  48 ## @item "EXP" specular exponent
  49 ## @end table
  50 ##
  51 ## The default lighting mode "cdata", changes the cdata property to give the
  52 ## impression
  53 ## of a lighted surface.  Please note: the alternative "light" mode, which
  54 ## creates a light
  55 ## object to illuminate the surface is not implemented (yet).
  56 ##
  57 ## Example:
  58 ##
  59 ## @example
  60 ## @group
  61 ## colormap (bone (64));
  62 ## surfl (peaks);
  63 ## shading interp;
  64 ## @end group
  65 ## @end example
  66 ## @seealso{surf, diffuse, specular, surface}
  67 ## @end deftypefn
  68
  69 ## Author: Kai Habel <kai.habel@gmx.de>
  70
  71 function retval = surfl (varargin)
  72
  73   [h, varargin] = __plt_get_axis_arg__ ("surfl", varargin{:});
  74
  75   oldh = gca ();
  76   unwind_protect
  77     axes (h);
  78     newplot ();
  79
  80     ## Check for lighting type.
  81     use_cdata = true;
  82     if (ischar (varargin{end}))
  83       lstr = varargin{end};
  84       if (strncmp (tolower (lstr), "light", 5))
  85         warning ("light method not supported (yet), using cdata method instead");
  86         ## This can be implemented when light objects are supported.
  87         use_cdata = false;
  88       elseif (strncmp (tolower (lstr), "cdata", 5))
  89         use_cdata = true;
  90       else
  91         error ("surfl: unknown lighting method");
  92       endif
  93       varargin(end) = [];
  94     endif
  95
  96     ## Check for reflection properties argument.
  97     ##
  98     ## r = [ambient light strength,
  99     ##      diffuse reflection strength,
 100     ##      specular reflection strength,
 101     ##      specular shine]
 102     if (length (varargin{end}) == 4 && isnumeric (varargin{end}))
 103       r = varargin{end};
 104       varargin(end) = [];
 105     else
 106       ## Default values.
 107       r = [0.55, 0.6, 0.4, 10];
 108     endif
 109
 110     ## Check for light vector (lv) argument.
 111     have_lv = false;
 112     if (isnumeric (varargin{end}))
 113       len = numel (varargin{end});
 114       lastarg = varargin{end};
 115       if (len == 3)
 116         lv = lastarg;
 117         varargin(end) = [];
 118         have_lv = true;
 119       elseif (len == 2)
 120         [lv(1), lv(2), lv(3)] = sph2cart ((lastarg(1) - 90) * pi/180, lastarg(2) * pi/180, 1.0);
 121         varargin(end) = [];
 122         have_lv = true;
 123       endif
 124     endif
 125
 126     tmp = surface (varargin{:});
 127     if (! ishold ())
 128       set (h, "view", [-37.5, 30],
 129            "xgrid", "on", "ygrid", "on", "zgrid", "on", "clim", [0 1]);
 130     endif
 131
 132     ## Get view vector (vv).
 133     a = axis;
 134     [az, el] = view;
 135     [vv(1), vv(2), vv(3)] = sph2cart ((az - 90) * pi/180.0, el * pi/180.0, 1.0);
 136     vv /= norm (vv);
 137
 138     if (!have_lv)
 139       ## Calculate light vector (lv) from view vector.
 140       Phi = 45.0 / 180.0 * pi;
 141       R = [cos(Phi), -sin(Phi), 0;
 142            sin(Phi),  cos(Phi), 0;
 143            0,          0,         1];
 144       lv = (R * vv.').';
 145     endif
 146
 147     vn = get (tmp, "vertexnormals");
 148     dar = get (h, "plotboxaspectratio");
 149     vn(:,:,1) *= dar(1);
 150     vn(:,:,2) *= dar(2);
 151     vn(:,:,3) *= dar(3);
 152
 153     ## Normalize vn.
 154     vn = vn ./ repmat (sqrt (sumsq (vn, 3)), [1, 1, 3]);
 155     [nr, nc] = size(get(tmp, "zdata"));
 156
 157     ## Ambient, diffuse, and specular term.
 158     cdata = (r(1) * ones (nr, nc)
 159              + r(2) * diffuse  (vn(:,:,1), vn(:,:,2), vn(:,:,3), lv)
 160              + r(3) * specular (vn(:,:,1), vn(:,:,2), vn(:,:,3), lv, vv, r(4)));
 161
 162     set (tmp, "cdata", cdata ./ sum (r(1:3)));
 163
 164   unwind_protect_cleanup
 165     axes (oldh);
 166   end_unwind_protect
 167
 168   if (nargout > 0)
 169     retval = tmp;
 170   endif
 171
 172 endfunction
 173
 174 %!demo
 175 %! clf
 176 %! [X,Y,Z]=sombrero;
 177 %! colormap(copper);
 178 %! surfl(X,Y,Z);
 179 %! shading interp;
 180
 181 %!demo
 182 %! clf
 183 %! [X,Y,Z]=sombrero;
 184 %! colormap(copper);
 185 %! [az, el] = view;
 186 %! surfl(X,Y,Z,[az+225,el],[0.2 0.6 0.4 25]);
 187 %! shading interp;
 188