Mercurial > octave-nkf
changeset 9110:22ae6b3411a7
Add isocolor, isonormals and isosurface functions (For Martin Helm). Add 3D filled triangular patches and the trisurf function
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Sat, 11 Apr 2009 16:26:01 +0200 |
| parents | 978c863bc8e5 |
| children | 96e7a72be5e7 |
| files | ChangeLog NEWS doc/ChangeLog doc/interpreter/contributors.in scripts/ChangeLog scripts/geometry/Makefile.in scripts/geometry/trimesh.m scripts/geometry/trisurf.m scripts/plot/Makefile.in scripts/plot/__go_draw_axes__.m scripts/plot/__interp_cube__.m scripts/plot/__marching_cube__.m scripts/plot/__patch__.m scripts/plot/isocolors.m scripts/plot/isonormals.m scripts/plot/isosurface.m scripts/plot/patch.m |
| diffstat | 17 files changed, 1688 insertions(+), 200 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog Sat Apr 11 16:05:16 2009 +0200 +++ b/ChangeLog Sat Apr 11 16:26:01 2009 +0200 @@ -1,3 +1,7 @@ +2009-04-11 David Bateman <dbateman@free.fr> + + * NEWS: Add new graphics functions. + 2009-04-05 John W. Eaton <jwe@octave.org> * configure.in: Use AC_USE_SYSTEM_EXTENSIONS instead of
--- a/NEWS Sat Apr 11 16:05:16 2009 +0200 +++ b/NEWS Sat Apr 11 16:26:01 2009 +0200 @@ -22,15 +22,16 @@ ** New graphics functions: - addlistener diffuse ezsurfc plotmatrix - addproperty ezcontour findall refresh - allchild ezcontourf gcbf refreshdata - available_backends ezmesh gcbo specular - backend ezmeshc ginput surfl - cla ezplot gtext waitforbuttonpress - clabel ezplot3 intwarning - comet ezpolar ishghandle - dellistener ezsurf linkprop + addlistener ezcontour gcbo refresh + addproperty ezcontourf ginput refreshdata + allchild ezmesh gtext specular + available_backends ezmeshc intwarning surfl + backend ezplot ishghandle trisurf + cla ezplot3 isocolors waitforbuttonpress + clabel ezpolar isonormals + comet ezsurf isosurface + dellistener findall linkprop + diffuse gcbf plotmatrix ** New experimental OpenGL/FLTK based plotting system.
--- a/doc/ChangeLog Sat Apr 11 16:05:16 2009 +0200 +++ b/doc/ChangeLog Sat Apr 11 16:26:01 2009 +0200 @@ -1,3 +1,7 @@ +2009-04-11 David Bateman <dbateman@free.fr> + + * interpreter/contributors.in: Add Martin Helm. + 2009-04-06 John W. Eaton <jwe@octave.org> * texinfo.tex: Prefer PDF image files if generating PDF output.
--- a/doc/interpreter/contributors.in Sat Apr 11 16:05:16 2009 +0200 +++ b/doc/interpreter/contributors.in Sat Apr 11 16:26:01 2009 +0200 @@ -70,6 +70,7 @@ Soren Hauberg Dave Hawthorne Daniel Heiserer +Martin Helm Yozo Hida Ryan Hinton Roman Hodek
--- a/scripts/ChangeLog Sat Apr 11 16:05:16 2009 +0200 +++ b/scripts/ChangeLog Sat Apr 11 16:26:01 2009 +0200 @@ -1,3 +1,19 @@ +2009-04-11 David Bateman <dbateman@free.fr> + + * geometry/trisurf.m: New file. + * geometry/Makefile.in (SOURCES): Add it here. + * geometry/trimesh.m: Convert to using 3D patches. + * plot/__go_draw_axes__.m: Allow 3D filled triangular patches. + * plot/__patch__.m: Rewrite to allow update of dependent variables + with listener functions amongst themselves. + * plot/patch.m: Add 3D demo. Update the documentation. + +2009-04-11 Martin Helm <martinh@sirius.mhelm.de> + + * plot/__interp_cube__.m, plot/__marching_cube__.m, isocolors.m, + isonnormals.m, isosurface.m: New files. + * plot/Makefile.in (SOURCES): Add them here. + 2009-04-11 Jaroslav Hajek <highegg@gmail.com> * set/intersect.m: Add missing branch.
--- a/scripts/geometry/Makefile.in Sat Apr 11 16:05:16 2009 +0200 +++ b/scripts/geometry/Makefile.in Sat Apr 11 16:26:01 2009 +0200 @@ -46,6 +46,7 @@ rectint.m \ trimesh.m \ triplot.m \ + trisurf.m \ tsearchn.m \ voronoi.m \ voronoin.m
--- a/scripts/geometry/trimesh.m Sat Apr 11 16:05:16 2009 +0200 +++ b/scripts/geometry/trimesh.m Sat Apr 11 16:26:01 2009 +0200 @@ -38,27 +38,29 @@ elseif (ischar (z)) triplot (tri, x, y, z, varargin{:}); else - idx = tri(:, [1,2,3,1]).'; - nt = size (tri, 1); - ## FIXME We should really use patch instead of plot3, but we don't - ## have a patch function, and probably won't in 3D that works with - ## gnuplot + newplot (); if (nargout > 0) - h = plot3 ([x(idx); NaN*ones(1, nt)](:), - [y(idx); NaN*ones(1, nt)](:), - [z(idx); NaN*ones(1, nt)](:), varargin{:}); + h = patch ("Vertices", [x(:), y(:), z(:)], "Faces", tri, + "FaceColor", "none", "EdgeColor", __next_line_color__(), + varargin{:}); else - plot3 ([x(idx); NaN*ones(1, nt)](:), - [y(idx); NaN*ones(1, nt)](:), - [z(idx); NaN*ones(1, nt)](:), varargin{:}); + patch ("Vertices", [x(:), y(:), z(:)], "Faces", tri, + "FaceColor", "none", "EdgeColor", __next_line_color__(), + varargin{:}); + endif + + if (! ishold ()) + set (gca(), "view", [-37.5, 30], + "xgrid", "on", "ygrid", "on", "zgrid", "on"); endif endif endfunction %!demo +%! N = 10; %! rand ('state', 10) -%! x = 3 - 6 * rand (1, 50); -%! y = 3 - 6 * rand (1, 50); +%! x = 3 - 6 * rand (N, N); +%! y = 3 - 6 * rand (N, N); %! z = peaks (x, y); -%! tri = delaunay (x, y); -%! trimesh (tri, x, y, z); +%! tri = delaunay (x(:), y(:)); +%! trimesh (tri, x(:), y(:), z(:));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/geometry/trisurf.m Sat Apr 11 16:26:01 2009 +0200 @@ -0,0 +1,75 @@ +## Copyright (C) 2007, 2008 David Bateman +## +## 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 +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} trisurf (@var{tri}, @var{x}, @var{y}, @var{z}) +## @deftypefnx {Function File} {@var{h} =} trisurf (@dots{}) +## Plot a triangular surface in 3D. The variable @var{tri} is the triangular +## meshing of the points @code{(@var{x}, @var{y})} which is returned +## from @code{delaunay}. The variable @var{z} is value at the point +## @code{(@var{x}, @var{y})}. The output argument @var{h} is the graphic +## handle to the plot. +## @seealso{triplot, delaunay3} +## @end deftypefn + +function h = trisurf (tri, x, y, z, varargin) + + if (nargin < 3) + print_usage (); + endif + + if (nargin == 3) + triplot (tri, x, y); + elseif (ischar (z)) + triplot (tri, x, y, z, varargin{:}); + else + if (nargin > 4 && isnumeric (varargin{1})) + c = varargin{1}; + varargin(1) = []; + else + c = z; + endif + + newplot (); + if (nargout > 0) + h = patch ("Faces", tri, "Vertices", [x(:), y(:), z(:)], + "FaceVertexCData", reshape (c, numel (c), 1), + "FaceColor", "flat", "EdgeColor", "none", + varargin{:}); + else + patch ("Faces", tri, "Vertices", [x(:), y(:), z(:)], + "FaceVertexCData", reshape (c, numel (c), 1), + "FaceColor", "flat", "EdgeColor", "none", + varargin{:}); + endif + + if (! ishold ()) + set (gca(), "view", [-37.5, 30], + "xgrid", "on", "ygrid", "on", "zgrid", "on"); + endif + endif +endfunction + +%!demo +%! N = 10; +%! rand ('state', 10) +%! x = 3 - 6 * rand (N, N); +%! y = 3 - 6 * rand (N, N); +%! z = peaks (x, y); +%! tri = delaunay (x(:), y(:)); +%! trisurf (tri, x(:), y(:), z(:));
--- a/scripts/plot/Makefile.in Sat Apr 11 16:05:16 2009 +0200 +++ b/scripts/plot/Makefile.in Sat Apr 11 16:26:01 2009 +0200 @@ -54,7 +54,9 @@ __go_draw_figure__.m \ __gnuplot_ginput__.m \ __gnuplot_version__.m \ + __interp_cube__.m \ __line__.m \ + __marching_cube__.m \ __next_line_color__.m \ __patch__.m \ __plr1__.m \ @@ -131,6 +133,9 @@ isfigure.m \ ishghandle.m \ ishold.m \ + isocolors.m \ + isonormals.m \ + isosurface.m \ legend.m \ line.m \ linkprop.m \
--- a/scripts/plot/__go_draw_axes__.m Sat Apr 11 16:05:16 2009 +0200 +++ b/scripts/plot/__go_draw_axes__.m Sat Apr 11 16:26:01 2009 +0200 @@ -318,6 +318,7 @@ if (! cautoscale && clim(1) == clim(2)) clim(2)++; endif + addedcmap = []; [view_cmd, view_fcn, view_zoom] = image_viewer (); use_gnuplot_for_images = (ischar (view_fcn) @@ -357,6 +358,7 @@ is_image_data(data_idx) = true; parametric(data_idx) = false; have_cdata(data_idx) = false; + have_3d_patch(data_idx) = false; [y_dim, x_dim] = size (img_data(:,:,1)); if (x_dim > 1) @@ -406,6 +408,8 @@ is_image_data(data_idx) = false; parametric(data_idx) = true; have_cdata(data_idx) = false; + have_3d_patch(data_idx) = false; + if (isempty (obj.keylabel)) titlespec{data_idx} = "title \"\""; else @@ -517,6 +521,7 @@ cdat = []; endif + data_3d_idx = NaN; for i = 1:nc xcol = obj.xdata(:,i); ycol = obj.ydata(:,i); @@ -533,22 +538,42 @@ if (strncmp (obj.facecolor, "none", 4)) hidden_removal = false; else + if (isnan (hidden_removal)) hidden_removal = true; endif if (nd == 3) - error ("gnuplot (as of v4.2) only supports 2D filled patches"); + if (numel (xcol) > 3) + error ("gnuplot (as of v4.2) only supports 3D filled triangular patches"); + else + if (isnan (data_3d_idx)) + data_idx++; + data_3d_idx = data_idx; + is_image_data(data_idx) = false; + parametric(data_idx) = false; + have_cdata(data_idx) = true; + have_3d_patch(data_idx) = true; + withclause{data_3d_idx} = sprintf ("with pm3d"); + usingclause{data_3d_idx} = "using 1:2:3:4"; + data{data_3d_idx} = []; + endif + local_idx = data_3d_idx; + ccdat = NaN; + endif + else + data_idx++; + local_idx = data_idx; + is_image_data(data_idx) = false; + parametric(data_idx) = false; + have_cdata(data_idx) = false; + have_3d_patch(data_idx) = false; endif - data_idx++; - is_image_data(data_idx) = false; - parametric(data_idx) = false; - have_cdata(data_idx) = false; if (i > 1 || isempty (obj.keylabel)) - titlespec{data_idx} = "title \"\""; + titlespec{local_idx} = "title \"\""; else tmp = undo_string_escapes (__maybe_munge_text__ (enhanced, obj, "keylabel")); - titlespec{data_idx} = cstrcat ("title \"", tmp, "\""); + titlespec{local_idx} = cstrcat ("title \"", tmp, "\""); endif if (isfield (obj, "facecolor")) if ((strncmp (obj.facecolor, "flat", 4) @@ -572,6 +597,8 @@ if (strncmp (obj.facecolor, "flat", 4)) if (numel(ccol) == 3) color = ccol; + elseif (nd == 3 && numel (xcol) == 3) + ccdat = ccol * ones (3,1); else r = 1 + round ((size (cmap, 1) - 1) * (ccol - clim(1))/(clim(2) - clim(1))); @@ -579,10 +606,22 @@ color = cmap(r, :); endif elseif (strncmp (obj.facecolor, "interp", 6)) - warning ("\"interp\" not supported, using 1st entry of cdata"); - r = 1 + round ((size (cmap, 1) - 1) * ccol(1)); - r = max (1, min (r, size (cmap, 1))); - color = cmap(r,:); + if (nd == 3 && numel (xcol) == 3) + ccdat = ccol; + if (! isvector (ccdat)) + tmp = rows(cmap) + rows(addedcmap) + ... + [1 : rows(ccdat)]; + addedcmap = [addedcmap; ccdat]; + ccdat = tmp(:); + else + ccdat = ccdat(:); + endif + else + warning ("\"interp\" not supported, using 1st entry of cdata"); + r = 1 + round ((size (cmap, 1) - 1) * ccol(1)); + r = max (1, min (r, size (cmap, 1))); + color = cmap(r,:); + endif endif elseif (isnumeric (obj.facecolor)) color = obj.facecolor; @@ -593,21 +632,32 @@ color = [0, 1, 0]; endif - if (mono) - colorspec = ""; - elseif (__gnuplot_has_feature__ ("transparent_patches") - && isscalar (obj.facealpha)) - colorspec = sprintf ("lc rgb \"#%02x%02x%02x\" fillstyle transparent solid %f", - round (255*color), obj.facealpha); + if (nd == 3 && numel (xcol) == 3) + if (isnan (ccdat)) + ccdat = (rows (cmap) + rows(addedcmap) + 1) * ones(3, 1); + addedcmap = [addedcmap; reshape(color, 1, 3)]; + endif + data{data_3d_idx} = [data{data_3d_idx}, ... + [[xcol; xcol(end)], [ycol; ycol(end)], ... + [zcol; zcol(end)], [ccdat; ccdat(end)]]']; else - colorspec = sprintf ("lc rgb \"#%02x%02x%02x\"", - round (255*color)); + if (mono) + colorspec = ""; + elseif (__gnuplot_has_feature__ ("transparent_patches") + && isscalar (obj.facealpha)) + colorspec = sprintf ("lc rgb \"#%02x%02x%02x\" fillstyle transparent solid %f", + round (255*color), obj.facealpha); + else + colorspec = sprintf ("lc rgb \"#%02x%02x%02x\"", + round (255*color)); + endif + + withclause{data_idx} = sprintf ("with filledcurve %s", + colorspec); + data{data_idx} = [xcol, ycol]'; + usingclause{data_idx} = sprintf ("record=%d using ($1):($2)", + numel (xcol)); endif - withclause{data_idx} = sprintf ("with filledcurve %s", - colorspec); - data{data_idx} = [xcol, ycol]'; - usingclause{data_idx} = sprintf ("record=%d using ($1):($2)", - numel (xcol)); endif endif @@ -618,6 +668,7 @@ is_image_data(data_idx) = false; parametric(data_idx) = false; have_cdata(data_idx) = false; + have_3d_patch(data_idx) = false; titlespec{data_idx} = "title \"\""; usingclause{data_idx} = sprintf ("record=%d", numel (obj.xdata)); @@ -793,6 +844,7 @@ is_image_data(data_idx) = false; parametric(data_idx) = false; have_cdata(data_idx) = true; + have_3d_patch(data_idx) = false; [style, typ, with] = do_linestyle_command (obj, data_idx, mono, plot_stream); if (isempty (obj.keylabel)) @@ -1037,12 +1089,22 @@ cmap = parent_figure_obj.colormap; cmap_sz = rows(cmap); - if (! any (isinf (clim))) if (truecolor || ! cdatadirect) - fprintf (plot_stream, "set cbrange [%g:%g];\n", clim); + if (rows(addedcmap) > 0) + for i = 1:data_idx + if (have_3d_patch(i)) + data{i}(end,:) = clim(2) * (data{i}(end, :) - 0.5) / cmap_sz; + endif + endfor + fprintf (plot_stream, "set cbrange [%g:%g];\n", clim(1), clim(2) * + (cmap_sz + rows(addedcmap)) / cmap_sz); + else + fprintf (plot_stream, "set cbrange [%g:%g];\n", clim); + endif else - fprintf (plot_stream, "set cbrange [1:%d];\n", cmap_sz); + fprintf (plot_stream, "set cbrange [1:%d];\n", cmap_sz + + rows (addedcmap)); endif endif @@ -1160,6 +1222,8 @@ endfor endif + cmap = [cmap; addedcmap]; + cmap_sz = cmap_sz + rows(addedcmap); if (length(cmap) > 0) fprintf (plot_stream, "set palette positive color model RGB maxcolors %i;\n", @@ -1190,7 +1254,11 @@ fprintf (plot_stream, "set view %.15g, %.15g;\n", rot_x, rot_z); endif endif - if (is_image_data (1)) + if (have_3d_patch (1)) + fputs (plot_stream, "set pm3d depthorder\n"); + fprintf (plot_stream, "%s \"-\" %s %s %s \\\n", plot_cmd, + usingclause{1}, titlespec{1}, withclause{1}); + elseif (is_image_data (1)) fprintf (plot_stream, "%s \"-\" %s %s %s \\\n", plot_cmd, usingclause{1}, titlespec{1}, withclause{1}); else @@ -1198,8 +1266,11 @@ usingclause{1}, titlespec{1}, withclause{1}); endif for i = 2:data_idx - if (is_image_data (i)) - fprintf (plot_stream, "%s \"-\" %s %s %s \\\n", plot_cmd, + if (have_3d_patch (i)) + fprintf (plot_stream, ", \"-\" %s %s %s \\\n", + usingclause{i}, titlespec{i}, withclause{i}); + elseif (is_image_data (i)) + fprintf (plot_stream, "%s \"-\" %s %s %s \\\n", plot_cmd, usingclause{i}, titlespec{i}, withclause{i}); else fprintf (plot_stream, ", \"-\" binary format='%%float64' %s %s %s \\\n", @@ -1208,7 +1279,21 @@ endfor fputs (plot_stream, ";\n"); for i = 1:data_idx - if (is_image_data(i)) + if (have_3d_patch (i)) + ## Can't write 3d patch data as binary as can't plot more than + ## a single patch at a time and have to plot all patches together + ## so that the gnuplot depth ordering is done correctly + for j = 1 : 4 : columns(data{i}) + if (j != 1) + fputs (plot_stream, "\n\n"); + endif + fprintf (plot_stream, "%.15g %.15g %.15g %.15g\n", data{i}(:,j).'); + fprintf (plot_stream, "%.15g %.15g %.15g %.15g\n\n", data{i}(:,j+1).'); + fprintf (plot_stream, "%.15g %.15g %.15g %.15g\n", data{i}(:,j+2).'); + fprintf (plot_stream, "%.15g %.15g %.15g %.15g\n", data{i}(:,j+3).'); + endfor + fputs (plot_stream, "e\n"); + elseif (is_image_data(i)) fwrite (plot_stream, data{i}, "float32"); else __gnuplot_write_data__ (plot_stream, data{i}, nd, parametric(i),
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/plot/__interp_cube__.m Sat Apr 11 16:26:01 2009 +0200 @@ -0,0 +1,181 @@ +## Copyright (C) 2009 Martin Helm +## +## This program 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. +## +## This program 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 this program; if not, see http://www.gnu.org/licenses/gpl.html. +## +## Author: Martin Helm <martin@mhelm.de> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{vxyz}, @var{idx}, @var{frac}] =} __interp_cube__ (@var{x}, @var{y}, @var{z}, @var{val}, @var{v}) +## Undocumented internal function. +## @end deftypefn + +function [Vxyz, idx, frac] = __interp_cube__(x, y, z, val, v, req = "values" ) + if (ismatrix (x) && ndims (x) == 3 && ismatrix (y) && ndims (y) == 3 ... + && ismatrix (z) && ndims (z) == 3 && size_equal (x, y, z, val)) + x = squeeze (x(1,:,1))(:); + y = squeeze (y(:,1,1))(:); + z = squeeze (z(1,1,:))(:); + elseif (isvector (x) && isvector (y) && isvector (z) ) + x = x(:); + y = y(:); + z = z(:); + else + error("x, y, z have wrong dimensions"); + endif + if (size (val) != [length(x), length(y), length(z)]) + error ("val has wrong dimensions"); + endif + if (size (v, 2) != 3) + error ( "v has to be N*3 matrix"); + endif + if (!ischar (req)) + error ("Invalid request parameter use 'values', 'normals' or 'normals8'"); + endif + if (isempty (v)) + Vxyz = idx = frac = []; + return + endif + + switch req + case "values" + [Vxyz, idx, frac] = interp_cube_trilin (x, y, z, val, v); + case "normals" + [idx, frac] = cube_idx (x, y, z, v); + + dx = x(2:end) - x(1:end-1); + dy = y(2:end) - y(1:end-1); + dz = z(2:end) - z(1:end-1); + dx = 0.5 .* [dx;dx(end)](idx(:,2)); + dy = 0.5 .* [dy;dy(end)](idx(:,1)); + dz = 0.5 .* [dz;dz(end)](idx(:,3)); + + p000 = [v(:, 1) - dx, v(:, 2) - dy, v(:, 3) - dz]; + p100 = [v(:, 1) + dx, v(:, 2) - dy, v(:, 3) - dz]; + p010 = [v(:, 1) - dx, v(:, 2) + dy, v(:, 3) - dz]; + p001 = [v(:, 1) - dx, v(:, 2) - dy, v(:, 3) + dz]; + p011 = [v(:, 1) - dx, v(:, 2) + dy, v(:, 3) + dz]; + p101 = [v(:, 1) + dx, v(:, 2) - dy, v(:, 3) + dz]; + p110 = [v(:, 1) + dx, v(:, 2) + dy, v(:, 3) - dz]; + p111 = [v(:, 1) + dx, v(:, 2) + dy, v(:, 3) + dz]; + + v000 = interp_cube_trilin (x, y, z, val, p000); + v100 = interp_cube_trilin (x, y, z, val, p100); + v010 = interp_cube_trilin (x, y, z, val, p010); + v001 = interp_cube_trilin (x, y, z, val, p001); + v011 = interp_cube_trilin (x, y, z, val, p011); + v101 = interp_cube_trilin (x, y, z, val, p101); + v110 = interp_cube_trilin (x, y, z, val, p110); + v111 = interp_cube_trilin (x, y, z, val, p111); + + Dx = -v000 .+ v100 .- v010 .- v001 .- v011 .+ v101 .+ v110 .+ v111; + Dy = -v000 .- v100 .+ v010 .- v001 .+ v011 .- v101 .+ v110 .+ v111; + Dz = -v000 .- v100 .- v010 .+ v001 .+ v011 .+ v101 .- v110 .+ v111; + Vxyz = 0.5 .* [Dx./dx, Dy./dy, Dz./dz]; + case "normals8" + [idx, frac] = cube_idx (x, y, z, v); + + dx = x(2:end) - x(1:end-1); + dy = y(2:end) - y(1:end-1); + dz = z(2:end) - z(1:end-1); + dx = [dx;dx(end)](idx(:,2)); + dy = [dy;dy(end)](idx(:,1)); + dz = [dz;dz(end)](idx(:,3)); + [Dx, Dy, Dz, idx, frac] = interp_cube_trilin_grad (x, y, z, val, v); + Vxyz = [Dx./dx, Dy./dy, Dz./dz]; + otherwise + error ("Invalid request type '%s', use 'values', 'normals' or 'normals8'", req); + endswitch +endfunction + +function [Vxyz, idx, frac] = interp_cube_trilin(x, y, z, val, v) + [idx, frac] = cube_idx (x(:), y(:), z(:), v); + sval = size (val); + i000 = sub2ind (sval, idx(:, 1), idx(:, 2), idx(:, 3)); + i100 = sub2ind (sval, idx(:, 1)+1, idx(:, 2), idx(:, 3)); + i010 = sub2ind (sval, idx(:, 1), idx(:, 2)+1, idx(:, 3)); + i001 = sub2ind (sval, idx(:, 1), idx(:, 2), idx(:, 3)+1); + i101 = sub2ind (sval, idx(:, 1)+1, idx(:, 2), idx(:, 3)+1); + i011 = sub2ind (sval, idx(:, 1), idx(:, 2)+1, idx(:, 3)+1); + i110 = sub2ind (sval, idx(:, 1)+1, idx(:, 2)+1, idx(:, 3)); + i111 = sub2ind (sval, idx(:, 1)+1, idx(:, 2)+1, idx(:, 3)+1 ); + Bx = frac(:, 1); + By = frac(:, 2); + Bz = frac(:, 3); + Vxyz = ... + val( i000 ) .* (1 .- Bx) .* (1 .- By) .* (1 .- Bz) .+ ... + val( i100 ) .* Bx .* (1 .- By) .* (1 .- Bz) .+ ... + val( i010 ) .* (1 .- Bx) .* By .* (1 .- Bz) .+ ... + val( i001 ) .* (1 .- Bx) .* (1 .- By) .* Bz .+ ... + val( i011 ) .* (1 .- Bx) .* By .* Bz .+ ... + val( i101 ) .* Bx .* (1 .- By) .* Bz .+ ... + val( i110 ) .* Bx .* By .* (1 .- Bz) .+ ... + val( i111 ) .* Bx .* By .* Bz; +endfunction + +function [Dx, Dy, Dz, idx, frac] = interp_cube_trilin_grad(x, y, z, val, v) + [idx, frac] = cube_idx (x(:), y(:), z(:), v); + sval = size (val); + i000 = sub2ind (sval, idx(:, 1), idx(:, 2), idx(:, 3)); + i100 = sub2ind (sval, idx(:, 1)+1, idx(:, 2), idx(:, 3)); + i010 = sub2ind (sval, idx(:, 1), idx(:, 2)+1, idx(:, 3)); + i001 = sub2ind (sval, idx(:, 1), idx(:, 2), idx(:, 3)+1); + i101 = sub2ind (sval, idx(:, 1)+1, idx(:, 2), idx(:, 3)+1); + i011 = sub2ind (sval, idx(:, 1), idx(:, 2)+1, idx(:, 3)+1); + i110 = sub2ind (sval, idx(:, 1)+1, idx(:, 2)+1, idx(:, 3)); + i111 = sub2ind (sval, idx(:, 1)+1, idx(:, 2)+1, idx(:, 3)+1 ); + Bx = frac(:, 1); + By = frac(:, 2); + Bz = frac(:, 3); + Dx = ... + val( i000 ) .* -1 .* (1 .- By) .* (1 .- Bz) .+ ... + val( i100 ) .* (1 .- By) .* (1 .- Bz) .+ ... + val( i010 ) .* -1 .* By .* (1 .- Bz) .+ ... + val( i001 ) .* -1 .* (1 .- By) .* Bz .+ ... + val( i011 ) .* -1 .* By .* Bz .+ ... + val( i101 ) .* (1 .- By) .* Bz .+ ... + val( i110 ) .* By .* (1 .- Bz) .+ ... + val( i111 ) .* By .* Bz; + Dy = ... + val( i000 ) .* (1 .- Bx) .* -1 .* (1 .- Bz) .+ ... + val( i100 ) .* Bx .* -1 .* (1 .- Bz) .+ ... + val( i010 ) .* (1 .- Bx) .* (1 .- Bz) .+ ... + val( i001 ) .* (1 .- Bx) .* -1 .* Bz .+ ... + val( i011 ) .* (1 .- Bx) .* Bz .+ ... + val( i101 ) .* Bx .* -1 .* Bz .+ ... + val( i110 ) .* Bx .* (1 .- Bz) .+ ... + val( i111 ) .* Bx .* Bz; + Dz = ... + val( i000 ) .* (1 .- Bx) .* (1 .- By) .* -1 .+ ... + val( i100 ) .* Bx .* (1 .- By) .* -1 .+ ... + val( i010 ) .* (1 .- Bx) .* By .* -1 .+ ... + val( i001 ) .* (1 .- Bx) .* (1 .- By) .+ ... + val( i011 ) .* (1 .- Bx) .* By + ... + val( i101 ) .* Bx .* (1 .- By) .+ ... + val( i110 ) .* Bx .* By .* -1 .+ ... + val( i111 ) .* Bx .* By; +endfunction + +function [idx, frac] = cube_idx(x, y, z, v) + idx = zeros (size (v)); + frac = zeros (size (v)); + idx(:, 2) = lookup (x(2:end-1), v(:, 1)) + 1; + frac(:, 2) = (v(:, 1) - x(idx(:, 2)) )... + ./ (x(idx(:, 2)+1) - x(idx(:, 2))); + idx(:, 1) = lookup (y(2:end-1), v(:, 2)) + 1; + frac(:, 1) = (v(:, 2) - y(idx(:, 1))) ... + ./ (y(idx(:, 1)+1) - y(idx(:, 1))); + idx(:, 3) = lookup (z(2:end-1), v(:, 3)) + 1; + frac(:, 3) = (v(:, 3) - z(idx(:, 3))) ... + ./ (z(idx(:, 3)+1) - z(idx(:, 3))); +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/plot/__marching_cube__.m Sat Apr 11 16:26:01 2009 +0200 @@ -0,0 +1,505 @@ +## Copyright (C) 2009 Martin Helm +## +## This program 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. +## +## This program 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 this program; if not, see http://www.gnu.org/licenses/gpl.html. +## +## Author: Martin Helm <martin@mhelm.de> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{t}, @var{p}, @var{col}] =} __marching_cube__ (@var{x}, @var{y}, @var{z}, @var{c}, @var{iso}, @var{color}) +## Undocumented internal function. +## @end deftypefn + +## usage: [T, P] = marching_cube( XX, YY, ZZ, C, ISO) +## usage: [T, P, COL] = marching_cube( XX, YY, ZZ, C, ISO, COLOR) +## +## Calculates the triangulation T with points P for the isosurface +## with level ISO. XX, YY, ZZ are meshgrid like values for the +## cube and C holds the scalar values of the field, +## COLOR holds additinal scalar values for coloring the surface. +## The orientation of the triangles is choosen such that the +## normals point from the lower values to the higher values. +## +## edgeTable and triTable are taken from Paul Bourke +## (http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/) +## Based on tables by Cory Gene Bloyd +## +## Example: +## +## x = linspace(0, 2, 20); +## y = linspace(0, 2, 20); +## z = linspace(0, 2, 20); +## +## [ xx, yy, zz ] = meshgrid(x, y, z); +## +## c = (xx-.5).^2 + (yy-.5).^2 + (zz-.5).^2; +## [T, p] = marching_cube(xx, yy, zz, c, 0.5); +## trimesh(T, p(:, 1), p(:, 2), p(:, 3)); +## +## with jhandles you can also use the patch function to visualize +## +## clf +## pa = patch('Faces',T,'Vertices',p,'FaceVertexCData',p, ... +## 'FaceColor','interp', 'EdgeColor', 'none'); +## set(pa, "VertexNormals", -get(pa, "VertexNormals")) # revert normals +## view(-30, 30) +## set(pa, "FaceLighting", "gouraud") +## light( "Position", [1 1 5]) +## + +function [T, p, col] = __marching_cube__( xx, yy, zz, c, iso, colors) + + persistent edge_table=[]; + persistent tri_table=[]; + + calc_cols = false; + lindex = 4; + + if (isempty (tri_table) || isempty (edge_table)) + [edge_table, tri_table] = init_mc (); + endif + + if ((nargin != 5 && nargin != 6) || (nargout != 2 && nargout != 3)) + print_usage (); + endif + + if (!ismatrix (xx) || !ismatrix (yy) || !ismatrix (zz) || !ismatrix (c) || ... + ndims (xx) != 3 || ndims (yy) != 3 || ndims (zz) != 3 || ndims (c) != 3) + error ("xx, yy, zz, c have to be matrizes of dim 3"); + endif + + if (!size_equal (xx, yy, zz, c)) + error ("xx, yy, zz, c are not the same size"); + endif + + if (any (size (xx) < [2 2 2])) + error ("grid size has to be at least 2x2x2"); + endif + + if (!isscalar (iso)) + error ("iso needs to be scalar value"); + endif + + if (nargin == 6) + if ( !ismatrix (colors) || ndims (colors) != 3 || size (colors) != size (c) ) + error ( "color has to be matrix of dim 3 and of same size as c" ); + endif + calc_cols = true; + lindex = 5; + endif + + n = size (c) - 1; + + ## phase I: assign information to each voxel which edges are intersected by + ## the isosurface + cc = zeros (n(1), n(2), n(3), "uint16"); + cedge = zeros (size (cc), "uint16"); + + vertex_idx = {1:n(1), 1:n(2), 1:n(3); ... + 2:n(1)+1, 1:n(2), 1:n(3); ... + 2:n(1)+1, 2:n(2)+1, 1:n(3); ... + 1:n(1), 2:n(2)+1, 1:n(3); ... + 1:n(1), 1:n(2), 2:n(3)+1; ... + 2:n(1)+1, 1:n(2), 2:n(3)+1; ... + 2:n(1)+1, 2:n(2)+1, 2:n(3)+1; ... + 1:n(1), 2:n(2)+1, 2:n(3)+1 }; + + ## calculate which vertices have values higher than iso + for ii=1:8 + idx = c(vertex_idx{ii, :}) > iso; + cc(idx) = bitset (cc(idx), ii); + endfor + + cedge = edge_table(cc+1); # assign the info about intersected edges + id = find (cedge); # select only voxels which are intersected + if (isempty (id)) + T = p = col = []; + return + endif + + ## phase II: calculate the list of intersection points + xyz_off = [1, 1, 1; 2, 1, 1; 2, 2, 1; 1, 2, 1; 1, 1, 2; 2, 1, 2; 2, 2, 2; 1, 2, 2]; + edges = [1 2; 2 3; 3 4; 4 1; 5 6; 6 7; 7 8; 8 5; 1 5; 2 6; 3 7; 4 8]; + offset = sub2ind (size (c), xyz_off(:, 1), xyz_off(:, 2), xyz_off(:, 3)) -1; + pp = zeros (length (id), lindex, 12); + ccedge = [vec(cedge(id)), id]; + ix_offset=0; + for jj=1:12 + id__ = bitget (ccedge(:, 1), jj); + id_ = ccedge(id__, 2); + [ix iy iz] = ind2sub (size (cc), id_); + id_c = sub2ind (size (c), ix, iy, iz); + id1 = id_c + offset(edges(jj, 1)); + id2 = id_c + offset(edges(jj, 2)); + if (calc_cols) + pp(id__, 1:5, jj) = [vertex_interp(iso, xx(id1), yy(id1), zz(id1), ... + xx(id2), yy(id2), zz(id2), c(id1), c(id2), colors(id1), colors(id2)), ... + (1:size (id_, 1))' + ix_offset ]; + else + pp(id__, 1:4, jj) = [vertex_interp(iso, xx(id1), yy(id1), zz(id1), ... + xx(id2), yy(id2), zz(id2), c(id1), c(id2)), ... + (1:size (id_, 1))' + ix_offset ]; + endif + ix_offset += size (id_, 1); + endfor + + ## phase III: calculate the triangulation from the point list + T = []; + tri = tri_table(cc(id)+1, :); + for jj=1:3:15 + id_ = find (tri(:, jj)>0); + p = [id_, lindex*ones(size (id_, 1), 1),tri(id_, jj:jj+2)]; + if (!isempty (p)) + p1 = sub2ind (size (pp), p(:,1), p(:,2), p(:,3)); + p2 = sub2ind (size (pp), p(:,1), p(:,2), p(:,4)); + p3 = sub2ind (size (pp), p(:,1), p(:,2), p(:,5)); + T = [T; pp(p1), pp(p2), pp(p3)]; + endif + endfor + + p = []; + col = []; + for jj = 1:12 + idp = pp(:, lindex, jj) > 0; + if (any (idp)) + p(pp(idp, lindex, jj), 1:3) = pp(idp, 1:3, jj); + if (calc_cols) + col(pp(idp, lindex, jj),1) = pp(idp, 4, jj); + endif + endif + endfor +endfunction + +function p = vertex_interp(isolevel,p1x, p1y, p1z,... + p2x, p2y, p2z,valp1,valp2, col1, col2) + + if (nargin == 9) + p = zeros (length (p1x), 3); + elseif (nargin == 11) + p = zeros (length (p1x), 4); + else + error ("Wrong number of arguments"); + endif + mu = zeros (length (p1x), 1); + id = abs (valp1-valp2) < (10*eps) .* (abs (valp1) .+ abs (valp2)); + if (any (id)) + p(id, 1:3) = [ p1x(id), p1y(id), p1z(id) ]; + if (nargin == 11) + p(id, 4) = col1(id); + endif + endif + nid = !id; + if (any (nid)) + mu(nid) = (isolevel - valp1(nid)) ./ (valp2(nid) - valp1(nid)); + p(nid, 1:3) = [p1x(nid) + mu(nid) .* (p2x(nid) - p1x(nid)), ... + p1y(nid) + mu(nid) .* (p2y(nid) - p1y(nid)), ... + p1z(nid) + mu(nid) .* (p2z(nid) - p1z(nid))]; + if (nargin == 11) + p(nid, 4) = col1(nid) + mu(nid) .* (col2(nid) - col1(nid)); + endif + endif +endfunction + +function [edge_table, tri_table] = init_mc() + edge_table = [ + 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, ... + 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, ... + 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, ... + 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, ... + 0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c, ... + 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, ... + 0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac, ... + 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, ... + 0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c, ... + 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60, ... + 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc, ... + 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, ... + 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c, ... + 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, ... + 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc , ... + 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, ... + 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, ... + 0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0, ... + 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, ... + 0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650, ... + 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, ... + 0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, ... + 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, ... + 0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460, ... + 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, ... + 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0, ... + 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, ... + 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230, ... + 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, ... + 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190, ... + 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, ... + 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 ]; + + tri_table =[ + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1; + 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1; + 3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1; + 3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1; + 9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1; + 1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1; + 9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1; + 2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1; + 8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1; + 9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1; + 4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1; + 3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1; + 1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1; + 4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1; + 4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1; + 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1; + 1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1; + 5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1; + 2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1; + 9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1; + 0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1; + 2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1; + 10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1; + 4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1; + 5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1; + 5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1; + 9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1; + 0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1; + 1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1; + 10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1; + 8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1; + 2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1; + 7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1; + 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1; + 2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1; + 11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1; + 9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1; + 5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1; + 11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1; + 11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1; + 1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1; + 9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1; + 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1; + 2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1; + 0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1; + 5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1; + 6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1; + 0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1; + 3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1; + 6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1; + 5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1; + 1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1; + 10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1; + 6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1; + 1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1; + 8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1; + 7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1; + 3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1; + 5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1; + 0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1; + 9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1; + 8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1; + 5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1; + 0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1; + 6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1; + 10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1; + 10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1; + 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1; + 1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1; + 3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1; + 0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1; + 10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1; + 0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1; + 3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1; + 6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1; + 9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1; + 8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1; + 3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1; + 6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1; + 0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1; + 10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1; + 10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1; + 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1; + 2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1; + 7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1; + 7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1; + 2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1; + 1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1; + 11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1; + 8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1; + 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1; + 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1; + 10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1; + 2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1; + 6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1; + 7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1; + 2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1; + 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1; + 10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1; + 10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1; + 0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1; + 7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1; + 6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1; + 8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1; + 9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1; + 6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1; + 1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1; + 4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1; + 10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1; + 8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1; + 0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1; + 1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1; + 8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1; + 10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1; + 4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1; + 10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1; + 5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1; + 11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1; + 9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1; + 6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1; + 7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1; + 3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1; + 7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1; + 9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1; + 3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1; + 6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1; + 9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1; + 1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1; + 4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1; + 7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1; + 6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1; + 3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1; + 0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1; + 6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1; + 1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1; + 0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1; + 11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1; + 6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1; + 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1; + 9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1; + 1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1; + 1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1; + 10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1; + 0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1; + 5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1; + 10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1; + 11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1; + 0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1; + 9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1; + 7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1; + 2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1; + 8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1; + 9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1; + 9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1; + 1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1; + 9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1; + 9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1; + 5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1; + 0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1; + 10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1; + 2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1; + 0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1; + 0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1; + 9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1; + 5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1; + 3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1; + 5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1; + 8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1; + 0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1; + 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1; + 0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1; + 1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1; + 3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1; + 4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1; + 9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1; + 11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1; + 11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1; + 2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1; + 9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1; + 3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1; + 1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1; + 4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1; + 4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1; + 0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1; + 3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1; + 3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1; + 0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1; + 9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1; + 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1; + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 ] + 1; +endfunction \ No newline at end of file
--- a/scripts/plot/__patch__.m Sat Apr 11 16:05:16 2009 +0200 +++ b/scripts/plot/__patch__.m Sat Apr 11 16:26:01 2009 +0200 @@ -27,165 +27,280 @@ ## Author: Kai Habel -function [h, fail] = __patch__ (p, varargin) - - fail = false; +function [h, failed] = __patch__ (p, varargin) - if (nargin < 3) - fail = true; - return; - endif - - iarg = 1; - have_x = have_z = have_c = have_faces = false; - if (isnumeric (varargin{1})) - if (! isnumeric (varargin{2})) - fail = true; - return; - endif + failed = false; - x = varargin{1}; - y = varargin{2}; - have_x = true; - iarg += 2; - - if (nargin > 3 && ndims (varargin{3}) == 2 && ndims (x) == 2 - && size_equal(x, varargin{3}) && !ischar(varargin{3})) - z = varargin{3}; - have_z = true; - iarg++; - endif - elseif (ischar (varargin{1}) - && (strcmpi (varargin{1}, "faces") - || strcmpi (varargin{1}, "vertices"))) - if (! isnumeric (varargin{2})) - fail = true; - return; + if (isstruct (varargin{1})) + if (isfield (varargin{1}, "vertices") && isfield (varargin{1}, "faces")) + args{1} = "faces"; + args{2} = field(varargin{1}, "faces"); + args{3} = "vertices"; + args{4} = field(varargin{1}, "vertices"); + args{5} = "facevertexcdata"; + if (isfield (varargin{1}, "facevertexcdata")) + args{6} = field(varargin{1}, "facevertexcdata"); + else + args{6} = []; + endif + args = [args; varargin(2:end)]; + args = setdata (args); + else + failed = true; endif - - if (strcmpi (varargin{1}, "faces")) - faces = varargin{2}; - if (strcmpi (varargin{3}, "vertices")) - vert = varargin{4}; - have_faces = true; - endif - elseif (strcmpi (varargin{1}, "vertices")) - vert = varargin{2}; - if (strcmpi (varargin{3}, "faces")) - faces = varargin{4}; - have_faces = true; - endif - endif - if (!have_faces) - fail = true; - return; + elseif (isnumeric (varargin{1})) + if (nargin < 3 || ! isnumeric (varargin{2})) + failed = true; else - iarg += 4; - endif - endif + x = varargin{1}; + y = varargin{2}; + iarg = 3; + + if (nargin > 3 && ndims (varargin{3}) == 2 && ndims (x) == 2 + && size_equal(x, varargin{3}) && !ischar(varargin{3})) + z = varargin{3}; + iarg++; + else + z = []; + endif - if ((have_x || have_faces) && nargin > iarg) - if (isnumeric (varargin{iarg})) - c = varargin{iarg}; - have_c = true; - iarg++; + if (isvector (x)) + x = x(:); + y = y(:); + z = z(:); + endif + args{1} = "xdata"; + args{2} = x; + args{3} = "ydata"; + args{4} = y; + args{5} = "zdata"; + args{6} = z; + + if (isnumeric (varargin{iarg})) + c = varargin{iarg}; + iarg++; + + if (ndims (c) == 3 && size (c, 2) == 1) + c = permute (c, [1, 3, 2]); + endif - if (ndims (c) == 3 && size (c, 2) == 1) - c = permute (c, [1, 3, 2]); + if (isvector (c) && numel (c) == columns (x)) + if (isnan (c)) + args{7} = "facecolor"; + args{8} = [1, 1, 1]; + args{9} = "cdata"; + args{10} = c; + elseif (isnumeric (c)) + args{7} = "facecolor"; + args{8} = "flat"; + args{9} = "cdata"; + args{10} = c; + else + error ("patch: color value not valid"); + endif + elseif (size (c, ndims (c)) == 3) + args{7} = "facecolor"; + args{8} = "flat"; + args{9} = "cdata"; + args{10} = c; + else + ## Color Vectors + if (rows (c) != rows (x) || rows (c) != length (y)) + error ("patch: size of x, y, and c must be equal") + else + args{7} = "facecolor"; + args{8} = "interp"; + args{9} = "cdata"; + args{10} = []; + endif + endif + elseif (ischar (varargin{iarg}) && rem (nargin - iarg, 2) != 0) + ## Assume that any additional argument over an even number is + ## color string. + args{7} = "facecolor"; + args{8} = tolower (varargin{iarg}); + args{9} = "cdata"; + args{10} = []; + iarg++; + else + args{7} = "facecolor"; + args{8} = [0, 1, 0]; + args{9} = "cdata"; + args{10} = []; endif - elseif (ischar (varargin{iarg}) && rem (nargin - iarg, 2) != 0) - ## Assume that any additional argument over an even number is - ## color string. - c = tolower (varargin{iarg}); - have_c = true; - iarg++; + + args = [args, varargin(iarg:end)]; + args = setvertexdata (args); + endif + else + args = varargin; + if (any(cellfun (@(x) strcmpi(x,"faces") || strcmpi(x, "vertices"), args))) + args = setdata (args); + else + args = setvertexdata (args); endif endif - if (rem (nargin - iarg, 2) != 0) - fail = true; - return; + if (!failed) + h = __go_patch__ (p, args {:}); + + ## Setup listener functions + addlistener (h, "xdata", @update_data); + addlistener (h, "ydata", @update_data); + addlistener (h, "zdata", @update_data); + addlistener (h, "cdata", @update_data); + + addlistener (h, "faces", @update_fvc); + addlistener (h, "vertices", @update_fvc); + addlistener (h, "facevertexcdata", @update_fvc); + endif +endfunction + +function args = delfields(args, flds) + idx = cellfun (@(x) any (strcmpi (x, flds)), args); + idx = idx | [false, idx(1:end-1)]; + args (idx) = []; +endfunction + +function args = setdata (args) + args = delfields (args, {"xdata", "ydata", "zdata", "cdata"}); + nargs = length (args); + idx = find (cellfun (@(x) strcmpi (x, "faces"), args)) + 1; + if (idx > nargs) + faces = []; + else + faces = args {idx}; + endif + idx = find (cellfun (@(x) strcmpi (x, "vertices"), args)) + 1; + if (idx > nargs) + vert = []; + else + vert = args {idx}; + endif + idx = find (cellfun (@(x) strcmpi (x, "facecolor"), args)) + 1; + if (isempty(idx) || idx > nargs) + fc = "flat"; + else + fc = args {idx}; + endif + idx = find (cellfun (@(x) strcmpi (x, "facevertexcdata"), args)) + 1; + if (isempty(idx) || idx > nargs) + fvc = []; + else + fvc = args {idx}; endif - if (have_x) - if (isvector (x)) - x = x(:); - y = y(:); - if (have_z) - z = z(:); - endif - endif - [nr, nc] = size (x); - if (have_z) - vert = [x(:), y(:), z(:)]; - else - vert = [x(:), y(:)]; - endif - faces = reshape (1:numel(x), rows (x), columns (x)); - faces = faces'; - elseif (have_faces) - nr = size (faces, 2); - nc = size (faces, 1); - idx = faces .'; - t1 = isnan (idx); - if (any (t1(:))) - t2 = find (t1 != t1([2:end,end],:)); - idx (t1) = idx (t2 (cell2mat (cellfun (@(x) x(1)*ones(1,x(2)), + nr = size (faces, 2); + nc = size (faces, 1); + idx = faces .'; + t1 = isnan (idx); + if (any (t1(:))) + t2 = find (t1 != t1([2:end,end],:)); + idx (t1) = idx (t2 (cell2mat (cellfun (@(x) x(1)*ones(1,x(2)), mat2cell ([1 : nc; sum(t1)], 2, ones(1,nc)), - "UniformOutput", false)))); - endif - x = reshape (vert(:,1)(idx), size (idx)); - y = reshape (vert(:,2)(idx), size (idx)); - if (size(vert,2) > 2) - have_z = true; - z = reshape (vert(:,3)(idx), size (idx)); - endif + "UniformOutput", false)))); + endif + x = reshape (vert(:,1)(idx), size (idx)); + y = reshape (vert(:,2)(idx), size (idx)); + if (size(vert,2) > 2) + z = reshape (vert(:,3)(idx), size (idx)); else - error ("patch: not supported"); + z = []; endif - cargs = {}; - if (have_c) - if (ischar (c)) - cargs{1} = "facecolor"; - cargs{2} = c; - elseif (isvector (c) && numel (c) == nc) - if (isnan (c)) - cargs{1} = "facecolor"; - cargs{2} = [1, 1, 1]; - cargs{3} = "cdata"; - cargs{4} = c; - elseif (isnumeric (c)) - cargs{1} = "facecolor"; - cargs{2} = "flat"; - cargs{3} = "cdata"; - cargs{4} = c; + if (ischar (fc) && (strcmpi (fc, "flat") || strcmpi (fc, "interp"))) + if (size(fvc, 1) == nc || size (fvc, 1) == 1) + c = reshape (fvc, [1, size(fvc)]); + else + if (size(fvc, 2) == 3) + c = cat(3, reshape (fvc(idx, 1), size(idx)), + reshape (fvc(idx, 2), size(idx)), + reshape (fvc(idx, 3), size(idx))); else - error ("patch: color value not valid"); - endif - elseif (size (c, ndims (c)) == 3) - cargs{1} = "facecolor"; - cargs{2} = "flat"; - cargs{3} = "cdata"; - cargs{4} = c; - else - ## Color Vectors - if (rows (c2) != rows (x) || rows (c2) != length (y)) - error ("patch: size of x, y, and c must be equal") - else - cargs{1} = "facecolor"; - cargs{2} = "interp"; + c = reshape (fvc(idx), size(idx)); endif endif else - cargs{1} = "facecolor"; - cargs{2} = [0, 1, 0]; + c = []; + endif + args = {"xdata", x, "ydata", y, "zdata", z, "cdata", c, args{:}}; +endfunction + +function args = setvertexdata (args) + args = delfields (args, {"vertices", "faces", "facevertexcdata"}); + nargs = length (args); + idx = find (cellfun (@(x) strcmpi (x, "xdata"), args)) + 1; + if (idx > nargs) + x = []; + else + x = args {idx}; + endif + idx = find (cellfun (@(x) strcmpi (x, "ydata"), args)) + 1; + if (idx > nargs) + y = []; + else + y = args {idx}; + endif + idx = find (cellfun (@(x) strcmpi (x, "zdata"), args)) + 1; + if (isempty(idx) || idx > nargs) + z = []; + else + z = args {idx}; + endif + idx = find (cellfun (@(x) strcmpi (x, "facecolor"), args)) + 1; + if (isempty(idx) || idx > nargs) + fc = "flat"; + else + fc = args {idx}; + endif + idx = find (cellfun (@(x) strcmpi (x, "cdata"), args)) + 1; + if (isempty(idx) || idx > nargs) + c = []; + else + c = args {idx}; endif - h = __go_patch__ (p, "xdata", x, "ydata", y, "faces", faces, - "vertices", vert, cargs{:}, varargin{iarg:end}); - if (have_z) - set (h, "zdata", z); + [nr, nc] = size (x); + if (!isempty (z)) + vert = [x(:), y(:), z(:)]; + else + vert = [x(:), y(:)]; + endif + faces = reshape (1:numel(x), rows (x), columns (x)); + faces = faces'; + + if (ischar (fc) && (strcmpi (fc, "flat") || strcmpi (fc, "interp"))) + if (ndims (c) == 3) + fvc = reshape (c, size (c, 1) * size (c, 2), size(c, 3)); + else + fvc = c(:); + endif + else + fvc = []; endif - + + args = {"faces", faces, "vertices", vert, "facevertexcdata", fvc, args{:}}; +endfunction + +function update_data (h, d) + update_handle (h, false); +endfunction + +function update_fvc (h, d) + update_handle (h, true); endfunction + +function update_handle (h, isfv) + persistent recursive = false; + + if (! recursive) + recursive = true; + f = get (h); + if (isfvc) + set (h, setvertexdata ([fieldnames(f), struct2cell(f)].'(:)){:}); + else + set (h, setdata ([fieldnames(f), struct2cell(f)].'(:)){:}); + endif + recursive = false; + endif +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/plot/isocolors.m Sat Apr 11 16:26:01 2009 +0200 @@ -0,0 +1,168 @@ +## Copyright (C) 2009 Martin Helm +## +## This program 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. +## +## This program 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 this program; if not, see http://www.gnu.org/licenses/gpl.html. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{cd}] =} isocolors (@var{c}, @var{v}) +## @deftypefnx {Function File} {[@var{cd}] =} isocolors (@var{x}, @var{y}, @var{z}, @var{c}, @var{v}) +## @deftypefnx {Function File} {[@var{cd}] =} isocolors (@var{x}, @var{y}, @var{z}, @var{r}, @var{g}, @var{b}, @var{v}) +## @deftypefnx {Function File} {[@var{cd}] =} isocolors (@var{r}, @var{g}, @var{b}, @var{v}) +## @deftypefnx {Function File} {[@var{cd}] =} isocolors (@dots{}, @var{p}) +## @deftypefnx {Function File} isocolors (@dots{}) +## +## If called with one output argument and the first input argument +## @var{c} is a three--dimensional array that contains color values and +## the second input argument @var{v} keeps the vertices of a geometry +## then return a matrix @var{cd} with color data information for the +## geometry at computed points +## @command{[x, y, z] = meshgrid (1:l, 1:m, 1:n)}. The output argument +## @var{cd} can be taken to manually set FaceVertexCData of a patch. +## +## If called with further input arguments @var{x}, @var{y} and @var{z} +## which are three--dimensional arrays of the same size than @var{c} +## then the color data is taken at those given points. Instead of the +## color data @var{c} this function can also be called with RGB values +## @var{r}, @var{g}, @var{b}. If input argumnets @var{x}, @var{y}, +## @var{z} are not given then again @command{meshgrid} computed values +## are taken. +## +## Optionally, the patch handle @var{p} can be given as the last input +## argument to all variations of function calls instead of the vertices +## data @var{v}. Finally, if no output argument is given then directly +## change the colors of a patch that is given by the patch handle +## @var{p}. +## +## For example, +## @example +## function [] = isofinish (p) +## set (gca, "DataAspectRatioMode", "manual", \ +## "DataAspectRatio", [1 1 1]); +## set (p, "FaceColor", "interp"); +## ## set (p, "FaceLighting", "flat"); +## ## light ("Position", [1 1 5]); ## Available with JHandles +## endfunction +## +## N = 15; ## Increase number of vertices in each direction +## iso = .4; ## Change isovalue to .1 to display a sphere +## lin = linspace (0, 2, N); +## [x, y, z] = meshgrid (lin, lin, lin); +## c = abs ((x-.5).^2 + (y-.5).^2 + (z-.5).^2); +## figure (); ## Open another figure window +## +## subplot (2, 2, 1); view (-38, 20); +## [f, v] = isosurface (x, y, z, c, iso); +## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "none"); +## cdat = rand (size (c)); ## Compute random patch color data +## isocolors (x, y, z, cdat, p); ## Directly set colors of patch +## isofinish (p); ## Call user function isofinish +## +## subplot (2, 2, 2); view (-38, 20); +## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "none"); +## [r, g, b] = meshgrid (lin, 2-lin, 2-lin); +## cdat = isocolors (x, y, z, c, v); ## Compute color data vertices +## set (p, "FaceVertexCData", cdat); ## Set color data manually +## isofinish (p); +## +## subplot (2, 2, 3); view (-38, 20); +## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "none"); +## cdat = isocolors (r, g, b, c, p); ## Compute color data patch +## set (p, "FaceVertexCData", cdat); ## Set color data manually +## isofinish (p); +## +## subplot (2, 2, 4); view (-38, 20); +## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "none"); +## r = g = b = repmat ([1:N] / N, [N, 1, N]); ## Black to white +## cdat = isocolors (x, y, z, r, g, b, v); +## set (p, "FaceVertexCData", cdat); +## isofinish (p); +## @end example +## +## @seealso{isosurface, isonormals, isocaps} +## +## @end deftypefn + +## Author: Martin Helm <martin@mhelm.de> + +function varargout = isocolors(varargin) + calc_rgb = false; + switch nargin + case 2 + c = varargin{1}; + vp = varargin{2}; + x = 1:size (c, 2); + y = 1:size (c, 1); + z = 1:size (c, 3); + case 4 + calc_rgb = true; + R = varargin{1}; + G = varargin{2}; + B = varargin{3}; + vp = varargin{4}; + x = 1:size (R, 1); + y = 1:size (R, 2); + z = 1:size (R, 3); + case 5 + x = varargin{1}; + y = varargin{2}; + z = varargin{3}; + c = varargin{4}; + vp = varargin{5}; + case 7 + calc_rgb = true; + x = varargin{1}; + y = varargin{2}; + z = varargin{3}; + R = varargin{4}; + G = varargin{5}; + B = varargin{6}; + vp = varargin{7}; + otherwise + print_usage (); + endswitch + if (ismatrix (vp) && size (vp,2) == 3) + pa = []; + v = vp; + elseif ( ishandle (vp) ) + pa = vp; + v = get (pa, "Vertices"); + else + error("Last argument is no vertex list and no patch handle"); + endif + if ( calc_rgb ) + new_col = zeros (size (v, 1), 3); + new_col(:, 1) = __interp_cube__ (x, y, z, R, v, "values" ); + new_col(:, 2) = __interp_cube__ (x, y, z, G, v, "values" ); + new_col(:, 3) = __interp_cube__ (x, y, z, B, v, "values" ); + else + new_col = __interp_cube__ (x, y, z, c, v, "values" ); + endif + switch nargout + case 0 + if (!isempty (pa)) + set (pa, "FaceVertexCData", new_col); + endif + case 1 + varargout = {new_col}; + otherwise + print_usage (); + endswitch +endfunction + +%!test +%! [x, y, z] = meshgrid (0:.5:2, 0:.5:2, 0:.5:2); +%! c = (x-.5).^2 + (y-.5).^2 + (z-.5).^2; +%! [f, v] = isosurface (x, y, z, c, .4); +%! cdat = isocolors (x, y, z, c, v); +%! assert (size (cdat, 1) == size (v, 1)); +## Can't create a patch handle for tests without a figure
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/plot/isonormals.m Sat Apr 11 16:26:01 2009 +0200 @@ -0,0 +1,78 @@ +## Copyright (C) 2009 Martin Helm +## +## This program 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. +## +## This program 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 this program; if not, see http://www.gnu.org/licenses/gpl.html. +## +## Author: Martin Helm <martin@mhelm.de> + +## usage: NORMALS = isonormals(X,Y,Z,V,VERT) +## usage: NORMALS = isonormals(V,VERT) +## usage: NORMALS = isonormals(V,P) +## usage: NORMALS = isonormals(X,Y,Z,V,P) +## usage: NORMALS = isonormals(...,'negate') +## usage: isonormals(V,P) +## usage: isonormals(X,Y,Z,V,P) +## + +function varargout = isonormals(varargin) + na = nargin; + negate = false; + if (ischar (varargin{nargin})) + na = nargin-1; + if (strcmp (lower (varargin{nargin}), "negate")) + negate = true; + else + error ("Unknown option '%s'", varargin{nargin}); + endif + endif + switch na + case 2 + c = varargin{1}; + vp = varargin{2}; + x = 1:size (c, 2); + y = 1:size (c, 1); + z = 1:size (c, 3); + case 5 + x = varargin{1}; + y = varargin{2}; + z = varargin{3}; + c = varargin{4}; + vp = varargin{5}; + otherwise + print_usage (); + endswitch + if (ismatrix (vp) && size (vp,2) == 3) + pa = []; + v = vp; + elseif (ishandle (vp)) + pa = vp; + v = get (pa, "Vertices"); + else + error ("Last argument is no vertex list and no patch handle"); + endif + if (negate) + normals = -__interp_cube__ (x, y, z, c, v, "normals"); + else + normals = __interp_cube__ (x, y, z, c, v, "normals"); + endif + switch nargout + case 0 + if (!isempty (pa)) + set (pa, "VertexNormals", normals); + endif + case 1 + varargout = {normals}; + otherwise + print_usage (); + endswitch +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/plot/isosurface.m Sat Apr 11 16:26:01 2009 +0200 @@ -0,0 +1,215 @@ +## Copyright (C) 2009 Martin Helm +## +## This program 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. +## +## This program 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 this program; if not, see http://www.gnu.org/licenses/gpl.html. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{fv}] =} isosurface (@var{val}, @var{iso}) +## @deftypefnx {Function File} {[@var{fv}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}) +## @deftypefnx {Function File} {[@var{fv}] =} isosurface (@dots{}, "noshare", "verbose") +## @deftypefnx {Function File} {[@var{fvc}] =} isosurface (@dots{}, @var{col}) +## @deftypefnx {Function File} {[@var{f}, @var{v}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}) +## @deftypefnx {Function File} {[@var{f}, @var{v}, @var{c}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}, @var{col}) +## @deftypefnx {Function File} {} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}, @var{col}, @var{opt}) +## +## If called with one output argument and the first input argument +## @var{val} is a three--dimensional array that contains the data of an +## isosurface geometry and the second input argument @var{iso} keeps the +## isovalue as a scalar value then return a structure array @var{fv} +## that contains the fields @var{Faces} and @var{Vertices} at computed +## points @command{[x, y, z] = meshgrid (1:l, 1:m, 1:n)}. The output +## argument @var{fv} can directly be taken as an input argument for the +## @command{patch} function. +## +## If called with further input arguments @var{x}, @var{y} and @var{z} +## which are three--dimensional arrays with the same size than @var{val} +## then the volume data is taken at those given points. +## +## The string input argument "noshare" is only for compatibility and +## has no effect. If given the string input argument +## "verbose" then print messages to the command line interface about the +## current progress. +## +## If called with the input argument @var{col} which is a +## three-dimensional array of the same size than @var{val} then take +## those values for the interpolation of coloring the isosurface +## geometry. Add the field @var{FaceVertexCData} to the structure +## array @var{fv}. +## +## If called with two or three output arguments then return the +## information about the faces @var{f}, vertices @var{v} and color data +## @var{c} as seperate arrays instead of a single structure array. +## +## If called with no output argument then directly process the +## isosurface geometry with the @command{patch} command. +## +## For example +## +## @example +## [x, y, z] = meshgrid (1:5, 1:5, 1:5); +## val = rand (5, 5, 5); +## isosurface (x, y, z, val, .5); +## @end example +## +## will directly draw a random isosurface geometry in a graphics window. +## Another example for an isosurface geometry with different additional +## coloring +## +## @example +## N = 15; ## Increase number of vertices in each direction +## iso = .4; ## Change isovalue to .1 to display a sphere +## lin = linspace (0, 2, N); +## [x, y, z] = meshgrid (lin, lin, lin); +## c = abs ((x-.5).^2 + (y-.5).^2 + (z-.5).^2); +## figure (); ## Open another figure window +## +## subplot (2, 2, 1); view (-38, 20); +## [f, v] = isosurface (x, y, z, c, iso); +## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "none"); +## set (gca, "DataAspectRatioMode","manual", "DataAspectRatio", [1 1 1]); +## # set (p, "FaceColor", "green", "FaceLighting", "phong"); +## # light ("Position", [1 1 5]); ## Available with the JHandles package +## +## subplot (2, 2, 2); view (-38, 20); +## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "blue"); +## set (gca, "DataAspectRatioMode","manual", "DataAspectRatio", [1 1 1]); +## # set (p, "FaceColor", "none", "FaceLighting", "phong"); +## # light ("Position", [1 1 5]); +## +## subplot (2, 2, 3); view (-38, 20); +## [f, v, c] = isosurface (x, y, z, c, iso, y); +## p = patch ("Faces", f, "Vertices", v, "FaceVertexCData", c, \ +## "FaceColor", "interp", "EdgeColor", "none"); +## set (gca, "DataAspectRatioMode","manual", "DataAspectRatio", [1 1 1]); +## # set (p, "FaceLighting", "phong"); +## # light ("Position", [1 1 5]); +## +## subplot (2, 2, 4); view (-38, 20); +## p = patch ("Faces", f, "Vertices", v, "FaceVertexCData", c, \ +## "FaceColor", "interp", "EdgeColor", "blue"); +## set (gca, "DataAspectRatioMode","manual", "DataAspectRatio", [1 1 1]); +## # set (p, "FaceLighting", "phong"); +## # light ("Position", [1 1 5]); +## @end example +## +## @seealso{isocolors, isonormals, isocaps} +## +## @end deftypefn + +## Author: Martin Helm <martin@mhelm.de> + +function varargout = isosurface(varargin) + + if (nargin < 2 || nargin > 8 || nargout > 3) + print_usage (); + endif + + calc_colors = false; + f = v = c = []; + verbose = false; + noshare = false; + if (nargin >= 5) + x = varargin{1}; + y = varargin{2}; + z = varargin{3}; + val = varargin{4}; + iso = varargin{5}; + if (nargin >= 6 && ismatrix (varargin{6})) + colors = varargin{6}; + calc_colors = true; + endif + else + val = varargin{1}; + [n1, n2, n3] = size (val); + [x, y, z] = meshgrid (1:n1, 1:n2, 1:n3); + iso = varargin{2}; + if (nargin >= 3 && ismatrix (varargin{3})) + colors = varargin{3}; + calc_colors = true; + endif + endif + if (calc_colors) + if (nargout == 2) + warning ( "Colors will be calculated, but you did not specify an output argument for it!" ); + endif + [fvc.faces, fvc.vertices, fvc.facevertexcdata] = __marching_cube__ (x, y, z, val, iso, colors); + else + [fvc.faces, fvc.vertices] = __marching_cube__ (x, y, z, val, iso); + endif + + if (isempty (fvc.vertices) || isempty (fvc.faces)) + warning ( "The resulting triangulation is empty" ); + endif + + switch (nargout) + case 0 + ## plot the calculated surface + newplot (); + if (calc_colors) + pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices, + "FaceVertexCData", fvc.facevertexcdata, + "FaceColor", "flat", "EdgeColor", "none"); + else + pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices, + "FaceColor", "g", "EdgeColor", "k"); + endif + if (! ishold ()) + set (gca(), "view", [-37.5, 30], + "xgrid", "on", "ygrid", "on", "zgrid", "on"); + endif + case 1 + varargout = {fvc}; + case 2 + varargout = {fvc.faces, fvc.vertices}; + case 3 + varargout = {fvc.faces, fvc.vertices, fvc.facevertexcdata}; + otherwise + print_usage (); + endswitch + +endfunction + + +%!shared x, y, z, val +%! [x, y, z] = meshgrid (0:1, 0:1, 0:1); ## Points for single +%! val = [0, 0; 0, 0]; ## cube and a 3--dim +%! val(:,:,2) = [0, 0; 1, 0]; ## array of values +%!test +%! fv = isosurface (x, y, z, val, 0.3); +%! assert (isfield (fv, "vertices"), true); +%! assert (isfield (fv, "faces"), true); +%! assert (size (fv.vertices), [3 3]); +%! assert (size (fv.faces), [1 3]); +%!test +%! fvc = isosurface (x, y, z, val, .3, y); +%! assert (isfield (fvc, "vertices"), true); +%! assert (isfield (fvc, "faces"), true); +%! assert (isfield (fvc, "facevertexcdata"), true); +%! assert (size (fvc.vertices), [3 3]); +%! assert (size (fvc.faces), [1 3]); +%! assert (size (fvc.facevertexcdata), [3 1]); +%!test +%! [f, v] = isosurface (x, y, z, val, .3); +%! assert (size (f), [1 3]); +%! assert (size (v), [3 3]); +%!test +%! [f, v, c] = isosurface (x, y, z, val, .3, y); +%! assert (size (f), [1 3]); +%! assert (size (v), [3 3]); +%! assert (size (c), [3 1]); + +%!demo +%! clf +%! [x,y,z] = meshgrid(-2:0.5:2, -2:0.5:2, -2:0.5:2); +%! v = x.^2 + y.^2 + z.^2; +%! isosurface (x, y, z, v, 1) \ No newline at end of file
--- a/scripts/plot/patch.m Sat Apr 11 16:05:16 2009 +0200 +++ b/scripts/plot/patch.m Sat Apr 11 16:26:01 2009 +0200 @@ -19,7 +19,8 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} patch () ## @deftypefnx {Function File} {} patch (@var{x}, @var{y}, @var{c}) -## @deftypefnx {Function File} {} patch (@var{x}, @var{y}, @var{c}, @var{opts}) +## @deftypefnx {Function File} {} patch (@var{x}, @var{y}, @var{z}, @var{c}) +## @deftypefnx {Function File} {} patch (@var{fv}) ## @deftypefnx {Function File} {} patch ('Faces', @var{f}, 'Vertices', @var{v}, @dots{}) ## @deftypefnx {Function File} {} patch (@dots{}, @var{prop}, @var{val}) ## @deftypefnx {Function File} {} patch (@var{h}, @dots{}) @@ -29,7 +30,11 @@ ## ## For a uniform colored patch, @var{c} can be given as an RGB vector, ## scalar value referring to the current colormap, or string value (for -## example, "r" or "red"). +## example, "r" or "red"). +## +## If passed a structure @var{fv} contain the fields "vertices", "faces" +## and optionally "facevertexcdata", create the patch based on these +## properties. ## @end deftypefn ## Author: jwe @@ -42,15 +47,14 @@ unwind_protect axes (h); - [tmp, fail] = __patch__ (h, varargin{:}); + [tmp, failed] = __patch__ (h, varargin{:}); + if (failed) + print_usage (); + endif unwind_protect_cleanup axes (oldh); end_unwind_protect - if (fail) - print_usage (); - endif - if (nargout > 0) retval = tmp; endif @@ -93,6 +97,19 @@ %! patch('Faces',fac,'Vertices',vert,'FaceColor','r'); %!demo +%! ## Specify vertices and faces separately +%! clf +%! t1 = (1/16:1/8:1)'*2*pi; +%! t2 = ((1/16:1/16:1)' + 1/32)*2*pi; +%! x1 = sin(t1) - 0.8; +%! y1 = cos(t1); +%! x2 = sin(t2) + 0.8; +%! y2 = cos(t2); +%! vert = [x1, y1; x2, y2]; +%! fac = [1:8,NaN(1,8);9:24]; +%! patch('Faces',fac,'Vertices',vert,'FaceVertexCData', [0, 1, 0; 0, 0, 1]); + +%!demo %! ## Property change on multiple patches %! clf %! t1 = (1/16:1/8:1)'*2*pi; @@ -104,3 +121,18 @@ %! h = patch([x1,x2],[y1,y2],cat (3,[0,0],[1,0],[0,1])); %! pause (1); %! set (h, 'FaceColor', 'r'); + +%!demo +%! clf +%! vertices = [0, 0, 0; +%! 1, 0, 0; +%! 1, 1, 0; +%! 0, 1, 0; +%! 0.5, 0.5, 1]; +%! faces = [1, 2, 5; +%! 2, 3, 5; +%! 3, 4, 5; +%! 4, 1, 5]; +%! patch('Vertices', vertices, 'Faces', faces, ... +%! 'FaceVertexCData', jet(4), 'FaceColor', 'flat') +%! view (-37.5, 30)