Updated OpenVDB library to 2.3.0;soc-2013-cycles_volume

1 files changed, 3432 insertions, 1443 deletions

diff --git a/extern/openvdb/internal/openvdb/tools/VolumeToMesh.h b/extern/openvdb/internal/openvdb/tools/VolumeToMesh.h
index 78b4d6c8066..ef0f39c7f17 100644
--- a/extern/openvdb/internal/openvdb/tools/VolumeToMesh.h
+++ b/extern/openvdb/internal/openvdb/tools/VolumeToMesh.h
@@ -31,25 +31,32 @@
 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
 
+#include <openvdb/Platform.h> // for OPENVDB_HAS_CXX11
 #include <openvdb/tree/ValueAccessor.h>
 #include <openvdb/util/Util.h> // for COORD_OFFSETS
 #include <openvdb/math/Operators.h> // for ISGradient
 #include <openvdb/tools/Morphology.h> // for dilateVoxels()
+#include <openvdb/tree/LeafManager.h>
 
 #include <boost/scoped_array.hpp>
 #include <boost/scoped_ptr.hpp>
-
 #include <tbb/blocked_range.h>
 #include <tbb/parallel_for.h>
 #include <tbb/parallel_reduce.h>
 
 #include <vector>
+#include <memory> // for auto_ptr/unique_ptr
+
+
+//////////
+
 
 namespace openvdb {
 OPENVDB_USE_VERSION_NAMESPACE
 namespace OPENVDB_VERSION_NAME {
 namespace tools {
 
+
 ////////////////////////////////////////
 
 
@@ -94,71 +101,60 @@ volumeToMesh(
 
 
 /// @brief Polygon flags, used for reference based meshing.
-enum {
-    POLYFLAG_EXTERIOR = 0x1, POLYFLAG_FRACTURE_SEAM = 0x2
-};
+enum { POLYFLAG_EXTERIOR = 0x1, POLYFLAG_FRACTURE_SEAM = 0x2,  POLYFLAG_SUBDIVIDED = 0x4};
 
 
 /// @brief Collection of quads and triangles
 class PolygonPool
 {
 public:
-    PolygonPool()
-    : mNumQuads(0)
-    , mNumTriangles(0)
-    , mQuads(NULL)
-    , mTriangles(NULL)
-    , mQuadFlags(NULL)
-    , mTriangleFlags(NULL)
-    {
-    }
 
-    void resetQuads(size_t size)
-    {
-        mNumQuads = size;
-        mQuads.reset(new openvdb::Vec4I[mNumQuads]);
-        mQuadFlags.reset(new char[mNumQuads]);
-    }
+    inline PolygonPool();
+    inline PolygonPool(const size_t numQuads, const size_t numTriangles);
 
-    void clearQuads()
-    {
-        mNumQuads = 0;
-        mQuads.reset(NULL);
-        mQuadFlags.reset(NULL);
-    }
 
-    void resetTriangles(size_t size)
-    {
-        mNumTriangles = size;
-        mTriangles.reset(new openvdb::Vec3I[mNumTriangles]);
-        mTriangleFlags.reset(new char[mNumTriangles]);
-    }
+    inline void copy(const PolygonPool& rhs);
 
-    void clearTriangles()
-    {
-        mNumTriangles = 0;
-        mTriangles.reset(NULL);
-        mTriangleFlags.reset(NULL);
-    }
+    inline void resetQuads(size_t size);
+    inline void clearQuads();
+
+    inline void resetTriangles(size_t size);
+    inline void clearTriangles();
 
-    const size_t& numQuads() const { return mNumQuads; }
-    const size_t& numTriangles() const { return mNumTriangles; }
 
     // polygon accessor methods
-    openvdb::Vec4I& quad(size_t n) { return mQuads[n]; }
-    const openvdb::Vec4I& quad(size_t n) const { return mQuads[n]; }
 
-    openvdb::Vec3I& triangle(size_t n) { return mTriangles[n]; }
-    const openvdb::Vec3I& triangle(size_t n) const { return mTriangles[n]; }
+    const size_t& numQuads() const                      { return mNumQuads; }
+
+    openvdb::Vec4I& quad(size_t n)                      { return mQuads[n]; }
+    const openvdb::Vec4I& quad(size_t n) const          { return mQuads[n]; }
+
+
+    const size_t& numTriangles() const                  { return mNumTriangles; }
+
+    openvdb::Vec3I& triangle(size_t n)                  { return mTriangles[n]; }
+    const openvdb::Vec3I& triangle(size_t n) const      { return mTriangles[n]; }
+
 
     // polygon flags accessor methods
-    char& quadFlags(size_t n) { return mQuadFlags[n]; }
-    const char& quadFlags(size_t n) const { return mQuadFlags[n]; }
 
-    char& triangleFlags(size_t n) { return mTriangleFlags[n]; }
-    const char& triangleFlags(size_t n) const { return mTriangleFlags[n]; }
+    char& quadFlags(size_t n)                           { return mQuadFlags[n]; }
+    const char& quadFlags(size_t n) const               { return mQuadFlags[n]; }
+
+    char& triangleFlags(size_t n)                       { return mTriangleFlags[n]; }
+    const char& triangleFlags(size_t n) const           { return mTriangleFlags[n]; }
+
+
+    // reduce the polygon containers, n has to
+    // be smaller than the current container size.
+
+    inline bool trimQuads(const size_t n, bool reallocate = false);
+    inline bool trimTrinagles(const size_t n, bool reallocate = false);
 
 private:
+    // disallow copy by assignment
+    void operator=(const PolygonPool&) {}
+
     size_t mNumQuads, mNumTriangles;
     boost::scoped_array<openvdb::Vec4I> mQuads;
     boost::scoped_array<openvdb::Vec3I> mTriangles;
@@ -180,23 +176,39 @@ typedef boost::scoped_array<PolygonPool> PolygonPoolList;
 class VolumeToMesh
 {
 public:
+
     /// @param isovalue         Determines which isosurface to mesh.
     /// @param adaptivity       Adaptivity threshold [0 to 1]
     VolumeToMesh(double isovalue = 0, double adaptivity = 0);
 
-    PointList& pointList();
+
+    //////////
+
+    // Mesh data accessors
+
     const size_t& pointListSize() const;
+    PointList& pointList();
 
+    const size_t& polygonPoolListSize() const;
     PolygonPoolList& polygonPoolList();
     const PolygonPoolList& polygonPoolList() const;
-    const size_t& polygonPoolListSize() const;
 
-    /// @brief main call
+    std::vector<unsigned char>& pointFlags();
+    const std::vector<unsigned char>& pointFlags() const;
+
+
+    //////////
+
+
+    /// @brief Main call
     /// @note Call with scalar typed grid.
     template<typename GridT>
     void operator()(const GridT&);
 
 
+    //////////
+
+
     /// @brief  When surfacing fractured SDF fragments, the original unfractured
     ///         SDF grid can be used to eliminate seam lines and tag polygons that are
     ///         coincident with the reference surface with the @c POLYFLAG_EXTERIOR
@@ -218,93 +230,452 @@ public:
     ///                         that do not exist in the reference grid. (Parts of the
     ///                         fragment surface that are not coincident with the
     ///                         reference surface.)
-    /// @param  smoothSeams     toggle to smooth seam line edges during mesh extraction,
-    ///                         removes staircase artifacts.
-    void setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity = 0, bool smoothSeams = true);
+    void setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity = 0);
+
+
+    /// @param mask A boolean grid whose active topology defines the region to mesh.
+    /// @param invertMask Toggle to mesh the complement of the mask.
+    /// @note The mask's tree configuration has to match @c GridT's tree configuration.
+    void setSurfaceMask(const GridBase::ConstPtr& mask, bool invertMask = false);
+
+    /// @param grid A scalar grid used as an spatial multiplier for the adaptivity threshold.
+    /// @note The grid's tree configuration has to match @c GridT's tree configuration.
+    void setSpatialAdaptivity(const GridBase::ConstPtr& grid);
+
+
+    /// @param tree A boolean tree whose active topology defines the adaptivity mask.
+    /// @note The tree configuration has to match @c GridT's tree configuration.
+    void setAdaptivityMask(const TreeBase::ConstPtr& tree);
+
+
+    /// @brief Subdivide volume and mesh into disjoint parts
+    /// @param partitions Number of partitions.
+    /// @param activePart Specific partition to mesh, 0 to @c partitions - 1.
+    void partition(unsigned partitions = 1, unsigned activePart = 0);
 
 private:
 
     PointList mPoints;
     PolygonPoolList mPolygons;
 
-    size_t mPointListSize, mPolygonPoolListSize;
+    size_t mPointListSize, mSeamPointListSize, mPolygonPoolListSize;
     double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
 
-    GridBase::ConstPtr mRefGrid;
-    TreeBase::Ptr mRefEdgeTree, mRefTopologyMaskTree, mSeamPointTree;
-    bool mSmoothSeams;
+    GridBase::ConstPtr mRefGrid, mSurfaceMaskGrid, mAdaptivityGrid;
+    TreeBase::ConstPtr mAdaptivityMaskTree;
+
+    TreeBase::Ptr mRefSignTree, mRefIdxTree;
+
+    bool mInvertSurfaceMask;
+    unsigned mPartitions, mActivePart;
+
+    boost::scoped_array<uint32_t> mQuantizedSeamPoints;
+
+    std::vector<unsigned char> mPointFlags;
 };
 
 
 ////////////////////////////////////////
 
 
-// Internal utility methods
+/// @brief  Given a set of tangent elements, @c points with corresponding @c normals,
+///         this method returns the intersection point of all tangent elements.
+///
+/// @note   Used to extract surfaces with sharp edges and corners from volume data,
+///         see the following paper for details: "Feature Sensitive Surface
+///         Extraction from Volume Data, Kobbelt et al. 2001".
+inline Vec3d findFeaturePoint(
+    const std::vector<Vec3d>& points,
+    const std::vector<Vec3d>& normals)
+{
+    typedef math::Mat3d Mat3d;
+
+    Vec3d avgPos(0.0);
+
+    if (points.empty()) return avgPos;
+
+    for (size_t n = 0, N = points.size(); n < N; ++n) {
+        avgPos += points[n];
+    }
+
+    avgPos /= double(points.size());
+
+    // Unique components of the 3x3 A^TA matrix, where A is
+    // the matrix of normals.
+    double m00=0,m01=0,m02=0,
+           m11=0,m12=0,
+           m22=0;
+
+    // The rhs vector, A^Tb, where b = n dot p
+    Vec3d rhs(0.0);
+
+    for (size_t n = 0, N = points.size(); n < N; ++n) {
+
+        const Vec3d& n_ref = normals[n];
+
+        // A^TA
+        m00 += n_ref[0] * n_ref[0]; // diagonal
+        m11 += n_ref[1] * n_ref[1];
+        m22 += n_ref[2] * n_ref[2];
+
+        m01 += n_ref[0] * n_ref[1]; // Upper-tri
+        m02 += n_ref[0] * n_ref[2];
+        m12 += n_ref[1] * n_ref[2];
+
+        // A^Tb (centered around the origin)
+        rhs += n_ref * n_ref.dot(points[n] - avgPos);
+    }
+
+    Mat3d A(m00,m01,m02,
+            m01,m11,m12,
+            m02,m12,m22);
 
+    /*
+    // Inverse
+    const double det = A.det();
+    if (det > 0.01) {
+        Mat3d A_inv = A.adjoint();
+        A_inv *= (1.0 / det);
 
+        return avgPos + A_inv * rhs;
+    }
+    */
+
+    // Compute the pseudo inverse
+
+    math::Mat3d eigenVectors;
+    Vec3d eigenValues;
+
+    diagonalizeSymmetricMatrix(A, eigenVectors, eigenValues, 300);
+
+    Mat3d D = Mat3d::identity();
+
+
+    double tolerance = std::max(std::abs(eigenValues[0]), std::abs(eigenValues[1]));
+    tolerance = std::max(tolerance, std::abs(eigenValues[2]));
+    tolerance *= 0.01;
+
+    int clamped = 0;
+    for (int i = 0; i < 3; ++i ) {
+        if (std::abs(eigenValues[i]) < tolerance) {
+            D[i][i] = 0.0;
+            ++clamped;
+        } else {
+            D[i][i] = 1.0 / eigenValues[i];
+        }
+    }
+
+    // Assemble the pseudo inverse and calc. the intersection point
+    if (clamped < 3) {
+        Mat3d pseudoInv = eigenVectors * D *  eigenVectors.transpose();
+        return avgPos + pseudoInv * rhs;
+    }
+
+    return avgPos;
+}
+
+
+////////////////////////////////////////
+
+
+// Internal utility methods
 namespace internal {
 
+template<typename T>
+struct UniquePtr
+{
+#ifdef OPENVDB_HAS_CXX11
+    typedef std::unique_ptr<T>  type;
+#else
+    typedef std::auto_ptr<T>    type;
+#endif
+};
+
 
-// Bit-flags
-enum { INSIDE = 0x1, XEDGE = 0x2, YEDGE = 0x4, ZEDGE = 0x8 };
-
-const bool sAmbiguous[256] =
-   {0,0,0,0,0,1,0,0,
-    0,0,1,0,0,0,0,0,
-    0,0,1,0,1,1,1,0,
-    1,0,1,0,1,0,1,0,
-    0,1,0,0,1,1,0,0,
-    1,1,1,0,1,1,0,0,
-    0,0,0,0,1,1,0,0,
-    1,0,1,0,1,1,1,0,
-    0,1,1,1,0,1,0,0,
-    1,1,1,1,0,0,0,0,
-    1,1,1,1,1,1,1,1,
-    1,1,1,1,1,1,1,1,
-    0,1,0,0,0,1,0,0,
-    1,1,1,1,0,1,0,0,
-    0,0,0,0,0,1,0,0,
-    1,1,1,1,1,1,1,0,
-    0,1,1,1,1,1,1,1,
-    0,0,1,0,0,0,0,0,
-    0,0,1,0,1,1,1,1,
-    0,0,1,0,0,0,1,0,
-    1,1,1,1,1,1,1,1,
-    1,1,1,1,1,1,1,1,
-    0,0,0,0,1,1,1,1,
-    0,0,1,0,1,1,1,0,
-    0,1,1,1,0,1,0,1,
-    0,0,1,1,0,0,0,0,
-    0,0,1,1,0,1,1,1,
-    0,0,1,1,0,0,1,0,
-    0,1,0,1,0,1,0,1,
-    0,1,1,1,0,1,0,0,
-    0,0,0,0,0,1,0,0,
-    0,0,1,0,0,0,0,0};
+/// @brief  Bit-flags used to classify cells.
+enum { SIGNS = 0xFF, EDGES = 0xE00, INSIDE = 0x100,
+       XEDGE = 0x200, YEDGE = 0x400, ZEDGE = 0x800, SEAM = 0x1000};
+
+
+/// @brief Used to quickly determine if a given cell is adaptable.
+const bool sAdaptable[256] = {
+    1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,0,1,0,0,0,1,0,1,0,1,0,1,0,1,
+    1,0,1,1,0,0,1,1,0,0,0,1,0,0,1,1,1,1,1,1,0,0,1,1,0,1,0,1,0,0,0,1,
+    1,0,0,0,1,0,1,1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+    1,0,1,1,1,0,1,1,0,0,0,0,1,0,1,1,1,1,1,1,1,0,1,1,0,0,0,0,0,0,0,1,
+    1,0,0,0,0,0,0,0,1,1,0,1,1,1,1,1,1,1,0,1,0,0,0,0,1,1,0,1,1,1,0,1,
+    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,1,1,0,1,0,0,0,1,
+    1,0,0,0,1,0,1,0,1,1,0,0,1,1,1,1,1,1,0,0,1,0,0,0,1,1,0,0,1,1,0,1,
+    1,0,1,0,1,0,1,0,1,0,0,0,1,0,1,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1};
+
+
+/// @brief  Contains the ambiguous face index for certain cell configuration.
+const unsigned char sAmbiguousFace[256] = {
+    0,0,0,0,0,5,0,0,0,0,5,0,0,0,0,0,0,0,1,0,0,5,1,0,4,0,0,0,4,0,0,0,
+    0,1,0,0,2,0,0,0,0,1,5,0,2,0,0,0,0,0,0,0,2,0,0,0,4,0,0,0,0,0,0,0,
+    0,0,2,2,0,5,0,0,3,3,0,0,0,0,0,0,6,6,0,0,6,0,0,0,0,0,0,0,0,0,0,0,
+    0,1,0,0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+    0,4,0,4,3,0,3,0,0,0,5,0,0,0,0,0,0,0,1,0,3,0,0,0,0,0,0,0,0,0,0,0,
+    6,0,6,0,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+    0,4,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+
+
+/// @brief  Lookup table for different cell sign configurations. The first entry specifies
+///         the total number of points that need to be generated inside a cell and the
+///         remaining 12 entries indicate different edge groups.
+const unsigned char sEdgeGroupTable[256][13] = {
+    {0,0,0,0,0,0,0,0,0,0,0,0,0},{1,1,0,0,1,0,0,0,0,1,0,0,0},{1,1,1,0,0,0,0,0,0,0,1,0,0},
+    {1,0,1,0,1,0,0,0,0,1,1,0,0},{1,0,1,1,0,0,0,0,0,0,0,1,0},{1,1,1,1,1,0,0,0,0,1,0,1,0},
+    {1,1,0,1,0,0,0,0,0,0,1,1,0},{1,0,0,1,1,0,0,0,0,1,1,1,0},{1,0,0,1,1,0,0,0,0,0,0,0,1},
+    {1,1,0,1,0,0,0,0,0,1,0,0,1},{1,1,1,1,1,0,0,0,0,0,1,0,1},{1,0,1,1,0,0,0,0,0,1,1,0,1},
+    {1,0,1,0,1,0,0,0,0,0,0,1,1},{1,1,1,0,0,0,0,0,0,1,0,1,1},{1,1,0,0,1,0,0,0,0,0,1,1,1},
+    {1,0,0,0,0,0,0,0,0,1,1,1,1},{1,0,0,0,0,1,0,0,1,1,0,0,0},{1,1,0,0,1,1,0,0,1,0,0,0,0},
+    {1,1,1,0,0,1,0,0,1,1,1,0,0},{1,0,1,0,1,1,0,0,1,0,1,0,0},{2,0,1,1,0,2,0,0,2,2,0,1,0},
+    {1,1,1,1,1,1,0,0,1,0,0,1,0},{1,1,0,1,0,1,0,0,1,1,1,1,0},{1,0,0,1,1,1,0,0,1,0,1,1,0},
+    {1,0,0,1,1,1,0,0,1,1,0,0,1},{1,1,0,1,0,1,0,0,1,0,0,0,1},{2,2,1,1,2,1,0,0,1,2,1,0,1},
+    {1,0,1,1,0,1,0,0,1,0,1,0,1},{1,0,1,0,1,1,0,0,1,1,0,1,1},{1,1,1,0,0,1,0,0,1,0,0,1,1},
+    {2,1,0,0,1,2,0,0,2,1,2,2,2},{1,0,0,0,0,1,0,0,1,0,1,1,1},{1,0,0,0,0,1,1,0,0,0,1,0,0},
+    {1,1,0,0,1,1,1,0,0,1,1,0,0},{1,1,1,0,0,1,1,0,0,0,0,0,0},{1,0,1,0,1,1,1,0,0,1,0,0,0},
+    {1,0,1,1,0,1,1,0,0,0,1,1,0},{2,2,2,1,1,1,1,0,0,1,2,1,0},{1,1,0,1,0,1,1,0,0,0,0,1,0},
+    {1,0,0,1,1,1,1,0,0,1,0,1,0},{2,0,0,2,2,1,1,0,0,0,1,0,2},{1,1,0,1,0,1,1,0,0,1,1,0,1},
+    {1,1,1,1,1,1,1,0,0,0,0,0,1},{1,0,1,1,0,1,1,0,0,1,0,0,1},{1,0,1,0,1,1,1,0,0,0,1,1,1},
+    {2,1,1,0,0,2,2,0,0,2,1,2,2},{1,1,0,0,1,1,1,0,0,0,0,1,1},{1,0,0,0,0,1,1,0,0,1,0,1,1},
+    {1,0,0,0,0,0,1,0,1,1,1,0,0},{1,1,0,0,1,0,1,0,1,0,1,0,0},{1,1,1,0,0,0,1,0,1,1,0,0,0},
+    {1,0,1,0,1,0,1,0,1,0,0,0,0},{1,0,1,1,0,0,1,0,1,1,1,1,0},{2,1,1,2,2,0,2,0,2,0,1,2,0},
+    {1,1,0,1,0,0,1,0,1,1,0,1,0},{1,0,0,1,1,0,1,0,1,0,0,1,0},{1,0,0,1,1,0,1,0,1,1,1,0,1},
+    {1,1,0,1,0,0,1,0,1,0,1,0,1},{2,1,2,2,1,0,2,0,2,1,0,0,2},{1,0,1,1,0,0,1,0,1,0,0,0,1},
+    {2,0,2,0,2,0,1,0,1,2,2,1,1},{2,2,2,0,0,0,1,0,1,0,2,1,1},{2,2,0,0,2,0,1,0,1,2,0,1,1},
+    {1,0,0,0,0,0,1,0,1,0,0,1,1},{1,0,0,0,0,0,1,1,0,0,0,1,0},{2,1,0,0,1,0,2,2,0,1,0,2,0},
+    {1,1,1,0,0,0,1,1,0,0,1,1,0},{1,0,1,0,1,0,1,1,0,1,1,1,0},{1,0,1,1,0,0,1,1,0,0,0,0,0},
+    {1,1,1,1,1,0,1,1,0,1,0,0,0},{1,1,0,1,0,0,1,1,0,0,1,0,0},{1,0,0,1,1,0,1,1,0,1,1,0,0},
+    {1,0,0,1,1,0,1,1,0,0,0,1,1},{1,1,0,1,0,0,1,1,0,1,0,1,1},{2,1,2,2,1,0,1,1,0,0,1,2,1},
+    {2,0,1,1,0,0,2,2,0,2,2,1,2},{1,0,1,0,1,0,1,1,0,0,0,0,1},{1,1,1,0,0,0,1,1,0,1,0,0,1},
+    {1,1,0,0,1,0,1,1,0,0,1,0,1},{1,0,0,0,0,0,1,1,0,1,1,0,1},{1,0,0,0,0,1,1,1,1,1,0,1,0},
+    {1,1,0,0,1,1,1,1,1,0,0,1,0},{2,1,1,0,0,2,2,1,1,1,2,1,0},{2,0,2,0,2,1,1,2,2,0,1,2,0},
+    {1,0,1,1,0,1,1,1,1,1,0,0,0},{2,2,2,1,1,2,2,1,1,0,0,0,0},{2,2,0,2,0,1,1,2,2,2,1,0,0},
+    {2,0,0,1,1,2,2,1,1,0,2,0,0},{2,0,0,1,1,1,1,2,2,1,0,1,2},{2,2,0,2,0,2,2,1,1,0,0,2,1},
+    {4,3,2,2,3,4,4,1,1,3,4,2,1},{3,0,2,2,0,1,1,3,3,0,1,2,3},{2,0,2,0,2,2,2,1,1,2,0,0,1},
+    {2,1,1,0,0,1,1,2,2,0,0,0,2},{3,1,0,0,1,2,2,3,3,1,2,0,3},{2,0,0,0,0,1,1,2,2,0,1,0,2},
+    {1,0,0,0,0,1,0,1,0,0,1,1,0},{1,1,0,0,1,1,0,1,0,1,1,1,0},{1,1,1,0,0,1,0,1,0,0,0,1,0},
+    {1,0,1,0,1,1,0,1,0,1,0,1,0},{1,0,1,1,0,1,0,1,0,0,1,0,0},{2,1,1,2,2,2,0,2,0,2,1,0,0},
+    {1,1,0,1,0,1,0,1,0,0,0,0,0},{1,0,0,1,1,1,0,1,0,1,0,0,0},{1,0,0,1,1,1,0,1,0,0,1,1,1},
+    {2,2,0,2,0,1,0,1,0,1,2,2,1},{2,2,1,1,2,2,0,2,0,0,0,1,2},{2,0,2,2,0,1,0,1,0,1,0,2,1},
+    {1,0,1,0,1,1,0,1,0,0,1,0,1},{2,2,2,0,0,1,0,1,0,1,2,0,1},{1,1,0,0,1,1,0,1,0,0,0,0,1},
+    {1,0,0,0,0,1,0,1,0,1,0,0,1},{1,0,0,0,0,0,0,1,1,1,1,1,0},{1,1,0,0,1,0,0,1,1,0,1,1,0},
+    {1,1,1,0,0,0,0,1,1,1,0,1,0},{1,0,1,0,1,0,0,1,1,0,0,1,0},{1,0,1,1,0,0,0,1,1,1,1,0,0},
+    {2,2,2,1,1,0,0,1,1,0,2,0,0},{1,1,0,1,0,0,0,1,1,1,0,0,0},{1,0,0,1,1,0,0,1,1,0,0,0,0},
+    {2,0,0,2,2,0,0,1,1,2,2,2,1},{2,1,0,1,0,0,0,2,2,0,1,1,2},{3,2,1,1,2,0,0,3,3,2,0,1,3},
+    {2,0,1,1,0,0,0,2,2,0,0,1,2},{2,0,1,0,1,0,0,2,2,1,1,0,2},{2,1,1,0,0,0,0,2,2,0,1,0,2},
+    {2,1,0,0,1,0,0,2,2,1,0,0,2},{1,0,0,0,0,0,0,1,1,0,0,0,1},{1,0,0,0,0,0,0,1,1,0,0,0,1},
+    {1,1,0,0,1,0,0,1,1,1,0,0,1},{2,1,1,0,0,0,0,2,2,0,1,0,2},{1,0,1,0,1,0,0,1,1,1,1,0,1},
+    {1,0,1,1,0,0,0,1,1,0,0,1,1},{2,1,1,2,2,0,0,1,1,1,0,1,2},{1,1,0,1,0,0,0,1,1,0,1,1,1},
+    {2,0,0,1,1,0,0,2,2,2,2,2,1},{1,0,0,1,1,0,0,1,1,0,0,0,0},{1,1,0,1,0,0,0,1,1,1,0,0,0},
+    {1,1,1,1,1,0,0,1,1,0,1,0,0},{1,0,1,1,0,0,0,1,1,1,1,0,0},{1,0,1,0,1,0,0,1,1,0,0,1,0},
+    {1,1,1,0,0,0,0,1,1,1,0,1,0},{1,1,0,0,1,0,0,1,1,0,1,1,0},{1,0,0,0,0,0,0,1,1,1,1,1,0},
+    {1,0,0,0,0,1,0,1,0,1,0,0,1},{1,1,0,0,1,1,0,1,0,0,0,0,1},{1,1,1,0,0,1,0,1,0,1,1,0,1},
+    {1,0,1,0,1,1,0,1,0,0,1,0,1},{1,0,1,1,0,1,0,1,0,1,0,1,1},{2,2,2,1,1,2,0,2,0,0,0,2,1},
+    {2,1,0,1,0,2,0,2,0,1,2,2,1},{2,0,0,2,2,1,0,1,0,0,1,1,2},{1,0,0,1,1,1,0,1,0,1,0,0,0},
+    {1,1,0,1,0,1,0,1,0,0,0,0,0},{2,1,2,2,1,2,0,2,0,1,2,0,0},{1,0,1,1,0,1,0,1,0,0,1,0,0},
+    {1,0,1,0,1,1,0,1,0,1,0,1,0},{1,1,1,0,0,1,0,1,0,0,0,1,0},{2,2,0,0,2,1,0,1,0,2,1,1,0},
+    {1,0,0,0,0,1,0,1,0,0,1,1,0},{1,0,0,0,0,1,1,1,1,0,1,0,1},{2,1,0,0,1,2,1,1,2,2,1,0,1},
+    {1,1,1,0,0,1,1,1,1,0,0,0,1},{2,0,2,0,2,1,2,2,1,1,0,0,2},{2,0,1,1,0,1,2,2,1,0,1,2,1},
+    {4,1,1,3,3,2,4,4,2,2,1,4,3},{2,2,0,2,0,2,1,1,2,0,0,1,2},{3,0,0,1,1,2,3,3,2,2,0,3,1},
+    {1,0,0,1,1,1,1,1,1,0,1,0,0},{2,2,0,2,0,1,2,2,1,1,2,0,0},{2,2,1,1,2,2,1,1,2,0,0,0,0},
+    {2,0,1,1,0,2,1,1,2,2,0,0,0},{2,0,2,0,2,2,1,1,2,0,2,1,0},{3,1,1,0,0,3,2,2,3,3,1,2,0},
+    {2,1,0,0,1,1,2,2,1,0,0,2,0},{2,0,0,0,0,2,1,1,2,2,0,1,0},{1,0,0,0,0,0,1,1,0,1,1,0,1},
+    {1,1,0,0,1,0,1,1,0,0,1,0,1},{1,1,1,0,0,0,1,1,0,1,0,0,1},{1,0,1,0,1,0,1,1,0,0,0,0,1},
+    {2,0,2,2,0,0,1,1,0,2,2,1,2},{3,1,1,2,2,0,3,3,0,0,1,3,2},{2,1,0,1,0,0,2,2,0,1,0,2,1},
+    {2,0,0,1,1,0,2,2,0,0,0,2,1},{1,0,0,1,1,0,1,1,0,1,1,0,0},{1,1,0,1,0,0,1,1,0,0,1,0,0},
+    {2,2,1,1,2,0,1,1,0,2,0,0,0},{1,0,1,1,0,0,1,1,0,0,0,0,0},{2,0,1,0,1,0,2,2,0,1,1,2,0},
+    {2,1,1,0,0,0,2,2,0,0,1,2,0},{2,1,0,0,1,0,2,2,0,1,0,2,0},{1,0,0,0,0,0,1,1,0,0,0,1,0},
+    {1,0,0,0,0,0,1,0,1,0,0,1,1},{1,1,0,0,1,0,1,0,1,1,0,1,1},{1,1,1,0,0,0,1,0,1,0,1,1,1},
+    {2,0,2,0,2,0,1,0,1,1,1,2,2},{1,0,1,1,0,0,1,0,1,0,0,0,1},{2,2,2,1,1,0,2,0,2,2,0,0,1},
+    {1,1,0,1,0,0,1,0,1,0,1,0,1},{2,0,0,2,2,0,1,0,1,1,1,0,2},{1,0,0,1,1,0,1,0,1,0,0,1,0},
+    {1,1,0,1,0,0,1,0,1,1,0,1,0},{2,2,1,1,2,0,2,0,2,0,2,1,0},{2,0,2,2,0,0,1,0,1,1,1,2,0},
+    {1,0,1,0,1,0,1,0,1,0,0,0,0},{1,1,1,0,0,0,1,0,1,1,0,0,0},{1,1,0,0,1,0,1,0,1,0,1,0,0},
+    {1,0,0,0,0,0,1,0,1,1,1,0,0},{1,0,0,0,0,1,1,0,0,1,0,1,1},{1,1,0,0,1,1,1,0,0,0,0,1,1},
+    {2,2,2,0,0,1,1,0,0,2,1,2,2},{2,0,1,0,1,2,2,0,0,0,2,1,1},{1,0,1,1,0,1,1,0,0,1,0,0,1},
+    {2,1,1,2,2,1,1,0,0,0,0,0,2},{2,1,0,1,0,2,2,0,0,1,2,0,1},{2,0,0,2,2,1,1,0,0,0,1,0,2},
+    {1,0,0,1,1,1,1,0,0,1,0,1,0},{1,1,0,1,0,1,1,0,0,0,0,1,0},{3,1,2,2,1,3,3,0,0,1,3,2,0},
+    {2,0,1,1,0,2,2,0,0,0,2,1,0},{1,0,1,0,1,1,1,0,0,1,0,0,0},{1,1,1,0,0,1,1,0,0,0,0,0,0},
+    {2,2,0,0,2,1,1,0,0,2,1,0,0},{1,0,0,0,0,1,1,0,0,0,1,0,0},{1,0,0,0,0,1,0,0,1,0,1,1,1},
+    {2,2,0,0,2,1,0,0,1,1,2,2,2},{1,1,1,0,0,1,0,0,1,0,0,1,1},{2,0,1,0,1,2,0,0,2,2,0,1,1},
+    {1,0,1,1,0,1,0,0,1,0,1,0,1},{3,1,1,3,3,2,0,0,2,2,1,0,3},{1,1,0,1,0,1,0,0,1,0,0,0,1},
+    {2,0,0,2,2,1,0,0,1,1,0,0,2},{1,0,0,1,1,1,0,0,1,0,1,1,0},{2,1,0,1,0,2,0,0,2,2,1,1,0},
+    {2,1,2,2,1,1,0,0,1,0,0,2,0},{2,0,1,1,0,2,0,0,2,2,0,1,0},{1,0,1,0,1,1,0,0,1,0,1,0,0},
+    {2,1,1,0,0,2,0,0,2,2,1,0,0},{1,1,0,0,1,1,0,0,1,0,0,0,0},{1,0,0,0,0,1,0,0,1,1,0,0,0},
+    {1,0,0,0,0,0,0,0,0,1,1,1,1},{1,1,0,0,1,0,0,0,0,0,1,1,1},{1,1,1,0,0,0,0,0,0,1,0,1,1},
+    {1,0,1,0,1,0,0,0,0,0,0,1,1},{1,0,1,1,0,0,0,0,0,1,1,0,1},{2,1,1,2,2,0,0,0,0,0,1,0,2},
+    {1,1,0,1,0,0,0,0,0,1,0,0,1},{1,0,0,1,1,0,0,0,0,0,0,0,1},{1,0,0,1,1,0,0,0,0,1,1,1,0},
+    {1,1,0,1,0,0,0,0,0,0,1,1,0},{2,1,2,2,1,0,0,0,0,1,0,2,0},{1,0,1,1,0,0,0,0,0,0,0,1,0},
+    {1,0,1,0,1,0,0,0,0,1,1,0,0},{1,1,1,0,0,0,0,0,0,0,1,0,0},{1,1,0,0,1,0,0,0,0,1,0,0,0},
+    {0,0,0,0,0,0,0,0,0,0,0,0,0}};
 
 
-template<class AccessorT>
-inline bool isAmbiguous(const AccessorT& accessor, const Coord& ijk,
-    typename AccessorT::ValueType isovalue, int dim)
+////////////////////////////////////////
+
+inline bool
+isPlanarQuad(
+    const Vec3d& p0, const Vec3d& p1,
+    const Vec3d& p2, const Vec3d& p3,
+    double epsilon = 0.001)
+{
+    // compute representative plane
+    Vec3d normal = (p2-p0).cross(p1-p3);
+    normal.normalize();
+    const Vec3d centroid = (p0 + p1 + p2 + p3);
+    const double d = centroid.dot(normal) * 0.25;
+
+
+    // test vertice distance to plane
+    double absDist = std::abs(p0.dot(normal) - d);
+    if (absDist > epsilon) return false;
+
+    absDist = std::abs(p1.dot(normal) - d);
+    if (absDist > epsilon) return false;
+
+    absDist = std::abs(p2.dot(normal) - d);
+    if (absDist > epsilon) return false;
+
+    absDist = std::abs(p3.dot(normal) - d);
+    if (absDist > epsilon) return false;
+
+    return true;
+}
+
+
+////////////////////////////////////////
+
+
+/// @{
+/// @brief  Utility methods for point quantization.
+
+enum { MASK_FIRST_10_BITS = 0x000003FF, MASK_DIRTY_BIT = 0x80000000, MASK_INVALID_BIT = 0x40000000 };
+
+inline uint32_t
+packPoint(const Vec3d& v)
+{
+    uint32_t data = 0;
+
+    // values are expected to be in the [0.0 to 1.0] range.
+    assert(!(v.x() > 1.0) && !(v.y() > 1.0) && !(v.z() > 1.0));
+    assert(!(v.x() < 0.0) && !(v.y() < 0.0) && !(v.z() < 0.0));
+
+    data |= (uint32_t(v.x() * 1023.0) & MASK_FIRST_10_BITS) << 20;
+    data |= (uint32_t(v.y() * 1023.0) & MASK_FIRST_10_BITS) << 10;
+    data |= (uint32_t(v.z() * 1023.0) & MASK_FIRST_10_BITS);
+
+    return data;
+}
+
+inline Vec3d
+unpackPoint(uint32_t data)
+{
+    Vec3d v;
+    v.z() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
+    data = data >> 10;
+    v.y() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
+    data = data >> 10;
+    v.x() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
+
+    return v;
+}
+
+/// @}
+
+////////////////////////////////////////
+
+
+/// @brief  General method that computes the cell-sign configuration at the given
+///         @c ijk coordinate.
+template<typename AccessorT>
+inline unsigned char
+evalCellSigns(const AccessorT& accessor, const Coord& ijk, typename AccessorT::ValueType iso)
 {
     unsigned signs = 0;
     Coord coord = ijk; // i, j, k
-    if (accessor.getValue(coord) < isovalue) signs |= 1u;
-    coord[0] += dim; // i+dim, j, k
-    if (accessor.getValue(coord) < isovalue) signs |= 2u;
-    coord[2] += dim; // i+dim, j, k+dim
-    if (accessor.getValue(coord) < isovalue) signs |= 4u;
-    coord[0] = ijk[0]; // i, j, k+dim
-    if (accessor.getValue(coord) < isovalue) signs |= 8u;
-    coord[1] += dim; coord[2] = ijk[2]; // i, j+dim, k
-    if (accessor.getValue(coord) < isovalue) signs |= 16u;
-    coord[0] += dim; // i+dim, j+dim, k
-    if (accessor.getValue(coord) < isovalue) signs |= 32u;
-    coord[2] += dim; // i+dim, j+dim, k+dim
-    if (accessor.getValue(coord) < isovalue) signs |= 64u;
-    coord[0] = ijk[0]; // i, j+dim, k+dim
-    if (accessor.getValue(coord) < isovalue) signs |= 128u;
-    return sAmbiguous[signs];
+    if (accessor.getValue(coord) < iso) signs |= 1u;
+    coord[0] += 1; // i+1, j, k
+    if (accessor.getValue(coord) < iso) signs |= 2u;
+    coord[2] += 1; // i+1, j, k+1
+    if (accessor.getValue(coord) < iso) signs |= 4u;
+    coord[0] = ijk[0]; // i, j, k+1
+    if (accessor.getValue(coord) < iso) signs |= 8u;
+    coord[1] += 1; coord[2] = ijk[2]; // i, j+1, k
+    if (accessor.getValue(coord) < iso) signs |= 16u;
+    coord[0] += 1; // i+1, j+1, k
+    if (accessor.getValue(coord) < iso) signs |= 32u;
+    coord[2] += 1; // i+1, j+1, k+1
+    if (accessor.getValue(coord) < iso) signs |= 64u;
+    coord[0] = ijk[0]; // i, j+1, k+1
+    if (accessor.getValue(coord) < iso) signs |= 128u;
+    return signs;
+}
+
+
+/// @brief  Leaf node optimized method that computes the cell-sign configuration
+///         at the given local @c offset
+template<typename LeafT>
+inline unsigned char
+evalCellSigns(const LeafT& leaf, const Index offset, typename LeafT::ValueType iso)
+{
+    unsigned char signs = 0;
+
+    // i, j, k
+    if (leaf.getValue(offset) < iso) signs |= 1u;
+
+    // i, j, k+1
+    if (leaf.getValue(offset + 1) < iso) signs |= 8u;
+
+    // i, j+1, k
+    if (leaf.getValue(offset + LeafT::DIM) < iso) signs |= 16u;
+
+    // i, j+1, k+1
+    if (leaf.getValue(offset + LeafT::DIM + 1) < iso) signs |= 128u;
+
+    // i+1, j, k
+    if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) ) < iso) signs |= 2u;
+
+    // i+1, j, k+1
+    if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1) < iso) signs |= 4u;
+
+    // i+1, j+1, k
+    if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM) < iso) signs |= 32u;
+
+    // i+1, j+1, k+1
+    if (leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1) < iso) signs |= 64u;
+
+    return signs;
+}
+
+
+/// @brief  Used to correct topological ambiguities related to two adjacent cells
+///         that share an ambiguous face.
+template<class AccessorT>
+inline void
+correctCellSigns(unsigned char& signs, unsigned char face,
+    const AccessorT& acc, Coord ijk, typename AccessorT::ValueType iso)
+{
+    if (face == 1) {
+        ijk[2] -= 1;
+        if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 3) signs = ~signs;
+    } else if (face == 3) {
+        ijk[2] += 1;
+        if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 1) signs = ~signs;
+    } else if (face == 2) {
+        ijk[0] += 1;
+        if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 4) signs = ~signs;
+    } else if (face == 4) {
+        ijk[0] -= 1;
+        if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 2) signs = ~signs;
+    } else if (face == 5) {
+        ijk[1] -= 1;
+        if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 6) signs = ~signs;
+    } else if (face == 6) {
+        ijk[1] += 1;
+        if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 5) signs = ~signs;
+    }
 }
 
 
