ENH: surfaceFeatureExtract: Use a dictionary instead of command line options

This commit is contained in:
laurence
2012-03-20 17:46:34 +00:00
parent c6a2710bbe
commit 73ec082aa7
2 changed files with 737 additions and 646 deletions

View File

@ -74,30 +74,6 @@ scalarField calcCurvature(const triSurface& surf)
// The rest of this function adapted from
// CGAL-3.7/examples/Jet_fitting_3/Mesh_estimation.cpp
// Licensed under CGAL-3.7/LICENSE.FREE_USE
// Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
// Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
// Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
// Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
// Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
// (Israel). All rights reserved.
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the
// "Software"), to deal in the Software without restriction, including
// without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit
// persons to whom the Software is furnished to do so, subject to the
// following conditions:
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT
// OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
// THE USE OR OTHER DEALINGS IN THE SOFTWARE.
//Vertex property map, with std::map
typedef std::map<Vertex*, int> Vertex2int_map_type;
@ -503,101 +479,77 @@ int main(int argc, char *argv[])
"extract and write surface features to file"
);
argList::noParallel();
argList::validArgs.append("surface");
argList::validArgs.append("output set");
argList::addOption
(
"includedAngle",
"degrees",
"construct feature set from included angle [0..180]"
"dict",
"word",
"name of dictionary to provide feature extraction information"
);
argList::addOption
(
"set",
"name",
"use existing feature set from file"
);
argList::addOption
(
"minLen",
"scalar",
"remove features shorter than the specified cumulative length"
);
argList::addOption
(
"minElem",
"int",
"remove features with fewer than the specified number of edges"
);
argList::addOption
(
"subsetBox",
"((x0 y0 z0)(x1 y1 z1))",
"remove edges outside specified bounding box"
);
argList::addOption
(
"deleteBox",
"((x0 y0 z0)(x1 y1 z1))",
"remove edges within specified bounding box"
);
argList::addBoolOption
(
"writeObj",
"write extendedFeatureEdgeMesh obj files"
);
argList::addBoolOption
(
"writeVTK",
"write extendedFeatureEdgeMesh vtk files"
);
argList::addBoolOption
(
"closeness",
"span to look for surface closeness"
);
argList::addOption
(
"featureProximity",
"scalar",
"distance to look for close features"
);
argList::addBoolOption
(
"writeVTK",
"write surface property VTK files"
);
argList::addBoolOption
(
"manifoldEdgesOnly",
"remove any non-manifold (open or more than two connected faces) edges"
);
argList::addOption
(
"plane",
"(nx ny nz)(z0 y0 z0)",
"use a plane to create feature edges for 2D meshing"
);
# ifdef ENABLE_CURVATURE
argList::addBoolOption
(
"curvature",
"calculate curvature and closeness fields"
);
# endif
# include "setRootCase.H"
# include "createTime.H"
bool writeVTK = args.optionFound("writeVTK");
word dictName
(
args.optionLookupOrDefault<word>("dict", "surfaceFeatureExtractDict")
);
bool writeObj = args.optionFound("writeObj");
Info<< "Reading " << dictName << nl << endl;
bool curvature = args.optionFound("curvature");
IOdictionary dict
(
IOobject
(
dictName,
runTime.system(),
runTime,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
);
forAllConstIter(dictionary, dict, iter)
{
dictionary surfaceDict = dict.subDict(iter().keyword());
const fileName surfFileName = iter().keyword();
const fileName sFeatFileName = surfFileName.lessExt().name();
Info<< "Surface : " << surfFileName << nl << endl;
const Switch writeVTK =
surfaceDict.lookupOrAddDefault<Switch>("writeVTK", "off");
const Switch writeObj =
surfaceDict.lookupOrAddDefault<Switch>("writeObj", "off");
const Switch writeFeatureEdgeMesh =
surfaceDict.lookupOrAddDefault<Switch>
(
"writeFeatureEdgeMesh",
"off"
);
const Switch curvature =
surfaceDict.lookupOrAddDefault<Switch>("curvature", "off");
const Switch featureProximity =
surfaceDict.lookupOrAddDefault<Switch>("featureProximity", "off");
const Switch closeness =
surfaceDict.lookupOrAddDefault<Switch>("closeness", "off");
const word extractionMethod = surfaceDict.lookup("extractionMethod");
#ifndef ENABLE_CURVATURE
if (curvature)
{
WarningIn(args.executable())
<< "Curvature calculation has been requested but "
<< args.executable() << " has not " << nl
<< " been compiled with CGAL. "
<< "Skipping the curvature calculation." << endl;
}
#else
if (curvature && env("FOAM_SIGFPE"))
{
WarningIn(args.executable())
@ -607,25 +559,16 @@ int main(int argc, char *argv[])
<< " (infinite curvature)" << nl
<< " Switch it off in case of problems." << endl;
}
#endif
Info<< "Feature line extraction is only valid on closed manifold surfaces."
<< endl;
const fileName surfFileName = args[1];
const fileName outFileName = args[2];
Info<< "Surface : " << surfFileName << nl
<< "Output feature set : " << outFileName << nl
<< endl;
fileName sFeatFileName = surfFileName.lessExt().name();
Info<< nl << "Feature line extraction is only valid on closed manifold "
<< "surfaces." << endl;
// Read
// ~~~~
triSurface surf(surfFileName);
triSurface surf("constant/triSurface/" + surfFileName);
Info<< "Statistics:" << endl;
surf.writeStats(Info);
@ -638,38 +581,52 @@ int main(int argc, char *argv[])
faces[fI] = surf[fI].triFaceFace();
}
// Either construct features from surface&featureangle or read set.
// Either construct features from surface & featureAngle or read set.
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
surfaceFeatures set(surf);
if (args.optionFound("set"))
if (extractionMethod == "extractFromFile")
{
const fileName setName = args["set"];
const fileName featureEdgeFile =
surfaceDict.subDict("extractFromFile").lookup
(
"featureEdgeFile"
);
Info<< "Reading existing feature set from file " << setName << endl;
edgeMesh eMesh(featureEdgeFile);
set = surfaceFeatures(surf, setName);
// Sometimes duplicate edges are present. Remove them.
eMesh.mergeEdges();
Info<< nl << "Reading existing feature edges from file "
<< featureEdgeFile << endl;
set = surfaceFeatures(surf, eMesh.points(), eMesh.edges());
}
else if (args.optionFound("includedAngle"))
else if (extractionMethod == "extractFromSurface")
{
const scalar includedAngle = args.optionRead<scalar>("includedAngle");
const scalar includedAngle =
readScalar
(
surfaceDict.subDict("extractFromSurface").lookup
(
"includedAngle"
)
);
Info<< "Constructing feature set from included angle " << includedAngle
<< endl;
Info<< nl << "Constructing feature set from included angle "
<< includedAngle << endl;
set = surfaceFeatures(surf, includedAngle);
// Info<< nl << "Writing initial features" << endl;
// set.write("initial.fSet");
// set.writeObj("initial");
}
else
{
FatalErrorIn(args.executable())
<< "No initial feature set. Provide either one"
<< " of -set (to read existing set)" << nl
<< " or -includedAngle (to new set construct from angle)"
<< " of extractFromFile (to read existing set)" << nl
<< " or extractFromSurface (to construct new set from angle)"
<< exit(FatalError);
}
@ -687,23 +644,25 @@ int main(int argc, char *argv[])
// Trim set
// ~~~~~~~~
scalar minLen = -GREAT;
if (args.optionReadIfPresent("minLen", minLen))
if (surfaceDict.isDict("trimFeatures"))
{
Info<< "Removing features of length < " << minLen << endl;
}
dictionary trimDict = surfaceDict.subDict("trimFeatures");
label minElem = 0;
if (args.optionReadIfPresent("minElem", minElem))
{
Info<< "Removing features with number of edges < " << minElem << endl;
}
scalar minLen =
trimDict.lookupOrAddDefault<scalar>("minLen", -GREAT);
label minElem = trimDict.lookupOrAddDefault<label>("minElem", 0);
// Trim away small groups of features
if (minElem > 0 || minLen > 0)
{
Info<< "Removing features of length < "
<< minLen << endl;
Info<< "Removing features with number of edges < "
<< minElem << endl;
set.trimFeatures(minLen, minElem);
Info<< endl << "Removed small features" << endl;
}
}
@ -713,12 +672,13 @@ int main(int argc, char *argv[])
// Convert to marked edges, points
List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
if (args.optionFound("subsetBox"))
if (surfaceDict.isDict("subsetFeatures"))
{
treeBoundBox bb
(
args.optionLookup("subsetBox")()
);
dictionary subsetDict = surfaceDict.subDict("subsetFeatures");
if (subsetDict.found("insideBox"))
{
treeBoundBox bb(subsetDict.lookup("insideBox")());
Info<< "Removing all edges outside bb " << bb << endl;
dumpBox(bb, "subsetBox.obj");
@ -731,12 +691,9 @@ int main(int argc, char *argv[])
edgeStat
);
}
else if (args.optionFound("deleteBox"))
else if (subsetDict.found("outsideBox"))
{
treeBoundBox bb
(
args.optionLookup("deleteBox")()
);
treeBoundBox bb(subsetDict.lookup("outsideBox")());
Info<< "Removing all edges inside bb " << bb << endl;
dumpBox(bb, "deleteBox.obj");
@ -750,7 +707,10 @@ int main(int argc, char *argv[])
);
}
if (args.optionFound("manifoldEdgesOnly"))
const Switch manifoldEdges =
subsetDict.lookupOrAddDefault<Switch>("manifoldEdges", "no");
if (manifoldEdges)
{
Info<< "Removing all non-manifold edges" << endl;
@ -763,9 +723,9 @@ int main(int argc, char *argv[])
}
}
if (args.optionFound("plane"))
if (subsetDict.found("plane"))
{
plane cutPlane(args.optionLookup("plane")());
plane cutPlane(subsetDict.lookup("plane")());
deleteEdges
(
@ -775,22 +735,23 @@ int main(int argc, char *argv[])
);
Info<< "Only edges that intersect the plane with normal "
<< cutPlane.normal() << " and base point " << cutPlane.refPoint()
<< cutPlane.normal()
<< " and base point " << cutPlane.refPoint()
<< " will be included as feature edges."<< endl;
}
}
surfaceFeatures newSet(surf);
newSet.setFromStatus(edgeStat);
//Info<< endl << "Writing trimmed features to " << outFileName << endl;
//newSet.write(outFileName);
// Info<< endl << "Writing edge objs." << endl;
// newSet.writeObj("final");
if (writeObj)
{
newSet.writeObj("final");
}
Info<< nl
<< "Final feature set:" << nl
<< "Final feature set after trimming and subsetting:" << nl
<< " feature points : " << newSet.featurePoints().size() << nl
<< " feature edges : " << newSet.featureEdges().size() << nl
<< " of which" << nl
@ -807,8 +768,8 @@ int main(int argc, char *argv[])
sFeatFileName + ".extendedFeatureEdgeMesh"
);
Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
<< endl;
Info<< nl << "Writing extendedFeatureEdgeMesh to "
<< feMesh.objectPath() << endl;
if (writeObj)
{
@ -818,6 +779,7 @@ int main(int argc, char *argv[])
feMesh.write();
// Write a featureEdgeMesh for backwards compatibility
if (writeFeatureEdgeMesh)
{
featureEdgeMesh bfeMesh
(
@ -922,10 +884,10 @@ int main(int argc, char *argv[])
// )
// );
if (args.optionFound("closeness"))
if (closeness)
{
Info<< nl << "Extracting internal and external closeness of surface."
<< endl;
Info<< nl << "Extracting internal and external closeness of "
<< "surface." << endl;
// Internal and external closeness
@ -938,18 +900,21 @@ int main(int argc, char *argv[])
//args.optionReadIfPresent("closeness", span);
scalar externalAngleTolerance = 10;
scalar externalToleranceCosAngle = Foam::cos
scalar externalToleranceCosAngle =
Foam::cos
(
degToRad(180 - externalAngleTolerance)
);
scalar internalAngleTolerance = 45;
scalar internalToleranceCosAngle = Foam::cos
scalar internalToleranceCosAngle =
Foam::cos
(
degToRad(180 - internalAngleTolerance)
);
Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle << nl
Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle
<< nl
<< "internalToleranceCosAngle: " << internalToleranceCosAngle
<< endl;
@ -989,7 +954,15 @@ int main(int argc, char *argv[])
}
else if (hitInfo[0].index() != fI)
{
drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
drawHitProblem
(
fI,
surf,
start,
faceCentres,
end,
hitInfo
);
// FatalErrorIn(args.executable())
// << "findLineAll did not hit its own face."
@ -1014,7 +987,15 @@ int main(int argc, char *argv[])
if (ownHitI < 0)
{
drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
drawHitProblem
(
fI,
surf,
start,
faceCentres,
end,
hitInfo
);
// FatalErrorIn(args.executable())
// << "findLineAll did not hit its own face."
@ -1022,35 +1003,45 @@ int main(int argc, char *argv[])
}
else if (ownHitI == 0)
{
// There are no internal hits, the first hit is the closest
// external hit
// There are no internal hits, the first hit is the
// closest external hit
if
(
(normals[fI] & normals[hitInfo[ownHitI + 1].index()])
(
normals[fI]
& normals[hitInfo[ownHitI + 1].index()]
)
< externalToleranceCosAngle
)
{
externalCloseness[fI] = mag
externalCloseness[fI] =
mag
(
faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
faceCentres[fI]
- hitInfo[ownHitI + 1].hitPoint()
);
}
}
else if (ownHitI == hitInfo.size() - 1)
{
// There are no external hits, the last but one hit is the
// closest internal hit
// There are no external hits, the last but one hit is
// the closest internal hit
if
(
(normals[fI] & normals[hitInfo[ownHitI - 1].index()])
(
normals[fI]
& normals[hitInfo[ownHitI - 1].index()]
)
< internalToleranceCosAngle
)
{
internalCloseness[fI] = mag
internalCloseness[fI] =
mag
(
faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
faceCentres[fI]
- hitInfo[ownHitI - 1].hitPoint()
);
}
}
@ -1058,25 +1049,35 @@ int main(int argc, char *argv[])
{
if
(
(normals[fI] & normals[hitInfo[ownHitI + 1].index()])
(
normals[fI]
& normals[hitInfo[ownHitI + 1].index()]
)
< externalToleranceCosAngle
)
{
externalCloseness[fI] = mag
externalCloseness[fI] =
mag
(
faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
faceCentres[fI]
- hitInfo[ownHitI + 1].hitPoint()
);
}
if
(
(normals[fI] & normals[hitInfo[ownHitI - 1].index()])
(
normals[fI]
& normals[hitInfo[ownHitI - 1].index()]
)
< internalToleranceCosAngle
)
{
internalCloseness[fI] = mag
internalCloseness[fI] =
mag
(
faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
faceCentres[fI]
- hitInfo[ownHitI - 1].hitPoint()
);
}
}
@ -1149,9 +1150,10 @@ int main(int argc, char *argv[])
#ifdef ENABLE_CURVATURE
if (args.optionFound("curvature"))
if (curvature)
{
Info<< nl << "Extracting curvature of surface at the points." << endl;
Info<< nl << "Extracting curvature of surface at the points."
<< endl;
scalarField k = calcCurvature(surf);
@ -1201,13 +1203,13 @@ int main(int argc, char *argv[])
#endif
if (args.optionFound("featureProximity"))
if (featureProximity)
{
Info<< nl << "Extracting proximity of close feature points and edges "
<< "to the surface" << endl;
Info<< nl << "Extracting proximity of close feature points and "
<< "edges to the surface" << endl;
const scalar searchDistance =
args.optionRead<scalar>("featureProximity");
readScalar(surfaceDict.lookup("maxFeatureProximity"));
const scalar radiusSqr = sqr(searchDistance);
@ -1274,6 +1276,9 @@ int main(int argc, char *argv[])
}
}
Info<< endl;
}
Info<< "End\n" << endl;
return 0;

View File

@ -0,0 +1,86 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object surfaceFeatureExtractDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
surface1.stl
{
// extractFromFile || extractFromSurface
extractionMethod extractFromFile;
extractFromFile
{
// Load from an existing feature edge file
featureEdgeFile "constant/triSurface/featureEdges.nas";
}
trimFeatures
{
// Remove features with fewer than the specified number of edges
minElem 0;
// Remove features shorter than the specified cumulative length
minLen 0.0;
}
subsetFeatures
{
// Use a plane to select feature edges
// (normal)(basePoint)
plane (1 0 0)(0 0 0);
// Select feature edges using a box
// (minPt)(maxPt)
insideBox (0 0 0)(1 1 1);
outsideBox (0 0 0)(1 1 1);
// Remove any non-manifold (open or > 2 connected faces) edges
manifoldEdges no;
}
// Output the curvature of the surface
curvature no;
// Output the proximity of feature points and edges to each other
featureProximity no;
// The maximum search distance to use when looking for other feature
// points and edges
maxFeatureProximity 1;
// Out put the closeness of surface elements to other surface elements.
closeness no;
// Write options
writeVTK no;
writeObj yes;
writeFeatureEdgeMesh no;
}
surface2.nas
{
extractionMethod extractFromSurface;
extractFromSurface
{
// Mark edges whose adjacent surface normals are at an angle less
// than includedAngle as features
// - 0 : selects no edges
// - 180: selects all edges
includedAngle 120;
}
}
// ************************************************************************* //