Geometry tools update finalizing with Unit tests

This commit is contained in:
Jacob Støren 2013-12-06 14:23:51 +01:00
parent 776ed7755e
commit d80611a61c
4 changed files with 312 additions and 70 deletions

View File

@ -29,6 +29,7 @@
#include "cvfGeometryTools.h"
#include "cvfBoundingBoxTree.h"
#include <array>
using namespace cvf;
@ -86,7 +87,7 @@ void ControlVolume::calculateCubeFaceStatus(const cvf::Vec3dArray& nodeCoords, d
template <typename NodeArrayType, typename NodeType, typename IndexType>
NodeType quadNormal (const ArrayWrapperConst<NodeArrayType, NodeType>& nodeCoords,
NodeType quadNormal (ArrayWrapperConst<NodeArrayType, NodeType> nodeCoords,
const IndexType cubeFaceIndices[4] )
{
return ( nodeCoords[cubeFaceIndices[2]] - nodeCoords[cubeFaceIndices[0]]) ^
@ -143,7 +144,9 @@ public:
private:
QuadFaceIntersectorImplHandle * m_implementation;
};
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> createVertices()
{
std::vector<cvf::Vec3d> vxs;
@ -166,15 +169,221 @@ std::vector<cvf::Vec3d> createVertices()
return vxs;
}
std::vector<caf::FixedArray<cvf::uint, 4> > getCubeFaces()
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<caf::UintArray4> getCubeFaces()
{
std::vector<caf::UintArray4 > cubeFaces;
cvf::uint faces[4*4] = {
0, 1, 2, 3,
4, 5, 6, 7,
5, 8, 9, 6,
10, 11, 12, 13
};
cubeFaces.resize(4);
cubeFaces[0] = &faces[0];
cubeFaces[1] = &faces[4];
cubeFaces[2] = &faces[8];
cubeFaces[3] = &faces[12];
cubeFaces[3] = cubeFaces[2];
return cubeFaces;
}
std::ostream& operator<< (std::ostream& stream, std::vector<uint> v)
{
for (size_t i = 0; i < v.size(); ++i)
{
stream << v[i] << " ";
}
return stream;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST(CellFaceIntersectionTst, Intersection1)
{
std::vector<cvf::Vec3d> nodes = createVertices();
std::vector<cvf::Vec3d> additionalVertices;
std::vector< std::vector<uint> > overlapPolygons;
std::vector<caf::UintArray4> faces = getCubeFaces();
EdgeIntersectStorage<uint> edgeIntersectionStorage;
edgeIntersectionStorage.setVertexCount(nodes.size());
{
std::vector<uint> polygon;
bool isOk = false;
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(
&polygon,
&additionalVertices,
edgeIntersectionStorage,
wrapArrayConst(&nodes),
faces[0].data(),
faces[1].data(),
1e-6);
EXPECT_EQ( 5, polygon.size());
EXPECT_EQ( (size_t)2, additionalVertices.size());
EXPECT_TRUE(isOk);
overlapPolygons.push_back(polygon);
std::cout << polygon << std::endl;
}
{
std::vector<uint> polygon;
bool isOk = false;
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(
&polygon,
&additionalVertices,
edgeIntersectionStorage,
wrapArrayConst(&nodes),
faces[0].data(),
faces[2].data(),
1e-6);
EXPECT_EQ( 5, polygon.size());
EXPECT_EQ( (size_t)4, additionalVertices.size());
EXPECT_TRUE(isOk);
overlapPolygons.push_back(polygon);
std::cout << polygon << std::endl;
}
{
std::vector<uint> polygon;
bool isOk = false;
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(
&polygon,
&additionalVertices,
edgeIntersectionStorage,
wrapArrayConst(&nodes),
faces[0].data(),
faces[3].data(),
1e-6);
EXPECT_EQ( 3, polygon.size());
EXPECT_EQ( (size_t)6, additionalVertices.size());
EXPECT_TRUE(isOk);
overlapPolygons.push_back(polygon);
std::cout << polygon << std::endl;
}
nodes.insert(nodes.end(), additionalVertices.begin(), additionalVertices.end());
std::vector<uint> basePolygon;
basePolygon.insert(basePolygon.begin(), faces[0].data(), &(faces[0].data()[4]));
for (uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx)
{
bool inserted = GeometryTools::insertVertexInPolygon(
&basePolygon,
wrapArrayConst(&nodes),
vxIdx,
1e-6
);
}
EXPECT_EQ( 8, basePolygon.size());
std::cout << "Bp: " << basePolygon << std::endl;
for (size_t pIdx = 0; pIdx < overlapPolygons.size(); ++pIdx)
{
for (uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx)
{
bool inserted = GeometryTools::insertVertexInPolygon(
&overlapPolygons[pIdx],
wrapArrayConst(&nodes),
vxIdx,
1e-6
);
}
if (pIdx == 0)
{
EXPECT_EQ(5, overlapPolygons[pIdx].size());
}
if (pIdx == 1)
{
EXPECT_EQ(5, overlapPolygons[pIdx].size());
}
if (pIdx == 2)
{
EXPECT_EQ(4, overlapPolygons[pIdx].size());
}
std::cout << "Op" << pIdx << ":" << overlapPolygons[pIdx] << std::endl;
}
Vec3d normal = quadNormal(wrapArrayConst(&nodes), faces[0].data());
std::vector<bool> faceOverlapPolygonWindingSameAsCubeFaceFlags;
faceOverlapPolygonWindingSameAsCubeFaceFlags.resize(overlapPolygons.size(), true);
{
std::vector<uint> freeFacePolygon;
bool hasHoles = false;
std::vector< std::vector<uint>* > overlapPolygonPtrs;
for (size_t pIdx = 0; pIdx < overlapPolygons.size(); ++pIdx)
{
overlapPolygonPtrs.push_back(&(overlapPolygons[pIdx]));
}
GeometryTools::calculatePartiallyFreeCubeFacePolygon(
wrapArrayConst(&nodes),
wrapArrayConst(&basePolygon),
normal,
overlapPolygonPtrs,
faceOverlapPolygonWindingSameAsCubeFaceFlags,
&freeFacePolygon,
&hasHoles
);
EXPECT_EQ( 4, freeFacePolygon.size());
EXPECT_FALSE(hasHoles);
std::cout << "FF1: " << freeFacePolygon << std::endl;
}
{
std::vector<uint> freeFacePolygon;
bool hasHoles = false;
std::vector< std::vector<uint>* > overlapPolygonPtrs;
for (size_t pIdx = 0; pIdx < 1; ++pIdx)
{
overlapPolygonPtrs.push_back(&(overlapPolygons[pIdx]));
}
GeometryTools::calculatePartiallyFreeCubeFacePolygon(
wrapArrayConst(&nodes),
wrapArrayConst(&basePolygon),
normal,
overlapPolygonPtrs,
faceOverlapPolygonWindingSameAsCubeFaceFlags,
&freeFacePolygon,
&hasHoles
);
EXPECT_EQ( 9, freeFacePolygon.size());
EXPECT_FALSE(hasHoles);
std::cout << "FF2: " << freeFacePolygon << std::endl;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST(CellFaceIntersectionTst, Intersection)
{

View File

@ -474,58 +474,6 @@ cvf::Vec3d GeometryTools::barycentricCoords(const cvf::Vec3d& t0, const cvf::Ve
return m;
}
//--------------------------------------------------------------------------------------------------
/// Inserts the vertex into the polygon if it fits along one of the edges within the tolerance.
/// The method returns true if it was inserted, or if it was already in the polygon, or if it was
/// within the tolerance of an existing vertex in the polygon.
/// In the latter situation it replaces the previous vertex in the polygon.
///
/// Todo: If a vertex is replaced, the VxToCv map in TimeStepGeometry should be updated
//--------------------------------------------------------------------------------------------------
bool GeometryTools::insertVertexInPolygon(std::list<std::pair<cvf::uint, bool> >* polygon, const cvf::Vec3dArray& nodeCoords, cvf::uint vertexIndex, double tolerance)
{
std::list<std::pair<cvf::uint, bool> >::iterator it;
for(it = polygon->begin(); it != polygon->end(); ++it)
{
if (it->first == vertexIndex) return true;
}
#if 1
bool existsOrInserted = false;
for(it = polygon->begin(); it != polygon->end(); ++it)
{
if ( (nodeCoords[it->first] - nodeCoords[vertexIndex]).length() < tolerance)
{
if (vertexIndex < it->first) it->first = vertexIndex;
existsOrInserted = true;
}
}
if (existsOrInserted) return true;
#endif
// Insert vertex in polygon if the distance to one of the edges is small enough
std::list<std::pair<cvf::uint, bool> >::iterator it2;
std::list<std::pair<cvf::uint, bool> >::iterator insertBefore;
for (it = polygon->begin(); it != polygon->end(); ++it)
{
it2 = it;
++it2; insertBefore = it2; if (it2 == polygon->end()) it2 = polygon->begin();
double sqDistToLine = GeometryTools::linePointSquareDist(nodeCoords[it->first], nodeCoords[it2->first], nodeCoords[vertexIndex]);
if (fabs(sqDistToLine) < tolerance*tolerance )
{
it = polygon->insert(insertBefore, std::make_pair(vertexIndex, false));
existsOrInserted = true;
}
}
return existsOrInserted;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

View File

@ -37,9 +37,14 @@ public:
LINES_OVERLAP
};
static bool insertVertexInPolygon(std::list<std::pair<cvf::uint, bool> >* polygon, const cvf::Vec3dArray& nodeCoords, cvf::uint vertexIndex, double tolerance);
static void addMidEdgeNodes(std::list<std::pair<cvf::uint, bool> >* polygon, const cvf::Vec3dArray& nodes, EdgeSplitStorage& edgeSplitStorage, std::vector<cvf::Vec3d>* createdVertexes);
template<typename VerticeArrayType, typename IndexType>
static bool insertVertexInPolygon( std::vector<IndexType> * polygon,
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodeCoords,
IndexType vertexIndex,
double tolerance);
static IntersectionStatus inPlaneLineIntersect3D(const cvf::Vec3d& planeNormal,
const cvf::Vec3d& p1, const cvf::Vec3d& p2, const cvf::Vec3d& p3, const cvf::Vec3d& p4,
cvf::Vec3d* intersectionPoint, double* fractionAlongLine1, double* fractionAlongLine2,
@ -55,7 +60,8 @@ public:
template<typename VerticeArrayType, typename IndexType>
static bool calculateOverlapPolygonOfTwoQuads( std::vector<IndexType> * polygon, std::vector<cvf::Vec3d>* createdVertexes,
static bool calculateOverlapPolygonOfTwoQuads( std::vector<IndexType> * polygon,
std::vector<cvf::Vec3d>* createdVertexes,
EdgeIntersectStorage<IndexType>& edgeIntersectionStorage,
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodes,
const IndexType cv1CubeFaceIndices[4],
@ -76,7 +82,7 @@ template <typename IndexType>
class EdgeIntersectStorage
{
public:
void setVertexCount(IndexType size);
void setVertexCount(size_t size);
bool findIntersection( IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2,
IndexType* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus,
double* fractionAlongEdge1, double* fractionAlongEdge2);

View File

@ -3,6 +3,83 @@
namespace cvf
{
//--------------------------------------------------------------------------------------------------
/// Inserts the vertex into the polygon if it fits along one of the edges within the tolerance.
/// The method returns true if it was inserted, or if it was already in the polygon, or if it was
/// within the tolerance of an existing vertex in the polygon.
/// In the latter situation it replaces the previous vertex in the polygon.
///
/// Todo: If a vertex is replaced, the VxToCv map in TimeStepGeometry should be updated
//--------------------------------------------------------------------------------------------------
template<typename VerticeArrayType, typename IndexType>
bool GeometryTools::insertVertexInPolygon( std::vector<IndexType> * polygon,
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodeCoords,
IndexType vertexIndex,
double tolerance)
{
CVF_ASSERT(polygon);
// Check if vertex is directly included already
for(std::vector<IndexType>::iterator it = polygon->begin(); it != polygon->end(); ++it)
{
if (*it == vertexIndex) return true;
}
#if 1
// Check if the new point is within tolerance of one of the polygon vertices
bool existsOrInserted = false;
for(std::vector<IndexType>::iterator it = polygon->begin(); it != polygon->end(); ++it)
{
if ( (nodeCoords[*it] - nodeCoords[vertexIndex]).length() < tolerance)
{
if (vertexIndex < *it) *it = vertexIndex;
existsOrInserted = true;
}
}
if (existsOrInserted) return true;
#endif
// Copy the start polygon to a list
std::list<IndexType> listPolygon;
for (size_t pcIdx = 0; pcIdx < polygon->size(); ++pcIdx)
{
listPolygon.push_back((*polygon)[pcIdx]);
}
// Insert vertex in polygon if the distance to one of the edges is small enough
std::list<IndexType >::iterator it2;
std::list<IndexType >::iterator insertBefore;
for (std::list<IndexType >::iterator it = listPolygon.begin(); it != listPolygon.end(); ++it)
{
it2 = it;
++it2; insertBefore = it2; if (it2 == listPolygon.end()) it2 = listPolygon.begin();
double sqDistToLine = GeometryTools::linePointSquareDist(nodeCoords[*it], nodeCoords[*it2], nodeCoords[vertexIndex]);
if (fabs(sqDistToLine) < tolerance*tolerance )
{
it = listPolygon.insert(insertBefore, vertexIndex);
existsOrInserted = true;
}
}
// Write polygon back into the vector
polygon->clear();
for (std::list<IndexType >::iterator it = listPolygon.begin(); it != listPolygon.end(); ++it)
{
polygon->push_back(*it);
}
return existsOrInserted;
}
//--------------------------------------------------------------------------------------------------
@ -116,6 +193,8 @@ bool GeometryTools::isPointTouchingIndexedPolygon( const cvf::Vec3d& polygonNor
//--------------------------------------------------------------------------------------------------
/// Returns true if we get an actual polygon
/// The returned polygon will keep the winding from the first face.
/// The second face must have opposite winding of the first
//--------------------------------------------------------------------------------------------------
template<typename VerticeArrayType, typename IndexType>
bool GeometryTools::calculateOverlapPolygonOfTwoQuads(std::vector<IndexType> * polygon,
@ -126,12 +205,12 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads(std::vector<IndexType> * p
const IndexType cv2CubeFaceIndices[4],
double tolerance)
{
CVF_ASSERT(polygon);
CVF_ASSERT(createdVertexes);
// Topology analysis
if (createdVertexes == NULL) return false;
size_t newVertexIndex = nodes.size() + createdVertexes->size();
IndexType newVertexIndex = static_cast<IndexType>(nodes.size() + createdVertexes->size());
bool cv1VxTouchCv2[4] = { false, false, false, false };
bool cv2VxTouchCv1[4] = { false, false, false, false };
@ -387,13 +466,13 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads(std::vector<IndexType> * p
// Insert indices to the new intersection vertices into the polygon of
// this connection according to fraction along edge
std::map<double,size_t> sortingMap;
std::map<double,IndexType> sortingMap;
for (i = 0; i < intersectionFractionsAlongEdge.size(); ++i)
{
sortingMap[intersectionFractionsAlongEdge[i]] = intersectionVxIndices[i];
}
std::map<double, size_t>::iterator it;
std::map<double, IndexType>::iterator it;
for (it = sortingMap.begin(); it != sortingMap.end(); ++it)
{
if (polygon->empty() || polygon->back() != it->second)
@ -478,7 +557,7 @@ void GeometryTools::calculatePartiallyFreeCubeFacePolygon(ArrayWrapperConst<Vert
bool* m_partiallyFreeCubeFaceHasHoles)
{
// Vertex Index to position in polygon
typedef std::map< size_t, std::vector<size_t >::const_iterator > VxIdxToPolygonPositionMap;
typedef std::map< IndexType, std::vector<IndexType>::const_iterator > VxIdxToPolygonPositionMap;
CVF_ASSERT(m_partiallyFreeCubeFaceHasHoles);
CVF_ASSERT(partialFacePolygon != NULL);
@ -745,7 +824,7 @@ void GeometryTools::calculatePartiallyFreeCubeFacePolygon(ArrayWrapperConst<Vert
///
//--------------------------------------------------------------------------------------------------
template <typename IndexType>
void EdgeIntersectStorage<IndexType>::setVertexCount(IndexType size)
void EdgeIntersectStorage<IndexType>::setVertexCount(size_t size)
{
m_edgeIntsectMap.resize(size);
}
@ -759,7 +838,7 @@ void EdgeIntersectStorage<IndexType>::canonizeAddress(IndexType& e1P1, IndexType
flipE1E2 = (flipE1 ? e1P2: e1P1) > (flipE2 ? e2P2: e2P1);
static size_t temp;
static IndexType temp;
if (flipE1)
{
temp = e1P1;
@ -835,15 +914,15 @@ bool EdgeIntersectStorage<IndexType>::findIntersection(IndexType e1P1, IndexType
if (!m_edgeIntsectMap[e1P1].size()) return false;
std::map<size_t, std::map<size_t, std::map<size_t, IntersectData > > >::iterator it;
std::map<IndexType, std::map<IndexType, std::map<IndexType, IntersectData > > >::iterator it;
it = m_edgeIntsectMap[e1P1].find(e1P2);
if (it == m_edgeIntsectMap[e1P1].end()) return false;
std::map<size_t, std::map<size_t, IntersectData > >::iterator it2;
std::map<IndexType, std::map<IndexType, IntersectData > >::iterator it2;
it2 = it->second.find(e2P1);
if (it2 == it->second.end()) return false;
std::map<size_t, IntersectData >::iterator it3;
std::map<IndexType, IntersectData >::iterator it3;
it3 = it2->second.find(e2P2);
if (it3 == it2->second.end()) return false;