@@ -343,7 +714,7 @@ isNonManifold(const AccessorT& accessor, const Coord& ijk,
     if (p[5]) signs |= 32u;
     if (p[6]) signs |= 64u;
     if (p[7]) signs |= 128u;
-    if (sAmbiguous[signs]) return true;
+    if (!sAdaptable[signs]) return true;
 
     // Manifold check
 
@@ -519,420 +890,1009 @@ isMergable(LeafType& leaf, const Coord& start, int dim,
 ////////////////////////////////////////
 
 
-template <class TreeT>
-class LeafCPtrList
+template<typename TreeT, typename LeafManagerT>
+class SignData
 {
 public:
-    typedef std::vector<const typename TreeT::LeafNodeType *> ListT;
+    typedef typename TreeT::ValueType ValueT;
+    typedef tree::ValueAccessor<const TreeT> AccessorT;
+
+    typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
+    typedef tree::ValueAccessor<IntTreeT> IntAccessorT;
+
+    typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
+    typedef tree::ValueAccessor<Int16TreeT> Int16AccessorT;
 
-    LeafCPtrList(const TreeT& tree)
+    //////////
+
+
+    SignData(const TreeT& distTree, const LeafManagerT& leafs, ValueT iso);
+
+    void run(bool threaded = true);
+
+    typename Int16TreeT::Ptr signTree() const { return mSignTree; }
+    typename IntTreeT::Ptr idxTree() const { return mIdxTree; }
+
+    //////////
+
+    SignData(SignData&, tbb::split);
+    void operator()(const tbb::blocked_range<size_t>&);
+    void join(const SignData& rhs)
     {
-        mLeafNodes.reserve(tree.leafCount());
-        typename TreeT::LeafCIter iter = tree.cbeginLeaf();
-        for ( ; iter; ++iter) mLeafNodes.push_back(iter.getLeaf());
+        mSignTree->merge(*rhs.mSignTree);
+        mIdxTree->merge(*rhs.mIdxTree);
     }
 
-    size_t size() const { return mLeafNodes.size(); }
+private:
 
-    const typename TreeT::LeafNodeType* operator[](size_t n) const
-        { return mLeafNodes[n]; }
+    const TreeT& mDistTree;
+    AccessorT mDistAcc;
 
-    tbb::blocked_range<size_t> getRange() const
-        { return tbb::blocked_range<size_t>(0, mLeafNodes.size()); }
+    const LeafManagerT& mLeafs;
+    ValueT mIsovalue;
 
-    const ListT& getList() const { return mLeafNodes; }
+    typename Int16TreeT::Ptr mSignTree;
+    Int16AccessorT mSignAcc;
+
+    typename IntTreeT::Ptr mIdxTree;
+    IntAccessorT mIdxAcc;
 
-private:
-    ListT mLeafNodes;
 };
 
 
-template <class TreeT>
-class LeafPtrList
+template<typename TreeT, typename LeafManagerT>
+SignData<TreeT, LeafManagerT>::SignData(const TreeT& distTree,
+    const LeafManagerT& leafs, ValueT iso)
+    : mDistTree(distTree)
+    , mDistAcc(mDistTree)
+    , mLeafs(leafs)
+    , mIsovalue(iso)
+    , mSignTree(new Int16TreeT(0))
+    , mSignAcc(*mSignTree)
+    , mIdxTree(new IntTreeT(int(util::INVALID_IDX)))
+    , mIdxAcc(*mIdxTree)
 {
-public:
-    typedef std::vector<typename TreeT::LeafNodeType *> ListT;
+}
 
-    LeafPtrList(TreeT& tree)
-    {
-        mLeafNodes.reserve(tree.leafCount());
-        typename TreeT::LeafIter iter = tree.beginLeaf();
-        for ( ; iter; ++iter) mLeafNodes.push_back(iter.getLeaf());
-    }
 
-    size_t size() const { return mLeafNodes.size(); }
+template<typename TreeT, typename LeafManagerT>
+SignData<TreeT, LeafManagerT>::SignData(SignData& rhs, tbb::split)
+    : mDistTree(rhs.mDistTree)
+    , mDistAcc(mDistTree)
+    , mLeafs(rhs.mLeafs)
+    , mIsovalue(rhs.mIsovalue)
+    , mSignTree(new Int16TreeT(0))
+    , mSignAcc(*mSignTree)
+    , mIdxTree(new IntTreeT(int(util::INVALID_IDX)))
+    , mIdxAcc(*mIdxTree)
+{
+}
+
+
+template<typename TreeT, typename LeafManagerT>
+void
+SignData<TreeT, LeafManagerT>::run(bool threaded)
+{
+    if (threaded) tbb::parallel_reduce(mLeafs.getRange(), *this);
+    else (*this)(mLeafs.getRange());
+}
+
+template<typename TreeT, typename LeafManagerT>
+void
+SignData<TreeT, LeafManagerT>::operator()(const tbb::blocked_range<size_t>& range)
+{
+    typedef typename Int16TreeT::LeafNodeType Int16LeafT;
+    typedef typename IntTreeT::LeafNodeType IntLeafT;
+    typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
+    unsigned char signs, face;
+    Coord ijk, coord;
+
+    typename internal::UniquePtr<Int16LeafT>::type signLeafPt(new Int16LeafT(ijk, 0));
 
-    typename TreeT::LeafNodeType* operator[](size_t n) const
-        { return mLeafNodes[n]; }
+    for (size_t n = range.begin(); n != range.end(); ++n) {
 
-    tbb::blocked_range<size_t> getRange() const
-        { return tbb::blocked_range<size_t>(0, mLeafNodes.size()); }
+        bool collectedData = false;
 
-    const ListT& getList() const { return mLeafNodes; }
+        coord = mLeafs.leaf(n).origin();
 
-private:
-    ListT mLeafNodes;
-};
+        if (!signLeafPt.get()) signLeafPt.reset(new Int16LeafT(coord, 0));
+        else signLeafPt->setOrigin(coord);
+
+        const typename TreeT::LeafNodeType *leafPt = mDistAcc.probeConstLeaf(coord);
+
+        coord.offset(TreeT::LeafNodeType::DIM - 1);
+
+        for (iter = mLeafs.leaf(n).cbeginValueOn(); iter; ++iter) {
+
+            ijk = iter.getCoord();
+
+            if (leafPt && ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2]) {
+                signs = evalCellSigns(*leafPt, iter.pos(), mIsovalue);
+            } else {
+                signs = evalCellSigns(mDistAcc, ijk, mIsovalue);
+            }
+
+            if (signs != 0 && signs != 0xFF) {
+                Int16 flags = (signs & 0x1) ? INSIDE : 0;
+
+                if (bool(signs & 0x1) != bool(signs & 0x2)) flags |= XEDGE;
+                if (bool(signs & 0x1) != bool(signs & 0x10)) flags |= YEDGE;
+                if (bool(signs & 0x1) != bool(signs & 0x8)) flags |= ZEDGE;
+
+                face = internal::sAmbiguousFace[signs];
+                if (face != 0) correctCellSigns(signs, face, mDistAcc, ijk, mIsovalue);
+
+                flags |= Int16(signs);
+
+                signLeafPt->setValue(ijk, flags);
+                collectedData = true;
+            }
+        }
+
+        if (collectedData) {
+
+            IntLeafT* idxLeaf = mIdxAcc.touchLeaf(coord);
+            idxLeaf->topologyUnion(*signLeafPt);
+            typename IntLeafT::ValueOnIter it = idxLeaf->beginValueOn();
+            for (; it; ++it) {
+                it.setValue(0);
+            }
+
+            mSignAcc.addLeaf(signLeafPt.release());
+        }
+    }
+}
 
 
 ////////////////////////////////////////
 
 
-template<typename DistTreeT>
-struct ReferenceData
+/// @brief Counts the total number of points per leaf, accounts for cells with multiple points.
+class CountPoints
 {
-    typedef typename DistTreeT::ValueType DistValueT;
-    typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
-    typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
-    typedef typename DistTreeT::template ValueConverter<Vec3s>::Type Vec3sTreeT;
-
-    ReferenceData()
-        : mDistTree(NULL)
-        , mEdgeTree(NULL)
-        , mTopologyMaskTree(NULL)
-        , mSeamPointTree(NULL)
-        , mSeamMaskTree(typename BoolTreeT::Ptr())
-        , mSmoothingMaskTree(typename BoolTreeT::Ptr())
-        , mInternalAdaptivity(DistValueT(0.0))
+public:
+    CountPoints(std::vector<size_t>& pointList) : mPointList(pointList) {}
+
+    template <typename LeafNodeType>
+    void operator()(LeafNodeType &leaf, size_t leafIndex) const
     {
+        size_t points = 0;
+
+        typename LeafNodeType::ValueOnCIter iter = leaf.cbeginValueOn();
+        for (; iter; ++iter) {
+            points += size_t(sEdgeGroupTable[(SIGNS & iter.getValue())][0]);
+        }
+
+        mPointList[leafIndex] = points;
     }
 
-    bool isValid() const
+private:
+    std::vector<size_t>& mPointList;
+};
+
+
+/// @brief Computes the point list indices for the index tree.
+template<typename Int16TreeT>
+class MapPoints
+{
+public:
+    typedef tree::ValueAccessor<const Int16TreeT> Int16AccessorT;
+
+    MapPoints(std::vector<size_t>& pointList, const Int16TreeT& signTree)
+        : mPointList(pointList)
+        , mSignAcc(signTree)
     {
-        return mDistTree && mEdgeTree && mTopologyMaskTree && mSeamMaskTree;
     }
 
-    const DistTreeT* mDistTree;
-    const CharTreeT* mEdgeTree;
-    BoolTreeT* mTopologyMaskTree;
-    Vec3sTreeT* mSeamPointTree;
-    typename BoolTreeT::Ptr mSeamMaskTree, mSmoothingMaskTree;
-    DistValueT mInternalAdaptivity;
-};
+    template <typename LeafNodeType>
+    void operator()(LeafNodeType &leaf, size_t leafIndex) const
+    {
+        size_t ptnIdx = mPointList[leafIndex];
+        typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
 
+        const typename Int16TreeT::LeafNodeType *signLeafPt =
+            mSignAcc.probeConstLeaf(leaf.origin());
 
-////////////////////////////////////////
+        for (; iter; ++iter) {
+            iter.setValue(ptnIdx);
+            unsigned signs = SIGNS & signLeafPt->getValue(iter.pos());
+            ptnIdx += size_t(sEdgeGroupTable[signs][0]);
+        }
+    }
+
+private:
+    std::vector<size_t>& mPointList;
+    Int16AccessorT mSignAcc;
+};
 
 
-template <class DistTreeT>
-class Count
+/// @brief Counts the total number of points per collapsed region
+template<typename IntTreeT>
+class CountRegions
 {
 public:
-    Count(const LeafPtrList<DistTreeT>&, std::vector<size_t>&);
-    inline Count(const Count<DistTreeT>&);
+    typedef tree::ValueAccessor<IntTreeT> IntAccessorT;
+    typedef typename IntTreeT::LeafNodeType IntLeafT;
+
+    CountRegions(IntTreeT& idxTree, std::vector<size_t>& regions)
+    : mIdxAcc(idxTree)
+    , mRegions(regions)
+    {
+    }
+
+    template <typename LeafNodeType>
+    void operator()(LeafNodeType &leaf, size_t leafIndex) const
+    {
 
-    void runParallel();
-    void runSerial();
+        size_t regions = 0;
+
+        IntLeafT tmpLeaf(*mIdxAcc.probeConstLeaf(leaf.origin()));
+
+        typename IntLeafT::ValueOnIter iter = tmpLeaf.beginValueOn();
+        for (; iter; ++iter) {
+            if(iter.getValue() == 0) {
+                iter.setValueOff();
+                regions += size_t(sEdgeGroupTable[(SIGNS & leaf.getValue(iter.pos()))][0]);
+            }
+        }
+
+        int onVoxelCount = int(tmpLeaf.onVoxelCount());
+        while (onVoxelCount > 0) {
+            ++regions;
+            iter = tmpLeaf.beginValueOn();
+            int regionId = iter.getValue();
+            for (; iter; ++iter) {
+                if (iter.getValue() == regionId) {
+                    iter.setValueOff();
+                    --onVoxelCount;
+                }
+            }
+        }
+
+        mRegions[leafIndex] = regions;
+    }
 
-    inline void operator()(const tbb::blocked_range<size_t>&) const;
 private:
-    const LeafPtrList<DistTreeT>& mLeafNodes;
-    std::vector<size_t>& mLeafRegionCount;
+    IntAccessorT mIdxAcc;
+    std::vector<size_t>& mRegions;
 };
 
 
-template <class DistTreeT>
-Count<DistTreeT>::Count(
-    const LeafPtrList<DistTreeT>& leafs,
-    std::vector<size_t>& leafRegionCount)
-    : mLeafNodes(leafs)
-    , mLeafRegionCount(leafRegionCount)
+////////////////////////////////////////
+
+
+// @brief linear interpolation.
+inline double evalRoot(double v0, double v1, double iso) { return (iso - v0) / (v1 - v0); }
+
+
+/// @brief Extracts the eight corner values for leaf inclusive cells.
+template<typename LeafT>
+inline void
+collectCornerValues(const LeafT& leaf, const Index offset, std::vector<double>& values)
 {
+    values[0] = double(leaf.getValue(offset)); // i, j, k
+    values[3] = double(leaf.getValue(offset + 1)); // i, j, k+1
+    values[4] = double(leaf.getValue(offset + LeafT::DIM)); // i, j+1, k
+    values[7] = double(leaf.getValue(offset + LeafT::DIM + 1)); // i, j+1, k+1
+    values[1] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM))); // i+1, j, k
+    values[2] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + 1)); // i+1, j, k+1
+    values[5] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM)); // i+1, j+1, k
+    values[6] = double(leaf.getValue(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1)); // i+1, j+1, k+1
 }
 
 
-template <class DistTreeT>
-inline
-Count<DistTreeT>::Count(const Count<DistTreeT>& rhs)
-    : mLeafNodes(rhs.mLeafNodes)
-    , mLeafRegionCount(rhs.mLeafRegionCount)
+/// @brief Extracts the eight corner values for a cell starting at the given @ijk coordinate.
+template<typename AccessorT>
+inline void
+collectCornerValues(const AccessorT& acc, const Coord& ijk, std::vector<double>& values)
 {
+    Coord coord = ijk;
+    values[0] = double(acc.getValue(coord)); // i, j, k
+
+    coord[0] += 1;
+    values[1] = double(acc.getValue(coord)); // i+1, j, k
+
+    coord[2] += 1;
+    values[2] = double(acc.getValue(coord)); // i+i, j, k+1
+
+    coord[0] = ijk[0];
+    values[3] = double(acc.getValue(coord)); // i, j, k+1
+
+    coord[1] += 1; coord[2] = ijk[2];
+    values[4] = double(acc.getValue(coord)); // i, j+1, k
+
+    coord[0] += 1;
+    values[5] = double(acc.getValue(coord)); // i+1, j+1, k
+
+    coord[2] += 1;
+    values[6] = double(acc.getValue(coord)); // i+1, j+1, k+1
+
+    coord[0] = ijk[0];
+    values[7] = double(acc.getValue(coord)); // i, j+1, k+1
 }
 
 
