Mercurial > octave-nkf

changeset 8450:6c4e2fc14beb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add new 3D plotting function surfl. Add light functions diffuse and specular
author Kai Habel
date Thu, 08 Jan 2009 15:51:50 +0100
parents 448188853722
children 6306ce1406f1
files scripts/ChangeLog scripts/plot/diffuse.m scripts/plot/specular.m scripts/plot/surfl.m
diffstat 4 files changed, 337 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Thu Jan 08 14:35:22 2009 +0100
+++ b/scripts/ChangeLog	Thu Jan 08 15:51:50 2009 +0100
@@ -25,6 +25,12 @@
 
 	* goemetry/voronoi.m: Speed up and handle dense grids.
 
+2009-01-08  Kai Habel <kai.habel@gmx.de>
+
+        * plot/surfl.m: New function
+        * plot/diffuse.m: dito
+        * plot/specular.m: dito
+
 2008-12-28  Jaroslav Hajek <highegg@gmail.com>
 
 	* miscellaneous/delete.m: Allow filename globs. Display warnings if
@@ -47,6 +53,10 @@
 
 	* plot/__contour__.m: __contour__.m: correct order of patches
 
+2008-12-30  Ben Abbott <bpabbott@mac.com>
+
+	* plot/__contour__.m: __contour__.m: correct order of patches
+
 2008-12-24  Doug Stewart  <dastew@sympatico.ca>
 
 	* plot/grid.m: Handle "minor" option.
--- /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 <kai.habel@gmx.de>
+
+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
--- /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 <kai.habel@gmx.de>
+
+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
--- /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} = [@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, diffuse, specular, surface}
+## @end deftypefn
+
+## Author: Kai Habel <kai.habel@gmx.de>
+
+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;
+