New plotting function surfl (again)

[Kai Habel]

Hello all,

this patch is also outstanding [1].

> the attached changeset adds the new plotting function surfl. To work
> properly it requires the previous changeset about surface normals to be
> applied as well.
>
> The function surfl plots lighted surfaces. For best visual experience
> colormaps like pink, copper, bone, or gray should be used.
>
> For example:
>
> colormap(bone);
> surfl(peaks);
> shading interp; 
I have used the time since then, to revise the patch. The new patch also
includes the functions specular and diffuse, which perform the light
calculation now. All three are matlab core functions.

Kai

[1] http://www.nabble.com/New-plotting-function-surfl-td20407867.html

# HG changeset patch
# User Kai Habel
# Date 1231426310 -3600
# Node ID 4c149a136ffb5508b96229dc3d0039dbcb83a7b7
# Parent  f53edfb28c12db747388c643b20bab1b334e7291
Add new 3D plotting function surfl. Add light functions diffuse and specular

diff -r f53edfb28c12 -r 4c149a136ffb scripts/ChangeLog
--- a/scripts/ChangeLog Thu Jan 08 14:35:22 2009 +0100
+++ b/scripts/ChangeLog Thu Jan 08 15:51:50 2009 +0100
@@ -1,3 +1,9 @@
+2009-01-08  Kai Habel <address@hidden>
+
+        * plot/surfl.m: New function
+        * plot/diffuse.m: dito
+        * plot/specular.m: dito
+
 2008-12-28  Jaroslav Hajek <address@hidden>
 
        * miscellaneous/delete.m: Allow filename globs. Display warnings if
@@ -15,6 +21,10 @@
 
        * plot/hist.m: Doc string now mentions matrix input argument.
        Correct error message.
+
+2008-12-30  Ben Abbott <address@hidden>
+
+       * plot/__contour__.m: __contour__.m: correct order of patches
 
 2008-12-24  Doug Stewart  <address@hidden>
 
