1########################################################################
2##
3## Copyright (C) 2007-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn  {} {} slice (@var{x}, @var{y}, @var{z}, @var{v}, @var{sx}, @var{sy}, @var{sz})
28## @deftypefnx {} {} slice (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi})
29## @deftypefnx {} {} slice (@var{v}, @var{sx}, @var{sy}, @var{sz})
30## @deftypefnx {} {} slice (@var{v}, @var{xi}, @var{yi}, @var{zi})
31## @deftypefnx {} {} slice (@dots{}, @var{method})
32## @deftypefnx {} {} slice (@var{hax}, @dots{})
33## @deftypefnx {} {@var{h} =} slice (@dots{})
34## Plot slices of 3-D data/scalar fields.
35##
36## Each element of the 3-dimensional array @var{v} represents a scalar value at
37## a location given by the parameters @var{x}, @var{y}, and @var{z}.  The
38## parameters @var{x}, @var{x}, and @var{z} are either 3-dimensional arrays of
39## the same size as the array @var{v} in the @qcode{"meshgrid"} format or
40## vectors.  The parameters @var{xi}, etc.@: respect a similar format to
41## @var{x}, etc., and they represent the points at which the array @var{vi}
42## is interpolated using interp3.  The vectors @var{sx}, @var{sy}, and
43## @var{sz} contain points of orthogonal slices of the respective axes.
44##
45## If @var{x}, @var{y}, @var{z} are omitted, they are assumed to be
46## @code{x = 1:size (@var{v}, 2)}, @code{y = 1:size (@var{v}, 1)} and
47## @code{z = 1:size (@var{v}, 3)}.
48##
49## @var{method} is one of:
50##
51## @table @asis
52## @item @qcode{"nearest"}
53## Return the nearest neighbor.
54##
55## @item @qcode{"linear"}
56## Linear interpolation from nearest neighbors.
57##
58## @item @qcode{"cubic"}
59## Cubic interpolation from four nearest neighbors (not implemented yet).
60##
61## @item @qcode{"spline"}
62## Cubic spline interpolation---smooth first and second derivatives
63## throughout the curve.
64## @end table
65##
66## The default method is @qcode{"linear"}.
67##
68## If the first argument @var{hax} is an axes handle, then plot into this axes,
69## rather than the current axes returned by @code{gca}.
70##
71## The optional return value @var{h} is a graphics handle to the created
72## surface object.
73##
74## Examples:
75##
76## @example
77## @group
78## [x, y, z] = meshgrid (linspace (-8, 8, 32));
79## v = sin (sqrt (x.^2 + y.^2 + z.^2)) ./ (sqrt (x.^2 + y.^2 + z.^2));
80## slice (x, y, z, v, [], 0, []);
81##
82## [xi, yi] = meshgrid (linspace (-7, 7));
83## zi = xi + yi;
84## slice (x, y, z, v, xi, yi, zi);
85## @end group
86## @end example
87## @seealso{interp3, surface, pcolor}
88## @end deftypefn
89
90function h = slice (varargin)
91
92  [hax, varargin, nargs] = __plt_get_axis_arg__ ("slice", varargin{:});
93
94  method = "linear";
95
96  if (ischar (varargin{end}))
97    method = varargin{end};
98    nargs -= 1;
99  endif
100
101  if (nargs == 4)
102    v = varargin{1};
103    if (ndims (v) != 3)
104      error ("slice: V must be a 3-dimensional array of values");
105    endif
106    [nx, ny, nz] = size (v);
107    [x, y, z] = meshgrid (1:nx, 1:ny, 1:nz);
108    sx = varargin{2};
109    sy = varargin{3};
110    sz = varargin{4};
111  elseif (nargs == 7)
112    v = varargin{4};
113    if (ndims (v) != 3)
114      error ("slice: V must be a 3-dimensional array of values");
115    endif
116    x = varargin{1};
117    y = varargin{2};
118    z = varargin{3};
119    if (isvector (x) && isvector (y) && isvector (z))
120      [x, y, z] = meshgrid (x, y, z);
121    elseif (ndims (x) == 3 && size_equal (x, y, z))
122      ## Do nothing.
123    else
124      error ("slice: X, Y, Z size mismatch");
125    endif
126    sx = varargin{5};
127    sy = varargin{6};
128    sz = varargin{7};
129  else
130    print_usage ();
131  endif
132
133  if (any ([isvector(sx), isvector(sy), isvector(sz)]))
134    have_sval = true;
135  elseif (ndims (sx) == 2 && size_equal (sx, sy, sz))
136    have_sval = false;
137  else
138    error ("slice: dimensional mismatch for (XI, YI, ZI) or (SX, SY, SZ)");
139  endif
140
141  oldfig = [];
142  if (! isempty (hax))
143    oldfig = get (0, "currentfigure");
144  endif
145  unwind_protect
146    hax = newplot (hax);
147
148    sidx = 1;
149    minv = min (v(:));
150    maxv = max (v(:));
151    set (hax, "clim", double ([minv, maxv]));
152
153    if (have_sval)
154      ns = length (sx) + length (sy) + length (sz);
155      hs = zeros (ns,1);
156      [ny, nx, nz] = size (v);
157      if (length (sz) > 0)
158        for i = 1:length (sz)
159          [xi, yi, zi] = meshgrid (squeeze (x(1,:,1)),
160                                   squeeze (y(:,1,1)), sz(i));
161          vz = squeeze (interp3 (x, y, z, v, xi, yi, zi, method));
162          htmp(sidx++) = surface (xi, yi, sz(i) * ones (size (yi)), vz);
163        endfor
164      endif
165
166      if (length (sy) > 0)
167        for i = length (sy):-1:1
168          [xi, yi, zi] = meshgrid (squeeze (x(1,:,1)),
169                                   sy(i),
170                                   squeeze (z(1,1,:)));
171          vy = squeeze (interp3 (x, y, z, v, xi, yi, zi, method));
172          htmp(sidx++) = surface (squeeze (xi),
173                                  squeeze (sy(i) * ones (size (zi))),
174                                  squeeze (zi), vy);
175        endfor
176      endif
177
178      if (length (sx) > 0)
179        for i = length (sx):-1:1
180          [xi, yi, zi] = meshgrid (sx(i), squeeze (y(:,1,1)), squeeze (z(1,1,:)));
181          vx = squeeze (interp3 (x, y, z, v, xi, yi, zi, method));
182          htmp(sidx++) = surface (squeeze (sx(i) * ones (size (zi))),
183                                  squeeze (yi), squeeze(zi), vx);
184        endfor
185      endif
186    else
187      vi = interp3 (x, y, z, v, sx, sy, sz);
188      htmp = surface (sx, sy, sz, vi);
189    endif
190
191    if (! ishold ())
192      set (hax, "view", [-37.5, 30.0],
193                "xgrid", "on", "ygrid", "on", "zgrid", "on");
194    endif
195
196  unwind_protect_cleanup
197    if (! isempty (oldfig))
198      set (0, "currentfigure", oldfig);
199    endif
200  end_unwind_protect
201
202  if (nargout > 0)
203    h = htmp;
204  endif
205
206endfunction
207
208
209%!demo
210%! clf;
211%! colormap ("default");
212%! [x, y, z] = meshgrid (linspace (-8, 8, 32));
213%! v = sin (sqrt (x.^2 + y.^2 + z.^2)) ./ (sqrt (x.^2 + y.^2 + z.^2));
214%! slice (x, y, z, v, [], 0, []);
215%! title ("slice() demo #1");
216
217%!demo
218%! clf;
219%! colormap ("default");
220%! [x, y, z] = meshgrid (linspace (-8, 8, 32));
221%! v = sin (sqrt (x.^2 + y.^2 + z.^2)) ./ (sqrt (x.^2 + y.^2 + z.^2));
222%! [xi, yi] = meshgrid (linspace (-7, 7));
223%! zi = xi + yi;
224%! slice (x, y, z, v, xi, yi, zi);
225%! title ("slice() demo #2");
226