mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
Moved GeometryTools w tests to ReservoirDataModel
This commit is contained in:
Jacob Støren
@@ -20,6 +20,8 @@ ${CEE_CURRENT_LIST_DIR}RigStatisticsMath.h
|
||||
${CEE_CURRENT_LIST_DIR}RigWellPath.h
|
||||
${CEE_CURRENT_LIST_DIR}RigFault.h
|
||||
${CEE_CURRENT_LIST_DIR}RigNNCData.h
|
||||
${CEE_CURRENT_LIST_DIR}cvfGeometryTools.h
|
||||
${CEE_CURRENT_LIST_DIR}cvfGeometryTools.inl
|
||||
)
|
||||
|
||||
set (SOURCE_GROUP_SOURCE_FILES
|
||||
@@ -38,6 +40,7 @@ ${CEE_CURRENT_LIST_DIR}RigStatisticsMath.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigWellPath.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigFault.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigNNCData.cpp
|
||||
${CEE_CURRENT_LIST_DIR}cvfGeometryTools.cpp
|
||||
)
|
||||
|
||||
list(APPEND CODE_HEADER_FILES
|
||||
|
||||
@@ -39,6 +39,7 @@ set( UNIT_TEST_CPP_SOURCES
|
||||
RigActiveCellInfo-Test.cpp
|
||||
RigReservoir-Test.cpp
|
||||
RigStatisticsMath-Test.cpp
|
||||
cvfGeometryTools-Test.cpp
|
||||
)
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,517 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2011-2012 Statoil ASA, Ceetron AS
|
||||
//
|
||||
// ResInsight 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.
|
||||
//
|
||||
// ResInsight 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 at <http://www.gnu.org/licenses/gpl.html>
|
||||
// for more details.
|
||||
//
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#include "gtest/gtest.h"
|
||||
|
||||
#include "cvfLibCore.h"
|
||||
#include "cvfLibViewing.h"
|
||||
#include "cvfLibRender.h"
|
||||
#include "cvfLibGeometry.h"
|
||||
#include "cafFixedArray.h"
|
||||
|
||||
#include "cvfArrayWrapperToEdit.h"
|
||||
#include "cvfArrayWrapperConst.h"
|
||||
|
||||
#include "cvfGeometryTools.h"
|
||||
#include "cvfBoundingBoxTree.h"
|
||||
|
||||
using namespace cvf;
|
||||
|
||||
#if 0
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void ControlVolume::calculateCubeFaceStatus(const cvf::Vec3dArray& nodeCoords, double areaTolerance)
|
||||
{
|
||||
int cubeFace;
|
||||
cvf::uint cubeFaceIndices[4];
|
||||
for (cubeFace = 0; cubeFace < 6; ++cubeFace)
|
||||
{
|
||||
surfaceNodeIndices(static_cast<Defines::CubeFace>(cubeFace), cubeFaceIndices);
|
||||
|
||||
std::vector<const brv::Connection*> conns;
|
||||
connections(static_cast<Defines::CubeFace>(cubeFace), &conns);
|
||||
|
||||
if (!conns.size())
|
||||
{
|
||||
m_cubeFaceStatus[cubeFace] = FREE_FACE;
|
||||
}
|
||||
else
|
||||
{
|
||||
double area = 0.5 * (nodeCoords[cubeFaceIndices[1]]-nodeCoords[cubeFaceIndices[0]] ^ nodeCoords[cubeFaceIndices[3]]-nodeCoords[cubeFaceIndices[0]]).length();
|
||||
area += 0.5 * (nodeCoords[cubeFaceIndices[3]]-nodeCoords[cubeFaceIndices[2]] ^ nodeCoords[cubeFaceIndices[1]]-nodeCoords[cubeFaceIndices[2]]).length();
|
||||
double totConnectionArea = 0;
|
||||
size_t i;
|
||||
for (i = 0; i < conns.size(); ++i)
|
||||
{
|
||||
totConnectionArea += conns[i]->brfArea();
|
||||
}
|
||||
|
||||
if ( totConnectionArea < area - areaTolerance )
|
||||
{
|
||||
m_cubeFaceStatus[cubeFace] = PARTIALLY_COVERED;
|
||||
}
|
||||
else
|
||||
{
|
||||
m_cubeFaceStatus[cubeFace] = COMPLETELY_COVERED;
|
||||
}
|
||||
}
|
||||
|
||||
// Create a polygon to store the complete polygon of the faces
|
||||
// not completely covered by connections
|
||||
// This polygon will be filled with nodes later
|
||||
|
||||
if (m_cubeFaceStatus[cubeFace] != COMPLETELY_COVERED )
|
||||
{
|
||||
m_freeFacePolygons[cubeFace] = new std::list<std::pair<cvf::uint, bool> >;
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
template <typename NodeArrayType, typename NodeType, typename IndexType>
|
||||
NodeType quadNormal (ArrayWrapperConst<NodeArrayType, NodeType> nodeCoords,
|
||||
const IndexType cubeFaceIndices[4] )
|
||||
{
|
||||
return ( nodeCoords[cubeFaceIndices[2]] - nodeCoords[cubeFaceIndices[0]]) ^
|
||||
( nodeCoords[cubeFaceIndices[3]] - nodeCoords[cubeFaceIndices[1]]);
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
class QuadFaceIntersectorImplHandle
|
||||
{
|
||||
public:
|
||||
virtual ~QuadFaceIntersectorImplHandle() {}
|
||||
virtual bool intersect() = 0;
|
||||
};
|
||||
|
||||
template < typename NodeArrayType, typename NodeType, typename IndicesArrayType, typename IndicesType>
|
||||
class QuadFaceIntersectorImpl : public QuadFaceIntersectorImplHandle
|
||||
{
|
||||
public:
|
||||
QuadFaceIntersectorImpl( ArrayWrapperToEdit<NodeArrayType, NodeType> nodeArray, ArrayWrapperToEdit<IndicesArrayType, IndicesType> indices)
|
||||
: m_nodeArray(nodeArray),
|
||||
m_indices(indices){}
|
||||
|
||||
|
||||
virtual bool intersect()
|
||||
{
|
||||
size_t nodeCount = m_nodeArray.size();
|
||||
NodeType a = m_nodeArray[0];
|
||||
IndicesType idx = m_indices[0];
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
|
||||
ArrayWrapperToEdit<NodeArrayType, NodeType> m_nodeArray;
|
||||
ArrayWrapperToEdit<IndicesArrayType, IndicesType> m_indices;
|
||||
};
|
||||
|
||||
|
||||
class QuadFaceIntersector
|
||||
{
|
||||
public:
|
||||
template <typename NodeArrayType, typename NodeType, typename IndicesArrayType, typename IndicesType>
|
||||
void setup( ArrayWrapperToEdit<NodeArrayType, NodeType> nodeArray, ArrayWrapperToEdit<IndicesArrayType, IndicesType> indices)
|
||||
{
|
||||
|
||||
m_implementation = new QuadFaceIntersectorImpl< NodeArrayType, NodeType, IndicesArrayType, IndicesType>( nodeArray, indices);
|
||||
}
|
||||
|
||||
bool intersect() { return m_implementation->intersect(); }
|
||||
private:
|
||||
QuadFaceIntersectorImplHandle * m_implementation;
|
||||
};
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<cvf::Vec3d> createVertices()
|
||||
{
|
||||
std::vector<cvf::Vec3d> vxs;
|
||||
vxs.resize(14, cvf::Vec3d::ZERO);
|
||||
|
||||
vxs[ 0]= cvf::Vec3d( 0 , 0 , 0 );
|
||||
vxs[ 1]= cvf::Vec3d( 1 , 0 , 0 );
|
||||
vxs[ 2]= cvf::Vec3d( 1 , 1 , 0 );
|
||||
vxs[ 3]= cvf::Vec3d( 0 , 1 , 0 );
|
||||
vxs[ 4]= cvf::Vec3d(-0.4 ,-0.2 , 0.0 );
|
||||
vxs[ 5]= cvf::Vec3d( 0.4 , 0.6 , 0.0 );
|
||||
vxs[ 6]= cvf::Vec3d( 0.8 , 0.2 , 0.0 );
|
||||
vxs[ 7]= cvf::Vec3d( 0.0 ,-0.6 , 0.0 );
|
||||
vxs[ 8]= cvf::Vec3d( 1.0 , 1.2 , 0.0 );
|
||||
vxs[ 9]= cvf::Vec3d( 1.4 , 0.8 , 0.0 );
|
||||
vxs[10]= cvf::Vec3d( 0.4 ,-0.2 , 0.0 );
|
||||
vxs[11]= cvf::Vec3d( 1.2 , 0.6 , 0.0 );
|
||||
vxs[12]= cvf::Vec3d( 1.6 , 0.2 , 0.0 );
|
||||
vxs[13]= cvf::Vec3d( 0.8 ,-0.6 , 0.0 );
|
||||
|
||||
return vxs;
|
||||
}
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
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];
|
||||
|
||||
return cubeFaces;
|
||||
}
|
||||
|
||||
std::ostream& operator<< (std::ostream& stream, std::vector<cvf::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<cvf::uint> > overlapPolygons;
|
||||
std::vector<caf::UintArray4> faces = getCubeFaces();
|
||||
|
||||
EdgeIntersectStorage<cvf::uint> edgeIntersectionStorage;
|
||||
edgeIntersectionStorage.setVertexCount(nodes.size());
|
||||
{
|
||||
std::vector<cvf::uint> polygon;
|
||||
bool isOk = false;
|
||||
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(
|
||||
&polygon,
|
||||
&additionalVertices,
|
||||
&edgeIntersectionStorage,
|
||||
wrapArrayConst(&nodes),
|
||||
faces[0].data(),
|
||||
faces[1].data(),
|
||||
1e-6);
|
||||
|
||||
EXPECT_EQ( (size_t)5, polygon.size());
|
||||
EXPECT_EQ( (size_t)2, additionalVertices.size());
|
||||
EXPECT_TRUE(isOk);
|
||||
overlapPolygons.push_back(polygon);
|
||||
std::cout << polygon << std::endl;
|
||||
|
||||
}
|
||||
|
||||
{
|
||||
std::vector<cvf::uint> polygon;
|
||||
bool isOk = false;
|
||||
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(
|
||||
&polygon,
|
||||
&additionalVertices,
|
||||
&edgeIntersectionStorage,
|
||||
wrapArrayConst(&nodes),
|
||||
faces[0].data(),
|
||||
faces[2].data(),
|
||||
1e-6);
|
||||
|
||||
EXPECT_EQ( (size_t)5, polygon.size());
|
||||
EXPECT_EQ( (size_t)4, additionalVertices.size());
|
||||
EXPECT_TRUE(isOk);
|
||||
overlapPolygons.push_back(polygon);
|
||||
std::cout << polygon << std::endl;
|
||||
|
||||
}
|
||||
|
||||
{
|
||||
std::vector<cvf::uint> polygon;
|
||||
bool isOk = false;
|
||||
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(
|
||||
&polygon,
|
||||
&additionalVertices,
|
||||
&edgeIntersectionStorage,
|
||||
wrapArrayConst(&nodes),
|
||||
faces[0].data(),
|
||||
faces[3].data(),
|
||||
1e-6);
|
||||
|
||||
EXPECT_EQ( (size_t)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<cvf::uint> basePolygon;
|
||||
basePolygon.insert(basePolygon.begin(), faces[0].data(), &(faces[0].data()[4]));
|
||||
|
||||
for (cvf::uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx)
|
||||
{
|
||||
bool inserted = GeometryTools::insertVertexInPolygon(
|
||||
&basePolygon,
|
||||
wrapArrayConst(&nodes),
|
||||
vxIdx,
|
||||
1e-6
|
||||
);
|
||||
}
|
||||
|
||||
EXPECT_EQ( (size_t)8, basePolygon.size());
|
||||
std::cout << "Bp: " << basePolygon << std::endl;
|
||||
|
||||
for (size_t pIdx = 0; pIdx < overlapPolygons.size(); ++pIdx)
|
||||
{
|
||||
for (cvf::uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx)
|
||||
{
|
||||
bool inserted = GeometryTools::insertVertexInPolygon(
|
||||
&overlapPolygons[pIdx],
|
||||
wrapArrayConst(&nodes),
|
||||
vxIdx,
|
||||
1e-6
|
||||
);
|
||||
}
|
||||
|
||||
if (pIdx == 0)
|
||||
{
|
||||
EXPECT_EQ((size_t)5, overlapPolygons[pIdx].size());
|
||||
}
|
||||
if (pIdx == 1)
|
||||
{
|
||||
EXPECT_EQ((size_t)5, overlapPolygons[pIdx].size());
|
||||
}
|
||||
if (pIdx == 2)
|
||||
{
|
||||
EXPECT_EQ((size_t)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<cvf::uint> freeFacePolygon;
|
||||
bool hasHoles = false;
|
||||
|
||||
std::vector< std::vector<cvf::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( (size_t)4, freeFacePolygon.size());
|
||||
EXPECT_FALSE(hasHoles);
|
||||
std::cout << "FF1: " << freeFacePolygon << std::endl;
|
||||
}
|
||||
|
||||
{
|
||||
std::vector<cvf::uint> freeFacePolygon;
|
||||
bool hasHoles = false;
|
||||
|
||||
std::vector< std::vector<cvf::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( (size_t)9, freeFacePolygon.size());
|
||||
EXPECT_FALSE(hasHoles);
|
||||
|
||||
std::cout << "FF2: " << freeFacePolygon << std::endl;
|
||||
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
TEST(CellFaceIntersectionTst, Intersection)
|
||||
{
|
||||
|
||||
std::vector<cvf::Vec3d> additionalVertices;
|
||||
cvf::Vec3dArray nodes;
|
||||
std::vector<size_t> polygon;
|
||||
|
||||
|
||||
cvf::Array<size_t> ids;
|
||||
size_t cv1CubeFaceIndices[4] = {0, 1, 2, 3};
|
||||
size_t cv2CubeFaceIndices[4] = {4, 5, 6, 7};
|
||||
|
||||
nodes.resize(8);
|
||||
nodes.setAll(cvf::Vec3d(0, 0, 0));
|
||||
EdgeIntersectStorage<size_t> edgeIntersectionStorage;
|
||||
edgeIntersectionStorage.setVertexCount(nodes.size());
|
||||
|
||||
// Face 1
|
||||
nodes[0] = cvf::Vec3d(0, 0, 0);
|
||||
nodes[1] = cvf::Vec3d(1, 0, 0);
|
||||
nodes[2] = cvf::Vec3d(1, 1, 0);
|
||||
nodes[3] = cvf::Vec3d(0, 1, 0);
|
||||
// Face 2
|
||||
nodes[4] = cvf::Vec3d(0, 0, 0);
|
||||
nodes[5] = cvf::Vec3d(1, 0, 0);
|
||||
nodes[6] = cvf::Vec3d(1, 1, 0);
|
||||
nodes[7] = cvf::Vec3d(0, 1, 0);
|
||||
|
||||
|
||||
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, &edgeIntersectionStorage,
|
||||
wrapArrayConst(&nodes), cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
|
||||
EXPECT_EQ( (size_t)4, polygon.size());
|
||||
EXPECT_EQ( (size_t)0, additionalVertices.size());
|
||||
EXPECT_TRUE(isOk);
|
||||
|
||||
// Face 1
|
||||
nodes[0] = cvf::Vec3d(0, 0, 0);
|
||||
nodes[1] = cvf::Vec3d(1, 0, 0);
|
||||
nodes[2] = cvf::Vec3d(1, 1, 0);
|
||||
nodes[3] = cvf::Vec3d(0, 1, 0);
|
||||
// Face 2
|
||||
nodes[4] = cvf::Vec3d(0.5, -0.25, 0);
|
||||
nodes[5] = cvf::Vec3d(1.25, 0.5, 0);
|
||||
nodes[6] = cvf::Vec3d(0.5, 1.25, 0);
|
||||
nodes[7] = cvf::Vec3d(-0.25, 0.5, 0);
|
||||
polygon.clear();
|
||||
|
||||
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, &edgeIntersectionStorage,
|
||||
wrapArrayConst(&nodes), cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
|
||||
EXPECT_EQ( (size_t)8, polygon.size());
|
||||
EXPECT_EQ( (size_t)8, additionalVertices.size());
|
||||
EXPECT_TRUE(isOk);
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST(CellFaceIntersectionTst, FreeFacePolygon)
|
||||
{
|
||||
|
||||
std::vector<cvf::Vec3d> additionalVertices;
|
||||
cvf::Vec3dArray nodes;
|
||||
std::vector<size_t> polygon;
|
||||
|
||||
|
||||
cvf::Array<size_t> ids;
|
||||
size_t cv1CubeFaceIndices[4] = {0, 1, 2, 3};
|
||||
size_t cv2CubeFaceIndices[4] = {4, 5, 6, 7};
|
||||
|
||||
nodes.resize(8);
|
||||
nodes.setAll(cvf::Vec3d(0, 0, 0));
|
||||
EdgeIntersectStorage<size_t> edgeIntersectionStorage;
|
||||
edgeIntersectionStorage.setVertexCount(nodes.size());
|
||||
|
||||
// Face 1
|
||||
nodes[0] = cvf::Vec3d(0, 0, 0);
|
||||
nodes[1] = cvf::Vec3d(1, 0, 0);
|
||||
nodes[2] = cvf::Vec3d(1, 1, 0);
|
||||
nodes[3] = cvf::Vec3d(0, 1, 0);
|
||||
// Face 2
|
||||
nodes[4] = cvf::Vec3d(0, 0, 0);
|
||||
nodes[5] = cvf::Vec3d(1, 0, 0);
|
||||
nodes[6] = cvf::Vec3d(1, 1, 0);
|
||||
nodes[7] = cvf::Vec3d(0, 1, 0);
|
||||
|
||||
|
||||
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, &edgeIntersectionStorage,
|
||||
wrapArrayConst(&nodes), cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
|
||||
EXPECT_EQ( (size_t)4, polygon.size());
|
||||
EXPECT_EQ( (size_t)0, additionalVertices.size());
|
||||
EXPECT_TRUE(isOk);
|
||||
|
||||
std::vector< bool > faceOverlapPolygonWinding;
|
||||
std::vector< std::vector<size_t>* > faceOverlapPolygons;
|
||||
faceOverlapPolygons.push_back(&polygon);
|
||||
faceOverlapPolygonWinding.push_back(true);
|
||||
|
||||
std::vector<size_t> partialFacePolygon;
|
||||
bool hasHoles = false;
|
||||
GeometryTools::calculatePartiallyFreeCubeFacePolygon(
|
||||
wrapArrayConst(&nodes),
|
||||
wrapArrayConst(cv1CubeFaceIndices, 4),
|
||||
Vec3d(0,0,1),
|
||||
faceOverlapPolygons,
|
||||
faceOverlapPolygonWinding,
|
||||
&partialFacePolygon,
|
||||
&hasHoles);
|
||||
|
||||
// Face 1
|
||||
nodes[0] = cvf::Vec3d(0, 0, 0);
|
||||
nodes[1] = cvf::Vec3d(1, 0, 0);
|
||||
nodes[2] = cvf::Vec3d(1, 1, 0);
|
||||
nodes[3] = cvf::Vec3d(0, 1, 0);
|
||||
// Face 2
|
||||
nodes[4] = cvf::Vec3d(0.5, -0.25, 0);
|
||||
nodes[5] = cvf::Vec3d(1.25, 0.5, 0);
|
||||
nodes[6] = cvf::Vec3d(0.5, 1.25, 0);
|
||||
nodes[7] = cvf::Vec3d(-0.25, 0.5, 0);
|
||||
polygon.clear();
|
||||
|
||||
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads(&polygon, &additionalVertices, &edgeIntersectionStorage,
|
||||
wrapArrayConst(&nodes), cv1CubeFaceIndices, cv2CubeFaceIndices, 1e-6);
|
||||
EXPECT_EQ( (size_t)8, polygon.size());
|
||||
EXPECT_EQ( (size_t)8, additionalVertices.size());
|
||||
EXPECT_TRUE(isOk);
|
||||
|
||||
|
||||
|
||||
}
|
||||
925
ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp
Normal file
925
ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp
Normal file
@@ -0,0 +1,925 @@
|
||||
#include "cvfGeometryTools.h"
|
||||
|
||||
#pragma warning (disable : 4503)
|
||||
namespace cvf
|
||||
{
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
cvf::Vec3d GeometryTools::computeFaceCenter(const cvf::Vec3d& v0, const cvf::Vec3d& v1, const cvf::Vec3d& v2, const cvf::Vec3d& v3)
|
||||
{
|
||||
cvf::Vec3d centerCoord = v0;
|
||||
centerCoord += v1;
|
||||
centerCoord += v2;
|
||||
centerCoord += v3;
|
||||
centerCoord *= 0.25;
|
||||
|
||||
return centerCoord;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
int GeometryTools::findClosestAxis(const cvf::Vec3d& vec )
|
||||
{
|
||||
int closestAxis = 0;
|
||||
double maxComponent = fabs(vec.x());
|
||||
|
||||
if (fabs(vec.y()) > maxComponent)
|
||||
{
|
||||
maxComponent = (float)fabs(vec.y());
|
||||
closestAxis = 1;
|
||||
}
|
||||
|
||||
if (fabs(vec.z()) > maxComponent)
|
||||
{
|
||||
closestAxis = 2;
|
||||
}
|
||||
|
||||
return closestAxis;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Return angle between vectors if v1 x v2 is same way as normal
|
||||
/// else return 2PI - angle
|
||||
/// This means if the angle is slightly "negative", using the right hand rule, this method will return
|
||||
/// nearly 2*PI
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
double const MY_PI = 4 * atan(1.0);
|
||||
|
||||
double GeometryTools::getAngle(const cvf::Vec3d& positiveNormalAxis, const cvf::Vec3d& v1, const cvf::Vec3d& v2)
|
||||
{
|
||||
bool isOk = false;
|
||||
cvf::Vec3d v1N = v1.getNormalized(&isOk);
|
||||
if (!isOk) return 0;
|
||||
cvf::Vec3d v2N = v2.getNormalized();
|
||||
if (!isOk) return 0;
|
||||
|
||||
double cosAng = v1N * v2N;
|
||||
// Guard acos against out-of-domain input
|
||||
if (cosAng <= -1.0)
|
||||
{
|
||||
cosAng = -1.0;
|
||||
}
|
||||
else if (cosAng >= 1.0)
|
||||
{
|
||||
cosAng = 1.0;
|
||||
}
|
||||
double angle = acos(cosAng);
|
||||
|
||||
cvf::Vec3d crossProd = v1N ^ v2N;
|
||||
double sign = positiveNormalAxis * crossProd;
|
||||
if (sign < 0)
|
||||
{
|
||||
angle = 2*MY_PI - angle;
|
||||
}
|
||||
return angle;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Return angle in radians between vectors [0, Pi]
|
||||
/// If v1 or v2 is zero, the method will return 0.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
double GeometryTools::getAngle(const cvf::Vec3d& v1, const cvf::Vec3d& v2)
|
||||
{
|
||||
bool isOk = false;
|
||||
cvf::Vec3d v1N = v1.getNormalized(&isOk);
|
||||
if (!isOk) return 0;
|
||||
cvf::Vec3d v2N = v2.getNormalized();
|
||||
if (!isOk) return 0;
|
||||
|
||||
double cosAng = v1N * v2N;
|
||||
// Guard acos against out-of-domain input
|
||||
if (cosAng <= -1.0)
|
||||
{
|
||||
cosAng = -1.0;
|
||||
}
|
||||
else if (cosAng >= 1.0)
|
||||
{
|
||||
cosAng = 1.0;
|
||||
}
|
||||
double angle = acos(cosAng);
|
||||
|
||||
return angle;
|
||||
}
|
||||
|
||||
/*
|
||||
Determine the intersection point of two line segments
|
||||
From Paul Bourke, but modified to really handle coincident lines
|
||||
and lines with touching vertexes.
|
||||
Returns an intersection status telling what kind of intersection it is (if any)
|
||||
*/
|
||||
|
||||
GeometryTools::IntersectionStatus inPlaneLineIntersect(
|
||||
double x1, double y1,
|
||||
double x2, double y2,
|
||||
double x3, double y3,
|
||||
double x4, double y4,
|
||||
double l1NormalizedTolerance, double l2NormalizedTolerance,
|
||||
double *x, double *y, double* fractionAlongLine1, double* fractionAlongLine2)
|
||||
{
|
||||
double mua, mub;
|
||||
double denom, numera, numerb;
|
||||
|
||||
denom = (y4-y3) * (x2-x1) - (x4-x3) * (y2-y1);
|
||||
numera = (x4-x3) * (y1-y3) - (y4-y3) * (x1-x3);
|
||||
numerb = (x2-x1) * (y1-y3) - (y2-y1) * (x1-x3);
|
||||
|
||||
double EPS = 1e-40;
|
||||
|
||||
// Are the line coincident?
|
||||
if (fabs(numera) < EPS && fabs(numerb) < EPS && fabs(denom) < EPS)
|
||||
{
|
||||
#if 0
|
||||
*x = 0;
|
||||
*y = 0;
|
||||
*fractionAlongLine1 = 0;
|
||||
*fractionAlongLine2 = 0;
|
||||
return GeometryTools::LINES_OVERLAP;
|
||||
#else
|
||||
cvf::Vec3d p12(x2-x1, y2-y1, 0);
|
||||
cvf::Vec3d p13(x3-x1, y3-y1, 0);
|
||||
cvf::Vec3d p34(x4-x3, y4-y3, 0);
|
||||
|
||||
double length12 = p12.length();
|
||||
double length34 = p34.length();
|
||||
|
||||
// Check if the p1 p2 line is a point
|
||||
|
||||
if (length12 < EPS )
|
||||
{
|
||||
cvf::Vec3d p34(x4-x3, y4-y3, 0);
|
||||
*x = x1;
|
||||
*y = y1;
|
||||
*fractionAlongLine1 = 1;
|
||||
*fractionAlongLine2 = p13.length()/p34.length();
|
||||
return GeometryTools::LINES_OVERLAP;
|
||||
}
|
||||
|
||||
cvf::Vec3d p14(x4-x1, y4-y1, 0);
|
||||
cvf::Vec3d p32(x2-x3, y2-y3, 0);
|
||||
|
||||
cvf::Vec3d e12 = p12.getNormalized();
|
||||
double normDist13 = e12*p13 / length12;
|
||||
double normDist14 = e12*p14 / length12;
|
||||
|
||||
// Check if both points on the p3 p4 line is outside line p1 p2.
|
||||
if( (normDist13 < 0 - l1NormalizedTolerance && normDist14 < 0-l1NormalizedTolerance )|| (normDist13 > 1 +l1NormalizedTolerance && normDist14 > 1+l1NormalizedTolerance ) )
|
||||
{
|
||||
*x = 0;
|
||||
*y = 0;
|
||||
*fractionAlongLine1 = 0;
|
||||
*fractionAlongLine2 = 0;
|
||||
return GeometryTools::NO_INTERSECTION;
|
||||
}
|
||||
|
||||
double normDist32 = e12*p32 / length34;
|
||||
double normDist31 = -e12*p13 / length34;
|
||||
|
||||
// Set up fractions along lines to the edge2 vertex actually touching edge 1.
|
||||
/// if two, select the one furthest from the start
|
||||
bool pt3IsInside = false;
|
||||
bool pt4IsInside = false;
|
||||
if ((0.0 - l1NormalizedTolerance) <= normDist13 && normDist13 <= (1.0 +l1NormalizedTolerance) ) pt3IsInside = true;
|
||||
if ((0.0 - l1NormalizedTolerance) <= normDist14 && normDist14 <= (1.0 +l1NormalizedTolerance) ) pt4IsInside = true;
|
||||
|
||||
if (pt3IsInside && !pt4IsInside)
|
||||
{
|
||||
*fractionAlongLine1 = normDist13;
|
||||
*fractionAlongLine2 = 0.0;
|
||||
*x = x3;
|
||||
*y = y3;
|
||||
}
|
||||
else if (pt4IsInside && !pt3IsInside)
|
||||
{
|
||||
*fractionAlongLine1 = normDist14;
|
||||
*fractionAlongLine2 = 1.0;
|
||||
*x = x4;
|
||||
*y = y4;
|
||||
}
|
||||
else if (pt3IsInside && pt4IsInside)
|
||||
{
|
||||
// Return edge 2 vertex furthest along edge 1
|
||||
if (normDist13 <= normDist14)
|
||||
{
|
||||
*fractionAlongLine1 = normDist14 ;
|
||||
*fractionAlongLine2 = 1.0;
|
||||
*x = x4;
|
||||
*y = y4;
|
||||
}
|
||||
else
|
||||
{
|
||||
*fractionAlongLine1 = normDist13;
|
||||
*fractionAlongLine2 = 0.0;
|
||||
*x = x3;
|
||||
*y = y3;
|
||||
}
|
||||
}
|
||||
else // both outside on each side
|
||||
{
|
||||
// Return End of edge 1
|
||||
*fractionAlongLine1 = 1.0;
|
||||
*fractionAlongLine2 = normDist32;
|
||||
*x = x2;
|
||||
*y = y2;
|
||||
}
|
||||
|
||||
return GeometryTools::LINES_OVERLAP;
|
||||
#endif
|
||||
}
|
||||
|
||||
/* Are the line parallel */
|
||||
if (fabs(denom) < EPS) {
|
||||
*x = 0;
|
||||
*y = 0;
|
||||
*fractionAlongLine1 = 0;
|
||||
*fractionAlongLine2 = 0;
|
||||
|
||||
return GeometryTools::NO_INTERSECTION;
|
||||
}
|
||||
|
||||
/* Is the intersection along the the segments */
|
||||
mua = numera / denom;
|
||||
mub = numerb / denom;
|
||||
|
||||
*x = x1 + mua * (x2 - x1);
|
||||
*y = y1 + mua * (y2 - y1);
|
||||
*fractionAlongLine1 = mua;
|
||||
*fractionAlongLine2 = mub;
|
||||
|
||||
if (mua < 0 - l1NormalizedTolerance || 1 + l1NormalizedTolerance < mua || mub < 0 - l2NormalizedTolerance || 1 + l2NormalizedTolerance < mub)
|
||||
{
|
||||
return GeometryTools::LINES_INTERSECT_OUTSIDE;
|
||||
}
|
||||
else if (fabs(mua) < l1NormalizedTolerance || fabs(1-mua) < l1NormalizedTolerance ||
|
||||
fabs(mub) < l2NormalizedTolerance || fabs(1-mub) < l2NormalizedTolerance )
|
||||
{
|
||||
if (fabs(mua) < l1NormalizedTolerance) *fractionAlongLine1 = 0;
|
||||
if (fabs(1-mua) < l1NormalizedTolerance) *fractionAlongLine1 = 1;
|
||||
if (fabs(mub) < l2NormalizedTolerance) *fractionAlongLine2 = 0;
|
||||
if (fabs(1-mub) < l2NormalizedTolerance) *fractionAlongLine2 = 1;
|
||||
return GeometryTools::LINES_TOUCH;
|
||||
}
|
||||
else
|
||||
{
|
||||
return GeometryTools::LINES_CROSSES;
|
||||
}
|
||||
}
|
||||
//----------------------------------------------------------------------------------------------------------
|
||||
/// Supposed to find the intersection point if lines intersect
|
||||
/// It returns the intersection status telling if the lines only touch or are overlapping
|
||||
//----------------------------------------------------------------------------------------------------------
|
||||
|
||||
GeometryTools::IntersectionStatus
|
||||
GeometryTools::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, double tolerance)
|
||||
{
|
||||
CVF_ASSERT (intersectionPoint != NULL);
|
||||
|
||||
int Z = findClosestAxis(planeNormal);
|
||||
int X = (Z + 1) % 3;
|
||||
int Y = (Z + 2) % 3;
|
||||
double x, y;
|
||||
|
||||
// Todo: handle zero length edges
|
||||
double l1NormTol = tolerance / (p2-p1).length();
|
||||
double l2NormTol = tolerance / (p4-p3).length();
|
||||
|
||||
IntersectionStatus intersectionStatus = inPlaneLineIntersect(p1[X], p1[Y], p2[X], p2[Y], p3[X], p3[Y], p4[X], p4[Y], l1NormTol, l2NormTol, &x, &y, fractionAlongLine1, fractionAlongLine2);
|
||||
|
||||
// Check if we have a valid intersection point
|
||||
if (intersectionStatus == NO_INTERSECTION || intersectionStatus == LINES_OVERLAP)
|
||||
{
|
||||
intersectionPoint->setZero();
|
||||
}
|
||||
else
|
||||
{
|
||||
*intersectionPoint = p1 + (*fractionAlongLine1)*(p2-p1);
|
||||
}
|
||||
|
||||
return intersectionStatus;
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double GeometryTools::linePointSquareDist(const cvf::Vec3d& p1, const cvf::Vec3d& p2, const cvf::Vec3d& p3)
|
||||
{
|
||||
cvf::Vec3d v31 = p3 - p1;
|
||||
cvf::Vec3d v21 = p2 - p1;
|
||||
|
||||
double geomTolerance = 1e-24;
|
||||
if (v21.lengthSquared() < geomTolerance)
|
||||
{
|
||||
// P2 and P1 coincide, use distance from P3 to P1
|
||||
return v31.lengthSquared();
|
||||
}
|
||||
|
||||
double u = (v31*v21)/(v21*v21);
|
||||
cvf::Vec3d pOnLine(0,0,0);
|
||||
if (0 < u && u < 1) pOnLine = p1 + u*v21;
|
||||
else if (u <= 0 ) pOnLine = p1;
|
||||
else pOnLine = p2;
|
||||
|
||||
return (p3-pOnLine).lengthSquared();
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
// Copyright 2001, softSurfer (www.softsurfer.com)
|
||||
// This code may be freely used and modified for any purpose
|
||||
// providing that this copyright notice is included with it.
|
||||
// SoftSurfer makes no warranty for this code, and cannot be held
|
||||
// liable for any real or imagined damage resulting from its use.
|
||||
// Users of this code must verify correctness for their application.
|
||||
// http://www.softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm
|
||||
//
|
||||
/// Intersect a line segment with a 3D triangle
|
||||
/// Input: A line segment p0, p1. A triangle t0, t1, t2.
|
||||
/// Output: *intersectionPoint = intersection point (when it exists)
|
||||
/// Return: -1 = triangle is degenerate (a segment or point)
|
||||
/// 0 = disjoint (no intersect)
|
||||
/// 1 = intersect in unique point I1
|
||||
/// 2 = are in the same plane
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
#define SMALL_NUM 0.00000001 // anything that avoids division overflow
|
||||
// dot product (3D) which allows vector operations in arguments
|
||||
#define dot(u,v) ((u).x() * (v).x() + (u).y() * (v).y() + (u).z() * (v).z())
|
||||
|
||||
int GeometryTools::intersectLineSegmentTriangle( const cvf::Vec3d p0, const cvf::Vec3d p1,
|
||||
const cvf::Vec3d t0, const cvf::Vec3d t1, const cvf::Vec3d t2,
|
||||
cvf::Vec3d* intersectionPoint )
|
||||
{
|
||||
CVF_ASSERT(intersectionPoint != NULL);
|
||||
cvf::Vec3d u, v, n; // triangle vectors
|
||||
cvf::Vec3d dir, w0, w; // ray vectors
|
||||
double r, a, b; // params to calc ray-plane intersect
|
||||
|
||||
// get triangle edge vectors and plane normal
|
||||
u = t1 - t0;
|
||||
v = t2 - t0;
|
||||
n = u ^ v; // cross product
|
||||
if (n == cvf::Vec3d::ZERO) // triangle is degenerate
|
||||
return -1; // do not deal with this case
|
||||
|
||||
dir = p1 - p0; // ray direction vector
|
||||
w0 = p0 - t0;
|
||||
a = -dot(n, w0);
|
||||
b = dot(n, dir);
|
||||
if (fabs(b) < SMALL_NUM) { // ray is parallel to triangle plane
|
||||
if (a == 0) // ray lies in triangle plane
|
||||
return 2;
|
||||
else return 0; // ray disjoint from plane
|
||||
}
|
||||
|
||||
// get intersect point of ray with triangle plane
|
||||
r = a / b;
|
||||
if (r < 0.0) // ray goes away from triangle
|
||||
return 0; // => no intersect
|
||||
|
||||
if (r > 1.0) // Line segment does not reach triangle
|
||||
return 0;
|
||||
|
||||
*intersectionPoint = p0 + r * dir; // intersect point of ray and plane
|
||||
|
||||
// is I inside T?
|
||||
double uu, uv, vv, wu, wv, D;
|
||||
uu = dot(u, u);
|
||||
uv = dot(u, v);
|
||||
vv = dot(v, v);
|
||||
w = *intersectionPoint - t0;
|
||||
wu = dot(w, u);
|
||||
wv = dot(w, v);
|
||||
D = uv * uv - uu * vv;
|
||||
|
||||
// get and test parametric coords
|
||||
double s, t;
|
||||
s = (uv * wv - vv * wu) / D;
|
||||
if (s < 0.0 || s > 1.0) // I is outside T
|
||||
return 0;
|
||||
|
||||
t = (uv * wu - uu * wv) / D;
|
||||
if (t < 0.0 || (s + t) > 1.0) // I is outside T
|
||||
return 0;
|
||||
|
||||
return 1; // I is in T
|
||||
}
|
||||
|
||||
/*
|
||||
// t0 = (x0, y0, z0)
|
||||
// t1 = (x1, y1, z1)
|
||||
// t2 = (x2, y2, z2)
|
||||
//
|
||||
// p = (xp, yp, zp)
|
||||
|
||||
cvf::Vec3d barycentricCoordsExperiment(const cvf::Vec3d& t0, const cvf::Vec3d& t1, const cvf::Vec3d& t2, const cvf::Vec3d& p)
|
||||
{
|
||||
det = x0(y1*z2 - y2*z1) + x1(y2*z0 - z2*y0) + x2(y0*z1 - y1*z0);
|
||||
|
||||
b0 = ((x1 * y2 - x2*y1)*zp + xp*(y1*z2-y2*z1) + yp*(x2*z1-x1*z2)) / det;
|
||||
b1 = ((x2 * y0 - x0*y2)*zp + xp*(y2*z0-y0*z2) + yp*(x0*z2-x2*z0)) / det;
|
||||
b2 = ((x0 * y1 - x1*y0)*zp + xp*(y0*z1-y1*z0) + yp*(x1*z0-x0*z1)) / det;
|
||||
}
|
||||
|
||||
*/
|
||||
|
||||
inline double TriArea2D(double x1, double y1, double x2, double y2, double x3, double y3)
|
||||
{
|
||||
return (x1-x2)*(y2-y3) - (x2-x3)*(y1-y2);
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
// Compute barycentric coordinates (area coordinates) (u, v, w) for
|
||||
// point p with respect to triangle (t0, t1, t2)
|
||||
// These can be used as weights for interpolating scalar values across the triangle
|
||||
// Based on section 3.4 in "Real Time collision detection" by Christer Ericson
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
cvf::Vec3d GeometryTools::barycentricCoords(const cvf::Vec3d& t0, const cvf::Vec3d& t1, const cvf::Vec3d& t2, const cvf::Vec3d& p)
|
||||
{
|
||||
// Unnormalized triangle normal
|
||||
cvf::Vec3d m = (t1 - t0 ^ t2 - t0);
|
||||
|
||||
// Absolute components for determining projection plane
|
||||
int X = 0, Y = 1, Z = 2;
|
||||
Z = findClosestAxis(m);
|
||||
switch (Z)
|
||||
{
|
||||
case 0: X = 1; Y = 2; break; // x is largest, project to the yz plane
|
||||
case 1: X = 0; Y = 2; break; // y is largest, project to the xz plane
|
||||
case 2: X = 0; Y = 1; break; // z is largest, project to the xy plane
|
||||
}
|
||||
|
||||
// Compute areas in plane of largest projection
|
||||
// Nominators and one-over-denominator for u and v ratios
|
||||
double nu, nv, ood;
|
||||
nu = TriArea2D(p[X], p[Y], t1[X], t1[Y], t2[X], t2[Y]); // Area of PBC in yz plane
|
||||
nv = TriArea2D(p[X], p[Y], t2[X], t2[Y], t0[X], t0[Y]); // Area of PCA in yz plane
|
||||
ood = 1.0f / m[Z]; // 1/(2*area of ABC in yz plane)
|
||||
|
||||
if (Z == 1) ood = -ood; // For some reason not explained
|
||||
|
||||
// Normalize
|
||||
|
||||
m[0] = nu * ood;
|
||||
m[1] = nv * ood;
|
||||
m[2] = 1.0f - m[0] - m[1];
|
||||
|
||||
return m;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void GeometryTools::addMidEdgeNodes(std::list<std::pair<cvf::uint, bool> >* polygon, const cvf::Vec3dArray& nodes, EdgeSplitStorage& edgeSplitStorage, std::vector<cvf::Vec3d>* createdVertexes)
|
||||
{
|
||||
size_t newVertexIndex = nodes.size() + createdVertexes->size();
|
||||
std::list<std::pair<cvf::uint, bool> >::iterator it;
|
||||
std::list<std::pair<cvf::uint, bool> >::iterator it2;
|
||||
|
||||
cvf::Vec3d midEdgeCoord(0,0,0);
|
||||
size_t midPointIndex = cvf::UNDEFINED_UINT;
|
||||
|
||||
for (it = polygon->begin(); it != polygon->end(); ++it)
|
||||
{
|
||||
it2 = it;
|
||||
++it2; if (it2 == polygon->end()) it2 = polygon->begin();
|
||||
|
||||
// Find or Create and add a mid-edge node
|
||||
|
||||
if (!edgeSplitStorage.findSplitPoint(it->first, it2->first, &midPointIndex))
|
||||
{
|
||||
|
||||
midEdgeCoord.setZero();
|
||||
midEdgeCoord += (it->first < nodes.size()) ? nodes[it->first] : (*createdVertexes)[it->first - nodes.size()];
|
||||
midEdgeCoord += (it2->first < nodes.size()) ? nodes[it2->first] : (*createdVertexes)[it2->first - nodes.size()];
|
||||
midEdgeCoord *= 0.5;
|
||||
|
||||
midPointIndex = newVertexIndex;
|
||||
createdVertexes->push_back(midEdgeCoord);
|
||||
++newVertexIndex;
|
||||
|
||||
edgeSplitStorage.addSplitPoint(it->first, it2->first, midPointIndex);
|
||||
}
|
||||
|
||||
if (it2 != polygon->begin())
|
||||
polygon->insert(it2, std::make_pair((cvf::uint)midPointIndex, true));
|
||||
else
|
||||
polygon->insert(polygon->end(), std::make_pair((cvf::uint)midPointIndex, true));
|
||||
|
||||
++it;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void EdgeSplitStorage::setVertexCount(size_t size)
|
||||
{
|
||||
m_edgeSplitMap.resize(size);
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
bool EdgeSplitStorage::findSplitPoint(size_t edgeP1Index, size_t edgeP2Index, size_t* splitPointIndex)
|
||||
{
|
||||
canonizeAddress(edgeP1Index, edgeP2Index);
|
||||
CVF_ASSERT(edgeP1Index < m_edgeSplitMap.size());
|
||||
|
||||
std::map< size_t, size_t >::iterator it;
|
||||
|
||||
it = m_edgeSplitMap[edgeP1Index].find(edgeP2Index);
|
||||
if (it == m_edgeSplitMap[edgeP1Index].end()) return false;
|
||||
|
||||
*splitPointIndex = it->second;
|
||||
return true;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void EdgeSplitStorage::addSplitPoint(size_t edgeP1Index, size_t edgeP2Index, size_t splitPointIndex)
|
||||
{
|
||||
canonizeAddress(edgeP1Index, edgeP2Index);
|
||||
CVF_ASSERT(edgeP1Index < m_edgeSplitMap.size());
|
||||
m_edgeSplitMap[edgeP1Index][edgeP2Index] = splitPointIndex;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void EdgeSplitStorage::canonizeAddress(size_t& edgeP1Index, size_t& edgeP2Index)
|
||||
{
|
||||
if (edgeP1Index > edgeP2Index)
|
||||
{
|
||||
size_t tmp;
|
||||
tmp = edgeP1Index;
|
||||
edgeP1Index = edgeP2Index;
|
||||
edgeP2Index = tmp;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
EarClipTesselator::EarClipTesselator():
|
||||
m_X(-1),
|
||||
m_Y(-1),
|
||||
m_areaTolerance(1e-12),
|
||||
m_nodeCoords(NULL)
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// \brief Do the main processing/actual triangulation
|
||||
/// \param triangleIndices Array that will receive the indices of the triangles resulting from the triangulation
|
||||
/// \return true when a tesselation was successully created
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
bool EarClipTesselator::calculateTriangles( std::vector<size_t>* triangleIndices )
|
||||
{
|
||||
CVF_ASSERT(m_nodeCoords != NULL);
|
||||
CVF_ASSERT(m_X > -1 && m_Y > -1);
|
||||
|
||||
size_t numVertices = m_polygonIndices.size();
|
||||
|
||||
if (numVertices < 3) return false;
|
||||
|
||||
// We want m_polygonIndices to be a counter-clockwise polygon to make the validation test work
|
||||
|
||||
if (calculatePolygonArea() < 0 )
|
||||
{
|
||||
m_polygonIndices.reverse();
|
||||
}
|
||||
|
||||
std::list<size_t>::iterator u, v, w;
|
||||
|
||||
// If we loop two times around polygon without clipping a single triangle we are toast.
|
||||
size_t count = 2*numVertices; // error detection
|
||||
|
||||
v = m_polygonIndices.end(); //nv - 1;
|
||||
--v;
|
||||
|
||||
while (numVertices > 2)
|
||||
{
|
||||
// if we loop, it is probably a non-simple polygon
|
||||
if (count <= 0 )
|
||||
{
|
||||
// Triangulate: ERROR - probable bad polygon!
|
||||
return false;
|
||||
}
|
||||
--count;
|
||||
|
||||
// Three consecutive vertices in current polygon, <u,v,w>
|
||||
// previous
|
||||
u = v;
|
||||
if (u == m_polygonIndices.end()) u = m_polygonIndices.begin(); // if (nv <= u) u = 0;
|
||||
|
||||
// new v
|
||||
v = u; ++v; //u + 1;
|
||||
if (v == m_polygonIndices.end()) v = m_polygonIndices.begin(); //if (nv <= v) v = 0;
|
||||
|
||||
// next
|
||||
w = v; ++w; //v + 1;
|
||||
if (w == m_polygonIndices.end()) w = m_polygonIndices.begin(); //if (nv <= w) w = 0;
|
||||
|
||||
|
||||
if ( isTriangleValid(u, v, w) )
|
||||
{
|
||||
// Indices of the vertices
|
||||
triangleIndices->push_back(*u);
|
||||
triangleIndices->push_back(*v);
|
||||
triangleIndices->push_back(*w);
|
||||
|
||||
// Remove v from remaining polygon
|
||||
m_polygonIndices.erase(v);
|
||||
v = w;
|
||||
numVertices--;
|
||||
|
||||
// Resets error detection counter
|
||||
count = 2*numVertices;
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Is this a valid triangle ? ( No points inside, and points not on a line. )
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
bool EarClipTesselator::isTriangleValid( std::list<size_t>::const_iterator u, std::list<size_t>::const_iterator v, std::list<size_t>::const_iterator w) const
|
||||
{
|
||||
CVF_ASSERT(m_X > -1 && m_Y > -1);
|
||||
|
||||
cvf::Vec3d A = (*m_nodeCoords)[*u];
|
||||
cvf::Vec3d B = (*m_nodeCoords)[*v];
|
||||
cvf::Vec3d C = (*m_nodeCoords)[*w];
|
||||
|
||||
if ( m_areaTolerance > (((B[m_X]-A[m_X])*(C[m_Y]-A[m_Y])) - ((B[m_Y]-A[m_Y])*(C[m_X]-A[m_X]))) ) return false;
|
||||
|
||||
std::list<size_t>::const_iterator c;
|
||||
std::list<size_t>::const_iterator outside;
|
||||
for (c = m_polygonIndices.begin(); c != m_polygonIndices.end(); ++c)
|
||||
{
|
||||
// The polygon points that actually make up the triangle candidate does not count
|
||||
// (but the same points on different positions in the polygon does!
|
||||
// Except those one off the triangle, that references the start or end of the triangle)
|
||||
|
||||
if ( (c == u) || (c == v) || (c == w)) continue;
|
||||
|
||||
// Originally the below tests was not included which resulted in missing triangles sometimes
|
||||
|
||||
outside = w; ++outside; if (outside == m_polygonIndices.end()) outside = m_polygonIndices.begin();
|
||||
if (c == outside && *c == *u)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
outside = u; if (outside == m_polygonIndices.begin()) outside = m_polygonIndices.end(); --outside;
|
||||
if (c == outside && *c == *w)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
cvf::Vec3d P = (*m_nodeCoords)[*c];
|
||||
|
||||
if (isPointInsideTriangle(A, B, C, P)) return false;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Decides if a point P is inside of the triangle defined by A, B, C.
|
||||
/// By calculating the "double area" (cross product) of Corner to corner x Corner to point vectors
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
bool EarClipTesselator::isPointInsideTriangle(const cvf::Vec3d& A, const cvf::Vec3d& B, const cvf::Vec3d& C, const cvf::Vec3d& P) const
|
||||
{
|
||||
CVF_ASSERT(m_X > -1 && m_Y > -1);
|
||||
|
||||
double ax = C[m_X] - B[m_X]; double ay = C[m_Y] - B[m_Y];
|
||||
double bx = A[m_X] - C[m_X]; double by = A[m_Y] - C[m_Y];
|
||||
double cx = B[m_X] - A[m_X]; double cy = B[m_Y] - A[m_Y];
|
||||
|
||||
double apx= P[m_X] - A[m_X]; double apy= P[m_Y] - A[m_Y];
|
||||
double bpx= P[m_X] - B[m_X]; double bpy= P[m_Y] - B[m_Y];
|
||||
double cpx= P[m_X] - C[m_X]; double cpy= P[m_Y] - C[m_Y];
|
||||
|
||||
double aCROSSbp = ax*bpy - ay*bpx;
|
||||
double cCROSSap = cx*apy - cy*apx;
|
||||
double bCROSScp = bx*cpy - by*cpx;
|
||||
double tol = 0;
|
||||
return ((aCROSSbp >= tol) && (bCROSScp >= tol) && (cCROSSap >= tol));
|
||||
};
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Computes area of the currently stored 2D polygon/contour
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
double EarClipTesselator::calculatePolygonArea() const
|
||||
{
|
||||
CVF_ASSERT(m_X > -1 && m_Y > -1);
|
||||
|
||||
double A = 0;
|
||||
|
||||
std::list<size_t>::const_iterator p = m_polygonIndices.end();
|
||||
--p;
|
||||
|
||||
std::list<size_t>::const_iterator q = m_polygonIndices.begin();
|
||||
while (q != m_polygonIndices.end())
|
||||
{
|
||||
A += (*m_nodeCoords)[*p][m_X] * (*m_nodeCoords)[*q][m_Y] - (*m_nodeCoords)[*q][m_X]*(*m_nodeCoords)[*p][m_Y];
|
||||
|
||||
p = q;
|
||||
q++;
|
||||
}
|
||||
|
||||
return A*0.5;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void EarClipTesselator::setNormal(const cvf::Vec3d& polygonNormal)
|
||||
{
|
||||
int Z = GeometryTools::findClosestAxis(polygonNormal);
|
||||
m_X = (Z + 1) % 3;
|
||||
m_Y = (Z + 2) % 3;
|
||||
m_polygonNormal = polygonNormal;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void EarClipTesselator::setPolygonIndices(const std::list<size_t>& polygon)
|
||||
{
|
||||
m_polygonIndices = polygon;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void EarClipTesselator::setPolygonIndices(const std::vector<size_t>& polygon)
|
||||
{
|
||||
size_t i;
|
||||
for (i = 0; i < polygon.size(); ++i)
|
||||
{
|
||||
m_polygonIndices.push_back(polygon[i]);
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void EarClipTesselator::setMinTriangleArea(double areaTolerance)
|
||||
{
|
||||
m_areaTolerance = 2*areaTolerance; // Convert to trapesoidal area
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void EarClipTesselator::setGlobalNodeArray(const cvf::Vec3dArray& nodeCoords)
|
||||
{
|
||||
m_nodeCoords = &nodeCoords;
|
||||
}
|
||||
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
FanEarClipTesselator::FanEarClipTesselator() :
|
||||
m_centerNodeIndex(std::numeric_limits<size_t>::max())
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void FanEarClipTesselator::setCenterNode(size_t centerNodeIndex)
|
||||
{
|
||||
m_centerNodeIndex = centerNodeIndex;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
bool FanEarClipTesselator::calculateTriangles(std::vector<size_t>* triangles)
|
||||
{
|
||||
CVF_ASSERT(m_centerNodeIndex != std::numeric_limits<size_t>::max());
|
||||
CVF_ASSERT(m_nodeCoords != NULL);
|
||||
CVF_ASSERT(m_X > -1 && m_Y > -1);
|
||||
|
||||
size_t nv = m_polygonIndices.size();
|
||||
|
||||
if (nv < 3) return false;
|
||||
|
||||
// We want m_polygonIndices to be a counter-clockwise polygon to make the validation test work
|
||||
|
||||
if (calculatePolygonArea() < 0 )
|
||||
{
|
||||
m_polygonIndices.reverse();
|
||||
}
|
||||
|
||||
std::list<size_t>::const_iterator it1;
|
||||
std::list<size_t>::const_iterator it2;
|
||||
|
||||
std::list< std::list<size_t> > restPolygons;
|
||||
bool wasPreviousTriangleValid = true;
|
||||
|
||||
for (it1 = m_polygonIndices.begin(); it1 != m_polygonIndices.end(); it1++)
|
||||
{
|
||||
it2 = it1;
|
||||
it2++;
|
||||
|
||||
if (it2 == m_polygonIndices.end()) it2 = m_polygonIndices.begin();
|
||||
|
||||
if (isTriangleValid(*it1, *it2, m_centerNodeIndex))
|
||||
{
|
||||
triangles->push_back(*it1);
|
||||
triangles->push_back(*it2);
|
||||
triangles->push_back(m_centerNodeIndex);
|
||||
wasPreviousTriangleValid = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (wasPreviousTriangleValid)
|
||||
{
|
||||
// Create new rest polygon.
|
||||
restPolygons.push_back(std::list<size_t>());
|
||||
restPolygons.back().push_back(m_centerNodeIndex);
|
||||
restPolygons.back().push_back(*it1);
|
||||
restPolygons.back().push_back(*it2);
|
||||
}
|
||||
else
|
||||
{
|
||||
restPolygons.back().push_back(*it2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
EarClipTesselator triMaker;
|
||||
triMaker.setNormal(m_polygonNormal);
|
||||
triMaker.setMinTriangleArea(m_areaTolerance);
|
||||
triMaker.setGlobalNodeArray(*m_nodeCoords);
|
||||
std::list< std::list<size_t> >::iterator rpIt;
|
||||
|
||||
for (rpIt = restPolygons.begin(); rpIt != restPolygons.end(); ++rpIt)
|
||||
{
|
||||
triMaker.setPolygonIndices(*rpIt);
|
||||
triMaker.calculateTriangles(triangles);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// This needs to be rewritten because we need to test for crossing edges, not only point inside.
|
||||
/// In addition the test for polygon
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
bool FanEarClipTesselator::isTriangleValid(size_t u, size_t v, size_t w)
|
||||
{
|
||||
CVF_ASSERT(m_X > -1 && m_Y > -1);
|
||||
|
||||
cvf::Vec3d A = (*m_nodeCoords)[u];
|
||||
cvf::Vec3d B = (*m_nodeCoords)[v];
|
||||
cvf::Vec3d C = (*m_nodeCoords)[w];
|
||||
|
||||
if ( m_areaTolerance > (((B[m_X]-A[m_X])*(C[m_Y]-A[m_Y])) - ((B[m_Y]-A[m_Y])*(C[m_X]-A[m_X]))) ) return false;
|
||||
|
||||
std::list<size_t>::const_iterator c;
|
||||
for (c = m_polygonIndices.begin(); c != m_polygonIndices.end(); ++c)
|
||||
{
|
||||
// The polygon points that actually make up the triangle candidate does not count
|
||||
// (but the same points on different positions in the polygon does! )
|
||||
// Todo so this test below is to accepting !! Bug !!
|
||||
if ( (*c == u) || (*c == v) || (*c == w)) continue;
|
||||
|
||||
cvf::Vec3d P = (*m_nodeCoords)[*c];
|
||||
|
||||
if (isPointInsideTriangle(A, B, C, P)) return false;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
167
ApplicationCode/ReservoirDataModel/cvfGeometryTools.h
Normal file
167
ApplicationCode/ReservoirDataModel/cvfGeometryTools.h
Normal file
@@ -0,0 +1,167 @@
|
||||
#pragma once
|
||||
#include "cvfBase.h"
|
||||
#include "cvfArray.h"
|
||||
#include <list>
|
||||
#include <map>
|
||||
#include "cvfArrayWrapperConst.h"
|
||||
|
||||
namespace cvf
|
||||
{
|
||||
|
||||
|
||||
class EdgeSplitStorage;
|
||||
template <typename IndexType> class EdgeIntersectStorage;
|
||||
|
||||
class GeometryTools
|
||||
{
|
||||
public:
|
||||
static cvf::Vec3d computeFaceCenter(const cvf::Vec3d& v0, const cvf::Vec3d& v1, const cvf::Vec3d& v2, const cvf::Vec3d& v3);
|
||||
|
||||
static double linePointSquareDist(const cvf::Vec3d& p1, const cvf::Vec3d& p2, const cvf::Vec3d& p3);
|
||||
static int intersectLineSegmentTriangle( const cvf::Vec3d p0, const cvf::Vec3d p1,
|
||||
const cvf::Vec3d t0, const cvf::Vec3d t1, const cvf::Vec3d t2,
|
||||
cvf::Vec3d* intersectionPoint );
|
||||
static cvf::Vec3d barycentricCoords(const cvf::Vec3d& t0, const cvf::Vec3d& t1, const cvf::Vec3d& t2, const cvf::Vec3d& p);
|
||||
|
||||
static int findClosestAxis(const cvf::Vec3d& vec );
|
||||
static double getAngle(const cvf::Vec3d& positiveNormalAxis, const cvf::Vec3d& v1, const cvf::Vec3d& v2);
|
||||
static double getAngle(const cvf::Vec3d& v1, const cvf::Vec3d& v2);
|
||||
|
||||
enum IntersectionStatus
|
||||
{
|
||||
NO_INTERSECTION,
|
||||
LINES_INTERSECT_OUTSIDE,
|
||||
LINES_TOUCH,
|
||||
LINES_CROSSES,
|
||||
LINES_OVERLAP
|
||||
};
|
||||
|
||||
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,
|
||||
double tolerance = 1e-6);
|
||||
|
||||
template<typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
|
||||
static bool isPointTouchingIndexedPolygon( const cvf::Vec3d& polygonNormal,
|
||||
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> vertices,
|
||||
ArrayWrapperConst<PolygonArrayType, IndexType> indices,
|
||||
const cvf::Vec3d& point,
|
||||
int* touchedEdgeIndex,
|
||||
double tolerance = 1e-6);
|
||||
|
||||
|
||||
template<typename VerticeArrayType, typename IndexType>
|
||||
static bool calculateOverlapPolygonOfTwoQuads( std::vector<IndexType> * polygon,
|
||||
std::vector<cvf::Vec3d>* createdVertexes,
|
||||
EdgeIntersectStorage<IndexType>* edgeIntersectionStorage,
|
||||
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodes,
|
||||
const IndexType cv1CubeFaceIndices[4],
|
||||
const IndexType cv2CubeFaceIndices[4],
|
||||
double tolerance);
|
||||
|
||||
template<typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
|
||||
static void calculatePartiallyFreeCubeFacePolygon(ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodeCoords,
|
||||
ArrayWrapperConst<PolygonArrayType, IndexType> completeFacePolygon,
|
||||
const cvf::Vec3d& faceNormal,
|
||||
const std::vector< std::vector<IndexType>* >& faceOverlapPolygons,
|
||||
const std::vector<bool> faceOverlapPolygonWindingSameAsCubeFaceFlags,
|
||||
std::vector<IndexType>* partialFacePolygon,
|
||||
bool* m_partiallyFreeCubeFaceHasHoles);
|
||||
};
|
||||
|
||||
template <typename IndexType>
|
||||
class EdgeIntersectStorage
|
||||
{
|
||||
public:
|
||||
void setVertexCount(size_t size);
|
||||
bool findIntersection( IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2,
|
||||
IndexType* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus,
|
||||
double* fractionAlongEdge1, double* fractionAlongEdge2);
|
||||
void addIntersection( IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2,
|
||||
IndexType vxIndexIntersectionPoint, GeometryTools::IntersectionStatus intersectionStatus,
|
||||
double fractionAlongEdge1, double fractionAlongEdge2);
|
||||
|
||||
private:
|
||||
struct IntersectData
|
||||
{
|
||||
IndexType intersectionPointIndex;
|
||||
GeometryTools::IntersectionStatus intersectionStatus;
|
||||
double fractionAlongEdge1;
|
||||
double fractionAlongEdge2;
|
||||
};
|
||||
|
||||
void canonizeAddress(IndexType& e1P1, IndexType& e1P2, IndexType& e2P1, IndexType& e2P2, bool& flipE1, bool& flipE2, bool& flipE1E2);
|
||||
|
||||
// A map containing the intersection data. The addressing is :
|
||||
// ( when leastVxIdxEdge1 < leastVxIdxEdge2 )
|
||||
// leastVxIdxEdge1, largestVxIdxEdge1, leastVxIdxEdge2, largestVxIdxEdge2, { vxIdxIntersection, fractionAlongEdg1, fractionAlonEdge2 }
|
||||
|
||||
std::vector< std::map<IndexType, std::map<IndexType, std::map<IndexType, IntersectData > > > > m_edgeIntsectMap;
|
||||
};
|
||||
|
||||
class EdgeSplitStorage
|
||||
{
|
||||
public:
|
||||
void setVertexCount(size_t size);
|
||||
bool findSplitPoint(size_t edgeP1Index, size_t edgeP2Index, size_t* splitPointIndex);
|
||||
void addSplitPoint(size_t edgeP1Index, size_t edgeP2Index, size_t splitPointIndex);
|
||||
|
||||
private:
|
||||
void canonizeAddress(size_t& e1P1, size_t& e1P2);
|
||||
|
||||
// Least VxIdx, LargestVxIdx, VertexIdx of splitpoint
|
||||
std::vector< std::map< size_t, size_t > > m_edgeSplitMap;
|
||||
};
|
||||
|
||||
|
||||
class EarClipTesselator
|
||||
{
|
||||
public:
|
||||
EarClipTesselator();
|
||||
void setNormal(const cvf::Vec3d& polygonNormal );
|
||||
void setMinTriangleArea(double areaTolerance);
|
||||
void setGlobalNodeArray(const cvf::Vec3dArray& nodeCoords);
|
||||
|
||||
void setPolygonIndices(const std::list<size_t>& polygon);
|
||||
void setPolygonIndices(const std::vector<size_t>& polygon);
|
||||
|
||||
virtual bool calculateTriangles(std::vector<size_t>* triangles);
|
||||
|
||||
protected:
|
||||
bool isTriangleValid( std::list<size_t>::const_iterator u, std::list<size_t>::const_iterator v, std::list<size_t>::const_iterator w) const;
|
||||
bool isPointInsideTriangle(const cvf::Vec3d& A, const cvf::Vec3d& B, const cvf::Vec3d& C, const cvf::Vec3d& P) const;
|
||||
double calculatePolygonArea() const;
|
||||
|
||||
protected:
|
||||
std::list<size_t> m_polygonIndices;
|
||||
const cvf::Vec3dArray* m_nodeCoords;
|
||||
int m_X, m_Y; // Index shift in vector to do simple 2D projection
|
||||
cvf::Vec3d m_polygonNormal;
|
||||
double m_areaTolerance;
|
||||
|
||||
};
|
||||
|
||||
|
||||
class FanEarClipTesselator : public EarClipTesselator
|
||||
{
|
||||
public:
|
||||
FanEarClipTesselator();
|
||||
void setCenterNode(size_t centerNodeIndex );
|
||||
|
||||
virtual bool calculateTriangles(std::vector<size_t>* triangles);
|
||||
private:
|
||||
bool isTriangleValid( size_t u, size_t v, size_t w);
|
||||
size_t m_centerNodeIndex;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#include "cvfGeometryTools.inl"
|
||||
957
ApplicationCode/ReservoirDataModel/cvfGeometryTools.inl
Normal file
957
ApplicationCode/ReservoirDataModel/cvfGeometryTools.inl
Normal file
@@ -0,0 +1,957 @@
|
||||
#include <cmath>
|
||||
|
||||
|
||||
#pragma warning (disable : 4503)
|
||||
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(typename 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(typename 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
|
||||
|
||||
typename std::list<IndexType >::iterator it2;
|
||||
typename std::list<IndexType >::iterator insertBefore;
|
||||
|
||||
for (typename 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 (typename std::list<IndexType >::iterator it = listPolygon.begin(); it != listPolygon.end(); ++it)
|
||||
{
|
||||
polygon->push_back(*it);
|
||||
}
|
||||
|
||||
return existsOrInserted;
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// \brief Test if a point touches a polygon within the specified tolerance
|
||||
///
|
||||
/// \param polygonNorm Polygon normal
|
||||
/// \param pPolygonVerts Array of polygon vertice coordinates
|
||||
/// \param piVertexIndices Array of integer node indices for this polygon
|
||||
/// \param iNumVerts Number of vertices in polygon
|
||||
/// \param point The point to be checked
|
||||
/// \param tolerance Tolerance in length
|
||||
/// \param touchedEdgeIndex returns -1 if point is inside, and edge index if point touches an edge.
|
||||
/// \return true if point lies inside or on the border of the polygon.
|
||||
///
|
||||
/// \assumpt Assumes that the polygon is planar
|
||||
/// \comment First check if point is on an edge, Then check if it is inside by
|
||||
/// counting the number of times a ray from point along positive X axis
|
||||
/// crosses an edge. Odd number says inside.
|
||||
/// \author SP (really by Eric Haines) and JJS
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
template<typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
|
||||
bool GeometryTools::isPointTouchingIndexedPolygon( const cvf::Vec3d& polygonNormal,
|
||||
cvf::ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> vertices,
|
||||
cvf::ArrayWrapperConst<PolygonArrayType, IndexType> indices,
|
||||
const cvf::Vec3d& point,
|
||||
int* touchedEdgeIndex,
|
||||
double tolerance)
|
||||
{
|
||||
size_t numIndices = indices.size();
|
||||
|
||||
int Z = findClosestAxis(polygonNormal);
|
||||
int X = (Z + 1) % 3;
|
||||
int Y = (Z + 2) % 3;
|
||||
|
||||
int crossings;
|
||||
|
||||
int xBelowVx0;
|
||||
int yBelowVx0;
|
||||
int yBelowVx1 = 0;
|
||||
|
||||
const double* vtx0;
|
||||
const double* vtx1 = NULL;
|
||||
|
||||
double dv0;
|
||||
|
||||
cvf::uint j;
|
||||
|
||||
// Check if point is on an edge or vertex
|
||||
size_t firstIdx;
|
||||
size_t secondIdx;
|
||||
|
||||
CVF_TIGHT_ASSERT(touchedEdgeIndex);
|
||||
|
||||
*touchedEdgeIndex = -1;
|
||||
for (firstIdx = 0, secondIdx = 1; firstIdx < numIndices; ++firstIdx, ++secondIdx)
|
||||
{
|
||||
if (secondIdx >= numIndices) secondIdx = 0;
|
||||
const cvf::Vec3d& vx0 = vertices[indices[firstIdx]];
|
||||
const cvf::Vec3d& vx1 = vertices[indices[secondIdx]];
|
||||
|
||||
double sqDist = GeometryTools::linePointSquareDist(vx0, vx1, point);
|
||||
if (sqDist < tolerance*tolerance)
|
||||
{
|
||||
*touchedEdgeIndex = static_cast<int>(firstIdx);
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
vtx0 = vertices[indices[numIndices-1]].ptr();
|
||||
|
||||
// get test bit for above/below Y axis. Y of Point is under Y of vtx0
|
||||
yBelowVx0 = ( dv0 = vtx0[Y] - point[Y] ) >= 0.0;
|
||||
|
||||
crossings = 0;
|
||||
for (j = 0; j < numIndices; j++)
|
||||
{
|
||||
// cleverness: bobble between filling endpoints of edges, so that the previous edge's shared endpoint is maintained.
|
||||
if (j & 0x1)
|
||||
{
|
||||
vtx0 = vertices[indices[j]].ptr();
|
||||
yBelowVx0 = (dv0 = vtx0[Y] - point[Y]) >= 0.0;
|
||||
}
|
||||
else
|
||||
{
|
||||
vtx1 = vertices[indices[j]].ptr();
|
||||
yBelowVx1 = (vtx1[Y] >= point[Y]);
|
||||
}
|
||||
|
||||
// check if Y of point is between Y of Vx0 and Vx1
|
||||
if (yBelowVx0 != yBelowVx1)
|
||||
{
|
||||
// check if X of point is not between X of Vx0 and Vx1
|
||||
if ( (xBelowVx0 = (vtx0[X] >= point[X])) == (vtx1[X] >= point[X]) )
|
||||
{
|
||||
if (xBelowVx0) crossings++;
|
||||
}
|
||||
else
|
||||
{
|
||||
// compute intersection of polygon segment with X ray, note if > point's X.
|
||||
crossings += (vtx0[X] - dv0*(vtx1[X] - vtx0[X])/(vtx1[Y] - vtx0[Y])) >= point[X];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// test if crossings is odd. If we care about its winding number > 0, then just: inside_flag = crossings > 0;
|
||||
if (crossings & 0x01) return true;
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// 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,
|
||||
std::vector<cvf::Vec3d>* createdVertexes,
|
||||
EdgeIntersectStorage<IndexType>* edgeIntersectionStorage,
|
||||
ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodes,
|
||||
const IndexType cv1CubeFaceIndices[4],
|
||||
const IndexType cv2CubeFaceIndices[4],
|
||||
double tolerance)
|
||||
{
|
||||
CVF_ASSERT(polygon);
|
||||
CVF_ASSERT(createdVertexes);
|
||||
|
||||
// Topology analysis
|
||||
|
||||
IndexType newVertexIndex = static_cast<IndexType>(nodes.size() + createdVertexes->size());
|
||||
|
||||
bool cv1VxTouchCv2[4] = { false, false, false, false };
|
||||
bool cv2VxTouchCv1[4] = { false, false, false, false };
|
||||
int cv1VxTouchCv2Edge[4] = { -1, -1, -1, -1 };
|
||||
int cv2VxTouchCv1Edge[4] = { -1, -1, -1, -1 };
|
||||
|
||||
int cv1Idx, cv2Idx;
|
||||
int numMatchedNodes = 0;
|
||||
|
||||
// First check for complete topological match.
|
||||
|
||||
for (cv1Idx = 0 ; cv1Idx < 4 ; ++cv1Idx)
|
||||
{
|
||||
bool found = false;
|
||||
for (cv2Idx = 0; cv2Idx < 4; ++cv2Idx)
|
||||
{
|
||||
if (cv1CubeFaceIndices[cv1Idx] == cv2CubeFaceIndices[cv2Idx])
|
||||
{
|
||||
cv1VxTouchCv2[cv1Idx] = true;
|
||||
cv2VxTouchCv1[cv2Idx] = true;
|
||||
found = true;
|
||||
++numMatchedNodes;
|
||||
continue;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (numMatchedNodes >= 4) // Todo: Handle collapsed cells
|
||||
{
|
||||
int k;
|
||||
for (k = 0; k < 4; ++k)
|
||||
{
|
||||
polygon->push_back(cv1CubeFaceIndices[k]);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
cvf::Vec3d diag1 = nodes[cv1CubeFaceIndices[2]] - nodes[cv1CubeFaceIndices[0]];
|
||||
cvf::Vec3d diag2 = nodes[cv1CubeFaceIndices[3]] - nodes[cv1CubeFaceIndices[1]];
|
||||
cvf::Vec3d normal = diag1^diag2;
|
||||
int numCv1VxesOnCv2 = numMatchedNodes;
|
||||
int numCv2VxesOnCv1 = numMatchedNodes;
|
||||
|
||||
for (cv1Idx = 0 ; cv1Idx < 4 ; ++cv1Idx)
|
||||
{
|
||||
if (!cv1VxTouchCv2[cv1Idx])
|
||||
{
|
||||
cv1VxTouchCv2[cv1Idx] = GeometryTools::isPointTouchingIndexedPolygon(normal,
|
||||
nodes,
|
||||
wrapArrayConst(cv2CubeFaceIndices, 4),
|
||||
nodes[cv1CubeFaceIndices[cv1Idx]],
|
||||
&(cv1VxTouchCv2Edge[cv1Idx]),
|
||||
tolerance);
|
||||
if (cv1VxTouchCv2[cv1Idx]) ++numCv1VxesOnCv2;
|
||||
}
|
||||
|
||||
if (!cv2VxTouchCv1[cv1Idx])
|
||||
{
|
||||
cv2VxTouchCv1[cv1Idx] = GeometryTools::isPointTouchingIndexedPolygon(normal,
|
||||
nodes,
|
||||
wrapArrayConst(cv1CubeFaceIndices, 4),
|
||||
nodes[cv2CubeFaceIndices[cv1Idx]],
|
||||
&(cv2VxTouchCv1Edge[cv1Idx]),
|
||||
tolerance);
|
||||
if (cv2VxTouchCv1[cv1Idx]) ++numCv2VxesOnCv1;
|
||||
}
|
||||
}
|
||||
|
||||
// Handle case where one of the faces are completely covered by the other
|
||||
|
||||
if (numCv1VxesOnCv2 >= 4)
|
||||
{
|
||||
int k;
|
||||
for (k = 0; k < 4; ++k)
|
||||
{
|
||||
polygon->push_back(cv1CubeFaceIndices[k]);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
if (numCv2VxesOnCv1 >= 4)
|
||||
{
|
||||
|
||||
int k;
|
||||
for (k = 0; k < 4; ++k)
|
||||
{
|
||||
polygon->push_back(cv2CubeFaceIndices[k]);
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
// Handle partial coverage
|
||||
// Algorithm outline as follows:
|
||||
|
||||
// Loop over edges in the face of Cv1. Intersect each one with all the edges of the Cv2 face.
|
||||
// Add first point of the cv1 edge to polygon if it really touches Cv2 ( touch of edge is considered as not touching)
|
||||
// Add each intersection point along the Cv1 edge if present
|
||||
// and finally: if the cv1 edge is going out of cv2, the add the cv2 vertexes from that intersection as long as they touch cv1.
|
||||
|
||||
int nextCv1Idx = 1;
|
||||
for (cv1Idx = 0 ; cv1Idx < 4 ; ++cv1Idx, ++nextCv1Idx)
|
||||
{
|
||||
if (nextCv1Idx > 3) nextCv1Idx = 0;
|
||||
|
||||
if (cv1VxTouchCv2[cv1Idx] && cv1VxTouchCv2Edge[cv1Idx] == -1) // Start of cv1 edge is touching inside the cv2 polygon (not on an cv2 edge)
|
||||
{
|
||||
if (polygon->empty() || polygon->back() != cv1CubeFaceIndices[cv1Idx])
|
||||
{
|
||||
polygon->push_back(cv1CubeFaceIndices[cv1Idx]);
|
||||
}
|
||||
|
||||
if (cv1VxTouchCv2[nextCv1Idx] && cv1VxTouchCv2Edge[nextCv1Idx] == -1)
|
||||
{
|
||||
// Both ends of this cv1 edge is touching inside(not on an edge) cv2 polygon, no intersections possible (assuming convex cube-face)
|
||||
// Continue with next Cv1-edge.
|
||||
continue;
|
||||
}
|
||||
}
|
||||
|
||||
// Find intersection(s) on this edge
|
||||
|
||||
std::vector<IndexType> intersectionVxIndices;
|
||||
std::vector<int> intersectedCv2EdgeIdxs;
|
||||
std::vector<double> intersectionFractionsAlongEdge;
|
||||
|
||||
int nextCv2Idx = 1;
|
||||
for (cv2Idx = 0; cv2Idx < 4; ++cv2Idx, ++nextCv2Idx)
|
||||
{
|
||||
if (nextCv2Idx > 3) nextCv2Idx = 0;
|
||||
|
||||
// Find a possible real intersection point.
|
||||
|
||||
cvf::Vec3d intersection(0,0,0);
|
||||
double fractionAlongEdge1;
|
||||
GeometryTools::IntersectionStatus intersectStatus = GeometryTools::NO_INTERSECTION;
|
||||
IndexType intersectionVxIndex = cvf::UNDEFINED_UINT;
|
||||
|
||||
// First handle some "trivial" ones to ease the burden for the real intersection calculation
|
||||
// It could be tested whether it really is necessary to do
|
||||
if (cv1VxTouchCv2Edge[cv1Idx] == cv2Idx && cv1VxTouchCv2Edge[nextCv1Idx] == cv2Idx )
|
||||
{
|
||||
intersectStatus = GeometryTools::LINES_OVERLAP;
|
||||
fractionAlongEdge1 = 1;
|
||||
intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
|
||||
}
|
||||
else if (cv1VxTouchCv2Edge[cv1Idx] == cv2Idx )
|
||||
{
|
||||
// When this happens, the cv1 vertex will already have been added to the polygon
|
||||
// by the statements in the top of the cv1 edge loop. Should it be treated specially ?
|
||||
intersectStatus = GeometryTools::LINES_TOUCH;
|
||||
fractionAlongEdge1 = 0;
|
||||
intersectionVxIndex = cv1CubeFaceIndices[cv1Idx];
|
||||
}
|
||||
else if (cv1VxTouchCv2Edge[nextCv1Idx] == cv2Idx )
|
||||
{
|
||||
intersectStatus = GeometryTools::LINES_TOUCH;
|
||||
fractionAlongEdge1 = 1;
|
||||
intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
|
||||
}
|
||||
else
|
||||
{
|
||||
double fractionAlongEdge2;
|
||||
bool found = false;
|
||||
if (edgeIntersectionStorage)
|
||||
found = edgeIntersectionStorage->findIntersection(
|
||||
cv1CubeFaceIndices[cv1Idx],
|
||||
cv1CubeFaceIndices[nextCv1Idx],
|
||||
cv2CubeFaceIndices[cv2Idx],
|
||||
cv2CubeFaceIndices[nextCv2Idx],
|
||||
&intersectionVxIndex, &intersectStatus,
|
||||
&fractionAlongEdge1, &fractionAlongEdge2);
|
||||
|
||||
if (!found)
|
||||
{
|
||||
|
||||
intersectStatus = GeometryTools::inPlaneLineIntersect3D(normal,
|
||||
nodes[cv1CubeFaceIndices[cv1Idx]],
|
||||
nodes[cv1CubeFaceIndices[nextCv1Idx]],
|
||||
nodes[cv2CubeFaceIndices[cv2Idx]],
|
||||
nodes[cv2CubeFaceIndices[nextCv2Idx]],
|
||||
&intersection, &fractionAlongEdge1, &fractionAlongEdge2,
|
||||
tolerance);
|
||||
|
||||
switch (intersectStatus)
|
||||
{
|
||||
case GeometryTools::LINES_CROSSES:
|
||||
{
|
||||
intersectionVxIndex = newVertexIndex;
|
||||
createdVertexes->push_back(intersection);
|
||||
++newVertexIndex;
|
||||
}
|
||||
break;
|
||||
case GeometryTools::LINES_TOUCH:
|
||||
{
|
||||
if (fractionAlongEdge1 <= 0.0) intersectionVxIndex = cv1CubeFaceIndices[cv1Idx];
|
||||
else if (fractionAlongEdge1 >= 1.0) intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
|
||||
else if (fractionAlongEdge2 <= 0.0) intersectionVxIndex = cv2CubeFaceIndices[cv2Idx];
|
||||
else if (fractionAlongEdge2 >= 1.0) intersectionVxIndex = cv2CubeFaceIndices[nextCv2Idx];
|
||||
else CVF_ASSERT(false); // Tolerance trouble
|
||||
}
|
||||
break;
|
||||
case GeometryTools::LINES_OVERLAP:
|
||||
{
|
||||
if (fractionAlongEdge1 <= 0.0) intersectionVxIndex = cv1CubeFaceIndices[cv1Idx];
|
||||
else if (fractionAlongEdge1 >= 1.0) intersectionVxIndex = cv1CubeFaceIndices[nextCv1Idx];
|
||||
else if (fractionAlongEdge2 <= 0.0) intersectionVxIndex = cv2CubeFaceIndices[cv2Idx];
|
||||
else if (fractionAlongEdge2 >= 1.0) intersectionVxIndex = cv2CubeFaceIndices[nextCv2Idx];
|
||||
else CVF_ASSERT(false); // Tolerance trouble
|
||||
}
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
if(edgeIntersectionStorage)
|
||||
edgeIntersectionStorage->addIntersection( cv1CubeFaceIndices[cv1Idx],
|
||||
cv1CubeFaceIndices[nextCv1Idx],
|
||||
cv2CubeFaceIndices[cv2Idx],
|
||||
cv2CubeFaceIndices[nextCv2Idx],
|
||||
intersectionVxIndex, intersectStatus,
|
||||
fractionAlongEdge1, fractionAlongEdge2);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
// Store data for each intersection along the Cv1-edge
|
||||
|
||||
if ( (intersectStatus == GeometryTools::LINES_CROSSES)
|
||||
|| (intersectStatus == GeometryTools::LINES_TOUCH)
|
||||
|| (intersectStatus == GeometryTools::LINES_OVERLAP) )
|
||||
{
|
||||
CVF_ASSERT(intersectionVxIndex != cvf::UNDEFINED_UINT);
|
||||
|
||||
intersectionFractionsAlongEdge.push_back(fractionAlongEdge1);
|
||||
intersectedCv2EdgeIdxs.push_back(cv2Idx);
|
||||
intersectionVxIndices.push_back(intersectionVxIndex);
|
||||
}
|
||||
}
|
||||
|
||||
// Insert the intersections into the polygon in the correct order along the Cv1 edge.
|
||||
// Find the last intersection in order to possibly continue the polygon along Cv2 into Cv1
|
||||
|
||||
size_t i;
|
||||
size_t lastIntersection = std::numeric_limits<size_t>::max();
|
||||
double largestFraction = -1;
|
||||
for (i = 0; i < intersectionFractionsAlongEdge.size(); ++i)
|
||||
{
|
||||
if (intersectionFractionsAlongEdge[i] > largestFraction)
|
||||
{
|
||||
lastIntersection = i;
|
||||
largestFraction = intersectionFractionsAlongEdge[i];
|
||||
}
|
||||
}
|
||||
|
||||
// Insert indices to the new intersection vertices into the polygon of
|
||||
// this connection according to fraction along edge
|
||||
|
||||
std::map<double,IndexType> sortingMap;
|
||||
for (i = 0; i < intersectionFractionsAlongEdge.size(); ++i)
|
||||
{
|
||||
sortingMap[intersectionFractionsAlongEdge[i]] = intersectionVxIndices[i];
|
||||
}
|
||||
|
||||
typename std::map<double, IndexType>::iterator it;
|
||||
for (it = sortingMap.begin(); it != sortingMap.end(); ++it)
|
||||
{
|
||||
if (polygon->empty() || polygon->back() != it->second)
|
||||
{
|
||||
polygon->push_back(it->second);
|
||||
}
|
||||
}
|
||||
|
||||
// Then if the Cv1 edge is going out of Cv2, add to the polygon, all the Cv2 face vertex-indices
|
||||
// from the intersection that touches Cv1.
|
||||
|
||||
// if cv1 edge in any way touches cv2 and ends up outside, it went out.
|
||||
|
||||
/*
|
||||
if cv1 edge is going out of cv2 then
|
||||
if intersected cv2 edge has endpoint touching cv1 then
|
||||
add endpoint to polygon. continue to add next endpoint until it does not touch Cv1
|
||||
*/
|
||||
if ( !cv1VxTouchCv2[nextCv1Idx]
|
||||
&& ( cv1VxTouchCv2[cv1Idx] || ( intersectedCv2EdgeIdxs.size() ) ) ) // Two touches along edge also qualifies
|
||||
{
|
||||
if(lastIntersection < intersectedCv2EdgeIdxs.size())
|
||||
{
|
||||
cv2Idx = intersectedCv2EdgeIdxs[lastIntersection];
|
||||
int count = 0;
|
||||
// Continue the polygon along the Cv2 edges as long as they touch cv1.
|
||||
// Depending on the faces having opposite winding, which is guaranteed as long as
|
||||
// no intersecting CVs share a connection
|
||||
while (cv2VxTouchCv1[cv2Idx] && count < 4 && (cv2VxTouchCv1Edge[cv2Idx] == -1)) // Touch of edge is regarded as being outside, so we must stop
|
||||
{
|
||||
if (polygon->empty() || polygon->back() != cv2CubeFaceIndices[cv2Idx])
|
||||
{
|
||||
polygon->push_back(cv2CubeFaceIndices[cv2Idx]);
|
||||
}
|
||||
--cv2Idx;
|
||||
if (cv2Idx < 0 ) cv2Idx = 3;
|
||||
++count;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
CVF_ASSERT(lastIntersection < intersectedCv2EdgeIdxs.size());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (polygon->size() > 2)
|
||||
{
|
||||
if (polygon->back() == polygon->front()) polygon->pop_back();
|
||||
}
|
||||
|
||||
// Sanity checks
|
||||
if (polygon->size() < 3)
|
||||
{
|
||||
// cvf::Trace::show(cvf::String("Degenerated connection polygon detected. (Less than 3 vertexes) Cv's probably not in contact: %1 , %2").arg(m_ownerCvId).arg(m_neighborCvId));
|
||||
polygon->clear();
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// This method assumes that all intersection and mid edge vertexes are created an are already
|
||||
/// merged into all the polygons. We can also assume that all the connection polygons are completely
|
||||
/// inside (or sharing edges with) the cube face polygon initially
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
#define DEBUG_PRINT 0
|
||||
|
||||
//template <typename NodeArrayType, typename NodeType, typename IndicesArrayType, typename IndicesType>
|
||||
//void setup( ArrayWrapper<NodeArrayType, NodeType> nodeArray, ArrayWrapper<IndicesArrayType, IndicesType> indices)
|
||||
template<typename VerticeArrayType, typename PolygonArrayType, typename IndexType>
|
||||
void GeometryTools::calculatePartiallyFreeCubeFacePolygon(ArrayWrapperConst<VerticeArrayType, cvf::Vec3d> nodeCoords,
|
||||
ArrayWrapperConst<PolygonArrayType, IndexType> completeFacePolygon,
|
||||
const cvf::Vec3d& faceNormal,
|
||||
const std::vector< std::vector<IndexType>* >& faceOverlapPolygons,
|
||||
const std::vector<bool> faceOverlapPolygonWindingSameAsCubeFaceFlags,
|
||||
std::vector<IndexType>* partialFacePolygon,
|
||||
bool* m_partiallyFreeCubeFaceHasHoles)
|
||||
{
|
||||
// Vertex Index to position in polygon
|
||||
typedef std::map< IndexType, typename std::vector<IndexType>::const_iterator > VxIdxToPolygonPositionMap;
|
||||
|
||||
CVF_ASSERT(m_partiallyFreeCubeFaceHasHoles);
|
||||
CVF_ASSERT(partialFacePolygon != NULL);
|
||||
|
||||
// Copy the start polygon
|
||||
std::list<IndexType> resultPolygon;
|
||||
for (size_t pcIdx = 0; pcIdx < completeFacePolygon.size(); ++pcIdx)
|
||||
{
|
||||
resultPolygon.push_back(completeFacePolygon[pcIdx]);
|
||||
}
|
||||
|
||||
// First build search maps to fast find whether and where an index is positioned in a polygon
|
||||
// Map from Vertex-index to position in polygon
|
||||
|
||||
std::vector< VxIdxToPolygonPositionMap > polygonSearchMaps;
|
||||
std::vector<bool> isConnectionPolygonMerged;
|
||||
|
||||
polygonSearchMaps.resize(faceOverlapPolygons.size());
|
||||
isConnectionPolygonMerged.resize(faceOverlapPolygons.size(), false);
|
||||
|
||||
// Build search maps
|
||||
{
|
||||
size_t count;
|
||||
for (size_t i = 0; i < faceOverlapPolygons.size(); ++i)
|
||||
{
|
||||
count = 0;
|
||||
for (typename std::vector<IndexType >::const_iterator pcIt = faceOverlapPolygons[i]->begin();
|
||||
pcIt != faceOverlapPolygons[i]->end();
|
||||
++pcIt)
|
||||
{
|
||||
polygonSearchMaps[i][*pcIt] = pcIt;
|
||||
++count;
|
||||
}
|
||||
|
||||
if (count < 3) isConnectionPolygonMerged[i] = true; // Ignore false polygons
|
||||
}
|
||||
}
|
||||
|
||||
#if DEBUG_PRINT
|
||||
{
|
||||
cvf::Trace::show("Circumference polygon: ");
|
||||
std::list<IndexType>::const_iterator polIt;
|
||||
for ( polIt = resultPolygon.begin(); polIt != resultPolygon.end(); ++polIt)
|
||||
{
|
||||
cvf::Trace::show(cvf::String("%1 \t%2 %3 %4").arg((int)(*polIt))
|
||||
.arg(nodeCoords[*polIt].x())
|
||||
.arg(nodeCoords[*polIt].y())
|
||||
.arg(nodeCoords[*polIt].z()));
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
#if DEBUG_PRINT
|
||||
{
|
||||
cvf::Trace::show("Connection polygons: ");
|
||||
for (size_t cIdx = 0; cIdx < faceOverlapPolygons.size(); cIdx++)
|
||||
{
|
||||
std::vector<IndexType >::const_iterator polIt;
|
||||
cvf::Trace::show("Connection " + cvf::String((long long)cIdx));
|
||||
for (polIt = faceOverlapPolygons[cIdx]->begin(); polIt != faceOverlapPolygons[cIdx]->end(); ++polIt)
|
||||
{
|
||||
cvf::Trace::show(cvf::String("%1 \t%2 %3 %4").arg((int)(*polIt))
|
||||
.arg(nodeCoords[*polIt].x())
|
||||
.arg(nodeCoords[*polIt].y())
|
||||
.arg(nodeCoords[*polIt].z()));
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
// Merge connection polygons with the main polygon as long as one of them has something in common.
|
||||
|
||||
// For each vx in the m_freeFacePolygons[cubeFace] polygon .
|
||||
// loop over all connections
|
||||
// if it has the node in common and that the edge angle will decrease if inserting
|
||||
// merge the connection polygon into the main polygon,
|
||||
// and remove the connection polygon from the merge able connection polygons.
|
||||
|
||||
|
||||
for (typename std::list<IndexType>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end(); ++pIt)
|
||||
{
|
||||
// Set iterator to previous node in polygon
|
||||
typename std::list<IndexType>::iterator prevPIt = pIt;
|
||||
if (prevPIt == resultPolygon.begin()) prevPIt = resultPolygon.end();
|
||||
--prevPIt;
|
||||
|
||||
cvf::Vec3d pToPrev = nodeCoords[*prevPIt] - nodeCoords[*pIt];
|
||||
|
||||
// Set iterator to next node in polygon. Used to insert before and as pointer to the next point
|
||||
typename std::list<IndexType>::iterator nextPIt = pIt;
|
||||
++nextPIt;
|
||||
typename std::list<IndexType>::iterator insertBeforePIt = nextPIt;
|
||||
if (nextPIt == resultPolygon.end()) nextPIt = resultPolygon.begin();
|
||||
|
||||
// Calculate existing edge to edge angle
|
||||
|
||||
cvf::Vec3d pToNext = nodeCoords[*nextPIt] - nodeCoords[*pIt];
|
||||
double mainPolygonEdgeAngle = GeometryTools::getAngle(faceNormal, pToNext , pToPrev);
|
||||
|
||||
// Find connections containing the pIt vertex index. Merge them into the main polygon
|
||||
|
||||
for (size_t opIdx = 0; opIdx < faceOverlapPolygons.size(); ++opIdx)
|
||||
{
|
||||
if (isConnectionPolygonMerged[opIdx]) continue; // Already merged
|
||||
|
||||
// Find position of pIt vertex index in the current connection polygon
|
||||
typename VxIdxToPolygonPositionMap::iterator vxIndexPositionInPolygonIt = polygonSearchMaps[opIdx].find(*pIt);
|
||||
|
||||
if (vxIndexPositionInPolygonIt != polygonSearchMaps[opIdx].end())
|
||||
{
|
||||
// Merge the connection polygon into the main polygon
|
||||
// if the angle prevPIt pIt nextPIt is larger than angle prevPIt pIt (startCPIt++)
|
||||
|
||||
typename std::vector<IndexType>::const_iterator startCPIt;
|
||||
startCPIt = vxIndexPositionInPolygonIt->second;
|
||||
|
||||
// First vx to insert is the one after the match
|
||||
|
||||
bool hasSameWinding = faceOverlapPolygonWindingSameAsCubeFaceFlags[opIdx];
|
||||
if (hasSameWinding)
|
||||
{
|
||||
// Same winding as main polygon. We need to go the opposite way
|
||||
if (startCPIt == faceOverlapPolygons[opIdx]->begin()) startCPIt = faceOverlapPolygons[opIdx]->end();
|
||||
--startCPIt;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Opposite winding. Go forward when merging
|
||||
++startCPIt; if (startCPIt == faceOverlapPolygons[opIdx]->end()) startCPIt = faceOverlapPolygons[opIdx]->begin();
|
||||
}
|
||||
|
||||
// Calculate possible new edge-to-edge angle and test against existing angle
|
||||
cvf::Vec3d pToStart = nodeCoords[*startCPIt] - nodeCoords[*pIt];
|
||||
double candidatePolygonEdgeAngle = GeometryTools::getAngle(faceNormal, pToStart , pToPrev);
|
||||
|
||||
if (candidatePolygonEdgeAngle < mainPolygonEdgeAngle )
|
||||
{
|
||||
// Merge ok
|
||||
typename std::vector<IndexType >::const_iterator pcIt = startCPIt;
|
||||
if (hasSameWinding)
|
||||
{
|
||||
do
|
||||
{
|
||||
resultPolygon.insert(insertBeforePIt, (*pcIt));
|
||||
|
||||
if (pcIt == faceOverlapPolygons[opIdx]->begin()) pcIt = faceOverlapPolygons[opIdx]->end();
|
||||
--pcIt;
|
||||
|
||||
} while (pcIt != startCPIt);
|
||||
}
|
||||
else
|
||||
{
|
||||
do
|
||||
{
|
||||
resultPolygon.insert(insertBeforePIt, (*pcIt));
|
||||
|
||||
++pcIt;
|
||||
if (pcIt == faceOverlapPolygons[opIdx]->end()) pcIt = faceOverlapPolygons[opIdx]->begin();
|
||||
|
||||
} while (pcIt != startCPIt);
|
||||
}
|
||||
|
||||
isConnectionPolygonMerged[opIdx] = true;
|
||||
|
||||
// Recalculate the next position to point into the new nodes
|
||||
// Set iterator in the main polygon to insert before and to the next point
|
||||
nextPIt = pIt;
|
||||
++nextPIt;
|
||||
insertBeforePIt = nextPIt;
|
||||
if (nextPIt == resultPolygon.end()) nextPIt = resultPolygon.begin();
|
||||
|
||||
// Recalculate the existing edge to edge angle
|
||||
pToNext = nodeCoords[*nextPIt] - nodeCoords[*pIt];
|
||||
mainPolygonEdgeAngle = GeometryTools::getAngle(faceNormal, pToNext , pToPrev);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Now remove all double edges
|
||||
|
||||
bool goneAround = false;
|
||||
for ( typename std::list<IndexType>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end() && !goneAround; ++pIt)
|
||||
{
|
||||
// Set iterator to next node in polygon.
|
||||
typename std::list<IndexType>::iterator nextPIt = pIt;
|
||||
++nextPIt;
|
||||
if (nextPIt == resultPolygon.end())
|
||||
{
|
||||
nextPIt = resultPolygon.begin();
|
||||
goneAround = true; // Gone around polygon. Stop even if pIt is jumping over end()
|
||||
}
|
||||
|
||||
// Set iterator to previous node in polygon
|
||||
|
||||
typename std::list<IndexType>::iterator prevPIt = pIt;
|
||||
|
||||
if (prevPIt == resultPolygon.begin()) prevPIt = resultPolygon.end();
|
||||
--prevPIt;
|
||||
|
||||
// If previous and next node are the same, erase
|
||||
while(*nextPIt == *prevPIt)
|
||||
{
|
||||
resultPolygon.erase(pIt);
|
||||
resultPolygon.erase(prevPIt);
|
||||
|
||||
if ( resultPolygon.begin() == resultPolygon.end()) break; // Polygon has been completely removed. Nothing left. Break out of while
|
||||
|
||||
pIt = nextPIt;
|
||||
++nextPIt;
|
||||
if (nextPIt == resultPolygon.end())
|
||||
{
|
||||
nextPIt = resultPolygon.begin();
|
||||
goneAround = true; // Gone around polygon pIt is jumping over end()
|
||||
}
|
||||
|
||||
prevPIt = pIt;
|
||||
if (prevPIt == resultPolygon.begin()) prevPIt = resultPolygon.end();
|
||||
--prevPIt;
|
||||
}
|
||||
|
||||
if ( resultPolygon.begin() == resultPolygon.end()) break; // Polygon has been completely removed. Nothing left. Break out of for loop
|
||||
|
||||
}
|
||||
|
||||
// Check for holes
|
||||
|
||||
bool hasHoles = false;
|
||||
for (size_t i = 0; i < isConnectionPolygonMerged.size(); ++i)
|
||||
{
|
||||
hasHoles = !isConnectionPolygonMerged[i];
|
||||
if(hasHoles)
|
||||
{
|
||||
*m_partiallyFreeCubeFaceHasHoles = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
#if DEBUG_PRINT
|
||||
{
|
||||
cvf::Trace::show("Polygon: ");
|
||||
for (std::list<IndexType>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end(); ++pIt)
|
||||
{
|
||||
cvf::Trace::show(cvf::String("%1 \t%2 %3 %4").arg((int)(*pIt))
|
||||
.arg(nodeCoords[*pIt].x())
|
||||
.arg(nodeCoords[*pIt].y())
|
||||
.arg(nodeCoords[*pIt].z()));
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
// Copy the result polygon to the output variable
|
||||
|
||||
partialFacePolygon->clear();
|
||||
for (typename std::list<IndexType>::iterator pIt = resultPolygon.begin(); pIt != resultPolygon.end(); ++pIt)
|
||||
{
|
||||
partialFacePolygon->push_back(*pIt);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
template <typename IndexType>
|
||||
void EdgeIntersectStorage<IndexType>::setVertexCount(size_t size)
|
||||
{
|
||||
m_edgeIntsectMap.resize(size);
|
||||
}
|
||||
|
||||
template <typename IndexType>
|
||||
void EdgeIntersectStorage<IndexType>::canonizeAddress(IndexType& e1P1, IndexType& e1P2, IndexType& e2P1, IndexType& e2P2,
|
||||
bool& flipE1, bool& flipE2, bool& flipE1E2)
|
||||
{
|
||||
flipE1 = e1P1 > e1P2;
|
||||
flipE2 = e2P1 > e2P2;
|
||||
|
||||
flipE1E2 = (flipE1 ? e1P2: e1P1) > (flipE2 ? e2P2: e2P1);
|
||||
|
||||
static IndexType temp;
|
||||
if (flipE1)
|
||||
{
|
||||
temp = e1P1;
|
||||
e1P1 = e1P2;
|
||||
e1P2 = temp;
|
||||
}
|
||||
|
||||
if (flipE2)
|
||||
{
|
||||
temp = e2P1;
|
||||
e2P1 = e2P2;
|
||||
e2P2 = temp;
|
||||
}
|
||||
|
||||
if (flipE1E2)
|
||||
{
|
||||
temp = e1P1;
|
||||
e1P1 = e2P1;
|
||||
e2P1 = temp;
|
||||
|
||||
temp = e1P2;
|
||||
e1P2 = e2P2;
|
||||
e2P2 = temp;
|
||||
}
|
||||
}
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
template <typename IndexType>
|
||||
void EdgeIntersectStorage<IndexType>::addIntersection(IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2,
|
||||
IndexType vxIndexIntersectionPoint, GeometryTools::IntersectionStatus intersectionStatus,
|
||||
double fractionAlongEdge1, double fractionAlongEdge2)
|
||||
{
|
||||
static bool flipE1 ;
|
||||
static bool flipE2 ;
|
||||
static bool flipE1E2;
|
||||
|
||||
canonizeAddress(e1P1, e1P2, e2P1, e2P2, flipE1, flipE2, flipE1E2);
|
||||
|
||||
static IntersectData iData;
|
||||
|
||||
iData.fractionAlongEdge1 = flipE1 ? 1 - fractionAlongEdge1 : fractionAlongEdge1;
|
||||
iData.fractionAlongEdge2 = flipE2 ? 1 - fractionAlongEdge2 : fractionAlongEdge2;
|
||||
iData.intersectionStatus = intersectionStatus;
|
||||
|
||||
if (flipE1E2)
|
||||
{
|
||||
double temp = iData.fractionAlongEdge1;
|
||||
iData.fractionAlongEdge1 = iData.fractionAlongEdge2;
|
||||
iData.fractionAlongEdge2 = temp;
|
||||
}
|
||||
|
||||
iData.intersectionPointIndex = vxIndexIntersectionPoint;
|
||||
CVF_ASSERT(e1P1 < m_edgeIntsectMap.size());
|
||||
m_edgeIntsectMap[e1P1][e1P2][e2P1][e2P2] = iData;
|
||||
}
|
||||
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
template <typename IndexType>
|
||||
bool EdgeIntersectStorage<IndexType>::findIntersection(IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2,
|
||||
IndexType* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus,
|
||||
double* fractionAlongEdge1, double* fractionAlongEdge2)
|
||||
{
|
||||
static bool flipE1 ;
|
||||
static bool flipE2 ;
|
||||
static bool flipE1E2;
|
||||
|
||||
canonizeAddress(e1P1, e1P2, e2P1, e2P2, flipE1, flipE2, flipE1E2);
|
||||
|
||||
if (!m_edgeIntsectMap[e1P1].size()) return false;
|
||||
|
||||
typename 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;
|
||||
|
||||
typename std::map<IndexType, std::map<IndexType, IntersectData > >::iterator it2;
|
||||
it2 = it->second.find(e2P1);
|
||||
if (it2 == it->second.end()) return false;
|
||||
|
||||
typename std::map<IndexType, IntersectData >::iterator it3;
|
||||
it3 = it2->second.find(e2P2);
|
||||
if (it3 == it2->second.end()) return false;
|
||||
|
||||
*vxIndexIntersectionPoint = it3->second.intersectionPointIndex;
|
||||
*intersectionStatus = it3->second.intersectionStatus;
|
||||
|
||||
if (flipE1E2)
|
||||
{
|
||||
*fractionAlongEdge1 = it3->second.fractionAlongEdge2;
|
||||
*fractionAlongEdge2 = it3->second.fractionAlongEdge1;
|
||||
}
|
||||
else
|
||||
{
|
||||
*fractionAlongEdge1 = it3->second.fractionAlongEdge1;
|
||||
*fractionAlongEdge2 = it3->second.fractionAlongEdge2;
|
||||
}
|
||||
|
||||
if (flipE1) *fractionAlongEdge1 = 1 - *fractionAlongEdge1;
|
||||
if (flipE2) *fractionAlongEdge2 = 1 - *fractionAlongEdge2;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
Reference in New Issue
Block a user