Dear Octave forge developers,
I am relatively new to octave (had a first look at it last year) and noticed
that in the field of 3d graphical data processing there are some missing
functions like isosurface and other related functions.
I made preliminary substitutes for some of them for my own use which I would
like to share if there is interest.
The isosurface function is based on a vectorized form of the marching cube
algorithm. The original marching cube algorithm was patented in the past, but
the patent expired a few years ago, so I can see no legal problem using it in
the open source community.
If someone sees a problem problem with that, I implement at the moment the
marching tetrahedron algorithm which can then be used instead.
For the next time I plan to implement isocaps too.
Code documentation and error handling needs to be improved, no plan when I
will find time, but I will do it.
Since it seems to be that jhandles is the only graphical backend which has all
necessary features for visualization of the corresponding 3d patches including
facelighting and gouraud shading, the first tests I made are exclusively done
with it (jhandles 0.3.4 and octave 3.0.3).
I also made some preliminary tests with fltk and octave 3.1.53+ built from the
mercurial repository that shows a noticable performance enhacement for the
algorithm (the subscripting with the development sources are much faster) and
the basic visualization features are working well without lighting (the light
function is missing).
Attached you can find the implementations to have a look at it. You can tell
me if they are appropriate or not. There is also a simple octave package which
is available on my private homepage. (I am not familiar with the packaging
system but it worked for me).
http://www.mhelm.de/octave/pkg/visualize3d-0.1.2.tar.gz
With some simple examples
http://www.mhelm.de/octave/index.html
Kind regards
Martin Helm
## 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 <[email protected]>
## usage: [Vxyz, idx, frac] = interp_cube(x, y, z, val, v)
##
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(x) == size(y) ...
&& size(y) == size(z) )
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' or 'normals'");
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);
Bx = By = Bz = 0.5;
Dx = ...
v000 .* -1 .* (1 .- By) .* (1 .- Bz) .+ ...
v100 .* (1 .- By) .* (1 .- Bz) .+ ...
v010 .* -1 .* By .* (1 .- Bz) .+ ...
v001 .* -1 .* (1 .- By) .* Bz .+ ...
v011 .* -1 .* By .* Bz .+ ...
v101 .* (1 .- By) .* Bz .+ ...
v110 .* By .* (1 .- Bz) .+ ...
v111 .* By .* Bz;
Dy = ...
v000 .* (1 .- Bx) .* -1 .* (1 .- Bz) .+ ...
v100 .* Bx .* -1 .* (1 .- Bz) .+ ...
v010 .* (1 .- Bx) .* (1 .- Bz) .+ ...
v001 .* (1 .- Bx) .* -1 .* Bz .+ ...
v011 .* (1 .- Bx) .* Bz .+ ...
v101 .* Bx .* -1 .* Bz .+ ...
v110 .* Bx .* (1 .- Bz) .+ ...
v111 .* Bx .* Bz;
Dz = ...
v000 .* (1 .- Bx) .* (1 .- By) .* -1 .+ ...
v100 .* Bx .* (1 .- By) .* -1 .+ ...
v010 .* (1 .- Bx) .* By .* -1 .+ ...
v001 .* (1 .- Bx) .* (1 .- By) .+ ...
v011 .* (1 .- Bx) .* By + ...
v101 .* Bx .* (1 .- By) .+ ...
v110 .* Bx .* By .* -1 .+ ...
v111 .* Bx .* By;
Vxyz = 2 .* [ 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' or 'normals'", 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## 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 <[email protected]>
## usage: NC = isocolors(X,Y,Z,C,VERT)
## usage: NC = isocolors(X,Y,Z,R,G,B,VERT)
## usage: NC = isocolors(C,VERT)
## usage: NC = isocolors(R,G,B,VERT)
## usage: NC = isocolors(...,PatchHandle)
## usage: isocolors(...,PatchHandle)
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## 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 <[email protected]>
## 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## 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 <[email protected]>
## usage: FV = isosurface(X,Y,Z,V,ISO)
## FV = isosurface(V,ISO)
## FVC = isosurface(...,C)
## FV = isosurface(...,'noshare')
## FV = isosurface(...,'verbose')
## [F,V] = isosurface(...)
## [F,V,C] = isosurface(...)
## isosurface(...)
##
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
[f, v, c] = marching_cube(x, y, z, val, iso, colors);
else
[f, v] = marching_cube(x, y, z, val, iso);
endif
if ( isempty(f) || isempty(v) )
warning( "The resulting triangulation is empty" );
endif
switch (nargout)
case 0
# plot the calculated surface
if ( calc_colors )
pa = patch('Faces', f,'Vertices',v ,'FaceVertexCData', c);
else
pa = patch('Faces', f,'Vertices', v);
endif
case 1
## this is a little bit special, we can return a struct but patch cannot
## use it, with this workaround a syntax like
## patch(isosurface(...){:}) is possible
if ( calc_colors )
varargout = {{"Faces", f, "Vertices", v, "FaceVertexCData", c }};
else
varargout = {{"Faces", f, "Vertices", v}};
endif
case 2
varargout = {f, v};
case 3
varargout = {f, v, c};
otherwise
print_usage();
endswitch
endfunction
## 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 <[email protected]>
## 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 edgeTable=[];
persistent triTable=[];
calc_cols = false;
lindex = 4;
if ( isempty(triTable) || isempty(edgeTable) )
[edgeTable, triTable] = 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(xx) != size(yy) || size(yy) != size(zz)|| size(zz) != size(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 less than iso
for ii=1:8
idx = c(vertex_idx{ii, :}) < iso;
cc(idx) = bitset(cc(idx), ii);
endfor
cedge = edgeTable(cc+1); # assign the info about intersected edges
id = find(cedge); # select only voxels which are intersected
## 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 = triTable(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 ( ! isempty(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));
p(id, 1:3) = [ p1x(id), p1y(id), p1z(id) ];
mu(!id) = (isolevel - valp1(!id)) ./ (valp2(!id) - valp1(!id));
p(!id, 1:3) = [p1x(!id) + mu(!id) .* (p2x(!id) - p1x(!id)), ...
p1y(!id) + mu(!id) .* (p2y(!id) - p1y(!id)), ...
p1z(!id) + mu(!id) .* (p2z(!id) - p1z(!id))];
if ( nargin == 11 )
p(!id, 4) = col1(!id) + mu(!id) .* (col2(!id) - col1(!id));
endif
endfunction
function [edgeTable, triTable] = init_mc()
edgeTable = [
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 ];
triTable =[
-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------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
