# HG changeset patch # User Kai Habel # Date 1231426310 -3600 # Node ID 6c4e2fc14beb66cc57a92a04cb7fc9650b2cd7bc # Parent 4481888537221cf8605db1873d0f4b35ee0aca82 Add new 3D plotting function surfl. Add light functions diffuse and specular diff -r 448188853722 -r 6c4e2fc14beb scripts/ChangeLog --- 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 + + * plot/surfl.m: New function + * plot/diffuse.m: dito + * plot/specular.m: dito + 2008-12-28 Jaroslav Hajek * 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 + + * plot/__contour__.m: __contour__.m: correct order of patches + 2008-12-24 Doug Stewart * plot/grid.m: Handle "minor" option. diff -r 448188853722 -r 6c4e2fc14beb 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 +## . + +## -*- 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 + +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 448188853722 -r 6c4e2fc14beb 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 +## . + +## -*- 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 + +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 448188853722 -r 6c4e2fc14beb 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 +## . + +## -*- 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 + +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; +