diff --git a/common/gal/graphics_abstraction_layer.cpp b/common/gal/graphics_abstraction_layer.cpp index 0282e69327..70e761f5c9 100644 --- a/common/gal/graphics_abstraction_layer.cpp +++ b/common/gal/graphics_abstraction_layer.cpp @@ -54,8 +54,8 @@ GAL::GAL() : // Initialize the cursor shape SetCursorColor( COLOR4D( 1.0, 1.0, 1.0, 1.0 ) ); - SetCursorSize( 15 ); - SetCursorEnabled( true ); + SetCursorSize( 80 ); + SetCursorEnabled( false ); strokeFont.LoadNewStrokeFont( newstroke_font, newstroke_font_bufsize ); } diff --git a/common/gal/opengl/opengl_gal.cpp b/common/gal/opengl/opengl_gal.cpp index 78cb44962e..e1859d9869 100644 --- a/common/gal/opengl/opengl_gal.cpp +++ b/common/gal/opengl/opengl_gal.cpp @@ -82,7 +82,7 @@ OPENGL_GAL::OPENGL_GAL( wxWindow* aParent, wxEvtHandler* aMouseListener, SetSize( aParent->GetSize() ); screenSize = VECTOR2D( aParent->GetSize() ); - initCursor( 20 ); + initCursor( 80 ); // Grid color settings are different in Cairo and OpenGL SetGridColor( COLOR4D( 0.8, 0.8, 0.8, 0.1 ) ); diff --git a/common/geometry/shape_line_chain.cpp b/common/geometry/shape_line_chain.cpp index 44c7bad7f4..6edc0d2b14 100644 --- a/common/geometry/shape_line_chain.cpp +++ b/common/geometry/shape_line_chain.cpp @@ -439,7 +439,7 @@ SHAPE_LINE_CHAIN& SHAPE_LINE_CHAIN::Simplify() const VECTOR2I SHAPE_LINE_CHAIN::NearestPoint(const VECTOR2I& aP) const { int min_d = INT_MAX; - int nearest; + int nearest = 0; for ( int i = 0; i < SegmentCount() ; i++ ) { int d = CSegment(i).Distance(aP); diff --git a/common/profile.h b/common/profile.h index 82ecdf199e..f779dfae59 100644 --- a/common/profile.h +++ b/common/profile.h @@ -27,8 +27,8 @@ * @brief Simple profiling functions for measuring code execution time. */ -#ifndef __PROFILE_H -#define __PROFILE_H +#ifndef __TPROFILE_H +#define __TPROFILE_H #include #include @@ -96,7 +96,6 @@ static inline uint64_t get_tics() gettimeofday( &tv, NULL ); return (uint64_t) tv.tv_sec * 1000000ULL + (uint64_t) tv.tv_usec; -#endif } @@ -144,3 +143,5 @@ static inline void prof_end( prof_counter* cnt ) else cnt->value = get_tics() - cnt->value; } + +#endif diff --git a/common/view/view.cpp b/common/view/view.cpp index 1aa45473d2..529b966d49 100644 --- a/common/view/view.cpp +++ b/common/view/view.cpp @@ -452,6 +452,14 @@ void VIEW::ChangeLayerDepth( int aLayer, int aDepth ) m_layers[aLayer].items->Query( r, visitor ); } +int VIEW::GetTopLayer( ) const +{ + if( m_topLayers.size() == 0 ) + return 0; + + return *m_topLayers.begin(); +} + void VIEW::SetTopLayer( int aLayer, bool aEnabled ) { diff --git a/common/view/wx_view_controls.cpp b/common/view/wx_view_controls.cpp index 9776710676..4aa69d2d61 100644 --- a/common/view/wx_view_controls.cpp +++ b/common/view/wx_view_controls.cpp @@ -60,11 +60,24 @@ WX_VIEW_CONTROLS::WX_VIEW_CONTROLS( VIEW* aView, wxWindow* aParentPanel ) : } +void VIEW_CONTROLS::ShowCursor( bool aEnabled ) +{ + m_view->GetGAL()->SetCursorEnabled( aEnabled ); +} + + void WX_VIEW_CONTROLS::onMotion( wxMouseEvent& aEvent ) { m_mousePosition.x = aEvent.GetX(); m_mousePosition.y = aEvent.GetY(); + if( m_forceCursorPosition ) + m_cursorPosition = m_view->ToScreen( m_forcedPosition ); + else if( m_snappingEnabled ) + m_cursorPosition = m_view->GetGAL()->GetGridPoint( m_mousePosition ); + else + m_cursorPosition = m_mousePosition; + bool isAutoPanning = false; if( m_autoPanEnabled ) diff --git a/include/geometry/rtree.h b/include/geometry/rtree.h new file mode 100644 index 0000000000..d64d7a074a --- /dev/null +++ b/include/geometry/rtree.h @@ -0,0 +1,1908 @@ +//TITLE +// +// R-TREES: A DYNAMIC INDEX STRUCTURE FOR SPATIAL SEARCHING +// +//DESCRIPTION +// +// A C++ templated version of the RTree algorithm. +// For more information please read the comments in RTree.h +// +//AUTHORS +// +// * 1983 Original algorithm and test code by Antonin Guttman and Michael Stonebraker, UC Berkely +// * 1994 ANCI C ported from original test code by Melinda Green - melinda@superliminal.com +// * 1995 Sphere volume fix for degeneracy problem submitted by Paul Brook +// * 2004 Templated C++ port by Greg Douglas +// * 2013 CERN (www.cern.ch) +// +//LICENSE: +// +// Entirely free for all uses. Enjoy! + +#ifndef RTREE_H +#define RTREE_H + +// NOTE This file compiles under MSVC 6 SP5 and MSVC .Net 2003 it may not work on other compilers without modification. + +// NOTE These next few lines may be win32 specific, you may need to modify them to compile on other platform +#include +#include +#include +#include + +#define ASSERT assert // RTree uses ASSERT( condition ) +#ifndef Min + #define Min std::min +#endif // Min +#ifndef Max + #define Max std::max +#endif // Max + +// +// RTree.h +// + +#define RTREE_TEMPLATE template +#define RTREE_SEARCH_TEMPLATE template +#define RTREE_QUAL RTree +#define RTREE_SEARCH_QUAL RTree + +#define RTREE_DONT_USE_MEMPOOLS // This version does not contain a fixed memory allocator, fill in lines with EXAMPLE to implement one. +#define RTREE_USE_SPHERICAL_VOLUME // Better split classification, may be slower on some systems + +// Fwd decl +class RTFileStream; // File I/O helper class, look below for implementation and notes. + + +/// \class RTree +/// Implementation of RTree, a multidimensional bounding rectangle tree. +/// Example usage: For a 3-dimensional tree use RTree myTree; +/// +/// This modified, templated C++ version by Greg Douglas at Auran (http://www.auran.com) +/// +/// DATATYPE Referenced data, should be int, void*, obj* etc. no larger than sizeof and simple type +/// ELEMTYPE Type of element such as int or float +/// NUMDIMS Number of dimensions such as 2 or 3 +/// ELEMTYPEREAL Type of element that allows fractional and large values such as float or double, for use in volume calcs +/// +/// NOTES: Inserting and removing data requires the knowledge of its constant Minimal Bounding Rectangle. +/// This version uses new/delete for nodes, I recommend using a fixed size allocator for efficiency. +/// Instead of using a callback function for returned results, I recommend and efficient pre-sized, grow-only memory +/// array similar to MFC CArray or STL Vector for returning search query result. +/// +template +class RTree +{ +protected: + + struct Node; // Fwd decl. Used by other internal structs and iterator +public: + + // These constant must be declared after Branch and before Node struct + // Stuck up here for MSVC 6 compiler. NSVC .NET 2003 is much happier. + enum { + MAXNODES = TMAXNODES, ///< Max elements in node + MINNODES = TMINNODES, ///< Min elements in node + }; + + struct Statistics { + int maxDepth; + int avgDepth; + + int maxNodeLoad; + int avgNodeLoad; + int totalItems; + }; + +public: + + RTree(); + virtual ~RTree(); + + /// Insert entry + /// \param a_min Min of bounding rect + /// \param a_max Max of bounding rect + /// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not allowed. + void Insert( const ELEMTYPE a_min[NUMDIMS], + const ELEMTYPE a_max[NUMDIMS], + const DATATYPE& a_dataId ); + + /// Remove entry + /// \param a_min Min of bounding rect + /// \param a_max Max of bounding rect + /// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not allowed. + void Remove( const ELEMTYPE a_min[NUMDIMS], + const ELEMTYPE a_max[NUMDIMS], + const DATATYPE& a_dataId ); + + /// Find all within search rectangle + /// \param a_min Min of search bounding rect + /// \param a_max Max of search bounding rect + /// \param a_searchResult Search result array. Caller should set grow size. Function will reset, not append to array. + /// \param a_resultCallback Callback function to return result. Callback should return 'true' to continue searching + /// \param a_context User context to pass as parameter to a_resultCallback + /// \return Returns the number of entries found + int Search( const ELEMTYPE a_min[NUMDIMS], + const ELEMTYPE a_max[NUMDIMS], + bool a_resultCallback( DATATYPE a_data, void* a_context ), + void* a_context ); + + template + int Search( const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], VISITOR& a_visitor ) + { + #ifdef _DEBUG + + for( int index = 0; index 0; } + + /// Access the current data element. Caller must be sure iterator is not NULL first. + DATATYPE& operator*() + { + ASSERT( IsNotNull() ); + StackElement& curTos = m_stack[m_tos - 1]; + return curTos.m_node->m_branch[curTos.m_branchIndex].m_data; + } + + /// Access the current data element. Caller must be sure iterator is not NULL first. + const DATATYPE& operator*() const + { + ASSERT( IsNotNull() ); + StackElement& curTos = m_stack[m_tos - 1]; + return curTos.m_node->m_branch[curTos.m_branchIndex].m_data; + } + + /// Find the next data element + bool operator++() { return FindNextData(); } + + /// Get the bounds for this node + void GetBounds( ELEMTYPE a_min[NUMDIMS], ELEMTYPE a_max[NUMDIMS] ) + { + ASSERT( IsNotNull() ); + StackElement& curTos = m_stack[m_tos - 1]; + Branch& curBranch = curTos.m_node->m_branch[curTos.m_branchIndex]; + + for( int index = 0; index < NUMDIMS; ++index ) + { + a_min[index] = curBranch.m_rect.m_min[index]; + a_max[index] = curBranch.m_rect.m_max[index]; + } + } + + private: + + /// Reset iterator + void Init() { m_tos = 0; } + + /// Find the next data element in the tree (For internal use only) + bool FindNextData() + { + for( ; ; ) + { + if( m_tos <= 0 ) + { + return false; + } + + StackElement curTos = Pop(); // Copy stack top cause it may change as we use it + + if( curTos.m_node->IsLeaf() ) + { + // Keep walking through data while we can + if( curTos.m_branchIndex + 1 < curTos.m_node->m_count ) + { + // There is more data, just point to the next one + Push( curTos.m_node, curTos.m_branchIndex + 1 ); + return true; + } + + // No more data, so it will fall back to previous level + } + else + { + if( curTos.m_branchIndex + 1 < curTos.m_node->m_count ) + { + // Push sibling on for future tree walk + // This is the 'fall back' node when we finish with the current level + Push( curTos.m_node, curTos.m_branchIndex + 1 ); + } + + // Since cur node is not a leaf, push first of next level to get deeper into the tree + Node* nextLevelnode = curTos.m_node->m_branch[curTos.m_branchIndex].m_child; + Push( nextLevelnode, 0 ); + + // If we pushed on a new leaf, exit as the data is ready at TOS + if( nextLevelnode->IsLeaf() ) + { + return true; + } + } + } + } + + /// Push node and branch onto iteration stack (For internal use only) + void Push( Node* a_node, int a_branchIndex ) + { + m_stack[m_tos].m_node = a_node; + m_stack[m_tos].m_branchIndex = a_branchIndex; + ++m_tos; + ASSERT( m_tos <= MAX_STACK ); + } + + /// Pop element off iteration stack (For internal use only) + StackElement& Pop() + { + ASSERT( m_tos > 0 ); + --m_tos; + return m_stack[m_tos]; + } + + StackElement m_stack[MAX_STACK]; ///< Stack as we are doing iteration instead of recursion + int m_tos; ///< Top Of Stack index + + friend class RTree; // Allow hiding of non-public functions while allowing manipulation by logical owner + }; + + + /// Get 'first' for iteration + void GetFirst( Iterator& a_it ) + { + a_it.Init(); + Node* first = m_root; + + while( first ) + { + if( first->IsInternalNode() && first->m_count > 1 ) + { + a_it.Push( first, 1 ); // Descend sibling branch later + } + else if( first->IsLeaf() ) + { + if( first->m_count ) + { + a_it.Push( first, 0 ); + } + + break; + } + + first = first->m_branch[0].m_child; + } + } + + /// Get Next for iteration + void GetNext( Iterator& a_it ) { ++a_it; } + + /// Is iterator NULL, or at end? + bool IsNull( Iterator& a_it ) { return a_it.IsNull(); } + + /// Get object at iterator position + DATATYPE& GetAt( Iterator& a_it ) { return *a_it; } +protected: + + /// Minimal bounding rectangle (n-dimensional) + struct Rect + { + ELEMTYPE m_min[NUMDIMS]; ///< Min dimensions of bounding box + ELEMTYPE m_max[NUMDIMS]; ///< Max dimensions of bounding box + }; + + /// May be data or may be another subtree + /// The parents level determines this. + /// If the parents level is 0, then this is data + struct Branch + { + Rect m_rect; ///< Bounds + union + { + Node* m_child; ///< Child node + DATATYPE m_data; ///< Data Id or Ptr + }; + }; + + /// Node for each branch level + struct Node + { + bool IsInternalNode() { return m_level > 0; } // Not a leaf, but a internal node + bool IsLeaf() { return m_level == 0; } // A leaf, contains data + + int m_count; ///< Count + int m_level; ///< Leaf is zero, others positive + Branch m_branch[MAXNODES]; ///< Branch + }; + + /// A link list of nodes for reinsertion after a delete operation + struct ListNode + { + ListNode* m_next; ///< Next in list + Node* m_node; ///< Node + }; + + /// Variables for finding a split partition + struct PartitionVars + { + int m_partition[MAXNODES + 1]; + int m_total; + int m_minFill; + int m_taken[MAXNODES + 1]; + int m_count[2]; + Rect m_cover[2]; + ELEMTYPEREAL m_area[2]; + + Branch m_branchBuf[MAXNODES + 1]; + int m_branchCount; + Rect m_coverSplit; + ELEMTYPEREAL m_coverSplitArea; + }; + + /// Data structure used for Nearest Neighbor search implementation + struct NNNode + { + Branch m_branch; + ELEMTYPE minDist; + bool isLeaf; + }; + + Node* AllocNode(); + void FreeNode( Node* a_node ); + void InitNode( Node* a_node ); + void InitRect( Rect* a_rect ); + bool InsertRectRec( Rect* a_rect, + const DATATYPE& a_id, + Node* a_node, + Node** a_newNode, + int a_level ); + bool InsertRect( Rect* a_rect, const DATATYPE& a_id, Node** a_root, int a_level ); + Rect NodeCover( Node* a_node ); + bool AddBranch( Branch* a_branch, Node* a_node, Node** a_newNode ); + void DisconnectBranch( Node* a_node, int a_index ); + int PickBranch( Rect* a_rect, Node* a_node ); + Rect CombineRect( Rect* a_rectA, Rect* a_rectB ); + void SplitNode( Node* a_node, Branch* a_branch, Node** a_newNode ); + ELEMTYPEREAL RectSphericalVolume( Rect* a_rect ); + ELEMTYPEREAL RectVolume( Rect* a_rect ); + ELEMTYPEREAL CalcRectVolume( Rect* a_rect ); + void GetBranches( Node* a_node, Branch* a_branch, PartitionVars* a_parVars ); + void ChoosePartition( PartitionVars* a_parVars, int a_minFill ); + void LoadNodes( Node* a_nodeA, Node* a_nodeB, PartitionVars* a_parVars ); + void InitParVars( PartitionVars* a_parVars, int a_maxRects, int a_minFill ); + void PickSeeds( PartitionVars* a_parVars ); + void Classify( int a_index, int a_group, PartitionVars* a_parVars ); + bool RemoveRect( Rect* a_rect, const DATATYPE& a_id, Node** a_root ); + bool RemoveRectRec( Rect* a_rect, + const DATATYPE& a_id, + Node* a_node, + ListNode** a_listNode ); + ListNode* AllocListNode(); + void FreeListNode( ListNode* a_listNode ); + bool Overlap( Rect* a_rectA, Rect* a_rectB ); + void ReInsert( Node* a_node, ListNode** a_listNode ); + ELEMTYPE MinDist( const ELEMTYPE a_point[NUMDIMS], Rect* a_rect ); + void InsertNNListSorted( std::vector* nodeList, NNNode* newNode ); + + bool Search( Node * a_node, Rect * a_rect, int& a_foundCount, bool a_resultCallback( + DATATYPE a_data, + void* a_context ), void* a_context ); + + template + bool Search( Node* a_node, Rect* a_rect, VISITOR& a_visitor, int& a_foundCount ) + { + ASSERT( a_node ); + ASSERT( a_node->m_level >= 0 ); + ASSERT( a_rect ); + + if( a_node->IsInternalNode() ) // This is an internal node in the tree + { + for( int index = 0; index < a_node->m_count; ++index ) + { + if( Overlap( a_rect, &a_node->m_branch[index].m_rect ) ) + { + if( !Search( a_node->m_branch[index].m_child, a_rect, a_visitor, a_foundCount ) ) + { + return false; // Don't continue searching + } + } + } + } + else // This is a leaf node + { + for( int index = 0; index < a_node->m_count; ++index ) + { + if( Overlap( a_rect, &a_node->m_branch[index].m_rect ) ) + { + DATATYPE& id = a_node->m_branch[index].m_data; + + if( !a_visitor( id ) ) + return false; + + a_foundCount++; + } + } + } + + return true; // Continue searching + } + + void RemoveAllRec( Node* a_node ); + void Reset(); + void CountRec( Node* a_node, int& a_count ); + + bool SaveRec( Node* a_node, RTFileStream& a_stream ); + bool LoadRec( Node* a_node, RTFileStream& a_stream ); + + Node* m_root; ///< Root of tree + ELEMTYPEREAL m_unitSphereVolume; ///< Unit sphere constant for required number of dimensions +}; + + +// Because there is not stream support, this is a quick and dirty file I/O helper. +// Users will likely replace its usage with a Stream implementation from their favorite API. +class RTFileStream +{ + FILE* m_file; +public: + + + RTFileStream() + { + m_file = NULL; + } + + ~RTFileStream() + { + Close(); + } + + bool OpenRead( const char* a_fileName ) + { + m_file = fopen( a_fileName, "rb" ); + + if( !m_file ) + { + return false; + } + + return true; + } + + bool OpenWrite( const char* a_fileName ) + { + m_file = fopen( a_fileName, "wb" ); + + if( !m_file ) + { + return false; + } + + return true; + } + + void Close() + { + if( m_file ) + { + fclose( m_file ); + m_file = NULL; + } + } + + template + size_t Write( const TYPE& a_value ) + { + ASSERT( m_file ); + return fwrite( (void*) &a_value, sizeof(a_value), 1, m_file ); + } + + template + size_t WriteArray( const TYPE* a_array, int a_count ) + { + ASSERT( m_file ); + return fwrite( (void*) a_array, sizeof(TYPE) * a_count, 1, m_file ); + } + + template + size_t Read( TYPE& a_value ) + { + ASSERT( m_file ); + return fread( (void*) &a_value, sizeof(a_value), 1, m_file ); + } + + template + size_t ReadArray( TYPE* a_array, int a_count ) + { + ASSERT( m_file ); + return fread( (void*) a_array, sizeof(TYPE) * a_count, 1, m_file ); + } +}; + + +RTREE_TEMPLATE RTREE_QUAL::RTree() +{ + ASSERT( MAXNODES > MINNODES ); + ASSERT( MINNODES > 0 ); + + + // We only support machine word size simple data type eg. integer index or object pointer. + // Since we are storing as union with non data branch + ASSERT( sizeof(DATATYPE) == sizeof(void*) || sizeof(DATATYPE) == sizeof(int) ); + + // Precomputed volumes of the unit spheres for the first few dimensions + const float UNIT_SPHERE_VOLUMES[] = + { + 0.000000f, 2.000000f, 3.141593f, // Dimension 0,1,2 + 4.188790f, 4.934802f, 5.263789f, // Dimension 3,4,5 + 5.167713f, 4.724766f, 4.058712f, // Dimension 6,7,8 + 3.298509f, 2.550164f, 1.884104f, // Dimension 9,10,11 + 1.335263f, 0.910629f, 0.599265f, // Dimension 12,13,14 + 0.381443f, 0.235331f, 0.140981f, // Dimension 15,16,17 + 0.082146f, 0.046622f, 0.025807f, // Dimension 18,19,20 + }; + + m_root = AllocNode(); + m_root->m_level = 0; + m_unitSphereVolume = (ELEMTYPEREAL) UNIT_SPHERE_VOLUMES[NUMDIMS]; +} + + +RTREE_TEMPLATE +RTREE_QUAL::~RTree() { + Reset(); // Free, or reset node memory +} + + +RTREE_TEMPLATE +void RTREE_QUAL::Insert( const ELEMTYPE a_min[NUMDIMS], + const ELEMTYPE a_max[NUMDIMS], + const DATATYPE& a_dataId ) +{ +#ifdef _DEBUG + + for( int index = 0; indexNearestNeighbor( a_point, 0, 0 ); +} + + +RTREE_TEMPLATE +DATATYPE RTREE_QUAL::NearestNeighbor( const ELEMTYPE a_point[NUMDIMS], + ELEMTYPE a_squareDistanceCallback( const ELEMTYPE a_point[NUMDIMS], DATATYPE a_data ), + ELEMTYPE* a_squareDistance ) +{ + typedef typename std::vector::iterator iterator; + std::vector nodeList; + Node* node = m_root; + NNNode* closestNode = 0; + while( !closestNode || !closestNode->isLeaf ) + { + //check every node on this level + for( int index = 0; index < node->m_count; ++index ) + { + NNNode* newNode = new NNNode; + newNode->isLeaf = node->IsLeaf(); + newNode->m_branch = node->m_branch[index]; + if( newNode->isLeaf && a_squareDistanceCallback ) + newNode->minDist = a_squareDistanceCallback( a_point, newNode->m_branch.m_data ); + else + newNode->minDist = this->MinDist( a_point, &(node->m_branch[index].m_rect) ); + + //TODO: a custom list could be more efficient than a vector + this->InsertNNListSorted( &nodeList, newNode ); + } + if( nodeList.size() == 0 ) + { + return 0; + } + closestNode = nodeList.back(); + node = closestNode->m_branch.m_child; + nodeList.pop_back(); + free(closestNode); + } + + // free memory used for remaining NNNodes in nodeList + for( iterator iter = nodeList.begin(); iter != nodeList.end(); ++iter ) + { + NNNode* node = *iter; + free(node); + } + + *a_squareDistance = closestNode->minDist; + return closestNode->m_branch.m_data; +} + + +RTREE_TEMPLATE +int RTREE_QUAL::Count() +{ + int count = 0; + + CountRec( m_root, count ); + + return count; +} + + +RTREE_TEMPLATE +void RTREE_QUAL::CountRec( Node* a_node, int& a_count ) +{ + if( a_node->IsInternalNode() ) // not a leaf node + { + for( int index = 0; index < a_node->m_count; ++index ) + { + CountRec( a_node->m_branch[index].m_child, a_count ); + } + } + else // A leaf node + { + a_count += a_node->m_count; + } +} + + +RTREE_TEMPLATE +bool RTREE_QUAL::Load( const char* a_fileName ) +{ + RemoveAll(); // Clear existing tree + + RTFileStream stream; + + if( !stream.OpenRead( a_fileName ) ) + { + return false; + } + + bool result = Load( stream ); + + stream.Close(); + + return result; +}; + + +RTREE_TEMPLATE +bool RTREE_QUAL::Load( RTFileStream& a_stream ) +{ + // Write some kind of header + int _dataFileId = ('R' << 0) | ('T' << 8) | ('R' << 16) | ('E' << 24); + int _dataSize = sizeof(DATATYPE); + int _dataNumDims = NUMDIMS; + int _dataElemSize = sizeof(ELEMTYPE); + int _dataElemRealSize = sizeof(ELEMTYPEREAL); + int _dataMaxNodes = TMAXNODES; + int _dataMinNodes = TMINNODES; + + int dataFileId = 0; + int dataSize = 0; + int dataNumDims = 0; + int dataElemSize = 0; + int dataElemRealSize = 0; + int dataMaxNodes = 0; + int dataMinNodes = 0; + + a_stream.Read( dataFileId ); + a_stream.Read( dataSize ); + a_stream.Read( dataNumDims ); + a_stream.Read( dataElemSize ); + a_stream.Read( dataElemRealSize ); + a_stream.Read( dataMaxNodes ); + a_stream.Read( dataMinNodes ); + + bool result = false; + + // Test if header was valid and compatible + if( (dataFileId == _dataFileId) + && (dataSize == _dataSize) + && (dataNumDims == _dataNumDims) + && (dataElemSize == _dataElemSize) + && (dataElemRealSize == _dataElemRealSize) + && (dataMaxNodes == _dataMaxNodes) + && (dataMinNodes == _dataMinNodes) + ) + { + // Recursively load tree + result = LoadRec( m_root, a_stream ); + } + + return result; +} + + +RTREE_TEMPLATE +bool RTREE_QUAL::LoadRec( Node* a_node, RTFileStream& a_stream ) +{ + a_stream.Read( a_node->m_level ); + a_stream.Read( a_node->m_count ); + + if( a_node->IsInternalNode() ) // not a leaf node + { + for( int index = 0; index < a_node->m_count; ++index ) + { + Branch* curBranch = &a_node->m_branch[index]; + + a_stream.ReadArray( curBranch->m_rect.m_min, NUMDIMS ); + a_stream.ReadArray( curBranch->m_rect.m_max, NUMDIMS ); + + curBranch->m_child = AllocNode(); + LoadRec( curBranch->m_child, a_stream ); + } + } + else // A leaf node + { + for( int index = 0; index < a_node->m_count; ++index ) + { + Branch* curBranch = &a_node->m_branch[index]; + + a_stream.ReadArray( curBranch->m_rect.m_min, NUMDIMS ); + a_stream.ReadArray( curBranch->m_rect.m_max, NUMDIMS ); + + a_stream.Read( curBranch->m_data ); + } + } + + return true; // Should do more error checking on I/O operations +} + + +RTREE_TEMPLATE +bool RTREE_QUAL::Save( const char* a_fileName ) +{ + RTFileStream stream; + + if( !stream.OpenWrite( a_fileName ) ) + { + return false; + } + + bool result = Save( stream ); + + stream.Close(); + + return result; +} + + +RTREE_TEMPLATE +bool RTREE_QUAL::Save( RTFileStream& a_stream ) +{ + // Write some kind of header + int dataFileId = ('R' << 0) | ('T' << 8) | ('R' << 16) | ('E' << 24); + int dataSize = sizeof(DATATYPE); + int dataNumDims = NUMDIMS; + int dataElemSize = sizeof(ELEMTYPE); + int dataElemRealSize = sizeof(ELEMTYPEREAL); + int dataMaxNodes = TMAXNODES; + int dataMinNodes = TMINNODES; + + a_stream.Write( dataFileId ); + a_stream.Write( dataSize ); + a_stream.Write( dataNumDims ); + a_stream.Write( dataElemSize ); + a_stream.Write( dataElemRealSize ); + a_stream.Write( dataMaxNodes ); + a_stream.Write( dataMinNodes ); + + // Recursively save tree + bool result = SaveRec( m_root, a_stream ); + + return result; +} + + +RTREE_TEMPLATE +bool RTREE_QUAL::SaveRec( Node* a_node, RTFileStream& a_stream ) +{ + a_stream.Write( a_node->m_level ); + a_stream.Write( a_node->m_count ); + + if( a_node->IsInternalNode() ) // not a leaf node + { + for( int index = 0; index < a_node->m_count; ++index ) + { + Branch* curBranch = &a_node->m_branch[index]; + + a_stream.WriteArray( curBranch->m_rect.m_min, NUMDIMS ); + a_stream.WriteArray( curBranch->m_rect.m_max, NUMDIMS ); + + SaveRec( curBranch->m_child, a_stream ); + } + } + else // A leaf node + { + for( int index = 0; index < a_node->m_count; ++index ) + { + Branch* curBranch = &a_node->m_branch[index]; + + a_stream.WriteArray( curBranch->m_rect.m_min, NUMDIMS ); + a_stream.WriteArray( curBranch->m_rect.m_max, NUMDIMS ); + + a_stream.Write( curBranch->m_data ); + } + } + + return true; // Should do more error checking on I/O operations +} + + +RTREE_TEMPLATE +void RTREE_QUAL::RemoveAll() +{ + // Delete all existing nodes + Reset(); + + m_root = AllocNode(); + m_root->m_level = 0; +} + + +RTREE_TEMPLATE +void RTREE_QUAL::Reset() +{ +#ifdef RTREE_DONT_USE_MEMPOOLS + // Delete all existing nodes + RemoveAllRec( m_root ); +#else // RTREE_DONT_USE_MEMPOOLS + // Just reset memory pools. We are not using complex types + // EXAMPLE +#endif // RTREE_DONT_USE_MEMPOOLS +} + + +RTREE_TEMPLATE +void RTREE_QUAL::RemoveAllRec( Node* a_node ) +{ + ASSERT( a_node ); + ASSERT( a_node->m_level >= 0 ); + + if( a_node->IsInternalNode() ) // This is an internal node in the tree + { + for( int index = 0; index < a_node->m_count; ++index ) + { + RemoveAllRec( a_node->m_branch[index].m_child ); + } + } + + FreeNode( a_node ); +} + + +RTREE_TEMPLATE +typename RTREE_QUAL::Node* RTREE_QUAL::AllocNode() +{ + Node* newNode; + +#ifdef RTREE_DONT_USE_MEMPOOLS + newNode = new Node; +#else // RTREE_DONT_USE_MEMPOOLS + // EXAMPLE +#endif // RTREE_DONT_USE_MEMPOOLS + InitNode( newNode ); + return newNode; +} + + +RTREE_TEMPLATE +void RTREE_QUAL::FreeNode( Node* a_node ) +{ + ASSERT( a_node ); + +#ifdef RTREE_DONT_USE_MEMPOOLS + delete a_node; +#else // RTREE_DONT_USE_MEMPOOLS + // EXAMPLE +#endif // RTREE_DONT_USE_MEMPOOLS +} + + +// Allocate space for a node in the list used in DeletRect to +// store Nodes that are too empty. +RTREE_TEMPLATE +typename RTREE_QUAL::ListNode* RTREE_QUAL::AllocListNode() +{ +#ifdef RTREE_DONT_USE_MEMPOOLS + return new ListNode; +#else // RTREE_DONT_USE_MEMPOOLS + // EXAMPLE +#endif // RTREE_DONT_USE_MEMPOOLS +} + + +RTREE_TEMPLATE +void RTREE_QUAL::FreeListNode( ListNode* a_listNode ) +{ +#ifdef RTREE_DONT_USE_MEMPOOLS + delete a_listNode; +#else // RTREE_DONT_USE_MEMPOOLS + // EXAMPLE +#endif // RTREE_DONT_USE_MEMPOOLS +} + + +RTREE_TEMPLATE +void RTREE_QUAL::InitNode( Node* a_node ) +{ + a_node->m_count = 0; + a_node->m_level = -1; +} + + +RTREE_TEMPLATE +void RTREE_QUAL::InitRect( Rect* a_rect ) +{ + for( int index = 0; index < NUMDIMS; ++index ) + { + a_rect->m_min[index] = (ELEMTYPE) 0; + a_rect->m_max[index] = (ELEMTYPE) 0; + } +} + + +// Inserts a new data rectangle into the index structure. +// Recursively descends tree, propagates splits back up. +// Returns 0 if node was not split. Old node updated. +// If node was split, returns 1 and sets the pointer pointed to by +// new_node to point to the new node. Old node updated to become one of two. +// The level argument specifies the number of steps up from the leaf +// level to insert; e.g. a data rectangle goes in at level = 0. +RTREE_TEMPLATE +bool RTREE_QUAL::InsertRectRec( Rect* a_rect, + const DATATYPE& a_id, + Node* a_node, + Node** a_newNode, + int a_level ) +{ + ASSERT( a_rect && a_node && a_newNode ); + ASSERT( a_level >= 0 && a_level <= a_node->m_level ); + + int index; + Branch branch; + Node* otherNode; + + // Still above level for insertion, go down tree recursively + if( a_node->m_level > a_level ) + { + index = PickBranch( a_rect, a_node ); + + if( !InsertRectRec( a_rect, a_id, a_node->m_branch[index].m_child, &otherNode, a_level ) ) + { + // Child was not split + a_node->m_branch[index].m_rect = + CombineRect( a_rect, &(a_node->m_branch[index].m_rect) ); + return false; + } + else // Child was split + { + a_node->m_branch[index].m_rect = NodeCover( a_node->m_branch[index].m_child ); + branch.m_child = otherNode; + branch.m_rect = NodeCover( otherNode ); + return AddBranch( &branch, a_node, a_newNode ); + } + } + else if( a_node->m_level == a_level ) // Have reached level for insertion. Add rect, split if necessary + { + branch.m_rect = *a_rect; + branch.m_child = (Node*) a_id; + // Child field of leaves contains id of data record + return AddBranch( &branch, a_node, a_newNode ); + } + else + { + // Should never occur + ASSERT( 0 ); + return false; + } +} + + +// Insert a data rectangle into an index structure. +// InsertRect provides for splitting the root; +// returns 1 if root was split, 0 if it was not. +// The level argument specifies the number of steps up from the leaf +// level to insert; e.g. a data rectangle goes in at level = 0. +// InsertRect2 does the recursion. +// +RTREE_TEMPLATE +bool RTREE_QUAL::InsertRect( Rect* a_rect, const DATATYPE& a_id, Node** a_root, int a_level ) +{ + ASSERT( a_rect && a_root ); + ASSERT( a_level >= 0 && a_level <= (*a_root)->m_level ); +#ifdef _DEBUG + + for( int index = 0; index < NUMDIMS; ++index ) + { + ASSERT( a_rect->m_min[index] <= a_rect->m_max[index] ); + } + +#endif // _DEBUG + + Node* newRoot; + Node* newNode; + Branch branch; + + if( InsertRectRec( a_rect, a_id, *a_root, &newNode, a_level ) ) // Root split + { + newRoot = AllocNode(); // Grow tree taller and new root + newRoot->m_level = (*a_root)->m_level + 1; + branch.m_rect = NodeCover( *a_root ); + branch.m_child = *a_root; + AddBranch( &branch, newRoot, NULL ); + branch.m_rect = NodeCover( newNode ); + branch.m_child = newNode; + AddBranch( &branch, newRoot, NULL ); + *a_root = newRoot; + return true; + } + + return false; +} + + +// Find the smallest rectangle that includes all rectangles in branches of a node. +RTREE_TEMPLATE +typename RTREE_QUAL::Rect RTREE_QUAL::NodeCover( Node* a_node ) +{ + ASSERT( a_node ); + + int firstTime = true; + Rect rect; + InitRect( &rect ); + + for( int index = 0; index < a_node->m_count; ++index ) + { + if( firstTime ) + { + rect = a_node->m_branch[index].m_rect; + firstTime = false; + } + else + { + rect = CombineRect( &rect, &(a_node->m_branch[index].m_rect) ); + } + } + + return rect; +} + + +// Add a branch to a node. Split the node if necessary. +// Returns 0 if node not split. Old node updated. +// Returns 1 if node split, sets *new_node to address of new node. +// Old node updated, becomes one of two. +RTREE_TEMPLATE +bool RTREE_QUAL::AddBranch( Branch* a_branch, Node* a_node, Node** a_newNode ) +{ + ASSERT( a_branch ); + ASSERT( a_node ); + + if( a_node->m_count < MAXNODES ) // Split won't be necessary + { + a_node->m_branch[a_node->m_count] = *a_branch; + ++a_node->m_count; + + return false; + } + else + { + ASSERT( a_newNode ); + + SplitNode( a_node, a_branch, a_newNode ); + return true; + } +} + + +// Disconnect a dependent node. +// Caller must return (or stop using iteration index) after this as count has changed +RTREE_TEMPLATE +void RTREE_QUAL::DisconnectBranch( Node* a_node, int a_index ) +{ + ASSERT( a_node && (a_index >= 0) && (a_index < MAXNODES) ); + ASSERT( a_node->m_count > 0 ); + + // Remove element by swapping with the last element to prevent gaps in array + a_node->m_branch[a_index] = a_node->m_branch[a_node->m_count - 1]; + + --a_node->m_count; +} + + +// Pick a branch. Pick the one that will need the smallest increase +// in area to accomodate the new rectangle. This will result in the +// least total area for the covering rectangles in the current node. +// In case of a tie, pick the one which was smaller before, to get +// the best resolution when searching. +RTREE_TEMPLATE +int RTREE_QUAL::PickBranch( Rect* a_rect, Node* a_node ) +{ + ASSERT( a_rect && a_node ); + + bool firstTime = true; + ELEMTYPEREAL increase; + ELEMTYPEREAL bestIncr = (ELEMTYPEREAL) -1; + ELEMTYPEREAL area; + ELEMTYPEREAL bestArea; + int best; + Rect tempRect; + + for( int index = 0; index < a_node->m_count; ++index ) + { + Rect* curRect = &a_node->m_branch[index].m_rect; + area = CalcRectVolume( curRect ); + tempRect = CombineRect( a_rect, curRect ); + increase = CalcRectVolume( &tempRect ) - area; + + if( (increase < bestIncr) || firstTime ) + { + best = index; + bestArea = area; + bestIncr = increase; + firstTime = false; + } + else if( (increase == bestIncr) && (area < bestArea) ) + { + best = index; + bestArea = area; + bestIncr = increase; + } + } + + return best; +} + + +// Combine two rectangles into larger one containing both +RTREE_TEMPLATE +typename RTREE_QUAL::Rect RTREE_QUAL::CombineRect( Rect* a_rectA, Rect* a_rectB ) +{ + ASSERT( a_rectA && a_rectB ); + + Rect newRect; + + for( int index = 0; index < NUMDIMS; ++index ) + { + newRect.m_min[index] = Min( a_rectA->m_min[index], a_rectB->m_min[index] ); + newRect.m_max[index] = Max( a_rectA->m_max[index], a_rectB->m_max[index] ); + } + + return newRect; +} + + +// Split a node. +// Divides the nodes branches and the extra one between two nodes. +// Old node is one of the new ones, and one really new one is created. +// Tries more than one method for choosing a partition, uses best result. +RTREE_TEMPLATE +void RTREE_QUAL::SplitNode( Node* a_node, Branch* a_branch, Node** a_newNode ) +{ + ASSERT( a_node ); + ASSERT( a_branch ); + + // Could just use local here, but member or external is faster since it is reused + PartitionVars localVars; + PartitionVars* parVars = &localVars; + int level; + + // Load all the branches into a buffer, initialize old node + level = a_node->m_level; + GetBranches( a_node, a_branch, parVars ); + + // Find partition + ChoosePartition( parVars, MINNODES ); + + // Put branches from buffer into 2 nodes according to chosen partition + *a_newNode = AllocNode(); + (*a_newNode)->m_level = a_node->m_level = level; + LoadNodes( a_node, *a_newNode, parVars ); + + ASSERT( (a_node->m_count + (*a_newNode)->m_count) == parVars->m_total ); +} + + +// Calculate the n-dimensional volume of a rectangle +RTREE_TEMPLATE +ELEMTYPEREAL RTREE_QUAL::RectVolume( Rect* a_rect ) +{ + ASSERT( a_rect ); + + ELEMTYPEREAL volume = (ELEMTYPEREAL) 1; + + for( int index = 0; indexm_max[index] - a_rect->m_min[index]; + } + + ASSERT( volume >= (ELEMTYPEREAL) 0 ); + + return volume; +} + + +// The exact volume of the bounding sphere for the given Rect +RTREE_TEMPLATE +ELEMTYPEREAL RTREE_QUAL::RectSphericalVolume( Rect* a_rect ) +{ + ASSERT( a_rect ); + + ELEMTYPEREAL sumOfSquares = (ELEMTYPEREAL) 0; + ELEMTYPEREAL radius; + + for( int index = 0; index < NUMDIMS; ++index ) + { + ELEMTYPEREAL halfExtent = + ( (ELEMTYPEREAL) a_rect->m_max[index] - (ELEMTYPEREAL) a_rect->m_min[index] ) * 0.5f; + sumOfSquares += halfExtent * halfExtent; + } + + radius = (ELEMTYPEREAL) sqrt( sumOfSquares ); + + // Pow maybe slow, so test for common dims like 2,3 and just use x*x, x*x*x. + if( NUMDIMS == 3 ) + { + return radius * radius * radius * m_unitSphereVolume; + } + else if( NUMDIMS == 2 ) + { + return radius * radius * m_unitSphereVolume; + } + else + { + return (ELEMTYPEREAL) (pow( radius, NUMDIMS ) * m_unitSphereVolume); + } +} + + +// Use one of the methods to calculate retangle volume +RTREE_TEMPLATE +ELEMTYPEREAL RTREE_QUAL::CalcRectVolume( Rect* a_rect ) +{ +#ifdef RTREE_USE_SPHERICAL_VOLUME + return RectSphericalVolume( a_rect ); // Slower but helps certain merge cases +#else // RTREE_USE_SPHERICAL_VOLUME + return RectVolume( a_rect ); // Faster but can cause poor merges +#endif // RTREE_USE_SPHERICAL_VOLUME +} + + +// Load branch buffer with branches from full node plus the extra branch. +RTREE_TEMPLATE +void RTREE_QUAL::GetBranches( Node* a_node, Branch* a_branch, PartitionVars* a_parVars ) +{ + ASSERT( a_node ); + ASSERT( a_branch ); + + ASSERT( a_node->m_count == MAXNODES ); + + // Load the branch buffer + for( int index = 0; index < MAXNODES; ++index ) + { + a_parVars->m_branchBuf[index] = a_node->m_branch[index]; + } + + a_parVars->m_branchBuf[MAXNODES] = *a_branch; + a_parVars->m_branchCount = MAXNODES + 1; + + // Calculate rect containing all in the set + a_parVars->m_coverSplit = a_parVars->m_branchBuf[0].m_rect; + + for( int index = 1; index < MAXNODES + 1; ++index ) + { + a_parVars->m_coverSplit = + CombineRect( &a_parVars->m_coverSplit, &a_parVars->m_branchBuf[index].m_rect ); + } + + a_parVars->m_coverSplitArea = CalcRectVolume( &a_parVars->m_coverSplit ); + + InitNode( a_node ); +} + + +// Method #0 for choosing a partition: +// As the seeds for the two groups, pick the two rects that would waste the +// most area if covered by a single rectangle, i.e. evidently the worst pair +// to have in the same group. +// Of the remaining, one at a time is chosen to be put in one of the two groups. +// The one chosen is the one with the greatest difference in area expansion +// depending on which group - the rect most strongly attracted to one group +// and repelled from the other. +// If one group gets too full (more would force other group to violate min +// fill requirement) then other group gets the rest. +// These last are the ones that can go in either group most easily. +RTREE_TEMPLATE +void RTREE_QUAL::ChoosePartition( PartitionVars* a_parVars, int a_minFill ) +{ + ASSERT( a_parVars ); + + ELEMTYPEREAL biggestDiff; + int group, chosen, betterGroup; + + InitParVars( a_parVars, a_parVars->m_branchCount, a_minFill ); + PickSeeds( a_parVars ); + + while( ( (a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total ) + && ( a_parVars->m_count[0] < (a_parVars->m_total - a_parVars->m_minFill) ) + && ( a_parVars->m_count[1] < (a_parVars->m_total - a_parVars->m_minFill) ) ) + { + biggestDiff = (ELEMTYPEREAL) -1; + + for( int index = 0; indexm_total; ++index ) + { + if( !a_parVars->m_taken[index] ) + { + Rect* curRect = &a_parVars->m_branchBuf[index].m_rect; + Rect rect0 = CombineRect( curRect, &a_parVars->m_cover[0] ); + Rect rect1 = CombineRect( curRect, &a_parVars->m_cover[1] ); + ELEMTYPEREAL growth0 = CalcRectVolume( &rect0 ) - a_parVars->m_area[0]; + ELEMTYPEREAL growth1 = CalcRectVolume( &rect1 ) - a_parVars->m_area[1]; + ELEMTYPEREAL diff = growth1 - growth0; + + if( diff >= 0 ) + { + group = 0; + } + else + { + group = 1; + diff = -diff; + } + + if( diff > biggestDiff ) + { + biggestDiff = diff; + chosen = index; + betterGroup = group; + } + else if( (diff == biggestDiff) + && (a_parVars->m_count[group] < a_parVars->m_count[betterGroup]) ) + { + chosen = index; + betterGroup = group; + } + } + } + + Classify( chosen, betterGroup, a_parVars ); + } + + // If one group too full, put remaining rects in the other + if( (a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total ) + { + if( a_parVars->m_count[0] >= a_parVars->m_total - a_parVars->m_minFill ) + { + group = 1; + } + else + { + group = 0; + } + + for( int index = 0; indexm_total; ++index ) + { + if( !a_parVars->m_taken[index] ) + { + Classify( index, group, a_parVars ); + } + } + } + + ASSERT( (a_parVars->m_count[0] + a_parVars->m_count[1]) == a_parVars->m_total ); + ASSERT( (a_parVars->m_count[0] >= a_parVars->m_minFill) + && (a_parVars->m_count[1] >= a_parVars->m_minFill) ); +} + + +// Copy branches from the buffer into two nodes according to the partition. +RTREE_TEMPLATE +void RTREE_QUAL::LoadNodes( Node* a_nodeA, Node* a_nodeB, PartitionVars* a_parVars ) +{ + ASSERT( a_nodeA ); + ASSERT( a_nodeB ); + ASSERT( a_parVars ); + + for( int index = 0; index < a_parVars->m_total; ++index ) + { + ASSERT( a_parVars->m_partition[index] == 0 || a_parVars->m_partition[index] == 1 ); + + if( a_parVars->m_partition[index] == 0 ) + { + AddBranch( &a_parVars->m_branchBuf[index], a_nodeA, NULL ); + } + else if( a_parVars->m_partition[index] == 1 ) + { + AddBranch( &a_parVars->m_branchBuf[index], a_nodeB, NULL ); + } + } +} + + +// Initialize a PartitionVars structure. +RTREE_TEMPLATE +void RTREE_QUAL::InitParVars( PartitionVars* a_parVars, int a_maxRects, int a_minFill ) +{ + ASSERT( a_parVars ); + + a_parVars->m_count[0] = a_parVars->m_count[1] = 0; + a_parVars->m_area[0] = a_parVars->m_area[1] = (ELEMTYPEREAL) 0; + a_parVars->m_total = a_maxRects; + a_parVars->m_minFill = a_minFill; + + for( int index = 0; index < a_maxRects; ++index ) + { + a_parVars->m_taken[index] = false; + a_parVars->m_partition[index] = -1; + } +} + + +RTREE_TEMPLATE +void RTREE_QUAL::PickSeeds( PartitionVars* a_parVars ) +{ + int seed0, seed1; + ELEMTYPEREAL worst, waste; + ELEMTYPEREAL area[MAXNODES + 1]; + + for( int index = 0; indexm_total; ++index ) + { + area[index] = CalcRectVolume( &a_parVars->m_branchBuf[index].m_rect ); + } + + worst = -a_parVars->m_coverSplitArea - 1; + + for( int indexA = 0; indexA < a_parVars->m_total - 1; ++indexA ) + { + for( int indexB = indexA + 1; indexB < a_parVars->m_total; ++indexB ) + { + Rect oneRect = CombineRect( &a_parVars->m_branchBuf[indexA].m_rect, + &a_parVars->m_branchBuf[indexB].m_rect ); + waste = CalcRectVolume( &oneRect ) - area[indexA] - area[indexB]; + + if( waste > worst ) + { + worst = waste; + seed0 = indexA; + seed1 = indexB; + } + } + } + + Classify( seed0, 0, a_parVars ); + Classify( seed1, 1, a_parVars ); +} + + +// Put a branch in one of the groups. +RTREE_TEMPLATE +void RTREE_QUAL::Classify( int a_index, int a_group, PartitionVars* a_parVars ) +{ + ASSERT( a_parVars ); + ASSERT( !a_parVars->m_taken[a_index] ); + + a_parVars->m_partition[a_index] = a_group; + a_parVars->m_taken[a_index] = true; + + if( a_parVars->m_count[a_group] == 0 ) + { + a_parVars->m_cover[a_group] = a_parVars->m_branchBuf[a_index].m_rect; + } + else + { + a_parVars->m_cover[a_group] = CombineRect( &a_parVars->m_branchBuf[a_index].m_rect, + &a_parVars->m_cover[a_group] ); + } + + a_parVars->m_area[a_group] = CalcRectVolume( &a_parVars->m_cover[a_group] ); + ++a_parVars->m_count[a_group]; +} + + +// Delete a data rectangle from an index structure. +// Pass in a pointer to a Rect, the tid of the record, ptr to ptr to root node. +// Returns 1 if record not found, 0 if success. +// RemoveRect provides for eliminating the root. +RTREE_TEMPLATE +bool RTREE_QUAL::RemoveRect( Rect* a_rect, const DATATYPE& a_id, Node** a_root ) +{ + ASSERT( a_rect && a_root ); + ASSERT( *a_root ); + + Node* tempNode; + ListNode* reInsertList = NULL; + + if( !RemoveRectRec( a_rect, a_id, *a_root, &reInsertList ) ) + { + // Found and deleted a data item + // Reinsert any branches from eliminated nodes + while( reInsertList ) + { + tempNode = reInsertList->m_node; + + for( int index = 0; index < tempNode->m_count; ++index ) + { + InsertRect( &(tempNode->m_branch[index].m_rect), + tempNode->m_branch[index].m_data, + a_root, + tempNode->m_level ); + } + + ListNode* remLNode = reInsertList; + reInsertList = reInsertList->m_next; + + FreeNode( remLNode->m_node ); + FreeListNode( remLNode ); + } + + // Check for redundant root (not leaf, 1 child) and eliminate + if( (*a_root)->m_count == 1 && (*a_root)->IsInternalNode() ) + { + tempNode = (*a_root)->m_branch[0].m_child; + + ASSERT( tempNode ); + FreeNode( *a_root ); + *a_root = tempNode; + } + + return false; + } + else + { + return true; + } +} + + +// Delete a rectangle from non-root part of an index structure. +// Called by RemoveRect. Descends tree recursively, +// merges branches on the way back up. +// Returns 1 if record not found, 0 if success. +RTREE_TEMPLATE +bool RTREE_QUAL::RemoveRectRec( Rect* a_rect, + const DATATYPE& a_id, + Node* a_node, + ListNode** a_listNode ) +{ + ASSERT( a_rect && a_node && a_listNode ); + ASSERT( a_node->m_level >= 0 ); + + if( a_node->IsInternalNode() ) // not a leaf node + { + for( int index = 0; index < a_node->m_count; ++index ) + { + if( Overlap( a_rect, &(a_node->m_branch[index].m_rect) ) ) + { + if( !RemoveRectRec( a_rect, a_id, a_node->m_branch[index].m_child, a_listNode ) ) + { + if( a_node->m_branch[index].m_child->m_count >= MINNODES ) + { + // child removed, just resize parent rect + a_node->m_branch[index].m_rect = + NodeCover( a_node->m_branch[index].m_child ); + } + else + { + // child removed, not enough entries in node, eliminate node + ReInsert( a_node->m_branch[index].m_child, a_listNode ); + DisconnectBranch( a_node, index ); // Must return after this call as count has changed + } + + return false; + } + } + } + + return true; + } + else // A leaf node + { + for( int index = 0; index < a_node->m_count; ++index ) + { + if( a_node->m_branch[index].m_child == (Node*) a_id ) + { + DisconnectBranch( a_node, index ); // Must return after this call as count has changed + return false; + } + } + + return true; + } +} + + +// Decide whether two rectangles overlap. +RTREE_TEMPLATE +bool RTREE_QUAL::Overlap( Rect* a_rectA, Rect* a_rectB ) +{ + ASSERT( a_rectA && a_rectB ); + + for( int index = 0; index < NUMDIMS; ++index ) + { + if( a_rectA->m_min[index] > a_rectB->m_max[index] + || a_rectB->m_min[index] > a_rectA->m_max[index] ) + { + return false; + } + } + + return true; +} + + +// Add a node to the reinsertion list. All its branches will later +// be reinserted into the index structure. +RTREE_TEMPLATE +void RTREE_QUAL::ReInsert( Node* a_node, ListNode** a_listNode ) +{ + ListNode* newListNode; + + newListNode = AllocListNode(); + newListNode->m_node = a_node; + newListNode->m_next = *a_listNode; + *a_listNode = newListNode; +} + + +// Search in an index tree or subtree for all data retangles that overlap the argument rectangle. +RTREE_TEMPLATE +bool RTREE_QUAL::Search( Node* a_node, Rect* a_rect, int& a_foundCount, bool a_resultCallback( + DATATYPE a_data, + void* a_context ), void* a_context ) +{ + ASSERT( a_node ); + ASSERT( a_node->m_level >= 0 ); + ASSERT( a_rect ); + + if( a_node->IsInternalNode() ) // This is an internal node in the tree + { + for( int index = 0; index < a_node->m_count; ++index ) + { + if( Overlap( a_rect, &a_node->m_branch[index].m_rect ) ) + { + if( !Search( a_node->m_branch[index].m_child, a_rect, a_foundCount, + a_resultCallback, a_context ) ) + { + return false; // Don't continue searching + } + } + } + } + else // This is a leaf node + { + for( int index = 0; index < a_node->m_count; ++index ) + { + if( Overlap( a_rect, &a_node->m_branch[index].m_rect ) ) + { + DATATYPE& id = a_node->m_branch[index].m_data; + + // NOTE: There are different ways to return results. Here's where to modify + if( &a_resultCallback ) + { + ++a_foundCount; + + if( !a_resultCallback( id, a_context ) ) + { + return false; // Don't continue searching + } + } + } + } + } + + return true; // Continue searching +} + + + + +//calculate the minimum distance between a point and a rectangle as defined by Manolopoulos et al. +//it uses the square distance to avoid the use of ELEMTYPEREAL values, which are slower. +RTREE_TEMPLATE +ELEMTYPE RTREE_QUAL::MinDist( const ELEMTYPE a_point[NUMDIMS], Rect* a_rect ) +{ + ELEMTYPE *q, *s, *t; + q = (ELEMTYPE*) a_point; + s = a_rect->m_min; + t = a_rect->m_max; + int minDist = 0; + for( int index = 0; index < NUMDIMS; index++ ) + { + int r = q[index]; + if( q[index] < s[index] ) + { + r = s[index]; + } + else if( q[index] >t[index] ) + { + r = t[index]; + } + int addend = q[index] - r; + minDist += addend * addend; + } + return minDist; +} + + +//insert a NNNode in a list sorted by its minDist (desc.) +RTREE_TEMPLATE +void RTREE_QUAL::InsertNNListSorted( std::vector* nodeList, NNNode* newNode ) +{ + typedef typename std::vector::iterator iterator; + iterator iter = nodeList->begin(); + while( iter != nodeList->end() && (*iter)->minDist > newNode->minDist ) + { + ++iter; + } + nodeList->insert(iter, newNode); +} + + +#undef RTREE_TEMPLATE +#undef RTREE_QUAL +#undef RTREE_SEARCH_TEMPLATE +#undef RTREE_SEARCH_QUAL + +#endif // RTREE_H diff --git a/include/geometry/seg.h b/include/geometry/seg.h index 9989c76be7..7b5f64ac2f 100644 --- a/include/geometry/seg.h +++ b/include/geometry/seg.h @@ -267,6 +267,12 @@ class SEG { { return (a - b).EuclideanNorm(); } + + ecoord SquaredLength() const + { + return (a - b).SquaredEuclideanNorm(); + } + /** * Function Index() @@ -301,18 +307,7 @@ inline VECTOR2I SEG::LineProject( const VECTOR2I& aP ) const { // fixme: numerical errors for large integers assert(false); - /*const VECTOR2I d = aB - aA; - ecoord det = d.Dot(d); - ecoord dxdy = (ecoord) d.x * d.y; - - ecoord qx = - ( (extended_type) aA.x * d.y * d.y + (extended_type) d.x * d.x * x - dxdy * - (aA.y - y) ) / det; - extended_type qy = - ( (extended_type) aA.y * d.x * d.x + (extended_type) d.y * d.y * y - dxdy * - (aA.x - x) ) / det; - - return VECTOR2 ( (T) qx, (T) qy );*/ + return VECTOR2I(0, 0); } diff --git a/include/geometry/shape.h b/include/geometry/shape.h index 6d56fef79a..c759b0bbd6 100644 --- a/include/geometry/shape.h +++ b/include/geometry/shape.h @@ -77,7 +77,10 @@ class SHAPE { * Returns a dynamically allocated copy of the shape * @retval copy of the shape */ - virtual SHAPE* Clone() const { assert(false); }; + virtual SHAPE* Clone() const { + assert(false); + return NULL; + }; /** * Function Collide() diff --git a/include/geometry/shape_index.h b/include/geometry/shape_index.h index 4040cc2549..a36ecff97d 100644 --- a/include/geometry/shape_index.h +++ b/include/geometry/shape_index.h @@ -2,6 +2,7 @@ * This program source code file is part of KiCad, a free EDA CAD application. * * Copyright (C) 2013 CERN + * @author Jacobo Aragunde Pérez * @author Tomasz Wlostowski * * This program is free software; you can redistribute it and/or @@ -25,266 +26,328 @@ #ifndef __SHAPE_INDEX_H #define __SHAPE_INDEX_H -#include +#include +#include +#include -template const SHAPE *defaultShapeFunctor( const T aItem ) + +/** + * shapeFunctor template function + * + * It is used by SHAPE_INDEX to get a SHAPE* from another type. + * By default relies on T::GetShape() method, should be specialized if the T object + * doesn't allow that method. + * @param object generic T object + * @return a SHAPE* object equivalent to object. + */ +template +static const SHAPE* shapeFunctor( T aItem ) { - return aItem->GetShape(); + return aItem->GetShape(); } -template > +/** + * shapeFunctor template function: specialization for T = SHAPE* + */ +template<> +const SHAPE* shapeFunctor( SHAPE* aItem ) +{ + return aItem; +} -class SHAPE_INDEX_LIST { - - struct ShapeEntry { - ShapeEntry(T aParent) - { - shape = ShapeFunctor(aParent); - bbox = shape->BBox(0); - parent = aParent; - } +/** + * boundingBox template method + * + * It is used by SHAPE_INDEX to get the bounding box of a generic T object. + * By default relies on T::BBox() method, should be specialized if the T object + * doesn't allow that method. + * @param object generic T object + * @return a BOX2I object containing the bounding box of the T object. + */ +template +BOX2I boundingBox( T object ) +{ + return shapeFunctor(object)->BBox(); +} - ~ShapeEntry() - { - - } +/** + * acceptVisitor template method + * + * It is used by SHAPE_INDEX to implement Accept(). + * By default relies on V::operation() redefinition, should be specialized if V class + * doesn't have its () operation defined to accept T objects. + * @param object generic T object + * @param visitor V visitor object + */ +template +void acceptVisitor( T object, V visitor ) +{ + visitor(object); +} - T parent; - const SHAPE *shape; - BOX2I bbox; - }; +/** + * collide template method + * + * It is used by SHAPE_INDEX to implement Query(). + * By default relies on T::Collide(U) method, should be specialized if the T object + * doesn't allow that method. + * @param object generic T object + * @param anotherObject generic U object + * @param minDistance minimum collision distance + * @return if object and anotherObject collide + */ +template +bool collide( T object, U anotherObject, int minDistance ) +{ + return shapeFunctor(object)->Collide( anotherObject, minDistance ); +} - typedef std::vector ShapeVec; - typedef typename std::vector::iterator ShapeVecIter; - -public: +template +bool queryCallback(T shape, void* context) { + V* visitor = (V*) context; + acceptVisitor(shape, *visitor); + return true; +} -// "Normal" iterator interface, for STL algorithms. - class iterator { +template +class SHAPE_INDEX { - public: - iterator() {}; + public: - iterator( ShapeVecIter aCurrent) - : m_current(aCurrent) {}; + SHAPE_INDEX(); - iterator(const iterator &b) : - m_current(b.m_current) {}; + ~SHAPE_INDEX(); - T operator*() const - { - return (*m_current).parent; - } + /** + * Function Add() + * + * Adds a SHAPE to the index. + * @param shape the new SHAPE + */ + void Add( T shape ); - void operator++() - { - ++m_current; - } + /** + * Function Remove() + * + * Removes a SHAPE to the index. + * @param shape the new SHAPE + */ + void Remove( T shape ); - iterator& operator++(int dummy) - { - ++m_current; - return *this; - } + /** + * Function RemoveAll() + * + * Removes all the contents of the index. + */ + void RemoveAll(); - bool operator ==( const iterator& rhs ) const - { - return m_current == rhs.m_current; - } + /** + * Function Accept() + * + * Accepts a visitor for every SHAPE object contained in this INDEX. + * @param visitor Visitor object to be run + */ + template + void Accept( V visitor ) + { + SHAPE_INDEX::Iterator iter = this->Begin(); + while(!iter.IsNull()) { + T shape = *iter; + acceptVisitor(shape, visitor); + iter++; + } + } - bool operator !=( const iterator& rhs ) const - { - return m_current != rhs.m_current; - } + /** + * Function Reindex() + * + * Rebuilds the index. This should be used if the geometry of the objects + * contained by the index has changed. + */ + void Reindex(); - const iterator& operator=(const iterator& rhs) - { - m_current = rhs.m_current; - return *this; - } + /** + * Function Query() + * + * Runs a callback on every SHAPE object contained in the bounding box of (shape). + * @param shape shape to search against + * @param minDistance distance threshold + * @param visitor object to be invoked on every object contained in the search area. + */ - private: - ShapeVecIter m_current; - }; + template + int Query( const SHAPE *shape, int minDistance, V& visitor, bool aExact ) + { + BOX2I box = shape->BBox(); + box.Inflate(minDistance); + + int min[2] = {box.GetX(), box.GetY()}; + int max[2] = {box.GetRight(), box.GetBottom()}; -// "Query" iterator, for iterating over a set of spatially matching shapes. - class query_iterator { - public: + return this->m_tree->Search(min, max, visitor); + } - query_iterator() - { + class Iterator + { + private: - } + typedef typename RTree::Iterator RTreeIterator; + RTreeIterator iterator; - query_iterator( ShapeVecIter aCurrent, ShapeVecIter aEnd, SHAPE *aShape, int aMinDistance, bool aExact) - : m_end(aEnd), - m_current(aCurrent), - m_shape(aShape), - m_minDistance(aMinDistance), - m_exact(aExact) - { - if(aShape) - { - m_refBBox = aShape->BBox(); - next(); - } - } + /** + * Function Init() + * + * Setup the internal tree iterator. + * @param tree pointer to a RTREE object + */ + void Init(RTree* tree) { + tree->GetFirst(iterator); + } - query_iterator(const query_iterator &b) - : m_end(b.m_end), - m_current(b.m_current), - m_shape(b.m_shape), - m_minDistance(b.m_minDistance), - m_exact(b.m_exact), - m_refBBox(b.m_refBBox) - { + public: - } + /** + * Iterator constructor + * + * Creates an iterator for the index object + * @param index SHAPE_INDEX object to iterate + */ + Iterator(SHAPE_INDEX* index) { + Init(index->m_tree); + } - - T operator*() const - { - return (*m_current).parent; - } + /** + * Operator * (prefix) + * + * Returns the next data element. + */ + T operator*() { + return *iterator; + } - query_iterator& operator++() - { - ++m_current; - next(); - return *this; - } + /** + * Operator ++ (prefix) + * + * Shifts the iterator to the next element. + */ + bool operator++() { + return ++iterator; + } - query_iterator& operator++(int dummy) - { - ++m_current; - next(); - return *this; - } + /** + * Operator ++ (postfix) + * + * Shifts the iterator to the next element. + */ + bool operator++(int) { + return ++iterator; + } - bool operator ==( const query_iterator& rhs ) const - { - return m_current == rhs.m_current; - } + /** + * Function IsNull() + * + * Checks if the iterator has reached the end. + * @return true if it is in an invalid position (data finished) + */ + bool IsNull() { + return iterator.IsNull(); + } - bool operator !=( const query_iterator& rhs ) const - { - return m_current != rhs.m_current; - } + /** + * Function IsNotNull() + * + * Checks if the iterator has not reached the end. + * @return true if it is in an valid position (data not finished) + */ + bool IsNotNull() { + return iterator.IsNotNull(); + } - const query_iterator& operator=(const query_iterator& rhs) - { - m_end = rhs.m_end; - m_current = rhs.m_current; - m_shape = rhs.m_shape; - m_minDistance = rhs.m_minDistance; - m_exact = rhs.m_exact; - m_refBBox = rhs.m_refBBox; - return *this; - } + /** + * Function Next() + * + * Returns the current element of the iterator and moves to the next + * position. + * @return SHAPE object pointed by the iterator before moving to the + * next position. + */ + T Next() { + T object = *iterator; + ++iterator; + return object; + } + }; - private: + /** + * Function Begin() + * + * Creates an iterator for the current index object + * @return iterator + */ + Iterator Begin(); - void next() - { - while(m_current != m_end) - { - if (m_refBBox.Distance(m_current->bbox) <= m_minDistance) - { - if(!m_exact || m_current->shape->Collide(m_shape, m_minDistance)) - return; - } - ++m_current; - } - } + private: - ShapeVecIter m_end; - ShapeVecIter m_current; - BOX2I m_refBBox; - bool m_exact; - SHAPE *m_shape; - int m_minDistance; - }; - - void Add(T aItem) - { - ShapeEntry s (aItem); - - m_shapes.push_back(s); - } - - void Remove(const T aItem) - { - ShapeVecIter i; - - for(i=m_shapes.begin(); i!=m_shapes.end();++i) - { - if(i->parent == aItem) - break; - } - - if(i == m_shapes.end()) - return; - - m_shapes.erase(i); - } - - int Size() const - { - return m_shapes.size(); - } - - template - int Query( const SHAPE *aShape, int aMinDistance, Visitor &v, bool aExact = true) //const - { - ShapeVecIter i; - int n = 0; - VECTOR2I::extended_type minDistSq = (VECTOR2I::extended_type) aMinDistance * aMinDistance; - - BOX2I refBBox = aShape->BBox(); - - for(i = m_shapes.begin(); i!=m_shapes.end(); ++i) - { - if (refBBox.SquaredDistance(i->bbox) <= minDistSq) - { - if(!aExact || i->shape->Collide(aShape, aMinDistance)) - { - n++; - if(!v( i->parent )) - return n; - } - } - } - return n; - } - - void Clear() - { - m_shapes.clear(); - } - - query_iterator qbegin( SHAPE *aShape, int aMinDistance, bool aExact ) - { - return query_iterator( m_shapes.begin(), m_shapes.end(), aShape, aMinDistance, aExact); - } - - const query_iterator qend() - { - return query_iterator( m_shapes.end(), m_shapes.end(), NULL, 0, false ); - } - - iterator begin() - { - return iterator( m_shapes.begin() ); - } - - iterator end() - { - return iterator( m_shapes.end() ); - } - -private: - - ShapeVec m_shapes; + RTree* m_tree; }; +/* + * Class members implementation + */ + +template +SHAPE_INDEX::SHAPE_INDEX() { + this->m_tree = new RTree(); +} + +template +SHAPE_INDEX::~SHAPE_INDEX() { + delete this->m_tree; +} + +template +void SHAPE_INDEX::Add(T shape) { + BOX2I box = boundingBox(shape); + int min[2]= {box.GetX(), box.GetY()}; + int max[2] = {box.GetRight(), box.GetBottom()}; + this->m_tree->Insert(min, max, shape); +} + +template +void SHAPE_INDEX::Remove(T shape) { + BOX2I box = boundingBox(shape); + int min[2]= {box.GetX(), box.GetY()}; + int max[2] = {box.GetRight(), box.GetBottom()}; + this->m_tree->Remove(min, max, shape); +} + +template +void SHAPE_INDEX::RemoveAll() { + this->m_tree->RemoveAll(); +} + +template +void SHAPE_INDEX::Reindex() { + RTree* newTree; + newTree = new RTree(); + + SHAPE_INDEX::Iterator iter = this->Begin(); + while(!iter.IsNull()) { + T shape = *iter; + BOX2I box = boundingBox(shape); + int min[2]= {box.GetX(), box.GetY()}; + int max[2] = {box.GetRight(), box.GetBottom()}; + newTree->Insert(min, max, shape); + iter++; + } + delete this->m_tree; + this->m_tree = newTree; +} + +template +typename SHAPE_INDEX::Iterator SHAPE_INDEX::Begin() { + return Iterator(this); +} + + #endif diff --git a/include/geometry/shape_index_list.h b/include/geometry/shape_index_list.h new file mode 100644 index 0000000000..47c5d98c97 --- /dev/null +++ b/include/geometry/shape_index_list.h @@ -0,0 +1,290 @@ +/* + * This program source code file is part of KiCad, a free EDA CAD application. + * + * Copyright (C) 2013 CERN + * @author Tomasz Wlostowski + * + * This program 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 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, you may find one here: + * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html + * or you may search the http://www.gnu.org website for the version 2 license, + * or you may write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA + */ + +#ifndef __SHAPE_INDEX_LIST_H +#define __SHAPE_INDEX_LIST_H + +#include + +template const SHAPE *defaultShapeFunctor( const T aItem ) +{ + return aItem->GetShape(); +} + +template > + +class SHAPE_INDEX_LIST { + + struct ShapeEntry { + ShapeEntry(T aParent) + { + shape = ShapeFunctor(aParent); + bbox = shape->BBox(0); + parent = aParent; + } + + ~ShapeEntry() + { + + } + + T parent; + const SHAPE *shape; + BOX2I bbox; + }; + + typedef std::vector ShapeVec; + typedef typename std::vector::iterator ShapeVecIter; + +public: + +// "Normal" iterator interface, for STL algorithms. + class iterator { + + public: + iterator() {}; + + iterator( ShapeVecIter aCurrent) + : m_current(aCurrent) {}; + + iterator(const iterator &b) : + m_current(b.m_current) {}; + + T operator*() const + { + return (*m_current).parent; + } + + void operator++() + { + ++m_current; + } + + iterator& operator++(int dummy) + { + ++m_current; + return *this; + } + + bool operator ==( const iterator& rhs ) const + { + return m_current == rhs.m_current; + } + + bool operator !=( const iterator& rhs ) const + { + return m_current != rhs.m_current; + } + + const iterator& operator=(const iterator& rhs) + { + m_current = rhs.m_current; + return *this; + } + + private: + ShapeVecIter m_current; + }; + +// "Query" iterator, for iterating over a set of spatially matching shapes. + class query_iterator { + public: + + query_iterator() + { + + } + + query_iterator( ShapeVecIter aCurrent, ShapeVecIter aEnd, SHAPE *aShape, int aMinDistance, bool aExact) + : m_end(aEnd), + m_current(aCurrent), + m_shape(aShape), + m_minDistance(aMinDistance), + m_exact(aExact) + { + if(aShape) + { + m_refBBox = aShape->BBox(); + next(); + } + } + + query_iterator(const query_iterator &b) + : m_end(b.m_end), + m_current(b.m_current), + m_shape(b.m_shape), + m_minDistance(b.m_minDistance), + m_exact(b.m_exact), + m_refBBox(b.m_refBBox) + { + + } + + + T operator*() const + { + return (*m_current).parent; + } + + query_iterator& operator++() + { + ++m_current; + next(); + return *this; + } + + query_iterator& operator++(int dummy) + { + ++m_current; + next(); + return *this; + } + + bool operator ==( const query_iterator& rhs ) const + { + return m_current == rhs.m_current; + } + + bool operator !=( const query_iterator& rhs ) const + { + return m_current != rhs.m_current; + } + + const query_iterator& operator=(const query_iterator& rhs) + { + m_end = rhs.m_end; + m_current = rhs.m_current; + m_shape = rhs.m_shape; + m_minDistance = rhs.m_minDistance; + m_exact = rhs.m_exact; + m_refBBox = rhs.m_refBBox; + return *this; + } + + private: + + void next() + { + while(m_current != m_end) + { + if (m_refBBox.Distance(m_current->bbox) <= m_minDistance) + { + if(!m_exact || m_current->shape->Collide(m_shape, m_minDistance)) + return; + } + ++m_current; + } + } + + ShapeVecIter m_end; + ShapeVecIter m_current; + BOX2I m_refBBox; + bool m_exact; + SHAPE *m_shape; + int m_minDistance; + }; + + void Add(T aItem) + { + ShapeEntry s (aItem); + + m_shapes.push_back(s); + } + + void Remove(const T aItem) + { + ShapeVecIter i; + + for(i=m_shapes.begin(); i!=m_shapes.end();++i) + { + if(i->parent == aItem) + break; + } + + if(i == m_shapes.end()) + return; + + m_shapes.erase(i); + } + + int Size() const + { + return m_shapes.size(); + } + + template + int Query( const SHAPE *aShape, int aMinDistance, Visitor &v, bool aExact = true) //const + { + ShapeVecIter i; + int n = 0; + VECTOR2I::extended_type minDistSq = (VECTOR2I::extended_type) aMinDistance * aMinDistance; + + BOX2I refBBox = aShape->BBox(); + + for(i = m_shapes.begin(); i!=m_shapes.end(); ++i) + { + if (refBBox.SquaredDistance(i->bbox) <= minDistSq) + { + if(!aExact || i->shape->Collide(aShape, aMinDistance)) + { + n++; + if(!v( i->parent )) + return n; + } + } + } + return n; + } + + void Clear() + { + m_shapes.clear(); + } + + query_iterator qbegin( SHAPE *aShape, int aMinDistance, bool aExact ) + { + return query_iterator( m_shapes.begin(), m_shapes.end(), aShape, aMinDistance, aExact); + } + + const query_iterator qend() + { + return query_iterator( m_shapes.end(), m_shapes.end(), NULL, 0, false ); + } + + iterator begin() + { + return iterator( m_shapes.begin() ); + } + + iterator end() + { + return iterator( m_shapes.end() ); + } + +private: + + ShapeVec m_shapes; +}; + +#endif diff --git a/include/view/view.h b/include/view/view.h index 5013fdaf91..714283ef44 100644 --- a/include/view/view.h +++ b/include/view/view.h @@ -362,6 +362,8 @@ public: */ void EnableTopLayer( bool aEnable ); + int GetTopLayer() const; + /** * Function ClearTopLayers() * Removes all layers from the on-the-top set (they are no longer displayed over the rest of diff --git a/include/view/view_controls.h b/include/view/view_controls.h index ca539f467c..7bf11860cb 100644 --- a/include/view/view_controls.h +++ b/include/view/view_controls.h @@ -46,9 +46,9 @@ class VIEW; class VIEW_CONTROLS { public: - VIEW_CONTROLS( VIEW* aView ) : m_view( aView ), m_snappingEnabled( false ), - m_grabMouse( false ), m_autoPanEnabled( false ), m_autoPanMargin( 0.1 ), - m_autoPanSpeed( 0.15 ) {}; + VIEW_CONTROLS( VIEW* aView ) : m_view( aView ), m_forceCursorPosition( false ), + m_snappingEnabled( false ), m_grabMouse( false ), m_autoPanEnabled( false ), + m_autoPanMargin( 0.1 ), m_autoPanSpeed( 0.15 ) {}; virtual ~VIEW_CONTROLS() {}; /** @@ -121,6 +121,21 @@ public: */ virtual const VECTOR2D GetCursorPosition() const = 0; + + /** + * Function ForceCursorPosition() + * Places the cursor immediately at a given point. Mouse movement is ignored. + * @param aEnabled enable forced cursor position + * @param aPosition the position + */ + virtual void ForceCursorPosition( bool aEnabled, const VECTOR2D& aPosition = VECTOR2D(0, 0) ) + { + m_forcedPosition = aPosition; + m_forceCursorPosition = aEnabled; + } + + virtual void ShowCursor( bool aEnabled ); + protected: /// Pointer to controlled VIEW. VIEW* m_view; @@ -128,6 +143,15 @@ protected: /// Current mouse position VECTOR2D m_mousePosition; + /// Current cursor position + VECTOR2D m_cursorPosition; + + /// Forced cursor position + VECTOR2D m_forcedPosition; + + /// Is the forced cursor position enabled + bool m_forceCursorPosition; + /// Should the cursor snap to grid or move freely bool m_snappingEnabled; diff --git a/include/view/view_group.h b/include/view/view_group.h index 18f5f5a9b6..7664060560 100644 --- a/include/view/view_group.h +++ b/include/view/view_group.h @@ -142,7 +142,7 @@ public: * * @return Pointer to the VIEW instance. */ - KiGfx::VIEW *GetView() const + KiGfx::VIEW* GetView() const { return m_view; } diff --git a/include/wxPcbStruct.h b/include/wxPcbStruct.h index 2110eb75d9..de4c55d0b2 100644 --- a/include/wxPcbStruct.h +++ b/include/wxPcbStruct.h @@ -1690,6 +1690,11 @@ public: */ void UpdateTitle(); + void SetTopLayer( LAYER_NUM aLayer ) + { + setTopLayer( aLayer ); + } + DECLARE_EVENT_TABLE() }; diff --git a/pcbnew/pcb_painter.cpp b/pcbnew/pcb_painter.cpp index a2cae1c427..5fa4c10e2e 100644 --- a/pcbnew/pcb_painter.cpp +++ b/pcbnew/pcb_painter.cpp @@ -206,8 +206,6 @@ PCB_PAINTER::PCB_PAINTER( GAL* aGal ) : bool PCB_PAINTER::Draw( const VIEW_ITEM* aItem, int aLayer ) { - const BOARD_ITEM* item = static_cast( aItem ); - // the "cast" applied in here clarifies which overloaded draw() is called switch( item->Type() ) { @@ -284,9 +282,9 @@ void PCB_PAINTER::draw( const TRACK* aTrack, int aLayer ) return; NETINFO_ITEM* net = ( (BOARD*) aTrack->GetParent() )->FindNet( netNumber ); - if( net == NULL ) + if( !net ) return; - + std::string netName = std::string( net->GetShortNetname().mb_str() ); VECTOR2D textPosition = start + line / 2.0; // center of the track double textOrientation = -atan( line.y / line.x ); diff --git a/polygon/CMakeLists.txt b/polygon/CMakeLists.txt index 52068b52af..482c3a1170 100644 --- a/polygon/CMakeLists.txt +++ b/polygon/CMakeLists.txt @@ -9,6 +9,12 @@ set(POLYGON_SRCS PolyLine.cpp polygon_test_point_inside.cpp clipper.cpp - ) + + poly2tri/common/shapes.cc + poly2tri/sweep/sweep.cc + poly2tri/sweep/cdt.cc + poly2tri/sweep/advancing_front.cc + poly2tri/sweep/sweep_context.cc +) add_library(polygon STATIC ${POLYGON_SRCS}) diff --git a/polygon/poly2tri/common/shapes.cc b/polygon/poly2tri/common/shapes.cc new file mode 100644 index 0000000000..77bafa1501 --- /dev/null +++ b/polygon/poly2tri/common/shapes.cc @@ -0,0 +1,367 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ +#include "shapes.h" +#include + +namespace p2t { + +Triangle::Triangle(Point& a, Point& b, Point& c) +{ + points_[0] = &a; points_[1] = &b; points_[2] = &c; + neighbors_[0] = NULL; neighbors_[1] = NULL; neighbors_[2] = NULL; + constrained_edge[0] = constrained_edge[1] = constrained_edge[2] = false; + delaunay_edge[0] = delaunay_edge[1] = delaunay_edge[2] = false; + interior_ = false; +} + +// Update neighbor pointers +void Triangle::MarkNeighbor(Point* p1, Point* p2, Triangle* t) +{ + if ((p1 == points_[2] && p2 == points_[1]) || (p1 == points_[1] && p2 == points_[2])) + neighbors_[0] = t; + else if ((p1 == points_[0] && p2 == points_[2]) || (p1 == points_[2] && p2 == points_[0])) + neighbors_[1] = t; + else if ((p1 == points_[0] && p2 == points_[1]) || (p1 == points_[1] && p2 == points_[0])) + neighbors_[2] = t; + else + assert(0); +} + +// Exhaustive search to update neighbor pointers +void Triangle::MarkNeighbor(Triangle& t) +{ + if (t.Contains(points_[1], points_[2])) { + neighbors_[0] = &t; + t.MarkNeighbor(points_[1], points_[2], this); + } else if (t.Contains(points_[0], points_[2])) { + neighbors_[1] = &t; + t.MarkNeighbor(points_[0], points_[2], this); + } else if (t.Contains(points_[0], points_[1])) { + neighbors_[2] = &t; + t.MarkNeighbor(points_[0], points_[1], this); + } +} + +/** + * Clears all references to all other triangles and points + */ +void Triangle::Clear() +{ + Triangle *t; + for( int i=0; i<3; i++ ) + { + t = neighbors_[i]; + if( t != NULL ) + { + t->ClearNeighbor( this ); + } + } + ClearNeighbors(); + points_[0]=points_[1]=points_[2] = NULL; +} + +void Triangle::ClearNeighbor(Triangle *triangle ) +{ + if( neighbors_[0] == triangle ) + { + neighbors_[0] = NULL; + } + else if( neighbors_[1] == triangle ) + { + neighbors_[1] = NULL; + } + else + { + neighbors_[2] = NULL; + } +} + +void Triangle::ClearNeighbors() +{ + neighbors_[0] = NULL; + neighbors_[1] = NULL; + neighbors_[2] = NULL; +} + +void Triangle::ClearDelunayEdges() +{ + delaunay_edge[0] = delaunay_edge[1] = delaunay_edge[2] = false; +} + +Point* Triangle::OppositePoint(Triangle& t, Point& p) +{ + Point *cw = t.PointCW(p); + double x = cw->x; + double y = cw->y; + x = p.x; + y = p.y; + return PointCW(*cw); +} + +// Legalized triangle by rotating clockwise around point(0) +void Triangle::Legalize(Point& point) +{ + points_[1] = points_[0]; + points_[0] = points_[2]; + points_[2] = &point; +} + +// Legalize triagnle by rotating clockwise around oPoint +void Triangle::Legalize(Point& opoint, Point& npoint) +{ + if (&opoint == points_[0]) { + points_[1] = points_[0]; + points_[0] = points_[2]; + points_[2] = &npoint; + } else if (&opoint == points_[1]) { + points_[2] = points_[1]; + points_[1] = points_[0]; + points_[0] = &npoint; + } else if (&opoint == points_[2]) { + points_[0] = points_[2]; + points_[2] = points_[1]; + points_[1] = &npoint; + } else { + assert(0); + } +} + +int Triangle::Index(const Point* p) +{ + if (p == points_[0]) { + return 0; + } else if (p == points_[1]) { + return 1; + } else if (p == points_[2]) { + return 2; + } + assert(0); +} + +int Triangle::EdgeIndex(const Point* p1, const Point* p2) +{ + if (points_[0] == p1) { + if (points_[1] == p2) { + return 2; + } else if (points_[2] == p2) { + return 1; + } + } else if (points_[1] == p1) { + if (points_[2] == p2) { + return 0; + } else if (points_[0] == p2) { + return 2; + } + } else if (points_[2] == p1) { + if (points_[0] == p2) { + return 1; + } else if (points_[1] == p2) { + return 0; + } + } + return -1; +} + +void Triangle::MarkConstrainedEdge(const int index) +{ + constrained_edge[index] = true; +} + +void Triangle::MarkConstrainedEdge(Edge& edge) +{ + MarkConstrainedEdge(edge.p, edge.q); +} + +// Mark edge as constrained +void Triangle::MarkConstrainedEdge(Point* p, Point* q) +{ + if ((q == points_[0] && p == points_[1]) || (q == points_[1] && p == points_[0])) { + constrained_edge[2] = true; + } else if ((q == points_[0] && p == points_[2]) || (q == points_[2] && p == points_[0])) { + constrained_edge[1] = true; + } else if ((q == points_[1] && p == points_[2]) || (q == points_[2] && p == points_[1])) { + constrained_edge[0] = true; + } +} + +// The point counter-clockwise to given point +Point* Triangle::PointCW(Point& point) +{ + if (&point == points_[0]) { + return points_[2]; + } else if (&point == points_[1]) { + return points_[0]; + } else if (&point == points_[2]) { + return points_[1]; + } + assert(0); +} + +// The point counter-clockwise to given point +Point* Triangle::PointCCW(Point& point) +{ + if (&point == points_[0]) { + return points_[1]; + } else if (&point == points_[1]) { + return points_[2]; + } else if (&point == points_[2]) { + return points_[0]; + } + assert(0); +} + +// The neighbor clockwise to given point +Triangle* Triangle::NeighborCW(Point& point) +{ + if (&point == points_[0]) { + return neighbors_[1]; + } else if (&point == points_[1]) { + return neighbors_[2]; + } + return neighbors_[0]; +} + +// The neighbor counter-clockwise to given point +Triangle* Triangle::NeighborCCW(Point& point) +{ + if (&point == points_[0]) { + return neighbors_[2]; + } else if (&point == points_[1]) { + return neighbors_[0]; + } + return neighbors_[1]; +} + +bool Triangle::GetConstrainedEdgeCCW(Point& p) +{ + if (&p == points_[0]) { + return constrained_edge[2]; + } else if (&p == points_[1]) { + return constrained_edge[0]; + } + return constrained_edge[1]; +} + +bool Triangle::GetConstrainedEdgeCW(Point& p) +{ + if (&p == points_[0]) { + return constrained_edge[1]; + } else if (&p == points_[1]) { + return constrained_edge[2]; + } + return constrained_edge[0]; +} + +void Triangle::SetConstrainedEdgeCCW(Point& p, bool ce) +{ + if (&p == points_[0]) { + constrained_edge[2] = ce; + } else if (&p == points_[1]) { + constrained_edge[0] = ce; + } else { + constrained_edge[1] = ce; + } +} + +void Triangle::SetConstrainedEdgeCW(Point& p, bool ce) +{ + if (&p == points_[0]) { + constrained_edge[1] = ce; + } else if (&p == points_[1]) { + constrained_edge[2] = ce; + } else { + constrained_edge[0] = ce; + } +} + +bool Triangle::GetDelunayEdgeCCW(Point& p) +{ + if (&p == points_[0]) { + return delaunay_edge[2]; + } else if (&p == points_[1]) { + return delaunay_edge[0]; + } + return delaunay_edge[1]; +} + +bool Triangle::GetDelunayEdgeCW(Point& p) +{ + if (&p == points_[0]) { + return delaunay_edge[1]; + } else if (&p == points_[1]) { + return delaunay_edge[2]; + } + return delaunay_edge[0]; +} + +void Triangle::SetDelunayEdgeCCW(Point& p, bool e) +{ + if (&p == points_[0]) { + delaunay_edge[2] = e; + } else if (&p == points_[1]) { + delaunay_edge[0] = e; + } else { + delaunay_edge[1] = e; + } +} + +void Triangle::SetDelunayEdgeCW(Point& p, bool e) +{ + if (&p == points_[0]) { + delaunay_edge[1] = e; + } else if (&p == points_[1]) { + delaunay_edge[2] = e; + } else { + delaunay_edge[0] = e; + } +} + +// The neighbor across to given point +Triangle& Triangle::NeighborAcross(Point& opoint) +{ + if (&opoint == points_[0]) { + return *neighbors_[0]; + } else if (&opoint == points_[1]) { + return *neighbors_[1]; + } + return *neighbors_[2]; +} + +void Triangle::DebugPrint() +{ + using namespace std; + cout << points_[0]->x << "," << points_[0]->y << " "; + cout << points_[1]->x << "," << points_[1]->y << " "; + cout << points_[2]->x << "," << points_[2]->y << endl; +} + +} + diff --git a/polygon/poly2tri/common/shapes.h b/polygon/poly2tri/common/shapes.h new file mode 100644 index 0000000000..8c19d91ef8 --- /dev/null +++ b/polygon/poly2tri/common/shapes.h @@ -0,0 +1,325 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +// Include guard +#ifndef SHAPES_H +#define SHAPES_H + +#include +#include +#include +#include + +namespace p2t { + +struct Edge; + +struct Point { + + double x, y; + + /// Default constructor does nothing (for performance). + Point() + { + x = 0.0; + y = 0.0; + } + + /// The edges this point constitutes an upper ending point + std::vector edge_list; + + /// Construct using coordinates. + Point(double x, double y) : x(x), y(y) {} + + /// Set this point to all zeros. + void set_zero() + { + x = 0.0; + y = 0.0; + } + + /// Set this point to some specified coordinates. + void set(double x_, double y_) + { + x = x_; + y = y_; + } + + /// Negate this point. + Point operator -() const + { + Point v; + v.set(-x, -y); + return v; + } + + /// Add a point to this point. + void operator +=(const Point& v) + { + x += v.x; + y += v.y; + } + + /// Subtract a point from this point. + void operator -=(const Point& v) + { + x -= v.x; + y -= v.y; + } + + /// Multiply this point by a scalar. + void operator *=(double a) + { + x *= a; + y *= a; + } + + /// Get the length of this point (the norm). + double Length() const + { + return sqrt(x * x + y * y); + } + + /// Convert this point into a unit point. Returns the Length. + double Normalize() + { + double len = Length(); + x /= len; + y /= len; + return len; + } + +}; + +// Represents a simple polygon's edge +struct Edge { + + Point* p, *q; + + /// Constructor + Edge(Point& p1, Point& p2) : p(&p1), q(&p2) + { + if (p1.y > p2.y) { + q = &p1; + p = &p2; + } else if (p1.y == p2.y) { + if (p1.x > p2.x) { + q = &p1; + p = &p2; + } else if (p1.x == p2.x) { + // Repeat points + assert(false); + } + } + + q->edge_list.push_back(this); + } +}; + +// Triangle-based data structures are know to have better performance than quad-edge structures +// See: J. Shewchuk, "Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator" +// "Triangulations in CGAL" +class Triangle { +public: + +/// Constructor +Triangle(Point& a, Point& b, Point& c); + +/// Flags to determine if an edge is a Constrained edge +bool constrained_edge[3]; +/// Flags to determine if an edge is a Delauney edge +bool delaunay_edge[3]; + +Point* GetPoint(const int& index); +Point* PointCW(Point& point); +Point* PointCCW(Point& point); +Point* OppositePoint(Triangle& t, Point& p); + +Triangle* GetNeighbor(const int& index); +void MarkNeighbor(Point* p1, Point* p2, Triangle* t); +void MarkNeighbor(Triangle& t); + +void MarkConstrainedEdge(const int index); +void MarkConstrainedEdge(Edge& edge); +void MarkConstrainedEdge(Point* p, Point* q); + +int Index(const Point* p); +int EdgeIndex(const Point* p1, const Point* p2); + +Triangle* NeighborCW(Point& point); +Triangle* NeighborCCW(Point& point); +bool GetConstrainedEdgeCCW(Point& p); +bool GetConstrainedEdgeCW(Point& p); +void SetConstrainedEdgeCCW(Point& p, bool ce); +void SetConstrainedEdgeCW(Point& p, bool ce); +bool GetDelunayEdgeCCW(Point& p); +bool GetDelunayEdgeCW(Point& p); +void SetDelunayEdgeCCW(Point& p, bool e); +void SetDelunayEdgeCW(Point& p, bool e); + +bool Contains(Point* p); +bool Contains(const Edge& e); +bool Contains(Point* p, Point* q); +void Legalize(Point& point); +void Legalize(Point& opoint, Point& npoint); +/** + * Clears all references to all other triangles and points + */ +void Clear(); +void ClearNeighbor(Triangle *triangle ); +void ClearNeighbors(); +void ClearDelunayEdges(); + +inline bool IsInterior(); +inline void IsInterior(bool b); + +Triangle& NeighborAcross(Point& opoint); + +void DebugPrint(); + +private: + +/// Triangle points +Point* points_[3]; +/// Neighbor list +Triangle* neighbors_[3]; + +/// Has this triangle been marked as an interior triangle? +bool interior_; +}; + +inline bool cmp(const Point* a, const Point* b) +{ + if (a->y < b->y) { + return true; + } else if (a->y == b->y) { + // Make sure q is point with greater x value + if (a->x < b->x) { + return true; + } + } + return false; +} + +/// Add two points_ component-wise. +inline Point operator +(const Point& a, const Point& b) +{ + return Point(a.x + b.x, a.y + b.y); +} + +/// Subtract two points_ component-wise. +inline Point operator -(const Point& a, const Point& b) +{ + return Point(a.x - b.x, a.y - b.y); +} + +/// Multiply point by scalar +inline Point operator *(double s, const Point& a) +{ + return Point(s * a.x, s * a.y); +} + +inline bool operator ==(const Point& a, const Point& b) +{ + return a.x == b.x && a.y == b.y; +} + +inline bool operator !=(const Point& a, const Point& b) +{ + return !(a.x == b.x) && !(a.y == b.y); +} + +/// Peform the dot product on two vectors. +inline double Dot(const Point& a, const Point& b) +{ + return a.x * b.x + a.y * b.y; +} + +/// Perform the cross product on two vectors. In 2D this produces a scalar. +inline double Cross(const Point& a, const Point& b) +{ + return a.x * b.y - a.y * b.x; +} + +/// Perform the cross product on a point and a scalar. In 2D this produces +/// a point. +inline Point Cross(const Point& a, double s) +{ + return Point(s * a.y, -s * a.x); +} + +/// Perform the cross product on a scalar and a point. In 2D this produces +/// a point. +inline Point Cross(const double s, const Point& a) +{ + return Point(-s * a.y, s * a.x); +} + +inline Point* Triangle::GetPoint(const int& index) +{ + return points_[index]; +} + +inline Triangle* Triangle::GetNeighbor(const int& index) +{ + return neighbors_[index]; +} + +inline bool Triangle::Contains(Point* p) +{ + return p == points_[0] || p == points_[1] || p == points_[2]; +} + +inline bool Triangle::Contains(const Edge& e) +{ + return Contains(e.p) && Contains(e.q); +} + +inline bool Triangle::Contains(Point* p, Point* q) +{ + return Contains(p) && Contains(q); +} + +inline bool Triangle::IsInterior() +{ + return interior_; +} + +inline void Triangle::IsInterior(bool b) +{ + interior_ = b; +} + +} + +#endif + + diff --git a/polygon/poly2tri/common/utils.h b/polygon/poly2tri/common/utils.h new file mode 100644 index 0000000000..78416f2f31 --- /dev/null +++ b/polygon/poly2tri/common/utils.h @@ -0,0 +1,123 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef UTILS_H +#define UTILS_H + +// Otherwise #defines like M_PI are undeclared under Visual Studio +#define _USE_MATH_DEFINES + +#include +#include + +namespace p2t { + +const double PI_3div4 = 3 * M_PI / 4; +const double PI_div2 = 1.57079632679489661923; +const double EPSILON = 1e-12; + +enum Orientation { CW, CCW, COLLINEAR }; + +/** + * Forumla to calculate signed area
+ * Positive if CCW
+ * Negative if CW
+ * 0 if collinear
+ * 
+ * A[P1,P2,P3]  =  (x1*y2 - y1*x2) + (x2*y3 - y2*x3) + (x3*y1 - y3*x1)
+ *              =  (x1-x3)*(y2-y3) - (y1-y3)*(x2-x3)
+ *
+ */ +Orientation Orient2d(Point& pa, Point& pb, Point& pc) +{ + double detleft = (pa.x - pc.x) * (pb.y - pc.y); + double detright = (pa.y - pc.y) * (pb.x - pc.x); + double val = detleft - detright; + if (val > -EPSILON && val < EPSILON) { + return COLLINEAR; + } else if (val > 0) { + return CCW; + } + return CW; +} + +/* +bool InScanArea(Point& pa, Point& pb, Point& pc, Point& pd) +{ + double pdx = pd.x; + double pdy = pd.y; + double adx = pa.x - pdx; + double ady = pa.y - pdy; + double bdx = pb.x - pdx; + double bdy = pb.y - pdy; + + double adxbdy = adx * bdy; + double bdxady = bdx * ady; + double oabd = adxbdy - bdxady; + + if (oabd <= EPSILON) { + return false; + } + + double cdx = pc.x - pdx; + double cdy = pc.y - pdy; + + double cdxady = cdx * ady; + double adxcdy = adx * cdy; + double ocad = cdxady - adxcdy; + + if (ocad <= EPSILON) { + return false; + } + + return true; +} + +*/ + +bool InScanArea(Point& pa, Point& pb, Point& pc, Point& pd) +{ + double oadb = (pa.x - pb.x)*(pd.y - pb.y) - (pd.x - pb.x)*(pa.y - pb.y); + if (oadb >= -EPSILON) { + return false; + } + + double oadc = (pa.x - pc.x)*(pd.y - pc.y) - (pd.x - pc.x)*(pa.y - pc.y); + if (oadc <= EPSILON) { + return false; + } + return true; +} + +} + +#endif + diff --git a/polygon/poly2tri/poly2tri.h b/polygon/poly2tri/poly2tri.h new file mode 100644 index 0000000000..487755e2e9 --- /dev/null +++ b/polygon/poly2tri/poly2tri.h @@ -0,0 +1,39 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef POLY2TRI_H +#define POLY2TRI_H + +#include "common/shapes.h" +#include "sweep/cdt.h" + +#endif + diff --git a/polygon/poly2tri/sweep/advancing_front.cc b/polygon/poly2tri/sweep/advancing_front.cc new file mode 100644 index 0000000000..019df4a6eb --- /dev/null +++ b/polygon/poly2tri/sweep/advancing_front.cc @@ -0,0 +1,109 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ +#include "advancing_front.h" + +namespace p2t { + +AdvancingFront::AdvancingFront(Node& head, Node& tail) +{ + head_ = &head; + tail_ = &tail; + search_node_ = &head; +} + +Node* AdvancingFront::LocateNode(const double& x) +{ + Node* node = search_node_; + + if (x < node->value) { + while ((node = node->prev) != NULL) { + if (x >= node->value) { + search_node_ = node; + return node; + } + } + } else { + while ((node = node->next) != NULL) { + if (x < node->value) { + search_node_ = node->prev; + return node->prev; + } + } + } + return NULL; +} + +Node* AdvancingFront::FindSearchNode(const double& x) +{ + (void)x; // suppress compiler warnings "unused parameter 'x'" + // TODO: implement BST index + return search_node_; +} + +Node* AdvancingFront::LocatePoint(const Point* point) +{ + const double px = point->x; + Node* node = FindSearchNode(px); + const double nx = node->point->x; + + if (px == nx) { + if (point != node->point) { + // We might have two nodes with same x value for a short time + if (point == node->prev->point) { + node = node->prev; + } else if (point == node->next->point) { + node = node->next; + } else { + assert(0); + } + } + } else if (px < nx) { + while ((node = node->prev) != NULL) { + if (point == node->point) { + break; + } + } + } else { + while ((node = node->next) != NULL) { + if (point == node->point) + break; + } + } + if(node) search_node_ = node; + return node; +} + +AdvancingFront::~AdvancingFront() +{ +} + +} + diff --git a/polygon/poly2tri/sweep/advancing_front.h b/polygon/poly2tri/sweep/advancing_front.h new file mode 100644 index 0000000000..bab73d449c --- /dev/null +++ b/polygon/poly2tri/sweep/advancing_front.h @@ -0,0 +1,118 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef ADVANCED_FRONT_H +#define ADVANCED_FRONT_H + +#include "../common/shapes.h" + +namespace p2t { + +struct Node; + +// Advancing front node +struct Node { + Point* point; + Triangle* triangle; + + Node* next; + Node* prev; + + double value; + + Node(Point& p) : point(&p), triangle(NULL), next(NULL), prev(NULL), value(p.x) + { + } + + Node(Point& p, Triangle& t) : point(&p), triangle(&t), next(NULL), prev(NULL), value(p.x) + { + } + +}; + +// Advancing front +class AdvancingFront { +public: + +AdvancingFront(Node& head, Node& tail); +// Destructor +~AdvancingFront(); + +Node* head(); +void set_head(Node* node); +Node* tail(); +void set_tail(Node* node); +Node* search(); +void set_search(Node* node); + +/// Locate insertion point along advancing front +Node* LocateNode(const double& x); + +Node* LocatePoint(const Point* point); + +private: + +Node* head_, *tail_, *search_node_; + +Node* FindSearchNode(const double& x); +}; + +inline Node* AdvancingFront::head() +{ + return head_; +} +inline void AdvancingFront::set_head(Node* node) +{ + head_ = node; +} + +inline Node* AdvancingFront::tail() +{ + return tail_; +} +inline void AdvancingFront::set_tail(Node* node) +{ + tail_ = node; +} + +inline Node* AdvancingFront::search() +{ + return search_node_; +} + +inline void AdvancingFront::set_search(Node* node) +{ + search_node_ = node; +} + +} + +#endif diff --git a/polygon/poly2tri/sweep/cdt.cc b/polygon/poly2tri/sweep/cdt.cc new file mode 100644 index 0000000000..d7838257c7 --- /dev/null +++ b/polygon/poly2tri/sweep/cdt.cc @@ -0,0 +1,72 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ +#include "cdt.h" + +namespace p2t { + +CDT::CDT(std::vector polyline) +{ + sweep_context_ = new SweepContext(polyline); + sweep_ = new Sweep; +} + +void CDT::AddHole(std::vector polyline) +{ + sweep_context_->AddHole(polyline); +} + +void CDT::AddPoint(Point* point) { + sweep_context_->AddPoint(point); +} + +void CDT::Triangulate() +{ + sweep_->Triangulate(*sweep_context_); +} + +std::vector CDT::GetTriangles() +{ + return sweep_context_->GetTriangles(); +} + +std::list CDT::GetMap() +{ + return sweep_context_->GetMap(); +} + +CDT::~CDT() +{ + delete sweep_context_; + delete sweep_; +} + +} + diff --git a/polygon/poly2tri/sweep/cdt.h b/polygon/poly2tri/sweep/cdt.h new file mode 100644 index 0000000000..3e6f024086 --- /dev/null +++ b/polygon/poly2tri/sweep/cdt.h @@ -0,0 +1,105 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef CDT_H +#define CDT_H + +#include "advancing_front.h" +#include "sweep_context.h" +#include "sweep.h" + +/** + * + * @author Mason Green + * + */ + +namespace p2t { + +class CDT +{ +public: + + /** + * Constructor - add polyline with non repeating points + * + * @param polyline + */ + CDT(std::vector polyline); + + /** + * Destructor - clean up memory + */ + ~CDT(); + + /** + * Add a hole + * + * @param polyline + */ + void AddHole(std::vector polyline); + + /** + * Add a steiner point + * + * @param point + */ + void AddPoint(Point* point); + + /** + * Triangulate - do this AFTER you've added the polyline, holes, and Steiner points + */ + void Triangulate(); + + /** + * Get CDT triangles + */ + std::vector GetTriangles(); + + /** + * Get triangle map + */ + std::list GetMap(); + + private: + + /** + * Internals + */ + + SweepContext* sweep_context_; + Sweep* sweep_; + +}; + +} + +#endif diff --git a/polygon/poly2tri/sweep/sweep.cc b/polygon/poly2tri/sweep/sweep.cc new file mode 100644 index 0000000000..258df5db08 --- /dev/null +++ b/polygon/poly2tri/sweep/sweep.cc @@ -0,0 +1,813 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ +#include +#include "sweep.h" +#include "sweep_context.h" +#include "advancing_front.h" +#include "../common/utils.h" + +namespace p2t { + +// Triangulate simple polygon with holes +void Sweep::Triangulate(SweepContext& tcx) +{ + tcx.InitTriangulation(); + tcx.CreateAdvancingFront(nodes_); + // Sweep points; build mesh + SweepPoints(tcx); + // Clean up + FinalizationPolygon(tcx); +} + +void Sweep::SweepPoints(SweepContext& tcx) +{ + for (int i = 1; i < tcx.point_count(); i++) { + Point& point = *tcx.GetPoint(i); + Node* node = &PointEvent(tcx, point); + for (unsigned int i = 0; i < point.edge_list.size(); i++) { + EdgeEvent(tcx, point.edge_list[i], node); + } + } +} + +void Sweep::FinalizationPolygon(SweepContext& tcx) +{ + // Get an Internal triangle to start with + Triangle* t = tcx.front()->head()->next->triangle; + Point* p = tcx.front()->head()->next->point; + while (!t->GetConstrainedEdgeCW(*p)) { + t = t->NeighborCCW(*p); + } + + // Collect interior triangles constrained by edges + tcx.MeshClean(*t); +} + +Node& Sweep::PointEvent(SweepContext& tcx, Point& point) +{ + Node& node = tcx.LocateNode(point); + Node& new_node = NewFrontTriangle(tcx, point, node); + + // Only need to check +epsilon since point never have smaller + // x value than node due to how we fetch nodes from the front + if (point.x <= node.point->x + EPSILON) { + Fill(tcx, node); + } + + //tcx.AddNode(new_node); + + FillAdvancingFront(tcx, new_node); + return new_node; +} + +void Sweep::EdgeEvent(SweepContext& tcx, Edge* edge, Node* node) +{ + tcx.edge_event.constrained_edge = edge; + tcx.edge_event.right = (edge->p->x > edge->q->x); + + if (IsEdgeSideOfTriangle(*node->triangle, *edge->p, *edge->q)) { + return; + } + + // For now we will do all needed filling + // TODO: integrate with flip process might give some better performance + // but for now this avoid the issue with cases that needs both flips and fills + FillEdgeEvent(tcx, edge, node); + EdgeEvent(tcx, *edge->p, *edge->q, node->triangle, *edge->q); +} + +void Sweep::EdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* triangle, Point& point) +{ + if (IsEdgeSideOfTriangle(*triangle, ep, eq)) { + return; + } + + Point* p1 = triangle->PointCCW(point); + Orientation o1 = Orient2d(eq, *p1, ep); + if (o1 == COLLINEAR) { + if( triangle->Contains(&eq, p1)) { + triangle->MarkConstrainedEdge(&eq, p1 ); + // We are modifying the constraint maybe it would be better to + // not change the given constraint and just keep a variable for the new constraint + tcx.edge_event.constrained_edge->q = p1; + triangle = &triangle->NeighborAcross(point); + EdgeEvent( tcx, ep, *p1, triangle, *p1 ); + } else { + std::runtime_error("EdgeEvent - collinear points not supported"); + assert(0); + } + return; + } + + Point* p2 = triangle->PointCW(point); + Orientation o2 = Orient2d(eq, *p2, ep); + if (o2 == COLLINEAR) { + if( triangle->Contains(&eq, p2)) { + triangle->MarkConstrainedEdge(&eq, p2 ); + // We are modifying the constraint maybe it would be better to + // not change the given constraint and just keep a variable for the new constraint + tcx.edge_event.constrained_edge->q = p2; + triangle = &triangle->NeighborAcross(point); + EdgeEvent( tcx, ep, *p2, triangle, *p2 ); + } else { + std::runtime_error("EdgeEvent - collinear points not supported"); + assert(0); + } + return; + } + + if (o1 == o2) { + // Need to decide if we are rotating CW or CCW to get to a triangle + // that will cross edge + if (o1 == CW) { + triangle = triangle->NeighborCCW(point); + } else{ + triangle = triangle->NeighborCW(point); + } + EdgeEvent(tcx, ep, eq, triangle, point); + } else { + // This triangle crosses constraint so lets flippin start! + FlipEdgeEvent(tcx, ep, eq, triangle, point); + } +} + +bool Sweep::IsEdgeSideOfTriangle(Triangle& triangle, Point& ep, Point& eq) +{ + int index = triangle.EdgeIndex(&ep, &eq); + + if (index != -1) { + triangle.MarkConstrainedEdge(index); + Triangle* t = triangle.GetNeighbor(index); + if (t) { + t->MarkConstrainedEdge(&ep, &eq); + } + return true; + } + return false; +} + +Node& Sweep::NewFrontTriangle(SweepContext& tcx, Point& point, Node& node) +{ + Triangle* triangle = new Triangle(point, *node.point, *node.next->point); + + triangle->MarkNeighbor(*node.triangle); + tcx.AddToMap(triangle); + + Node* new_node = new Node(point); + nodes_.push_back(new_node); + + new_node->next = node.next; + new_node->prev = &node; + node.next->prev = new_node; + node.next = new_node; + + if (!Legalize(tcx, *triangle)) { + tcx.MapTriangleToNodes(*triangle); + } + + return *new_node; +} + +void Sweep::Fill(SweepContext& tcx, Node& node) +{ + Triangle* triangle = new Triangle(*node.prev->point, *node.point, *node.next->point); + + // TODO: should copy the constrained_edge value from neighbor triangles + // for now constrained_edge values are copied during the legalize + triangle->MarkNeighbor(*node.prev->triangle); + triangle->MarkNeighbor(*node.triangle); + + tcx.AddToMap(triangle); + + // Update the advancing front + node.prev->next = node.next; + node.next->prev = node.prev; + + // If it was legalized the triangle has already been mapped + if (!Legalize(tcx, *triangle)) { + tcx.MapTriangleToNodes(*triangle); + } + +} + +void Sweep::FillAdvancingFront(SweepContext& tcx, Node& n) +{ + + // Fill right holes + Node* node = n.next; + + while (node->next) { + // if HoleAngle exceeds 90 degrees then break. + if (LargeHole_DontFill(node)) break; + Fill(tcx, *node); + node = node->next; + } + + // Fill left holes + node = n.prev; + + while (node->prev) { + // if HoleAngle exceeds 90 degrees then break. + if (LargeHole_DontFill(node)) break; + Fill(tcx, *node); + node = node->prev; + } + + // Fill right basins + if (n.next && n.next->next) { + double angle = BasinAngle(n); + if (angle < PI_3div4) { + FillBasin(tcx, n); + } + } +} + +// True if HoleAngle exceeds 90 degrees. +bool Sweep::LargeHole_DontFill(Node* node) { + + Node* nextNode = node->next; + Node* prevNode = node->prev; + if (!AngleExceeds90Degrees(node->point, nextNode->point, prevNode->point)) + return false; + + // Check additional points on front. + Node* next2Node = nextNode->next; + // "..Plus.." because only want angles on same side as point being added. + if ((next2Node != NULL) && !AngleExceedsPlus90DegreesOrIsNegative(node->point, next2Node->point, prevNode->point)) + return false; + + Node* prev2Node = prevNode->prev; + // "..Plus.." because only want angles on same side as point being added. + if ((prev2Node != NULL) && !AngleExceedsPlus90DegreesOrIsNegative(node->point, nextNode->point, prev2Node->point)) + return false; + + return true; +} + +bool Sweep::AngleExceeds90Degrees(Point* origin, Point* pa, Point* pb) { + double angle = Angle(*origin, *pa, *pb); + bool exceeds90Degrees = ((angle > PI_div2) || (angle < -PI_div2)); + return exceeds90Degrees; +} + +bool Sweep::AngleExceedsPlus90DegreesOrIsNegative(Point* origin, Point* pa, Point* pb) { + double angle = Angle(*origin, *pa, *pb); + bool exceedsPlus90DegreesOrIsNegative = (angle > PI_div2) || (angle < 0); + return exceedsPlus90DegreesOrIsNegative; +} + +double Sweep::Angle(Point& origin, Point& pa, Point& pb) { + /* Complex plane + * ab = cosA +i*sinA + * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx) + * atan2(y,x) computes the principal value of the argument function + * applied to the complex number x+iy + * Where x = ax*bx + ay*by + * y = ax*by - ay*bx + */ + double px = origin.x; + double py = origin.y; + double ax = pa.x- px; + double ay = pa.y - py; + double bx = pb.x - px; + double by = pb.y - py; + double x = ax * by - ay * bx; + double y = ax * bx + ay * by; + double angle = atan2(x, y); + return angle; +} + +double Sweep::BasinAngle(Node& node) +{ + double ax = node.point->x - node.next->next->point->x; + double ay = node.point->y - node.next->next->point->y; + return atan2(ay, ax); +} + +double Sweep::HoleAngle(Node& node) +{ + /* Complex plane + * ab = cosA +i*sinA + * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx) + * atan2(y,x) computes the principal value of the argument function + * applied to the complex number x+iy + * Where x = ax*bx + ay*by + * y = ax*by - ay*bx + */ + double ax = node.next->point->x - node.point->x; + double ay = node.next->point->y - node.point->y; + double bx = node.prev->point->x - node.point->x; + double by = node.prev->point->y - node.point->y; + return atan2(ax * by - ay * bx, ax * bx + ay * by); +} + +bool Sweep::Legalize(SweepContext& tcx, Triangle& t) +{ + // To legalize a triangle we start by finding if any of the three edges + // violate the Delaunay condition + for (int i = 0; i < 3; i++) { + if (t.delaunay_edge[i]) + continue; + + Triangle* ot = t.GetNeighbor(i); + + if (ot) { + Point* p = t.GetPoint(i); + Point* op = ot->OppositePoint(t, *p); + int oi = ot->Index(op); + + // If this is a Constrained Edge or a Delaunay Edge(only during recursive legalization) + // then we should not try to legalize + if (ot->constrained_edge[oi] || ot->delaunay_edge[oi]) { + t.constrained_edge[i] = ot->constrained_edge[oi]; + continue; + } + + bool inside = Incircle(*p, *t.PointCCW(*p), *t.PointCW(*p), *op); + + if (inside) { + // Lets mark this shared edge as Delaunay + t.delaunay_edge[i] = true; + ot->delaunay_edge[oi] = true; + + // Lets rotate shared edge one vertex CW to legalize it + RotateTrianglePair(t, *p, *ot, *op); + + // We now got one valid Delaunay Edge shared by two triangles + // This gives us 4 new edges to check for Delaunay + + // Make sure that triangle to node mapping is done only one time for a specific triangle + bool not_legalized = !Legalize(tcx, t); + if (not_legalized) { + tcx.MapTriangleToNodes(t); + } + + not_legalized = !Legalize(tcx, *ot); + if (not_legalized) + tcx.MapTriangleToNodes(*ot); + + // Reset the Delaunay edges, since they only are valid Delaunay edges + // until we add a new triangle or point. + // XXX: need to think about this. Can these edges be tried after we + // return to previous recursive level? + t.delaunay_edge[i] = false; + ot->delaunay_edge[oi] = false; + + // If triangle have been legalized no need to check the other edges since + // the recursive legalization will handles those so we can end here. + return true; + } + } + } + return false; +} + +bool Sweep::Incircle(Point& pa, Point& pb, Point& pc, Point& pd) +{ + double adx = pa.x - pd.x; + double ady = pa.y - pd.y; + double bdx = pb.x - pd.x; + double bdy = pb.y - pd.y; + + double adxbdy = adx * bdy; + double bdxady = bdx * ady; + double oabd = adxbdy - bdxady; + + if (oabd <= 0) + return false; + + double cdx = pc.x - pd.x; + double cdy = pc.y - pd.y; + + double cdxady = cdx * ady; + double adxcdy = adx * cdy; + double ocad = cdxady - adxcdy; + + if (ocad <= 0) + return false; + + double bdxcdy = bdx * cdy; + double cdxbdy = cdx * bdy; + + double alift = adx * adx + ady * ady; + double blift = bdx * bdx + bdy * bdy; + double clift = cdx * cdx + cdy * cdy; + + double det = alift * (bdxcdy - cdxbdy) + blift * ocad + clift * oabd; + + return det > 0; +} + +void Sweep::RotateTrianglePair(Triangle& t, Point& p, Triangle& ot, Point& op) +{ + Triangle* n1, *n2, *n3, *n4; + n1 = t.NeighborCCW(p); + n2 = t.NeighborCW(p); + n3 = ot.NeighborCCW(op); + n4 = ot.NeighborCW(op); + + bool ce1, ce2, ce3, ce4; + ce1 = t.GetConstrainedEdgeCCW(p); + ce2 = t.GetConstrainedEdgeCW(p); + ce3 = ot.GetConstrainedEdgeCCW(op); + ce4 = ot.GetConstrainedEdgeCW(op); + + bool de1, de2, de3, de4; + de1 = t.GetDelunayEdgeCCW(p); + de2 = t.GetDelunayEdgeCW(p); + de3 = ot.GetDelunayEdgeCCW(op); + de4 = ot.GetDelunayEdgeCW(op); + + t.Legalize(p, op); + ot.Legalize(op, p); + + // Remap delaunay_edge + ot.SetDelunayEdgeCCW(p, de1); + t.SetDelunayEdgeCW(p, de2); + t.SetDelunayEdgeCCW(op, de3); + ot.SetDelunayEdgeCW(op, de4); + + // Remap constrained_edge + ot.SetConstrainedEdgeCCW(p, ce1); + t.SetConstrainedEdgeCW(p, ce2); + t.SetConstrainedEdgeCCW(op, ce3); + ot.SetConstrainedEdgeCW(op, ce4); + + // Remap neighbors + // XXX: might optimize the markNeighbor by keeping track of + // what side should be assigned to what neighbor after the + // rotation. Now mark neighbor does lots of testing to find + // the right side. + t.ClearNeighbors(); + ot.ClearNeighbors(); + if (n1) ot.MarkNeighbor(*n1); + if (n2) t.MarkNeighbor(*n2); + if (n3) t.MarkNeighbor(*n3); + if (n4) ot.MarkNeighbor(*n4); + t.MarkNeighbor(ot); +} + +void Sweep::FillBasin(SweepContext& tcx, Node& node) +{ + if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) { + tcx.basin.left_node = node.next->next; + } else { + tcx.basin.left_node = node.next; + } + + // Find the bottom and right node + tcx.basin.bottom_node = tcx.basin.left_node; + while (tcx.basin.bottom_node->next + && tcx.basin.bottom_node->point->y >= tcx.basin.bottom_node->next->point->y) { + tcx.basin.bottom_node = tcx.basin.bottom_node->next; + } + if (tcx.basin.bottom_node == tcx.basin.left_node) { + // No valid basin + return; + } + + tcx.basin.right_node = tcx.basin.bottom_node; + while (tcx.basin.right_node->next + && tcx.basin.right_node->point->y < tcx.basin.right_node->next->point->y) { + tcx.basin.right_node = tcx.basin.right_node->next; + } + if (tcx.basin.right_node == tcx.basin.bottom_node) { + // No valid basins + return; + } + + tcx.basin.width = tcx.basin.right_node->point->x - tcx.basin.left_node->point->x; + tcx.basin.left_highest = tcx.basin.left_node->point->y > tcx.basin.right_node->point->y; + + FillBasinReq(tcx, tcx.basin.bottom_node); +} + +void Sweep::FillBasinReq(SweepContext& tcx, Node* node) +{ + // if shallow stop filling + if (IsShallow(tcx, *node)) { + return; + } + + Fill(tcx, *node); + + if (node->prev == tcx.basin.left_node && node->next == tcx.basin.right_node) { + return; + } else if (node->prev == tcx.basin.left_node) { + Orientation o = Orient2d(*node->point, *node->next->point, *node->next->next->point); + if (o == CW) { + return; + } + node = node->next; + } else if (node->next == tcx.basin.right_node) { + Orientation o = Orient2d(*node->point, *node->prev->point, *node->prev->prev->point); + if (o == CCW) { + return; + } + node = node->prev; + } else { + // Continue with the neighbor node with lowest Y value + if (node->prev->point->y < node->next->point->y) { + node = node->prev; + } else { + node = node->next; + } + } + + FillBasinReq(tcx, node); +} + +bool Sweep::IsShallow(SweepContext& tcx, Node& node) +{ + double height; + + if (tcx.basin.left_highest) { + height = tcx.basin.left_node->point->y - node.point->y; + } else { + height = tcx.basin.right_node->point->y - node.point->y; + } + + // if shallow stop filling + if (tcx.basin.width > height) { + return true; + } + return false; +} + +void Sweep::FillEdgeEvent(SweepContext& tcx, Edge* edge, Node* node) +{ + if (tcx.edge_event.right) { + FillRightAboveEdgeEvent(tcx, edge, node); + } else { + FillLeftAboveEdgeEvent(tcx, edge, node); + } +} + +void Sweep::FillRightAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node) +{ + while (node->next->point->x < edge->p->x) { + // Check if next node is below the edge + if (Orient2d(*edge->q, *node->next->point, *edge->p) == CCW) { + FillRightBelowEdgeEvent(tcx, edge, *node); + } else { + node = node->next; + } + } +} + +void Sweep::FillRightBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) +{ + if (node.point->x < edge->p->x) { + if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) { + // Concave + FillRightConcaveEdgeEvent(tcx, edge, node); + } else{ + // Convex + FillRightConvexEdgeEvent(tcx, edge, node); + // Retry this one + FillRightBelowEdgeEvent(tcx, edge, node); + } + } +} + +void Sweep::FillRightConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) +{ + Fill(tcx, *node.next); + if (node.next->point != edge->p) { + // Next above or below edge? + if (Orient2d(*edge->q, *node.next->point, *edge->p) == CCW) { + // Below + if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) { + // Next is concave + FillRightConcaveEdgeEvent(tcx, edge, node); + } else { + // Next is convex + } + } + } + +} + +void Sweep::FillRightConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) +{ + // Next concave or convex? + if (Orient2d(*node.next->point, *node.next->next->point, *node.next->next->next->point) == CCW) { + // Concave + FillRightConcaveEdgeEvent(tcx, edge, *node.next); + } else{ + // Convex + // Next above or below edge? + if (Orient2d(*edge->q, *node.next->next->point, *edge->p) == CCW) { + // Below + FillRightConvexEdgeEvent(tcx, edge, *node.next); + } else{ + // Above + } + } +} + +void Sweep::FillLeftAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node) +{ + while (node->prev->point->x > edge->p->x) { + // Check if next node is below the edge + if (Orient2d(*edge->q, *node->prev->point, *edge->p) == CW) { + FillLeftBelowEdgeEvent(tcx, edge, *node); + } else { + node = node->prev; + } + } +} + +void Sweep::FillLeftBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) +{ + if (node.point->x > edge->p->x) { + if (Orient2d(*node.point, *node.prev->point, *node.prev->prev->point) == CW) { + // Concave + FillLeftConcaveEdgeEvent(tcx, edge, node); + } else { + // Convex + FillLeftConvexEdgeEvent(tcx, edge, node); + // Retry this one + FillLeftBelowEdgeEvent(tcx, edge, node); + } + } +} + +void Sweep::FillLeftConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) +{ + // Next concave or convex? + if (Orient2d(*node.prev->point, *node.prev->prev->point, *node.prev->prev->prev->point) == CW) { + // Concave + FillLeftConcaveEdgeEvent(tcx, edge, *node.prev); + } else{ + // Convex + // Next above or below edge? + if (Orient2d(*edge->q, *node.prev->prev->point, *edge->p) == CW) { + // Below + FillLeftConvexEdgeEvent(tcx, edge, *node.prev); + } else{ + // Above + } + } +} + +void Sweep::FillLeftConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node) +{ + Fill(tcx, *node.prev); + if (node.prev->point != edge->p) { + // Next above or below edge? + if (Orient2d(*edge->q, *node.prev->point, *edge->p) == CW) { + // Below + if (Orient2d(*node.point, *node.prev->point, *node.prev->prev->point) == CW) { + // Next is concave + FillLeftConcaveEdgeEvent(tcx, edge, node); + } else{ + // Next is convex + } + } + } + +} + +void Sweep::FlipEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* t, Point& p) +{ + Triangle& ot = t->NeighborAcross(p); + Point& op = *ot.OppositePoint(*t, p); + + if (&ot == NULL) { + // If we want to integrate the fillEdgeEvent do it here + // With current implementation we should never get here + //throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle"); + assert(0); + } + + if (InScanArea(p, *t->PointCCW(p), *t->PointCW(p), op)) { + // Lets rotate shared edge one vertex CW + RotateTrianglePair(*t, p, ot, op); + tcx.MapTriangleToNodes(*t); + tcx.MapTriangleToNodes(ot); + + if (p == eq && op == ep) { + if (eq == *tcx.edge_event.constrained_edge->q && ep == *tcx.edge_event.constrained_edge->p) { + t->MarkConstrainedEdge(&ep, &eq); + ot.MarkConstrainedEdge(&ep, &eq); + Legalize(tcx, *t); + Legalize(tcx, ot); + } else { + // XXX: I think one of the triangles should be legalized here? + } + } else { + Orientation o = Orient2d(eq, op, ep); + t = &NextFlipTriangle(tcx, (int)o, *t, ot, p, op); + FlipEdgeEvent(tcx, ep, eq, t, p); + } + } else { + Point& newP = NextFlipPoint(ep, eq, ot, op); + FlipScanEdgeEvent(tcx, ep, eq, *t, ot, newP); + EdgeEvent(tcx, ep, eq, t, p); + } +} + +Triangle& Sweep::NextFlipTriangle(SweepContext& tcx, int o, Triangle& t, Triangle& ot, Point& p, Point& op) +{ + if (o == CCW) { + // ot is not crossing edge after flip + int edge_index = ot.EdgeIndex(&p, &op); + ot.delaunay_edge[edge_index] = true; + Legalize(tcx, ot); + ot.ClearDelunayEdges(); + return t; + } + + // t is not crossing edge after flip + int edge_index = t.EdgeIndex(&p, &op); + + t.delaunay_edge[edge_index] = true; + Legalize(tcx, t); + t.ClearDelunayEdges(); + return ot; +} + +Point& Sweep::NextFlipPoint(Point& ep, Point& eq, Triangle& ot, Point& op) +{ + Orientation o2d = Orient2d(eq, op, ep); + if (o2d == CW) { + // Right + return *ot.PointCCW(op); + } else if (o2d == CCW) { + // Left + return *ot.PointCW(op); + } else{ + //throw new RuntimeException("[Unsupported] Opposing point on constrained edge"); + assert(0); + } +} + +void Sweep::FlipScanEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle& flip_triangle, + Triangle& t, Point& p) +{ + Triangle& ot = t.NeighborAcross(p); + Point& op = *ot.OppositePoint(t, p); + + if (&t.NeighborAcross(p) == NULL) { + // If we want to integrate the fillEdgeEvent do it here + // With current implementation we should never get here + //throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle"); + assert(0); + } + + if (InScanArea(eq, *flip_triangle.PointCCW(eq), *flip_triangle.PointCW(eq), op)) { + // flip with new edge op->eq + FlipEdgeEvent(tcx, eq, op, &ot, op); + // TODO: Actually I just figured out that it should be possible to + // improve this by getting the next ot and op before the the above + // flip and continue the flipScanEdgeEvent here + // set new ot and op here and loop back to inScanArea test + // also need to set a new flip_triangle first + // Turns out at first glance that this is somewhat complicated + // so it will have to wait. + } else{ + Point& newP = NextFlipPoint(ep, eq, ot, op); + FlipScanEdgeEvent(tcx, ep, eq, flip_triangle, ot, newP); + } +} + +Sweep::~Sweep() { + + // Clean up memory + for(int i = 0; i < nodes_.size(); i++) { + delete nodes_[i]; + } + +} + +} + diff --git a/polygon/poly2tri/sweep/sweep.h b/polygon/poly2tri/sweep/sweep.h new file mode 100644 index 0000000000..f62c4cc3f2 --- /dev/null +++ b/polygon/poly2tri/sweep/sweep.h @@ -0,0 +1,285 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ +/** + * Sweep-line, Constrained Delauney Triangulation (CDT) See: Domiter, V. and + * Zalik, B.(2008)'Sweep-line algorithm for constrained Delaunay triangulation', + * International Journal of Geographical Information Science + * + * "FlipScan" Constrained Edge Algorithm invented by Thomas Åhlén, thahlen@gmail.com + */ + +#ifndef SWEEP_H +#define SWEEP_H + +#include + +namespace p2t { + +class SweepContext; +struct Node; +struct Point; +struct Edge; +class Triangle; + +class Sweep +{ +public: + + /** + * Triangulate + * + * @param tcx + */ + void Triangulate(SweepContext& tcx); + + /** + * Destructor - clean up memory + */ + ~Sweep(); + +private: + + /** + * Start sweeping the Y-sorted point set from bottom to top + * + * @param tcx + */ + void SweepPoints(SweepContext& tcx); + + /** + * Find closes node to the left of the new point and + * create a new triangle. If needed new holes and basins + * will be filled to. + * + * @param tcx + * @param point + * @return + */ + Node& PointEvent(SweepContext& tcx, Point& point); + + /** + * + * + * @param tcx + * @param edge + * @param node + */ + void EdgeEvent(SweepContext& tcx, Edge* edge, Node* node); + + void EdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* triangle, Point& point); + + /** + * Creates a new front triangle and legalize it + * + * @param tcx + * @param point + * @param node + * @return + */ + Node& NewFrontTriangle(SweepContext& tcx, Point& point, Node& node); + + /** + * Adds a triangle to the advancing front to fill a hole. + * @param tcx + * @param node - middle node, that is the bottom of the hole + */ + void Fill(SweepContext& tcx, Node& node); + + /** + * Returns true if triangle was legalized + */ + bool Legalize(SweepContext& tcx, Triangle& t); + + /** + * Requirement:
+ * 1. a,b and c form a triangle.
+ * 2. a and d is know to be on opposite side of bc
+ * 
+   *                a
+   *                +
+   *               / \
+   *              /   \
+   *            b/     \c
+   *            +-------+
+   *           /    d    \
+   *          /           \
+   *
+ * Fact: d has to be in area B to have a chance to be inside the circle formed by + * a,b and c
+ * d is outside B if orient2d(a,b,d) or orient2d(c,a,d) is CW
+ * This preknowledge gives us a way to optimize the incircle test + * @param a - triangle point, opposite d + * @param b - triangle point + * @param c - triangle point + * @param d - point opposite a + * @return true if d is inside circle, false if on circle edge + */ + bool Incircle(Point& pa, Point& pb, Point& pc, Point& pd); + + /** + * Rotates a triangle pair one vertex CW + *
+   *       n2                    n2
+   *  P +-----+             P +-----+
+   *    | t  /|               |\  t |
+   *    |   / |               | \   |
+   *  n1|  /  |n3           n1|  \  |n3
+   *    | /   |    after CW   |   \ |
+   *    |/ oT |               | oT \|
+   *    +-----+ oP            +-----+
+   *       n4                    n4
+   *
+ */ + void RotateTrianglePair(Triangle& t, Point& p, Triangle& ot, Point& op); + + /** + * Fills holes in the Advancing Front + * + * + * @param tcx + * @param n + */ + void FillAdvancingFront(SweepContext& tcx, Node& n); + + // Decision-making about when to Fill hole. + // Contributed by ToolmakerSteve2 + bool LargeHole_DontFill(Node* node); + bool AngleExceeds90Degrees(Point* origin, Point* pa, Point* pb); + bool AngleExceedsPlus90DegreesOrIsNegative(Point* origin, Point* pa, Point* pb); + double Angle(Point& origin, Point& pa, Point& pb); + + /** + * + * @param node - middle node + * @return the angle between 3 front nodes + */ + double HoleAngle(Node& node); + + /** + * The basin angle is decided against the horizontal line [1,0] + */ + double BasinAngle(Node& node); + + /** + * Fills a basin that has formed on the Advancing Front to the right + * of given node.
+ * First we decide a left,bottom and right node that forms the + * boundaries of the basin. Then we do a reqursive fill. + * + * @param tcx + * @param node - starting node, this or next node will be left node + */ + void FillBasin(SweepContext& tcx, Node& node); + + /** + * Recursive algorithm to fill a Basin with triangles + * + * @param tcx + * @param node - bottom_node + * @param cnt - counter used to alternate on even and odd numbers + */ + void FillBasinReq(SweepContext& tcx, Node* node); + + bool IsShallow(SweepContext& tcx, Node& node); + + bool IsEdgeSideOfTriangle(Triangle& triangle, Point& ep, Point& eq); + + void FillEdgeEvent(SweepContext& tcx, Edge* edge, Node* node); + + void FillRightAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node); + + void FillRightBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node); + + void FillRightConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node); + + void FillRightConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node); + + void FillLeftAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node); + + void FillLeftBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node); + + void FillLeftConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node); + + void FillLeftConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node); + + void FlipEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* t, Point& p); + + /** + * After a flip we have two triangles and know that only one will still be + * intersecting the edge. So decide which to contiune with and legalize the other + * + * @param tcx + * @param o - should be the result of an orient2d( eq, op, ep ) + * @param t - triangle 1 + * @param ot - triangle 2 + * @param p - a point shared by both triangles + * @param op - another point shared by both triangles + * @return returns the triangle still intersecting the edge + */ + Triangle& NextFlipTriangle(SweepContext& tcx, int o, Triangle& t, Triangle& ot, Point& p, Point& op); + + /** + * When we need to traverse from one triangle to the next we need + * the point in current triangle that is the opposite point to the next + * triangle. + * + * @param ep + * @param eq + * @param ot + * @param op + * @return + */ + Point& NextFlipPoint(Point& ep, Point& eq, Triangle& ot, Point& op); + + /** + * Scan part of the FlipScan algorithm
+ * When a triangle pair isn't flippable we will scan for the next + * point that is inside the flip triangle scan area. When found + * we generate a new flipEdgeEvent + * + * @param tcx + * @param ep - last point on the edge we are traversing + * @param eq - first point on the edge we are traversing + * @param flipTriangle - the current triangle sharing the point eq with edge + * @param t + * @param p + */ + void FlipScanEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle& flip_triangle, Triangle& t, Point& p); + + void FinalizationPolygon(SweepContext& tcx); + + std::vector nodes_; + +}; + +} + +#endif diff --git a/polygon/poly2tri/sweep/sweep_context.cc b/polygon/poly2tri/sweep/sweep_context.cc new file mode 100644 index 0000000000..6c0b0447dd --- /dev/null +++ b/polygon/poly2tri/sweep/sweep_context.cc @@ -0,0 +1,216 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ +#include "sweep_context.h" +#include +#include "advancing_front.h" + +namespace p2t { + +SweepContext::SweepContext(std::vector polyline) : + front_(0), + head_(0), + tail_(0), + af_head_(0), + af_middle_(0), + af_tail_(0) +{ + basin = Basin(); + edge_event = EdgeEvent(); + + points_ = polyline; + + InitEdges(points_); +} + +void SweepContext::AddHole(std::vector polyline) +{ + InitEdges(polyline); + for(unsigned int i = 0; i < polyline.size(); i++) { + points_.push_back(polyline[i]); + } +} + +void SweepContext::AddPoint(Point* point) { + points_.push_back(point); +} + +std::vector SweepContext::GetTriangles() +{ + return triangles_; +} + +std::list SweepContext::GetMap() +{ + return map_; +} + +void SweepContext::InitTriangulation() +{ + double xmax(points_[0]->x), xmin(points_[0]->x); + double ymax(points_[0]->y), ymin(points_[0]->y); + + // Calculate bounds. + for (unsigned int i = 0; i < points_.size(); i++) { + Point& p = *points_[i]; + if (p.x > xmax) + xmax = p.x; + if (p.x < xmin) + xmin = p.x; + if (p.y > ymax) + ymax = p.y; + if (p.y < ymin) + ymin = p.y; + } + + double dx = kAlpha * (xmax - xmin); + double dy = kAlpha * (ymax - ymin); + head_ = new Point(xmax + dx, ymin - dy); + tail_ = new Point(xmin - dx, ymin - dy); + + // Sort points along y-axis + std::sort(points_.begin(), points_.end(), cmp); + +} + +void SweepContext::InitEdges(std::vector polyline) +{ + int num_points = polyline.size(); + for (int i = 0; i < num_points; i++) { + int j = i < num_points - 1 ? i + 1 : 0; + edge_list.push_back(new Edge(*polyline[i], *polyline[j])); + } +} + +Point* SweepContext::GetPoint(const int& index) +{ + return points_[index]; +} + +void SweepContext::AddToMap(Triangle* triangle) +{ + map_.push_back(triangle); +} + +Node& SweepContext::LocateNode(Point& point) +{ + // TODO implement search tree + return *front_->LocateNode(point.x); +} + +void SweepContext::CreateAdvancingFront(std::vector nodes) +{ + + (void) nodes; + // Initial triangle + Triangle* triangle = new Triangle(*points_[0], *tail_, *head_); + + map_.push_back(triangle); + + af_head_ = new Node(*triangle->GetPoint(1), *triangle); + af_middle_ = new Node(*triangle->GetPoint(0), *triangle); + af_tail_ = new Node(*triangle->GetPoint(2)); + front_ = new AdvancingFront(*af_head_, *af_tail_); + + // TODO: More intuitive if head is middles next and not previous? + // so swap head and tail + af_head_->next = af_middle_; + af_middle_->next = af_tail_; + af_middle_->prev = af_head_; + af_tail_->prev = af_middle_; +} + +void SweepContext::RemoveNode(Node* node) +{ + delete node; +} + +void SweepContext::MapTriangleToNodes(Triangle& t) +{ + for (int i = 0; i < 3; i++) { + if (!t.GetNeighbor(i)) { + Node* n = front_->LocatePoint(t.PointCW(*t.GetPoint(i))); + if (n) + n->triangle = &t; + } + } +} + +void SweepContext::RemoveFromMap(Triangle* triangle) +{ + map_.remove(triangle); +} + +void SweepContext::MeshClean(Triangle& triangle) +{ + std::vector triangles; + triangles.push_back(&triangle); + + while(!triangles.empty()){ + Triangle *t = triangles.back(); + triangles.pop_back(); + + if (t != NULL && !t->IsInterior()) { + t->IsInterior(true); + triangles_.push_back(t); + for (int i = 0; i < 3; i++) { + if (!t->constrained_edge[i]) + triangles.push_back(t->GetNeighbor(i)); + } + } + } +} + +SweepContext::~SweepContext() +{ + + // Clean up memory + + delete head_; + delete tail_; + delete front_; + delete af_head_; + delete af_middle_; + delete af_tail_; + + typedef std::list type_list; + + for(type_list::iterator iter = map_.begin(); iter != map_.end(); ++iter) { + Triangle* ptr = *iter; + delete ptr; + } + + for(unsigned int i = 0; i < edge_list.size(); i++) { + delete edge_list[i]; + } + +} + +} diff --git a/polygon/poly2tri/sweep/sweep_context.h b/polygon/poly2tri/sweep/sweep_context.h new file mode 100644 index 0000000000..1010c0e8a8 --- /dev/null +++ b/polygon/poly2tri/sweep/sweep_context.h @@ -0,0 +1,186 @@ +/* + * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors + * http://code.google.com/p/poly2tri/ + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of Poly2Tri nor the names of its contributors may be + * used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef SWEEP_CONTEXT_H +#define SWEEP_CONTEXT_H + +#include +#include +#include + +namespace p2t { + +// Inital triangle factor, seed triangle will extend 30% of +// PointSet width to both left and right. +const double kAlpha = 0.3; + +struct Point; +class Triangle; +struct Node; +struct Edge; +class AdvancingFront; + +class SweepContext { +public: + +/// Constructor +SweepContext(std::vector polyline); +/// Destructor +~SweepContext(); + +void set_head(Point* p1); + +Point* head(); + +void set_tail(Point* p1); + +Point* tail(); + +int point_count(); + +Node& LocateNode(Point& point); + +void RemoveNode(Node* node); + +void CreateAdvancingFront(std::vector nodes); + +/// Try to map a node to all sides of this triangle that don't have a neighbor +void MapTriangleToNodes(Triangle& t); + +void AddToMap(Triangle* triangle); + +Point* GetPoint(const int& index); + +Point* GetPoints(); + +void RemoveFromMap(Triangle* triangle); + +void AddHole(std::vector polyline); + +void AddPoint(Point* point); + +AdvancingFront* front(); + +void MeshClean(Triangle& triangle); + +std::vector GetTriangles(); +std::list GetMap(); + +std::vector edge_list; + +struct Basin { + Node* left_node; + Node* bottom_node; + Node* right_node; + double width; + bool left_highest; + + Basin() : left_node(NULL), bottom_node(NULL), right_node(NULL), width(0.0), left_highest(false) + { + } + + void Clear() + { + left_node = NULL; + bottom_node = NULL; + right_node = NULL; + width = 0.0; + left_highest = false; + } +}; + +struct EdgeEvent { + Edge* constrained_edge; + bool right; + + EdgeEvent() : constrained_edge(NULL), right(false) + { + } +}; + +Basin basin; +EdgeEvent edge_event; + +private: + +friend class Sweep; + +std::vector triangles_; +std::list map_; +std::vector points_; + +// Advancing front +AdvancingFront* front_; +// head point used with advancing front +Point* head_; +// tail point used with advancing front +Point* tail_; + +Node *af_head_, *af_middle_, *af_tail_; + +void InitTriangulation(); +void InitEdges(std::vector polyline); + +}; + +inline AdvancingFront* SweepContext::front() +{ + return front_; +} + +inline int SweepContext::point_count() +{ + return points_.size(); +} + +inline void SweepContext::set_head(Point* p1) +{ + head_ = p1; +} + +inline Point* SweepContext::head() +{ + return head_; +} + +inline void SweepContext::set_tail(Point* p1) +{ + tail_ = p1; +} + +inline Point* SweepContext::tail() +{ + return tail_; +} + +} + +#endif