-template <class DistTreeT>
-void
-Count<DistTreeT>::runParallel()
+/// @brief Computes the average cell point for a given edge group.
+inline Vec3d
+computePoint(const std::vector<double>& values, unsigned char signs,
+    unsigned char edgeGroup, double iso)
 {
-    tbb::parallel_for(mLeafNodes.getRange(), *this);
+    Vec3d avg(0.0, 0.0, 0.0);
+    int samples = 0;
+
+    if (sEdgeGroupTable[signs][1] == edgeGroup) { // Edged: 0 - 1
+        avg[0] += evalRoot(values[0], values[1], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][2] == edgeGroup) { // Edged: 1 - 2
+        avg[0] += 1.0;
+        avg[2] += evalRoot(values[1], values[2], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][3] == edgeGroup) { // Edged: 3 - 2
+        avg[0] += evalRoot(values[3], values[2], iso);
+        avg[2] += 1.0;
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][4] == edgeGroup) { // Edged: 0 - 3
+        avg[2] += evalRoot(values[0], values[3], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][5] == edgeGroup) { // Edged: 4 - 5
+        avg[0] += evalRoot(values[4], values[5], iso);
+        avg[1] += 1.0;
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][6] == edgeGroup) { // Edged: 5 - 6
+        avg[0] += 1.0;
+        avg[1] += 1.0;
+        avg[2] += evalRoot(values[5], values[6], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][7] == edgeGroup) { // Edged: 7 - 6
+        avg[0] += evalRoot(values[7], values[6], iso);
+        avg[1] += 1.0;
+        avg[2] += 1.0;
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][8] == edgeGroup) { // Edged: 4 - 7
+        avg[1] += 1.0;
+        avg[2] += evalRoot(values[4], values[7], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][9] == edgeGroup) { // Edged: 0 - 4
+        avg[1] += evalRoot(values[0], values[4], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][10] == edgeGroup) { // Edged: 1 - 5
+        avg[0] += 1.0;
+        avg[1] += evalRoot(values[1], values[5], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][11] == edgeGroup) { // Edged: 2 - 6
+        avg[0] += 1.0;
+        avg[1] += evalRoot(values[2], values[6], iso);
+        avg[2] += 1.0;
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][12] == edgeGroup) { // Edged: 3 - 7
+        avg[1] += evalRoot(values[3], values[7], iso);
+        avg[2] += 1.0;
+        ++samples;
+    }
+
+    if (samples > 1) {
+        double w = 1.0 / double(samples);
+        avg[0] *= w;
+        avg[1] *= w;
+        avg[2] *= w;
+    }
+
+    return avg;
 }
 
 
-template <class DistTreeT>
-void
-Count<DistTreeT>::runSerial()
+/// @brief  Computes the average cell point for a given edge group, ignoring edge
+///         samples present in the @c signsMask configuration.
+inline int
+computeMaskedPoint(Vec3d& avg, const std::vector<double>& values, unsigned char signs,
+    unsigned char signsMask, unsigned char edgeGroup, double iso)
+{
+    avg = Vec3d(0.0, 0.0, 0.0);
+    int samples = 0;
+
+    if (sEdgeGroupTable[signs][1] == edgeGroup
+        && sEdgeGroupTable[signsMask][1] == 0) { // Edged: 0 - 1
+        avg[0] += evalRoot(values[0], values[1], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][2] == edgeGroup
+        && sEdgeGroupTable[signsMask][2] == 0) { // Edged: 1 - 2
+        avg[0] += 1.0;
+        avg[2] += evalRoot(values[1], values[2], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][3] == edgeGroup
+        && sEdgeGroupTable[signsMask][3] == 0) { // Edged: 3 - 2
+        avg[0] += evalRoot(values[3], values[2], iso);
+        avg[2] += 1.0;
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][4] == edgeGroup
+        && sEdgeGroupTable[signsMask][4] == 0) { // Edged: 0 - 3
+        avg[2] += evalRoot(values[0], values[3], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][5] == edgeGroup
+        && sEdgeGroupTable[signsMask][5] == 0) { // Edged: 4 - 5
+        avg[0] += evalRoot(values[4], values[5], iso);
+        avg[1] += 1.0;
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][6] == edgeGroup
+        && sEdgeGroupTable[signsMask][6] == 0) { // Edged: 5 - 6
+        avg[0] += 1.0;
+        avg[1] += 1.0;
+        avg[2] += evalRoot(values[5], values[6], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][7] == edgeGroup
+        && sEdgeGroupTable[signsMask][7] == 0) { // Edged: 7 - 6
+        avg[0] += evalRoot(values[7], values[6], iso);
+        avg[1] += 1.0;
+        avg[2] += 1.0;
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][8] == edgeGroup
+        && sEdgeGroupTable[signsMask][8] == 0) { // Edged: 4 - 7
+        avg[1] += 1.0;
+        avg[2] += evalRoot(values[4], values[7], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][9] == edgeGroup
+        && sEdgeGroupTable[signsMask][9] == 0) { // Edged: 0 - 4
+        avg[1] += evalRoot(values[0], values[4], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][10] == edgeGroup
+        && sEdgeGroupTable[signsMask][10] == 0) { // Edged: 1 - 5
+        avg[0] += 1.0;
+        avg[1] += evalRoot(values[1], values[5], iso);
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][11] == edgeGroup
+        && sEdgeGroupTable[signsMask][11] == 0) { // Edged: 2 - 6
+        avg[0] += 1.0;
+        avg[1] += evalRoot(values[2], values[6], iso);
+        avg[2] += 1.0;
+        ++samples;
+    }
+
+    if (sEdgeGroupTable[signs][12] == edgeGroup
+        && sEdgeGroupTable[signsMask][12] == 0) { // Edged: 3 - 7
+        avg[1] += evalRoot(values[3], values[7], iso);
+        avg[2] += 1.0;
+        ++samples;
+    }
+
+    if (samples > 1) {
+        double w = 1.0 / double(samples);
+        avg[0] *= w;
+        avg[1] *= w;
+        avg[2] *= w;
+    }
+
+    return samples;
+}
+
+
+/// @brief  Computes the average cell point for a given edge group, by computing
+///         convex weights based on the distance from the sample point @c p.
+inline Vec3d
+computeWeightedPoint(const Vec3d& p, const std::vector<double>& values,
+    unsigned char signs, unsigned char edgeGroup, double iso)
 {
-    (*this)(mLeafNodes.getRange());
+    std::vector<Vec3d> samples;
+    samples.reserve(8);
+
+    std::vector<double> weights;
+    weights.reserve(8);
+
+    Vec3d avg(0.0, 0.0, 0.0);
+
+    if (sEdgeGroupTable[signs][1] == edgeGroup) { // Edged: 0 - 1
+        avg[0] = evalRoot(values[0], values[1], iso);
+        avg[1] = 0.0;
+        avg[2] = 0.0;
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][2] == edgeGroup) { // Edged: 1 - 2
+        avg[0] = 1.0;
+        avg[1] = 0.0;
+        avg[2] = evalRoot(values[1], values[2], iso);
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][3] == edgeGroup) { // Edged: 3 - 2
+        avg[0] = evalRoot(values[3], values[2], iso);
+        avg[1] = 0.0;
+        avg[2] = 1.0;
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][4] == edgeGroup) { // Edged: 0 - 3
+        avg[0] = 0.0;
+        avg[1] = 0.0;
+        avg[2] = evalRoot(values[0], values[3], iso);
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][5] == edgeGroup) { // Edged: 4 - 5
+        avg[0] = evalRoot(values[4], values[5], iso);
+        avg[1] = 1.0;
+        avg[2] = 0.0;
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][6] == edgeGroup) { // Edged: 5 - 6
+        avg[0] = 1.0;
+        avg[1] = 1.0;
+        avg[2] = evalRoot(values[5], values[6], iso);
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][7] == edgeGroup) { // Edged: 7 - 6
+        avg[0] = evalRoot(values[7], values[6], iso);
+        avg[1] = 1.0;
+        avg[2] = 1.0;
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][8] == edgeGroup) { // Edged: 4 - 7
+        avg[0] = 0.0;
+        avg[1] = 1.0;
+        avg[2] = evalRoot(values[4], values[7], iso);
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][9] == edgeGroup) { // Edged: 0 - 4
+        avg[0] = 0.0;
+        avg[1] = evalRoot(values[0], values[4], iso);
+        avg[2] = 0.0;
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][10] == edgeGroup) { // Edged: 1 - 5
+        avg[0] = 1.0;
+        avg[1] = evalRoot(values[1], values[5], iso);
+        avg[2] = 0.0;
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][11] == edgeGroup) { // Edged: 2 - 6
+        avg[0] = 1.0;
+        avg[1] = evalRoot(values[2], values[6], iso);
+        avg[2] = 1.0;
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+    if (sEdgeGroupTable[signs][12] == edgeGroup) { // Edged: 3 - 7
+        avg[0] = 0.0;
+        avg[1] = evalRoot(values[3], values[7], iso);
+        avg[2] = 1.0;
+
+        samples.push_back(avg);
+        weights.push_back((avg-p).lengthSqr());
+    }
+
+
+    double minWeight = std::numeric_limits<double>::max();
+    double maxWeight = -std::numeric_limits<double>::max();
+
+    for (size_t i = 0, I = weights.size(); i < I; ++i) {
+        minWeight = std::min(minWeight, weights[i]);
+        maxWeight = std::max(maxWeight, weights[i]);
+    }
+
+    const double offset = maxWeight + minWeight * 0.1;
+    for (size_t i = 0, I = weights.size(); i < I; ++i) {
+        weights[i] = offset - weights[i];
+    }
+
+
+    double weightSum = 0.0;
+    for (size_t i = 0, I = weights.size(); i < I; ++i) {
+        weightSum += weights[i];
+    }
+
+    avg[0] = 0.0;
+    avg[1] = 0.0;
+    avg[2] = 0.0;
+
+    if (samples.size() > 1) {
+        for (size_t i = 0, I = samples.size(); i < I; ++i) {
+            avg += samples[i] * (weights[i] / weightSum);
+        }
+    } else {
+        avg = samples.front();
+    }
+
+    return avg;
 }
 
 
-template <class DistTreeT>
+/// @brief  Computes the average cell points defined by the sign configuration
+///         @c signs and the given corner values @c values.
 inline void
-Count<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
+computeCellPoints(std::vector<Vec3d>& points,
+    const std::vector<double>& values, unsigned char signs, double iso)
 {
-    for (size_t n = range.begin(); n != range.end(); ++n) {
-        mLeafRegionCount[n] = size_t(mLeafNodes[n]->onVoxelCount());
+    for (size_t n = 1, N = sEdgeGroupTable[signs][0] + 1; n < N; ++n) {
+        points.push_back(computePoint(values, signs, n, iso));
     }
 }
 
 
-////////////////////////////////////////
+/// @brief  Given a sign configuration @c lhsSigns and an edge group @c groupId,
+///         finds the corresponding edge group in a different sign configuration
+///         @c rhsSigns. Returns -1 if no match is found.
+inline int
+matchEdgeGroup(unsigned char groupId, unsigned char lhsSigns, unsigned char rhsSigns)
+{
+    int id = -1;
+    for (size_t i = 1; i <= 12; ++i) {
+        if (sEdgeGroupTable[lhsSigns][i] == groupId && sEdgeGroupTable[rhsSigns][i] != 0) {
+            id = sEdgeGroupTable[rhsSigns][i];
+            break;
+        }
+    }
+    return id;
+}
 
 
-template <class DistTreeT>
-class Merge
+/// @brief  Computes the average cell points defined by the sign configuration
+///         @c signs and the given corner values @c values. Combines data from
+///         two different level sets to eliminate seam lines when meshing
+///         fractured segments.
+inline void
+computeCellPoints(std::vector<Vec3d>& points, std::vector<bool>& weightedPointMask,
+    const std::vector<double>& lhsValues, const std::vector<double>& rhsValues,
+    unsigned char lhsSigns, unsigned char rhsSigns,
+    double iso, size_t pointIdx, const boost::scoped_array<uint32_t>& seamPoints)
+{
+    for (size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n) {
+
+        int id = matchEdgeGroup(n, lhsSigns, rhsSigns);
+
+        if (id != -1) {
+
+            const unsigned char e(id);
+            uint32_t& quantizedPoint = seamPoints[pointIdx + (id - 1)];
+
+            if ((quantizedPoint & MASK_DIRTY_BIT) && !(quantizedPoint & MASK_INVALID_BIT)) {
+                Vec3d p = unpackPoint(quantizedPoint);
+                points.push_back(computeWeightedPoint(p, rhsValues, rhsSigns, e, iso));
+                weightedPointMask.push_back(true);
+            } else {
+                points.push_back(computePoint(rhsValues, rhsSigns, e, iso));
+                weightedPointMask.push_back(false);
+            }
+
+        } else {
+            points.push_back(computePoint(lhsValues, lhsSigns, n, iso));
+            weightedPointMask.push_back(false);
+        }
+    }
+}
+
+
+template <typename TreeT, typename LeafManagerT>
+class GenPoints
 {
 public:
-    typedef typename DistTreeT::ValueType DistValueT;
-    typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
-    typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
+    typedef tree::ValueAccessor<const TreeT> AccessorT;
+
+    typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
+    typedef tree::ValueAccessor<IntTreeT> IntAccessorT;
+    typedef tree::ValueAccessor<const IntTreeT> IntCAccessorT;
 
-    Merge(
-        const DistTreeT& distTree,
-        LeafPtrList<IntTreeT>& auxLeafs,
-        std::vector<size_t>& leafRegionCount,
-        const DistValueT iso,
-        const DistValueT adaptivity);
+    typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
+    typedef tree::ValueAccessor<const Int16TreeT> Int16CAccessorT;
 
-    inline Merge(const Merge<DistTreeT>&);
+    typedef boost::scoped_array<uint32_t> QuantizedPointList;
 
-    void setRefData(const ReferenceData<DistTreeT>&);
+    //////////
+
+
+    GenPoints(const LeafManagerT& signLeafs, const TreeT& distTree,
+        IntTreeT& idxTree, PointList& points, std::vector<size_t>& indices,
+        const math::Transform& xform, double iso);
+
+    void run(bool threaded = true);
+
+    void setRefData(const Int16TreeT* refSignTree = NULL, const TreeT* refDistTree = NULL,
+        IntTreeT* refIdxTree = NULL, const QuantizedPointList* seamPoints = NULL,
+        std::vector<unsigned char>* mSeamPointMaskPt = NULL);
+
+    //////////
 
-    void runParallel();
-    void runSerial();
 
     void operator()(const tbb::blocked_range<size_t>&) const;
 
 private:
-    const DistTreeT& mDistTree;
-    LeafPtrList<IntTreeT>& mAuxLeafs;
-    std::vector<size_t>& mLeafRegionCount;
-    const DistValueT mIsovalue, mAdaptivity;
-    const ReferenceData<DistTreeT>* mRefData;
-};
+    const LeafManagerT& mSignLeafs;
 
+    AccessorT mDistAcc;
+    IntTreeT& mIdxTree;
 
-template <class DistTreeT>
-Merge<DistTreeT>::Merge(
-    const DistTreeT& distTree,
-    LeafPtrList<IntTreeT>& auxLeafs,
-    std::vector<size_t>& leafRegionCount,
-    const DistValueT iso,
-    const DistValueT adaptivity)
-    : mDistTree(distTree)
-    , mAuxLeafs(auxLeafs)
-    , mLeafRegionCount(leafRegionCount)
-    , mIsovalue(iso)
-    , mAdaptivity(adaptivity)
-    , mRefData(NULL)
-{
-}
+    PointList& mPoints;
+    std::vector<size_t>& mIndices;
+    const math::Transform& mTransform;
+    const double mIsovalue;
 
+    // reference data
+    const Int16TreeT *mRefSignTreePt;
+    const TreeT* mRefDistTreePt;
+    const IntTreeT* mRefIdxTreePt;
+    const QuantizedPointList* mSeamPointsPt;
+    std::vector<unsigned char>* mSeamPointMaskPt;
+};
 
-template <class DistTreeT>
-inline
-Merge<DistTreeT>::Merge(const Merge<DistTreeT>& rhs)
-    : mDistTree(rhs.mDistTree)
-    , mAuxLeafs(rhs.mAuxLeafs)
-    , mLeafRegionCount(rhs.mLeafRegionCount)
-    , mIsovalue(rhs.mIsovalue)
-    , mAdaptivity(rhs.mAdaptivity)
-    , mRefData(rhs.mRefData)
+
+template <typename TreeT, typename LeafManagerT>
+GenPoints<TreeT, LeafManagerT>::GenPoints(const LeafManagerT& signLeafs,
+    const TreeT& distTree, IntTreeT& idxTree, PointList& points,
+    std::vector<size_t>& indices, const math::Transform& xform, double iso)
+    : mSignLeafs(signLeafs)
+    , mDistAcc(distTree)
+    , mIdxTree(idxTree)
+    , mPoints(points)
+    , mIndices(indices)
+    , mTransform(xform)
+    , mIsovalue(iso)
+    , mRefSignTreePt(NULL)
+    , mRefDistTreePt(NULL)
+    , mRefIdxTreePt(NULL)
+    , mSeamPointsPt(NULL)
+    , mSeamPointMaskPt(NULL)
 {
 }
 
 
-template <class DistTreeT>
+template <typename TreeT, typename LeafManagerT>
 void
-Merge<DistTreeT>::runParallel()
+GenPoints<TreeT, LeafManagerT>::run(bool threaded)
 {
-    tbb::parallel_for(mAuxLeafs.getRange(), *this);
+    if (threaded) tbb::parallel_for(mSignLeafs.getRange(), *this);
+    else (*this)(mSignLeafs.getRange());
 }
 
 
-template <class DistTreeT>
+template <typename TreeT, typename LeafManagerT>
 void
-Merge<DistTreeT>::runSerial()
+GenPoints<TreeT, LeafManagerT>::setRefData(const Int16TreeT *refSignTree, const TreeT *refDistTree,
+    IntTreeT* refIdxTree, const QuantizedPointList* seamPoints, std::vector<unsigned char>* seamPointMask)
 {
-    (*this)(mAuxLeafs.getRange());
+    mRefSignTreePt = refSignTree;
+    mRefDistTreePt = refDistTree;
+    mRefIdxTreePt = refIdxTree;
+    mSeamPointsPt = seamPoints;
+    mSeamPointMaskPt = seamPointMask;
 }
 
-template <class DistTreeT>
-void
-Merge<DistTreeT>::setRefData(const ReferenceData<DistTreeT>& refData)
-{
-    mRefData = &refData;
-}
 
-template <class DistTreeT>
+template <typename TreeT, typename LeafManagerT>
 void
-Merge<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
+GenPoints<TreeT, LeafManagerT>::operator()(
+    const tbb::blocked_range<size_t>& range) const
 {
-    typedef math::Vec3<DistValueT> Vec3T;
-    typedef typename BoolTreeT::LeafNodeType BoolLeafT;
-    typedef typename IntTreeT::LeafNodeType IntLeafT;
-    typedef typename BoolLeafT::template ValueConverter<Vec3T>::Type Vec3LeafT;
+    typename IntTreeT::LeafNodeType::ValueOnIter iter;
+    unsigned char signs, refSigns;
+    Index offset;
+    Coord ijk, coord;
+    std::vector<Vec3d> points(4);
+    std::vector<bool> weightedPointMask(4);
+    std::vector<double> values(8), refValues(8);
 
-    typedef typename IntLeafT::ValueOnIter IntIterT;
-    typedef typename BoolLeafT::ValueOnCIter BoolCIterT;
 
-    typedef typename tree::ValueAccessor<BoolTreeT> BoolTreeAccessorT;
-    typedef typename tree::ValueAccessor<const BoolTreeT> BoolTreeCAccessorT;
+    IntAccessorT idxAcc(mIdxTree);
 
-    boost::scoped_ptr<BoolTreeAccessorT> seamMaskAcc;
-    boost::scoped_ptr<BoolTreeCAccessorT> topologyMaskAcc;
-    if (mRefData && mRefData->isValid()) {
-        seamMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mSeamMaskTree.get()));
-        topologyMaskAcc.reset(new BoolTreeCAccessorT(*mRefData->mTopologyMaskTree));
-    }
-    const bool hasRefData = seamMaskAcc && topologyMaskAcc;
+    // reference data accessors
+    boost::scoped_ptr<Int16CAccessorT> refSignAcc;
+    if (mRefSignTreePt) refSignAcc.reset(new Int16CAccessorT(*mRefSignTreePt));
 
-    const int LeafDim = BoolLeafT::DIM;
-    tree::ValueAccessor<const DistTreeT> distAcc(mDistTree);
+    boost::scoped_ptr<IntCAccessorT> refIdxAcc;
+    if (mRefIdxTreePt) refIdxAcc.reset(new IntCAccessorT(*mRefIdxTreePt));
+
+    boost::scoped_ptr<AccessorT> refDistAcc;
+    if (mRefDistTreePt) refDistAcc.reset(new AccessorT(*mRefDistTreePt));
 
-    // Allocate reusable leaf buffers
-    BoolLeafT mask;
-    Vec3LeafT gradientBuffer;
-    Coord ijk, coord, end;
 
     for (size_t n = range.begin(); n != range.end(); ++n) {
 
-        DistValueT adaptivity = mAdaptivity;
-        IntLeafT& auxLeaf = *mAuxLeafs[n];
+        coord = mSignLeafs.leaf(n).origin();
 
-        const Coord& origin = auxLeaf.getOrigin();
-        end[0] = origin[0] + LeafDim;
-        end[1] = origin[1] + LeafDim;
-        end[2] = origin[2] + LeafDim;
+        const typename TreeT::LeafNodeType *leafPt = mDistAcc.probeConstLeaf(coord);
+        typename IntTreeT::LeafNodeType *idxLeafPt = idxAcc.probeLeaf(coord);
 
-        mask.setValuesOff();
 
-        // Mask off seam line adjacent voxels
-        if (hasRefData) {
-            const BoolLeafT* seamMask = seamMaskAcc->probeConstLeaf(origin);
-            if (seamMask != NULL) {
-                for (BoolCIterT it = seamMask->cbeginValueOn(); it; ++it) {
-                    ijk = it.getCoord();
-                    coord[0] = ijk[0] - (ijk[0] % 2);
-                    coord[1] = ijk[1] - (ijk[1] % 2);
-                    coord[2] = ijk[2] - (ijk[2] % 2);
-                    mask.setActiveState(coord, true);
+        // reference data leafs
+        const typename Int16TreeT::LeafNodeType *refSignLeafPt = NULL;
+        if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(coord);
+
+        const typename IntTreeT::LeafNodeType *refIdxLeafPt = NULL;
+        if (refIdxAcc) refIdxLeafPt = refIdxAcc->probeConstLeaf(coord);
+
+        const typename TreeT::LeafNodeType *refDistLeafPt = NULL;
+        if (refDistAcc) refDistLeafPt = refDistAcc->probeConstLeaf(coord);
+
+
+        // generate cell points
+        size_t ptnIdx = mIndices[n];
+        coord.offset(TreeT::LeafNodeType::DIM - 1);
+
+
+
+        for (iter = idxLeafPt->beginValueOn(); iter; ++iter) {
+
+            if(iter.getValue() != 0) continue;
+
+            iter.setValue(ptnIdx);
+            iter.setValueOff();
+            offset = iter.pos();
+            ijk = iter.getCoord();
+
+            const bool inclusiveCell = ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2];
+
+            const Int16& flags = mSignLeafs.leaf(n).getValue(offset);
+            signs    = SIGNS & flags;
+            refSigns = 0;
+
+            if ((flags & SEAM) && refSignLeafPt && refIdxLeafPt) {
+                if (refSignLeafPt->isValueOn(offset)) {
+                    refSigns = (SIGNS & refSignLeafPt->getValue(offset));
                 }
             }
-            if (topologyMaskAcc->probeConstLeaf(origin) == NULL) {
-                adaptivity = mRefData->mInternalAdaptivity;
+
+
+            if (inclusiveCell) collectCornerValues(*leafPt, offset, values);
+            else collectCornerValues(mDistAcc, ijk, values);
+
+
+            points.clear();
+            weightedPointMask.clear();
+
+            if (refSigns == 0) {
+                computeCellPoints(points, values, signs, mIsovalue);
+            } else {
+
+                if (inclusiveCell) collectCornerValues(*refDistLeafPt, offset, refValues);
+                else collectCornerValues(*refDistAcc, ijk, refValues);
+
+                computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
+                    mIsovalue, refIdxLeafPt->getValue(offset), *mSeamPointsPt);
+            }
+
+
+            for (size_t i = 0, I = points.size(); i < I; ++i) {
+
+                // offset by cell-origin
+                points[i][0] += double(ijk[0]);
+                points[i][1] += double(ijk[1]);
+                points[i][2] += double(ijk[2]);
+
+
+                points[i] = mTransform.indexToWorld(points[i]);
+
+                mPoints[ptnIdx][0] = float(points[i][0]);
+                mPoints[ptnIdx][1] = float(points[i][1]);
+                mPoints[ptnIdx][2] = float(points[i][2]);
+
+                if (mSeamPointMaskPt && !weightedPointMask.empty() && weightedPointMask[i]) {
+                    (*mSeamPointMaskPt)[ptnIdx] = 1;
+                }
+
+                ++ptnIdx;
             }
         }
 
-        // Mask off ambiguous voxels
-        for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
-            ijk = it.getCoord();
-            coord[0] = ijk[0] - (ijk[0] % 2);
-            coord[1] = ijk[1] - (ijk[1] % 2);
-            coord[2] = ijk[2] - (ijk[2] % 2);
-            if(mask.isValueOn(coord)) continue;
-            mask.setActiveState(coord, isAmbiguous(distAcc, ijk, mIsovalue, 1));
-        }
-
-       int dim = 2;
-       // Mask off topologically ambiguous 2x2x2 voxel sub-blocks
-       for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
-           for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
-               for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
-                   if (isNonManifold(distAcc, ijk, mIsovalue, dim)) {
-                       mask.setActiveState(ijk, true);
-                   }
-               }
-           }
-       }
-
-       // Compute the gradient for the remaining voxels
-       gradientBuffer.setValuesOff();
-
-       for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
-
-           ijk = it.getCoord();
-           coord[0] = ijk[0] - (ijk[0] % dim);
-           coord[1] = ijk[1] - (ijk[1] % dim);
-           coord[2] = ijk[2] - (ijk[2] % dim);
-           if(mask.isValueOn(coord)) continue;
-
-           Vec3T norm(math::ISGradient<math::CD_2ND>::result(distAcc, ijk));
-           // Normalize (Vec3's normalize uses isApproxEqual, which uses abs and does more work)
-           DistValueT length = norm.length();
-           if (length > DistValueT(1.0e-7)) {
-               norm *= DistValueT(1.0) / length;
-           }
-           gradientBuffer.setValue(ijk, norm);
-       }
-
-       int regionId = 1, next_dim = dim << 1;
-
-       // Process the first adaptivity level.
-       for (ijk[0] = 0; ijk[0] < LeafDim; ijk[0] += dim) {
-           coord[0] = ijk[0] - (ijk[0] % next_dim);
-           for (ijk[1] = 0; ijk[1] < LeafDim; ijk[1] += dim) {
-               coord[1] = ijk[1] - (ijk[1] % next_dim);
-               for (ijk[2] = 0; ijk[2] < LeafDim; ijk[2] += dim) {
-                   coord[2] = ijk[2] - (ijk[2] % next_dim);
-                   if(mask.isValueOn(ijk) || !isMergable(gradientBuffer, ijk, dim, adaptivity)) {
-                       mask.setActiveState(coord, true);
-                       continue;
-                   }
-                   mergeVoxels(auxLeaf, ijk, dim, regionId++);
-               }
-           }
-       }
-
-
-       // Process remaining adaptivity levels
-       for (dim = 4; dim < LeafDim; dim = dim << 1) {
-           next_dim = dim << 1;
-           coord[0] = ijk[0] - (ijk[0] % next_dim);
-           for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
-               coord[1] = ijk[1] - (ijk[1] % next_dim);
-               for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
-                   coord[2] = ijk[2] - (ijk[2] % next_dim);
-                   for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
-
-                       if (mask.isValueOn(ijk) || isNonManifold(distAcc, ijk, mIsovalue, dim) ||
-                           !isMergable(gradientBuffer, ijk, dim, adaptivity)) {
-                           mask.setActiveState(coord, true);
-                           continue;
-                       }
-                       mergeVoxels(auxLeaf, ijk, dim, regionId++);
-                   }
-               }
-           }
-       }
-
-
-       if (!(mask.isValueOn(origin) || isNonManifold(distAcc, origin, mIsovalue, LeafDim))
-           && isMergable(gradientBuffer, origin, LeafDim, adaptivity)) {
-           mergeVoxels(auxLeaf, origin, LeafDim, regionId++);
-       }
-
-
-       // Count unique regions
-       size_t numVoxels = 0;
-       IntLeafT tmpLeaf(auxLeaf);
-       for (IntIterT it = tmpLeaf.beginValueOn(); it; ++it) {
-           if(it.getValue() == 0) {
-               it.setValueOff();
-               ++numVoxels;
-           }
-       }
-
-       while (tmpLeaf.onVoxelCount() > 0) {
-           ++numVoxels;
-           IntIterT it = tmpLeaf.beginValueOn();
-           regionId = it.getValue();
-           for (; it; ++it) {
-               if (it.getValue() == regionId) it.setValueOff();
-           }
-       }
-
-       mLeafRegionCount[n] = numVoxels;
+        // generate collapsed region points
+        int onVoxelCount = int(idxLeafPt->onVoxelCount());
+        while (onVoxelCount > 0) {
+
+            iter = idxLeafPt->beginValueOn();
+            int regionId = iter.getValue(), count = 0;
+
+            Vec3d avg(0.0), point;
+
+            for (; iter; ++iter) {
+                if (iter.getValue() != regionId) continue;
+
+                iter.setValue(ptnIdx);
+                iter.setValueOff();
+                --onVoxelCount;
+
+                ijk = iter.getCoord();
+                offset = iter.pos();
+
+                signs = (SIGNS & mSignLeafs.leaf(n).getValue(offset));
+
+                if (ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2]) {
+                    collectCornerValues(*leafPt, offset, values);
+                } else {
+                    collectCornerValues(mDistAcc, ijk, values);
+                }
+
+                points.clear();
+                computeCellPoints(points, values, signs, mIsovalue);
+
+                avg[0] += double(ijk[0]) + points[0][0];
+                avg[1] += double(ijk[1]) + points[0][1];
+                avg[2] += double(ijk[2]) + points[0][2];
+
+                ++count;
+            }
+
+
+            if (count > 1) {
+                double w = 1.0 / double(count);
+                avg[0] *= w;
+                avg[1] *= w;
+                avg[2] *= w;
+            }
+
+            avg = mTransform.indexToWorld(avg);
+
+            mPoints[ptnIdx][0] = float(avg[0]);
+            mPoints[ptnIdx][1] = float(avg[1]);
+            mPoints[ptnIdx][2] = float(avg[2]);
+
+            ++ptnIdx;
+        }
     }
 }
 
