New plotting function surfl

[kahacjde]

Hello all,

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;

The function should work pretty well, so please give it try and report
bugs.

Kai

# HG changeset patch
# User Kai Habel
# Date 1226245327 -3600
# Node ID cacabddbb04affd8ea987eb6528813f8f47a4891
# Parent  36c59eb57b3764c762d273499c2730c47093ab04
Add new 3D plotting function surfl

diff -r 36c59eb57b37 -r cacabddbb04a scripts/ChangeLog
--- a/scripts/ChangeLog Sun Nov 09 15:58:24 2008 +0100
+++ b/scripts/ChangeLog Sun Nov 09 16:42:07 2008 +0100
@@ -1,3 +1,7 @@
[phone]  Kai Habel <kai.habel@...>
+
+ * plot/surfl.m: New function file
+
 2008-09-28  Jaroslav Hajek <highegg@...>
 
  * optimization/__fdjac__.m: New function file.
diff -r 36c59eb57b37 -r cacabddbb04a scripts/plot/surfl.m
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/surfl.m Sun Nov 09 16:42:07 2008 +0100
@@ -0,0 +1,192 @@
+## Copyright (C) 2008 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 45° counter-clockwise (using the z-axis) from the current view.
+##
+## The material properties of the surface can specified using a 4-element vector
+## @var{P} = [@var{AM} @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, surface}
+## @end deftypefn
+
+## Author: Kai Habel <kai.habel@...>
+
+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
+    # refl_prop = [ambient light,
+    #              diffuse reflection,
+    #              specular reflection,
+    #              specular shine]
+    if ((length(varargin{end}) == 4) && isnumeric(varargin{end}))
+      refl_prop = varargin{end};
+      varargin(end) = [];
+    else
+      # default values
+      refl_prop = [0.55 0.6 0.4 10];
+    endif
+    refl_norm = 1 / sum(refl_prop(1:3));
+
+    # 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)
+        [lx ly lz] = sph2cart(lastarg(1) * pi/180, lastarg(2) * pi/180, 1.0);
+        lv = [lx ly lz];
+        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;
+    [vx vy vz] = sph2cart((az-90.0) * pi/180.0, el * pi/180.0, 1.0);
+    vv = [vx vy vz];
+    vv = vv / norm(vv);
+
+    if (!have_lv)
+      # calculate light vector (lv) from view vector
+      phi = 45.0;
+      Phi = phi/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);
+
+    vn = vn ./ repmat(sqrt(sumsq(vn,3)),[1 1 3]);
+
+    # ambient term
+    [nr, nc] = size(get(tmp,'zdata'));
+    cdata = ones(nr,nc) * refl_prop(1);
+
+    # diffuse term
+    lmat = repmat(lv(:),[1 nr nc]);
+    lmat = shiftdim(lmat,1);
+    diffterm = dot (lmat, vn, 3);
+    diffterm(diffterm <= 0) = 0;
+    cdata += refl_prop(2) * diffterm;
+
+    # specular term
+    tmp_vec = (vv + lv) / norm(vv + lv);
+    tmp_mat = repmat(tmp_vec(:), [1 nr nc]);
+    tmp_mat = shiftdim(tmp_mat, 1);
+    specterm = dot(tmp_mat, vn, 3);
+    specterm(specterm <= 0) = 0;
+    specterm = specterm .^ refl_prop(4);
+    cdata += refl_prop(3) * specterm;
+
+    set(tmp, "cdata", cdata .* refl_norm);
+    
+  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);
+
+%!demo
+%! [X,Y,Z]=sombrero;
+%! colormap(copper);
+%! [az, el] = view;
+%! surfl(X,Y,Z,[az+225,el],[0.2 0.6 0.4 25]);
+

[Ben Abbott]

On Nov 9, 2008, at 10:54 AM, Kai Habel wrote:

> Hello all,
>
> 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;
>
> The function should work pretty well, so please give it try and report
> bugs.

Very cool!

It works for me on Mac OSX with Aquaterm as well as with x11. I'm  
unable to try wxt at the moment ... I broke it when I started running  
the developers sources for gnuplot :-(

To be sure, what *I* see is what *you* see, I've attached a pdf.

Ben

[kahacjde]

...

>
> Very cool!
>
> It works for me on Mac OSX with Aquaterm as well as with x11. I'm
> unable to try wxt at the moment ... I broke it when I started running
> the developers sources for gnuplot :-(
>
> To be sure, what *I* see is what *you* see, I've attached a pdf.
>
> Ben
>
>

Thanks for trying. Yes, that's how it should look like. If you have not
tried already, an 'interp shading' removes the grid lines and gives even
better graphics.

For the interested, the implementation was not that complicated. The
algorithms are well documented in the OpenGL literature (e.g. "Red Book"
OpenGL [1]).

Kai

[1] http://fly.srk.fer.hr/~unreal/theredbook/chapter06.html