ENH: readSTLASCII: single precision point merging

src/OpenFOAM/meshes/meshTools/mergePoints.C

@@ -23,25 +23,25 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "mergePoints.H"
-#include "SortableList.H"
 #include "ListOps.H"
+#include "point.H"
+#include "Field.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::mergePoints
+template<class Type>
+Foam::label Foam::mergePoints
 (
-    const UList<point>& points,
+    const UList<Type>& points,
     const scalar mergeTol,
     const bool verbose,
     labelList& pointMap,
-    List<point>& newPoints,
-    const point& origin
+    const Type& origin
 )
 {
-    point compareOrigin = origin;
+    Type compareOrigin = origin;
 
-    if (origin == point(VGREAT, VGREAT, VGREAT))
+    if (origin == Type::max)
     {
         if (points.size())
         {
@@ -53,12 +53,9 @@ bool Foam::mergePoints
     pointMap.setSize(points.size());
     pointMap = -1;
 
-    // Storage for merged points
-    newPoints.setSize(points.size());
-
     if (points.empty())
     {
-        return false;
+        return points.size();
     }
 
     // We're comparing distance squared to origin first.
@@ -70,33 +67,56 @@ bool Foam::mergePoints
     //     x^2+y^2+z^2 + 2*mergeTol*(x+z+y) + mergeTol^2*...
     // so the difference will be 2*mergeTol*(x+y+z)
 
-    const scalar mergeTolSqr = sqr(mergeTol);
+    const scalar mergeTolSqr = Foam::sqr(scalar(mergeTol));
 
     // Sort points by magSqr
-    const pointField d(points - compareOrigin);
-    SortableList<scalar> sortedMagSqr(magSqr(d));
-    scalarField sortedTol(points.size());
-    forAll(sortedMagSqr.indices(), sortI)
+    const Field<Type> d(points - compareOrigin);
+
+    List<scalar> magSqrD(d.size());
+    forAll(d, pointI)
     {
-        const point& pt = d[sortedMagSqr.indices()[sortI]];
-        sortedTol[sortI] = 2*mergeTol*(mag(pt.x())+mag(pt.y())+mag(pt.z()));
+        magSqrD[pointI] = magSqr(d[pointI]);
     }
+    labelList order;
+    sortedOrder(magSqrD, order);
+
+
+    Field<scalar> sortedTol(points.size());
+    forAll(order, sortI)
+    {
+        label pointI = order[sortI];
 
-    bool hasMerged = false;
+        // Convert to scalar precision
+        const point pt
+        (
+            scalar(d[pointI].x()),
+            scalar(d[pointI].y()),
+            scalar(d[pointI].z())
+        );
+        sortedTol[sortI] = 2*mergeTol*(mag(pt.x())+mag(pt.y())+mag(pt.z()));
+    }
 
     label newPointI = 0;
 
 
     // Handle 0th point separately (is always unique)
-    label pointI = sortedMagSqr.indices()[0];
-    pointMap[pointI] = newPointI;
-    newPoints[newPointI++] = points[pointI];
+    label pointI = order[0];
+    pointMap[pointI] = newPointI++;
 
 
-    for (label sortI = 1; sortI < sortedMagSqr.size(); sortI++)
+    for (label sortI = 1; sortI < order.size(); sortI++)
     {
         // Get original point index
-        label pointI = sortedMagSqr.indices()[sortI];
+        label pointI = order[sortI];
+        const scalar mag2 = magSqrD[order[sortI]];
+        // Convert to scalar precision
+        const point pt
+        (
+            scalar(points[pointI].x()),
+            scalar(points[pointI].y()),
+            scalar(points[pointI].z())
+        );
+
 
         // Compare to previous points to find equal one.
         label equalPointI = -1;
@@ -105,17 +125,19 @@ bool Foam::mergePoints
         (
             label prevSortI = sortI - 1;
             prevSortI >= 0
-         && mag
-            (
-                sortedMagSqr[prevSortI]
-             -  sortedMagSqr[sortI]
-            ) <= sortedTol[sortI];
+         && (mag(magSqrD[order[prevSortI]] - mag2) <= sortedTol[sortI]);
             prevSortI--
         )
         {
-            label prevPointI = sortedMagSqr.indices()[prevSortI];
+            label prevPointI = order[prevSortI];
+            const point prevPt
+            (
+                scalar(points[prevPointI].x()),
+                scalar(points[prevPointI].y()),
+                scalar(points[prevPointI].z())
+            );
 
-            if (magSqr(points[pointI] - points[prevPointI]) <= mergeTolSqr)
+            if (magSqr(pt - prevPt) <= mergeTolSqr)
             {
                 // Found match.
                 equalPointI = prevPointI;
@@ -130,8 +152,6 @@ bool Foam::mergePoints
             // Same coordinate as equalPointI. Map to same new point.
             pointMap[pointI] = pointMap[equalPointI];
 
-            hasMerged = true;
-
             if (verbose)
             {
                 Pout<< "Foam::mergePoints : Merging points "
@@ -144,14 +164,41 @@ bool Foam::mergePoints
         else
         {
             // Differs. Store new point.
-            pointMap[pointI] = newPointI;
-            newPoints[newPointI++] = points[pointI];
+            pointMap[pointI] = newPointI++;
         }
     }
 
-    newPoints.setSize(newPointI);
+    return newPointI;
+}
+
+
+template<class Type>
+bool Foam::mergePoints
+(
+    const UList<Type>& points,
+    const scalar mergeTol,
+    const bool verbose,
+    labelList& pointMap,
+    List<Type>& newPoints,
+    const Type& origin = Type::zero
+)
+{
+    label nUnique = mergePoints
+    (
+        points,
+        mergeTol,
+        verbose,
+        pointMap,
+        origin
+    );
+
+    newPoints.setSize(nUnique);
+    forAll(pointMap, pointI)
+    {
+        newPoints[pointMap[pointI]] = points[pointI];
+    }
 
-    return hasMerged;
+    return (nUnique != points.size());
 }
 
 

src/OpenFOAM/meshes/meshTools/mergePoints.H

@@ -32,8 +32,8 @@ SourceFiles
 #ifndef mergePoints_H
 #define mergePoints_H
 
-#include "scalarField.H"
-#include "pointField.H"
+#include "scalar.H"
+#include "labelList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,23 +44,41 @@ namespace Foam
                            Function mergePoints Declaration
 \*---------------------------------------------------------------------------*/
 
-//- Sort & merge points. All points closer than/equal mergeTol get merged.
-//  Outputs the new unique points and a map from old to new. Returns
-//  true if anything merged, false otherwise.
+//- Sorts and merges points. All points closer than/equal mergeTol get merged.
+//  Returns the number of unique points and a map from old to new.
+template<class Type>
+label mergePoints
+(
+    const UList<Type>& points,
+    const scalar mergeTol,
+    const bool verbose,
+    labelList& pointMap,
+    const Type& origin = Type::zero
+);
+
+//- Sorts and merges points. Determines new points. Returns true if anything
+//  merged (though newPoints still sorted even if not merged).
+template<class Type>
 bool mergePoints
 (
-    const UList<point>& points,
+    const UList<Type>& points,
     const scalar mergeTol,
     const bool verbose,
     labelList& pointMap,
-    List<point>& newPoints,
-    const point& origin = point::zero
+    List<Type>& newPoints,
+    const Type& origin = Type::zero
 );
 
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "mergePoints.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/primitives/Vector/floatVector/floatVector.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Vector of floats.
+
+\*---------------------------------------------------------------------------*/
+
+#include "floatVector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<>
+const char* const floatVector::typeName = "floatVector";
+
+template<>
+const char* floatVector::componentNames[] = {"x", "y", "z"};
+
+template<>
+const floatVector floatVector::zero(0, 0, 0);
+
+template<>
+const floatVector floatVector::one(1, 1, 1);
+
+template<>
+const floatVector floatVector::max
+(
+    floatScalarVGREAT,
+    floatScalarVGREAT,
+    floatScalarVGREAT
+);
+
+template<>
+const floatVector floatVector::min
+(
+    -floatScalarVGREAT,
+    -floatScalarVGREAT,
+    -floatScalarVGREAT
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/primitives/Vector/floatVector/floatVector.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::floatVector
+
+Description
+    A float version of vector
+
+SourceFiles
+    floatVector.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef floatVector_H
+#define floatVector_H
+
+#include "Vector.H"
+#include "contiguous.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef Vector<float> floatVector;
+
+
+//- Data associated with floatVector type are contiguous
+template<>
+inline bool contiguous<floatVector>() {return true;}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/triSurface/triSurface/interfaces/STL/readSTLASCII.L

@@ -33,8 +33,10 @@ License
 
 #include "IFstream.H"
 #include "triSurface.H"
-#include "STLpoint.H"
+#include "floatVector.H"
 #include "OSspecific.H"
+#include "mergePoints.H"
+//#include "memInfo.H"
 
 using namespace Foam;
 
@@ -77,8 +79,8 @@ class STLLexer
         label lineNo_;
         word startError_;
 
-        DynamicList<STLpoint> STLpoints_;
-        //DynamicList<STLpoint> STLnormals_;
+        DynamicList<floatVector> STLpoints_;
+        //DynamicList<floatVector > STLnormals_;
         DynamicList<label> STLlabels_;
         HashTable<label, word> STLsolidNames_;
 
@@ -103,12 +105,12 @@ public:
             return nTriangles_;
         }
 
-        DynamicList<STLpoint>& STLpoints()
+        DynamicList<floatVector>& STLpoints()
         {
             return STLpoints_;
         }
 
-        //DynamicList<STLpoint>& STLnormals()
+        //DynamicList<floatVector>& STLnormals()
         //{
         //    return STLnormals_;
         //}
@@ -202,8 +204,8 @@ endsolid              {space}("endsolid"|"ENDSOLID")({some_space}{word})*
     // End of read character pointer returned by strtof
     // char* endPtr;
 
-    STLpoint normal;
-    STLpoint vertex;
+    floatVector normal;
+    floatVector vertex;
     label cmpt = 0;   // component index used for reading vertex
 
     static const char* stateNames[7] =
@@ -387,16 +389,25 @@ bool triSurface::readSTLASCII(const fileName& STLfileName)
             << exit(FatalError);
     }
 
+    //memInfo memStat;
+
+    //memStat.update();
+    //Pout<< "At start:" << memStat.rss() << endl;
+
     // Create the lexer obtaining the approximate number of vertices in the STL
     // from the file size
     STLLexer lexer(&STLstream.stdStream(), Foam::fileSize(STLfileName)/400);
     while (lexer.lex() != 0)
     {}
 
-    DynamicList<STLpoint>& STLpoints = lexer.STLpoints();
+    //memStat.update();
+    //Pout<< "After lexing:" << memStat.rss() << endl;
+
+    DynamicList<floatVector>& STLpoints = lexer.STLpoints();
+    DynamicList<label>& STLlabels = lexer.STLlabels();
 
     /*
-    DynamicList<STLpoint>& STLnormals = lexer.STLnormals();
+    DynamicList<floatVector>& STLnormals = lexer.STLnormals();
 
     if (STLpoints.size() != 3*STLnormals.size())
     {
@@ -410,35 +421,51 @@ bool triSurface::readSTLASCII(const fileName& STLfileName)
     }
     */
 
-    pointField rawPoints(STLpoints.size());
+    labelList pointMap;
+    label nUniquePoints = mergePoints
+    (
+        STLpoints,
+        100*SMALL,              // merge distance
+        false,                  // verbose
+        pointMap
+    );
 
-    forAll(rawPoints, i)
-    {
-        rawPoints[i] = STLpoints[i];
-    }
+    //memStat.update();
+    //Pout<< "After merging:" << memStat.rss() << endl;
 
-    STLpoints.clear();
+
+    pointField& sp = storedPoints();
 
     setSize(lexer.nTriangles());
-    DynamicList<label>& STLlabels = lexer.STLlabels();
+    sp.setSize(nUniquePoints);
+    forAll(STLpoints, pointI)
+    {
+        const floatVector& pt = STLpoints[pointI];
+        sp[pointMap[pointI]] = vector
+        (
+            scalar(pt.x()),
+            scalar(pt.y()),
+            scalar(pt.z())
+        );
+    }
 
     // Assign triangles
-    label pointi = 0;
+    label pointI = 0;
     forAll(*this, i)
     {
-        operator[](i)[0] = pointi++;
-        operator[](i)[1] = pointi++;
-        operator[](i)[2] = pointi++;
+        operator[](i)[0] = pointMap[pointI++];
+        operator[](i)[1] = pointMap[pointI++];
+        operator[](i)[2] = pointMap[pointI++];
         operator[](i).region() = STLlabels[i];
     }
 
-    STLlabels.clear();
+    //memStat.update();
+    //Pout<< "After assigning:" << memStat.rss() << endl;
 
-    // Assign coordinates
-    storedPoints().transfer(rawPoints);
 
-    // Stitch all points within SMALL meters.
-    stitchTriangles();
+    STLpoints.clear();
+    STLlabels.clear();
+
 
     // Convert solidNames into regionNames
     patches_.setSize(lexer.STLsolidNames().size());
@@ -457,6 +484,9 @@ bool triSurface::readSTLASCII(const fileName& STLfileName)
     // Fill in the missing information in the patches
     setDefaultPatches();
 
+    //memStat.update();
+    //Pout<< "After patchifying:" << memStat.rss() << endl;
+
     return true;
 }
 

src/triSurface/triSurface/interfaces/STL/readSTLBINARY.C

@@ -28,7 +28,8 @@ License
 #include "IFstream.H"
 #include "OSspecific.H"
 #include "gzstream.h"
-
+#include "floatVector.H"
+#include "mergePoints.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -94,40 +95,55 @@ bool Foam::triSurface::readSTLBINARY(const fileName& STLfileName)
     // Everything OK so go ahead and read the triangles.
 
     // Allocate storage for raw points
-    pointField rawPoints(3*nTris);
-
-    // Allocate storage for triangles
+    List<floatVector> STLpoints(3*nTris);
     setSize(nTris);
 
-    label rawPointI = 0;
+    label pointI = 0;
 
-    // Read the triangles
-    forAll(*this, i)
+    for (label i = 0; i < nTris; i++)
     {
         // Read an STL triangle
         STLtriangle stlTri(STLfile);
 
-        // Set the rawPoints to the vertices of the STL triangle
-        // and set the point labels of the labelledTri
-        rawPoints[rawPointI] = stlTri.a();
-        operator[](i)[0] = rawPointI++;
-
-        rawPoints[rawPointI] = stlTri.b();
-        operator[](i)[1] = rawPointI++;
-
-        rawPoints[rawPointI] = stlTri.c();
-        operator[](i)[2] = rawPointI++;
-
+        // Set the STLpoints to the vertices of the STL triangle
+        STLpoints[pointI++] = stlTri.a();
+        STLpoints[pointI++] = stlTri.b();
+        STLpoints[pointI++] = stlTri.c();
         operator[](i).region() = stlTri.region();
     }
 
-    //STLfile.close();
+    // Stitch points
+    labelList pointMap;
+    label nUniquePoints = mergePoints
+    (
+        STLpoints,
+        10*SMALL,               // merge distance
+        false,                  // verbose
+        pointMap                // old to new
+    );
+
+    pointField& sp = storedPoints();
 
-    // Assign coordinates
-    storedPoints().transfer(rawPoints);
+    sp.setSize(nUniquePoints);
+    forAll(STLpoints, pointI)
+    {
+        const floatVector& pt = STLpoints[pointI];
+        sp[pointMap[pointI]] = vector
+        (
+            scalar(pt.x()),
+            scalar(pt.y()),
+            scalar(pt.z())
+        );
+    }
 
-    // Stitch all points within SMALL meters.
-    stitchTriangles();
+    // Assign triangles
+    pointI = 0;
+    forAll(*this, i)
+    {
+        operator[](i)[0] = pointMap[pointI++];
+        operator[](i)[1] = pointMap[pointI++];
+        operator[](i)[2] = pointMap[pointI++];
+    }
 
     return true;
 }