@@ -940,374 +1900,436 @@ Merge<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
 ////////////////////////////////////////
 
 
-template <class DistTreeT>
-class PointGen
+template<typename TreeT>
+class SeamWeights
 {
 public:
-    typedef typename DistTreeT::ValueType DistValueT;
-    typedef tree::ValueAccessor<const DistTreeT> DistTreeAccessorT;
-    typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
+    typedef tree::ValueAccessor<const TreeT> AccessorT;
 
-    PointGen(
-        const DistTreeT& distTree,
-        const LeafPtrList<IntTreeT>& auxLeafs,
-        std::vector<size_t>& leafIndices,
-        const openvdb::math::Transform& xform,
-        PointList& points,
-        double iso = 0.0);
+    typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
+    typedef tree::ValueAccessor<const IntTreeT> IntAccessorT;
 
-    PointGen(const PointGen<DistTreeT>&);
+    typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
+    typedef tree::ValueAccessor<const Int16TreeT> Int16AccessorT;
 
-    void setRefData(const ReferenceData<DistTreeT>&);
+    typedef boost::scoped_array<uint32_t> QuantizedPointList;
 
-    void runParallel();
-    void runSerial();
+    //////////
 
-    void operator()(const tbb::blocked_range<size_t>&) const;
+    SeamWeights(const TreeT& distTree, const Int16TreeT& refSignTree,
+        IntTreeT& refIdxTree, QuantizedPointList& points, double iso);
+
+    template <typename LeafNodeType>
+    void operator()(LeafNodeType &signLeaf, size_t leafIndex) const;
 
 private:
-    const DistTreeT& mDistTree;
-    const LeafPtrList<IntTreeT>& mAuxLeafs;
-    const std::vector<size_t>& mLeafIndices;
-    const openvdb::math::Transform& mTransform;
-    const PointList& mPoints;
-    const double mIsovalue;
-    const ReferenceData<DistTreeT>* mRefData;
+    AccessorT mDistAcc;
+    Int16AccessorT mRefSignAcc;
+    IntAccessorT mRefIdxAcc;
 
-    double root(double v0, double v1) const { return (mIsovalue - v0) / (v1 - v0); }
-    int calcAvgPoint(DistTreeAccessorT&, const Coord&, openvdb::Vec3d&) const;
+    QuantizedPointList& mPoints;
+    const double mIsovalue;
 };
 
 
-template <class DistTreeT>
-PointGen<DistTreeT>::PointGen(
-    const DistTreeT& distTree,
-    const LeafPtrList<IntTreeT>& auxLeafs,
-    std::vector<size_t>& leafIndices,
-    const openvdb::math::Transform& xform,
-    PointList& points,
-    double iso)
-    : mDistTree(distTree)
-    , mAuxLeafs(auxLeafs)
-    , mLeafIndices(leafIndices)
-    , mTransform(xform)
+template<typename TreeT>
+SeamWeights<TreeT>::SeamWeights(const TreeT& distTree, const Int16TreeT& refSignTree,
+    IntTreeT& refIdxTree, QuantizedPointList& points, double iso)
+    : mDistAcc(distTree)
+    , mRefSignAcc(refSignTree)
+    , mRefIdxAcc(refIdxTree)
     , mPoints(points)
     , mIsovalue(iso)
-    , mRefData(NULL)
 {
 }
 
 
-template <class DistTreeT>
-PointGen<DistTreeT>::PointGen(const PointGen<DistTreeT>& rhs)
-    : mDistTree(rhs.mDistTree)
-    , mAuxLeafs(rhs.mAuxLeafs)
-    , mLeafIndices(rhs.mLeafIndices)
-    , mTransform(rhs.mTransform)
-    , mPoints(rhs.mPoints)
-    , mIsovalue(rhs.mIsovalue)
-    , mRefData(rhs.mRefData)
+template<typename TreeT>
+template <typename LeafNodeType>
+void
+SeamWeights<TreeT>::operator()(LeafNodeType &signLeaf, size_t /*leafIndex*/) const
 {
-}
+    Coord coord = signLeaf.origin();
+    const typename Int16TreeT::LeafNodeType *refSignLeafPt = mRefSignAcc.probeConstLeaf(coord);
 
+    if (!refSignLeafPt) return;
 
-template <class DistTreeT>
-void
-PointGen<DistTreeT>::setRefData(
-    const ReferenceData<DistTreeT>& refData)
-{
-    mRefData = &refData;
-}
+    const typename TreeT::LeafNodeType *distLeafPt = mDistAcc.probeConstLeaf(coord);
+    const typename IntTreeT::LeafNodeType *refIdxLeafPt = mRefIdxAcc.probeConstLeaf(coord);
 
-template <class DistTreeT>
-void
-PointGen<DistTreeT>::runParallel()
-{
-    tbb::parallel_for(mAuxLeafs.getRange(), *this);
-}
+    std::vector<double> values(8);
+    unsigned char lhsSigns, rhsSigns;
+    Vec3d point;
+    Index offset;
 
+    Coord ijk;
+    coord.offset(TreeT::LeafNodeType::DIM - 1);
 
-template <class DistTreeT>
-void
-PointGen<DistTreeT>::runSerial()
-{
-    (*this)(mAuxLeafs.getRange());
-}
+    typename LeafNodeType::ValueOnCIter iter = signLeaf.cbeginValueOn();
+    for (; iter; ++iter) {
 
+        offset = iter.pos();
+        ijk = iter.getCoord();
 
-template <class DistTreeT>
-void
-PointGen<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
-{
-    typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
-    typedef typename DistTreeT::template ValueConverter<Vec3s>::Type Vec3sTreeT;
+        const bool inclusiveCell = ijk[0] < coord[0] && ijk[1] < coord[1] && ijk[2] < coord[2];
 
-    typedef tree::ValueAccessor<BoolTreeT> BoolTreeAccessorT;
-    typedef tree::ValueAccessor<Vec3sTreeT> Vec3sTreeAccessorT;
+        if ((iter.getValue() & SEAM) && refSignLeafPt->isValueOn(offset)) {
 
-    typedef typename BoolTreeT::LeafNodeType BoolLeafT;
-    typedef typename IntTreeT::LeafNodeType IntLeafT;
-    typedef typename Vec3sTreeT::LeafNodeType Vec3sLeafT;
+            lhsSigns = SIGNS & iter.getValue();
+            rhsSigns = SIGNS & refSignLeafPt->getValue(offset);
 
-    boost::scoped_ptr<DistTreeAccessorT> refDistAcc;
-    boost::scoped_ptr<BoolTreeAccessorT> refMaskAcc, refSmoothMaskAcc;
-    boost::scoped_ptr<Vec3sTreeAccessorT> refPtnAcc;
 
-    if (mRefData && mRefData->isValid()) {
-        refDistAcc.reset(new DistTreeAccessorT(*mRefData->mDistTree));
-        refMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mTopologyMaskTree));
-        refSmoothMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mSmoothingMaskTree));
-        refPtnAcc.reset(new Vec3sTreeAccessorT(*mRefData->mSeamPointTree));
-    }
+            if (inclusiveCell) {
+                collectCornerValues(*distLeafPt, offset, values);
+            } else {
+                collectCornerValues(mDistAcc, ijk, values);
+            }
 
 
-    const bool hasRefData = refDistAcc && refMaskAcc;
-    typename IntTreeT::LeafNodeType::ValueOnIter auxIter;
-    DistTreeAccessorT distAcc(mDistTree);
+            for (size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n) {
 
-    Coord ijk;
-    openvdb::Vec3d avg, tmp;
+                int id = matchEdgeGroup(n, lhsSigns, rhsSigns);
 
-    for (size_t n = range.begin(); n != range.end(); ++n) {
+                if (id != -1) {
+
+                    uint32_t& data = mPoints[refIdxLeafPt->getValue(offset) + (id - 1)];
+
+                    if (!(data & MASK_DIRTY_BIT)) {
 
-        size_t idx = mLeafIndices[n];
-        IntLeafT& auxLeaf = *mAuxLeafs[n];
+                        int smaples = computeMaskedPoint(point, values, lhsSigns, rhsSigns, n, mIsovalue);
 
-        BoolLeafT* maskLeaf = NULL;
-        BoolLeafT* smoothMaskLeaf = NULL;
-        Vec3sLeafT* ptnLeaf = NULL;
-        if (hasRefData) {
-            maskLeaf = refMaskAcc->probeLeaf(auxLeaf.getOrigin());
-            smoothMaskLeaf = refSmoothMaskAcc->probeLeaf(auxLeaf.getOrigin());
-            ptnLeaf = refPtnAcc->probeLeaf(auxLeaf.getOrigin());
+                        if (smaples > 0) data = packPoint(point);
+                        else data = MASK_INVALID_BIT;
+
+                        data |= MASK_DIRTY_BIT;
+                    }
+                }
+            }
         }
+    }
+}
 
-        for (auxIter = auxLeaf.beginValueOn(); auxIter; ++auxIter) {
 
-            if(auxIter.getValue() == 0) {
+////////////////////////////////////////
 
-                auxIter.setValue(idx);
-                auxIter.setValueOff();
-                ijk = auxIter.getCoord();
 
-                if (hasRefData && maskLeaf && maskLeaf->isValueOn(ijk)) {
+template <typename TreeT, typename LeafManagerT>
+class MergeVoxelRegions
+{
+public:
+    typedef typename TreeT::ValueType ValueT;
+    typedef tree::ValueAccessor<const TreeT> AccessorT;
 
-                    if (ptnLeaf && ptnLeaf->isValueOn(ijk)) {
-                        avg = ptnLeaf->getValue(ijk);
-                    } else {
-                        int e1 = calcAvgPoint(*refDistAcc.get(), ijk, avg);
+    typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
+    typedef tree::ValueAccessor<IntTreeT> IntAccessorT;
 
-                        if (e1 != (XEDGE|YEDGE|ZEDGE)) {
-                            int e2 = calcAvgPoint(distAcc, ijk, tmp);
-                            if((e2 & (~e1)) != 0) smoothMaskLeaf->setValueOn(ijk);
-                        }
-                    }
-                } else {
-                    calcAvgPoint(distAcc, ijk, avg);
-                }
+    typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
 
-                openvdb::Vec3s& ptn = mPoints[idx];
-                ptn[0] = float(avg[0]);
-                ptn[1] = float(avg[1]);
-                ptn[2] = float(avg[2]);
+    typedef typename LeafManagerT::TreeType::template ValueConverter<Int16>::Type Int16TreeT;
+    typedef tree::ValueAccessor<const Int16TreeT> Int16AccessorT;
 
-                ++idx;
-            }
-        }
+    typedef typename TreeT::template ValueConverter<float>::Type FloatTreeT;
+    typedef Grid<FloatTreeT> FloatGridT;
 
-        while(auxLeaf.onVoxelCount() > 0) {
 
-            avg[0] = 0;
-            avg[1] = 0;
-            avg[2] = 0;
+    //////////
 
-            auxIter = auxLeaf.beginValueOn();
-            int regionId = auxIter.getValue(), points = 0;
+    MergeVoxelRegions(const LeafManagerT& signLeafs, const Int16TreeT& signTree,
+        const TreeT& distTree, IntTreeT& idxTree, ValueT iso, ValueT adaptivity);
 
-            for (; auxIter; ++auxIter) {
-                if(auxIter.getValue() == regionId) {
+    void run(bool threaded = true);
 
-                    auxIter.setValue(idx);
-                    auxIter.setValueOff();
-                    ijk = auxIter.getCoord();
+    void setSpatialAdaptivity(
+        const math::Transform& distGridXForm, const FloatGridT& adaptivityField);
 
-                    if (hasRefData && maskLeaf && maskLeaf->isValueOn(ijk)) {
-                        calcAvgPoint(*refDistAcc.get(), ijk, tmp);
-                    } else {
-                        calcAvgPoint(distAcc, ijk, tmp);
-                    }
+    void setAdaptivityMask(const BoolTreeT* mask);
 
-                    avg += tmp;
-                    ++points;
-                }
-            }
+    void setRefData(const Int16TreeT* signTree, ValueT adaptivity);
 
-            if (points > 1) {
-                double w = 1.0 / double(points);
-                avg[0] *= w;
-                avg[1] *= w;
-                avg[2] *= w;
-            }
+    //////////
 
-            openvdb::Vec3s& ptn = mPoints[idx];
-            ptn[0] = float(avg[0]);
-            ptn[1] = float(avg[1]);
-            ptn[2] = float(avg[2]);
-            ++idx;
-        }
-    }
+
+    void operator()(const tbb::blocked_range<size_t>&) const;
+
+private:
+
+    const LeafManagerT& mSignLeafs;
+
+    const Int16TreeT& mSignTree;
+    Int16AccessorT mSignAcc;
+
+    const TreeT& mDistTree;
+    AccessorT mDistAcc;
+
+    IntTreeT& mIdxTree;
+    ValueT mIsovalue, mSurfaceAdaptivity, mInternalAdaptivity;
+
+    const math::Transform* mTransform;
+    const FloatGridT* mAdaptivityGrid;
+    const BoolTreeT* mMask;
+
+    const Int16TreeT* mRefSignTree;
+};
+
+template <typename TreeT, typename LeafManagerT>
+MergeVoxelRegions<TreeT, LeafManagerT>::MergeVoxelRegions(
+    const LeafManagerT& signLeafs, const Int16TreeT& signTree,
+    const TreeT& distTree, IntTreeT& idxTree, ValueT iso, ValueT adaptivity)
+    : mSignLeafs(signLeafs)
+    , mSignTree(signTree)
+    , mSignAcc(mSignTree)
+    , mDistTree(distTree)
+    , mDistAcc(mDistTree)
+    , mIdxTree(idxTree)
+    , mIsovalue(iso)
+    , mSurfaceAdaptivity(adaptivity)
+    , mInternalAdaptivity(adaptivity)
+    , mTransform(NULL)
+    , mAdaptivityGrid(NULL)
+    , mMask(NULL)
+    , mRefSignTree(NULL)
+{
 }
 
-template <class DistTreeT>
-int
-PointGen<DistTreeT>::calcAvgPoint(DistTreeAccessorT& acc,
-    const Coord& ijk, openvdb::Vec3d& avg) const
+
+template <typename TreeT, typename LeafManagerT>
+void
+MergeVoxelRegions<TreeT, LeafManagerT>::run(bool threaded)
 {
-    double values[8];
-    bool signMask[8];
-    Coord coord;
+    if (threaded) tbb::parallel_for(mSignLeafs.getRange(), *this);
+    else (*this)(mSignLeafs.getRange());
+}
 
-    // Sample corner values
-    coord = ijk;
-    values[0] = double(acc.getValue(coord)); // i, j, k
 
-    coord[0] += 1;
-    values[1] = double(acc.getValue(coord)); // i+1, j, k
+template <typename TreeT, typename LeafManagerT>
+void
+MergeVoxelRegions<TreeT, LeafManagerT>::setSpatialAdaptivity(
+    const math::Transform& distGridXForm, const FloatGridT& adaptivityField)
+{
+    mTransform = &distGridXForm;
+    mAdaptivityGrid = &adaptivityField;
+}
 
-    coord[2] += 1;
-    values[2] = double(acc.getValue(coord)); // i+i, j, k+1
 
-    coord[0] = ijk[0];
-    values[3] = double(acc.getValue(coord)); // i, j, k+1
+template <typename TreeT, typename LeafManagerT>
+void
+MergeVoxelRegions<TreeT, LeafManagerT>::setAdaptivityMask(const BoolTreeT* mask)
+{
+    mMask = mask;
+}
 
-    coord[1] += 1; coord[2] = ijk[2];
-    values[4] = double(acc.getValue(coord)); // i, j+1, k
+template <typename TreeT, typename LeafManagerT>
+void
+MergeVoxelRegions<TreeT, LeafManagerT>::setRefData(const Int16TreeT* signTree, ValueT adaptivity)
+{
+    mRefSignTree = signTree;
+    mInternalAdaptivity = adaptivity;
+}
 
-    coord[0] += 1;
-    values[5] = double(acc.getValue(coord)); // i+1, j+1, k
 
-    coord[2] += 1;
-    values[6] = double(acc.getValue(coord)); // i+1, j+1, k+1
+template <typename TreeT, typename LeafManagerT>
+void
+MergeVoxelRegions<TreeT, LeafManagerT>::operator()(const tbb::blocked_range<size_t>& range) const
+{
+    typedef math::Vec3<ValueT> Vec3T;
 
-    coord[0] = ijk[0];
-    values[7] = double(acc.getValue(coord)); // i, j+1, k+1
+    typedef typename TreeT::LeafNodeType LeafT;
+    typedef typename IntTreeT::LeafNodeType IntLeafT;
+    typedef typename BoolTreeT::LeafNodeType BoolLeafT;
+    typedef typename LeafT::template ValueConverter<Vec3T>::Type Vec3LeafT;
 
-    // init sign mask
-    for (int n = 0; n < 8; ++n) signMask[n] = (values[n] < mIsovalue);
+    const int LeafDim = LeafT::DIM;
 
-    int samples = 0, edgeFlags = 0;
-    avg[0] = 0.0;
-    avg[1] = 0.0;
-    avg[2] = 0.0;
+    IntAccessorT idxAcc(mIdxTree);
 
-    if (signMask[0] != signMask[1]) { // Edged: 0 - 1
-        avg[0] += root(values[0], values[1]);
-        ++samples;
-        edgeFlags |= XEDGE;
-    }
+    typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
 
-    if (signMask[1] != signMask[2]) { // Edged: 1 - 2
-        avg[0] += 1.0;
-        avg[2] += root(values[1], values[2]);
-        ++samples;
-        edgeFlags |= ZEDGE;
+    typedef typename tree::ValueAccessor<const FloatTreeT> FloatTreeCAccessorT;
+    boost::scoped_ptr<FloatTreeCAccessorT> adaptivityAcc;
+    if (mAdaptivityGrid) {
+        adaptivityAcc.reset(new FloatTreeCAccessorT(mAdaptivityGrid->tree()));
     }
 
-    if (signMask[3] != signMask[2]) { // Edged: 3 - 2
-        avg[0] += root(values[3], values[2]);
-        avg[2] += 1.0;
-        ++samples;
-        edgeFlags |= XEDGE;
+    typedef typename tree::ValueAccessor<const Int16TreeT> Int16TreeCAccessorT;
+    boost::scoped_ptr<Int16TreeCAccessorT> refAcc;
+    if (mRefSignTree) {
+        refAcc.reset(new Int16TreeCAccessorT(*mRefSignTree));
     }
 
-    if (signMask[0] != signMask[3]) { // Edged: 0 - 3
-        avg[2] += root(values[0], values[3]);
-        ++samples;
-        edgeFlags |= ZEDGE;
+    typedef typename tree::ValueAccessor<const BoolTreeT> BoolTreeCAccessorT;
+    boost::scoped_ptr<BoolTreeCAccessorT> maskAcc;
+    if (mMask) {
+        maskAcc.reset(new BoolTreeCAccessorT(*mMask));
     }
 
-    if (signMask[4] != signMask[5]) { // Edged: 4 - 5
-        avg[0] += root(values[4], values[5]);
-        avg[1] += 1.0;
-        ++samples;
-        edgeFlags |= XEDGE;
-    }
+    // Allocate reusable leaf buffers
+    BoolLeafT mask;
+    Vec3LeafT gradientBuffer;
+    Coord ijk, nijk, coord, end;
 
-    if (signMask[5] != signMask[6]) { // Edged: 5 - 6
-        avg[0] += 1.0;
-        avg[1] += 1.0;
-        avg[2] += root(values[5], values[6]);
-        ++samples;
-        edgeFlags |= ZEDGE;
-    }
+    for (size_t n = range.begin(); n != range.end(); ++n) {
 
-    if (signMask[7] != signMask[6]) { // Edged: 7 - 6
-        avg[0] += root(values[7], values[6]);
-        avg[1] += 1.0;
-        avg[2] += 1.0;
-        ++samples;
-        edgeFlags |= XEDGE;
-    }
+        const Coord& origin = mSignLeafs.leaf(n).origin();
 
-    if (signMask[4] != signMask[7]) { // Edged: 4 - 7
-        avg[1] += 1.0;
-        avg[2] += root(values[4], values[7]);
-        ++samples;
-        edgeFlags |= ZEDGE;
-    }
+        ValueT adaptivity = mSurfaceAdaptivity;
 
-    if (signMask[0] != signMask[4]) { // Edged: 0 - 4
-        avg[1] += root(values[0], values[4]);
-        ++samples;
-        edgeFlags |= YEDGE;
-    }
+        if (refAcc && refAcc->probeConstLeaf(origin) == NULL) {
+            adaptivity = mInternalAdaptivity;
+        }
 
-    if (signMask[1] != signMask[5]) { // Edged: 1 - 5
-        avg[0] += 1.0;
-        avg[1] += root(values[1], values[5]);
-        ++samples;
-        edgeFlags |= YEDGE;
-    }
+        IntLeafT& idxLeaf = *idxAcc.probeLeaf(origin);
 
-    if (signMask[2] != signMask[6]) { // Edged: 2 - 6
-        avg[0] += 1.0;
-        avg[1] += root(values[2], values[6]);
-        avg[2] += 1.0;
-        ++samples;
-        edgeFlags |= YEDGE;
-    }
+        end[0] = origin[0] + LeafDim;
+        end[1] = origin[1] + LeafDim;
+        end[2] = origin[2] + LeafDim;
 
-    if (signMask[3] != signMask[7]) { // Edged: 3 - 7
-        avg[1] += root(values[3], values[7]);
-        avg[2] += 1.0;
-        ++samples;
-        edgeFlags |= YEDGE;
-    }
+        mask.setValuesOff();
+
+        // Mask off seam line adjacent voxels
+        if (maskAcc) {
+            const BoolLeafT* maskLeaf = maskAcc->probeConstLeaf(origin);
+            if (maskLeaf != NULL) {
+                typename BoolLeafT::ValueOnCIter it;
+                for (it = maskLeaf->cbeginValueOn(); it; ++it) {
+                    ijk = it.getCoord();
+                    coord[0] = ijk[0] - (ijk[0] % 2);
+                    coord[1] = ijk[1] - (ijk[1] % 2);
+                    coord[2] = ijk[2] - (ijk[2] % 2);
+                    mask.setActiveState(coord, true);
+                }
+            }
+        }
 
-    if (samples > 1) {
-        double w = 1.0 / double(samples);
-        avg[0] *= w;
-        avg[1] *= w;
-        avg[2] *= w;
-    }
 
-    // offset by cell-origin
-    avg[0] += double(ijk[0]);
-    avg[1] += double(ijk[1]);
-    avg[2] += double(ijk[2]);
+        LeafT adaptivityLeaf(origin, adaptivity);
 
-    avg = mTransform.indexToWorld(avg);
+        if (mAdaptivityGrid) {
+            for (Index offset = 0; offset < LeafT::NUM_VALUES; ++offset) {
+                ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
+                Vec3d xyz = mAdaptivityGrid->transform().worldToIndex(
+                    mTransform->indexToWorld(ijk));
+                ValueT tmpA = ValueT(adaptivityAcc->getValue(util::nearestCoord(xyz)));
+                adaptivityLeaf.setValueOnly(offset, tmpA * adaptivity);
+            }
+        }
 
-    return edgeFlags;
+        // Mask off ambiguous voxels
+        for (iter = mSignLeafs.leaf(n).cbeginValueOn(); iter; ++iter) {
+            ijk = iter.getCoord();
+            coord[0] = ijk[0] - (ijk[0] % 2);
+            coord[1] = ijk[1] - (ijk[1] % 2);
+            coord[2] = ijk[2] - (ijk[2] % 2);
+            if(mask.isValueOn(coord)) continue;
+
+
+
+            int flags = int(iter.getValue());
+            unsigned char signs = SIGNS & flags;
+            if ((flags & SEAM) || !sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
+                mask.setActiveState(coord, true);
+                continue;
+            }
+
+            for (int i = 0; i < 26; ++i) {
+                nijk = ijk + util::COORD_OFFSETS[i];
+                signs = SIGNS & mSignAcc.getValue(nijk);
+                if (!sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
+                    mask.setActiveState(coord, true);
+                    break;
+                }
+            }
+        }
+
+        int dim = 2;
+        // Mask off topologically ambiguous 2x2x2 voxel sub-blocks
+        for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
+            for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
+                for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
+                    if (isNonManifold(mDistAcc, ijk, mIsovalue, dim)) {
+                        mask.setActiveState(ijk, true);
+                    }
+                }
+            }
+        }
+
+        // Compute the gradient for the remaining voxels
+        gradientBuffer.setValuesOff();
+        for (iter = mSignLeafs.leaf(n).cbeginValueOn(); iter; ++iter) {
+
+            ijk = iter.getCoord();
+            coord[0] = ijk[0] - (ijk[0] % dim);
+            coord[1] = ijk[1] - (ijk[1] % dim);
+            coord[2] = ijk[2] - (ijk[2] % dim);
+            if(mask.isValueOn(coord)) continue;
+
+            Vec3T norm(math::ISGradient<math::CD_2ND>::result(mDistAcc, ijk));
+            // Normalize (Vec3's normalize uses isApproxEqual, which uses abs and does more work)
+            ValueT length = norm.length();
+            if (length > ValueT(1.0e-7)) {
+                norm *= ValueT(1.0) / length;
+            }
+            gradientBuffer.setValue(ijk, norm);
+        }
+
+        int regionId = 1, next_dim = dim << 1;
+
+        // Process the first adaptivity level.
+        for (ijk[0] = 0; ijk[0] < LeafDim; ijk[0] += dim) {
+            coord[0] = ijk[0] - (ijk[0] % next_dim);
+            for (ijk[1] = 0; ijk[1] < LeafDim; ijk[1] += dim) {
+                coord[1] = ijk[1] - (ijk[1] % next_dim);
+                for (ijk[2] = 0; ijk[2] < LeafDim; ijk[2] += dim) {
+                    coord[2] = ijk[2] - (ijk[2] % next_dim);
+                    adaptivity = adaptivityLeaf.getValue(ijk);
+                    if (mask.isValueOn(ijk) || !isMergable(gradientBuffer, ijk, dim, adaptivity)) {
+                        mask.setActiveState(coord, true);
+                        continue;
+                    }
+                    mergeVoxels(idxLeaf, ijk, dim, regionId++);
+                }
+            }
+        }
+
+
+        // Process remaining adaptivity levels
+        for (dim = 4; dim < LeafDim; dim = dim << 1) {
+            next_dim = dim << 1;
+            coord[0] = ijk[0] - (ijk[0] % next_dim);
+            for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
+                coord[1] = ijk[1] - (ijk[1] % next_dim);
+                for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
+                    coord[2] = ijk[2] - (ijk[2] % next_dim);
+                    for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
+                        adaptivity = adaptivityLeaf.getValue(ijk);
+                        if (mask.isValueOn(ijk) || isNonManifold(mDistAcc, ijk, mIsovalue, dim) ||
+                            !isMergable(gradientBuffer, ijk, dim, adaptivity))
+                        {
+                            mask.setActiveState(coord, true);
+                            continue;
+                        }
+                        mergeVoxels(idxLeaf, ijk, dim, regionId++);
+                    }
+                }
+            }
+        }
+
+        adaptivity = adaptivityLeaf.getValue(origin);
+        if (!(mask.isValueOn(origin) || isNonManifold(mDistAcc, origin, mIsovalue, LeafDim))
+            && isMergable(gradientBuffer, origin, LeafDim, adaptivity))
+        {
+            mergeVoxels(idxLeaf, origin, LeafDim, regionId++);
+        }
+    }
 }
 
 
 ////////////////////////////////////////
 
 
