HiveBrain v1.2.0
Get Started
← Back to all entries
patternMinor

Interpolation on marching cubes algorithm

Submitted by: @import:stackexchange-cs··
0
Viewed 0 times
algorithmcubesinterpolationmarching

Problem

I'm currently learning how isosurfaces are extrated from volumetric data. I already understood the original marching cubes algorithm which is based on 3D-voxel data which stores only values of either 1 or 0. For learning purposes I implemented the easier version of it for 2D-volumetric data called marching squares algorithm. Results are shown below:

After all I read that the marching cubes algorithm can benefit of float values in the voxels instead of the binary values in terms of interpolating the vertices. However I do not understand how this should work simply because I'm missing the rules of which values should only be assumed for the voxels and also how the interpolation based on them would be done.

Could someone give me a fundamental view of how interpolation with the marching cubes algorithm is asumed to work? An example with 2D would be fine but showing how it works in 3D would be even nicer.

PS:
For references, I read such stuff here.

Solution

I worked it out on my own. I'm trying to write it down correctly so it's readable for newcomers.

Introduction

First of all, interpolation itself requires that there are float values and not binary values. This gives the ability to create different edge-crossing-vertices. How it's going I'm explaining now.

The interpolation which is needed here is linear interpolation.
For understanding purposes I'll explain it here too.

General Interpolation

After all, what is interpolation? Well good question.
Interpolation is described as name for algorithms which smooth a given input of data e.g. points (with n-dimensions). Simpler there are points searched between the points for a given interpolation value. First it's required to give the points a base so that's clear where points are searched for in between. After that a specific algorithm for the interpolation has to be chosen. There are different kinds of interpolation algorithm out there each for a specific application.
Best known is linear interpolation which results for a value one point on a line which crosses two input points.

The formula for linear interpolation (2D here) is as following (simply maths):
$$
y = y_{0} + (value - x_{0} ) * \frac{y_{1} - y_{0} }{x_{1} - x_{0}}
$$
with $$ value \in \mathbb{R} \wedge 0

  • The value of the voxels must be only between zero and one



  • There must be a boundary value be declared which describes whether a voxel is supposed to be in surface.



Looking back to the blog the boundary value is called the iso surface value.
This value describes as explained whether a voxel is being asumed to be in the surface so for acquiring the edge crossing. This is important to know because without it you wouldn't find out where even a edge-crossing computation is needed.

The interpolation is done on each coordinate of the two voxels position where and crossing-edge vertex is.
The implementation is quite easy and looks like this (extracted from here):

/*
   Linearly interpolate the position where an isosurface cuts
   an edge between two vertices, each with their own scalar value
*/
XYZ VertexInterp(isolevel,p1,p2,valp1,valp2)
double isolevel;
XYZ p1,p2;
double valp1,valp2;
{
   double mu;
   XYZ p;

   if (ABS(isolevel-valp1) < 0.00001)
      return(p1);
   if (ABS(isolevel-valp2) < 0.00001)
      return(p2);
   if (ABS(valp1-valp2) < 0.00001)
      return(p1);
   mu = (isolevel - valp1) / (valp2 - valp1);
   p.x = p1.x + mu * (p2.x - p1.x);
   p.y = p1.y + mu * (p2.y - p1.y);
   p.z = p1.z + mu * (p2.z - p1.z);

   return(p);
}


Conclusion

The interpolation of the edge-crossing vertices is a nice feature for smoothing the presentation of voxels.However the interpolation relies on that there are already floating values which means that it cannot smooth your existing binary data. For such cases different algorithms and approaches are needed.

Code Snippets

/*
   Linearly interpolate the position where an isosurface cuts
   an edge between two vertices, each with their own scalar value
*/
XYZ VertexInterp(isolevel,p1,p2,valp1,valp2)
double isolevel;
XYZ p1,p2;
double valp1,valp2;
{
   double mu;
   XYZ p;

   if (ABS(isolevel-valp1) < 0.00001)
      return(p1);
   if (ABS(isolevel-valp2) < 0.00001)
      return(p2);
   if (ABS(valp1-valp2) < 0.00001)
      return(p1);
   mu = (isolevel - valp1) / (valp2 - valp1);
   p.x = p1.x + mu * (p2.x - p1.x);
   p.y = p1.y + mu * (p2.y - p1.y);
   p.z = p1.z + mu * (p2.z - p1.z);

   return(p);
}

Context

StackExchange Computer Science Q#71102, answer score: 3

Revisions (0)

No revisions yet.