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 {} {} surfnorm (@var{x}, @var{y}, @var{z}) 28## @deftypefnx {} {} surfnorm (@var{z}) 29## @deftypefnx {} {} surfnorm (@dots{}, @var{prop}, @var{val}, @dots{}) 30## @deftypefnx {} {} surfnorm (@var{hax}, @dots{}) 31## @deftypefnx {} {[@var{nx}, @var{ny}, @var{nz}] =} surfnorm (@dots{}) 32## Find the vectors normal to a meshgridded surface. 33## 34## If @var{x} and @var{y} are vectors, then a typical vertex is 35## (@var{x}(j), @var{y}(i), @var{z}(i,j)). Thus, columns of @var{z} correspond 36## to different @var{x} values and rows of @var{z} correspond to different 37## @var{y} values. If only a single input @var{z} is given then @var{x} is 38## taken to be @code{1:columns (@var{z})} and @var{y} is 39## @code{1:rows (@var{z})}. 40## 41## If no return arguments are requested, a surface plot with the normal 42## vectors to the surface is plotted. 43## 44## Any property/value input pairs are assigned to the surface object. The full 45## list of properties is documented at @ref{Surface Properties}. 46## 47## If the first argument @var{hax} is an axes handle, then plot into this axes, 48## rather than the current axes returned by @code{gca}. 49## 50## If output arguments are requested then the components of the normal 51## vectors are returned in @var{nx}, @var{ny}, and @var{nz} and no plot is 52## made. The normal vectors are unnormalized (magnitude != 1). To normalize, 53## use 54## 55## @example 56## @group 57## len = sqrt (nx.^2 + ny.^2 + nz.^2); 58## nx ./= len; ny ./= len; nz ./= len; 59## @end group 60## @end example 61## 62## An example of the use of @code{surfnorm} is 63## 64## @example 65## surfnorm (peaks (25)); 66## @end example 67## 68## Algorithm: The normal vectors are calculated by taking the cross product 69## of the diagonals of each of the quadrilateral faces in the meshgrid to find 70## the normal vectors at the center of each face. Next, for each meshgrid 71## point the four nearest normal vectors are averaged to obtain the final 72## normal to the surface at the meshgrid point. 73## 74## For surface objects, the @qcode{"VertexNormals"} property contains 75## equivalent information, except possibly near the boundary of the surface 76## where different interpolation schemes may yield slightly different values. 77## 78## @seealso{isonormals, quiver3, surf, meshgrid} 79## @end deftypefn 80 81function [Nx, Ny, Nz] = surfnorm (varargin) 82 83 [hax, varargin, nargin] = __plt_get_axis_arg__ ("surfnorm", varargin{:}); 84 85 if (nargin == 0 || nargin == 2) 86 print_usage (); 87 endif 88 89 if (nargin == 1) 90 z = varargin{1}; 91 [x, y] = meshgrid (1:columns (z), 1:rows (z)); 92 ioff = 2; 93 else 94 x = varargin{1}; 95 y = varargin{2}; 96 z = varargin{3}; 97 ioff = 4; 98 endif 99 100 if (iscomplex (z) || iscomplex (x) || iscomplex (y)) 101 error ("surfnorm: X, Y, and Z must be 2-D real matrices"); 102 endif 103 if (! size_equal (x, y, z)) 104 error ("surfnorm: X, Y, and Z must have the same dimensions"); 105 endif 106 107 ## FIXME: Matlab uses a bicubic interpolation, not linear, along the boundary. 108 ## Do a linear extrapolation for mesh points on the boundary so that the mesh 109 ## is increased by 1 on each side. This allows each original meshgrid point 110 ## to be surrounded by four quadrilaterals and the same calculation can be 111 ## used for interior and boundary points. The extrapolation works badly for 112 ## closed surfaces like spheres. 113 xx = [2 * x(:,1) - x(:,2), x, 2 * x(:,end) - x(:,end-1)]; 114 xx = [2 * xx(1,:) - xx(2,:); xx; 2 * xx(end,:) - xx(end-1,:)]; 115 yy = [2 * y(:,1) - y(:,2), y, 2 * y(:,end) - y(:,end-1)]; 116 yy = [2 * yy(1,:) - yy(2,:); yy; 2 * yy(end,:) - yy(end-1,:)]; 117 zz = [2 * z(:,1) - z(:,2), z, 2 * z(:,end) - z(:,end-1)]; 118 zz = [2 * zz(1,:) - zz(2,:); zz; 2 * zz(end,:) - zz(end-1,:)]; 119 120 u.x = xx(1:end-1,1:end-1) - xx(2:end,2:end); 121 u.y = yy(1:end-1,1:end-1) - yy(2:end,2:end); 122 u.z = zz(1:end-1,1:end-1) - zz(2:end,2:end); 123 v.x = xx(1:end-1,2:end) - xx(2:end,1:end-1); 124 v.y = yy(1:end-1,2:end) - yy(2:end,1:end-1); 125 v.z = zz(1:end-1,2:end) - zz(2:end,1:end-1); 126 127 c = cross ([u.x(:), u.y(:), u.z(:)], [v.x(:), v.y(:), v.z(:)]); 128 w.x = reshape (c(:,1), size (u.x)); 129 w.y = reshape (c(:,2), size (u.y)); 130 w.z = reshape (c(:,3), size (u.z)); 131 132 ## Create normal vectors as mesh vectices from normals at mesh centers 133 nx = (w.x(1:end-1,1:end-1) + w.x(1:end-1,2:end) + 134 w.x(2:end,1:end-1) + w.x(2:end,2:end)) / 4; 135 ny = (w.y(1:end-1,1:end-1) + w.y(1:end-1,2:end) + 136 w.y(2:end,1:end-1) + w.y(2:end,2:end)) / 4; 137 nz = (w.z(1:end-1,1:end-1) + w.z(1:end-1,2:end) + 138 w.z(2:end,1:end-1) + w.z(2:end,2:end)) / 4; 139 140 if (nargout == 0) 141 oldfig = []; 142 if (! isempty (hax)) 143 oldfig = get (0, "currentfigure"); 144 endif 145 unwind_protect 146 hax = newplot (hax); 147 148 surf (x, y, z, varargin{ioff:end}); 149 old_hold_state = get (hax, "nextplot"); 150 unwind_protect 151 set (hax, "nextplot", "add"); 152 153 ## Normalize the normal vectors 154 nmag = sqrt (nx.^2 + ny.^2 + nz.^2); 155 156 ## And correct for the aspect ratio of the display 157 daratio = daspect (hax); 158 damag = sqrt (sumsq (daratio)); 159 160 ## FIXME: May also want to normalize the vectors relative to the size 161 ## of the diagonal. 162 163 nx ./= nmag / (daratio(1)^2 / damag); 164 ny ./= nmag / (daratio(2)^2 / damag); 165 nz ./= nmag / (daratio(3)^2 / damag); 166 167 plot3 ([x(:).'; x(:).' + nx(:).' ; NaN(size(x(:).'))](:), 168 [y(:).'; y(:).' + ny(:).' ; NaN(size(y(:).'))](:), 169 [z(:).'; z(:).' + nz(:).' ; NaN(size(z(:).'))](:), 170 "r"); 171 unwind_protect_cleanup 172 set (hax, "nextplot", old_hold_state); 173 end_unwind_protect 174 175 unwind_protect_cleanup 176 if (! isempty (oldfig)) 177 set (0, "currentfigure", oldfig); 178 endif 179 end_unwind_protect 180 else 181 Nx = nx; 182 Ny = ny; 183 Nz = nz; 184 endif 185 186endfunction 187 188 189%!demo 190%! clf; 191%! colormap ("default"); 192%! surfnorm (peaks (19)); 193%! shading faceted; 194%! title ({"surfnorm() shows surface and normals at each vertex", ... 195%! "peaks() function with 19 faces"}); 196 197%!demo 198%! clf; 199%! colormap ("default"); 200%! [x, y, z] = sombrero (10); 201%! surfnorm (x, y, z); 202%! title ({"surfnorm() shows surface and normals at each vertex", ... 203%! "sombrero() function with 10 faces"}); 204 205## Test input validation 206%!error surfnorm () 207%!error surfnorm (1,2) 208%!error <X, Y, and Z must be 2-D real matrices> surfnorm (i) 209%!error <X, Y, and Z must be 2-D real matrices> surfnorm (i, 1, 1) 210%!error <X, Y, and Z must be 2-D real matrices> surfnorm (1, i, 1) 211%!error <X, Y, and Z must be 2-D real matrices> surfnorm (1, 1, i) 212%!error <X, Y, and Z must have the same dimensions> surfnorm ([1 2], 1, 1) 213%!error <X, Y, and Z must have the same dimensions> surfnorm (1, [1 2], 1) 214%!error <X, Y, and Z must have the same dimensions> surfnorm (1, 1, [1 2]) 215