-struct QuadMeshOp
+// Constructs qudas
+struct UniformPrimBuilder
 {
-    QuadMeshOp(): mIdx(0), mPolygonPool(NULL) {}
+    UniformPrimBuilder(): mIdx(0), mPolygonPool(NULL) {}
 
     void init(const size_t upperBound, PolygonPool& quadPool)
     {
@@ -1331,7 +2353,10 @@ struct QuadMeshOp
         ++mIdx;
     }
 
-    void done() {}
+    void done()
+    {
+        mPolygonPool->trimQuads(mIdx);
+    }
 
 private:
     size_t mIdx;
@@ -1339,16 +2364,16 @@ private:
 };
 
 
-struct AdaptiveMeshOp
+// Constructs qudas and triangles
+struct AdaptivePrimBuilder
 {
-    AdaptiveMeshOp(): mQuadIdx(0), mTriangleIdx(0), mPolygonPool(NULL), mTmpPolygonPool() {}
+    AdaptivePrimBuilder() : mQuadIdx(0), mTriangleIdx(0), mPolygonPool(NULL) {}
 
     void init(const size_t upperBound, PolygonPool& polygonPool)
     {
         mPolygonPool = &polygonPool;
-
-        mTmpPolygonPool.resetQuads(upperBound);
-        mTmpPolygonPool.resetTriangles(upperBound);
+        mPolygonPool->resetQuads(upperBound);
+        mPolygonPool->resetTriangles(upperBound);
 
         mQuadIdx = 0;
         mTriangleIdx = 0;
@@ -1358,35 +2383,35 @@ struct AdaptiveMeshOp
     {
         if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
             && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
-            mTmpPolygonPool.quadFlags(mQuadIdx) = flags;
+            mPolygonPool->quadFlags(mQuadIdx) = flags;
             addQuad(verts, reverse);
         } else if (
             verts[0] == verts[3] &&
             verts[1] != verts[2] &&
             verts[1] != verts[0] &&
             verts[2] != verts[0]) {
-            mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
+            mPolygonPool->triangleFlags(mTriangleIdx) = flags;
             addTriangle(verts[0], verts[1], verts[2], reverse);
         } else if (
             verts[1] == verts[2] &&
             verts[0] != verts[3] &&
             verts[0] != verts[1] &&
             verts[3] != verts[1]) {
-            mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
+            mPolygonPool->triangleFlags(mTriangleIdx) = flags;
             addTriangle(verts[0], verts[1], verts[3], reverse);
         } else if (
             verts[0] == verts[1] &&
             verts[2] != verts[3] &&
             verts[2] != verts[0] &&
             verts[3] != verts[0]) {
-            mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
+            mPolygonPool->triangleFlags(mTriangleIdx) = flags;
             addTriangle(verts[0], verts[2], verts[3], reverse);
         } else if (
             verts[2] == verts[3] &&
             verts[0] != verts[1] &&
             verts[0] != verts[2] &&
             verts[1] != verts[2]) {
-            mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
+            mPolygonPool->triangleFlags(mTriangleIdx) = flags;
             addTriangle(verts[0], verts[1], verts[2], reverse);
         }
     }
@@ -1394,19 +2419,8 @@ struct AdaptiveMeshOp
 
     void done()
     {
-        mPolygonPool->resetQuads(mQuadIdx);
-        for (size_t  i = 0; i < mQuadIdx; ++i) {
-            mPolygonPool->quad(i) = mTmpPolygonPool.quad(i);
-            mPolygonPool->quadFlags(i) = mTmpPolygonPool.quadFlags(i);
-        }
-        mTmpPolygonPool.clearQuads();
-
-        mPolygonPool->resetTriangles(mTriangleIdx);
-        for (size_t  i = 0; i < mTriangleIdx; ++i) {
-            mPolygonPool->triangle(i) = mTmpPolygonPool.triangle(i);
-            mPolygonPool->triangleFlags(i) = mTmpPolygonPool.triangleFlags(i);
-        }
-        mTmpPolygonPool.clearTriangles();
+        mPolygonPool->trimQuads(mQuadIdx, /*reallocate=*/true);
+        mPolygonPool->trimTrinagles(mTriangleIdx, /*reallocate=*/true);
     }
 
 private:
@@ -1414,9 +2428,9 @@ private:
     void addQuad(const Vec4I& verts, bool reverse)
     {
         if (!reverse) {
-            mTmpPolygonPool.quad(mQuadIdx) = verts;
+            mPolygonPool->quad(mQuadIdx) = verts;
         } else {
-            Vec4I& quad = mTmpPolygonPool.quad(mQuadIdx);
+            Vec4I& quad = mPolygonPool->quad(mQuadIdx);
             quad[0] = verts[3];
             quad[1] = verts[2];
             quad[2] = verts[1];
@@ -1427,7 +2441,7 @@ private:
 
     void addTriangle(unsigned v0, unsigned v1, unsigned v2, bool reverse)
     {
-        Vec3I& prim = mTmpPolygonPool.triangle(mTriangleIdx);
+        Vec3I& prim = mPolygonPool->triangle(mTriangleIdx);
 
         prim[1] = v1;
 
@@ -1443,206 +2457,295 @@ private:
 
     size_t mQuadIdx, mTriangleIdx;
     PolygonPool *mPolygonPool;
-    PolygonPool mTmpPolygonPool;
 };
 
 
-////////////////////////////////////////
+template<typename SignAccT, typename IdxAccT, typename PrimBuilder>
+inline void
+constructPolygons(Int16 flags, Int16 refFlags, const Vec4i& offsets, const Coord& ijk,
+    const SignAccT& signAcc, const IdxAccT& idxAcc, PrimBuilder& mesher, Index32 pointListSize)
+{
+    const Index32 v0 = idxAcc.getValue(ijk);
+    if (v0 == util::INVALID_IDX) return;
 
+    char tag[2];
+    tag[0] = (flags & SEAM) ? POLYFLAG_FRACTURE_SEAM : 0;
+    tag[1] = tag[0] | char(POLYFLAG_EXTERIOR);
 
-template<class DistTreeT, class MeshingOp>
-class MeshGen
-{
-public:
-    typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
-    typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
+    const bool isInside = flags & INSIDE;
 
-    MeshGen(const LeafCPtrList<CharTreeT>& edgeLeafs, const IntTreeT& auxTree, PolygonPoolList&);
-    MeshGen(const MeshGen<DistTreeT, MeshingOp>&);
+    Coord coord;
+    openvdb::Vec4I quad;
+    unsigned char cell;
+    Index32 tmpIdx = 0;
 
-    void setRefData(const ReferenceData<DistTreeT>&);
+    if (flags & XEDGE) {
 
-    void runParallel();
-    void runSerial();
+        quad[0] = v0 + offsets[0];
 
-    void operator()(const tbb::blocked_range<size_t>&) const;
+        // i, j-1, k
+        coord[0] = ijk[0];
+        coord[1] = ijk[1] - 1;
+        coord[2] = ijk[2];
 
-private:
-    const LeafCPtrList<CharTreeT>& mEdgeLeafs;
-    const IntTreeT& mAuxTree;
-    const PolygonPoolList& mPolygonPoolList;
-    size_t mID;
-    const ReferenceData<DistTreeT>* mRefData;
-};
+        quad[1] = idxAcc.getValue(coord);
+        cell = SIGNS & signAcc.getValue(coord);
+        if (sEdgeGroupTable[cell][0] > 1) {
+            tmpIdx = quad[1] + Index32(sEdgeGroupTable[cell][5] - 1);
+            if (tmpIdx < pointListSize) quad[1] = tmpIdx;
+        }
 
+        // i, j-1, k-1
+        coord[2] -= 1;
 
-template<class DistTreeT, class MeshingOp>
-MeshGen<DistTreeT, MeshingOp>::MeshGen(const LeafCPtrList<CharTreeT>& edgeLeafs,
-    const IntTreeT& auxTree, PolygonPoolList& polygonPoolList)
-    : mEdgeLeafs(edgeLeafs)
-    , mAuxTree(auxTree)
-    , mPolygonPoolList(polygonPoolList)
-    , mRefData(NULL)
-{
-}
+        quad[2] = idxAcc.getValue(coord);
+        cell = SIGNS & signAcc.getValue(coord);
+        if (sEdgeGroupTable[cell][0] > 1) {
+            tmpIdx = quad[2] + Index32(sEdgeGroupTable[cell][7] - 1);
+            if (tmpIdx < pointListSize) quad[2] = tmpIdx;
+        }
 
+        // i, j, k-1
+        coord[1] = ijk[1];
 
-template<class DistTreeT, class MeshingOp>
-MeshGen<DistTreeT, MeshingOp>::MeshGen(const MeshGen<DistTreeT, MeshingOp>& rhs)
-    : mEdgeLeafs(rhs.mEdgeLeafs)
-    , mAuxTree(rhs.mAuxTree)
-    , mPolygonPoolList(rhs.mPolygonPoolList)
-    , mRefData(rhs.mRefData)
-{
+        quad[3] = idxAcc.getValue(coord);
+        cell = SIGNS & signAcc.getValue(coord);
+        if (sEdgeGroupTable[cell][0] > 1) {
+            tmpIdx = quad[3] + Index32(sEdgeGroupTable[cell][3] - 1);
+            if (tmpIdx < pointListSize) quad[3] = tmpIdx;
+        }
+
+        if (quad[1] != util::INVALID_IDX &&
+            quad[2] != util::INVALID_IDX && quad[3] != util::INVALID_IDX) {
+            mesher.addPrim(quad, isInside, tag[bool(refFlags & XEDGE)]);
+        }
+    }
+
+
+    if (flags & YEDGE) {
+
+        quad[0] = v0 + offsets[1];
+
+        // i, j, k-1
+        coord[0] = ijk[0];
+        coord[1] = ijk[1];
+        coord[2] = ijk[2] - 1;
+
+        quad[1] = idxAcc.getValue(coord);
+        cell = SIGNS & signAcc.getValue(coord);
+        if (sEdgeGroupTable[cell][0] > 1) {
+            tmpIdx = quad[1] + Index32(sEdgeGroupTable[cell][12] - 1);
+            if (tmpIdx < pointListSize) quad[1] = tmpIdx;
+        }
+
+        // i-1, j, k-1
+        coord[0] -= 1;
+
+        quad[2] = idxAcc.getValue(coord);
+        cell = SIGNS & signAcc.getValue(coord);
+        if (sEdgeGroupTable[cell][0] > 1) {
+            tmpIdx = quad[2] + Index32(sEdgeGroupTable[cell][11] - 1);
+            if (tmpIdx < pointListSize) quad[2] = tmpIdx;
+        }
+
+        // i-1, j, k
+        coord[2] = ijk[2];
+
+        quad[3] = idxAcc.getValue(coord);
+        cell = SIGNS & signAcc.getValue(coord);
+        if (sEdgeGroupTable[cell][0] > 1) {
+            tmpIdx = quad[3] + Index32(sEdgeGroupTable[cell][10] - 1);
+            if (tmpIdx < pointListSize) quad[3] = tmpIdx;
+        }
+
+        if (quad[1] != util::INVALID_IDX &&
+            quad[2] != util::INVALID_IDX && quad[3] != util::INVALID_IDX) {
+            mesher.addPrim(quad, isInside, tag[bool(refFlags & YEDGE)]);
+        }
+    }
+
+    if (flags & ZEDGE) {
+
+        quad[0] = v0 + offsets[2];
+
+        // i, j-1, k
+        coord[0] = ijk[0];
+        coord[1] = ijk[1] - 1;
+        coord[2] = ijk[2];
+
+        quad[1] = idxAcc.getValue(coord);
+        cell = SIGNS & signAcc.getValue(coord);
+        if (sEdgeGroupTable[cell][0] > 1) {
+            tmpIdx = quad[1] + Index32(sEdgeGroupTable[cell][8] - 1);
+            if (tmpIdx < pointListSize) quad[1] = tmpIdx;
+        }
+
+        // i-1, j-1, k
+        coord[0] -= 1;
+
+        quad[2] = idxAcc.getValue(coord);
+        cell = SIGNS & signAcc.getValue(coord);
+        if (sEdgeGroupTable[cell][0] > 1) {
+            tmpIdx = quad[2] + Index32(sEdgeGroupTable[cell][6] - 1);
+            if (tmpIdx < pointListSize) quad[2] = tmpIdx;
+        }
+
+        // i-1, j, k
+        coord[1] = ijk[1];
+
+        quad[3] = idxAcc.getValue(coord);
+        cell = SIGNS & signAcc.getValue(coord);
+        if (sEdgeGroupTable[cell][0] > 1) {
+            tmpIdx = quad[3] + Index32(sEdgeGroupTable[cell][2] - 1);
+            if (tmpIdx < pointListSize) quad[3] = tmpIdx;
+        }
+
+        if (quad[1] != util::INVALID_IDX &&
+            quad[2] != util::INVALID_IDX && quad[3] != util::INVALID_IDX) {
+            mesher.addPrim(quad, !isInside, tag[bool(refFlags & ZEDGE)]);
+        }
+    }
 }
 
 
-template<class DistTreeT, class MeshingOp>
-void
-MeshGen<DistTreeT, MeshingOp>::setRefData(
-    const ReferenceData<DistTreeT>& refData)
+////////////////////////////////////////
+
+
+template<typename LeafManagerT, typename PrimBuilder>
+class GenPolygons
 {
-    mRefData = &refData;
-}
+public:
+    typedef typename LeafManagerT::TreeType::template ValueConverter<int>::Type IntTreeT;
+    typedef typename LeafManagerT::TreeType::template ValueConverter<Int16>::Type Int16TreeT;
 
+    typedef tree::ValueAccessor<const IntTreeT> IntAccessorT;
+    typedef tree::ValueAccessor<const Int16TreeT> Int16AccessorT;
 
-template<class DistTreeT, class MeshingOp>
-void
-MeshGen<DistTreeT, MeshingOp>::runParallel()
+    //////////
+
+
+    GenPolygons(const LeafManagerT& signLeafs, const Int16TreeT& signTree,
+        const IntTreeT& idxTree, PolygonPoolList& polygons, Index32 pointListSize);
+
+    void run(bool threaded = true);
+
+
+    void setRefSignTree(const Int16TreeT *r) { mRefSignTree = r; }
+
+    //////////
+
+
+    void operator()(const tbb::blocked_range<size_t>&) const;
+
+private:
+    const LeafManagerT& mSignLeafs;
+    const Int16TreeT& mSignTree;
+    const IntTreeT& mIdxTree;
+    const PolygonPoolList& mPolygonPoolList;
+    const Index32 mPointListSize;
+
+    const Int16TreeT *mRefSignTree;
+ };
+
+
+template<typename LeafManagerT, typename PrimBuilder>
+GenPolygons<LeafManagerT, PrimBuilder>::GenPolygons(const LeafManagerT& signLeafs,
+    const Int16TreeT& signTree, const IntTreeT& idxTree, PolygonPoolList& polygons,
+    Index32 pointListSize)
+    : mSignLeafs(signLeafs)
+    , mSignTree(signTree)
+    , mIdxTree(idxTree)
+    , mPolygonPoolList(polygons)
+    , mPointListSize(pointListSize)
+    , mRefSignTree(NULL)
 {
-    tbb::parallel_for(mEdgeLeafs.getRange(), *this);
 }
 
-
-template<class DistTreeT, class MeshingOp>
+template<typename LeafManagerT, typename PrimBuilder>
 void
-MeshGen<DistTreeT, MeshingOp>::runSerial()
+GenPolygons<LeafManagerT, PrimBuilder>::run(bool threaded)
 {
-    (*this)(mEdgeLeafs.getRange());
+    if (threaded) tbb::parallel_for(mSignLeafs.getRange(), *this);
+    else (*this)(mSignLeafs.getRange());
 }
 
-
-template<class DistTreeT, class MeshingOp>
+template<typename LeafManagerT, typename PrimBuilder>
 void
-MeshGen<DistTreeT, MeshingOp>::operator()(
+GenPolygons<LeafManagerT, PrimBuilder>::operator()(
     const tbb::blocked_range<size_t>& range) const
 {
-    typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
-    typedef typename BoolTreeT::LeafNodeType BoolLeafT;
-    typedef typename CharTreeT::LeafNodeType CharLeafT;
+    typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
+    IntAccessorT idxAcc(mIdxTree);
+    Int16AccessorT signAcc(mSignTree);
 
-    typedef openvdb::tree::ValueAccessor<const CharTreeT> CharTreeAccessorT;
-    typedef openvdb::tree::ValueAccessor<const IntTreeT> IntTreeAccessorT;
-    typedef openvdb::tree::ValueAccessor<const BoolTreeT> BoolTreeAccessorT;
 
-    boost::scoped_ptr<CharTreeAccessorT> refEdgeAcc;
-    boost::scoped_ptr<BoolTreeAccessorT> refMaskAcc;
-    if (mRefData && mRefData->isValid()) {
-        refEdgeAcc.reset(new CharTreeAccessorT(*mRefData->mEdgeTree));
-        refMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mSeamMaskTree.get()));
-    }
-    const bool hasRefData = refEdgeAcc && refMaskAcc;
+    PrimBuilder mesher;
+    size_t edgeCount;
+    Coord ijk, origin;
 
 
-    typename CharTreeT::LeafNodeType::ValueOnCIter iter;
-    IntTreeAccessorT auxAcc(mAuxTree);
+    // reference data
+    boost::scoped_ptr<Int16AccessorT> refSignAcc;
+    if (mRefSignTree) refSignAcc.reset(new Int16AccessorT(*mRefSignTree));
 
-    Coord ijk, coord;
-    char refEdgeFlags, isSemLinePoly;
-    const char isExteriorPoly[2] = {0, char(POLYFLAG_EXTERIOR)};
-    openvdb::Vec4I quad;
-    size_t edgeCount;
-
-    MeshingOp mesher;
 
     for (size_t n = range.begin(); n != range.end(); ++n) {
 
-        const Coord origin = mEdgeLeafs[n]->getOrigin();
+        origin = mSignLeafs.leaf(n).origin();
 
         // Get an upper bound on the number of primitives.
         edgeCount = 0;
-        iter = mEdgeLeafs[n]->cbeginValueOn();
+        iter = mSignLeafs.leaf(n).cbeginValueOn();
         for (; iter; ++iter) {
-            char edgeFlags = iter.getValue() >> 1;
-            edgeCount += edgeFlags & 0x1;
-
-            edgeFlags = edgeFlags >> 1;
-            edgeCount += edgeFlags & 0x1;
-
-            edgeFlags = edgeFlags >> 1;
-            edgeCount += edgeFlags & 0x1;
+            if (iter.getValue() & XEDGE) ++edgeCount;
+            if (iter.getValue() & YEDGE) ++edgeCount;
+            if (iter.getValue() & ZEDGE) ++edgeCount;
         }
 
+        if(edgeCount == 0) continue;
+
         mesher.init(edgeCount, mPolygonPoolList[n]);
 
-        const CharLeafT* refEdgeLeaf = NULL;
-        const BoolLeafT* refMaskLeaf = NULL;
+        const typename Int16TreeT::LeafNodeType *signleafPt = signAcc.probeConstLeaf(origin);
+        const typename IntTreeT::LeafNodeType *idxLeafPt = idxAcc.probeConstLeaf(origin);
 
-        if (hasRefData) {
-            refEdgeLeaf = refEdgeAcc->probeConstLeaf(origin);
-            refMaskLeaf = refMaskAcc->probeConstLeaf(origin);
-        }
+        if (!signleafPt || !idxLeafPt) continue;
 
-        iter = mEdgeLeafs[n]->cbeginValueOn();
-        for (; iter; ++iter) {
-            ijk = iter.getCoord();
-            const char& edgeFlags = iter.getValue();
 
-            const bool isInside = edgeFlags & INSIDE;
-
-            refEdgeFlags = 0;
-            isSemLinePoly = 0;
-            if (hasRefData) {
-                if(refEdgeLeaf) refEdgeFlags = refEdgeLeaf->getValue(ijk);
-                if (refMaskLeaf && refMaskLeaf->isValueOn(ijk)) {
-                    isSemLinePoly = char(POLYFLAG_FRACTURE_SEAM);
-                }
-            }
+        const typename Int16TreeT::LeafNodeType *refSignLeafPt = NULL;
+        if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
 
+        Vec4i offsets;
 
-            int v0 = auxAcc.getValue(ijk);
+        iter = mSignLeafs.leaf(n).cbeginValueOn();
+        for (; iter; ++iter) {
+            ijk = iter.getCoord();
 
-            if (edgeFlags & XEDGE) {
+            Int16 flags = iter.getValue();
 
-                quad[0] = v0;
-                coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]; // i, j-1, k
-                quad[1] = auxAcc.getValue(coord);
-                coord[2] -= 1; // i, j-1, k-1
-                quad[2] = auxAcc.getValue(coord);
-                coord[1] = ijk[1]; // i, j, k-1
-                quad[3] = auxAcc.getValue(coord);
+            if (!(flags & 0xE00)) continue;
 
-                mesher.addPrim(quad, isInside,
-                    (isSemLinePoly | isExteriorPoly[bool(refEdgeFlags & XEDGE)]));
+            Int16 refFlags = 0;
+            if (refSignLeafPt) {
+                refFlags = refSignLeafPt->getValue(iter.pos());
             }
 
+            offsets[0] = 0;
+            offsets[1] = 0;
+            offsets[2] = 0;
 
-            if (edgeFlags & YEDGE) {
-
-                quad[0] = v0;
-                coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1; // i, j, k-1
-                quad[1] = auxAcc.getValue(coord);
-                coord[0] -= 1; // i-1, j, k-1
-                quad[2] = auxAcc.getValue(coord);
-                coord[2] = ijk[2]; // i-1, j, k
-                quad[3] = auxAcc.getValue(coord);
+            const unsigned char cell = (SIGNS & flags);
 
-                mesher.addPrim(quad, isInside,
-                    (isSemLinePoly | isExteriorPoly[bool(refEdgeFlags & YEDGE)]));
+            if (sEdgeGroupTable[cell][0] > 1) {
+                offsets[0] = (sEdgeGroupTable[cell][1] - 1);
+                offsets[1] = (sEdgeGroupTable[cell][9] - 1);
+                offsets[2] = (sEdgeGroupTable[cell][4] - 1);
             }
 
-            if (edgeFlags & ZEDGE) {
-
-                quad[0] = v0;
-                coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]; // i, j-1, k
-                quad[1] = auxAcc.getValue(coord);
-                coord[0] -= 1; // i-1, j-1, k
-                quad[2] = auxAcc.getValue(coord);
-                coord[1] = ijk[1]; // i, j, k
-                quad[3] = auxAcc.getValue(coord);
-
-                mesher.addPrim(quad, !isInside,
-                     (isSemLinePoly | isExteriorPoly[bool(refEdgeFlags & ZEDGE)]));
+            if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
+                constructPolygons(flags, refFlags, offsets, ijk, *signleafPt, *idxLeafPt, mesher, mPointListSize);
+            } else {
+                constructPolygons(flags, refFlags, offsets, ijk, signAcc, idxAcc, mesher, mPointListSize);
             }
         }
 
@@ -1653,474 +2756,699 @@ MeshGen<DistTreeT, MeshingOp>::operator()(
 
 ////////////////////////////////////////
 
+// Masking and mesh partitioning
 
-template<class DistTreeT, class AuxDataT = int>
-class AuxiliaryData
+struct PartOp
 {
-public:
-    typedef openvdb::tree::ValueAccessor<const DistTreeT> SourceAccessorT;
-    typedef typename DistTreeT::ValueType ValueT;
 
-    typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
-    typedef openvdb::tree::ValueAccessor<CharTreeT> EdgeAccessorT;
+    PartOp(size_t leafCount, size_t partitions, size_t activePart)
+    {
+        size_t leafSegments = leafCount / partitions;
+        mStart = leafSegments * activePart;
+        mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
+    }
 
-    typedef typename DistTreeT::template ValueConverter<AuxDataT>::Type AuxTreeT;
-    typedef openvdb::tree::ValueAccessor<AuxTreeT> AuxAccessorT;
+    template <typename LeafNodeType>
+    void operator()(LeafNodeType &leaf, size_t leafIndex) const
+    {
+        if (leafIndex < mStart || leafIndex >= mEnd) leaf.setValuesOff();
+    }
 
-    AuxiliaryData(const DistTreeT&, const LeafCPtrList<DistTreeT>&,
-        double iso = 0.0, bool extraCheck = false);
-    AuxiliaryData(AuxiliaryData&, tbb::split);
+private:
+    size_t mStart, mEnd;
+};
 
-    void runParallel();
-    void runSerial();
 
-    typename CharTreeT::Ptr edgeTree() const { return mEdgeTree; }
-    typename AuxTreeT::Ptr auxTree() const { return mAuxTree; }
+////////////////////////////////////////
 
-    void operator()(const tbb::blocked_range<size_t>&);
 
-    void join(const AuxiliaryData& rhs)
-    {
-        mEdgeTree->merge(*rhs.mEdgeTree);
-        mAuxTree->merge(*rhs.mAuxTree);
-    }
+template<typename SrcTreeT>
+class PartGen
+{
+public:
+    typedef tree::LeafManager<const SrcTreeT> LeafManagerT;
+    typedef typename SrcTreeT::template ValueConverter<bool>::Type BoolTreeT;
+    typedef tree::ValueAccessor<BoolTreeT> BoolAccessorT;
 
-private:
-    const LeafCPtrList<DistTreeT>& mLeafNodes;
-    const DistTreeT& mSourceTree;
-    SourceAccessorT mSourceAccessor;
+    //////////
 
-    typename CharTreeT::Ptr mEdgeTree;
-    EdgeAccessorT mEdgeAccessor;
 
-    typename AuxTreeT::Ptr mAuxTree;
-    AuxAccessorT mAuxAccessor;
+    PartGen(const LeafManagerT& leafs, size_t partitions, size_t activePart);
 
-    const double mIsovalue;
-    const bool mExtraCheck;
+    void run(bool threaded = true);
+
+    BoolTreeT&  tree() { return mTree; }
 
-    int edgeCheck(const Coord& ijk, const bool thisInside);
+
+    //////////
+
+    PartGen(PartGen&, tbb::split);
+    void operator()(const tbb::blocked_range<size_t>&);
+    void join(PartGen& rhs) { mTree.merge(rhs.mTree); }
+
+private:
+    const LeafManagerT& mLeafManager;
+    BoolTreeT mTree;
+    size_t mStart, mEnd;
 };
 
-template<class DistTreeT, class AuxDataT>
-AuxiliaryData<DistTreeT, AuxDataT>::AuxiliaryData(const DistTreeT& tree,
-    const LeafCPtrList<DistTreeT>& leafNodes, double iso, bool extraCheck)
-    : mLeafNodes(leafNodes)
-    , mSourceTree(tree)
-    , mSourceAccessor(mSourceTree)
-    , mEdgeTree(new CharTreeT(0))
-    , mEdgeAccessor(*mEdgeTree)
-    , mAuxTree(new AuxTreeT(AuxDataT(0)))
-    , mAuxAccessor(*mAuxTree)
-    , mIsovalue(iso)
-    , mExtraCheck(extraCheck)
+template<typename SrcTreeT>
+PartGen<SrcTreeT>::PartGen(const LeafManagerT& leafs, size_t partitions, size_t activePart)
+    : mLeafManager(leafs)
+    , mTree(false)
+    , mStart(0)
+    , mEnd(0)
 {
+    size_t leafCount = leafs.leafCount();
+    size_t leafSegments = leafCount / partitions;
+    mStart = leafSegments * activePart;
+    mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
 }
 
-template<class DistTreeT, class AuxDataT>
-AuxiliaryData<DistTreeT, AuxDataT>::AuxiliaryData(AuxiliaryData& rhs, tbb::split)
-    : mLeafNodes(rhs.mLeafNodes)
-    , mSourceTree(rhs.mSourceTree)
-    , mSourceAccessor(mSourceTree)
-    , mEdgeTree(new CharTreeT(0))
-    , mEdgeAccessor(*mEdgeTree)
-    , mAuxTree(new AuxTreeT(AuxDataT(0)))
-    , mAuxAccessor(*mAuxTree)
-    , mIsovalue(rhs.mIsovalue)
-    , mExtraCheck(rhs.mExtraCheck)
+template<typename SrcTreeT>
+PartGen<SrcTreeT>::PartGen(PartGen& rhs, tbb::split)
+    : mLeafManager(rhs.mLeafManager)
+    , mTree(false)
+    , mStart(rhs.mStart)
+    , mEnd(rhs.mEnd)
 {
 }
 
 
+template<typename SrcTreeT>
+void
+PartGen<SrcTreeT>::run(bool threaded)
+{
+    if (threaded) tbb::parallel_reduce(mLeafManager.getRange(), *this);
+    else (*this)(mLeafManager.getRange());
+}
+
 
-template<class DistTreeT, typename AuxDataT>
+template<typename SrcTreeT>
 void
-AuxiliaryData<DistTreeT, AuxDataT>::runParallel()
+PartGen<SrcTreeT>::operator()(const tbb::blocked_range<size_t>& range)
 {
-    tbb::parallel_reduce(mLeafNodes.getRange(), *this);
+    Coord ijk;
+    BoolAccessorT acc(mTree);
+
+    typedef typename BoolTreeT::LeafNodeType BoolLeafT;
+    typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
+
+    for (size_t n = range.begin(); n != range.end(); ++n) {
+        if (n < mStart || n >= mEnd) continue;
+        BoolLeafT* leaf = acc.touchLeaf(mLeafManager.leaf(n).origin());
+        leaf->topologyUnion(mLeafManager.leaf(n));
+    }
 }
 
-template<class DistTreeT, typename AuxDataT>
+
+////////////////////////////////////////
+
+
+template<typename TreeT, typename LeafManagerT>
+class GenSeamMask
+{
+public:
+    typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
+
+    //////////
+
+    GenSeamMask(const LeafManagerT& leafs, const TreeT& tree);
+
+    void run(bool threaded = true);
+
+    BoolTreeT& mask() { return mMaskTree; }
+
+    //////////
+
+    GenSeamMask(GenSeamMask&, tbb::split);
+    void operator()(const tbb::blocked_range<size_t>&);
+    void join(GenSeamMask& rhs) { mMaskTree.merge(rhs.mMaskTree); }
+
+private:
+
+    const LeafManagerT& mLeafManager;
+    const TreeT& mTree;
+
+    BoolTreeT mMaskTree;
+};
+
+
+template<typename TreeT, typename LeafManagerT>
+GenSeamMask<TreeT, LeafManagerT>::GenSeamMask(const LeafManagerT& leafs, const TreeT& tree)
+    : mLeafManager(leafs)
+    , mTree(tree)
+    , mMaskTree(false)
+{
+}
+
+
+template<typename TreeT, typename LeafManagerT>
+GenSeamMask<TreeT, LeafManagerT>::GenSeamMask(GenSeamMask& rhs, tbb::split)
+    : mLeafManager(rhs.mLeafManager)
+    , mTree(rhs.mTree)
+    , mMaskTree(false)
+{
+}
+
+
+template<typename TreeT, typename LeafManagerT>
 void
-AuxiliaryData<DistTreeT, AuxDataT>::runSerial()
+GenSeamMask<TreeT, LeafManagerT>::run(bool threaded)
 {
-    (*this)(mLeafNodes.getRange());
+    if (threaded) tbb::parallel_reduce(mLeafManager.getRange(), *this);
+    else (*this)(mLeafManager.getRange());
 }
 
-template<class DistTreeT, typename AuxDataT>
+
+template<typename TreeT, typename LeafManagerT>
 void
-AuxiliaryData<DistTreeT, AuxDataT>::operator()(const tbb::blocked_range<size_t>& range)
+GenSeamMask<TreeT, LeafManagerT>::operator()(const tbb::blocked_range<size_t>& range)
 {
-    typename DistTreeT::LeafNodeType::ValueOnCIter iter;
     Coord ijk;
-    bool thisInside;
-    int edgeFlags;
-    ValueT val;
+    tree::ValueAccessor<const TreeT> acc(mTree);
+    tree::ValueAccessor<BoolTreeT> maskAcc(mMaskTree);
 
-    if (!mExtraCheck) {
-        for (size_t n = range.begin(); n != range.end(); ++n) {
-            for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
-                ijk = iter.getCoord();
-                thisInside = iter.getValue() < mIsovalue;
-                edgeFlags = edgeCheck(ijk, thisInside);
+    typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter it;
+
+    for (size_t n = range.begin(); n != range.end(); ++n) {
+
+        it = mLeafManager.leaf(n).cbeginValueOn();
+
+        for (; it; ++it) {
+
+            ijk = it.getCoord();
 
-                if (edgeFlags != 0) {
-                    edgeFlags |= int(thisInside);
-                    mEdgeAccessor.setValue(ijk, char(edgeFlags));
+            unsigned char rhsSigns = acc.getValue(ijk) & SIGNS;
+
+            if (sEdgeGroupTable[rhsSigns][0] > 0) {
+                unsigned char lhsSigns = it.getValue() & SIGNS;
+                if (rhsSigns != lhsSigns) {
+                    maskAcc.setValueOn(ijk);
                 }
             }
         }
-    } else {
-        for (size_t n = range.begin(); n != range.end(); ++n) {
-            for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
+    }
+}
 
-                ijk = iter.getCoord();
-                thisInside = iter.getValue() < mIsovalue;
-                edgeFlags = edgeCheck(ijk, thisInside);
 
-                if (edgeFlags != 0) {
-                    edgeFlags |= int(thisInside);
-                    mEdgeAccessor.setValue(ijk, char(edgeFlags));
-                }
+////////////////////////////////////////
 
-                --ijk[0];
-                if (!mSourceAccessor.probeValue(ijk, val)) {
-                    thisInside = val < mIsovalue;
-                    edgeFlags = edgeCheck(ijk, thisInside);
 
-                    if (edgeFlags != 0) {
-                        edgeFlags |= int(thisInside);
-                        mEdgeAccessor.setValue(ijk, char(edgeFlags));
-                    }
-                }
+template<typename TreeT>
+class TagSeamEdges
+{
+public:
+    typedef tree::ValueAccessor<const TreeT> AccessorT;
 
-                ++ijk[0];
-                --ijk[1];
-                if (!mSourceAccessor.probeValue(ijk, val)) {
-                    thisInside = val < mIsovalue;
-                    edgeFlags = edgeCheck(ijk, thisInside);
+    TagSeamEdges(const TreeT& tree) : mAcc(tree) {}
 
-                    if (edgeFlags != 0) {
-                        edgeFlags |= int(thisInside);
-                        mEdgeAccessor.setValue(ijk, char(edgeFlags));
-                    }
-                }
+    template <typename LeafNodeType>
+    void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
+    {
+        const typename TreeT::LeafNodeType *maskLeaf =
+            mAcc.probeConstLeaf(leaf.origin());
 
-                ++ijk[1];
-                --ijk[2];
-                if (!mSourceAccessor.probeValue(ijk, val)) {
-                    thisInside = val < mIsovalue;
-                    edgeFlags = edgeCheck(ijk, thisInside);
+        if (!maskLeaf) return;
 
-                    if (edgeFlags != 0) {
-                        edgeFlags |= int(thisInside);
-                        mEdgeAccessor.setValue(ijk, char(edgeFlags));
-                    }
-                }
+        typename LeafNodeType::ValueOnIter it = leaf.beginValueOn();
+
+        for (; it; ++it) {
+
+            if (maskLeaf->isValueOn(it.pos())) {
+                it.setValue(it.getValue() | SEAM);
             }
         }
     }
-}
 
-template<class DistTreeT, typename AuxDataT>
-inline int
-AuxiliaryData<DistTreeT, AuxDataT>::edgeCheck(const Coord& ijk, const bool thisInside)
+private:
+    AccessorT mAcc;
+};
+
+
+
+template<typename BoolTreeT>
+struct MaskEdges
 {
-    int edgeFlags = 0;
-    Coord n_ijk, coord;
+    typedef tree::ValueAccessor<const BoolTreeT> BoolAccessorT;
+
+    MaskEdges(const BoolTreeT& valueMask) : mMaskAcc(valueMask) {}
+
+    template <typename LeafNodeType>
+    void operator()(LeafNodeType &leaf, size_t /*leafIndex*/) const
+    {
+        typename LeafNodeType::ValueOnIter it = leaf.beginValueOn();
+
+        const typename BoolTreeT::LeafNodeType * maskLeaf =
+            mMaskAcc.probeConstLeaf(leaf.origin());
 
-    // Eval upwind x-edge
-    n_ijk = ijk; ++n_ijk[0];
-    bool otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
-    if (otherInside != thisInside) {
+        if (maskLeaf) {
+            for (; it; ++it) {
+                if (!maskLeaf->isValueOn(it.pos())) {
+                    it.setValue(0x1FF & it.getValue());
+                }
+            }
+        } else {
+            for (; it; ++it) {
+                it.setValue(0x1FF & it.getValue());
+            }
+        }
+    }
 
-        edgeFlags = XEDGE;
+private:
+    BoolAccessorT mMaskAcc;
+};
 
-        mAuxAccessor.setActiveState(ijk, true);
 
-        coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
-        mAuxAccessor.setActiveState(coord, true);
+class FlagUsedPoints
+{
+public:
+    //////////
 
-        coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]-1;
-        mAuxAccessor.setActiveState(coord, true);
+    FlagUsedPoints(const PolygonPoolList& polygons, size_t polyListCount,
+        std::vector<unsigned char>& usedPointMask)
+        : mPolygons(polygons)
+        , mPolyListCount(polyListCount)
+        , mUsedPointMask(usedPointMask)
+    {
+    }
 
-        coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
-        mAuxAccessor.setActiveState(coord, true);
+    void run(bool threaded = true)
+    {
+        if (threaded) {
+            tbb::parallel_for(tbb::blocked_range<size_t>(0, mPolyListCount), *this);
+        } else {
+            (*this)(tbb::blocked_range<size_t>(0, mPolyListCount));
+        }
     }
 
-    // Eval upwind y-edge
-    n_ijk[0] = ijk[0]; ++n_ijk[1];
-    otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
-    if (otherInside != thisInside) {
+    //////////
 
-        edgeFlags |= YEDGE;
+    void operator()(const tbb::blocked_range<size_t>& range) const
+    {
+        // Concurrent writes to same memory address can occur, but
+        // all threads are writing the same value and char is atomic.
+        for (size_t n = range.begin(); n != range.end(); ++n) {
+            const PolygonPool& polygons = mPolygons[n];
+            for (size_t i = 0; i < polygons.numQuads(); ++i) {
+                const Vec4I& quad = polygons.quad(i);
+                mUsedPointMask[quad[0]] = 1;
+                mUsedPointMask[quad[1]] = 1;
+                mUsedPointMask[quad[2]] = 1;
+                mUsedPointMask[quad[3]] = 1;
+            }
 
-        mAuxAccessor.setActiveState(ijk, true);
+            for (size_t i = 0; i < polygons.numTriangles(); ++i) {
+                const Vec3I& triangle = polygons.triangle(i);
+                mUsedPointMask[triangle[0]] = 1;
+                mUsedPointMask[triangle[1]] = 1;
+                mUsedPointMask[triangle[2]] = 1;
+            }
+        }
+    }
 
-        coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
-        mAuxAccessor.setActiveState(coord, true);
 
-        coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
-        mAuxAccessor.setActiveState(coord, true);
+private:
+    const PolygonPoolList& mPolygons;
+    size_t mPolyListCount;
+    std::vector<unsigned char>& mUsedPointMask;
+};
+
+class RemapIndices
+{
+public:
+    //////////
 
-        coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
-        mAuxAccessor.setActiveState(coord, true);
+    RemapIndices(PolygonPoolList& polygons,
+        size_t polyListCount, const std::vector<unsigned>& indexMap)
+        : mPolygons(polygons)
+        , mPolyListCount(polyListCount)
+        , mIndexMap(indexMap)
+    {
     }
 
-    // Eval upwind z-edge
-    n_ijk[1] = ijk[1]; ++n_ijk[2];
-    otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
-    if (otherInside != thisInside) {
+    void run(bool threaded = true)
+    {
+        if (threaded) {
+            tbb::parallel_for(tbb::blocked_range<size_t>(0, mPolyListCount), *this);
+        } else {
+            (*this)(tbb::blocked_range<size_t>(0, mPolyListCount));
+        }
+    }
 
-        edgeFlags |= ZEDGE;
+    //////////
 
-        mAuxAccessor.setActiveState(ijk, true);
+    void operator()(const tbb::blocked_range<size_t>& range) const
+    {
+        for (size_t n = range.begin(); n != range.end(); ++n) {
+            PolygonPool& polygons = mPolygons[n];
+            for (size_t i = 0; i < polygons.numQuads(); ++i) {
+                Vec4I& quad = polygons.quad(i);
+                quad[0] = mIndexMap[quad[0]];
+                quad[1] = mIndexMap[quad[1]];
+                quad[2] = mIndexMap[quad[2]];
+                quad[3] = mIndexMap[quad[3]];
+            }
 
-        coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
-        mAuxAccessor.setActiveState(coord, true);
+            for (size_t i = 0; i < polygons.numTriangles(); ++i) {
+                Vec3I& triangle = polygons.triangle(i);
+                triangle[0] = mIndexMap[triangle[0]];
+                triangle[1] = mIndexMap[triangle[1]];
+                triangle[2] = mIndexMap[triangle[2]];
+            }
+        }
+    }
+
+
+private:
+    PolygonPoolList& mPolygons;
+    size_t mPolyListCount;
+    const std::vector<unsigned>& mIndexMap;
+};
 
-        coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
-        mAuxAccessor.setActiveState(coord, true);
 
-        coord[0] = ijk[0]-1; coord[1] = ijk[1]-1; coord[2] = ijk[2];
-        mAuxAccessor.setActiveState(coord, true);
+class MovePoints
+{
+public:
+    //////////
+
+    MovePoints(
+        internal::UniquePtr<openvdb::Vec3s>::type& newPointList,
+        const PointList& oldPointList,
+        const std::vector<unsigned>& indexMap,
+        const std::vector<unsigned char>& usedPointMask)
+        : mNewPointList(newPointList)
+        , mOldPointList(oldPointList)
+        , mIndexMap(indexMap)
+        , mUsedPointMask(usedPointMask)
+    {
+    }
+
+    void run(bool threaded = true)
+    {
+        if (threaded) {
+            tbb::parallel_for(tbb::blocked_range<size_t>(0, mIndexMap.size()), *this);
+        } else {
+            (*this)(tbb::blocked_range<size_t>(0, mIndexMap.size()));
+        }
     }
-    return edgeFlags;
-}
+
+    //////////
+
+    void operator()(const tbb::blocked_range<size_t>& range) const
+    {
+        for (size_t n = range.begin(); n != range.end(); ++n) {
+            if (mUsedPointMask[n]) {
+                const size_t index = mIndexMap[n];
+                mNewPointList.get()[index] = mOldPointList[n];
+            }
+        }
+    }
+
+private:
+    internal::UniquePtr<openvdb::Vec3s>::type& mNewPointList;
+    const PointList& mOldPointList;
+    const std::vector<unsigned>& mIndexMap;
+    const std::vector<unsigned char>& mUsedPointMask;
+};
 
 
 ////////////////////////////////////////
 
 
-template <class DistTreeT>
-class SeamMaskGen
+template<typename SrcTreeT>
+class GenTopologyMask
 {
 public:
-    typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
-    typedef typename DistTreeT::template ValueConverter<int>::Type  IntTreeT;
-    typedef tree::ValueAccessor<const BoolTreeT> BoolTreeAccessorT;
-    typedef tree::ValueAccessor<const IntTreeT> IntTreeAccessorT;
+    typedef tree::LeafManager<const SrcTreeT> LeafManagerT;
+    typedef typename SrcTreeT::template ValueConverter<bool>::Type BoolTreeT;
+    typedef tree::ValueAccessor<const SrcTreeT> SrcAccessorT;
+    typedef tree::ValueAccessor<BoolTreeT> BoolAccessorT;
+    typedef Grid<BoolTreeT> BoolGridT;
 
-    SeamMaskGen(LeafPtrList<BoolTreeT>& seamMaskLeafs,
-        const BoolTreeT& topologyMaskTree, const IntTreeT& auxTree);
 
-    SeamMaskGen(const SeamMaskGen<DistTreeT>&);
+    //////////
 
-    void runParallel();
-    void runSerial();
 
-    void operator()(const tbb::blocked_range<size_t>&) const;
+    GenTopologyMask(const BoolGridT& mask, const LeafManagerT& srcLeafs,
+        const math::Transform& srcXForm, bool invertMask);
+
+    void run(bool threaded = true);
+
+    BoolTreeT& tree() { return mTree; }
+
+
+    //////////
+
+    GenTopologyMask(GenTopologyMask&, tbb::split);
+
+    void operator()(const tbb::blocked_range<size_t>&);
+
+    void join(GenTopologyMask& rhs) { mTree.merge(rhs.mTree); }
 
 private:
-    LeafPtrList<BoolTreeT>& mSeamMaskLeafs;
-    const BoolTreeT& mTopologyMaskTree;
-    BoolTreeAccessorT mTopologyMaskAcc;
-    const IntTreeT& mAuxTree;
-    IntTreeAccessorT mAuxAcc;
+
+    const BoolGridT& mMask;
+    const LeafManagerT& mLeafManager;
+    const math::Transform& mSrcXForm;
+    bool mInvertMask;
+    BoolTreeT mTree;
 };
 
 
-template <class DistTreeT>
-SeamMaskGen<DistTreeT>::SeamMaskGen(LeafPtrList<BoolTreeT>& seamMaskLeafs,
-    const BoolTreeT& topologyMaskTree, const IntTreeT& auxTree)
-    : mSeamMaskLeafs(seamMaskLeafs)
-    , mTopologyMaskTree(topologyMaskTree)
-    , mTopologyMaskAcc(mTopologyMaskTree)
-    , mAuxTree(auxTree)
-    , mAuxAcc(mAuxTree)
+template<typename SrcTreeT>
+GenTopologyMask<SrcTreeT>::GenTopologyMask(const BoolGridT& mask, const LeafManagerT& srcLeafs,
+    const math::Transform& srcXForm, bool invertMask)
+    : mMask(mask)
+    , mLeafManager(srcLeafs)
+    , mSrcXForm(srcXForm)
+    , mInvertMask(invertMask)
+    , mTree(false)
 {
 }
 
-template <class DistTreeT>
-SeamMaskGen<DistTreeT>::SeamMaskGen(const SeamMaskGen<DistTreeT>& rhs)
-    : mSeamMaskLeafs(rhs.mSeamMaskLeafs)
-    , mTopologyMaskTree(rhs.mTopologyMaskTree)
-    , mTopologyMaskAcc(mTopologyMaskTree)
-    , mAuxTree(rhs.mAuxTree)
-    , mAuxAcc(mAuxTree)
-{
-}
 
-template <class DistTreeT>
-void
-SeamMaskGen<DistTreeT>::runParallel()
+template<typename SrcTreeT>
+GenTopologyMask<SrcTreeT>::GenTopologyMask(GenTopologyMask& rhs, tbb::split)
+    : mMask(rhs.mMask)
+    , mLeafManager(rhs.mLeafManager)
+    , mSrcXForm(rhs.mSrcXForm)
+    , mInvertMask(rhs.mInvertMask)
+    , mTree(false)
 {
-    tbb::parallel_for(mSeamMaskLeafs.getRange(), *this);
 }
 
-template <class DistTreeT>
+
+template<typename SrcTreeT>
 void
-SeamMaskGen<DistTreeT>::runSerial()
+GenTopologyMask<SrcTreeT>::run(bool threaded)
 {
-    (*this)(mSeamMaskLeafs.getRange());
+    if (threaded) {
+        tbb::parallel_reduce(mLeafManager.getRange(), *this);
+    } else {
+        (*this)(mLeafManager.getRange());
+    }
 }
 
-template <class DistTreeT>
+
+template<typename SrcTreeT>
 void
-SeamMaskGen<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
+GenTopologyMask<SrcTreeT>::operator()(const tbb::blocked_range<size_t>& range)
 {
-    typedef typename BoolTreeT::LeafNodeType::ValueOnIter ValueOnIterT;
-    Coord ijk, n_ijk;
-    for (size_t leafIdx = range.begin(); leafIdx != range.end(); ++leafIdx) {
-        ValueOnIterT it = mSeamMaskLeafs[leafIdx]->beginValueOn();
-        for (; it; ++it) {
-            ijk = it.getCoord();
-            if (!mTopologyMaskAcc.isValueOn(ijk)) {
-                it.setValueOff();
-            } else {
-                bool turnOff = true;
-                for (size_t n = 0; n < 6; ++n) {
-                    n_ijk = ijk + util::COORD_OFFSETS[n];
-                    if (!mAuxTree.isValueOn(n_ijk) && mTopologyMaskAcc.isValueOn(n_ijk)) {
-                        turnOff = false;
-                        break;
+    Coord ijk;
+    Vec3d xyz;
+    typedef typename BoolTreeT::LeafNodeType BoolLeafT;
+    const math::Transform& maskXForm = mMask.transform();
+    tree::ValueAccessor<const BoolTreeT> maskAcc(mMask.tree());
+    tree::ValueAccessor<BoolTreeT> acc(mTree);
+
+    typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
+    for (size_t n = range.begin(); n != range.end(); ++n) {
+
+        ijk = mLeafManager.leaf(n).origin();
+        BoolLeafT* leaf = new BoolLeafT(ijk, false);
+        bool addLeaf = false;
+
+        if (maskXForm == mSrcXForm) {
+
+            const BoolLeafT* maskLeaf = maskAcc.probeConstLeaf(ijk);
+
+            if (maskLeaf) {
+
+                for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
+                    Index pos = iter.pos();
+                    if(maskLeaf->isValueOn(pos) != mInvertMask) {
+                        leaf->setValueOn(pos);
+                        addLeaf = true;
                     }
                 }
-                if (turnOff) it.setValueOff();
+
+            } else if (maskAcc.isValueOn(ijk) != mInvertMask) {
+                leaf->topologyUnion(mLeafManager.leaf(n));
+                addLeaf = true;
+            }
+
+        } else {
+            for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
+                ijk = iter.getCoord();
+                xyz = maskXForm.worldToIndex(mSrcXForm.indexToWorld(ijk));
+                if(maskAcc.isValueOn(util::nearestCoord(xyz)) != mInvertMask) {
+                    leaf->setValueOn(iter.pos());
+                    addLeaf = true;
+                }
             }
         }
+
+        if (addLeaf) acc.addLeaf(leaf);
+        else delete leaf;
     }
 }
 
+
 ////////////////////////////////////////
 
 
-template <class DistTreeT>
-class EdgeSmooth
+template<typename SrcTreeT>
+class GenBoundaryMask
 {
 public:
-    typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
-    typedef typename DistTreeT::template ValueConverter<int>::Type  IntTreeT;
-    typedef tree::ValueAccessor<const BoolTreeT> BoolTreeAccessorT;
-    typedef tree::ValueAccessor<const IntTreeT> IntTreeAccessorT;
+    typedef typename SrcTreeT::template ValueConverter<int>::Type IntTreeT;
+    typedef typename SrcTreeT::template ValueConverter<bool>::Type BoolTreeT;
+    typedef tree::LeafManager<const SrcTreeT> LeafManagerT;
 
-    EdgeSmooth(
-        LeafPtrList<BoolTreeT>& leafs,
-        const BoolTreeT& edgeMaskTree,
-        const IntTreeT& auxTree,
-        PointList& points,
-        const math::Transform& xform);
+    //////////
 
-    EdgeSmooth(const EdgeSmooth<DistTreeT>&);
+    GenBoundaryMask(const LeafManagerT& leafs, const BoolTreeT&, const IntTreeT&);
 
-    void runParallel(const size_t iterations);
-    void runSerial(const size_t iterations);
+    void run(bool threaded = true);
 
-    void operator()(const tbb::blocked_range<size_t>&) const;
+    BoolTreeT&  tree() { return mTree; }
+
+    //////////
+
+    GenBoundaryMask(GenBoundaryMask&, tbb::split);
+    void operator()(const tbb::blocked_range<size_t>&);
+    void join(GenBoundaryMask& rhs) { mTree.merge(rhs.mTree); }
 
 private:
-    LeafPtrList<BoolTreeT>& mLeafs;
-    const BoolTreeT& mEdgeMaskTree;
-    const IntTreeT& mAuxTree;
-    PointList& mPoints;
-    const math::Transform& mTransform;
+    // This typedef is needed for Windows
+    typedef tree::ValueAccessor<const IntTreeT> IntTreeAccessorT;
 
-    bool pointInAABB(const Vec3s& p, const Vec3s& bmin, const Vec3s& bmax) const
-    {
-        for (int i = 0; i < 3; ++i) {
-            if (p[i] < bmin[i] || p[i] > bmax[i]) {
-                return false;
-            }
-        }
-        return true;
-    }
+    bool neighboringLeaf(const Coord&, const IntTreeAccessorT&) const;
 
+    const LeafManagerT& mLeafManager;
+    const BoolTreeT& mMaskTree;
+    const IntTreeT& mIdxTree;
+    BoolTreeT mTree;
+    CoordBBox mLeafBBox;
 };
 
 
-template <class DistTreeT>
-EdgeSmooth<DistTreeT>::EdgeSmooth(
-    LeafPtrList<BoolTreeT>& leafs,
-    const BoolTreeT& edgeMaskTree,
-    const IntTreeT& auxTree,
-    PointList& points,
-    const math::Transform& xform)
-    : mLeafs(leafs)
-    , mEdgeMaskTree(edgeMaskTree)
-    , mAuxTree(auxTree)
-    , mPoints(points)
-    , mTransform(xform)
+template<typename SrcTreeT>
+GenBoundaryMask<SrcTreeT>::GenBoundaryMask(const LeafManagerT& leafs,
+    const BoolTreeT& maskTree, const IntTreeT& auxTree)
+    : mLeafManager(leafs)
+    , mMaskTree(maskTree)
+    , mIdxTree(auxTree)
+    , mTree(false)
 {
+    mIdxTree.evalLeafBoundingBox(mLeafBBox);
+    mLeafBBox.expand(IntTreeT::LeafNodeType::DIM);
 }
 
-template <class DistTreeT>
-EdgeSmooth<DistTreeT>::EdgeSmooth(const EdgeSmooth<DistTreeT>& rhs)
-    : mLeafs(rhs.mLeafs)
-    , mEdgeMaskTree(rhs.mEdgeMaskTree)
-    , mAuxTree(rhs.mAuxTree)
-    , mPoints(rhs.mPoints)
-    , mTransform(rhs.mTransform)
+
+template<typename SrcTreeT>
+GenBoundaryMask<SrcTreeT>::GenBoundaryMask(GenBoundaryMask& rhs, tbb::split)
+    : mLeafManager(rhs.mLeafManager)
+    , mMaskTree(rhs.mMaskTree)
+    , mIdxTree(rhs.mIdxTree)
+    , mTree(false)
+    , mLeafBBox(rhs.mLeafBBox)
 {
 }
 
-template <class DistTreeT>
+
+template<typename SrcTreeT>
 void
-EdgeSmooth<DistTreeT>::runParallel(const size_t iterations)
+GenBoundaryMask<SrcTreeT>::run(bool threaded)
 {
-    for (size_t i = 0; i < iterations; ++i) {
-        tbb::parallel_for(mLeafs.getRange(), *this);
+    if (threaded) {
+        tbb::parallel_reduce(mLeafManager.getRange(), *this);
+    } else {
+        (*this)(mLeafManager.getRange());
     }
 }
 
-template <class DistTreeT>
-void
-EdgeSmooth<DistTreeT>::runSerial(const size_t iterations)
+
+template<typename SrcTreeT>
+bool
+GenBoundaryMask<SrcTreeT>::neighboringLeaf(const Coord& ijk, const IntTreeAccessorT& acc) const
 {
-    for (size_t i = 0; i < iterations; ++i) {
-        (*this)(mLeafs.getRange());
-    }
+    if (acc.probeConstLeaf(ijk)) return true;
+
+    const int dim = IntTreeT::LeafNodeType::DIM;
+
+    // face adjacent neghbours
+    if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1], ijk[2]))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1], ijk[2]))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] + dim, ijk[2]))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] - dim, ijk[2]))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0], ijk[1], ijk[2] + dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0], ijk[1], ijk[2] - dim))) return true;
+
+    // edge adjacent neighbors
+    if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1], ijk[2] - dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1], ijk[2] - dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1], ijk[2] + dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1], ijk[2] + dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] + dim, ijk[2]))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] + dim, ijk[2]))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] - dim, ijk[2]))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] - dim, ijk[2]))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] - dim, ijk[2] + dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] - dim, ijk[2] - dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] + dim, ijk[2] + dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0], ijk[1] + dim, ijk[2] - dim))) return true;
+
+    // corner adjacent neighbors
+    if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] - dim, ijk[2] - dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] - dim, ijk[2] + dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] - dim, ijk[2] + dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] - dim, ijk[2] - dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] + dim, ijk[2] - dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] - dim, ijk[1] + dim, ijk[2] + dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] + dim, ijk[2] + dim))) return true;
+    if (acc.probeConstLeaf(Coord(ijk[0] + dim, ijk[1] + dim, ijk[2] - dim))) return true;
+
+    return false;
 }
 
