File: surfl.m

package info (click to toggle)
octave 9.4.0-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 144,300 kB
  • sloc: cpp: 332,784; ansic: 77,239; fortran: 20,963; objc: 9,396; sh: 8,213; yacc: 4,925; lex: 4,389; perl: 1,544; java: 1,366; awk: 1,259; makefile: 648; xml: 189
file content (232 lines) | stat: -rw-r--r-- 7,083 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
########################################################################
##
## Copyright (C) 2009-2024 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or <https://octave.org/copyright/>.
##
## 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
## <https://www.gnu.org/licenses/>.
##
########################################################################

## -*- texinfo -*-
## @deftypefn  {} {} surfl (@var{z})
## @deftypefnx {} {} surfl (@var{x}, @var{y}, @var{z})
## @deftypefnx {} {} surfl (@dots{}, @var{lsrc})
## @deftypefnx {} {} surfl (@var{x}, @var{y}, @var{z}, @var{lsrc}, @var{P})
## @deftypefnx {} {} surfl (@dots{}, "cdata")
## @deftypefnx {} {} surfl (@dots{}, "light")
## @deftypefnx {} {} surfl (@var{hax}, @dots{})
## @deftypefnx {} {@var{h} =} surfl (@dots{})
## Plot a 3-D surface using shading based on various lighting models.
##
## The surface mesh is plotted using shaded rectangles.  The vertices of the
## rectangles [@var{x}, @var{y}] are typically the output of @code{meshgrid}.
## over a 2-D rectangular region in the x-y plane.  @var{z} determines the
## height above the plane of each vertex.  If only a single @var{z} matrix is
## given, then it is plotted over the meshgrid
## @code{@var{x} = 1:columns (@var{z}), @var{y} = 1:rows (@var{z})}.
## Thus, columns of @var{z} correspond to different @var{x} values and rows
## of @var{z} correspond to different @var{y} values.
##
## The default lighting mode @qcode{"cdata"}, changes the cdata property of the
## surface object to give the impression of a lighted surface.
##
## The alternate mode @qcode{"light"} creates a light object to illuminate the
## surface.
##
## The light source location may be specified using @var{lsrc} which can be
## a 2-element vector [azimuth, elevation] in degrees, or a 3-element vector
## [lx, ly, lz].  The default value is rotated 45 degrees counterclockwise to
## 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 @asis
## @item @qcode{"AM"} strength of ambient light
##
## @item @qcode{"D"} strength of diffuse reflection
##
## @item @qcode{"SP"} strength of specular reflection
##
## @item @qcode{"EXP"} specular exponent
## @end table
##
## If the first argument @var{hax} is an axes handle, then plot into this axes,
## rather than the current axes returned by @code{gca}.
##
## The optional return value @var{h} is a graphics handle to the created
## surface object.
##
## Example:
##
## @example
## @group
## colormap (bone (64));
## surfl (peaks);
## shading interp;
## @end group
## @end example
## @seealso{diffuse, specular, surf, shading, colormap, caxis}
## @end deftypefn

function h = surfl (varargin)

  [hax, varargin, nargin] = __plt_get_axis_arg__ ("surfl", varargin{:});

  if (nargin == 0)
    print_usage ();
  endif

  ## Check for lighting type.
  use_cdata = true;
  if (ischar (varargin{end}))
    switch (lower (varargin{end}))
      case "light"
        use_cdata = false;
      case "cdata"
        use_cdata = true;
      otherwise
        error ("surfl: unknown lighting method");
    endswitch
    varargin(end) = [];
  endif

  ## Check for reflection properties argument.
  ##
  ## r = [ambient light strength,
  ##      diffuse reflection strength,
  ##      specular reflection strength,
  ##      specular shine]
  if (isnumeric (varargin{end}) && isvector (varargin{end})
      && (numel (varargin{end}) == 4))
    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}) && isvector (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

  oldfig = [];
  if (! isempty (hax))
    oldfig = get (0, "currentfigure");
  endif
  unwind_protect
    hax = newplot (hax);

    htmp = surface (varargin{:});
    if (! ishold ())
      set (hax, "view", [-37.5, 30],
                "xgrid", "on", "ygrid", "on", "zgrid", "on");
    endif

    ## Get view vector (vv).
    [az, el] = view ();
    [vv(1), vv(2), vv(3)] = sph2cart ((az - 90) * pi/180.0, el * pi/180.0, 1.0);

    if (! have_lv)
      ## Calculate light vector (lv) from view vector.
      phi = pi / 4;  # 45 degrees
      R = [cos(phi), -sin(phi), 0;
           sin(phi),  cos(phi), 0;
           0,         0,        1];
      lv = (R * vv.').';
    endif

    if (use_cdata)
      set (hax, "clim", [0 1]);

      __update_normals__ (htmp);
      vn = get (htmp, "vertexnormals");
      dar = get (hax, "dataaspectratio");
      vn(:,:,1) *= dar(1);
      vn(:,:,2) *= dar(2);
      vn(:,:,3) *= dar(3);

      ## Normalize vn.
      vn ./= repmat (sqrt (sumsq (vn, 3)), [1, 1, 3]);
      [nr, nc] = size (get (htmp, "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)));
      cdata ./= sum (r(1:3));

      set (htmp, "cdata", cdata);
    else
      light (hax, "position", lv);
      set (htmp, "ambientstrength", r(1), "diffusestrength", r(2), ...
                 "specularstrength", r(3), "specularexponent", r(4));
    endif

  unwind_protect_cleanup
    if (! isempty (oldfig))
      set (0, "currentfigure", oldfig);
    endif
  end_unwind_protect

  if (nargout > 0)
    h = htmp;
  endif

endfunction


%!demo
%! clf;
%! [X,Y,Z] = sombrero ();
%! colormap (copper (64));
%! surfl (X,Y,Z);
%! shading interp;
%! title ("surfl() with defaults");

%!demo
%! clf;
%! [X,Y,Z] = sombrero ();
%! colormap (copper (64));
%! surfl (X,Y,Z, [62.50,30], [0.2 0.6 0.4 25]);
%! shading interp;
%! title ("surfl() with lighting vector and material properties");

%!demo
%! clf;
%! [X, Y] = meshgrid (-3:1/8:3);
%! Z = peaks (X, Y);
%! surfl (X, Y, Z, "light");
%! shading interp;
%! title ("surfl() with light object");