diff -r f53edfb28c12 -r 4c149a136ffb scripts/plot/diffuse.m
--- /dev/null   Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/diffuse.m    Thu Jan 08 15:51:50 2009 +0100
@@ -0,0 +1,59 @@
+## Copyright (C) 2009 Kai Habel
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} diffuse (@var{sx}, @var{sy}, @var{sz}, 
@var{l})
+## Calculate diffuse reflection strength of a surface defined by the normal
+## vector elements @var{sx}, @var{sy}, @var{sz}. 
+## The light vector can be specified using parameter @var{L}. It can be
+## given as 2-element vector [azimuth, elevation] in degrees or as 3-element
+## vector [lx, ly, lz]. 
+## @seealso{specular, surfl}
+## @end deftypefn
+
+## Author: Kai Habel <address@hidden>
+
+function retval = diffuse (sx, sy, sz, lv)
+
+  ## general checks
+  if (nargin != 4)
+    usage ("number of arguments must be 4")
+  endif
+
+  ## check for normal vector
+  if (!size_equal (sx, sy, sz))
+    usage ("SX, SY, and SZ must have same size")
+  endif
+  
+  ## check for light vector (lv) argument
+  if (length (lv) < 2 || length (lv) > 3)
+    usage ("light vector LV must be a 2- or 3-element vector");
+  elseif (length (lv) == 2)
+    [lv(1), lv(2), lv(3)] = sph2cart (lv(1) * pi/180, lv(2) * pi/180, 1.0);
+  endif
+
+  ## normalize view and light vector
+  if (sum (abs (lv)) > 0)
+    lv  /= norm (lv);
+  endif
+
+  ns = sqrt (sx.^2 + sy.^2 + sz.^2);
+  retval = (sx * lv(1) + sy * lv(2) + sz * lv(3)) ./ ns;
+  retval(retval < 0) = 0;
+  
+endfunction
diff -r f53edfb28c12 -r 4c149a136ffb scripts/plot/specular.m
--- /dev/null   Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/specular.m   Thu Jan 08 15:51:50 2009 +0100
@@ -0,0 +1,90 @@
+## Copyright (C) 2009 Kai Habel
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} specular (@var{sx}, @var{sy}, @var{sz}, 
@var{l}, @var{v})
+## @deftypefnx {Function File} {} specular (@var{sx}, @var{sy}, @var{sz}, 
@var{l}, @var{v}, @var{se})
+## Calculate specular reflection strength of a surface defined by the normal
+## vector elements @var{sx}, @var{sy}, @var{sz} using Phong's approximation. 
+## The light and view vectors can be specified using parameter @var{L} and 
@var{V} respectively.
+## Both can be given as 2-element vectors [azimuth, elevation] in degrees or 
as 3-element
+## vector [x, y, z]. An optional 6th argument describes the specular exponent 
(spread) @var{se}.
+## @seealso{surfl, diffuse}
+## @end deftypefn
+
+## Author: Kai Habel <address@hidden>
+
+function retval = specular (sx, sy, sz, lv, vv, se)
+
+  ## general checks
+  if ((nargin < 5) || (nargin > 6))
+    usage ("number of arguments must be 5 or 6")
+  endif
+
+  ## checks for specular exponent (se)
+  if (nargin < 6)
+    se = 10;
+  else
+    if (!isnumeric (se) || (numel (se) != 1) || (se <= 0))
+      usage ("specular exponent must be positive scalar");
+    endif
+  endif
+
+  ## checks for normal vector
+  if (!size_equal (sx, sy, sz))
+    usage ("SX, SY, and SZ must have same size")
+  endif
+  
+  ## check for light vector (lv) argument
+  if (length (lv) < 2 || length (lv) > 3)
+    usage ("light vector LV must be a 2- or 3-element vector");
+  elseif (length (lv) == 2)
+    [lv(1), lv(2), lv(3)] = sph2cart (lv(1) * pi/180, lv(2) * pi/180, 1.0);
+  endif
+
+  ## check for view vector (vv) argument
+  if ((length (vv) < 2) || (length (lv) > 3))
+    error ("view vector VV must be a 2- or 3-element vector");
+  elseif (length (vv) == 2)
+    [vv(1), vv(2), vv(3)] = sph2cart (vv(1) * pi / 180, vv(2) * pi / 180, 1.0);
+  endif
+
+  ## normalize view and light vector
+  if (sum (abs (lv)) > 0)
+    lv  /= norm (lv);
+  endif
+  if (sum (abs (vv)) > 0)
+    vv  /= norm (vv);
+  endif
+
+  ## calculate normal vector lengths and dot-products
+  ns = sqrt (sx.^2 + sy.^2 + sz.^2);
+  l_dot_n = (sx * lv(1) + sy * lv(2) + sz * lv(3)) ./ ns;
+  v_dot_n = (sx * vv(1) + sy * vv(2) + sz * vv(3)) ./ ns;
+
+  ## calculate specular reflection using Phong's approximation
+  retval = 2 * l_dot_n .* v_dot_n - dot (lv, vv);
+  
+  ## set zero if light is on the other side
+  retval(l_dot_n < 0) = 0;
+
+  ## allow postive values only
+  retval(retval < 0) = 0;
+  retval = retval .^ se;
+  
+endfunction
diff -r f53edfb28c12 -r 4c149a136ffb scripts/plot/surfl.m
--- /dev/null   Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/surfl.m      Thu Jan 08 15:51:50 2009 +0100
@@ -0,0 +1,178 @@
+## Copyright (C) 2009 Kai Habel
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} surfl (@var{x}, @var{y}, @var{z})
+## @deftypefnx {Function File} {} surfl (@var{z})
+## @deftypefnx {Function File} {} surfl (@var{x}, @var{y}, @var{z}, @var{L})
+## @deftypefnx {Function File} {} surfl (@var{x}, @var{y}, @var{z}, @var{L}, 
@var{P})
+## @deftypefnx {Function File} {} surfl (...,"light")
+## Plot a lighted surface given matrices @var{x}, and @var{y} from 
@code{meshgrid} and
+## a matrix @var{z} corresponding to the @var{x} and @var{y} coordinates of
+## the mesh.  If @var{x} and @var{y} are vectors, then a typical vertex
+## is (@var{x}(j), @var{y}(i), @var{z}(i,j)).  Thus, columns of @var{z}
+## correspond to different @var{x} values and rows of @var{z} correspond
+## to different @var{y} values.
+##
+## The light direction can be specified using @var{L}. It can be
+## given as 2-element vector [azimuth, elevation] in degrees or as 3-element 
vector [lx, ly, lz].
+## The default value is rotated 45° counter-clockwise from the current view.
+##
+## The material properties of the surface can specified using a 4-element 
vector
+## @var{P} = address@hidden @var{D} @var{SP} @var{exp}] which defaults to
+## @var{p} = [0.55 0.6 0.4 10]. 
+## @table @code
+## @item "AM" strength of ambient light
+## @item "D" strength of diffuse reflection
+## @item "SP" strength of specular reflection
+## @item "EXP" specular exponent
+## @end table
+## 
+## The default lighting mode "cdata", changes the cdata property to give the 
impression
+## of a lighted surface. Please note: the alternative "light" mode, which 
creates a light
+## object to iluminate the the surface is not implemented (yet).
+##
+## Example:
+##
+## @example
+## colormap(bone);
+## surfl(peaks);
+## shading interp;
+## @end example
+## @seealso{surf, diffuse, specular, surface}
+## @end deftypefn
+
+## Author: Kai Habel <address@hidden>
+
+function retval = surfl (varargin)
+
+  [h, varargin] = __plt_get_axis_arg__ ("surfl", varargin{:});
+
+  oldh = gca ();
+  unwind_protect
+    axes (h);
+    newplot ();
+
+    ## check for lighting type
+    use_cdata = true;
+    if (ischar (varargin{end}))
+      lstr = varargin{end};
+      if strncmp (tolower (lstr), "light", 5)
+        warning ("light method not supported (yet), using cdata method 
instead");
+        # this can be implemented when light objects are being
+        # supported.
+        use_cdata = false;
+      elseif (strncmp (tolower (lstr), "cdata", 5))
+        use_cdata = true;
+      else
+        usage ("unknown lighting method");
+      endif
+      varargin(end) = [];
+    endif
+
+    ## check for reflection properties argument
+    ## r = [ambient light strength,
+    ##      diffuse reflection strength,
+    ##      specular reflection strength,
+    ##      specular shine] 
+    if ((length (varargin{end}) == 4) && isnumeric (varargin{end}))
+      r = varargin{end};
+      varargin(end) = [];
+    else
+      ## default values
+      r = [0.55, 0.6, 0.4, 10];
+    endif
+    
+
+    ## check for light vector (lv) argument
+    have_lv = false;
+    if (isnumeric (varargin{end}))
+      len = numel (varargin{end});
+      lastarg = varargin{end};
+      if (len == 3)
+        lv = lastarg;
+        varargin(end) = [];
+        have_lv = true;
+      elseif (len == 2)
+        [lv(1), lv(2), lv(3)] = sph2cart ((lastarg(1) - 90) * pi/180, 
lastarg(2) * pi/180, 1.0);
+        varargin(end) = [];
+        have_lv = true;
+      endif
+    endif
+    
+    tmp = surface (varargin{:});
+    if (! ishold ())
+      set (h, "view", [-37.5, 30],
+          "xgrid", "on", "ygrid", "on", "zgrid", "on", "clim", [0 1]);
+    endif
+
+    ## get view vector (vv)
+    a = axis;
+    [az, el] = view;
+    [vv(1), vv(2), vv(3)] = sph2cart ((az - 90) * pi/180.0, el * pi/180.0, 
1.0);
+    vv /= norm (vv);
+
+    if (!have_lv)
+      ## calculate light vector (lv) from view vector
+      Phi = 45.0 / 180.0 * pi;
+      R = [cos(Phi), -sin(Phi), 0;\
+           sin(Phi),  cos(Phi), 0;\
+           0,          0,         1];
+      lv = (R * vv.').';
+    endif
+
+    vn = get (tmp, "vertexnormals");
+    dar = get (h, "dataaspectratio");
+    vn(:, :, 1) *= dar(1);
+    vn(:, :, 2) *= dar(2);
+    vn(:, :, 3) *= dar(3);
+
+    ## normalize vn
+    vn = vn ./ repmat (sqrt (sumsq (vn, 3)), [1, 1, 3]);
+    [nr, nc] = size(get(tmp, "zdata"));
+
+    ## ambient, diffuse, and specular term
+    cdata = r(1) * ones (nr, nc) \
+          + r(2) * diffuse  (vn(:, :, 1), vn(:, :, 2), vn(:, :, 3), lv) \
+          + r(3) * specular (vn(:, :, 1), vn(:, :, 2), vn(:, :, 3), lv, vv, 
r(4));
+
+    set (tmp, "cdata", cdata ./ sum (r(1:3)));
+    
+  unwind_protect_cleanup
+    axes (oldh);
+  end_unwind_protect
+
+  if (nargout > 0)
+    retval = tmp;
+  endif
+
+endfunction
+
+%!demo
+%! [X,Y,Z]=sombrero;
+%! colormap(copper);
+%! surfl(X,Y,Z);
+%! shading interp;
+
+%!demo
+%! [X,Y,Z]=sombrero;
+%! colormap(copper);
+%! [az, el] = view;
+%! surfl(X,Y,Z,[az+225,el],[0.2 0.6 0.4 25]);
+%! shading interp;
+