-template <class DistTreeT>
+
+template<typename SrcTreeT>
 void
-EdgeSmooth<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
+GenBoundaryMask<SrcTreeT>::operator()(const tbb::blocked_range<size_t>& range)
 {
-    typedef typename BoolTreeT::LeafNodeType::ValueOnIter ValueOnIterT;
-    BoolTreeAccessorT maskAcc(mEdgeMaskTree);
-    IntTreeAccessorT auxAcc(mAuxTree);
-
-    float dx = float(mTransform.voxelSize()[0]);
-    openvdb::Vec3s avg, bmin, bmax;
-    Coord ijk, nijk;
-    int count;
-
-    for (size_t leafIdx = range.begin(); leafIdx != range.end(); ++leafIdx) {
-        typename BoolTreeT::LeafNodeType::ValueOnIter valueIt = mLeafs[leafIdx]->beginValueOn();
-        for ( ; valueIt; ++valueIt) {
-
-            ijk = valueIt.getCoord();
-            openvdb::Vec3s& ptn = mPoints[auxAcc.getValue(ijk)];
-
-            avg = ptn;
-            count = 1;
-            for (int n = 0; n < 26; ++n) {
-                nijk = ijk + util::COORD_OFFSETS[n];
-                if (maskAcc.isValueOn(nijk)) {
-                    avg += mPoints[auxAcc.getValue(nijk)];
-                    ++count;
-                }
-            }
+    Coord ijk;
+    tree::ValueAccessor<const BoolTreeT> maskAcc(mMaskTree);
+    tree::ValueAccessor<const IntTreeT> idxAcc(mIdxTree);
+    tree::ValueAccessor<BoolTreeT> acc(mTree);
 
-            if (count > 1) {
-                avg *= (1.0 / float(count));
+    typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
 
-                // Constrain to current cell
-                bmin = openvdb::Vec3s(mTransform.indexToWorld(ijk));
-                bmax = bmin + dx;
+    for (size_t n = range.begin(); n != range.end(); ++n) {
 
-                bool inCell = true;
-                for (int i = 0; i < 10; ++i) {
+        const typename SrcTreeT::LeafNodeType&
+            leaf = mLeafManager.leaf(n);
 
-                    inCell = pointInAABB(avg, bmin, bmax);
+        ijk = leaf.origin();
 
-                    if (inCell) break;
+        if (!mLeafBBox.isInside(ijk) || !neighboringLeaf(ijk, idxAcc)) continue;
 
-                    avg += ptn;
-                    avg *= 0.5;
-                }
+        const typename BoolTreeT::LeafNodeType*
+            maskLeaf = maskAcc.probeConstLeaf(ijk);
 
-                if (inCell) ptn = avg;
-            }
+        if (!maskLeaf || !leaf.hasSameTopology(maskLeaf)) {
+            acc.touchLeaf(ijk)->topologyUnion(leaf);
         }
     }
 }
@@ -2129,235 +3457,239 @@ EdgeSmooth<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
 ////////////////////////////////////////
 
 
-template<class CharAccessorT, typename AuxAccessorT>
-class AuxDataGenerator
+template<typename TreeT>
+class GenTileMask
 {
 public:
-    AuxDataGenerator(CharAccessorT& edgeAcc, AuxAccessorT& auxAcc)
-        : mEdgeAcc(edgeAcc), mAuxAcc(auxAcc) {}
+    typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
 
+    typedef typename TreeT::ValueType ValueT;
 
-    void setXEdge(char edgeFlags, const Coord& ijk)
-    {
-        mEdgeAcc.setValue(ijk, edgeFlags | XEDGE);
+    //////////
 
-        mAuxAcc.setActiveState(ijk, true);
+    GenTileMask(const std::vector<Vec4i>& tiles, const TreeT& distTree, ValueT iso);
 
-        Coord coord = ijk;
-        coord[1] = ijk[1]-1;
-        mAuxAcc.setActiveState(coord, true);
+    void run(bool threaded = true);
 
-        coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]-1;
-        mAuxAcc.setActiveState(coord, true);
+    BoolTreeT& tree() { return mTree; }
 
-        coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
-        mAuxAcc.setActiveState(coord, true);
-    }
+    //////////
 
-    void setYEdge(char edgeFlags, const Coord& ijk)
-    {
-        mEdgeAcc.setValue(ijk, edgeFlags | YEDGE);
-
-        mAuxAcc.setActiveState(ijk, true);
-
-        Coord coord = ijk;
-        coord[2] = ijk[2]-1;
-        mAuxAcc.setActiveState(coord, true);
-
-        coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
-        mAuxAcc.setActiveState(coord, true);
+    GenTileMask(GenTileMask&, tbb::split);
+    void operator()(const tbb::blocked_range<size_t>&);
+    void join(GenTileMask& rhs) { mTree.merge(rhs.mTree); }
 
-        coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
-        mAuxAcc.setActiveState(coord, true);
-    }
+private:
 
-    void setZEdge(char edgeFlags, const Coord& ijk)
-    {
-        mEdgeAcc.setValue(ijk, edgeFlags | ZEDGE);
+    const std::vector<Vec4i>& mTiles;
+    const TreeT& mDistTree;
+    ValueT mIsovalue;
 
-        mAuxAcc.setActiveState(ijk, true);
+    BoolTreeT mTree;
+};
 
-        Coord coord = ijk;
-        coord[1] = ijk[1]-1;
-        mAuxAcc.setActiveState(coord, true);
 
-        coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
-        mAuxAcc.setActiveState(coord, true);
+template<typename TreeT>
+GenTileMask<TreeT>::GenTileMask(
+    const std::vector<Vec4i>& tiles, const TreeT& distTree, ValueT iso)
+    : mTiles(tiles)
+    , mDistTree(distTree)
+    , mIsovalue(iso)
+    , mTree(false)
+{
+}
 
-        coord[0] = ijk[0]-1; coord[1] = ijk[1]-1; coord[2] = ijk[2];
-        mAuxAcc.setActiveState(coord, true);
-    }
 
-private:
-    CharAccessorT& mEdgeAcc;
-    AuxAccessorT& mAuxAcc;
-};
+template<typename TreeT>
+GenTileMask<TreeT>::GenTileMask(GenTileMask& rhs, tbb::split)
+    : mTiles(rhs.mTiles)
+    , mDistTree(rhs.mDistTree)
+    , mIsovalue(rhs.mIsovalue)
+    , mTree(false)
+{
+}
 
 
-////////////////////////////////////////
+template<typename TreeT>
+void
+GenTileMask<TreeT>::run(bool threaded)
+{
+    if (threaded) tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mTiles.size()), *this);
+    else (*this)(tbb::blocked_range<size_t>(0, mTiles.size()));
+}
 
 
-template<class DistTreeT, class AuxTreeT, class CharTreeT>
-inline void
-tileAuxiliaryData(
-    const DistTreeT& distTree, CharTreeT& edgeTree, AuxTreeT& auxTree,
-    double iso)
+template<typename TreeT>
+void
+GenTileMask<TreeT>::operator()(const tbb::blocked_range<size_t>& range)
 {
-    typedef tree::ValueAccessor<const DistTreeT> DistAccessorT;
-    typedef tree::ValueAccessor<AuxTreeT> AuxTreeAccessorT;
-    typedef tree::ValueAccessor<CharTreeT> CharTreeAccessorT;
-
-    typename DistTreeT::ValueType isoValue = typename DistTreeT::ValueType(iso);
+    tree::ValueAccessor<const TreeT> distAcc(mDistTree);
+    CoordBBox region, bbox;
+    Coord ijk, nijk;
+    bool processRegion = true;
+    ValueT value;
 
-    DistAccessorT distAcc(distTree);
-    CharTreeAccessorT edgeAcc(edgeTree);
-    AuxTreeAccessorT auxAcc(auxTree);
 
-    AuxDataGenerator<CharTreeAccessorT, AuxTreeAccessorT> auxData(edgeAcc, auxAcc);
+    for (size_t n = range.begin(); n != range.end(); ++n) {
 
-    Coord ijk, nijk;
-    typename DistTreeT::ValueType value;
-    CoordBBox bbox;
-    bool processTileFace;
+        const Vec4i& tile = mTiles[n];
 
-    typename DistTreeT::ValueOnCIter tileIter(distTree);
-    tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
+        bbox.min()[0] = tile[0];
+        bbox.min()[1] = tile[1];
+        bbox.min()[2] = tile[2];
 
-    for ( ; tileIter; ++tileIter) {
-        tileIter.getBoundingBox(bbox);
+        bbox.max() = bbox.min();
+        bbox.max().offset(tile[3]);
 
-        const bool thisInside = (distAcc.getValue(bbox.min()) < isoValue);
+        const bool thisInside = (distAcc.getValue(bbox.min()) < mIsovalue);
         const int thisDepth = distAcc.getValueDepth(bbox.min());
 
-
-        // Eval x-edges
+        // eval x-edges
 
         ijk = bbox.max();
         nijk = ijk;
         ++nijk[0];
 
-        processTileFace = true;
+        processRegion = true;
         if (thisDepth >= distAcc.getValueDepth(nijk)) {
-            processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
+            processRegion = thisInside != (distAcc.getValue(nijk) < mIsovalue);
         }
 
-        if (processTileFace) {
-            for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
-                nijk[1] = ijk[1];
-                for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
-                    nijk[2] = ijk[2];
-                    if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
-                        auxData.setXEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
-                    }
-                }
-            }
+
+        if (processRegion) {
+            region = bbox;
+            region.min()[0] = region.max()[0] = ijk[0];
+            mTree.fill(region, true);
         }
 
+
         ijk = bbox.min();
         --ijk[0];
 
-        processTileFace = true;
+        processRegion = true;
         if (thisDepth >= distAcc.getValueDepth(ijk)) {
-            processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
+            processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
         }
 
-        if (processTileFace) {
-            for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
-                for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
-                    if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
-                         auxData.setXEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
-                    }
-                }
-            }
+        if (processRegion) {
+            region = bbox;
+            region.min()[0] = region.max()[0] = ijk[0];
+            mTree.fill(region, true);
         }
 
 
