ENH: surfaceClean: have minimum quality

applications/utilities/surface/surfaceClean/collapseBase.C

@@ -66,8 +66,8 @@ static void writeRegionOBJ
 
     triSurface regionSurf(surf.subsetMesh(include, pointMap, faceMap));
 
-    //Pout<< "Region " << regionI << " surface:" << nl;
+    Pout<< "Region " << regionI << " surface:" << nl;
-    //regionSurf.writeStats(Pout);
+    regionSurf.writeStats(Pout);
 
     regionSurf.write(regionName);
 
@@ -97,7 +97,7 @@ static void splitTri
     DynamicList<labelledTri>& tris
 )
 {
-    label oldNTris = tris.size();
+    //label oldNTris = tris.size();
 
     label fp = findIndex(f, e[0]);
     label fp1 = f.fcIndex(fp);
@@ -175,7 +175,7 @@ static void splitTri
     }
     else
     {
-        FatalErrorIn("splitTri")
+        FatalErrorIn("splitTri(..)")
             << "Edge " << e << " not part of triangle " << f
             << " fp:" << fp
             << " fp1:" << fp1
@@ -183,13 +183,13 @@ static void splitTri
             << abort(FatalError);
     }
 
-    Pout<< "Split face " << f << " along edge " << e
+    //Pout<< "Split face " << f << " along edge " << e
-        << " into triangles:" << endl;
+    //    << " into triangles:" << endl;
-
+    //
-    for (label i = oldNTris; i < tris.size(); i++)
+    //for (label i = oldNTris; i < tris.size(); i++)
-    {
+    //{
-        Pout<< "   " << tris[i] << nl;
+    //    Pout<< "   " << tris[i] << nl;
-    }
+    //}
 }
 
 
