Monday, 12 November 2012

V-HACD: Replacing Triangle's Constrained Delaunay Triangulation

As mentioned in my previous post, the first version of V-HACD relies on the library Triangle to compute 2D constrained Delaunay triangulations. In fact, the V-HACD algorithm involves clipping the mesh against a plane and filling the produced holes (cf. Figure 1).

The Triangle library produces excellent results and is very stable. However, as pointed out by Erwin Coumans it has a non permissive license. Another alternative to Triangle is the poly2tri, which is released under a New BSD License. Poly2tri turned out to be non-usable in my case because it supports only simple polygons as constraints. After searching the internet for a "good" and stable implementation with a permissive license, I ended up developing a simple triangulation algorithm (it is not a Constrained Delaunay Triangulation), which is enough for V-HACD needs. My implementation is a very simplified version of the Incremental Delaunay algorithm 

The new triangulation algorithm was uploaded to the git repository and could be compared to the Triangle library by adding #define USE_TRIANGLE to vhacdMesh.cpp

The code is still ugly and not optimized at all. For now, I am trying to have some first results and understand better the algorithm limitations. Hopefully, I'll have soon the time to clean up the code.

Wednesday, 7 November 2012


This post is out-of-date. Check out V-HACD 2.0!



V-HACD: Hierarchical Approximate Convex Decomposition Revisited

Lately, I have found time to work on improving the HACD algorithm. The new V-HACD library tries to tackle the problem of convex-hulls inter-penetration usually encountered when using HACD. Figure 1 illustrates this limitation by comparing V-HACD results to those generated by using HACD.



Figure 1. V-HACD vs. HACD: V-HACD generates non-overlapping convex-hulls.

V-HACD works only with manifold closed meshes of arbitrary genus, which makes it less general than HACD algorithm (which handles open and non-manifold meshes). The V-HACD code  is available under New BSD License. However, the current version relies on the following other libraries:

  • Triangle to compute constrained delaunay triangulation which has a non permissive license, which should be replaced [Thanks to Erwin Coumans for noticing that]
  • Ole Kniemeyer's implementation of the convex-hull algorithm provided with bullet
  • John Tsiombikas's kd-tree algorithm

The first results are encouraging. However, the code is still buggy and not optimized. I hope that people would be interested in helping me improve this first version.

To play with the algorithm, a pre-compiled win32 executable is available here
Pre-computed decomposition results are available here


V-HACD parameters are the following:
testVHACD fileName.off depth maxConcavity invertInputFaces posSampling angleSampling posRefine angleRefine alpha targetNTrianglesDecimatedMesh
  • fileName.off: 3D mesh in off format (type: string, example: block.off )
  • depth: maximum number of decomposition stages (type: integer, default: 10)
  • maxConcavity: maximum allowed concavity (type: float, default 0.01)
  • invertInputFaces: indicates whether mesh normals should be inverted or not (type: boolean, default 0)
  • posSampling: clipping plane position sampling resolution for coarse search (type: int, default 10)
  • angleSampling: clipping plane orientation sampling resolution for coarse search (type: int, default 10) 
  • posRefine: clipping plane position sampling resolution for refined search (type: int, default 5) 
  • angleRefine: clipping plane orientation sampling resolution for refined search (type: int, default 5)
  • alpha: parameter controlling the compromise between concavity and balance between convex-hulls. (type: float, default: 0.01)
  • targetNTrianglesDecimatedMesh: number of triangles in the decimated mesh. V-HACD decimates the input mesh to reduce computation times. (type: integer, default: 1000)

To apply V-HACD to the input mesh "input.off" use the following command line:

testVHACD.exe input.off 30 0.01 0 64 32 8 64 0.001 2000

Below, some screen-shots of the obtained results.





The current version of V-HACD provides poor results or errors for the following three meshes. As soon as I get some free time I'll try to fix these bugs.

  • dancer2
  • elk
  • hand1
  • hand2
  • hero
  • octopus
  • polygirl
  • shark_b
  • Sketched-Brunnen
  • torus
  • tstTorusModel
  • tstTorusModel3

