Mercurial > octave-libgccjit
comparison scripts/plot/draw/surfnorm.m @ 18947:91f626902d17
surfnorm.m: Overhaul for Matlab compatibility.
* surfnorm.m: Redo docstring. Use better input validation and accept optional
property/value pairs for surface object. Fix incorrect meshgrid call reversing
x and y. Color the surface normals red for compatibility with Matlab. Use
in-place division operation for better performance. Add tests for input
validation.
author | Rik <rik@octave.org> |
---|---|
date | Sun, 20 Jul 2014 13:26:10 -0700 |
parents | d63878346099 |
children |
comparison
equal
deleted
inserted
replaced
18943:714ce8ca71ea | 18947:91f626902d17 |
---|---|
17 ## <http://www.gnu.org/licenses/>. | 17 ## <http://www.gnu.org/licenses/>. |
18 | 18 |
19 ## -*- texinfo -*- | 19 ## -*- texinfo -*- |
20 ## @deftypefn {Function File} {} surfnorm (@var{x}, @var{y}, @var{z}) | 20 ## @deftypefn {Function File} {} surfnorm (@var{x}, @var{y}, @var{z}) |
21 ## @deftypefnx {Function File} {} surfnorm (@var{z}) | 21 ## @deftypefnx {Function File} {} surfnorm (@var{z}) |
22 ## @deftypefnx {Function File} {} surfnorm (@dots{}, @var{prop}, @var{val}, @dots{}) | |
23 ## @deftypefnx {Function File} {} surfnorm (@var{hax}, @dots{}) | |
22 ## @deftypefnx {Function File} {[@var{nx}, @var{ny}, @var{nz}] =} surfnorm (@dots{}) | 24 ## @deftypefnx {Function File} {[@var{nx}, @var{ny}, @var{nz}] =} surfnorm (@dots{}) |
23 ## @deftypefnx {Function File} {} surfnorm (@var{h}, @dots{}) | 25 ## Find the vectors normal to a meshgridded surface. |
24 ## Find the vectors normal to a meshgridded surface. The meshed gridded | 26 ## |
25 ## surface is defined by @var{x}, @var{y}, and @var{z}. If @var{x} and | 27 ## If @var{x} and @var{y} are vectors, then a typical vertex is |
26 ## @var{y} are not defined, then it is assumed that they are given by | 28 ## (@var{x}(j), @var{y}(i), @var{z}(i,j)). Thus, columns of @var{z} correspond |
29 ## to different @var{x} values and rows of @var{z} correspond to different | |
30 ## @var{y} values. If only a single input @var{z} is given then @var{x} is | |
31 ## taken to be @code{1:rows (@var{z})} and @var{y} is | |
32 ## @code{1:columns (@var{z})}. | |
33 ## | |
34 ## If no return arguments are requested, a surface plot with the normal | |
35 ## vectors to the surface is plotted. | |
36 ## | |
37 ## Any property/value input pairs are assigned to the surface object. | |
38 ## | |
39 ## If the first argument @var{hax} is an axes handle, then plot into this axis, | |
40 ## rather than the current axes returned by @code{gca}. | |
41 ## | |
42 ## If output arguments are requested then the components of the normal | |
43 ## vectors are returne in @var{nx}, @var{ny}, and @var{nz} and no plot is | |
44 ## made. | |
45 ## | |
46 ## An example of the use of @code{surfnorm} is | |
27 ## | 47 ## |
28 ## @example | 48 ## @example |
29 ## @group | 49 ## surfnorm (peaks (25)); |
30 ## [@var{x}, @var{y}] = meshgrid (1:rows (@var{z}), | |
31 ## 1:columns (@var{z})); | |
32 ## @end group | |
33 ## @end example | 50 ## @end example |
34 ## | 51 ## |
35 ## If no return arguments are requested, a surface plot with the normal | 52 ## Algorithm: The normal vectors are calculated by taking the cross product |
36 ## vectors to the surface is plotted. Otherwise the components of the normal | 53 ## of the diagonals of each of the quadrilaterals in the meshgrid to find the |
37 ## vectors at the mesh gridded points are returned in @var{nx}, @var{ny}, | |
38 ## and @var{nz}. | |
39 ## | |
40 ## The normal vectors are calculated by taking the cross product of the | |
41 ## diagonals of each of the quadrilaterals in the meshgrid to find the | |
42 ## normal vectors of the centers of these quadrilaterals. The four nearest | 54 ## normal vectors of the centers of these quadrilaterals. The four nearest |
43 ## normal vectors to the meshgrid points are then averaged to obtain the | 55 ## normal vectors to the meshgrid points are then averaged to obtain the |
44 ## normal to the surface at the meshgridded points. | 56 ## normal to the surface at the meshgridded points. |
45 ## | 57 ## |
46 ## An example of the use of @code{surfnorm} is | 58 ## @seealso{isonormals, quiver3, surf, meshgrid} |
47 ## | |
48 ## @example | |
49 ## surfnorm (peaks (25)); | |
50 ## @end example | |
51 ## @seealso{surf, quiver3} | |
52 ## @end deftypefn | 59 ## @end deftypefn |
53 | 60 |
54 function [Nx, Ny, Nz] = surfnorm (varargin) | 61 function [Nx, Ny, Nz] = surfnorm (varargin) |
55 | 62 |
56 [hax, varargin, nargin] = __plt_get_axis_arg__ ("surfnorm", varargin{:}); | 63 [hax, varargin, nargin] = __plt_get_axis_arg__ ("surfnorm", varargin{:}); |
57 | 64 |
58 if (nargin != 1 && nargin != 3) | 65 if (nargin == 0 || nargin == 2) |
59 print_usage (); | 66 print_usage (); |
60 endif | 67 endif |
61 | 68 |
62 if (nargin == 1) | 69 if (nargin == 1) |
63 z = varargin{1}; | 70 z = varargin{1}; |
64 [x, y] = meshgrid (1:rows (z), 1:columns (z)); | 71 [x, y] = meshgrid (1:columns (z), 1:rows (z)); |
65 ioff = 2; | 72 ioff = 2; |
66 else | 73 else |
67 x = varargin{1}; | 74 x = varargin{1}; |
68 y = varargin{2}; | 75 y = varargin{2}; |
69 z = varargin{3}; | 76 z = varargin{3}; |
70 ioff = 4; | 77 ioff = 4; |
71 endif | 78 endif |
72 | 79 |
73 if (!ismatrix (z) || isvector (z) || isscalar (z)) | 80 if (iscomplex (z) || iscomplex (x) || iscomplex (y)) |
74 error ("surfnorm: Z argument must be a matrix"); | 81 error ("surfnorm: X, Y, and Z must be 2-D real matrices"); |
75 endif | 82 endif |
76 if (! size_equal (x, y, z)) | 83 if (! size_equal (x, y, z)) |
77 error ("surfnorm: X, Y, and Z must have the same dimensions"); | 84 error ("surfnorm: X, Y, and Z must have the same dimensions"); |
78 endif | 85 endif |
79 | 86 |
80 ## Make life easier, and avoid having to do the extrapolation later, do | 87 ## Do a linear extrapolation for mesh points on the boundary so that the mesh |
81 ## a simpler linear extrapolation here. This is approximative, and works | 88 ## is increased by 1 on each side. This allows each original meshgrid point |
82 ## badly for closed surfaces like spheres. | 89 ## to be surrounded by four quadrilaterals and the same calculation can be |
83 xx = [2 .* x(:,1) - x(:,2), x, 2 .* x(:,end) - x(:,end-1)]; | 90 ## used for interior and boundary points. The extrapolation works badly for |
84 xx = [2 .* xx(1,:) - xx(2,:); xx; 2 .* xx(end,:) - xx(end-1,:)]; | 91 ## closed surfaces like spheres. |
85 yy = [2 .* y(:,1) - y(:,2), y, 2 .* y(:,end) - y(:,end-1)]; | 92 xx = [2 * x(:,1) - x(:,2), x, 2 * x(:,end) - x(:,end-1)]; |
86 yy = [2 .* yy(1,:) - yy(2,:); yy; 2 .* yy(end,:) - yy(end-1,:)]; | 93 xx = [2 * xx(1,:) - xx(2,:); xx; 2 * xx(end,:) - xx(end-1,:)]; |
87 zz = [2 .* z(:,1) - z(:,2), z, 2 .* z(:,end) - z(:,end-1)]; | 94 yy = [2 * y(:,1) - y(:,2), y, 2 * y(:,end) - y(:,end-1)]; |
88 zz = [2 .* zz(1,:) - zz(2,:); zz; 2 .* zz(end,:) - zz(end-1,:)]; | 95 yy = [2 * yy(1,:) - yy(2,:); yy; 2 * yy(end,:) - yy(end-1,:)]; |
96 zz = [2 * z(:,1) - z(:,2), z, 2 * z(:,end) - z(:,end-1)]; | |
97 zz = [2 * zz(1,:) - zz(2,:); zz; 2 * zz(end,:) - zz(end-1,:)]; | |
89 | 98 |
90 u.x = xx(1:end-1,1:end-1) - xx(2:end,2:end); | 99 u.x = xx(1:end-1,1:end-1) - xx(2:end,2:end); |
91 u.y = yy(1:end-1,1:end-1) - yy(2:end,2:end); | 100 u.y = yy(1:end-1,1:end-1) - yy(2:end,2:end); |
92 u.z = zz(1:end-1,1:end-1) - zz(2:end,2:end); | 101 u.z = zz(1:end-1,1:end-1) - zz(2:end,2:end); |
93 v.x = xx(1:end-1,2:end) - xx(2:end,1:end-1); | 102 v.x = xx(1:end-1,2:end) - xx(2:end,1:end-1); |
99 w.y = reshape (c(:,2), size (u.y)); | 108 w.y = reshape (c(:,2), size (u.y)); |
100 w.z = reshape (c(:,3), size (u.z)); | 109 w.z = reshape (c(:,3), size (u.z)); |
101 | 110 |
102 ## Create normal vectors as mesh vectices from normals at mesh centers | 111 ## Create normal vectors as mesh vectices from normals at mesh centers |
103 nx = (w.x(1:end-1,1:end-1) + w.x(1:end-1,2:end) + | 112 nx = (w.x(1:end-1,1:end-1) + w.x(1:end-1,2:end) + |
104 w.x(2:end,1:end-1) + w.x(2:end,2:end)) ./ 4; | 113 w.x(2:end,1:end-1) + w.x(2:end,2:end)) / 4; |
105 ny = (w.y(1:end-1,1:end-1) + w.y(1:end-1,2:end) + | 114 ny = (w.y(1:end-1,1:end-1) + w.y(1:end-1,2:end) + |
106 w.y(2:end,1:end-1) + w.y(2:end,2:end)) ./ 4; | 115 w.y(2:end,1:end-1) + w.y(2:end,2:end)) / 4; |
107 nz = (w.z(1:end-1,1:end-1) + w.z(1:end-1,2:end) + | 116 nz = (w.z(1:end-1,1:end-1) + w.z(1:end-1,2:end) + |
108 w.z(2:end,1:end-1) + w.z(2:end,2:end)) ./ 4; | 117 w.z(2:end,1:end-1) + w.z(2:end,2:end)) / 4; |
109 | 118 |
119 ## FIXME: According to Matlab documentation the vertex normals | |
120 ## returned are not normalized. | |
110 ## Normalize the normal vectors | 121 ## Normalize the normal vectors |
111 len = sqrt (nx.^2 + ny.^2 + nz.^2); | 122 len = sqrt (nx.^2 + ny.^2 + nz.^2); |
112 nx = nx ./ len; | 123 nx ./= len; |
113 ny = ny ./ len; | 124 ny ./= len; |
114 nz = nz ./ len; | 125 nz ./= len; |
115 | 126 |
116 if (nargout == 0) | 127 if (nargout == 0) |
117 oldfig = []; | 128 oldfig = []; |
118 if (! isempty (hax)) | 129 if (! isempty (hax)) |
119 oldfig = get (0, "currentfigure"); | 130 oldfig = get (0, "currentfigure"); |
123 | 134 |
124 surf (x, y, z, varargin{ioff:end}); | 135 surf (x, y, z, varargin{ioff:end}); |
125 old_hold_state = get (hax, "nextplot"); | 136 old_hold_state = get (hax, "nextplot"); |
126 unwind_protect | 137 unwind_protect |
127 set (hax, "nextplot", "add"); | 138 set (hax, "nextplot", "add"); |
128 plot3 ([x(:)'; x(:).' + nx(:).' ; NaN(size(x(:).'))](:), | 139 |
129 [y(:)'; y(:).' + ny(:).' ; NaN(size(y(:).'))](:), | 140 ## FIXME: Scale unit normals by data aspect ratio in order for |
130 [z(:)'; z(:).' + nz(:).' ; NaN(size(z(:).'))](:), | 141 ## normals to appear correct. |
131 varargin{ioff:end}); | 142 ##daratio = daspect (hax); |
143 ##daspect ("manual"); | |
144 ##len = norm (daratio); | |
145 ## This assumes an even meshgrid which isn't a great assumption | |
146 ##dx = x(1,2) - x(1,1); | |
147 ##dy = y(2,1) - y(1,1); | |
148 ##nx *= daratio(1); | |
149 ##ny *= daratio(2); | |
150 ##nz *= daratio(3); | |
151 ##len = sqrt (nx.^2 + ny.^2 + nz.^2); | |
152 ##nx ./= len; | |
153 ##ny ./= len; | |
154 ##nz ./= len; | |
155 plot3 ([x(:).'; x(:).' + nx(:).' ; NaN(size(x(:).'))](:), | |
156 [y(:).'; y(:).' + ny(:).' ; NaN(size(y(:).'))](:), | |
157 [z(:).'; z(:).' + nz(:).' ; NaN(size(z(:).'))](:), | |
158 "r"); | |
132 unwind_protect_cleanup | 159 unwind_protect_cleanup |
133 set (hax, "nextplot", old_hold_state); | 160 set (hax, "nextplot", old_hold_state); |
134 end_unwind_protect | 161 end_unwind_protect |
135 | 162 |
136 unwind_protect_cleanup | 163 unwind_protect_cleanup |
148 | 175 |
149 | 176 |
150 %!demo | 177 %!demo |
151 %! clf; | 178 %! clf; |
152 %! colormap ('default'); | 179 %! colormap ('default'); |
153 %! [x, y, z] = peaks (10); | 180 %! surfnorm (peaks (32)); |
154 %! surfnorm (x, y, z); | 181 %! shading interp; |
182 %! title ({'surfnorm() shows surface and normals at each vertex', ... | |
183 %! 'peaks() function with 32 faces'}); | |
155 | 184 |
156 %!demo | 185 %!demo |
157 %! clf; | 186 %! clf; |
158 %! colormap ('default'); | 187 %! colormap ('default'); |
159 %! surfnorm (peaks (10)); | 188 %! [x, y, z] = sombrero (10); |
160 | 189 %! surfnorm (x, y, z); |
161 %!demo | 190 |
162 %! clf; | 191 %% Test input validation |
163 %! colormap ('default'); | 192 %!error surfnorm () |
164 %! surfnorm (peaks (32)); | 193 %!error surfnorm (1,2) |
165 %! shading interp; | 194 %!error <X, Y, and Z must be 2-D real matrices> surfnorm (i) |
166 | 195 %!error <X, Y, and Z must be 2-D real matrices> surfnorm (i, 1, 1) |
196 %!error <X, Y, and Z must be 2-D real matrices> surfnorm (1, i, 1) | |
197 %!error <X, Y, and Z must be 2-D real matrices> surfnorm (1, 1, i) | |
198 %!error <X, Y, and Z must have the same dimensions> surfnorm ([1 2], 1, 1) | |
199 %!error <X, Y, and Z must have the same dimensions> surfnorm (1, [1 2], 1) | |
200 %!error <X, Y, and Z must have the same dimensions> surfnorm (1, 1, [1 2]) | |
201 |