HUMIM Tracker: /home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/core/containers/mesh/common/util/mesh_smooth

Go to the documentation of this file.
00001 #ifndef OPENTISSUE_CORE_CONTAINERS_MESH_COMMON_UTIL_MESH_SMOOTH_ISOSURFACE_H
00002 #define OPENTISSUE_CORE_CONTAINERS_MESH_COMMON_UTIL_MESH_SMOOTH_ISOSURFACE_H
00003 //
00004 // OpenTissue Template Library
00005 // - A generic toolbox for physics-based modeling and simulation.
00006 // Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
00007 //
00008 // OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
00009 //
00010 #include <OpenTissue/configuration.h>
00011 
00012 #include <OpenTissue/core/containers/mesh/common/util/mesh_compute_mesh_center.h>
00013 #include <OpenTissue/core/containers/mesh/common/util/mesh_compute_mesh_minimum_coord.h>
00014 #include <OpenTissue/core/containers/mesh/common/util/mesh_compute_mesh_maximum_coord.h>
00015 #include <OpenTissue/core/containers/mesh/common/util/mesh_deformation_modifiers.h>
00016 
00017 #include <cassert>
00018 #include <vector>
00019 
00020 namespace OpenTissue
00021 {
00022   namespace mesh
00023   {
00024 
00025     /*
00026     First, I will assume you are more or less familiar with the Marching Cubes algorithm.
00027 
00028     The cube description I use looks like this:
00029 
00030     Y
00031     |
00032     3----2----2
00033     /|        /|
00034     11 |      10 |
00035     /  3      /  1
00036     7----6----6   |
00037     |   |     |   |
00038     |   0----0|---1---X
00039     7  /      5  /
00040     | 8       | 9
00041     |/        |/
00042     4----4----5
00043     /
00044     Z
00045 
00046     Given that, here's the basic procedure I have developed:
00047 
00048     You know that Marching Cubes treats a volume as a regular 3D grid of "cubes", with
00049     8 voxel values of the volume at the corners. For each of the 12 edges in any given cube,
00050     the isovalue can fall above, below, or between the values at the endpoints. This means
00051     that any edge can have at most one interpolated vertex "on" it.
00052 
00053     The main data structure for my extraction algorithm is a "Cell" illustrated by the
00054     following diagram:
00055 
00056     Y
00057     |
00058     |
00059     Vy
00060     |
00061     |
00062     Vs---Vx----X
00063     /
00064     Vz
00065     /
00066     Z
00067 
00068     The "CellInfo" struct below represents this structure. If you imagine putting lots
00069     of these cells next to each other in space, you can quickly see they form a grid similar
00070     to a 3D grid of cubes, but with no duplicated edges:
00071 
00072     |   |   |
00073     +-- +-- +--
00074     /   /   /
00075     |   |   |
00076     +-- +-- +--
00077     /   /   /
00078 
00079     Two 2D arrays of cells, each the same size as a "slice" in the 3D volume, will between them contain
00080     all necessary information for any vertices and triangles generated by the first set of cubes
00081     in the volume. Each 2D array of cells I'll from hereon refer to as a "grid".
00082 
00083     Looking at the Cell diagram above and CellInfo struct below:
00084     Vx is the integer index into a list of vertices for the vertex interpolated along that edge - the x-axis.
00085     That is, if the isovalue falls between the values of the endpoints of that edge, then a vertex is
00086     interpolated and stored in a list, and its index in that list is stored in Vx.
00087     It is -1 initially, or if no vertex has yet been interpolated. Likewise for the other 2 axes.
00088 
00089     During extraction, each "cube" of 8 values are visited, and used to construct the index into
00090     'triangleEdgeTable', which has the list of triplets of edges forming at most 5 triangles for
00091     the cube. As each of the 3 edges forming the triangle are visited, the "EdgeCellInfo" struct
00092     is used to see which cell that edge maps to in one of the 2 grids.
00093 
00094     For example, edge 9 in one cell maps to Vz in the cell next to it, in the same grid; i.e. to
00095     the right of it in the cube diagram above. Edge 5 in one cell maps to Vy in the cell to the
00096     right of it and in the next grid.
00097 
00098     Now that we know which cell the edge maps to, we can see if the vertex on that edge has been
00099     interpolated by checking to see if the corresponding 'vertIndex' field is -1. If it is, we
00100     interpolate the vertex, add it to the list, and store its index. If not, then the vertex has
00101     already been interpolated and we just use it.
00102 
00103     We now have the 3 indexes of the vertices that form a triangle, and can add that triangle to the
00104     list of triangles. The normal of the triangle is computed and stored as well as accumulated
00105     within each vertex. Later, the per-vertex normals are averaged and normalized.
00106 
00107     After processing 2 slices of the volume, the "top" grid (i.e. the 2nd grid) can be reused as the
00108     "bottom" grid for the next set of slices. It is swapped with the 1st grid, and the new 2nd grid is
00109     reset and reused.
00110 
00111     I have also implemented a vertex "smoothing" operation, which basically collects vertices that
00112     are within 1/2 unit of a "corner" of the cube and averages them, forming as few as half as many
00113     vertices and triangles and preserving the overall shape and morphology of the surface.
00114 
00115     The 'smoothedIndex' field holds a pointer to the smoothed vertex in the list, and the 'vertIndex'
00116     holds pointers to which smoothed vertex is closest to the corresponding edge. Everything else
00117     works the same, except degenerate triangles (triangles where two or more of the vertices are
00118     the same) are not added to the triangle list.
00119 
00120     At the end, smoothed vertices are averaged, and all vertices are pre-scaled according to a caller-specified
00121     scaling factor (to account for isotropism) and translated so as to be in a default [-0.5,0.5] space
00122     centered at the origin.
00123 
00124     */
00125 
00126     namespace detail
00127     {
00128 
00129       /* Marching Cubes Isosurface extraction routines.
00130       The Marching Cubes algorithm is patented and was originally developed
00131       by William Lorenson and Harvey Cline:
00132 
00133       "Marching Cubes: A High Resolution 3D Surface Construction Algorithm",
00134       William E. Lorensen and Harvey E. Cline,
00135       Computer Graphics (Proceedings of SIGGRAPH '87), Vol. 21, No. 4, pp. 163-169.
00136       */
00137 
00138       /* The triangleEdgeTable was obtained from open-source code written by Cory Bloyd
00139       which was based on a marching cubes description by Paul Bourke:
00140       http://www.swin.edu.au/astronomy/pbourke/modelling/polygonise/
00141 
00142       Everything else in here is mine! :-)
00143       */
00144 
00145 
00146 
00151       template<typename mesh_type>
00152       class CellInfo
00153       {
00154       public:
00155 
00156         typedef typename mesh_type::vertex_handle vertex_handle;
00157 
00158       public:
00159 
00160         vertex_handle  m_vertex_handle[3];  
00161         vertex_handle  m_smoothed_handle;   
00162 
00163       public:
00164 
00165         void clear()
00166         {
00167           vertex_handle null_handle;
00168           m_vertex_handle[0]  = null_handle;
00169           m_vertex_handle[1]  = null_handle;
00170           m_vertex_handle[2]  = null_handle;
00171           m_smoothed_handle = null_handle;
00172         }
00173       };
00174 
00179       class EdgeCellInfo
00180       {
00181       public:
00182         int m_cell;
00183         int m_grid;
00184         int m_axis;
00185         int m_offset[3];
00186       };
00187 
00188       template<typename cell_info_iterator>
00189       void clear_cell_grid(cell_info_iterator begin,cell_info_iterator end)
00190       {
00191         for(cell_info_iterator c=begin;c!=end;++c)
00192           c->clear();
00193       }
00194 
00195     }// namespace detail
00196 
00206     template<typename grid_type,typename mesh_type>
00207     void smooth_isosurface(
00208       grid_type const & phi
00209       , typename grid_type::value_type const & isolevel
00210       , mesh_type & mesh
00211       , size_t sub_sample  = 1u
00212       , bool smoothing   = true
00213       )
00214     {
00215       typedef detail::CellInfo<mesh_type>           cell_info_type;
00216       typedef std::vector<cell_info_type>           cell_grid_type;
00217       typedef typename grid_type::value_type         value_type;
00218 
00219       typedef typename mesh_type::math_types        math_types;
00220       typedef typename math_types::value_traits     value_traits;
00221       typedef typename math_types::vector3_type     vector3_type;
00222       typedef typename math_types::real_type        real_type;
00223 
00224       typedef typename mesh_type::vertex_handle     vertex_handle;
00225       typedef typename mesh_type::vertex_iterator   vertex_iterator;
00226 
00227       static int triangle_edge_table[256][16] =
00228       {
00229         {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00230         {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00231         {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00232         {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00233         {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00234         {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00235         {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00236         {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
00237         {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00238         {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00239         {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00240         {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
00241         {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00242         {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
00243         {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
00244         {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00245         {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00246         {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00247         {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00248         {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
00249         {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00250         {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
00251         {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
00252         {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
00253         {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00254         {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
00255         {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
00256         {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
00257         {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
00258         {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
00259         {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
00260         {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
00261         {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00262         {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00263         {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00264         {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
00265         {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00266         {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
00267         {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
00268         {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
00269         {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00270         {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
00271         {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
00272         {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
00273         {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
00274         {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
00275         {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
00276         {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
00277         {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00278         {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
00279         {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
00280         {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00281         {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
00282         {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
00283         {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
00284         {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
00285         {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
00286         {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
00287         {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
00288         {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
00289         {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
00290         {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
00291         {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
00292         {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00293         {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00294         {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00295         {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00296         {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
00297         {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00298         {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
00299         {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
00300         {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
00301         {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00302         {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
00303         {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
00304         {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
00305         {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
00306         {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
00307         {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
00308         {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
00309         {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00310         {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
00311         {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
00312         {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
00313         {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
00314         {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
00315         {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
00316         {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
00317         {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
00318         {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
00319         {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
00320         {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
00321         {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
00322         {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
00323         {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
00324         {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
00325         {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00326         {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
00327         {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
00328         {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
00329         {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
00330         {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
00331         {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00332         {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
00333         {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
00334         {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
00335         {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
00336         {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
00337         {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
00338         {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
00339         {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
00340         {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00341         {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
00342         {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
00343         {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
00344         {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
00345         {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
00346         {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
00347         {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
00348         {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00349         {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
00350         {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
00351         {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
00352         {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
00353         {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
00354         {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00355         {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
00356         {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00357         {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00358         {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00359         {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00360         {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
00361         {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00362         {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
00363         {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
00364         {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
00365         {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00366         {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
00367         {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
00368         {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
00369         {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
00370         {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
00371         {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
00372         {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
00373         {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00374         {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
00375         {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
00376         {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
00377         {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
00378         {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
00379         {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
00380         {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
00381         {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
00382         {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00383         {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
00384         {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
00385         {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
00386         {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
00387         {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
00388         {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00389         {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00390         {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
00391         {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
00392         {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
00393         {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
00394         {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
00395         {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
00396         {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
00397         {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
00398         {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
00399         {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
00400         {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
00401         {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
00402         {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
00403         {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
00404         {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
00405         {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
00406         {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
00407         {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
00408         {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
00409         {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
00410         {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
00411         {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
00412         {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
00413         {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
00414         {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
00415         {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
00416         {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00417         {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
00418         {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
00419         {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00420         {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00421         {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00422         {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
00423         {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
00424         {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
00425         {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
00426         {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
00427         {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
00428         {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
00429         {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
00430         {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
00431         {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
00432         {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
00433         {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00434         {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
00435         {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
00436         {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00437         {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
00438         {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
00439         {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
00440         {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
00441         {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
00442         {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
00443         {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
00444         {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00445         {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
00446         {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
00447         {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
00448         {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
00449         {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
00450         {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00451         {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
00452         {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00453         {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
00454         {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
00455         {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
00456         {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
00457         {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
00458         {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
00459         {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
00460         {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
00461         {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
00462         {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
00463         {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
00464         {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00465         {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
00466         {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
00467         {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00468         {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00469         {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00470         {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
00471         {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
00472         {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00473         {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
00474         {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
00475         {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00476         {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00477         {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
00478         {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00479         {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
00480         {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00481         {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00482         {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00483         {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
00484         {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
00485       };
00486 
00487 
00491       static int edge_vertex_table[12][2] =
00492       {
00493         {0,1}, {1,2}, {3,2}, {0,3},
00494         {4,5}, {5,6}, {7,6}, {4,7},
00495         {0,4}, {1,5}, {2,6}, {3,7}
00496       };
00497 
00498       //--- Clear any data stored in mesh.
00499       mesh.clear();
00500 
00501       //--- Determine constants used by the algorithm
00502       if( sub_sample <= 0 )  sub_sample = 1;
00503       if( sub_sample > 16 )  sub_sample = 16;
00504       size_t  sampled_I   = phi.I()/sub_sample;
00505       size_t  sampled_J   = phi.J()/sub_sample;
00506       //size_t  sampled_K   = phi.K()/sub_sample;
00507       size_t  grid_width  = sampled_I-1;
00508       size_t  grid_height = sampled_J-1;
00509       // NOTE: What about some sort of warning here for the user? An assertion perhaps.
00510       if( grid_width == 0 || grid_height == 0 )
00511         return;
00512 
00513       size_t cnt_cells = sampled_I * sampled_J;
00514 
00515       //--- Determine offsets needed by the algorithm...
00516       detail::EdgeCellInfo     edge_cell_lut[12];
00517       {
00518         //--- compute the edge-cell lookup: given any edge in the cube (there are 12),
00519         //--- this table indicates which cell the interpolated vertex for that edge is
00520         //--- (as indexed from the current cell), which grid that cell is in (either the
00521         //--- first or the second), and which axis the edge maps to (either x, y, or z).
00522         //--- grid is which grid of CellInfos to use: the first or the second
00523         //--  cell is how to index from the current cell to the proper cell for that edge
00524         //--- axis is which axis the needed vertex is on: 0=>x, 1=>y, 2=>z
00525         edge_cell_lut[ 0].m_cell = 0;          edge_cell_lut[ 0].m_grid = 0;    edge_cell_lut[ 0].m_axis = 0;
00526         edge_cell_lut[ 1].m_cell = 1;          edge_cell_lut[ 1].m_grid = 0;    edge_cell_lut[ 1].m_axis = 1;
00527         edge_cell_lut[ 2].m_cell = sampled_I;   edge_cell_lut[ 2].m_grid = 0;    edge_cell_lut[ 2].m_axis = 0;
00528         edge_cell_lut[ 3].m_cell = 0;          edge_cell_lut[ 3].m_grid = 0;    edge_cell_lut[ 3].m_axis = 1;
00529         edge_cell_lut[ 4].m_cell = 0;          edge_cell_lut[ 4].m_grid = 1;    edge_cell_lut[ 4].m_axis = 0;
00530         edge_cell_lut[ 5].m_cell = 1;          edge_cell_lut[ 5].m_grid = 1;    edge_cell_lut[ 5].m_axis = 1;
00531         edge_cell_lut[ 6].m_cell = sampled_I;   edge_cell_lut[ 6].m_grid = 1;    edge_cell_lut[ 6].m_axis = 0;
00532         edge_cell_lut[ 7].m_cell = 0;          edge_cell_lut[ 7].m_grid = 1;    edge_cell_lut[ 7].m_axis = 1;
00533         edge_cell_lut[ 8].m_cell = 0;          edge_cell_lut[ 8].m_grid = 0;    edge_cell_lut[ 8].m_axis = 2;
00534         edge_cell_lut[ 9].m_cell = 1;          edge_cell_lut[ 9].m_grid = 0;    edge_cell_lut[ 9].m_axis = 2;
00535         edge_cell_lut[10].m_cell = 1+sampled_I; edge_cell_lut[10].m_grid = 0;    edge_cell_lut[10].m_axis = 2;
00536         edge_cell_lut[11].m_cell = sampled_I;   edge_cell_lut[11].m_grid = 0;    edge_cell_lut[11].m_axis = 2;
00537 
00538         //--- The offset is how much to add to the coordinates of the current
00539         //--- cell to get the right interpolated  vertex coordinate
00540         edge_cell_lut[ 0].m_offset[0] = 0; edge_cell_lut[ 0].m_offset[1] = 0; edge_cell_lut[ 0].m_offset[2] = 0;
00541         edge_cell_lut[ 1].m_offset[0] = 1; edge_cell_lut[ 1].m_offset[1] = 0; edge_cell_lut[ 1].m_offset[2] = 0;
00542         edge_cell_lut[ 2].m_offset[0] = 0; edge_cell_lut[ 2].m_offset[1] = 1; edge_cell_lut[ 2].m_offset[2] = 0;
00543         edge_cell_lut[ 3].m_offset[0] = 0; edge_cell_lut[ 3].m_offset[1] = 0; edge_cell_lut[ 3].m_offset[2] = 0;
00544         edge_cell_lut[ 4].m_offset[0] = 0; edge_cell_lut[ 4].m_offset[1] = 0; edge_cell_lut[ 4].m_offset[2] = 1;
00545         edge_cell_lut[ 5].m_offset[0] = 1; edge_cell_lut[ 5].m_offset[1] = 0; edge_cell_lut[ 5].m_offset[2] = 1;
00546         edge_cell_lut[ 6].m_offset[0] = 0; edge_cell_lut[ 6].m_offset[1] = 1; edge_cell_lut[ 6].m_offset[2] = 1;
00547         edge_cell_lut[ 7].m_offset[0] = 0; edge_cell_lut[ 7].m_offset[1] = 0; edge_cell_lut[ 7].m_offset[2] = 1;
00548         edge_cell_lut[ 8].m_offset[0] = 0; edge_cell_lut[ 8].m_offset[1] = 0; edge_cell_lut[ 8].m_offset[2] = 0;
00549         edge_cell_lut[ 9].m_offset[0] = 1; edge_cell_lut[ 9].m_offset[1] = 0; edge_cell_lut[ 9].m_offset[2] = 0;
00550         edge_cell_lut[10].m_offset[0] = 1; edge_cell_lut[10].m_offset[1] = 1; edge_cell_lut[10].m_offset[2] = 0;
00551         edge_cell_lut[11].m_offset[0] = 0; edge_cell_lut[11].m_offset[1] = 1; edge_cell_lut[11].m_offset[2] = 0;
00552       }
00553       //--- Notice the convention that lowercase i,j,k are local subsampled
00554       //--- indices, uppercasre I,J,K are the non-sampled indices
00555       //--- (ie. the ones you should use to access the value phi(i,j,k))
00556       //--- The conversion formulas are:
00557       //---
00558       //---  I = i*sub_sample
00559       //---  J = j*sub_sample
00560       //---  K = k*sub_sample
00561       //---
00562       //--- Compute offsets to access values of nodes in phi. That is:
00563       //---
00564       //---  phi(           I,           J,K)
00565       //---  phi(I+sub_sample,           J,K)
00566       //---  phi(I+sub_sample,J+sub_sample,K)
00567       //---  phi(           I,J+sub_sample,K)
00568       //---
00569       size_t offset[4];
00570       {
00571         offset[0] = 0;                       //--- offset to subsampled location (  i,j,k)
00572         offset[1] = sub_sample;              //--- offset to subsampled location (i+1,j,k)
00573         offset[2] = (1+phi.I())*sub_sample;  //--- offset to subsampled location (i+1,j+1,k)
00574         offset[3] = phi.I()*sub_sample;      //--- offset to subsampled location (  i,j+1,k)
00575       }
00576 
00577       //--- Allocate memory for internal data structures and initialize any internal data
00578       static cell_grid_type   grid_storage[2];     //--- The actual memory used to store the grid cells!!!
00579       cell_info_type * cell_grid[2];        //--- Used to swap cell grids when moving on to next z-plane.
00580       cell_info_type * grid[2];             //--- Used to traverse cells in the grids!!!
00581       value_type  const * data = phi.data();   //--- pointer to location phi(0,0,0)
00582       {
00583         grid_storage[0].resize( cnt_cells );
00584         grid_storage[1].resize( cnt_cells );
00585         cell_grid[0] = &(grid_storage[0][0]);
00586         cell_grid[1] = &(grid_storage[1][0]);
00587         clear_cell_grid(cell_grid[0],cell_grid[0]+cnt_cells);
00588       }
00589 
00590       size_t k,K;
00591       for( k = 0, K = 0; K < (phi.K()-sub_sample) ; ++k, K+=sub_sample ) //--- run over sub-sampled z-planes
00592       {
00593         grid[0] = cell_grid[0];
00594         grid[1] = cell_grid[1];
00595         clear_cell_grid(cell_grid[1],cell_grid[1]+cnt_cells);
00596         for(size_t j = 0; j < grid_height; ++j )
00597         {
00598           value_type const * vox0 = data + j*sub_sample*phi.I() + K*phi.I()*phi.J();                //--- pointer to sub sampled location phi(i,j,k),    with i=0
00599           value_type const * vox1 = data + j*sub_sample*phi.I() + (K+sub_sample)*phi.I()*phi.J();   //--- pointer to sub sampled location phi(i,j,k+1),  with i=0
00600 
00601           for(size_t i = 0; i < grid_width; ++i )
00602           {
00603             //--- Compute the index for which set of triangles to use
00604             int index = 0;
00605             if( (float)(vox0[offset[0]]) <= isolevel )      index |= 1;    //--- phi(  i,  j,k  )
00606             if( (float)(vox0[offset[1]]) <= isolevel )      index |= 2;    //--- phi(i+1,  j,k  )
00607             if( (float)(vox0[offset[2]]) <= isolevel )      index |= 4;    //--- phi(i+1,j+1,k  )
00608             if( (float)(vox0[offset[3]]) <= isolevel )      index |= 8;    //--- phi(  i,j+1,k  )
00609             if( (float)(vox1[offset[0]]) <= isolevel )      index |= 16;   //--- phi(  i,  j,k+1)
00610             if( (float)(vox1[offset[1]]) <= isolevel )      index |= 32;   //--- phi(i+1,  j,k+1)
00611             if( (float)(vox1[offset[2]]) <= isolevel )      index |= 64;   //--- phi(i+1,j+1,k+1)
00612             if( (float)(vox1[offset[3]]) <= isolevel )      index |= 128;  //--- phi(  i,j+1,k+1)
00613             //--- Quick check for the most common indexes of all (i.e. no triangles formed)
00614             if( index==0 || index==255 )
00615             {
00616               //--- Advance to next pixel/cell
00617               vox0 += sub_sample;
00618               vox1 += sub_sample;
00619               grid[0]++;
00620               grid[1]++;
00621               continue;
00622             }
00623             int * triangle_edge = triangle_edge_table[index];
00624             for(int t = 0; ; t += 3 )
00625             {
00626               if( triangle_edge[t] == - 1 )
00627                 break;
00628               std::vector<vertex_handle> vertex_handle(3);
00629               //--- Examine each of the three edges of the cell
00630               for( int e = 0; e < 3; ++e )
00631               {
00632                 int E = triangle_edge[t+e];                                   //--- E???
00633                 detail::EdgeCellInfo * edge = edge_cell_lut + E;              //--- this is just a faster version of: edge = &(edge_cell_lut[E]);
00634                 int axis = edge->m_axis;                                      //--- the axis on where the intersection lies
00635                 cell_info_type * cell = &(grid[edge->m_grid][edge->m_cell]);  //--- get cell corresponding to the intersected edge
00636                 if( cell->m_vertex_handle[axis].is_null() )                   //--- see if we already got a vertex on this edge
00637                 {
00638                   //--- Get the 2 vertices that make up the edge.
00639                   //--- These values are on [0,7], where values [0,3] refer
00640                   //--- to grid 0, and [4,7] refer to grid 1
00641                   int v0 = edge_vertex_table[E][0];
00642                   int v1 = edge_vertex_table[E][1];
00643                   //--- get the corresponding phi values of the 2 vertices
00644                   value_type val0 = ( v0 < 4 ) ? vox0[offset[v0]] : vox1[offset[v0-4]];
00645                   value_type val1 = ( v1 < 4 ) ? vox0[offset[v1]] : vox1[offset[v1-4]];
00646                   //--- Compute the interpolant (i.e. the fraction the new vertex is between the 2 others)
00647                   value_type frac = (isolevel-val0)/(val1-val0);
00648                   //--- Add a new vertex to the list of vertices
00649 
00650                   vector3_type coord(0,0,0);
00651                   cell_info_type * new_cell = 0;
00652 
00653                   switch( axis )
00654                   {
00655                   case 0:
00656                     coord(0) = static_cast<real_type>((i + edge->m_offset[0]) + frac);
00657                     coord(1) = static_cast<real_type> (j + edge->m_offset[1]);
00658                     coord(2) = static_cast<real_type> (k + edge->m_offset[2]);
00659                     new_cell = &(grid[edge->m_grid][edge->m_cell+1]);
00660                     break;
00661                   case 1:
00662                     coord(0) = static_cast<real_type> (i + edge->m_offset[0]);
00663                     coord(1) = static_cast<real_type>((j + edge->m_offset[1]) + frac);
00664                     coord(2) = static_cast<real_type> (k + edge->m_offset[2]);
00665                     new_cell = &(grid[edge->m_grid][edge->m_cell+sampled_I]);
00666                     break;
00667                   default:
00668                     coord(0) = static_cast<real_type> (i + edge->m_offset[0]);
00669                     coord(1) = static_cast<real_type> (j + edge->m_offset[1]);
00670                     coord(2) = static_cast<real_type>((k + edge->m_offset[2]) + frac);
00671                     new_cell = &(grid[edge->m_grid+1][edge->m_cell]);
00672                     break;
00673                   }
00674 
00675                   //--- if we're doing vertex smoothing, we need to accumulate
00676                   //--- all potential verts into a single "smoothed" vert
00677                   if( smoothing )
00678                   {
00679                     if( frac <= static_cast<value_type>(0.5) )//--- hmmm, why this?
00680                     {
00681                       new_cell = cell;
00682                     }
00683                     if( new_cell->m_smoothed_handle.is_null() )
00684                     {
00685                       //--- the smoothed vertex has not been started yet. Make the interpolated vertex the first
00686                       new_cell->m_smoothed_handle = mesh.add_vertex( coord);
00687                       vertex_iterator v = mesh.get_vertex_iterator(new_cell->m_smoothed_handle);
00688                       v->m_tag = 1;  //--- use tag to remember the  number of vertices that have been smoothed
00689                       cell->m_vertex_handle[axis] = new_cell->m_smoothed_handle;
00690                       vertex_handle[e] = new_cell->m_smoothed_handle;
00691                     }
00692                     else
00693                     {
00694                       //--- we have already started smoothing. Accumulate the interpolated vertex
00695                       vertex_iterator v = mesh.get_vertex_iterator(new_cell->m_smoothed_handle);
00696                       v->m_tag++;               //--- increment the number of smoothed vertices...
00697                       v->m_coord += coord;
00698                       vertex_handle[e] = new_cell->m_smoothed_handle;
00699                       cell->m_vertex_handle[axis] = new_cell->m_smoothed_handle;
00700                     }
00701                   }
00702                   else
00703                   {
00704                     //--- store the index of the new vertex into the cell, and remember it for this edge
00705                     vertex_handle[e] = mesh.add_vertex(coord);
00706                     vertex_iterator v = mesh.get_vertex_iterator(vertex_handle[e]);
00707                     v->m_tag = 1;
00708                     cell->m_vertex_handle[axis] = vertex_handle[e];
00709                   }
00710                 }
00711                 else
00712                 {
00713                   //--- This edge has already been interpolated, so just use it
00714                   vertex_handle[e] = cell->m_vertex_handle[axis];
00715                 }
00716               }
00717               //--- Add a new triangle
00718               if( vertex_handle[0] == vertex_handle[1] || vertex_handle[1] == vertex_handle[2] || vertex_handle[0] == vertex_handle[2] )
00719               {
00720                 //--- It is weird to not have anything in here, but degenerate triangles are discovered faster  due to lazy evaluation.
00721               }
00722               else
00723               {
00724                 mesh.add_face( vertex_handle[0], vertex_handle[1], vertex_handle[2] );
00725               }
00726             }
00727             //--- advance to next pixel/cell
00728             vox0 += sub_sample;
00729             vox1 += sub_sample;
00730             grid[0]++;
00731             grid[1]++;
00732           }
00733           //--- Advance to next row (have to increment once more because loop stops at sampled_I-1)
00734           grid[0]++;
00735           grid[1]++;
00736         }
00737         //--- Swap the two grids as we advance to the next slice
00738         cell_info_type * tmp = cell_grid[0];
00739         cell_grid[0] = cell_grid[1];
00740         cell_grid[1] = tmp;
00741       }
00742 
00743       {
00744         vertex_iterator  v = mesh.vertex_begin();
00745         vertex_iterator  vend = mesh.vertex_end();
00746         for(;v!=vend;++v)
00747           if(v->m_tag>1)
00748           {
00749             real_type N = static_cast<real_type>(1.0/v->m_tag);
00750             v->m_coord *= N;
00751           }
00752       }
00753 
00754       mesh::scale(mesh,sub_sample*phi.dx(),sub_sample*phi.dy(),sub_sample*phi.dz());
00755       mesh::translate(mesh,phi.min_coord());
00756 
00757       mesh::compute_angle_weighted_vertex_normals(mesh);
00758     }
00759 
00760   } // namespace mesh
00761 } // namespace OpenTissue
00762 
00763 //OPENTISSUE_CORE_CONTAINERS_MESH_COMMON_UTIL_MESH_SMOOTH_ISOSURFACE_H
00764 #endif
/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/core/containers/mesh/common/util/mesh_smooth_isosurface.h