Wednesday, 23 November 2011

HACD parameters

In this post, I'll to give an overview of the HACD parameters and explain their meaning and how they should be set. The text will be improved over time. My main concern is to have things written down for reference...


  1. Parameters overview
Parameter
Description
Default
NTargetTrianglesDecimatedMesh
Target number of triangles in the decimated mesh. The decimation stage was added mainly to decrease the computation costs for dense meshes (refer to Section 2.1).
1000
NVerticesPerCH
Maximum number of vertices in the generated convex-hulls.
100
ScaleFactor
Normalization factor used to ensure that the other parameters (e.g. concavity) are expressed w.r.t. a fixed size. Refer to Section 2.3 for details

1000
 SmallClusterThreshold
Threshold on the clusters area (expressed as a percentage of the entire mesh area) under which the cluster is considered small and it is forced to be merged with other clusters at the price of a high concavity.
0.25
AddFacesPoints
If enabled an additional ray located at the center of each triangle pointing toward its normal is considered when computing the concavity of a non-flat cluster. The parameter was added to handle coarse meshes (i.e. with a low number of vertices)

ON
AddExtraDistPoitns
If enabled additional rays are considered to handle bowl-like shapes.
ON
NClusters 
Minimum number of convex-hulls to be generated
1
Concavity
Maximum allowed concavity
100
ConnectDist
If the distance between two triangles, each belonging to a different connected components (CCs), is lower than the ConnectDist threshold an additional edge connecting them is added to the dual graph. This parameter was added to handle meshes with multiple CCs.
30
VolumeWeight
Weight controlling the contribution of the volume related cost to the global edgecollapse cost (refer to XXX for details).
0.0
(not used)
CompacityWeight
Weight controlling the contribution of the shape factor related cost to the global edgecollapse cost (refer to XXX for details).  
0.0001

FlatRegionThreshold
Threshold expressed a percentage of ScaleFactor under which a cluster is considered flat.
1
ComputationWeight

Weight controlling the contribution of the computation related cost to the global edgecollapse cost (refer to XXX for details).
0.01
2. Detailed description
  • NTargetTrianglesDecimatedMesh
In order to reduce the computations times for dense meshes, the HACD library makes it possible to decimate the original mesh before running the decomposition process. More details about the decimation algorithm are provided here http://kmamou.blogspot.com/2011/10/hacd-optimization.html

  • NVerticesPerCH
this parameter was introduced in order to comply with the constraints that most physics engines put on the number of vertices/triangles per convex-hull (CH). If the function HACD::Compute() is called with the parameter fullCH=false, then the generated CHs will have a number of vertices lower than NVerticesPerCH.

In order to optimally choose the best vertices to keep in the final CH, the ICHull::Process(unsigned long nPointsCH) function implements a slightly different version of the Incremental Convex Hull algorithm (cf. demo code ). Here, at each step, the point with the highest volume increment  is chosen, until all points are processed or the CH has exactly NVerticesPerCH points.

The code looks like this:

        while (!vertices.GetData().m_tag && addedPoints < nPointsCH) // not processed
        {
            if (!FindMaxVolumePoint())
            {
                break;
            }                 
            vertices.GetData().m_tag = true;                     
            if (ProcessPoint())
            {
                addedPoints++;
                CleanUp(addedPoints);
                vertices.Next();
            }
        }
        // delete remaining points
        while (!vertices.GetData().m_tag)
        {
            vertices.Delete();
        }

  • ScaleFactor
A normalization process is applied to the input mesh in order to ensure that the other parameters (e.g. concavity) are expressed w.r.t. a fixed size. This process is inverted before producing the final CHs.

The HACD::Compute() function follows the following main steps:

    bool HACD::Compute(bool fullCH, bool exportDistPoints)
    {
        if (m_targetNTrianglesDecimatedMesh > 0)
        {
           DecimateMesh(targetNTrianglesDecimatedMesh);
        }
        NormalizeData();
        CreateGraph();
        InitializeDualGraph();
        InitializePriorityQueue();
        Simplify();
        DenormalizeData();
        CreateFinalCH();
        return true;
    }