-        // Eval y-edges
+        // eval y-edges
 
         ijk = bbox.max();
         nijk = ijk;
         ++nijk[1];
 
-        processTileFace = true;
+        processRegion = true;
         if (thisDepth >= distAcc.getValueDepth(nijk)) {
-            processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
+            processRegion = thisInside != (distAcc.getValue(nijk) < mIsovalue);
         }
 
-        if (processTileFace) {
-            for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
-                nijk[0] = ijk[0];
-                for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
-                    nijk[2] = ijk[2];
-
-                    if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
-                        auxData.setYEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
-                    }
-                }
-            }
+        if (processRegion) {
+            region = bbox;
+            region.min()[1] = region.max()[1] = ijk[1];
+            mTree.fill(region, true);
         }
 
 
         ijk = bbox.min();
         --ijk[1];
 
-        processTileFace = true;
+        processRegion = true;
         if (thisDepth >= distAcc.getValueDepth(ijk)) {
-            processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
+            processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
         }
 
-        if (processTileFace) {
-            for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
-                for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
-
-                    if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
-                        auxData.setYEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
-                    }
-                }
-            }
+        if (processRegion) {
+            region = bbox;
+            region.min()[1] = region.max()[1] = ijk[1];
+            mTree.fill(region, true);
         }
 
 
-        // Eval z-edges
+        // eval z-edges
 
         ijk = bbox.max();
         nijk = ijk;
         ++nijk[2];
 
-        processTileFace = true;
+        processRegion = true;
         if (thisDepth >= distAcc.getValueDepth(nijk)) {
-            processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
+            processRegion = thisInside != (distAcc.getValue(nijk) < mIsovalue);
         }
 
-        if (processTileFace) {
-            for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
-                nijk[0] = ijk[0];
-                for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
-                    nijk[1] = ijk[1];
+        if (processRegion) {
+            region = bbox;
+            region.min()[2] = region.max()[2] = ijk[2];
+            mTree.fill(region, true);
+        }
 
-                    if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
-                        auxData.setZEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
-                    }
-                }
-            }
+        ijk = bbox.min();
+        --ijk[2];
+
+        processRegion = true;
+        if (thisDepth >= distAcc.getValueDepth(ijk)) {
+            processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
         }
 
+        if (processRegion) {
+            region = bbox;
+            region.min()[2] = region.max()[2] = ijk[2];
+            mTree.fill(region, true);
+        }
+
+
         ijk = bbox.min();
+        --ijk[1];
         --ijk[2];
 
-        processTileFace = true;
+        processRegion = true;
         if (thisDepth >= distAcc.getValueDepth(ijk)) {
-            processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
+            processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
         }
 
-        if (processTileFace) {
-            for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
-                for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
+        if (processRegion) {
+            region = bbox;
+            region.min()[1] = region.max()[1] = ijk[1];
+            region.min()[2] = region.max()[2] = ijk[2];
+            mTree.fill(region, true);
+        }
 
-                    if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
-                        auxData.setZEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
-                    }
-                }
-            }
+
+        ijk = bbox.min();
+        --ijk[0];
+        --ijk[1];
+
+        processRegion = true;
+        if (thisDepth >= distAcc.getValueDepth(ijk)) {
+            processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
+        }
+
+        if (processRegion) {
+            region = bbox;
+            region.min()[1] = region.max()[1] = ijk[1];
+            region.min()[0] = region.max()[0] = ijk[0];
+            mTree.fill(region, true);
+        }
+
+        ijk = bbox.min();
+        --ijk[0];
+        --ijk[2];
+
+        processRegion = true;
+        if (thisDepth >= distAcc.getValueDepth(ijk)) {
+            processRegion = !distAcc.probeValue(ijk, value) && thisInside != (value < mIsovalue);
+        }
+
+        if (processRegion) {
+            region = bbox;
+            region.min()[2] = region.max()[2] = ijk[2];
+            region.min()[0] = region.max()[0] = ijk[0];
+            mTree.fill(region, true);
         }
     }
 }
@@ -2366,6 +3698,58 @@ tileAuxiliaryData(
 ////////////////////////////////////////
 
 
+template<class DistTreeT, class SignTreeT, class IdxTreeT>
+inline void
+tileData(const DistTreeT& distTree, SignTreeT& signTree, IdxTreeT& idxTree, double iso)
+{
+    typename DistTreeT::ValueOnCIter tileIter(distTree);
+    tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
+
+    if (!tileIter) return; // volume has no active tiles.
+
+    size_t tileCount = 0;
+    for ( ; tileIter; ++tileIter) {
+        ++tileCount;
+    }
+
+    std::vector<Vec4i> tiles(tileCount);
+
+    tileCount = 0;
+    tileIter = distTree.cbeginValueOn();
+    tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
+
+    CoordBBox bbox;
+    for (; tileIter; ++tileIter) {
+        Vec4i& tile = tiles[tileCount++];
+        tileIter.getBoundingBox(bbox);
+        tile[0] = bbox.min()[0];
+        tile[1] = bbox.min()[1];
+        tile[2] = bbox.min()[2];
+        tile[3] = bbox.max()[0] - bbox.min()[0];
+    }
+
+    typename DistTreeT::ValueType isovalue = typename DistTreeT::ValueType(iso);
+
+    GenTileMask<DistTreeT> tileMask(tiles, distTree, isovalue);
+    tileMask.run();
+
+    typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
+    typedef tree::LeafManager<BoolTreeT> BoolLeafManagerT;
+
+    BoolLeafManagerT leafs(tileMask.tree());
+
+
+    internal::SignData<DistTreeT, BoolLeafManagerT> op(distTree, leafs, isovalue);
+    op.run();
+
+    signTree.merge(*op.signTree());
+    idxTree.merge(*op.idxTree());
+}
+
+
+////////////////////////////////////////
+
+
 // Utility class for the volumeToMesh wrapper
 class PointListCopy
 {
@@ -2388,30 +3772,216 @@ private:
 };
 
 
-} // namespace internal
+// Checks if the isovalue is in proximity to the active voxel boundary.
+template <typename LeafManagerT>
+inline bool
+needsActiveVoxePadding(const LeafManagerT& leafs, double iso, double voxelSize)
+{
+    double interiorWidth = 0.0, exteriorWidth = 0.0;
+    {
+        typename LeafManagerT::TreeType::LeafNodeType::ValueOffCIter it;
+        bool foundInterior = false, foundExterior = false;
+        for (size_t n = 0, N = leafs.leafCount(); n < N; ++n) {
+
+            for (it = leafs.leaf(n).cbeginValueOff(); it; ++it) {
+                double value = double(it.getValue());
+                if (value < 0.0) {
+                    interiorWidth = value;
+                    foundInterior = true;
+                } else if (value > 0.0) {
+                    exteriorWidth = value;
+                    foundExterior = true;
+                }
+
+                if (foundInterior && foundExterior) break;
+            }
+
+            if (foundInterior && foundExterior) break;
+        }
+
+    }
+
+    double minDist = std::min(std::abs(interiorWidth - iso), std::abs(exteriorWidth - iso));
+    return !(minDist > (2.0 * voxelSize));
+}
+
+
+} // end namespace internal
+
+
+////////////////////////////////////////
+
+
+inline
+PolygonPool::PolygonPool()
+    : mNumQuads(0)
+    , mNumTriangles(0)
+    , mQuads(NULL)
+    , mTriangles(NULL)
+    , mQuadFlags(NULL)
+    , mTriangleFlags(NULL)
+{
+}
+
+
+inline
+PolygonPool::PolygonPool(const size_t numQuads, const size_t numTriangles)
+    : mNumQuads(numQuads)
+    , mNumTriangles(numTriangles)
+    , mQuads(new openvdb::Vec4I[mNumQuads])
+    , mTriangles(new openvdb::Vec3I[mNumTriangles])
+    , mQuadFlags(new char[mNumQuads])
+    , mTriangleFlags(new char[mNumTriangles])
+{
+}
+
+
+inline void
+PolygonPool::copy(const PolygonPool& rhs)
+{
+    resetQuads(rhs.numQuads());
+    resetTriangles(rhs.numTriangles());
+
+    for (size_t i = 0; i < mNumQuads; ++i) {
+        mQuads[i] = rhs.mQuads[i];
+        mQuadFlags[i] = rhs.mQuadFlags[i];
+    }
+
+    for (size_t i = 0; i < mNumTriangles; ++i) {
+        mTriangles[i] = rhs.mTriangles[i];
+        mTriangleFlags[i] = rhs.mTriangleFlags[i];
+    }
+}
+
+
+inline void
+PolygonPool::resetQuads(size_t size)
+{
+    mNumQuads = size;
+    mQuads.reset(new openvdb::Vec4I[mNumQuads]);
+    mQuadFlags.reset(new char[mNumQuads]);
+}
+
+
+inline void
+PolygonPool::clearQuads()
+{
+    mNumQuads = 0;
+    mQuads.reset(NULL);
+    mQuadFlags.reset(NULL);
+}
+
+
+inline void
+PolygonPool::resetTriangles(size_t size)
+{
+    mNumTriangles = size;
+    mTriangles.reset(new openvdb::Vec3I[mNumTriangles]);
+    mTriangleFlags.reset(new char[mNumTriangles]);
+}
+
+
+inline void
+PolygonPool::clearTriangles()
+{
+    mNumTriangles = 0;
+    mTriangles.reset(NULL);
+    mTriangleFlags.reset(NULL);
+}
+
+
+inline bool
+PolygonPool::trimQuads(const size_t n, bool reallocate)
+{
+    if (!(n < mNumQuads)) return false;
+
+    if (reallocate) {
+
+        if (n == 0) {
+            mQuads.reset(NULL);
+        } else {
+
+            boost::scoped_array<openvdb::Vec4I> quads(new openvdb::Vec4I[n]);
+            boost::scoped_array<char> flags(new char[n]);
+
+            for (size_t i = 0; i < n; ++i) {
+                quads[i] = mQuads[i];
+                flags[i] = mQuadFlags[i];
+            }
+
+            mQuads.swap(quads);
+            mQuadFlags.swap(flags);
+        }
+    }
+
+    mNumQuads = n;
+    return true;
+}
+
+
+inline bool
+PolygonPool::trimTrinagles(const size_t n, bool reallocate)
+{
+    if (!(n < mNumTriangles)) return false;
+
+    if (reallocate) {
+
+        if (n == 0) {
+            mTriangles.reset(NULL);
+        } else {
+
+            boost::scoped_array<openvdb::Vec3I> triangles(new openvdb::Vec3I[n]);
+            boost::scoped_array<char> flags(new char[n]);
+
+            for (size_t i = 0; i < n; ++i) {
+                triangles[i] = mTriangles[i];
+                flags[i] = mTriangleFlags[i];
+            }
+
+            mTriangles.swap(triangles);
+            mTriangleFlags.swap(flags);
+        }
+    }
+
+    mNumTriangles = n;
+    return true;
+}
 
 
 ////////////////////////////////////////
 
 
 inline VolumeToMesh::VolumeToMesh(double isovalue, double adaptivity)
-    : mPointListSize(0)
+    : mPoints(NULL)
+    , mPolygons()
+    , mPointListSize(0)
+    , mSeamPointListSize(0)
     , mPolygonPoolListSize(0)
     , mIsovalue(isovalue)
     , mPrimAdaptivity(adaptivity)
     , mSecAdaptivity(0.0)
     , mRefGrid(GridBase::ConstPtr())
-    , mRefEdgeTree(TreeBase::Ptr())
-    , mRefTopologyMaskTree(TreeBase::Ptr())
+    , mSurfaceMaskGrid(GridBase::ConstPtr())
+    , mAdaptivityGrid(GridBase::ConstPtr())
+    , mAdaptivityMaskTree(TreeBase::ConstPtr())
+    , mRefSignTree(TreeBase::Ptr())
+    , mRefIdxTree(TreeBase::Ptr())
+    , mInvertSurfaceMask(false)
+    , mPartitions(1)
+    , mActivePart(0)
+    , mQuantizedSeamPoints(NULL)
+    , mPointFlags(0)
 {
 }
 
+
 inline PointList&
 VolumeToMesh::pointList()
 {
     return mPoints;
 }
 
