Some GNU Octave scripts

  1. Isosurface
  2. Marching Cube Algorithm
  3. A simple way to create patches with lighting
  4. Marching Tetrahedron Algorithm (in progress)
  5. 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.

Isosurface function

top
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).

Example 1: isosurface without additional coloring

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
isosurface

Example 2: isosurface with additional coloring

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
isosurface

Example 3: a realistic dataset

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])
aneurism

A simple way to create patches with lighting

top

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

calc_shading.m

The demo calc_shading will result in the following images when using backend("fltk"):

isosurface

isosurface

Marching Cube Algorithm

top

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)

isosurface

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).

isosurface

Marching Tetrahedron Algorithm

top
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