The HACD::NormalizeData() function centers the mesh and scale it so its coordinates are in the interval [-m_sacle, m_sclae]x[-m_sacle, m_sclae]x[-m_sacle, m_sclae]. The code proceeds as follows:

       void HACD::NormalizeData()
       {
              const Real invDiag = static_cast<Real>(2.0 * m_scale / m_diag);
              for (size_t v = 0; v < m_nPoints ; v++)
              {
                     m_points[v] = (m_points[v] - m_barycenter) * invDiag;
              }
       }

The HACD::DenormalizeData() function invert the normalization operated by HACD::NormalizeData():
       void HACD::DenormalizeData()
       {
              const Real diag = static_cast<Real>(m_diag / (2.0 * m_scale));
              for (size_t v = 0; v < m_nPoints ; v++)
              {
                     m_points[v] = m_points[v] * diag + m_barycenter;
              }
      }
  • SmallClusterThreshold
Due to numerical stability issues (or maybe some bugs I haven't spotted yet :) ) the HACD algorithm may produce small clusters. In order to detect them and make sure they will be merged, the SmallClusterThreshold was introduced. A cluster is considered to be small if its area is smaller than SmallClusterThreshold% of the entire mesh area.

The condition2 in the HACD::Simplify() function forces small clusters to be merged (m_area is the area of the entire mesh):
void HACD::Simplify()
       {
              double areaThreshold = m_area * m_smallClusterThreshold / 100.0;
              while ( !m_pqueue.empty() )
              {
               currentEdge = m_pqueue.top();
               m_pqueue.pop();
               v1 = m_graph.m_edges[currentEdge.m_name].m_v1;
               v2 = m_graph.m_edges[currentEdge.m_name].m_v2;
               bool condition1 = (m_graph.m_edges[currentEdge.m_name].m_concavity <  
                                  m_concavity) && (globalConcavity < m_concavity) && 
                                  (m_graph.GetNVertices() > m_nMinClusters) && 
                                  (m_graph.GetNEdges() > 0);
              bool condition2 = (m_graph.m_vertices[v1].m_surf < areaThreshold || 
                                 m_graph.m_vertices[v2].m_surf < areaThreshold);
              if (condition1 || condition2)
              {
                   m_graph.EdgeCollapse(v1, v2);
                   long idEdge;
                   for(size_t itE = 0; itE < m_graph.m_vertices[v1].m_edges.Size(); ++itE)
                   {
                      idEdge = m_graph.m_vertices[v1].m_edges[itE];
                      ComputeEdgeCost(idEdge);                                            
                   }                 
              }

       }

  • AddFacesPoints and AddExtraDistPoitns
The parameter AddFacesPoints was introduced to improve the precision of the concavity computation for meshes with a low number of vertices. The idea is to add additional rays located each at the center  of a triangle and pointing to the same direction as its normal.

The Parameter AddExtraDistPoints was added to handle bowl-like shapes. As illustrated below, if only the rays located on the current cluster are considered when computing its concavity, you may end up with a big cluster corresponding to the external surface (which is convex) of the bowl containing a lot of small clusters located on the internal part (which is concave). 

Bad decomposition for bowl-like shapes if AddExtraDistPoints is not enabled
In order to avoid such a bad decomposition, the idea consists in introducing new rays that would constrain the propagation of the cluster corresponding to the external surface of the bowl by taking into account rays located on the concave part. More precisely, during the initialization stage, an additional ray (the yellow ray in the figure below) is associated with each triangle (colored in red). 

The additional ray, denoted R (the yellow arrow), is defined as follows. Let N be the normal (the blue arrow) to the current triangle T (colored in red) and X be the ray starting at the center of T and with direction (-N) (the dotted green arrow). The starting point of R, denoted P0 (the yellow point), is defined a the nearest intersection point of X and the mesh. Moreover, the normal of the surface at P0 should point to the same direction as X. R has the direction of the normal of the surface at the P0.

Additional ray (yellow) is associated with the red triangle when AddExtraDistPoints is activated


The code looks like this:

    void HACD::InitializeDualGraph()
    {
         for(unsigned long f = 0; f < m_nTriangles; f++)
        {
            i = m_triangles[f].X();
            j = m_triangles[f].Y();
            k = m_triangles[f].Z();
            m_graph.m_vertices[f].m_distPoints.PushBack(DPoint(i, 0, false, false));
            m_graph.m_vertices[f].m_distPoints.PushBack(DPoint(j, 0, false, false));
            m_graph.m_vertices[f].m_distPoints.PushBack(DPoint(k, 0, false, false));
           
            u = m_points[j] - m_points[i];
            v = m_points[k] - m_points[i];
            w = m_points[k] - m_points[j];
            normal = u ^ v;

            m_normals[i] += normal;
            m_normals[j] += normal;
            m_normals[k] += normal;

            m_graph.m_vertices[f].m_surf = normal.GetNorm();
            m_area += m_graph.m_vertices[f].m_surf;
            normal.Normalize();
            if(m_addFacesPoints)
            {
                m_faceNormals[f] = normal;
                m_facePoints[f] = (m_points[i] + m_points[j] + m_points[k]) / 3.0;
            }
            if (m_addExtraDistPoints)   
            {
                Vec3<Real> seedPoint((m_points[i] + m_points[j] + m_points[k]) / 3.0);
                Vec3<Real> hitPoint;
                Vec3<Real> hitNormal;
                normal = -normal;
                if (rm.Raycast(seedPoint,normal,hitTriangle,dist, hitPoint, hitNormal))
                {
                     faceIndex = hitTriangle;
                }  
                if (faceIndex < m_nTriangles )
                {
                     m_extraDistPoints[f] = hitPoint;
                     m_extraDistNormals[f] = hitNormal;                
                     m_graph.m_vertices[f].m_distPoints.PushBack(DPoint(m_nPoints+f, 0, false, true));
              }
            }
        }
        for (size_t v = 0; v < m_nPoints; v++)
        {
              m_normals[v].Normalize();
        }
    }
  • Concavity
This parameter specifies the maximum allowed concavity for each cluster. As discussed in http://kmamou.blogspot.com/2011/10/hacd-hierarchical-approximate-convex.html different concavity measures are considered for flat (i.e. 2D) surfaces and for non-flat surfaces.
  • ConnectDist
In order to handle meshes with different connected components. The idea consists  in adding "virtual edges" between triangles belonging to different CCs. More precisely, if the distance between two triangles T1 and T2 belonging each to a different CC is lower than a threshold distConnect then an edge connecting T1 to T2 is added to the dual graph.

The code looks like this:
       void HACD::CreateGraph()
    {
        …
        if (m_ccConnectDist >= 0.0)
        {
            m_graph.ExtractCCs();
            if (m_graph.m_nCCs > 1)
            {
                std::vector< std::set<long> > cc2V;
                cc2V.resize(m_graph.m_nCCs);
                long cc;
                for(size_t t = 0; t < m_nTriangles; ++t)
                {
                    cc = m_graph.m_vertices[t].m_cc;
                    cc2V[cc].insert(m_triangles[t].X());
                    cc2V[cc].insert(m_triangles[t].Y());
                    cc2V[cc].insert(m_triangles[t].Z());
                }
               
                for(size_t cc1 = 0; cc1 < m_graph.m_nCCs; ++cc1)
                {
                    for(size_t cc2 = cc1+1; cc2 < m_graph.m_nCCs; ++cc2)
                    {
                        std::set<long>::const_iterator itV1(cc2V[cc1].begin()), itVEnd1(cc2V[cc1].end());
                        for(; itV1 != itVEnd1; ++itV1)
                        {
                                                double distC1C2 = std::numeric_limits<double>::max();
                            double dist;
                            t1 = -1;
                            t2 = -1;
                            std::set<long>::const_iterator itV2(cc2V[cc2].begin()), itVEnd2(cc2V[cc2].end());
                            for(; itV2 != itVEnd2; ++itV2)
                            {
                                dist = (m_points[*itV1] - m_points[*itV2]).GetNorm();
                                if (dist < distC1C2)
                                {
                                    distC1C2 = dist;
                                    t1 = *vertexToTriangles[*itV1].begin();
                                   
                                                              std::set<long>::const_iterator it2(vertexToTriangles[*itV2].begin()),
                                                                                                                 it2End(vertexToTriangles[*itV2].end());
                                                              t2 = -1;
                                                              for(; it2 != it2End; ++it2)
                                                              {
                                                                     if (*it2 != t1)
                                                                     {
                                                                           t2 = *it2;
                                                                           break;
                                                                     }
                                                              }
                                }
                            }
                            if (distC1C2 <= m_ccConnectDist && t1 >= 0 && t2 >= 0)
                            {
                                m_graph.AddEdge(t1, t2);                   
                            }
                        }
                    }
                }
            }
        }
    }
  •  FlatRegionThreshold
When computing concavity we need to distinguish between flat surfaces and non-float surfaces. The measure of flatness considered in HACD is related to the ratio between the convex-hull volume and the area of its boundary. If this later ration is small compared to m_scale then the mesh is considered flat. Otherwise it is considered non-flat. The parameter FlatRegionThreshold is the threshold which separate flat region from non flat region. It is expressed as percentage of m_scale.

In practice, the final concavity is computed as weighted sum of the flat region concavity and the 3D concavity. The weight is a function of the flatness of the cluster.

The code is as follows:

void HACD::ComputeEdgeCost(size_t e)
{
double surfCH = ch->ComputeArea() / 2.0;
double volumeCH = ch->ComputeVolume();
double vol2Surf = volumeCH / surfCH;
double concavity_flat = sqrt(fabs(surfCH-surf));
double weightFlat = std::max(0.0, 1.0 - pow(- vol2Surf * 100.0 / (m_scale * m_flatRegionThreshold), 2.0));
concavity_flat *= weightFlat;
if(!ch->IsFlat())
{
   concavity = Concavity(*ch, distPoints);
}
concavity += concavity_flat;
}

Wednesday, 5 October 2011

HACD optimization

Today, I have updated the HACD library in order to reduce both memory usage and computation complexity (cf. http://sourceforge.net/projects/hacd/). 

The new version:

  • uses John's Micro Allocator to avoid intensive dynamic memory allocation (thanks John for the great work),
  • exploits a simplified version of  an axis-aligned-bounding-volume AABB tree to accelerate dual graph generation (the code is inspired by John's post on this subject, thanks again John :) )
  • has an integrated mesh simplification pre-processing step, which makes it possible to handle dense meshes.
In this post, I'll give more details about this last feature.

To enable mesh decimation, the user specifies the target number of triangles that should be produced before running the convex decomposition. HACD will handle automatically the simplification and the decomposition processes. For the details of the mesh decimation algorithm, have a look at Michael Garland's  webpage http://mgarland.org/research/thesis.html. To turn this feature off just set the parameter targetNTrianglesDecimatedMesh=0.

I have tested the updated HACD algorithm on the 3D model "Bunny", which has 70K triangles. The HACD's parameters were set as follows:
  • minNClusters = 12
  • maxConcavity = 1000
  • addExtraDistPoints = 1
  • addFacesPoints = 1
  • ConnectedComponentsDist = 30
  • targetNTrianglesDecimatedMesh = 500, 1000, 2000, 4000, 8000 and 16000.

The snapshots of the produced convex decompositions are reported below. The computation times on my machine (Mac Intel Core 2 Duo, 4 GB RAM DDR3) ranged between 3 sec. for 500 triangles and 200 sec. for 16000 triangles (cf. Table 1 for the details).

----------------------------------------
# triangles  Time (sec.)
----------------------------------------
500 2.8
1000 4.6
2000 10
4000 21
8000 70
16000 192
----------------------------------------
Table 1. Computation times for targetNTrianglesDecimatedMesh= 500, 1000, 2000, 4000, 8000 and 16000 
# triangles after simplification - 500

# triangles after simplification - 1000

# triangles after simplification - 2000

# triangles after simplification - 4000

# triangles after simplification - 8000

# triangles after simplification - 16000