+
 inline const size_t&
 VolumeToMesh::pointListSize() const
 {
@@ -2441,15 +4011,60 @@ VolumeToMesh::polygonPoolListSize() const
 
 
 inline void
-VolumeToMesh::setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity, bool smoothSeams)
+VolumeToMesh::setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity)
 {
     mRefGrid = grid;
     mSecAdaptivity = secAdaptivity;
+
     // Clear out old auxiliary data
-    mRefEdgeTree = TreeBase::Ptr();
-    mRefTopologyMaskTree = TreeBase::Ptr();
-    mSeamPointTree = TreeBase::Ptr();
-    mSmoothSeams = smoothSeams;
+    mRefSignTree = TreeBase::Ptr();
+    mRefIdxTree = TreeBase::Ptr();
+    mSeamPointListSize = 0;
+    mQuantizedSeamPoints.reset(NULL);
+}
+
+
+inline void
+VolumeToMesh::setSurfaceMask(const GridBase::ConstPtr& mask, bool invertMask)
+{
+    mSurfaceMaskGrid = mask;
+    mInvertSurfaceMask = invertMask;
+}
+
+
+inline void
+VolumeToMesh::setSpatialAdaptivity(const GridBase::ConstPtr& grid)
+{
+    mAdaptivityGrid = grid;
+}
+
+
+inline void
+VolumeToMesh::setAdaptivityMask(const TreeBase::ConstPtr& tree)
+{
+   mAdaptivityMaskTree = tree;
+}
+
+
+inline void
+VolumeToMesh::partition(unsigned partitions, unsigned activePart)
+{
+    mPartitions = std::max(partitions, unsigned(1));
+    mActivePart = std::min(activePart, mPartitions-1);
+}
+
+
+inline std::vector<unsigned char>&
+VolumeToMesh::pointFlags()
+{
+    return mPointFlags;
+}
+
+
+inline const std::vector<unsigned char>&
+VolumeToMesh::pointFlags() const
+{
+    return mPointFlags;
 }
 
 
@@ -2458,172 +4073,542 @@ inline void
 VolumeToMesh::operator()(const GridT& distGrid)
 {
     typedef typename GridT::TreeType DistTreeT;
+    typedef tree::LeafManager<const DistTreeT> DistLeafManagerT;
     typedef typename DistTreeT::ValueType DistValueT;
 
-    typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
-    typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
-    typedef typename DistTreeT::template ValueConverter<Vec3s>::Type Vec3sTreeT;
     typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
+    typedef tree::LeafManager<BoolTreeT> BoolLeafManagerT;
+    typedef Grid<BoolTreeT> BoolGridT;
+
+    typedef typename DistTreeT::template ValueConverter<Int16>::Type Int16TreeT;
+    typedef tree::LeafManager<Int16TreeT> Int16LeafManagerT;
+
+    typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
+    typedef typename DistTreeT::template ValueConverter<float>::Type FloatTreeT;
+    typedef Grid<FloatTreeT> FloatGridT;
 
-    const bool noAdaptivity = mPrimAdaptivity < 1e-6 && mSecAdaptivity < 1e-6;
 
     const openvdb::math::Transform& transform = distGrid.transform();
     const DistTreeT& distTree = distGrid.tree();
-    typename CharTreeT::Ptr edgeTree; // edge flags
-    typename IntTreeT::Ptr auxTree; // auxiliary data
+    const DistValueT isovalue = DistValueT(mIsovalue);
+
+    typename Int16TreeT::Ptr signTreePt;
+    typename IntTreeT::Ptr idxTreePt;
+    typename BoolTreeT::Ptr pointMask;
 
-    const bool nonLevelSetGrid = distGrid.getGridClass() != GRID_LEVEL_SET;
+    BoolTreeT valueMask(false), seamMask(false);
+    const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
+    bool maskEdges = false;
 
-    const bool extraSignCheck = nonLevelSetGrid ||
-        (std::abs(mIsovalue - double(distGrid.background())) < (1.5 * transform.voxelSize()[0]));
 
+    const BoolGridT * surfaceMask = NULL;
+    if (mSurfaceMaskGrid && mSurfaceMaskGrid->type() == BoolGridT::gridType()) {
+        surfaceMask = static_cast<const BoolGridT*>(mSurfaceMaskGrid.get());
+    }
+
+    const FloatGridT * adaptivityField = NULL;
+    if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridT::gridType()) {
+        adaptivityField = static_cast<const FloatGridT*>(mAdaptivityGrid.get());
+    }
+
+    if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeT::treeType()) {
+        const BoolTreeT *adaptivityMaskPt =
+            static_cast<const BoolTreeT*>(mAdaptivityMaskTree.get());
+        seamMask.topologyUnion(*adaptivityMaskPt);
+    }
 
 
     // Collect auxiliary data
     {
-        internal::LeafCPtrList<DistTreeT> sourceLeafs(distTree);
-        internal::AuxiliaryData<DistTreeT> op(distTree, sourceLeafs, mIsovalue, extraSignCheck);
-        op.runParallel();
-        edgeTree = op.edgeTree();
-        auxTree = op.auxTree();
+        DistLeafManagerT distLeafs(distTree);
 
-        // Collect auxiliary data from active tiles
-        if (nonLevelSetGrid) {
-            internal::tileAuxiliaryData(distTree, *edgeTree, *auxTree, mIsovalue);
+        // Check if the isovalue is in proximity to the active voxel boundary.
+        bool padActiveVoxels = false;
+        int padVoxels = 3;
+
+        if (distGrid.getGridClass() != GRID_LEVEL_SET) {
+            padActiveVoxels = true;
+        } else {
+            padActiveVoxels = internal::needsActiveVoxePadding(distLeafs,
+                mIsovalue, transform.voxelSize()[0]);
         }
-    }
 
+        // always pad the active region for small volumes (the performance hit is neglectable).
+        if (!padActiveVoxels) {
+            Coord dim;
+            distTree.evalActiveVoxelDim(dim);
+            int maxDim = std::max(std::max(dim[0], dim[1]), dim[2]);
+            if (maxDim < 1000) {
+                padActiveVoxels = true;
+                padVoxels = 1;
+            }
+        }
 
-    // Optionally collect auxiliary data from a reference surface.
-    internal::ReferenceData<DistTreeT> refData;
-    if (mRefGrid) {
-        const GridT* refGrid = static_cast<const GridT*>(mRefGrid.get());
-        if (refGrid && refGrid->activeVoxelCount() > 0) {
-
-            refData.mDistTree = &refGrid->tree();
-            refData.mInternalAdaptivity = DistValueT(mSecAdaptivity);
-
-            // Cache reference data for reuse.
-            if (!mRefEdgeTree && !mRefTopologyMaskTree) {
-                internal::LeafCPtrList<DistTreeT> leafs(*refData.mDistTree);
-                internal::AuxiliaryData<DistTreeT, bool> op(
-                    *refData.mDistTree, leafs, mIsovalue, extraSignCheck);
-                op.runParallel();
-                mRefEdgeTree = op.edgeTree();
-                mRefTopologyMaskTree = op.auxTree();
-                mSeamPointTree = typename Vec3sTreeT::Ptr(new Vec3sTreeT(Vec3s(0.0)));
+        if (surfaceMask || mPartitions > 1) {
+
+            maskEdges = true;
+
+            if (surfaceMask) {
+
+                { // Mask
+                    internal::GenTopologyMask<DistTreeT> masking(
+                        *surfaceMask, distLeafs, transform, mInvertSurfaceMask);
+                    masking.run();
+                    valueMask.merge(masking.tree());
+                }
+
+                if (mPartitions > 1) { // Partition
+                    tree::LeafManager<BoolTreeT> leafs(valueMask);
+                    leafs.foreach(internal::PartOp(leafs.leafCount() , mPartitions, mActivePart));
+                    valueMask.pruneInactive();
+                }
+
+            } else { // Partition
+
+                internal::PartGen<DistTreeT> partitioner(distLeafs, mPartitions, mActivePart);
+                partitioner.run();
+                valueMask.merge(partitioner.tree());
+            }
+
+            {
+                if (padActiveVoxels) tools::dilateVoxels(valueMask, padVoxels);
+                BoolLeafManagerT leafs(valueMask);
+
+                internal::SignData<DistTreeT, BoolLeafManagerT>
+                    signDataOp(distTree, leafs, isovalue);
+                signDataOp.run();
+
+                signTreePt = signDataOp.signTree();
+                idxTreePt = signDataOp.idxTree();
             }
 
-            if (mRefEdgeTree && mRefTopologyMaskTree) {
-                refData.mEdgeTree = static_cast<CharTreeT*>(mRefEdgeTree.get());
-                refData.mTopologyMaskTree =  static_cast<BoolTreeT*>(mRefTopologyMaskTree.get());
-                refData.mSeamPointTree = static_cast<Vec3sTreeT*>(mSeamPointTree.get());
-                refData.mSeamMaskTree = typename BoolTreeT::Ptr(new BoolTreeT(false));
-                refData.mSmoothingMaskTree = typename BoolTreeT::Ptr(new BoolTreeT(false));
+            {
+                internal::GenBoundaryMask<DistTreeT> boundary(distLeafs, valueMask, *idxTreePt);
+                boundary.run();
+
+                BoolLeafManagerT bleafs(boundary.tree());
+
+                internal::SignData<DistTreeT, BoolLeafManagerT>
+                    signDataOp(distTree, bleafs, isovalue);
+                signDataOp.run();
+
+                signTreePt->merge(*signDataOp.signTree());
+                idxTreePt->merge(*signDataOp.idxTree());
+            }
+
+        } else {
+
+            // Collect voxel-sign configurations
+            if (padActiveVoxels) {
+
+                BoolTreeT regionMask(false);
+                regionMask.topologyUnion(distTree);
+                tools::dilateVoxels(regionMask, padVoxels);
+
+                BoolLeafManagerT leafs(regionMask);
+
+                internal::SignData<DistTreeT, BoolLeafManagerT>
+                    signDataOp(distTree, leafs, isovalue);
+                signDataOp.run();
+
+                signTreePt = signDataOp.signTree();
+                idxTreePt = signDataOp.idxTree();
+            } else {
+
+                internal::SignData<DistTreeT, DistLeafManagerT>
+                    signDataOp(distTree, distLeafs, isovalue);
+                signDataOp.run();
+
+                signTreePt = signDataOp.signTree();
+                idxTreePt = signDataOp.idxTree();
             }
         }
+
     }
 
-    BoolTreeT edgeMaskTree(0.0);
 
-    // Generate the seamline mask
-    if (refData.mSeamMaskTree) {
-        refData.mSeamMaskTree->topologyUnion(*auxTree.get());
+    // Collect auxiliary data from active tiles
+    internal::tileData(distTree, *signTreePt, *idxTreePt, isovalue);
+
+    // Optionally collect auxiliary data from a reference level set.
+    Int16TreeT *refSignTreePt = NULL;
+    IntTreeT *refIdxTreePt = NULL;
+    const DistTreeT *refDistTreePt = NULL;
+
+    if (mRefGrid && mRefGrid->type() == GridT::gridType()) {
+
+        const GridT* refGrid = static_cast<const GridT*>(mRefGrid.get());
+        refDistTreePt = &refGrid->tree();
+
+        // Collect and cache auxiliary data from the reference grid.
+        if (!mRefSignTree && !mRefIdxTree) {
+
+            DistLeafManagerT refDistLeafs(*refDistTreePt);
+            internal::SignData<DistTreeT, DistLeafManagerT>
+                signDataOp(*refDistTreePt, refDistLeafs, isovalue);
 
-        internal::LeafPtrList<BoolTreeT> leafs(*refData.mSeamMaskTree.get());
-        internal::SeamMaskGen<DistTreeT> op(leafs, *refData.mTopologyMaskTree, *auxTree.get());
-        op.runParallel();
+            signDataOp.run();
 
-        refData.mSeamMaskTree->pruneInactive();
-        edgeMaskTree.topologyUnion(*refData.mSeamMaskTree);
-        dilateVoxels(*refData.mSeamMaskTree);
-        dilateVoxels(*refData.mSeamMaskTree);
-        dilateVoxels(*refData.mSeamMaskTree);
+            mRefSignTree = signDataOp.signTree();
+            mRefIdxTree = signDataOp.idxTree();
+        }
+
+        // Get cached auxiliary data
+        if (mRefSignTree && mRefIdxTree) {
+            refSignTreePt = static_cast<Int16TreeT*>(mRefSignTree.get());
+            refIdxTreePt = static_cast<IntTreeT*>(mRefIdxTree.get());
+        }
     }
 
 
     // Process auxiliary data
+    Int16LeafManagerT signLeafs(*signTreePt);
+
+    if (maskEdges) {
+        signLeafs.foreach(internal::MaskEdges<BoolTreeT>(valueMask));
+        valueMask.clear();
+    }
+
+
+    // Generate the seamline mask
+    if (refSignTreePt) {
+        internal::GenSeamMask<Int16TreeT, Int16LeafManagerT> seamOp(signLeafs, *refSignTreePt);
+        seamOp.run();
+
+        tools::dilateVoxels(seamOp.mask(), 3);
+        signLeafs.foreach(internal::TagSeamEdges<BoolTreeT>(seamOp.mask()));
+
+        seamMask.merge(seamOp.mask());
+    }
+
+
+    std::vector<size_t> regions(signLeafs.leafCount(), 0);
+    if (regions.empty()) return;
+
+    if (adaptive) {
+
+        internal::MergeVoxelRegions<DistTreeT, Int16LeafManagerT> merge(
+            signLeafs, *signTreePt, distTree, *idxTreePt, isovalue, DistValueT(mPrimAdaptivity));
+
+        if (adaptivityField) {
+            merge.setSpatialAdaptivity(transform, *adaptivityField);
+        }
+
+        if (refSignTreePt || mAdaptivityMaskTree) {
+            merge.setAdaptivityMask(&seamMask);
+        }
+
+        if (refSignTreePt) {
+            merge.setRefData(refSignTreePt, DistValueT(mSecAdaptivity));
+        }
+
+        merge.run();
+
+        signLeafs.foreach(internal::CountRegions<IntTreeT>(*idxTreePt, regions));
+
+    } else {
+        signLeafs.foreach(internal::CountPoints(regions));
+    }
+
+
     {
-        internal::LeafPtrList<IntTreeT> auxLeafs(*auxTree);
-        std::vector<size_t> regions(auxLeafs.size(), 0);
+        mPointListSize = 0;
+        size_t tmp = 0;
+        for (size_t n = 0, N = regions.size(); n < N; ++n) {
+            tmp = regions[n];
+            regions[n] = mPointListSize;
+            mPointListSize += tmp;
+        }
+    }
 
-        {
-            if (noAdaptivity) {
-                internal::Count<IntTreeT> count(auxLeafs, regions);
-                count.runParallel();
-            } else {
-                internal::Merge<DistTreeT> merge(distTree, auxLeafs,
-                    regions, DistValueT(mIsovalue), DistValueT(mPrimAdaptivity));
-                merge.setRefData(refData);
-                merge.runParallel();
+
+    // Generate the unique point list
+    mPoints.reset(new openvdb::Vec3s[mPointListSize]);
+    mPointFlags.clear();
+
+    // Generate seam line sample points
+    if (refSignTreePt && refIdxTreePt) {
+
+        if (mSeamPointListSize == 0) {
+
+            std::vector<size_t> pointMap;
+
+            {
+                Int16LeafManagerT refSignLeafs(*refSignTreePt);
+                pointMap.resize(refSignLeafs.leafCount(), 0);
+
+                refSignLeafs.foreach(internal::CountPoints(pointMap));
+
+                size_t tmp = 0;
+                for (size_t n = 0, N = pointMap.size(); n < N; ++n) {
+                    tmp = pointMap[n];
+                    pointMap[n] = mSeamPointListSize;
+                    mSeamPointListSize += tmp;
+                }
             }
 
-            mPointListSize = 0;
-            size_t tmp = 0;
-            for (size_t n = 0, N = regions.size(); n < N; ++n) {
-                tmp = regions[n];
-                regions[n] = mPointListSize;
-                mPointListSize += tmp;
+            if (!pointMap.empty() && mSeamPointListSize != 0) {
+
+                mQuantizedSeamPoints.reset(new uint32_t[mSeamPointListSize]);
+                memset(mQuantizedSeamPoints.get(), 0, sizeof(uint32_t) * mSeamPointListSize);
+
+                typedef tree::LeafManager<IntTreeT> IntLeafManagerT;
+
+                IntLeafManagerT refIdxLeafs(*refIdxTreePt);
+                refIdxLeafs.foreach(internal::MapPoints<Int16TreeT>(pointMap, *refSignTreePt));
             }
         }
 
-        if (refData.isValid()) { // match leaf topology
-            tree::ValueAccessor<BoolTreeT> acc(*refData.mSmoothingMaskTree);
-            for (size_t n = 0, N = auxLeafs.size(); n < N; ++n) {
-                acc.touchLeaf(auxLeafs[n]->getOrigin());
-            }
+        if (mSeamPointListSize != 0) {
+            signLeafs.foreach(internal::SeamWeights<DistTreeT>(
+                distTree, *refSignTreePt, *refIdxTreePt, mQuantizedSeamPoints, mIsovalue));
+        }
+    }
+
+
+    internal::GenPoints<DistTreeT, Int16LeafManagerT>
+        pointOp(signLeafs, distTree, *idxTreePt, mPoints, regions, transform, mIsovalue);
+
+
+    if (mSeamPointListSize != 0) {
+        mPointFlags.resize(mPointListSize);
+        pointOp.setRefData(refSignTreePt, refDistTreePt, refIdxTreePt,
+            &mQuantizedSeamPoints, &mPointFlags);
+    }
+
+    pointOp.run();
+
+
+    mPolygonPoolListSize = signLeafs.leafCount();
+    mPolygons.reset(new PolygonPool[mPolygonPoolListSize]);
+
+
+    if (adaptive) {
+
+        internal::GenPolygons<Int16LeafManagerT, internal::AdaptivePrimBuilder>
+            mesher(signLeafs, *signTreePt, *idxTreePt, mPolygons, Index32(mPointListSize));
+
+        mesher.setRefSignTree(refSignTreePt);
+        mesher.run();
+
+    } else {
+
+        internal::GenPolygons<Int16LeafManagerT, internal::UniformPrimBuilder>
+            mesher(signLeafs, *signTreePt, *idxTreePt, mPolygons, Index32(mPointListSize));
+
+        mesher.setRefSignTree(refSignTreePt);
+        mesher.run();
+    }
+
+    // Clean up unused points, only necessary if masking and/or
+    // automatic mesh partitioning is enabled.
+    if ((surfaceMask || mPartitions > 1) && mPointListSize > 0) {
+
+        // Flag used points
+        std::vector<unsigned char> usedPointMask(mPointListSize, 0);
+
+        internal::FlagUsedPoints flagPoints(mPolygons, mPolygonPoolListSize, usedPointMask);
+        flagPoints.run();
+
+        // Create index map
+        std::vector<unsigned> indexMap(mPointListSize);
+        size_t usedPointCount = 0;
+        for (size_t p = 0; p < mPointListSize; ++p) {
+            if (usedPointMask[p]) indexMap[p] = usedPointCount++;
         }
 
-        // Generate the unique point list
-        mPoints.reset(new openvdb::Vec3s[mPointListSize]);
+        if (usedPointCount < mPointListSize) {
+
+            // move points
+            internal::UniquePtr<openvdb::Vec3s>::type
+                newPointList(new openvdb::Vec3s[usedPointCount]);
 
-        internal::PointGen<DistTreeT>
-            pointGen(distTree, auxLeafs, regions, transform, mPoints, mIsovalue);
-        pointGen.setRefData(refData);
-        pointGen.runParallel();
+            internal::MovePoints movePoints(newPointList, mPoints, indexMap, usedPointMask);
+            movePoints.run();
+
+            mPointListSize = usedPointCount;
+            mPoints.reset(newPointList.release());
+
+            // update primitives
+            internal::RemapIndices remap(mPolygons, mPolygonPoolListSize, indexMap);
+            remap.run();
+        }
     }
 
-    // Smooth seam line edges
-    if (mSmoothSeams && refData.isValid()) {
-        refData.mSmoothingMaskTree->pruneInactive();
-        internal::LeafPtrList<BoolTreeT> leafs(*refData.mSmoothingMaskTree);
-        internal::EdgeSmooth<DistTreeT> op(leafs, edgeMaskTree, *auxTree, mPoints, transform);
-        op.runParallel(3);
 
-        // Cache shared points
-        tree::ValueAccessor<Vec3sTreeT> ptnAcc(*refData.mSeamPointTree);
-        tree::ValueAccessor<IntTreeT> auxAcc(*auxTree);
-        Coord ijk;
+    // Subdivide nonplanar quads near the seamline edges
+    // todo: thread and clean up
+    if (refSignTreePt || refIdxTreePt || refDistTreePt) {
+        std::vector<Vec3s> newPoints;
+
+        for (size_t n = 0; n <  mPolygonPoolListSize; ++n) {
+
+            PolygonPool& polygons = mPolygons[n];
+
+            std::vector<size_t> nonPlanarQuads;
+            nonPlanarQuads.reserve(polygons.numQuads());
+
+            for (size_t i = 0; i < polygons.numQuads(); ++i) {
+
+                char& flags = polygons.quadFlags(i);
+
+                if ((flags & POLYFLAG_FRACTURE_SEAM) && !(flags & POLYFLAG_EXTERIOR)) {
 
-        typename BoolTreeT::LeafIter leafIt = refData.mSeamMaskTree->beginLeaf();
-        for ( ; leafIt; ++leafIt) {
-            typename BoolTreeT::LeafNodeType::ValueOnIter valueIt = leafIt->beginValueOn();
-            for ( ; valueIt; ++valueIt) {
-                ijk = valueIt.getCoord();
-                const int idx = auxAcc.getValue(ijk);
-                if (idx != 0 && !ptnAcc.isValueOn(ijk)) {
-                   ptnAcc.setValue(ijk, mPoints[idx]);
+                    openvdb::Vec4I& quad = polygons.quad(i);
+
+                    const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
+                        || mPointFlags[quad[2]] || mPointFlags[quad[3]];
+
+                    if (!edgePoly) continue;
+
+                    const Vec3s& p0 = mPoints[quad[0]];
+                    const Vec3s& p1 = mPoints[quad[1]];
+                    const Vec3s& p2 = mPoints[quad[2]];
+                    const Vec3s& p3 = mPoints[quad[3]];
+
+                    if (!internal::isPlanarQuad(p0, p1, p2, p3, 1e-6f)) {
+                        nonPlanarQuads.push_back(i);
+                    }
                 }
             }
+
+
+            if (!nonPlanarQuads.empty()) {
+
+                PolygonPool tmpPolygons;
+
+                tmpPolygons.resetQuads(polygons.numQuads() - nonPlanarQuads.size());
+                tmpPolygons.resetTriangles(polygons.numTriangles() + 4 * nonPlanarQuads.size());
+
+                size_t triangleIdx = 0;
+                for (size_t i = 0; i < nonPlanarQuads.size(); ++i) {
+
+                    size_t& quadIdx = nonPlanarQuads[i];
+
+                    openvdb::Vec4I& quad = polygons.quad(quadIdx);
+                    char& quadFlags = polygons.quadFlags(quadIdx);
+                    //quadFlags |= POLYFLAG_SUBDIVIDED;
+
+                    Vec3s centroid = (mPoints[quad[0]] + mPoints[quad[1]] +
+                        mPoints[quad[2]] + mPoints[quad[3]]) * 0.25;
+
+                    size_t pointIdx = newPoints.size() + mPointListSize;
+
+                    newPoints.push_back(centroid);
+
+
+                    {
+                        Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
+
+                        triangle[0] = quad[0];
+                        triangle[1] = pointIdx;
+                        triangle[2] = quad[3];
+
+                        tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
+
+                        if (mPointFlags[triangle[0]] || mPointFlags[triangle[2]]) {
+                            tmpPolygons.triangleFlags(triangleIdx) |= POLYFLAG_SUBDIVIDED;
+                        }
+                    }
+
+                    ++triangleIdx;
+
+                    {
+                        Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
+
+                        triangle[0] = quad[0];
+                        triangle[1] = quad[1];
+                        triangle[2] = pointIdx;
+
+                        tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
+
+                        if (mPointFlags[triangle[0]] || mPointFlags[triangle[1]]) {
+                            tmpPolygons.triangleFlags(triangleIdx) |= POLYFLAG_SUBDIVIDED;
+                        }
+                    }
+
+                    ++triangleIdx;
+
+                    {
+                        Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
+
+                        triangle[0] = quad[1];
+                        triangle[1] = quad[2];
+                        triangle[2] = pointIdx;
+
+                        tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
+
+                        if (mPointFlags[triangle[0]] || mPointFlags[triangle[1]]) {
+                            tmpPolygons.triangleFlags(triangleIdx) |= POLYFLAG_SUBDIVIDED;
+                        }
+                    }
+
+
+                    ++triangleIdx;
+
+                    {
+                        Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
+
+                        triangle[0] = quad[2];
+                        triangle[1] = quad[3];
+                        triangle[2] = pointIdx;
+
+                        tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
+
+                        if (mPointFlags[triangle[0]] || mPointFlags[triangle[1]]) {
+                            tmpPolygons.triangleFlags(triangleIdx) |= POLYFLAG_SUBDIVIDED;
+                        }
+                    }
+
+                    ++triangleIdx;
+
+                    quad[0] = util::INVALID_IDX;
+                }
+
+
+                for (size_t i = 0; i < polygons.numTriangles(); ++i) {
+                    tmpPolygons.triangle(triangleIdx) = polygons.triangle(i);
+                    tmpPolygons.triangleFlags(triangleIdx) = polygons.triangleFlags(i);
+                    ++triangleIdx;
+                }
+
+
+                size_t quadIdx = 0;
+                for (size_t i = 0; i < polygons.numQuads(); ++i) {
+                    openvdb::Vec4I& quad = polygons.quad(i);
+
+                    if (quad[0] != util::INVALID_IDX) {
+                        tmpPolygons.quad(quadIdx) = quad;
+                        tmpPolygons.quadFlags(quadIdx) = polygons.quadFlags(i);
+                        ++quadIdx;
+                    }
+                }
+
+
+                polygons.copy(tmpPolygons);
+            }
+
         }
-    }
 
-    // Generate mesh
-    {
-        internal::LeafCPtrList<CharTreeT> edgeLeafs(*edgeTree);
-        mPolygonPoolListSize = edgeLeafs.size();
-        mPolygons.reset(new PolygonPool[mPolygonPoolListSize]);
-
-        if (noAdaptivity) {
-            internal::MeshGen<DistTreeT, internal::QuadMeshOp>
-                meshGen(edgeLeafs, *auxTree, mPolygons);
-            meshGen.setRefData(refData);
-            meshGen.runParallel();
-        } else {
-            internal::MeshGen<DistTreeT, internal::AdaptiveMeshOp>
-                meshGen(edgeLeafs, *auxTree, mPolygons);
-            meshGen.setRefData(refData);
-            meshGen.runParallel();
+
+        if (!newPoints.empty()) {
+
+            size_t newPointCount = newPoints.size() + mPointListSize;
+
+            internal::UniquePtr<openvdb::Vec3s>::type
+                newPointList(new openvdb::Vec3s[newPointCount]);
+
+            for (size_t i = 0; i < mPointListSize; ++i) {
+                newPointList.get()[i] = mPoints[i];
+            }
+
+            for (size_t i = mPointListSize; i < newPointCount; ++i) {
+                newPointList.get()[i] = newPoints[i - mPointListSize];
+            }
+
+            mPointListSize = newPointCount;
+            mPoints.reset(newPointList.release());
+            mPointFlags.resize(mPointListSize, 0);
         }
     }
 }
@@ -2699,6 +4684,10 @@ volumeToMesh(
     volumeToMesh(grid,points, triangles, quads, isovalue, 0.0);
 }
 
+
+////////////////////////////////////////
+
+
 } // namespace tools
 } // namespace OPENVDB_VERSION_NAME
 } // namespace openvdb