- Isosurface
- Marching Cube Algorithm
- A simple way to create patches with lighting
- Marching Tetrahedron Algorithm (in progress)
- Finite Elements

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