The functions are available as an octave package visualize3d-0.1.5.tar.gz. Includes experimental isocolors and isonormals. The use of this package is deprecated because the iso* functions (newer versions of them) are now included in the octave core since version 3.2.
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(...)
isosurface.m (isosurface.html)
This is an very early draft of an isosurface function for Octave (based on the marching_cube function), which will be as Matlab compatible as possible. To show some possibilities some examples of its usage are shown below. In order to be able to use the graphical possibilities, you need a graphics backend which is able to work with filled 3d patches (JHandles or fltk will do it).
N=50;
iso=.5;
x = linspace(0, 2, N);
y = linspace(0, 2, N);
z = linspace(0, 2, N);
[ xx, yy, zz ] = meshgrid(x, y, z);
c = abs((xx-.5).^2 + (yy-.5).^2 + (zz-.5).^2);
[T, p] = isosurface(xx, yy, zz, c, iso);
figure();
pa = patch('Faces',T,'Vertices',p); # draws the surface with standard settings
view(-30, 30)
set (pa, "EdgeColor", "none")
set(gca, "DataAspectRatioMode", "manual")
set(gca, "DataAspectRatio", [1 1 1])
set(pa, "FaceLighting", "phong")
light( "Position", [1 1 5]) # available with jhandles
N=50;
iso=.5;
x = linspace(0, 2, N);
y = linspace(0, 2, N);
z = linspace(0, 2, N);
[ xx, yy, zz ] = meshgrid(x, y, z);
c = abs((xx-.5).^2 + (yy-.5).^2 + (zz-.5).^2);
[T, p, col] = isosurface(xx, yy, zz, c, iso, yy);
figure();
pa = patch('Faces',T,'Vertices',p,'FaceVertexCData',col, ...
'FaceColor','interp', 'EdgeColor', 'none');
view(-30, 30)
set (pa, "EdgeColor", "none")
set(gca, "DataAspectRatioMode", "manual")
set(gca, "DataAspectRatio", [1 1 1])
set(pa, "FaceLighting", "phong")
light( "Position", [1 1 5]) # available with jhandles
Uses aneurism dataset from http://www.cb.uu.se/~tc18/code_data_set/3D_greyscale/aneurism.raw.gz
Results in 175060 triangles with iso value 100 (with jhandles you need to enhance the default java heap space, I use -Xmx512M)
fid = fopen("aneurism.raw", "r");
[c, cnt] = fread(fid, Inf, "uchar");
fclose(fid);
c = reshape(c, 256, 256, 256);
[xx, yy, zz] = meshgrid(1:256);
tic;
[T, p, col] = isosurface(xx, yy, zz, c, 100, zz);
toc
clf
pa = patch('Faces',T,'Vertices',p,'FaceVertexCData',col, 'FaceColor','interp', 'EdgeColor', 'none');
view(-30, 30)
colormap(jet(256))
set(gca, "DataAspectRatioMode", "manual")
set(gca, "DataAspectRatio", [1 1 1])
set(pa, "FaceLighting", "phong")
light( "Position", [0 0 500])
Note: Gnuplot and the fltk backend in octave have no standard way to define lights at the moment. The following script is a small workaround to create some graphics with an improved 3d feeling.
Function File: [CDATA] = calc_shading (AMB, DIF, SPEC, SHINE,
COLORS, NORMALS, LVEC, VVEC)
If called with one output argument returns calculated color data
for use with the patch function.
AMB, DIF, SPEC, SHINE are the ambient, diffuse, specular and
specular exponent properties for the lighting.
COLORS can be a one or three column matrix with the same length as
the NORMALS matrix of normal vectors.
A scalar or a RGB row vector can also be used for COLORS. The
object has thean a single color.
LVEC: The light direction
VVEC: The view direction
Type
backend ("fltk")
demo calc_shading
to see how it works.
Without 'backend ("fltk")' it will also work but is very slow.
See also: isocolors, isonormals
The demo calc_shading will result in the following images when using backend("fltk"):
Note: For a long time it was not possible to use the marching cube algorithm in open source software because it is patented. A few years ago the patent expired, so there is now no longer a restriction to use it.
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');
view(-30, 30)
set(pa, "FaceLighting", "gouraud")
light( "Position", [1 1 5])
marching_cubes.m (marching_cubes.html)
If you have octaviz installed you can also use vtk_trisurf(T, p(:, 1), p(:, 2), p(:, 3)) to show the
triangulation (much faster and better results).
usage: [T, P] = polygonize_tetra( TH, X, Y, Z, C, ISO) Calculates the triangulation T with points P for the isosurface with level ISO. TH, X, Y, Z are vectors like for the delaunay3 function and C holds the scalar values of the field. based on a c program from Paul Bourke (http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/) Example: x = rand(1000, 1); y = rand(1000, 1); z = rand(1000, 1); TH = delaunay3(x, y, z); c = (xx-.5).^2 + (yy-.5).^2 + (zz-.5).^2; [T, p] = polygonize_tetra(TH, x, y, z, c, 0.5); trimesh(T, p(:, 1), p(:, 2), p(:, 3));top