@@ -206,14 +206,14 @@ static bool insertSorted
 {
     if (findIndex(sortedVerts, vertI) != -1)
     {
-        FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
+        FatalErrorIn("insertSorted(..)") << "Trying to insert vertex " << vertI
             << " which is already in list of sorted vertices "
             << sortedVerts << abort(FatalError);
     }
 
     if (weight <= 0 || weight >= 1)
     {
-        FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
+        FatalErrorIn("insertSorted(..)") << "Trying to insert vertex " << vertI
             << " with illegal weight " << weight
             << " into list of sorted vertices "
             << sortedVerts << abort(FatalError);
@@ -228,7 +228,7 @@ static bool insertSorted
 
         if (mag(w - weight) < SMALL)
         {
-            WarningIn("insertSorted")
+            WarningIn("insertSorted(..)")
                 << "Trying to insert weight " << weight << " which is close to"
                 << " existing weight " << w << " in " << sortedWeights
                 << endl;
@@ -263,64 +263,103 @@ static bool insertSorted
 }
 
 
+// Is triangle candidate for collapse? Small height or small quality
+bool isSliver
+(
+    const triSurface& surf,
+    const scalar minLen,
+    const scalar minQuality,
+    const label faceI,
+    const label edgeI
+)
+{
+    const pointField& localPoints = surf.localPoints();
+
+    // Check
+    // - opposite vertex projects onto base edge
+    // - normal distance is small
+    // - or triangle quality is small
+
+    label opposite0 =
+        triSurfaceTools::oppositeVertex
+        (
+            surf,
+            faceI,
+            edgeI
+        );
+
+    const edge& e = surf.edges()[edgeI];
+    const labelledTri& f = surf[faceI];
+
+    pointHit pHit =
+        e.line(localPoints).nearestDist
+        (
+            localPoints[opposite0]
+        );
+
+    if
+    (
+        pHit.hit()
+     && (
+            pHit.distance() < minLen
+         || f.tri(surf.points()).quality() < minQuality
+        )
+    )
+    {
+        // Remove faceI and split all other faces using this
+        // edge. This is done by 'replacing' the edgeI with the
+        // opposite0 vertex
+        //Pout<< "Splitting face " << faceI << " since distance "
+        //    << pHit.distance()
+        //    << " from vertex " << opposite0
+        //    << " to edge " << edgeI
+        //    << "  points "
+        //    << localPoints[e[0]]
+        //    << localPoints[e[1]]
+        //    << " is too small or triangle quality "
+        //    << f.tri(surf.points()).quality()
+        //    << " too small." << endl;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 // Mark all faces that are going to be collapsed.
 // faceToEdge: per face -1 or the base edge of the face.
 static void markCollapsedFaces
 (
     const triSurface& surf,
     const scalar minLen,
+    const scalar minQuality,
     labelList& faceToEdge
 )
 {
     faceToEdge.setSize(surf.size());
     faceToEdge = -1;
 
-    const pointField& localPoints = surf.localPoints();
     const labelListList& edgeFaces = surf.edgeFaces();
 
     forAll(edgeFaces, edgeI)
     {
-        const edge& e = surf.edges()[edgeI];
-
         const labelList& eFaces = surf.edgeFaces()[edgeI];
 
         forAll(eFaces, i)
         {
             label faceI = eFaces[i];
 
-            // Check distance of vertex to edge.
+            bool isCandidate = isSliver(surf, minLen, minQuality, faceI, edgeI);
-            label opposite0 =
-                triSurfaceTools::oppositeVertex
-                (
-                    surf,
-                    faceI,
-                    edgeI
-                );
 
-            pointHit pHit =
+            if (isCandidate)
-                e.line(localPoints).nearestDist
-                (
-                    localPoints[opposite0]
-                );
-
-            if (pHit.hit() && pHit.distance() < minLen)
             {
-                // Remove faceI and split all other faces using this
-                // edge. This is done by 'replacing' the edgeI with the
-                // opposite0 vertex
-                Pout<< "Splitting face " << faceI << " since distance "
-                    << pHit.distance()
-                    << " from vertex " << opposite0
-                    << " to edge " << edgeI
-                    << "  points "
-                    << localPoints[e[0]]
-                    << localPoints[e[1]]
-                    << " is too small" << endl;
-
                 // Mark face as being collapsed
                 if (faceToEdge[faceI] != -1)
                 {
-                    FatalErrorIn("markCollapsedFaces")
+                    FatalErrorIn("markCollapsedFaces(..)")
                         << "Cannot collapse face " << faceI << " since "
                         << " is marked to be collapsed both to edge "
                         << faceToEdge[faceI] << " and " << edgeI
@@ -347,7 +386,7 @@ static void markRegion
 {
     if (faceToEdge[faceI] == -1 || collapseRegion[faceI] != -1)
     {
-        FatalErrorIn("markRegion")
+        FatalErrorIn("markRegion(..)")
             << "Problem : crossed into uncollapsed/regionized face"
             << abort(FatalError);
     }
@@ -383,7 +422,7 @@ static void markRegion
                 }
                 else if (collapseRegion[nbrFaceI] != regionI)
                 {
-                    FatalErrorIn("markRegion")
+                    FatalErrorIn("markRegion(..)")
                         << "Edge:" << edgeI << " between face " << faceI
                         << " with region " << regionI
                         << " and face " << nbrFaceI
@@ -411,8 +450,8 @@ static label markRegions
     {
         if (collapseRegion[faceI] == -1 && faceToEdge[faceI] != -1)
         {
-            Pout<< "markRegions : Marking region:" << regionI
+            //Pout<< "markRegions : Marking region:" << regionI
-                << " starting from face " << faceI << endl;
+            //    << " starting from face " << faceI << endl;
 
             // Collapsed face. Mark connected region with current region number
             markRegion(surf, faceToEdge, regionI++, faceI, collapseRegion);
@@ -728,7 +767,12 @@ static void getSplitVerts
 }
 
 
-label collapseBase(triSurface& surf, const scalar minLen)
+label collapseBase
+(
+    triSurface& surf,
+    const scalar minLen,
+    const scalar minQuality
+)
 {
     label nTotalSplit = 0;
 
@@ -743,7 +787,7 @@ label collapseBase(triSurface& surf, const scalar minLen)
         labelList faceToEdge(surf.size(), -1);
 
         // Calculate faceToEdge (face collapses)
-        markCollapsedFaces(surf, minLen, faceToEdge);
+        markCollapsedFaces(surf, minLen, minQuality, faceToEdge);
 
 
         // Find regions of connected collapsed faces
@@ -754,8 +798,8 @@ label collapseBase(triSurface& surf, const scalar minLen)
 
         label nRegions = markRegions(surf, faceToEdge, collapseRegion);
 
-        Pout<< "Detected " << nRegions << " regions of faces to be collapsed"
+        //Pout<< "Detected " << nRegions << " regions of faces to be collapsed"
-            << nl << endl;
+        //    << nl << endl;
 
         // Pick up all vertices on outside of region
         labelListList outsideVerts
@@ -772,10 +816,10 @@ label collapseBase(triSurface& surf, const scalar minLen)
         {
             spanPoints[regionI] = getSpanPoints(surf, outsideVerts[regionI]);
 
-            Pout<< "For region " << regionI << " found extrema at points "
+            //Pout<< "For region " << regionI << " found extrema at points "
-                << surf.localPoints()[spanPoints[regionI][0]]
+            //    << surf.localPoints()[spanPoints[regionI][0]]
-                << surf.localPoints()[spanPoints[regionI][1]]
+            //    << surf.localPoints()[spanPoints[regionI][1]]
-                << endl;
+            //    << endl;
 
             // Project all non-span points onto the span edge.
             projectNonSpanPoints
@@ -787,21 +831,21 @@ label collapseBase(triSurface& surf, const scalar minLen)
                 orderedWeights[regionI]
             );
 
-            Pout<< "For region:" << regionI
+            //Pout<< "For region:" << regionI
-                << " span:" << spanPoints[regionI]
+            //    << " span:" << spanPoints[regionI]
-                << " orderedVerts:" << orderedVertices[regionI]
+            //    << " orderedVerts:" << orderedVertices[regionI]
-                << " orderedWeights:" << orderedWeights[regionI]
+            //    << " orderedWeights:" << orderedWeights[regionI]
-                << endl;
+            //    << endl;
 
-            writeRegionOBJ
+            //writeRegionOBJ
-            (
+            //(
-                surf,
+            //    surf,
-                regionI,
+            //    regionI,
-                collapseRegion,
+            //    collapseRegion,
-                outsideVerts[regionI]
+            //    outsideVerts[regionI]
-            );
+            //);
 
-            Pout<< endl;
+            //Pout<< endl;
         }
 
 
@@ -864,20 +908,19 @@ label collapseBase(triSurface& surf, const scalar minLen)
                     // Split edge using splitVerts. All non-collapsed triangles
                     // using edge will get split.
 
-
+                    //{
-                    {
+                    //    const pointField& localPoints = surf.localPoints();
-                        const pointField& localPoints = surf.localPoints();
+                    //    Pout<< "edge " << edgeI << ' ' << e
-                        Pout<< "edge " << edgeI << ' ' << e
+                    //        << "  points "
-                            << "  points "
+                    //        << localPoints[e[0]] << ' ' << localPoints[e[1]]
-                            << localPoints[e[0]] << ' ' << localPoints[e[1]]
+                    //        << " split into edges with extra points:"
-                            << " split into edges with extra points:"
+                    //        << endl;
-                            << endl;
+                    //    forAll(splitVerts, i)
-                        forAll(splitVerts, i)
+                    //    {
-                        {
+                    //        Pout<< "    " << splitVerts[i] << " weight "
-                            Pout<< "    " << splitVerts[i] << " weight "
+                    //            << splitWeights[i] << nl;
-                                << splitWeights[i] << nl;
+                    //    }
-                        }
+                    //}
-                    }
 
                     const labelList& eFaces = surf.edgeFaces()[edgeI];
 
@@ -914,7 +957,8 @@ label collapseBase(triSurface& surf, const scalar minLen)
             }
         }
 
-        Pout<< "collapseBase : splitting " << nSplit << " triangles"
+        Info<< "collapseBase : collapsing " << nSplit
+            << " triangles by splitting their base edge."
             << endl;
 
         nTotalSplit += nSplit;
@@ -927,15 +971,15 @@ label collapseBase(triSurface& surf, const scalar minLen)
         // Pack the triangles
         newTris.shrink();
 
-        Pout<< "Resetting surface from " << surf.size() << " to "
+        //Pout<< "Resetting surface from " << surf.size() << " to "
-            << newTris.size() << " triangles" << endl;
+        //    << newTris.size() << " triangles" << endl;
         surf = triSurface(newTris, surf.patches(), surf.localPoints());
 
-        {
+        //{
-            fileName fName("bla" + name(iter) + ".obj");
+        //    fileName fName("bla" + name(iter) + ".obj");
-            Pout<< "Writing surf to " << fName << endl;
+        //    Pout<< "Writing surf to " << fName << endl;
-            surf.write(fName);
+        //    surf.write(fName);
-        }
+        //}
 
         iter++;
     }

applications/utilities/surface/surfaceClean/collapseBase.H

@@ -41,9 +41,14 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//- Keep collapsing all triangles whose height is < minLen.
+//- Keep collapsing all triangles whose height is < minLen or quality < minQ.
 //  Returns number of triangles collapsed.
-label collapseBase(triSurface& surf, const scalar minLen);
+label collapseBase
+(
+    triSurface& surf,
+    const scalar minLen,
+    const scalar minQuality
+);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/surface/surfaceClean/collapseEdge.C

@@ -121,8 +121,8 @@ label collapseEdge(triSurface& surf, const scalar minLen)
                         pointMap[v1] = v;
                         newPoints[v] = 0.5*(localPoints[v1] + localPoints[v]);
 
-                        Pout<< "Collapsing triange " << faceI << " to edge mid "
+                        //Pout<< "Collapsing triange " << faceI
-                            << newPoints[v] << endl;
+                        //    << " to edge mid " << newPoints[v] << endl;
 
                         nCollapsed++;
                         okToCollapse[faceI] = false;
@@ -136,7 +136,8 @@ label collapseEdge(triSurface& surf, const scalar minLen)
             }
         }
 
-        Pout<< "collapseEdge : collapsing " << nCollapsed << " triangles"
+        Info<< "collapseEdge : collapsing " << nCollapsed
+            << " triangles to a single edge."
             << endl;
 
         nTotalCollapsed += nCollapsed;

applications/utilities/surface/surfaceClean/surfaceClean.C

@@ -50,6 +50,7 @@ int main(int argc, char *argv[])
     argList::noParallel();
     argList::validArgs.append("surfaceFile");
     argList::validArgs.append("min length");
+    argList::validArgs.append("min quality");
     argList::validArgs.append("output surfaceFile");
     argList::addBoolOption
     (
@@ -60,10 +61,13 @@ int main(int argc, char *argv[])
 
     const fileName inFileName = args[1];
     const scalar minLen = args.argRead<scalar>(2);
-    const fileName outFileName = args[3];
+    const scalar minQuality = args.argRead<scalar>(3);
+    const fileName outFileName = args[4];
 
     Info<< "Reading surface " << inFileName << nl
-        << "Collapsing all triangles with edges or heights < " << minLen << nl
+        << "Collapsing all triangles with" << nl
+        << "    edges or heights < " << minLen << nl
+        << "    quality          < " << minQuality << nl
         << "Writing result to " << outFileName << nl << endl;
 
 
@@ -90,7 +94,7 @@ int main(int argc, char *argv[])
     }
     while (true)
     {
-        label nSplitEdge = collapseBase(surf, minLen);
+        label nSplitEdge = collapseBase(surf, minLen, minQuality);
 
         if (nSplitEdge == 0)
         {