From: duv Date: Fri, 5 Feb 2016 11:32:52 +0000 (+0300) Subject: Creaform restoreorient version for 6.9.0 X-Git-Url: http://git.dev.opencascade.org/gitweb/?a=commitdiff_plain;h=a83f359622d147df5f7706802885f17f37483d82;p=occt-copy.git Creaform restoreorient version for 6.9.0 --- diff --git a/src/BRepMesh/BRepMesh_Block.hxx b/src/BRepMesh/BRepMesh_Block.hxx new file mode 100644 index 0000000000..53c63caf20 --- /dev/null +++ b/src/BRepMesh/BRepMesh_Block.hxx @@ -0,0 +1,191 @@ +// Created on: 2015-11-27 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef __BLOCK_H__ +#define __BLOCK_H__ + +#include + +template class Block +{ +public: + + Block (int theSize, void (*theErrFunction) (char*) = NULL) + { + myFirst = myLast = NULL; + myBlockSize = theSize; + myErrorFun = theErrFunction; + } + + ~Block() + { + while (myFirst) + { + BlockStruct* aNext = myFirst->Next; + delete myFirst; + myFirst = aNext; + } + } + + Type* New (int theNum = 1) + { + Type* aLast; + + if (!myLast || myLast->Current + theNum > myLast->Last) + { + if (myLast && myLast->Next) { myLast = myLast->Next; } + else + { + BlockStruct* aNext = (BlockStruct*) new char [sizeof (BlockStruct) + (myBlockSize-1) *sizeof (Type)]; + + if (!aNext) { if (myErrorFun) { (*myErrorFun) ("Not enough memory!"); } exit (1); } + + if (myLast) { myLast->Next = aNext; } + else { myFirst = aNext; } + + myLast = aNext; + myLast->Current = & (myLast->Data[0]); + myLast->Last = myLast->Current + myBlockSize; + myLast->Next = NULL; + } + } + + aLast = myLast->Current; + myLast->Current += theNum; + return aLast; + } + + Type* ScanFirst() + { + myScanCurrentBlock = myFirst; + + if (!myScanCurrentBlock) { return NULL; } + + myScanCurrentData = & (myScanCurrentBlock->Data[0]); + return myScanCurrentData++; + } + + Type* ScanNext() + { + if (myScanCurrentData >= myScanCurrentBlock->Current) + { + myScanCurrentBlock = myScanCurrentBlock->Next; + + if (!myScanCurrentBlock) { return NULL; } + + myScanCurrentData = & (myScanCurrentBlock->Data[0]); + } + + return myScanCurrentData++; + } + + void Reset() + { + BlockStruct* aBlock; + + if (!myFirst) { return; } + + for (aBlock = myFirst; ; aBlock = aBlock->Next) + { + aBlock->Current = & (aBlock->Data[0]); + + if (aBlock == myLast) { break; } + } + + myLast = myFirst; + } + +private: + + typedef struct BlockSt + { + Type* Current; + Type* Last; + struct BlockSt* Next; + Type Data[1]; + } BlockStruct; + + int myBlockSize; + BlockStruct* myFirst; + BlockStruct* myLast; + + BlockStruct* myScanCurrentBlock; + Type* myScanCurrentData; + + void (*myErrorFun) (char*); +}; + +template class DBlock +{ +public: + + DBlock (int theSize, void (*theErrFunction) (char*) = NULL) { myFirst = NULL; myFirstFree = NULL; myBLockSize = theSize; myErrFun = theErrFunction; } + + ~DBlock() { while (myFirst) { BlockStruct* aNext = myFirst->Next; delete myFirst; myFirst = aNext; } } + + Type* New() + { + BlockItem* anItem; + + if (!myFirstFree) + { + BlockStruct* next = myFirst; + myFirst = (BlockStruct*) new char [sizeof (BlockStruct) + (myBLockSize-1) *sizeof (BlockItem)]; + + if (!myFirst) { if (myErrFun) { (*myErrFun) ("Not enough memory!"); } exit (1); } + + myFirstFree = & (myFirst->Data[0]); + + for (anItem=myFirstFree; anItemNextFree = anItem + 1; } + + anItem->NextFree = NULL; + myFirst->Next = next; + } + + anItem = myFirstFree; + myFirstFree = anItem->NextFree; + return (Type*) anItem; + } + + void Delete (Type* t) + { + ( (BlockItem*) t)->NextFree = myFirstFree; + myFirstFree = (BlockItem*) t; + } + +private: + + typedef union BlockItemSt + { + Type Item; + BlockItemSt* NextFree; + } BlockItem; + + typedef struct BlockSt + { + struct BlockSt* Next; + BlockItem Data[1]; + } BlockStruct; + + int myBLockSize; + BlockStruct* myFirst; + BlockItem* myFirstFree; + + void (*myErrFun) (char*); +}; + +#endif + diff --git a/src/BRepMesh/BRepMesh_EdgeSet.hxx b/src/BRepMesh/BRepMesh_EdgeSet.hxx new file mode 100644 index 0000000000..5927203508 --- /dev/null +++ b/src/BRepMesh/BRepMesh_EdgeSet.hxx @@ -0,0 +1,98 @@ +// Created on: 2015-08-28 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _BRepMesh_EdgeSet_Header +#define _BRepMesh_EdgeSet_Header + +#include +#include + +//! Class representing set of boundary edges. +template +class BRepMesh_EdgeSet : public BVH_PrimitiveSet +{ +public: + + typedef typename BVH::VectorType::Type BVH_VecNt; + typedef typename BVH::ArrayType::Type BVH_ArrayNt; + +public: + + //! Creates empty edge set. + BRepMesh_EdgeSet (const BVH_ArrayNt& theTriangulationVertices); + + //! Releases resources of edge set. + virtual ~BRepMesh_EdgeSet(); + +private: + + BRepMesh_EdgeSet (const BRepMesh_EdgeSet& theOther); + + BRepMesh_EdgeSet& operator= (const BRepMesh_EdgeSet& theOther); + +public: + + //! Array of indices of edge vertices. + BVH_Array4i Elements; + +private: + + const BVH_ArrayNt& Vertices; + +public: + + //! Returns index of edge which is closest to the given one. + Standard_Integer ClosestEdge (const BVH_VecNt& theEdgeV1, + const BVH_VecNt& theEdgeV2, + T* theDist = NULL); + + //! Returns index of edge which is closest to the given one (accounting for edge orientation). + Standard_Integer CoherentEdge (const BVH_VecNt& theEdgeV1, + const BVH_VecNt& theEdgeV2, + T* theDist = NULL); + + //! Returns the distance between to edges. + static T EdgeToEdgeDistance (const BVH_VecNt& theE1V1, + const BVH_VecNt& theE1V2, + const BVH_VecNt& theE2V1, + const BVH_VecNt& theE2V2); + + //! Returns the distance between to edges. + static T EdgeToBoxDistance (const BVH_VecNt& theEV1, + const BVH_VecNt& theEV2, + const BVH_VecNt& theMinPoint, + const BVH_VecNt& theMaxPoint); + + //! Returns total number of edges. + virtual Standard_Integer Size() const; + + //! Returns AABB of entire set of objects. + using BVH_PrimitiveSet::Box; + + //! Returns AABB of the given edge. + virtual BVH_Box Box (const Standard_Integer theIndex) const; + + //! Returns centroid position along the given axis. + virtual T Center (const Standard_Integer theIndex, + const Standard_Integer theAxis) const; + + //! Performs transposing the two given edges in the set. + virtual void Swap (const Standard_Integer theIndex1, + const Standard_Integer theIndex2); + +}; + +#include +#endif // _BRepMesh_EdgeSet_Header diff --git a/src/BRepMesh/BRepMesh_EdgeSet.lxx b/src/BRepMesh/BRepMesh_EdgeSet.lxx new file mode 100644 index 0000000000..ab5324eb9d --- /dev/null +++ b/src/BRepMesh/BRepMesh_EdgeSet.lxx @@ -0,0 +1,452 @@ +// Created on: 2015-08-28 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +// ======================================================================= +// function : BRepMesh_EdgeSet +// purpose : +// ======================================================================= +template +BRepMesh_EdgeSet::BRepMesh_EdgeSet (const BVH_ArrayNt& theTriangulationVertices) + : Vertices (theTriangulationVertices) +{ + // +} + +// ======================================================================= +// function : ~BRepMesh_EdgeSet +// purpose : +// ======================================================================= +template +BRepMesh_EdgeSet::~BRepMesh_EdgeSet() +{ + // +} + +// ======================================================================= +// function : Size +// purpose : +// ======================================================================= +template +Standard_Integer BRepMesh_EdgeSet::Size() const +{ + return BVH::Array::Size (Elements); +} + +// ======================================================================= +// function : Box +// purpose : +// ======================================================================= +template +BVH_Box BRepMesh_EdgeSet::Box (const Standard_Integer theIndex) const +{ + const BVH_Vec4i& anIndex = BVH::Array::Value (Elements, theIndex); + + const BVH_VecNt& aPoint0 = BVH::Array::Value (Vertices, anIndex.x()); + const BVH_VecNt& aPoint1 = BVH::Array::Value (Vertices, anIndex.y()); + + BVH_VecNt aMinPoint = aPoint0; + BVH_VecNt aMaxPoint = aPoint0; + + BVH::BoxMinMax::CwiseMin (aMinPoint, aPoint1); + BVH::BoxMinMax::CwiseMax (aMaxPoint, aPoint1); + + return BVH_Box (aMinPoint, aMaxPoint); +} + +// ======================================================================= +// function : Center +// purpose : +// ======================================================================= +template +T BRepMesh_EdgeSet::Center (const Standard_Integer theIndex, + const Standard_Integer theAxis) const +{ + const BVH_Vec4i& anIndex = BVH::Array::Value (Elements, theIndex); + + const BVH_VecNt& aPoint0 = BVH::Array::Value (Vertices, anIndex.x()); + const BVH_VecNt& aPoint1 = BVH::Array::Value (Vertices, anIndex.y()); + + return ( BVH::VecComp::Get (aPoint0, theAxis) + + BVH::VecComp::Get (aPoint1, theAxis) ) * static_cast (1.0 / 2.0); +} + +// ======================================================================= +// function : Swap +// purpose : +// ======================================================================= +template +void BRepMesh_EdgeSet::Swap (const Standard_Integer theIndex1, + const Standard_Integer theIndex2) +{ + BVH_Vec4i& anIndices1 = BVH::Array::ChangeValue (Elements, theIndex1); + BVH_Vec4i& anIndices2 = BVH::Array::ChangeValue (Elements, theIndex2); + + std::swap (anIndices1, anIndices2); +} + +// ======================================================================= +// function : ClosestEdge +// purpose : +// ======================================================================= +template +Standard_Integer BRepMesh_EdgeSet::ClosestEdge (const BVH_VecNt& theEdgeV1, + const BVH_VecNt& theEdgeV2, + T* theDist) +{ + Standard_Integer aStack[32]; + + const NCollection_Handle >& aBVH = BVH(); + + if (aBVH.IsNull()) + { + return -1; + } + + Standard_Integer aNode = 0; + + T aMinDist = std::numeric_limits::max(); + Standard_Integer aClosestEdge = 0; + + Standard_Integer aHead = -1; // stack head position + + for (;;) + { + if (aBVH->IsOuter (aNode)) + { + for (Standard_Integer anElemNb = aBVH->BegPrimitive (aNode); anElemNb <= aBVH->EndPrimitive (aNode); ++anElemNb) + { + const BVH_Vec4i& anEdge = Elements[anElemNb]; + const BVH_VecNt& aVrt0 = Vertices[anEdge.x()]; + const BVH_VecNt& aVrt1 = Vertices[anEdge.y()]; + T aDist = EdgeToEdgeDistance (theEdgeV1, theEdgeV2, aVrt0, aVrt1); + + if (aDist < aMinDist) + { + aMinDist = aDist; + aClosestEdge = anElemNb; + } + } + + if (aHead < 0) + { + break; + } + + aNode = aStack[aHead--]; + } + else + { + T aLftDist = EdgeToBoxDistance (theEdgeV1, + theEdgeV2, + aBVH->MinPoint (aBVH->LeftChild (aNode)), + aBVH->MaxPoint (aBVH->LeftChild (aNode))); + + T aRghDist = EdgeToBoxDistance (theEdgeV1, + theEdgeV2, + aBVH->MinPoint (aBVH->RightChild (aNode)), + aBVH->MaxPoint (aBVH->RightChild (aNode))); + + const bool aProceedLft = aLftDist < aMinDist; + const bool aProceedRgh = aRghDist < aMinDist; + + if (aProceedLft && aProceedRgh) + { + int aLeft = aBVH->LeftChild (aNode); + int aRight = aBVH->RightChild (aNode); + + aNode = (aLftDist < aRghDist) ? aLeft : aRight; + aStack[++aHead] = (aRghDist <= aLftDist) ? aLeft : aRight; + } + else + { + if (aProceedLft || aProceedRgh) + { + aNode = aProceedLft ? aBVH->LeftChild (aNode) : aBVH->RightChild (aNode); + } + else + { + if (aHead < 0) + { + break; + } + + aNode = aStack[aHead--]; + } + } + } + } + + if (theDist != NULL) + { + *theDist = aMinDist; + } + + return aClosestEdge; +} + +// ======================================================================= +// function : CoherentEdge +// purpose : +// ======================================================================= +template +Standard_Integer BRepMesh_EdgeSet::CoherentEdge (const BVH_VecNt& theEdgeV1, + const BVH_VecNt& theEdgeV2, + T* theDist) +{ + Standard_Integer aStack[32]; + + const NCollection_Handle >& aBVH = BVH(); + + if (aBVH.IsNull()) + { + return -1; + } + + Standard_Integer aNode = 0; + + T aMinDist = std::numeric_limits::max(); + Standard_Integer aClosestEdge = 0; + + Standard_Integer aHead = -1; // stack head position + + for (;;) + { + if (aBVH->IsOuter (aNode)) + { + for (Standard_Integer anElemNb = aBVH->BegPrimitive (aNode); anElemNb <= aBVH->EndPrimitive (aNode); ++anElemNb) + { + const BVH_Vec4i& anEdge = Elements[anElemNb]; + const BVH_VecNt& aVrt0 = Vertices[anEdge.x()]; + const BVH_VecNt& aVrt1 = Vertices[anEdge.y()]; + T aDist = EdgeToEdgeDistance (theEdgeV1, theEdgeV2, aVrt0, aVrt1); + + // add scaled vertex distances + aDist += (1e-6 * (theEdgeV1 - aVrt0).Modulus()); + aDist += (1e-6 * (theEdgeV1 - aVrt1).Modulus()); + aDist += (1e-6 * (theEdgeV2 - aVrt0).Modulus()); + aDist += (1e-6 * (theEdgeV2 - aVrt1).Modulus()); + + if (aDist < aMinDist) + { + aMinDist = aDist; + aClosestEdge = anElemNb; + } + } + + if (aHead < 0) + { + break; + } + + aNode = aStack[aHead--]; + } + else + { + T aLftDist = EdgeToBoxDistance (theEdgeV1, + theEdgeV2, + aBVH->MinPoint (aBVH->LeftChild (aNode)), + aBVH->MaxPoint (aBVH->LeftChild (aNode))); + + T aRghDist = EdgeToBoxDistance (theEdgeV1, + theEdgeV2, + aBVH->MinPoint (aBVH->RightChild (aNode)), + aBVH->MaxPoint (aBVH->RightChild (aNode))); + + const bool aProceedLft = aLftDist < aMinDist; + const bool aProceedRgh = aRghDist < aMinDist; + + if (aProceedLft && aProceedRgh) + { + int aLeft = aBVH->LeftChild (aNode); + int aRight = aBVH->RightChild (aNode); + + aNode = (aLftDist < aRghDist) ? aLeft : aRight; + aStack[++aHead] = (aRghDist <= aLftDist) ? aLeft : aRight; + } + else + { + if (aProceedLft || aProceedRgh) + { + aNode = aProceedLft ? aBVH->LeftChild (aNode) : aBVH->RightChild (aNode); + } + else + { + if (aHead < 0) + { + break; + } + + aNode = aStack[aHead--]; + } + } + } + } + + if (theDist != NULL) + { + *theDist = aMinDist; + } + + return aClosestEdge; +} + +#define PARALLEL_THRESOLD 1e-8 + +// ======================================================================= +// function : EdgeToEdgeDistance +// purpose : +// ======================================================================= +template +T BRepMesh_EdgeSet::EdgeToEdgeDistance (const BVH_VecNt& theE1V1, + const BVH_VecNt& theE1V2, + const BVH_VecNt& theE2V1, + const BVH_VecNt& theE2V2) +{ + BVH_VecNt u = theE1V2 - theE1V1; + BVH_VecNt v = theE2V2 - theE2V1; + BVH_VecNt w = theE1V1 - theE2V1; + Standard_Real a = u.Dot (u); // always >= 0 + Standard_Real b = u.Dot (v); + Standard_Real c = v.Dot (v); // always >= 0 + Standard_Real d = u.Dot (w); + Standard_Real e = v.Dot (w); + Standard_Real D = a * c - b * b; // always >= 0 + Standard_Real sc, sN, sD = D; // sc = sN / sD, default sD = D >= 0 + Standard_Real tc, tN, tD = D; // tc = tN / tD, default tD = D >= 0 + + // compute the line parameters of the two closest points + if (D < PARALLEL_THRESOLD) // the lines are almost parallel + { + sN = 0.0; // force using point P0 on segment S1 + sD = 1.0; // to prevent possible division by 0.0 later + tN = e; + tD = c; + } + else // get the closest points on the infinite lines + { + sN = (b * e - c * d); + tN = (a * e - b * d); + + if (sN < 0.0) // sc < 0 => the s=0 edge is visible + { + sN = 0.0; + tN = e; + tD = c; + } + else if (sN > sD) // sc > 1 => the s=1 edge is visible + { + sN = sD; + tN = e + b; + tD = c; + } + } + + if (tN < 0.0) // tc < 0 => the t=0 edge is visible + { + tN = 0.0; + + // recompute sc for this edge + if (-d < 0.0) + { + sN = 0.0; + } + else if (-d > a) + { + sN = sD; + } + else + { + sN = -d; + sD = a; + } + } + else if (tN > tD) // tc > 1 => the t=1 edge is visible + { + tN = tD; + + // recompute sc for this edge + if ((-d + b) < 0.0) + { + sN = 0; + } + else if ((-d + b) > a) + { + sN = sD; + } + else + { + sN = (-d + b); + sD = a; + } + } + + // finally do the division to get sc and tc + sc = (abs (sN) < PARALLEL_THRESOLD ? 0.0 : sN / sD); + tc = (abs (tN) < PARALLEL_THRESOLD ? 0.0 : tN / tD); + + // get the difference of the two closest points + BVH_VecNt dP = w + (u * sc) - (v * tc); // = S1(sc) - S2(tc) + + return dP.Modulus(); // return the closest distance +} + +#undef PARALLEL_THRESOLD + +// ======================================================================= +// function : EdgeToBoxDistance +// purpose : +// ======================================================================= +template +T BRepMesh_EdgeSet::EdgeToBoxDistance (const BVH_VecNt& theEV1, + const BVH_VecNt& theEV2, + const BVH_VecNt& theMinPoint, + const BVH_VecNt& theMaxPoint) +{ + T anSqrDist = std::numeric_limits::max(); + + BVH_VecNt aVClosest = theEV1; + aVClosest = aVClosest.cwiseMax (theMinPoint); + aVClosest = aVClosest.cwiseMin (theMaxPoint); + + anSqrDist = (theEV1 - aVClosest).SquareModulus(); + + if (anSqrDist == (T) 0.0) + { + return anSqrDist; + } + + aVClosest = theEV2; + aVClosest = aVClosest.cwiseMax (theMinPoint); + aVClosest = aVClosest.cwiseMin (theMaxPoint); + + anSqrDist = std::min (anSqrDist, (theEV2 - aVClosest).SquareModulus()); + + if (anSqrDist == (T) 0.0) + { + return anSqrDist; + } + + LineToBoxDistance aDist; + LineToBoxDistance::Result aResult = aDist (theEV1, theEV2 - theEV1, theMinPoint, theMaxPoint); + + if (aResult.lineParameter > (T) 0.0 && aResult.lineParameter < (T) 1.0) + { + anSqrDist = std::min (anSqrDist, aResult.sqrDistance); + } + + return sqrt (anSqrDist); +} diff --git a/src/BRepMesh/BRepMesh_LineBoxDistance.hxx b/src/BRepMesh/BRepMesh_LineBoxDistance.hxx new file mode 100644 index 0000000000..a0a70d0388 --- /dev/null +++ b/src/BRepMesh/BRepMesh_LineBoxDistance.hxx @@ -0,0 +1,605 @@ +// Created on: 2015-11-27 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _BRepMesh_LineBoxDistance_Header +#define _BRepMesh_LineBoxDistance_Header + +#include + +template +class LineToBoxDistance +{ +public: + + typedef typename BVH::VectorType::Type BVH_VecNt; + +public: + struct Result + { + T distance, sqrDistance; + T lineParameter; + BVH_VecNt closestPoint[2]; + }; + + Result operator() (const BVH_VecNt& theOrigin, + const BVH_VecNt& theDirect, + const BVH_VecNt& theMinPoint, + const BVH_VecNt& theMaxPoint); + +protected: + // Compute the distance and closest point between a line and an + // axis-aligned box whose center is the origin. On input, 'point' is the + // line origin and 'direction' is the line direction. On output, 'point' + // is the point on the box closest to the line. The 'direction' is + // non-const to allow transforming the problem into the first octant. + void DoQuery (BVH_VecNt& point, BVH_VecNt& direction, + BVH_VecNt const& boxExtent, Result& result); + +private: + void Face (int i0, int i1, int i2, BVH_VecNt& pnt, + BVH_VecNt const& dir, BVH_VecNt const& PmE, + BVH_VecNt const& boxExtent, Result& result); + + void CaseNoZeros (BVH_VecNt& pnt, BVH_VecNt const& dir, + BVH_VecNt const& boxExtent, Result& result); + + void Case0 (int i0, int i1, int i2, BVH_VecNt& pnt, + BVH_VecNt const& dir, BVH_VecNt const& boxExtent, + Result& result); + + void Case00 (int i0, int i1, int i2, BVH_VecNt& pnt, + BVH_VecNt const& dir, BVH_VecNt const& boxExtent, + Result& result); + + void Case000 (BVH_VecNt& pnt, BVH_VecNt const& boxExtent, + Result& result); +}; + +// ======================================================================= +// function : +// purpose : +// ======================================================================= +template +typename LineToBoxDistance::Result +LineToBoxDistance::operator() (const BVH_VecNt& theOrigin, + const BVH_VecNt& theDirect, + const BVH_VecNt& theMinPoint, + const BVH_VecNt& theMaxPoint) +{ + // Translate the line and box so that the box has center at the origin. + BVH_VecNt boxCenter, boxExtent; + boxCenter = (theMaxPoint + theMinPoint) * (T) 0.5; + boxExtent = (theMaxPoint - theMinPoint) * (T) 0.5; + BVH_VecNt point = theOrigin - boxCenter; + BVH_VecNt direction = theDirect; + + Result result; + DoQuery (point, direction, boxExtent, result); + + // Compute the closest point on the line. + result.closestPoint[0] = + theOrigin + theDirect * result.lineParameter; + + // Compute the closest point on the box. + result.closestPoint[1] = boxCenter + point; + return result; +} + +// ======================================================================= +// function : DoQuery +// purpose : +// ======================================================================= +template +void LineToBoxDistance::DoQuery ( + BVH_VecNt& point, BVH_VecNt& direction, + BVH_VecNt const& boxExtent, Result& result) +{ + result.sqrDistance = (T)0; + result.lineParameter = (T)0; + + // Apply reflections so that direction vector has nonnegative components. + bool reflect[3]; + + for (int i = 0; i < 3; ++i) + { + if (direction[i] < (T)0) + { + point[i] = -point[i]; + direction[i] = -direction[i]; + reflect[i] = true; + } + else + { + reflect[i] = false; + } + } + + if (direction[0] > (T)0) + { + if (direction[1] > (T)0) + { + if (direction[2] > (T)0) // (+,+,+) + { + CaseNoZeros (point, direction, boxExtent, result); + } + else // (+,+,0) + { + Case0 (0, 1, 2, point, direction, boxExtent, result); + } + } + else + { + if (direction[2] > (T)0) // (+,0,+) + { + Case0 (0, 2, 1, point, direction, boxExtent, result); + } + else // (+,0,0) + { + Case00 (0, 1, 2, point, direction, boxExtent, result); + } + } + } + else + { + if (direction[1] > (T)0) + { + if (direction[2] > (T)0) // (0,+,+) + { + Case0 (1, 2, 0, point, direction, boxExtent, result); + } + else // (0,+,0) + { + Case00 (1, 0, 2, point, direction, boxExtent, result); + } + } + else + { + if (direction[2] > (T)0) // (0,0,+) + { + Case00 (2, 0, 1, point, direction, boxExtent, result); + } + else // (0,0,0) + { + Case000 (point, boxExtent, result); + } + } + } + + // Undo the reflections applied previously. + for (int i = 0; i < 3; ++i) + { + if (reflect[i]) + { + point[i] = -point[i]; + } + } + + result.distance = sqrt (result.sqrDistance); +} + +// ======================================================================= +// function : Face +// purpose : +// ======================================================================= +template +void LineToBoxDistance::Face (int i0, int i1, + int i2, BVH_VecNt& pnt, BVH_VecNt const& dir, + BVH_VecNt const& PmE, BVH_VecNt const& boxExtent, Result& result) +{ + BVH_VecNt PpE; + T lenSqr, inv, tmp, param, t, delta; + + PpE[i1] = pnt[i1] + boxExtent[i1]; + PpE[i2] = pnt[i2] + boxExtent[i2]; + + if (dir[i0] * PpE[i1] >= dir[i1] * PmE[i0]) + { + if (dir[i0] * PpE[i2] >= dir[i2] * PmE[i0]) + { + // v[i1] >= -e[i1], v[i2] >= -e[i2] (distance = 0) + pnt[i0] = boxExtent[i0]; + inv = ((T)1) / dir[i0]; + pnt[i1] -= dir[i1] * PmE[i0] * inv; + pnt[i2] -= dir[i2] * PmE[i0] * inv; + result.lineParameter = -PmE[i0] * inv; + } + else + { + // v[i1] >= -e[i1], v[i2] < -e[i2] + lenSqr = dir[i0] * dir[i0] + dir[i2] * dir[i2]; + tmp = lenSqr * PpE[i1] - dir[i1] * (dir[i0] * PmE[i0] + + dir[i2] * PpE[i2]); + + if (tmp <= ((T)2)*lenSqr * boxExtent[i1]) + { + t = tmp / lenSqr; + lenSqr += dir[i1] * dir[i1]; + tmp = PpE[i1] - t; + delta = dir[i0] * PmE[i0] + dir[i1] * tmp + dir[i2] * PpE[i2]; + param = -delta / lenSqr; + result.sqrDistance += PmE[i0] * PmE[i0] + tmp * tmp + + PpE[i2] * PpE[i2] + delta * param; + + result.lineParameter = param; + pnt[i0] = boxExtent[i0]; + pnt[i1] = t - boxExtent[i1]; + pnt[i2] = -boxExtent[i2]; + } + else + { + lenSqr += dir[i1] * dir[i1]; + delta = dir[i0] * PmE[i0] + dir[i1] * PmE[i1] + dir[i2] * PpE[i2]; + param = -delta / lenSqr; + result.sqrDistance += PmE[i0] * PmE[i0] + PmE[i1] * PmE[i1] + + PpE[i2] * PpE[i2] + delta * param; + + result.lineParameter = param; + pnt[i0] = boxExtent[i0]; + pnt[i1] = boxExtent[i1]; + pnt[i2] = -boxExtent[i2]; + } + } + } + else + { + if (dir[i0] * PpE[i2] >= dir[i2] * PmE[i0]) + { + // v[i1] < -e[i1], v[i2] >= -e[i2] + lenSqr = dir[i0] * dir[i0] + dir[i1] * dir[i1]; + tmp = lenSqr * PpE[i2] - dir[i2] * (dir[i0] * PmE[i0] + + dir[i1] * PpE[i1]); + + if (tmp <= ((T)2)*lenSqr * boxExtent[i2]) + { + t = tmp / lenSqr; + lenSqr += dir[i2] * dir[i2]; + tmp = PpE[i2] - t; + delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * tmp; + param = -delta / lenSqr; + result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] + + tmp * tmp + delta * param; + + result.lineParameter = param; + pnt[i0] = boxExtent[i0]; + pnt[i1] = -boxExtent[i1]; + pnt[i2] = t - boxExtent[i2]; + } + else + { + lenSqr += dir[i2] * dir[i2]; + delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * PmE[i2]; + param = -delta / lenSqr; + result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] + + PmE[i2] * PmE[i2] + delta * param; + + result.lineParameter = param; + pnt[i0] = boxExtent[i0]; + pnt[i1] = -boxExtent[i1]; + pnt[i2] = boxExtent[i2]; + } + } + else + { + // v[i1] < -e[i1], v[i2] < -e[i2] + lenSqr = dir[i0] * dir[i0] + dir[i2] * dir[i2]; + tmp = lenSqr * PpE[i1] - dir[i1] * (dir[i0] * PmE[i0] + + dir[i2] * PpE[i2]); + + if (tmp >= (T)0) + { + // v[i1]-edge is closest + if (tmp <= ((T)2)*lenSqr * boxExtent[i1]) + { + t = tmp / lenSqr; + lenSqr += dir[i1] * dir[i1]; + tmp = PpE[i1] - t; + delta = dir[i0] * PmE[i0] + dir[i1] * tmp + dir[i2] * PpE[i2]; + param = -delta / lenSqr; + result.sqrDistance += PmE[i0] * PmE[i0] + tmp * tmp + + PpE[i2] * PpE[i2] + delta * param; + + result.lineParameter = param; + pnt[i0] = boxExtent[i0]; + pnt[i1] = t - boxExtent[i1]; + pnt[i2] = -boxExtent[i2]; + } + else + { + lenSqr += dir[i1] * dir[i1]; + delta = dir[i0] * PmE[i0] + dir[i1] * PmE[i1] + + dir[i2] * PpE[i2]; + param = -delta / lenSqr; + result.sqrDistance += PmE[i0] * PmE[i0] + PmE[i1] * PmE[i1] + + PpE[i2] * PpE[i2] + delta * param; + + result.lineParameter = param; + pnt[i0] = boxExtent[i0]; + pnt[i1] = boxExtent[i1]; + pnt[i2] = -boxExtent[i2]; + } + + return; + } + + lenSqr = dir[i0] * dir[i0] + dir[i1] * dir[i1]; + tmp = lenSqr * PpE[i2] - dir[i2] * (dir[i0] * PmE[i0] + + dir[i1] * PpE[i1]); + + if (tmp >= (T)0) + { + // v[i2]-edge is closest + if (tmp <= ((T)2)*lenSqr * boxExtent[i2]) + { + t = tmp / lenSqr; + lenSqr += dir[i2] * dir[i2]; + tmp = PpE[i2] - t; + delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * tmp; + param = -delta / lenSqr; + result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] + + tmp * tmp + delta * param; + + result.lineParameter = param; + pnt[i0] = boxExtent[i0]; + pnt[i1] = -boxExtent[i1]; + pnt[i2] = t - boxExtent[i2]; + } + else + { + lenSqr += dir[i2] * dir[i2]; + delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + + dir[i2] * PmE[i2]; + param = -delta / lenSqr; + result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] + + PmE[i2] * PmE[i2] + delta * param; + + result.lineParameter = param; + pnt[i0] = boxExtent[i0]; + pnt[i1] = -boxExtent[i1]; + pnt[i2] = boxExtent[i2]; + } + + return; + } + + // (v[i1],v[i2])-corner is closest + lenSqr += dir[i2] * dir[i2]; + delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * PpE[i2]; + param = -delta / lenSqr; + result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] + + PpE[i2] * PpE[i2] + delta * param; + + result.lineParameter = param; + pnt[i0] = boxExtent[i0]; + pnt[i1] = -boxExtent[i1]; + pnt[i2] = -boxExtent[i2]; + } + } +} + +// ======================================================================= +// function : CaseNoZeros +// purpose : +// ======================================================================= +template +void LineToBoxDistance::CaseNoZeros ( + BVH_VecNt& pnt, BVH_VecNt const& dir, + BVH_VecNt const& boxExtent, Result& result) +{ + BVH_VecNt PmE = pnt - boxExtent; + T prodDxPy = dir[0] * PmE[1]; + T prodDyPx = dir[1] * PmE[0]; + T prodDzPx, prodDxPz, prodDzPy, prodDyPz; + + if (prodDyPx >= prodDxPy) + { + prodDzPx = dir[2] * PmE[0]; + prodDxPz = dir[0] * PmE[2]; + + if (prodDzPx >= prodDxPz) + { + // line intersects x = e0 + Face (0, 1, 2, pnt, dir, PmE, boxExtent, result); + } + else + { + // line intersects z = e2 + Face (2, 0, 1, pnt, dir, PmE, boxExtent, result); + } + } + else + { + prodDzPy = dir[2] * PmE[1]; + prodDyPz = dir[1] * PmE[2]; + + if (prodDzPy >= prodDyPz) + { + // line intersects y = e1 + Face (1, 2, 0, pnt, dir, PmE, boxExtent, result); + } + else + { + // line intersects z = e2 + Face (2, 0, 1, pnt, dir, PmE, boxExtent, result); + } + } +} + +// ======================================================================= +// function : Case0 +// purpose : +// ======================================================================= +template +void LineToBoxDistance::Case0 (int i0, int i1, + int i2, BVH_VecNt& pnt, BVH_VecNt const& dir, + BVH_VecNt const& boxExtent, Result& result) +{ + T PmE0 = pnt[i0] - boxExtent[i0]; + T PmE1 = pnt[i1] - boxExtent[i1]; + T prod0 = dir[i1] * PmE0; + T prod1 = dir[i0] * PmE1; + T delta, invLSqr, inv; + + if (prod0 >= prod1) + { + // line intersects P[i0] = e[i0] + pnt[i0] = boxExtent[i0]; + + T PpE1 = pnt[i1] + boxExtent[i1]; + delta = prod0 - dir[i0] * PpE1; + + if (delta >= (T)0) + { + invLSqr = ((T)1) / (dir[i0] * dir[i0] + dir[i1] * dir[i1]); + result.sqrDistance += delta * delta * invLSqr; + pnt[i1] = -boxExtent[i1]; + result.lineParameter = - (dir[i0] * PmE0 + dir[i1] * PpE1) * invLSqr; + } + else + { + inv = ((T)1) / dir[i0]; + pnt[i1] -= prod0 * inv; + result.lineParameter = -PmE0 * inv; + } + } + else + { + // line intersects P[i1] = e[i1] + pnt[i1] = boxExtent[i1]; + + T PpE0 = pnt[i0] + boxExtent[i0]; + delta = prod1 - dir[i1] * PpE0; + + if (delta >= (T)0) + { + invLSqr = ((T)1) / (dir[i0] * dir[i0] + dir[i1] * dir[i1]); + result.sqrDistance += delta * delta * invLSqr; + pnt[i0] = -boxExtent[i0]; + result.lineParameter = - (dir[i0] * PpE0 + dir[i1] * PmE1) * invLSqr; + } + else + { + inv = ((T)1) / dir[i1]; + pnt[i0] -= prod1 * inv; + result.lineParameter = -PmE1 * inv; + } + } + + if (pnt[i2] < -boxExtent[i2]) + { + delta = pnt[i2] + boxExtent[i2]; + result.sqrDistance += delta * delta; + pnt[i2] = -boxExtent[i2]; + } + else if (pnt[i2] > boxExtent[i2]) + { + delta = pnt[i2] - boxExtent[i2]; + result.sqrDistance += delta * delta; + pnt[i2] = boxExtent[i2]; + } +} + +// ======================================================================= +// function : Case00 +// purpose : +// ======================================================================= +template +void LineToBoxDistance::Case00 (int i0, int i1, + int i2, BVH_VecNt& pnt, BVH_VecNt const& dir, + BVH_VecNt const& boxExtent, Result& result) +{ + T delta; + + result.lineParameter = (boxExtent[i0] - pnt[i0]) / dir[i0]; + + pnt[i0] = boxExtent[i0]; + + if (pnt[i1] < -boxExtent[i1]) + { + delta = pnt[i1] + boxExtent[i1]; + result.sqrDistance += delta * delta; + pnt[i1] = -boxExtent[i1]; + } + else if (pnt[i1] > boxExtent[i1]) + { + delta = pnt[i1] - boxExtent[i1]; + result.sqrDistance += delta * delta; + pnt[i1] = boxExtent[i1]; + } + + if (pnt[i2] < -boxExtent[i2]) + { + delta = pnt[i2] + boxExtent[i2]; + result.sqrDistance += delta * delta; + pnt[i2] = -boxExtent[i2]; + } + else if (pnt[i2] > boxExtent[i2]) + { + delta = pnt[i2] - boxExtent[i2]; + result.sqrDistance += delta * delta; + pnt[i2] = boxExtent[i2]; + } +} + +// ======================================================================= +// function : Case000 +// purpose : +// ======================================================================= +template +void LineToBoxDistance::Case000 ( + BVH_VecNt& pnt, BVH_VecNt const& boxExtent, Result& result) +{ + T delta; + + if (pnt[0] < -boxExtent[0]) + { + delta = pnt[0] + boxExtent[0]; + result.sqrDistance += delta * delta; + pnt[0] = -boxExtent[0]; + } + else if (pnt[0] > boxExtent[0]) + { + delta = pnt[0] - boxExtent[0]; + result.sqrDistance += delta * delta; + pnt[0] = boxExtent[0]; + } + + if (pnt[1] < -boxExtent[1]) + { + delta = pnt[1] + boxExtent[1]; + result.sqrDistance += delta * delta; + pnt[1] = -boxExtent[1]; + } + else if (pnt[1] > boxExtent[1]) + { + delta = pnt[1] - boxExtent[1]; + result.sqrDistance += delta * delta; + pnt[1] = boxExtent[1]; + } + + if (pnt[2] < -boxExtent[2]) + { + delta = pnt[2] + boxExtent[2]; + result.sqrDistance += delta * delta; + pnt[2] = -boxExtent[2]; + } + else if (pnt[2] > boxExtent[2]) + { + delta = pnt[2] - boxExtent[2]; + result.sqrDistance += delta * delta; + pnt[2] = boxExtent[2]; + } +} + +#endif // _BRepMesh_LineBoxDistance_Header diff --git a/src/BRepMesh/BRepMesh_MinStCut.cxx b/src/BRepMesh/BRepMesh_MinStCut.cxx new file mode 100644 index 0000000000..b2ef61f0af --- /dev/null +++ b/src/BRepMesh/BRepMesh_MinStCut.cxx @@ -0,0 +1,1254 @@ +// Created on: 2015-11-27 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#pragma warning (disable: 4706) +#pragma warning (disable: 4701) +#pragma warning (disable: 4127) + +#include +#include + +Graph::Graph (void (*theErrFun) (char*)) +{ + myErrorFun = theErrFun; + myNodeBlockFirst = NULL; + myArcForBlockFirst = NULL; + myArcREvBlockFirst = NULL; + myFlow = 0; +} + +Graph::~Graph() +{ + while (myNodeBlockFirst) + { + NodeBlock* aNext = myNodeBlockFirst->Next; + delete myNodeBlockFirst; + myNodeBlockFirst = aNext; + } + + while (myArcForBlockFirst) + { + ArcForBlock* next = myArcForBlockFirst->Next; + delete myArcForBlockFirst->Start; + myArcForBlockFirst = next; + } + + while (myArcREvBlockFirst) + { + ArcRevBlock* next = myArcREvBlockFirst->Next; + delete myArcREvBlockFirst->Start; + myArcREvBlockFirst = next; + } +} + +Graph::NodeId Graph::AddNode() +{ + Node* aNode; + + if (!myNodeBlockFirst || myNodeBlockFirst->Current+1 > &myNodeBlockFirst->Nodes[NODE_BLOCK_SIZE-1]) + { + NodeBlock* next = myNodeBlockFirst; + myNodeBlockFirst = (NodeBlock*) new NodeBlock; + + if (!myNodeBlockFirst) { if (myErrorFun) { (*myErrorFun) ("Not enough memory!"); } exit (1); } + + myNodeBlockFirst->Current = & (myNodeBlockFirst->Nodes[0]); + myNodeBlockFirst->Next = next; + } + + aNode = myNodeBlockFirst->Current ++; + aNode->FirstOut = (ArcForward*) 0; + aNode->FirstIn = (ArcReverse*) 0; + + aNode->TrCapacity = 0; + + return (NodeId) aNode; +} + +void Graph::AddEdge (NodeId theFromNode, NodeId theToNode, CapacityType theCapacity, CapacityType theReverseCapacity) +{ + ArcForward* anArcFor; + ArcReverse* anArcRev; + + if (!myArcForBlockFirst || myArcForBlockFirst->Current+1 > &myArcForBlockFirst->ArcsFor[ARC_BLOCK_SIZE]) + { + ArcForBlock* aNext = myArcForBlockFirst; + char* aPtr = new char[sizeof (ArcForBlock) +1]; + + if (!aPtr) { if (myErrorFun) { (*myErrorFun) ("Not enough memory!"); } exit (1); } + + if ( (int) aPtr & 1) { myArcForBlockFirst = (ArcForBlock*) (aPtr + 1); } + else { myArcForBlockFirst = (ArcForBlock*) aPtr; } + + myArcForBlockFirst->Start = aPtr; + myArcForBlockFirst->Current = & (myArcForBlockFirst->ArcsFor[0]); + myArcForBlockFirst->Next = aNext; + } + + if (!myArcREvBlockFirst || myArcREvBlockFirst->Current+1 > &myArcREvBlockFirst->ArcsRev[ARC_BLOCK_SIZE]) + { + ArcRevBlock* aNext = myArcREvBlockFirst; + char* aPter = new char[sizeof (ArcRevBlock) +1]; + + if (!aPter) { if (myErrorFun) { (*myErrorFun) ("Not enough memory!"); } exit (1); } + + if ( (int) aPter & 1) { myArcREvBlockFirst = (ArcRevBlock*) (aPter + 1); } + else { myArcREvBlockFirst = (ArcRevBlock*) aPter; } + + myArcREvBlockFirst->Start = aPter; + myArcREvBlockFirst->Current = & (myArcREvBlockFirst->ArcsRev[0]); + myArcREvBlockFirst->Next = aNext; + } + + anArcFor = myArcForBlockFirst->Current ++; + anArcRev = myArcREvBlockFirst->Current ++; + + anArcRev->Sister = (ArcForward*) theFromNode; + anArcFor->Shift = (int) theToNode; + anArcFor->ResidualCap = theCapacity; + anArcFor->ReverseResidualCap = theReverseCapacity; + + ( (Node*) theFromNode)->FirstOut = + (ArcForward*) ( (int) ( ( (Node*) theFromNode)->FirstOut) + 1); + ( (Node*) theToNode)->FirstIn = + (ArcReverse*) ( (int) ( ( (Node*) theToNode)->FirstIn) + 1); +} + +void Graph::SetTWeights (NodeId theNode, CapacityType theCapacityToSource, CapacityType theCapacityToSink) +{ + myFlow += (theCapacityToSource < theCapacityToSink) ? theCapacityToSource : theCapacityToSink; + ( (Node*) theNode)->TrCapacity = theCapacityToSource - theCapacityToSink; +} + +void Graph::AddTWeights (NodeId theNode, CapacityType theCapacityToSource, CapacityType theCapacityToSink) +{ + register CapacityType delta = ( (Node*) theNode)->TrCapacity; + + if (delta > 0) { theCapacityToSource += delta; } + else { theCapacityToSink -= delta; } + + myFlow += (theCapacityToSource < theCapacityToSink) ? theCapacityToSource : theCapacityToSink; + ( (Node*) theNode)->TrCapacity = theCapacityToSource - theCapacityToSink; +} + +void Graph::prepareGraph() +{ + Node* aNode; + ArcForBlock* anArcBlockFor, *anArcBlockForFirst; + ArcRevBlock* anArcBlockRev, *anArcBlockRevFirst, *anArcBlockREvScan; + ArcForward* anArcFor; + ArcReverse* anArcRev, *anArcRevSCan, anArcREvTmp; + NodeBlock* aNodeBlock; + bool aForwardFlag = false, aReverseFlag = false; + int aK; + + if (!myArcREvBlockFirst) + { + NodeId aNodeFrom = AddNode(), aNodeTo = AddNode(); + AddEdge (aNodeFrom, aNodeTo, 1, 0); + } + + /* FIRST STAGE */ + anArcREvTmp.Sister = NULL; + + for (anArcRev=myArcREvBlockFirst->Current; anArcRev<&myArcREvBlockFirst->ArcsRev[ARC_BLOCK_SIZE]; anArcRev++) + { + anArcRev->Sister = NULL; + } + + anArcBlockFor = anArcBlockForFirst = myArcForBlockFirst; + anArcBlockRev = anArcBlockRevFirst = anArcBlockREvScan = myArcREvBlockFirst; + anArcFor = &anArcBlockFor->ArcsFor[0]; + anArcRev = anArcRevSCan = &anArcBlockRev->ArcsRev[0]; + + for (aNodeBlock=myNodeBlockFirst; aNodeBlock; aNodeBlock=aNodeBlock->Next) + { + for (aNode=&aNodeBlock->Nodes[0]; aNodeCurrent; aNode++) + { + /* outgoing arcs */ + aK = (int) aNode->FirstOut; + + if (anArcFor + aK > &anArcBlockFor->ArcsFor[ARC_BLOCK_SIZE]) + { + if (aK > ARC_BLOCK_SIZE) { if (myErrorFun) { (*myErrorFun) ("# of arcs per node exceeds block size!"); } exit (1); } + + if (aForwardFlag) { anArcBlockFor = NULL; } + else { anArcBlockFor = anArcBlockFor->Next; anArcBlockREvScan = anArcBlockREvScan->Next; } + + if (anArcBlockFor == NULL) + { + ArcForBlock* next = myArcForBlockFirst; + char* ptr = new char[sizeof (ArcForBlock) +1]; + + if (!ptr) { if (myErrorFun) { (*myErrorFun) ("Not enough memory!"); } exit (1); } + + if ( (int) ptr & 1) { myArcForBlockFirst = (ArcForBlock*) (ptr + 1); } + else { myArcForBlockFirst = (ArcForBlock*) ptr; } + + myArcForBlockFirst->Start = ptr; + myArcForBlockFirst->Current = & (myArcForBlockFirst->ArcsFor[0]); + myArcForBlockFirst->Next = next; + anArcBlockFor = myArcForBlockFirst; + aForwardFlag = true; + } + else { anArcRevSCan = &anArcBlockREvScan->ArcsRev[0]; } + + anArcFor = &anArcBlockFor->ArcsFor[0]; + } + + if (anArcBlockREvScan) + { + anArcRevSCan += aK; + aNode->Parent = (ArcForward*) anArcRevSCan; + } + else { aNode->Parent = (ArcForward*) &anArcREvTmp; } + + anArcFor += aK; + aNode->FirstOut = anArcFor; + anArcBlockFor->last_node = aNode; + + /* incoming arcs */ + aK = (int) aNode->FirstIn; + + if (anArcRev + aK > &anArcBlockRev->ArcsRev[ARC_BLOCK_SIZE]) + { + if (aK > ARC_BLOCK_SIZE) { if (myErrorFun) { (*myErrorFun) ("# of arcs per node exceeds block size!"); } exit (1); } + + if (aReverseFlag) { anArcBlockRev = NULL; } + else { anArcBlockRev = anArcBlockRev->Next; } + + if (anArcBlockRev == NULL) + { + ArcRevBlock* next = myArcREvBlockFirst; + char* ptr = new char[sizeof (ArcRevBlock) +1]; + + if (!ptr) { if (myErrorFun) { (*myErrorFun) ("Not enough memory!"); } exit (1); } + + if ( (int) ptr & 1) { myArcREvBlockFirst = (ArcRevBlock*) (ptr + 1); } + else { myArcREvBlockFirst = (ArcRevBlock*) ptr; } + + myArcREvBlockFirst->Start = ptr; + myArcREvBlockFirst->Current = & (myArcREvBlockFirst->ArcsRev[0]); + myArcREvBlockFirst->Next = next; + anArcBlockRev = myArcREvBlockFirst; + aReverseFlag = true; + } + + anArcRev = &anArcBlockRev->ArcsRev[0]; + } + + anArcRev += aK; + aNode->FirstIn = anArcRev; + anArcBlockRev->last_node = aNode; + } + + aNode->FirstOut = anArcFor; + aNode->FirstIn = anArcRev; + } + + for (anArcBlockFor=myArcForBlockFirst; anArcBlockFor; anArcBlockFor=anArcBlockFor->Next) + { + anArcBlockFor->Current = anArcBlockFor->last_node->FirstOut; + } + + for (anArcBlockFor=anArcBlockForFirst, anArcBlockRev=anArcBlockRevFirst; + anArcBlockFor; + anArcBlockFor=anArcBlockFor->Next, anArcBlockRev=anArcBlockRev->Next) + for (anArcFor=&anArcBlockFor->ArcsFor[0], anArcRev=&anArcBlockRev->ArcsRev[0]; + anArcFor<&anArcBlockFor->ArcsFor[ARC_BLOCK_SIZE]; + anArcFor++, anArcRev++) + { + ArcForward* anArcForward; + ArcReverse* anArcReverse; + Node* aNodeFrom; + int aShift = 0, aShiftNew; + CapacityType aRCap, aReverseRCap, aRCapNew, aReverseRCapNew; + + if (! (aNodeFrom= (Node*) (anArcRev->Sister))) { continue; } + + anArcForward = anArcFor; + anArcReverse = anArcRev; + + do + { + anArcReverse->Sister = NULL; + + aShiftNew = (int)(((char*) (anArcForward->Shift)) - (char*) aNodeFrom); + aRCapNew = anArcForward->ResidualCap; + aReverseRCapNew = anArcForward->ReverseResidualCap; + + if (aShift) + { + anArcForward->Shift = aShift; + anArcForward->ResidualCap = aRCap; + anArcForward->ReverseResidualCap = aReverseRCap; + } + + aShift = aShiftNew; + aRCap = aRCapNew; + aReverseRCap = aReverseRCapNew; + + anArcForward = -- aNodeFrom->FirstOut; + + if ( (ArcReverse*) (aNodeFrom->Parent) != &anArcREvTmp) + { + aNodeFrom->Parent = (ArcForward*) ( ( (ArcReverse*) (aNodeFrom->Parent)) - 1); + anArcReverse = (ArcReverse*) (aNodeFrom->Parent); + } + } + while (aNodeFrom= (Node*) (anArcReverse->Sister)); + + anArcForward->Shift = aShift; + anArcForward->ResidualCap = aRCap; + anArcForward->ReverseResidualCap = aReverseRCap; + } + + for (anArcBlockFor=myArcForBlockFirst; anArcBlockFor; anArcBlockFor=anArcBlockFor->Next) + { + aNode = anArcBlockFor->last_node; + anArcFor = aNode->FirstOut; + anArcBlockFor->Current->Shift = anArcFor->Shift; + anArcBlockFor->Current->ResidualCap = anArcFor->ResidualCap; + anArcBlockFor->Current->ReverseResidualCap = anArcFor->ReverseResidualCap; + anArcFor->Shift = (int) (anArcBlockFor->Current + 1); + aNode->FirstOut = (ArcForward*) ( ( (char*) anArcFor) - 1); + } + + /* THIRD STAGE */ + for (anArcBlockRev=myArcREvBlockFirst; anArcBlockRev; anArcBlockRev=anArcBlockRev->Next) + { + anArcBlockRev->Current = anArcBlockRev->last_node->FirstIn; + } + + for (aNodeBlock = myNodeBlockFirst; aNodeBlock; aNodeBlock=aNodeBlock->Next) + for (aNode = &aNodeBlock->Nodes[0]; aNodeCurrent; aNode++) + { + ArcForward* aForwardFirst, *aForwardLAst; + + aForwardFirst = aNode->FirstOut; + + if (IS_ODD (aForwardFirst)) + { + aForwardFirst = (ArcForward*) ( ( (char*) aForwardFirst) + 1); + aForwardLAst = (ArcForward*) ( (aForwardFirst ++)->Shift); + } + else { aForwardLAst = (aNode + 1)->FirstOut; } + + for (anArcFor=aForwardFirst; anArcForShift); + anArcRev = -- to->FirstIn; + anArcRev->Sister = anArcFor; + } + } + + for (anArcBlockRev=myArcREvBlockFirst; anArcBlockRev; anArcBlockRev=anArcBlockRev->Next) + { + aNode = anArcBlockRev->last_node; + anArcRev = aNode->FirstIn; + anArcBlockRev->Current->Sister = anArcRev->Sister; + anArcRev->Sister = (ArcForward*) (anArcBlockRev->Current + 1); + aNode->FirstIn = (ArcReverse*) ( ( (char*) anArcRev) - 1); + } +} + +#define TERMINAL ( (ArcForward *) 1 ) +#define ORPHAN ( (ArcForward *) 2 ) + +#define INFINITE_D 1000000000 + +inline void Graph::setActive (Node* theNode) +{ + if (!theNode->Next) + { + if (myQueueLast[1]) { myQueueLast[1] -> Next = theNode; } + else { myQueueFirst[1] = theNode; } + + myQueueLast[1] = theNode; + theNode -> Next = theNode; + } +} + +inline Graph::Node* Graph::nextActive() +{ + Node* aNode; + + while (true) + { + if (! (aNode=myQueueFirst[0])) + { + myQueueFirst[0] = aNode = myQueueFirst[1]; + myQueueLast[0] = myQueueLast[1]; + myQueueFirst[1] = NULL; + myQueueLast[1] = NULL; + + if (!aNode) { return NULL; } + } + + if (aNode->Next == aNode) { myQueueFirst[0] = myQueueLast[0] = NULL; } + else { myQueueFirst[0] = aNode -> Next; } + + aNode -> Next = NULL; + + if (aNode->Parent) { return aNode; } + } +} + +void Graph::maxflowInit() +{ + Node* aNode; + NodeBlock* aNodeBlock; + + myQueueFirst[0] = myQueueLast[0] = NULL; + myQueueFirst[1] = myQueueLast[1] = NULL; + myOrphanFirst = NULL; + + for (aNodeBlock=myNodeBlockFirst; aNodeBlock; aNodeBlock=aNodeBlock->Next) + { + for (aNode=&aNodeBlock->Nodes[0]; aNodeCurrent; aNode++) + { + aNode -> Next = NULL; + aNode -> TimeStamp = 0; + + if (aNode->TrCapacity > 0) + { + aNode -> IsSink = 0; + aNode -> Parent = TERMINAL; + setActive (aNode); + aNode -> TimeStamp = 0; + aNode -> Distance = 1; + } + else if (aNode->TrCapacity < 0) + { + aNode -> IsSink = 1; + aNode -> Parent = TERMINAL; + setActive (aNode); + aNode -> TimeStamp = 0; + aNode -> Distance = 1; + } + else + { + aNode -> Parent = NULL; + } + } + } + + myTime = 0; +} + +void Graph::augment (Node* theSStart, Node* theTStart, CapacityType* theCapMiddle, CapacityType* theRevCapMiddle) +{ + Node* aNode; + ArcForward* anArc; + CapacityType aBottleneck; + NodePtr* aNodePtr; + + aBottleneck = *theCapMiddle; + + for (aNode=theSStart; ;) + { + anArc = aNode -> Parent; + + if (anArc == TERMINAL) { break; } + + if (IS_ODD (anArc)) + { + anArc = MAKE_EVEN (anArc); + + if (aBottleneck > anArc->ResidualCap) { aBottleneck = anArc -> ResidualCap; } + + aNode = NEIGHBOR_NODE_REV (aNode, anArc -> Shift); + } + else + { + if (aBottleneck > anArc->ReverseResidualCap) { aBottleneck = anArc -> ReverseResidualCap; } + + aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); + } + } + + if (aBottleneck > aNode->TrCapacity) { aBottleneck = aNode -> TrCapacity; } + + for (aNode=theTStart; ;) + { + anArc = aNode -> Parent; + + if (anArc == TERMINAL) { break; } + + if (IS_ODD (anArc)) + { + anArc = MAKE_EVEN (anArc); + + if (aBottleneck > anArc->ReverseResidualCap) { aBottleneck = anArc -> ReverseResidualCap; } + + aNode = NEIGHBOR_NODE_REV (aNode, anArc -> Shift); + } + else + { + if (aBottleneck > anArc->ResidualCap) { aBottleneck = anArc -> ResidualCap; } + + aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); + } + } + + if (aBottleneck > - aNode->TrCapacity) { aBottleneck = - aNode -> TrCapacity; } + + + *theRevCapMiddle += aBottleneck; + *theCapMiddle -= aBottleneck; + + for (aNode=theSStart; ;) + { + anArc = aNode -> Parent; + + if (anArc == TERMINAL) { break; } + + if (IS_ODD (anArc)) + { + anArc = MAKE_EVEN (anArc); + anArc -> ReverseResidualCap += aBottleneck; + anArc -> ResidualCap -= aBottleneck; + + if (!anArc->ResidualCap) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + aNodePtr -> Next = myOrphanFirst; + myOrphanFirst = aNodePtr; + } + + aNode = NEIGHBOR_NODE_REV (aNode, anArc -> Shift); + } + else + { + anArc -> ResidualCap += aBottleneck; + anArc -> ReverseResidualCap -= aBottleneck; + + if (!anArc->ReverseResidualCap) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + aNodePtr -> Next = myOrphanFirst; + myOrphanFirst = aNodePtr; + } + + aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); + } + } + + aNode -> TrCapacity -= aBottleneck; + + if (!aNode->TrCapacity) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + aNodePtr -> Next = myOrphanFirst; + myOrphanFirst = aNodePtr; + } + + for (aNode=theTStart; ;) + { + anArc = aNode -> Parent; + + if (anArc == TERMINAL) { break; } + + if (IS_ODD (anArc)) + { + anArc = MAKE_EVEN (anArc); + anArc -> ResidualCap += aBottleneck; + anArc -> ReverseResidualCap -= aBottleneck; + + if (!anArc->ReverseResidualCap) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + aNodePtr -> Next = myOrphanFirst; + myOrphanFirst = aNodePtr; + } + + aNode = NEIGHBOR_NODE_REV (aNode, anArc -> Shift); + } + else + { + anArc -> ReverseResidualCap += aBottleneck; + anArc -> ResidualCap -= aBottleneck; + + if (!anArc->ResidualCap) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + aNodePtr -> Next = myOrphanFirst; + myOrphanFirst = aNodePtr; + } + + aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); + } + } + + aNode -> TrCapacity += aBottleneck; + + if (!aNode->TrCapacity) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + aNodePtr -> Next = myOrphanFirst; + myOrphanFirst = aNodePtr; + } + + myFlow += aBottleneck; +} + +void Graph::processSourceOrphan (Node* theNode) +{ + Node* aNode; + ArcForward* anArc0For, *anArc0ForFirst, *anArc0ForLast; + ArcReverse* anArc0Rev, *anArc0RevFirst, *anArc0RevLast; + ArcForward* anArc0Min = NULL, *anArc; + NodePtr* aNodePtr; + int aDist, aDistMin = INFINITE_D; + + anArc0ForFirst = theNode -> FirstOut; + + if (IS_ODD (anArc0ForFirst)) + { + anArc0ForFirst = (ArcForward*) ( ( (char*) anArc0ForFirst) + 1); + anArc0ForLast = (ArcForward*) ( (anArc0ForFirst ++) -> Shift); + } + else { anArc0ForLast = (theNode + 1) -> FirstOut; } + + anArc0RevFirst = theNode -> FirstIn; + + if (IS_ODD (anArc0RevFirst)) + { + anArc0RevFirst = (ArcReverse*) ( ( (char*) anArc0RevFirst) + 1); + anArc0RevLast = (ArcReverse*) ( (anArc0RevFirst ++) -> Sister); + } + else { anArc0RevLast = (theNode + 1) -> FirstIn; } + + + for (anArc0For=anArc0ForFirst; anArc0ForReverseResidualCap) + { + aNode = NEIGHBOR_NODE (theNode, anArc0For -> Shift); + + if (!aNode->IsSink && (anArc=aNode->Parent)) + { + /* checking the origin of j */ + aDist = 0; + + while (true) + { + if (aNode->TimeStamp == myTime) + { + aDist += aNode -> Distance; + break; + } + + anArc = aNode -> Parent; + aDist ++; + + if (anArc==TERMINAL) + { + aNode -> TimeStamp = myTime; + aNode -> Distance = 1; + break; + } + + if (anArc==ORPHAN) { aDist = INFINITE_D; break; } + + if (IS_ODD (anArc)) + { aNode = NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc) -> Shift); } + else + { aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); } + } + + if (aDistShift); aNode->TimeStamp!=myTime;) + { + aNode -> TimeStamp = myTime; + aNode -> Distance = aDist --; + anArc = aNode->Parent; + + if (IS_ODD (anArc)) + { aNode = NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc) -> Shift); } + else + { aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); } + } + } + } + } + + for (anArc0Rev=anArc0RevFirst; anArc0Rev Sister; + + if (anArc0For->ResidualCap) + { + aNode = NEIGHBOR_NODE_REV (theNode, anArc0For -> Shift); + + if (!aNode->IsSink && (anArc=aNode->Parent)) + { + aDist = 0; + + while (true) + { + if (aNode->TimeStamp == myTime) + { + aDist += aNode -> Distance; + break; + } + + anArc = aNode -> Parent; + aDist ++; + + if (anArc==TERMINAL) + { + aNode -> TimeStamp = myTime; + aNode -> Distance = 1; + break; + } + + if (anArc==ORPHAN) { aDist = INFINITE_D; break; } + + if (IS_ODD (anArc)) + { aNode = NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc) -> Shift); } + else + { aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); } + } + + if (aDistShift); aNode->TimeStamp!=myTime;) + { + aNode -> TimeStamp = myTime; + aNode -> Distance = aDist --; + anArc = aNode->Parent; + + if (IS_ODD (anArc)) + { aNode = NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc) -> Shift); } + else + { aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); } + } + } + } + } + } + + if (theNode->Parent = anArc0Min) + { + theNode -> TimeStamp = myTime; + theNode -> Distance = aDistMin + 1; + } + else + { + theNode -> TimeStamp = 0; + + for (anArc0For=anArc0ForFirst; anArc0For Shift); + + if (!aNode->IsSink && (anArc=aNode->Parent)) + { + if (anArc0For->ReverseResidualCap) { setActive (aNode); } + + if (anArc!=TERMINAL && anArc!=ORPHAN && IS_ODD (anArc) && NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc)->Shift) ==theNode) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + + if (myOrphanLast) { myOrphanLast -> Next = aNodePtr; } + else { myOrphanFirst = aNodePtr; } + + myOrphanLast = aNodePtr; + aNodePtr -> Next = NULL; + } + } + } + + for (anArc0Rev=anArc0RevFirst; anArc0Rev Sister; + aNode = NEIGHBOR_NODE_REV (theNode, anArc0For -> Shift); + + if (!aNode->IsSink && (anArc=aNode->Parent)) + { + if (anArc0For->ResidualCap) { setActive (aNode); } + + if (anArc!=TERMINAL && anArc!=ORPHAN && !IS_ODD (anArc) && NEIGHBOR_NODE (aNode, anArc->Shift) ==theNode) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + + if (myOrphanLast) { myOrphanLast -> Next = aNodePtr; } + else { myOrphanFirst = aNodePtr; } + + myOrphanLast = aNodePtr; + aNodePtr -> Next = NULL; + } + } + } + } +} + +void Graph::processSinkOrphan (Node* theNode) +{ + Node* aNode; + ArcForward* anArc0For, *anArc0ForFirst, *anArc0ForLast; + ArcReverse* anArc0Rev, *anArc0RevFirst, *anArc0RevLast; + ArcForward* anArc0Min = NULL, *anArc; + NodePtr* aNodePtr; + int aDist, aDistMin = INFINITE_D; + + anArc0ForFirst = theNode -> FirstOut; + + if (IS_ODD (anArc0ForFirst)) + { + anArc0ForFirst = (ArcForward*) ( ( (char*) anArc0ForFirst) + 1); + anArc0ForLast = (ArcForward*) ( (anArc0ForFirst ++) -> Shift); + } + else { anArc0ForLast = (theNode + 1) -> FirstOut; } + + anArc0RevFirst = theNode -> FirstIn; + + if (IS_ODD (anArc0RevFirst)) + { + anArc0RevFirst = (ArcReverse*) ( ( (char*) anArc0RevFirst) + 1); + anArc0RevLast = (ArcReverse*) ( (anArc0RevFirst ++) -> Sister); + } + else { anArc0RevLast = (theNode + 1) -> FirstIn; } + + + for (anArc0For=anArc0ForFirst; anArc0ForResidualCap) + { + aNode = NEIGHBOR_NODE (theNode, anArc0For -> Shift); + + if (aNode->IsSink && (anArc=aNode->Parent)) + { + aDist = 0; + + while (true) + { + if (aNode->TimeStamp == myTime) + { + aDist += aNode -> Distance; + break; + } + + anArc = aNode -> Parent; + aDist ++; + + if (anArc==TERMINAL) + { + aNode -> TimeStamp = myTime; + aNode -> Distance = 1; + break; + } + + if (anArc==ORPHAN) { aDist = INFINITE_D; break; } + + if (IS_ODD (anArc)) + { aNode = NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc) -> Shift); } + else + { aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); } + } + + if (aDistShift); aNode->TimeStamp!=myTime;) + { + aNode -> TimeStamp = myTime; + aNode -> Distance = aDist --; + anArc = aNode->Parent; + + if (IS_ODD (anArc)) + { aNode = NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc) -> Shift); } + else + { aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); } + } + } + } + } + + for (anArc0Rev=anArc0RevFirst; anArc0Rev Sister; + + if (anArc0For->ReverseResidualCap) + { + aNode = NEIGHBOR_NODE_REV (theNode, anArc0For -> Shift); + + if (aNode->IsSink && (anArc=aNode->Parent)) + { + aDist = 0; + + while (true) + { + if (aNode->TimeStamp == myTime) + { + aDist += aNode -> Distance; + break; + } + + anArc = aNode -> Parent; + aDist ++; + + if (anArc==TERMINAL) + { + aNode -> TimeStamp = myTime; + aNode -> Distance = 1; + break; + } + + if (anArc==ORPHAN) { aDist = INFINITE_D; break; } + + if (IS_ODD (anArc)) + { aNode = NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc) -> Shift); } + else + { aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); } + } + + if (aDistShift); aNode->TimeStamp!=myTime;) + { + aNode -> TimeStamp = myTime; + aNode -> Distance = aDist --; + anArc = aNode->Parent; + + if (IS_ODD (anArc)) + { aNode = NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc) -> Shift); } + else + { aNode = NEIGHBOR_NODE (aNode, anArc -> Shift); } + } + } + } + } + } + + if (theNode->Parent = anArc0Min) + { + theNode -> TimeStamp = myTime; + theNode -> Distance = aDistMin + 1; + } + else + { + theNode -> TimeStamp = 0; + + for (anArc0For=anArc0ForFirst; anArc0For Shift); + + if (aNode->IsSink && (anArc=aNode->Parent)) + { + if (anArc0For->ResidualCap) { setActive (aNode); } + + if (anArc!=TERMINAL && anArc!=ORPHAN && IS_ODD (anArc) && NEIGHBOR_NODE_REV (aNode, MAKE_EVEN (anArc)->Shift) ==theNode) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + + if (myOrphanLast) { myOrphanLast -> Next = aNodePtr; } + else { myOrphanFirst = aNodePtr; } + + myOrphanLast = aNodePtr; + aNodePtr -> Next = NULL; + } + } + } + + for (anArc0Rev=anArc0RevFirst; anArc0Rev Sister; + aNode = NEIGHBOR_NODE_REV (theNode, anArc0For -> Shift); + + if (aNode->IsSink && (anArc=aNode->Parent)) + { + if (anArc0For->ReverseResidualCap) { setActive (aNode); } + + if (anArc!=TERMINAL && anArc!=ORPHAN && !IS_ODD (anArc) && NEIGHBOR_NODE (aNode, anArc->Shift) ==theNode) + { + aNode -> Parent = ORPHAN; + aNodePtr = myNodePtrBLock -> New(); + aNodePtr -> Ptr = aNode; + + if (myOrphanLast) { myOrphanLast -> Next = aNodePtr; } + else { myOrphanFirst = aNodePtr; } + + myOrphanLast = aNodePtr; + aNodePtr -> Next = NULL; + } + } + } + } +} + +Graph::FlowType Graph::MaximumFlow() +{ + Node* aNode0, *aNode1, *aCurrentNode = NULL, *aSStart, *aTStart; + CapacityType* aCapacityMiddle, *aRevCapacityMiddle; + ArcForward* anArcFor, *anArcForFirst, *anArcForLast; + ArcReverse* anArcRev, *anArcRevFirst, *anArcRevLAst; + NodePtr* aNodePtr, *aNodePtrNext; + + prepareGraph(); + maxflowInit(); + myNodePtrBLock = new DBlock (NODEPTR_BLOCK_SIZE, myErrorFun); + + while (true) + { + if (aNode0=aCurrentNode) + { + aNode0 -> Next = NULL; + + if (!aNode0->Parent) { aNode0 = NULL; } + } + + if (!aNode0) + { + if (! (aNode0 = nextActive())) { break; } + } + + aSStart = NULL; + + anArcForFirst = aNode0 -> FirstOut; + + if (IS_ODD (anArcForFirst)) + { + anArcForFirst = (ArcForward*) ( ( (char*) anArcForFirst) + 1); + anArcForLast = (ArcForward*) ( (anArcForFirst ++) -> Shift); + } + else { anArcForLast = (aNode0 + 1) -> FirstOut; } + + anArcRevFirst = aNode0 -> FirstIn; + + if (IS_ODD (anArcRevFirst)) + { + anArcRevFirst = (ArcReverse*) ( ( (char*) anArcRevFirst) + 1); + anArcRevLAst = (ArcReverse*) ( (anArcRevFirst ++) -> Sister); + } + else { anArcRevLAst = (aNode0 + 1) -> FirstIn; } + + if (!aNode0->IsSink) + { + for (anArcFor=anArcForFirst; anArcForResidualCap) + { + aNode1 = NEIGHBOR_NODE (aNode0, anArcFor -> Shift); + + if (!aNode1->Parent) + { + aNode1 -> IsSink = 0; + aNode1 -> Parent = MAKE_ODD (anArcFor); + aNode1 -> TimeStamp = aNode0 -> TimeStamp; + aNode1 -> Distance = aNode0 -> Distance + 1; + setActive (aNode1); + } + else if (aNode1->IsSink) + { + aSStart = aNode0; + aTStart = aNode1; + aCapacityMiddle = & (anArcFor -> ResidualCap); + aRevCapacityMiddle = & (anArcFor -> ReverseResidualCap); + break; + } + else if (aNode1->TimeStamp <= aNode0->TimeStamp && + aNode1->Distance > aNode0->Distance) + { + aNode1 -> Parent = MAKE_ODD (anArcFor); + aNode1 -> TimeStamp = aNode0 -> TimeStamp; + aNode1 -> Distance = aNode0 -> Distance + 1; + } + } + + if (!aSStart) + for (anArcRev=anArcRevFirst; anArcRev Sister; + + if (anArcFor->ReverseResidualCap) + { + aNode1 = NEIGHBOR_NODE_REV (aNode0, anArcFor -> Shift); + + if (!aNode1->Parent) + { + aNode1 -> IsSink = 0; + aNode1 -> Parent = anArcFor; + aNode1 -> TimeStamp = aNode0 -> TimeStamp; + aNode1 -> Distance = aNode0 -> Distance + 1; + setActive (aNode1); + } + else if (aNode1->IsSink) + { + aSStart = aNode0; + aTStart = aNode1; + aCapacityMiddle = & (anArcFor -> ReverseResidualCap); + aRevCapacityMiddle = & (anArcFor -> ResidualCap); + break; + } + else if (aNode1->TimeStamp <= aNode0->TimeStamp && + aNode1->Distance > aNode0->Distance) + { + aNode1 -> Parent = anArcFor; + aNode1 -> TimeStamp = aNode0 -> TimeStamp; + aNode1 -> Distance = aNode0 -> Distance + 1; + } + } + } + } + else + { + for (anArcFor=anArcForFirst; anArcForReverseResidualCap) + { + aNode1 = NEIGHBOR_NODE (aNode0, anArcFor -> Shift); + + if (!aNode1->Parent) + { + aNode1 -> IsSink = 1; + aNode1 -> Parent = MAKE_ODD (anArcFor); + aNode1 -> TimeStamp = aNode0 -> TimeStamp; + aNode1 -> Distance = aNode0 -> Distance + 1; + setActive (aNode1); + } + else if (!aNode1->IsSink) + { + aSStart = aNode1; + aTStart = aNode0; + aCapacityMiddle = & (anArcFor -> ReverseResidualCap); + aRevCapacityMiddle = & (anArcFor -> ResidualCap); + break; + } + else if (aNode1->TimeStamp <= aNode0->TimeStamp && + aNode1->Distance > aNode0->Distance) + { + aNode1 -> Parent = MAKE_ODD (anArcFor); + aNode1 -> TimeStamp = aNode0 -> TimeStamp; + aNode1 -> Distance = aNode0 -> Distance + 1; + } + } + + for (anArcRev=anArcRevFirst; anArcRev Sister; + + if (anArcFor->ResidualCap) + { + aNode1 = NEIGHBOR_NODE_REV (aNode0, anArcFor -> Shift); + + if (!aNode1->Parent) + { + aNode1 -> IsSink = 1; + aNode1 -> Parent = anArcFor; + aNode1 -> TimeStamp = aNode0 -> TimeStamp; + aNode1 -> Distance = aNode0 -> Distance + 1; + setActive (aNode1); + } + else if (!aNode1->IsSink) + { + aSStart = aNode1; + aTStart = aNode0; + aCapacityMiddle = & (anArcFor -> ResidualCap); + aRevCapacityMiddle = & (anArcFor -> ReverseResidualCap); + break; + } + else if (aNode1->TimeStamp <= aNode0->TimeStamp && + aNode1->Distance > aNode0->Distance) + { + aNode1 -> Parent = anArcFor; + aNode1 -> TimeStamp = aNode0 -> TimeStamp; + aNode1 -> Distance = aNode0 -> Distance + 1; + } + } + } + } + + myTime ++; + + if (aSStart) + { + aNode0 -> Next = aNode0; + aCurrentNode = aNode0; + + augment (aSStart, aTStart, aCapacityMiddle, aRevCapacityMiddle); + + while (aNodePtr=myOrphanFirst) + { + aNodePtrNext = aNodePtr -> Next; + aNodePtr -> Next = NULL; + + while (aNodePtr=myOrphanFirst) + { + myOrphanFirst = aNodePtr -> Next; + aNode0 = aNodePtr -> Ptr; + myNodePtrBLock -> Delete (aNodePtr); + + if (!myOrphanFirst) { myOrphanLast = NULL; } + + if (aNode0->IsSink) { processSinkOrphan (aNode0); } + else { processSourceOrphan (aNode0); } + } + + myOrphanFirst = aNodePtrNext; + } + } + else { aCurrentNode = NULL; } + } + + delete myNodePtrBLock; + + return myFlow; +} + +Graph::TerminalType Graph::Label (NodeId i) +{ + if ( ( (Node*) i)->Parent && ! ( (Node*) i)->IsSink) { return SOURCE; } + + return SINK; +} diff --git a/src/BRepMesh/BRepMesh_MinStCut.hxx b/src/BRepMesh/BRepMesh_MinStCut.hxx new file mode 100644 index 0000000000..3657a8ffad --- /dev/null +++ b/src/BRepMesh/BRepMesh_MinStCut.hxx @@ -0,0 +1,243 @@ +// Created on: 2015-11-27 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef __GRAPH_H__ +#define __GRAPH_H__ + +#include +#include + +#define NODE_BLOCK_SIZE 512 +#define ARC_BLOCK_SIZE 1024 +#define NODEPTR_BLOCK_SIZE 128 + +class Graph +{ +public: + typedef enum + { + SOURCE = 0, + SINK = 1 + } TerminalType; + + typedef double CapacityType; + typedef double FlowType; + typedef void* NodeId; + + Graph (void (*theErrFun) (char*) = NULL); + + ~Graph(); + + NodeId AddNode(); + + void AddEdge (NodeId theFromNode, NodeId theToNode, CapacityType theCapNode, CapacityType theRevCapacity); + + void SetTWeights (NodeId theNode, CapacityType theCapacityToSource, CapacityType theCapacityToSink); + + void AddTWeights (NodeId theNode, CapacityType theCapacityToSource, CapacityType theCapacityToSink); + + TerminalType Label (NodeId i); + + FlowType MaximumFlow(); + +private: + + struct ArcForwardSt; + struct ArcReverseSt; + +#define IS_ODD(a) ((int)(a) & 1) +#define MAKE_ODD(a) ((ArcForward *) ((int)(a) | 1)) +#define MAKE_EVEN(a) ((ArcForward *) ((int)(a) & (~1))) +#define MAKE_ODD_REV(a) ((ArcReverse *) ((int)(a) | 1)) +#define MAKE_EVEN_REV(a) ((ArcReverse *) ((int)(a) & (~1))) + + typedef struct NodeSt + { + ArcForwardSt* FirstOut; + ArcReverseSt* FirstIn; + + ArcForwardSt* Parent; + + NodeSt* Next; + + int TimeStamp; + int Distance; + short IsSink; + + CapacityType TrCapacity; + } Node; + +#define NEIGHBOR_NODE(i, shift) ((Node *) ((char *)(i) + (shift))) +#define NEIGHBOR_NODE_REV(i, shift) ((Node *) ((char *)(i) - (shift))) + typedef struct ArcForwardSt + { + int Shift; + CapacityType ResidualCap; + CapacityType ReverseResidualCap; + } ArcForward; + + typedef struct ArcReverseSt + { + ArcForward* Sister; + } ArcReverse; + + typedef struct NodePtrSt + { + NodeSt* Ptr; + NodePtrSt* Next; + } NodePtr; + + typedef struct NodeBLockSt + { + Node* Current; + struct NodeBLockSt* Next; + Node Nodes[NODE_BLOCK_SIZE]; + } NodeBlock; + +#define last_node LAST_NODE.LAST_NODE + + typedef struct ArcForBlockSt + { + char* Start; + ArcForward* Current; + struct ArcForBlockSt* Next; + ArcForward ArcsFor[ARC_BLOCK_SIZE]; + union + { + ArcForward Dummy; + Node* LAST_NODE; + } LAST_NODE; + } ArcForBlock; + + typedef struct ArcRevBlockSt + { + char* Start; + ArcReverse* Current; + struct ArcRevBlockSt* Next; + ArcReverse ArcsRev[ARC_BLOCK_SIZE]; + union + { + ArcReverse Dummy; + Node* LAST_NODE; + } LAST_NODE; + } ArcRevBlock; + + NodeBlock* myNodeBlockFirst; + ArcForBlock* myArcForBlockFirst; + ArcRevBlock* myArcREvBlockFirst; + DBlock* myNodePtrBLock; + + void (*myErrorFun) (char*); + + FlowType myFlow; + + Node* myQueueFirst[2], *myQueueLast[2]; + NodePtr* myOrphanFirst, *myOrphanLast; + int myTime; + + void setActive (Node* theNode); + Node* nextActive(); + + void prepareGraph(); + void maxflowInit(); + void augment (Node* theSStart, Node* theTStart, CapacityType* theCapMiddle, CapacityType* theRevCapMiddle); + void processSourceOrphan (Node* theNode); + void processSinkOrphan (Node* theNode); +}; + +class Energy : Graph +{ +public: + typedef NodeId Var; + + typedef CapacityType Value; + typedef FlowType TotalValue; + + Energy (void (*theErrFun) (char*) = NULL); + + ~Energy(); + + Var AddVariable(); + + void AddConstant (Value E); + + void AddTerm (Var x, + Value E0, Value E1); + + void AddPairwise (Var x, Var y, + Value E00, Value E01, + Value E10, Value E11); + + TotalValue Minimize(); + + int Label (Var x); + +private: + + TotalValue myEConst; + void (*myErrorFun) (char*); +}; + +inline Energy::Energy (void (*err_function) (char*)) : Graph (err_function) +{ + myEConst = 0; + myErrorFun = err_function; +} + +inline Energy::~Energy() {} + +inline Energy::Var Energy::AddVariable() { return AddNode(); } + +inline void Energy::AddConstant (Value A) { myEConst += A; } + +inline void Energy::AddTerm (Var x, + Value A, Value B) +{ + AddTWeights (x, B, A); +} + +inline void Energy::AddPairwise (Var x, Var y, + Value A, Value B, + Value C, Value D) +{ + AddTWeights (x, D, A); + B -= A; + C -= D; + + assert (B + C >= 0); + + if (B < 0) + { + AddTWeights (x, 0, B); + AddTWeights (y, 0, -B); + AddEdge (x, y, 0, B+C); + } + else if (C < 0) + { + AddTWeights (x, 0, -C); + AddTWeights (y, 0, C); + AddEdge (x, y, B+C, 0); + } + else + { + AddEdge (x, y, B, C); + } +} + +inline Energy::TotalValue Energy::Minimize() { return myEConst + MaximumFlow(); } + +inline int Energy::Label (Var x) { return (int) Graph::Label (x); } + +#endif diff --git a/src/BRepMesh/BRepMesh_RestoreOrientationTool.cxx b/src/BRepMesh/BRepMesh_RestoreOrientationTool.cxx new file mode 100644 index 0000000000..0866fa41dc --- /dev/null +++ b/src/BRepMesh/BRepMesh_RestoreOrientationTool.cxx @@ -0,0 +1,881 @@ +// Created on: 2015-08-27 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include + +#include + +#include +#include + +#include + +IMPLEMENT_STANDARD_HANDLE (BRepMesh_RestoreOrientationTool, Standard_Transient) +IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_RestoreOrientationTool, Standard_Transient) + +namespace +{ + +// The maximum positive double value. +const Standard_Real THE_MAX_VAL = std::numeric_limits::max(); + +// The 'small' (safe for inversion) 4D vector of doubles. +const BVH_Vec3d THE_MIN_VEC (1e-100, 1e-100, 1e-100); + +// Relative weight of coherency component for optimization. +const Standard_Real COHERENCY_WEIGHT = 50.0; + +// Minimum number of traced rays per face. +const Standard_Integer MIN_FACE_RAYS = 50; + +// Maximum number of traced rays per face. +const Standard_Integer MAX_FACE_RAYS = 2000; + +// Indicates that ray does not hit any geometry. +const Standard_Integer INVALID_HIT = -1; + +//! Samples value from array according to probabilities. +Standard_Integer sampleValue (Standard_Real theKsi, const std::vector& theCDF) +{ + return static_cast (std::lower_bound (theCDF.begin(), theCDF.end(), theKsi) - theCDF.begin()); +} + +//! Utility class to manage coherence table. +class CoherencyTable +{ +public: + + //! Allocates square coherency table with theMaxSize * theMaxSize maximum number of elements. + CoherencyTable (Standard_Size theMaxSize) + : myMaxSize (theMaxSize), + mySize (0), + myTable (theMaxSize * theMaxSize, 0.0) + { + // + } + + //! Sets coherence value. + void setCoherence (Standard_Size theI, Standard_Size theJ, Standard_Real theValue) + { + myTable[theI * myMaxSize + theJ] = static_cast (theValue); + } + + //! Returns coherence value. + Standard_Real getCoherence (size_t theI, size_t theJ) + { + return static_cast (myTable[theI * myMaxSize + theJ]); + } + + //! Returns actual table size. + Standard_Size Size() + { + return mySize; + } + + //! Sets actual table size. + void SetSize (Standard_Size theSize) + { + mySize = theSize; + } + +private: + + Standard_Size myMaxSize; + Standard_Size mySize; + + std::vector myTable; +}; + +//! Returns time of intersection of ray and triangle. +Standard_Real rayTriangleIntersection (const BVH_Vec3d& theOrigin, + const BVH_Vec3d& theDirect, + const BVH_Vec3d& thePoint0, + const BVH_Vec3d& thePoint1, + const BVH_Vec3d& thePoint2) +{ + const BVH_Vec3d aEdge0 = thePoint1 - thePoint0; + const BVH_Vec3d aEdge1 = thePoint0 - thePoint2; + + BVH_Vec3d aNormal = BVH_Vec3d::Cross (aEdge1, aEdge0); + const Standard_Real aDot = aNormal.Dot (theDirect); + if (fabs (aDot) < 1.0e-16) + { + return THE_MAX_VAL; // ray in the plane of triangle + } + + const BVH_Vec3d aEdge2 = (thePoint0 - theOrigin) * (1.0 / aDot); + + const Standard_Real aTime = aNormal.Dot (aEdge2); + if (aTime < 0.0) + { + return THE_MAX_VAL; + } + + const BVH_Vec3d theInc = BVH_Vec3d::Cross (theDirect, aEdge2); + + Standard_Real U = theInc.Dot (aEdge1); + Standard_Real V = theInc.Dot (aEdge0); + + if ((U >= 0.0) && (V >= 0.0) && (U + V <= 1.0)) + { + return aTime; + } + + return THE_MAX_VAL; +} + +//! Checks if two boxes (AABBs) are overlapped. +bool overlapBoxes (const BVH_Vec3d& theBoxMin1, + const BVH_Vec3d& theBoxMax1, + const BVH_Vec3d& theBoxMin2, + const BVH_Vec3d& theBoxMax2) { + + return !(theBoxMin1.x() > theBoxMax2.x() + || theBoxMax1.x() < theBoxMin2.x() + || theBoxMin1.y() > theBoxMax2.y() + || theBoxMax1.y() < theBoxMin2.y() + || theBoxMin1.z() > theBoxMax2.z() + || theBoxMax1.z() < theBoxMin2.z()); +} + +} + +// ======================================================================= +// function : BRepMesh_RestoreOrientationTool +// purpose : +// ======================================================================= +BRepMesh_RestoreOrientationTool::BRepMesh_RestoreOrientationTool (const bool theVisibilityOnly /* false */) + : myIsDone (false), + myVisibilityOnly (theVisibilityOnly) +{ + // +} + +// ======================================================================= +// function : BRepMesh_RestoreOrientationTool +// purpose : +// ======================================================================= +BRepMesh_RestoreOrientationTool::BRepMesh_RestoreOrientationTool (const TopoDS_Shape& theShape, + const bool theVisibilityOnly /* false */) + : myIsDone (false), + myVisibilityOnly (theVisibilityOnly) +{ + Init (theShape); +} + +// ======================================================================= +// function : Init +// purpose : +// ======================================================================= +void BRepMesh_RestoreOrientationTool::Init (const TopoDS_Shape& theShape) +{ + if (theShape.IsNull()) + { + Standard_NullObject::Raise(); + } + + myFaceList.clear(); + + for (TopExp_Explorer anIter (theShape, TopAbs_FACE); anIter.More(); anIter.Next()) + { + myFaceList.push_back (static_cast (anIter.Current())); + } + + myTriangulation.Elements.clear(); + myTriangulation.Vertices.clear(); + + Standard_Integer aCurrentTriangle = 0; + + for (BRepMesh_FaceList::iterator anIter = myFaceList.begin(); anIter != myFaceList.end(); ++anIter) + { + Handle (BRepMesh_TriangulatedPatch) aNewPatch (new BRepMesh_TriangulatedPatch(*anIter)); + + Standard_Integer aTriangleCount = extractTriangulation (*anIter, *aNewPatch); + + if (aTriangleCount == 0) + { + continue; + } + + aCurrentTriangle += aTriangleCount; + myFaceIdIntervals.push_back (aCurrentTriangle - 1); + aNewPatch->PreprocessData(); + + myPatches.push_back (aNewPatch); + + extractTriangulation (*anIter, myTriangulation); + } +} + +// ======================================================================= +// function : extractTriangulation +// purpose : +// ======================================================================= +Standard_Integer BRepMesh_RestoreOrientationTool::extractTriangulation (const TopoDS_Face& theFace, + BVH_Triangulation& theTriangulation) +{ + TopLoc_Location aLocation; + + Handle (Poly_Triangulation) aPolyTriangulation = BRep_Tool::Triangulation (theFace, aLocation); + + if (aPolyTriangulation.IsNull()) + { + return 0; + } + + bool isMirrored = aLocation.Transformation().VectorialPart().Determinant() < 0; + bool needToFlip = static_cast (isMirrored ^ (theFace.Orientation() == TopAbs_REVERSED)); + + BRepAdaptor_Surface aFaceAdaptor (theFace, Standard_False); + + const Standard_Integer aVertOffset = + static_cast (theTriangulation.Vertices.size()) - 1; + + for (Standard_Integer aVertIdx = 1; aVertIdx <= aPolyTriangulation->NbNodes(); ++aVertIdx) + { + gp_Pnt aVertex = aPolyTriangulation->Nodes().Value (aVertIdx); + + aVertex.Transform (aLocation.Transformation()); + + theTriangulation.Vertices.push_back (BVH_Vec3d (aVertex.X(), + aVertex.Y(), + aVertex.Z())); + } + + for (Standard_Integer aTriIdx = 1; aTriIdx <= aPolyTriangulation->NbTriangles(); ++aTriIdx) + { + Standard_Integer aVertex1; + Standard_Integer aVertex2; + Standard_Integer aVertex3; + + if (!needToFlip) + { + aPolyTriangulation->Triangles().Value (aTriIdx).Get (aVertex1, + aVertex2, + aVertex3); + } + else + { + aPolyTriangulation->Triangles().Value (aTriIdx).Get (aVertex3, + aVertex2, + aVertex1); + } + + theTriangulation.Elements.push_back (BVH_Vec4i (aVertex1 + aVertOffset, + aVertex2 + aVertOffset, + aVertex3 + aVertOffset, + 1)); + } + + theTriangulation.MarkDirty(); // triangulation needs BVH rebuilding + + return aPolyTriangulation->NbTriangles(); +} + +// ======================================================================= +// function : computeVisibility +// purpose : +// ======================================================================= +void BRepMesh_RestoreOrientationTool::computeVisibility (BVH_Triangulation& theTriangulation, + Handle (BRepMesh_TriangulatedPatch)& thePatch, + Standard_Integer theRaysNb) +{ + // deterministic sequence + math_BullardGenerator anRNG (0); + + Standard_Real aFrontVisibility = 0.0; + Standard_Real aBackVisibility = 0.0; + + Standard_Real aFrontDistance = 0.0; + Standard_Real aBackDistance = 0.0; + + for (Standard_Integer anIdx = 0; anIdx < theRaysNb; ++anIdx) + { + // select triangle with probability proportional to triangle area + Standard_Integer aTrgIdx = sampleValue (anRNG.NextReal(), thePatch->TriangleCDF()); + + const BVH_Vec4i& aTriangle = thePatch->Elements[aTrgIdx]; + + // sample point in triangle uniformly + Standard_Real aKsi, aPsi; + do + { + aKsi = anRNG.NextReal(); + aPsi = anRNG.NextReal(); + } while (aKsi + aPsi > 1.0); + + const BVH_Vec3d& aPoint0 = thePatch->Vertices[aTriangle.x()]; + const BVH_Vec3d& aPoint1 = thePatch->Vertices[aTriangle.y()]; + const BVH_Vec3d& aPoint2 = thePatch->Vertices[aTriangle.z()]; + + BVH_Vec3d aPoint (aPoint0 * (1 - aKsi - aPsi) + + aPoint1 * aKsi + + aPoint2 * aPsi); + + const BVH_Vec3d aEdge0 = aPoint1 - aPoint0; + const BVH_Vec3d aEdge1 = aPoint0 - aPoint2; + + const BVH_Vec3d aNormal = (BVH_Vec3d::Cross (aEdge1, aEdge0).Normalized()); + + // hemisphere direction + BVH_Vec3d aDirection; + Standard_Real aDot = 0.0; + + do + { + aKsi = anRNG.NextReal(); + aPsi = anRNG.NextReal(); + + Standard_Real anSqrtPsi = sqrt (aPsi); + + aDirection = BVH_Vec3d (anSqrtPsi * cos (2.f * M_PI * aKsi), + anSqrtPsi * sin (2.f * M_PI * aKsi), + sqrt (1.f - aPsi)); + + aDot = aDirection.Dot (aNormal); + + } while (std::abs (aDot) < 0.2); + + if (aDot < 0.0) + { + aDirection = -aDirection; + } + + const Standard_Real anEpsilon = std::max (1.0e-6, 1.0e-4 * theTriangulation.Box().Size().Modulus()); + + Standard_Real aDist = 0.0; + if (traceRay (aPoint + aNormal * anEpsilon, aDirection, theTriangulation, aDist) == INVALID_HIT) + { + aFrontVisibility += 1.0; + } + else + { + aFrontDistance += aDist; + } + + if (traceRay (aPoint - aNormal * anEpsilon, -aDirection, theTriangulation, aDist) == INVALID_HIT) + { + aBackVisibility += 1.0; + } + else + { + aBackDistance += aDist; + } + } + + Standard_Real aInvRayNb = 1.0 / theRaysNb; + thePatch->SetFrontVisibility (aFrontVisibility * aInvRayNb); + thePatch->SetBackVisibility (aBackVisibility * aInvRayNb); + thePatch->SetFrontDistance (aFrontDistance * aInvRayNb); + thePatch->SetBackDistance (aBackDistance * aInvRayNb); +} + +// ======================================================================= +// function : traceRay +// purpose : +// ======================================================================= +Standard_Integer BRepMesh_RestoreOrientationTool::traceRay (const BVH_Vec3d& theOrigin, + const BVH_Vec3d& theDirect, + BVH_Triangulation& theTriangulation, + Standard_Real& theDistance) +{ + Standard_Integer aStack[32]; + + const NCollection_Handle >& aBVH = theTriangulation.BVH(); + + if (aBVH.IsNull()) + { + return INVALID_HIT; + } + + Standard_Integer aNode = 0; + + Standard_Real aMinDist = std::numeric_limits::max(); + Standard_Integer aHitTriangle = -1; + + Standard_Integer aHead = -1; // stack head position + + BVH_Vec3d anInvDirect = theDirect.cwiseAbs().cwiseMax (THE_MIN_VEC); + anInvDirect = BVH_Vec3d (1.0 / anInvDirect.x(), 1.0 / anInvDirect.y(), 1.0 / anInvDirect.z()); + + anInvDirect = BVH_Vec3d (_copysign (anInvDirect.x(), theDirect.x()), + _copysign (anInvDirect.y(), theDirect.y()), + _copysign (anInvDirect.z(), theDirect.z())); + + for (;;) + { + if (aBVH->IsOuter (aNode)) + { + for (Standard_Integer anElemNb = aBVH->BegPrimitive (aNode); anElemNb <= aBVH->EndPrimitive (aNode); ++anElemNb) + { + const BVH_Vec4i& aTriangle = theTriangulation.Elements[anElemNb]; + const BVH_Vec3d& aVrt0 = theTriangulation.Vertices[aTriangle.x()]; + const BVH_Vec3d& aVrt1 = theTriangulation.Vertices[aTriangle.y()]; + const BVH_Vec3d& aVrt2 = theTriangulation.Vertices[aTriangle.z()]; + + Standard_Real aDist = rayTriangleIntersection (theOrigin, theDirect, aVrt0, aVrt1, aVrt2); + + + if (aDist < aMinDist) + { + aMinDist = aDist; + aHitTriangle = anElemNb; + } + } + + if (aHead < 0) + { + break; + } + + aNode = aStack[aHead--]; + } + else + { + BVH_Vec3d aTime0 = (aBVH->MinPoint (aBVH->LeftChild (aNode)) - theOrigin) * anInvDirect; + BVH_Vec3d aTime1 = (aBVH->MaxPoint (aBVH->LeftChild (aNode)) - theOrigin) * anInvDirect; + + BVH_Vec3d aTimeMax = aTime0.cwiseMax (aTime1); + BVH_Vec3d aTimeMin = aTime0.cwiseMin (aTime1); + + aTime0 = (aBVH->MinPoint (aBVH->RightChild (aNode)) - theOrigin) * anInvDirect; + aTime1 = (aBVH->MaxPoint (aBVH->RightChild (aNode)) - theOrigin) * anInvDirect; + + Standard_Real aTimeFinal = aTimeMax.minComp(); + Standard_Real aTimeStart = aTimeMin.maxComp(); + + aTimeMax = aTime0.cwiseMax (aTime1); + aTimeMin = aTime0.cwiseMin (aTime1); + + Standard_Real aTimeMin1 = (aTimeStart <= aTimeFinal) && (aTimeFinal >= 0.0) + ? aTimeStart : THE_MAX_VAL; + + aTimeFinal = aTimeMax.minComp(); + aTimeStart = aTimeMin.maxComp(); + + Standard_Real aTimeMin2 = (aTimeStart <= aTimeFinal) && (aTimeFinal >= 0.f) + ? aTimeStart : THE_MAX_VAL; + + const bool aHitLft = (aTimeMin1 != THE_MAX_VAL); + const bool aHitRgh = (aTimeMin2 != THE_MAX_VAL); + + if (aHitLft && aHitRgh) + { + int aLeft = aBVH->LeftChild (aNode); + int aRight = aBVH->RightChild (aNode); + + aNode = (aTimeMin1 < aTimeMin2) ? aLeft : aRight; + aStack[++aHead] = (aTimeMin2 <= aTimeMin1) ? aLeft : aRight; + } + else + { + if (aHitLft || aHitRgh) + { + aNode = aHitLft ? aBVH->LeftChild (aNode) : aBVH->RightChild (aNode); + } + else + { + if (aHead < 0) + { + break; + } + + aNode = aStack[aHead--]; + } + } + } + } + + theDistance = aMinDist; + + return aHitTriangle; +} + +// ======================================================================= +// function : edgeCoherence +// purpose : +// ======================================================================= +Standard_Real BRepMesh_RestoreOrientationTool::edgeCoherence (const BRepMesh_OrientedEdge& theEdge1, + const BRepMesh_TriangulatedPatch& thePatch1, + const BRepMesh_OrientedEdge& theEdge2, + const BRepMesh_TriangulatedPatch& thePatch2) +{ + BVH_Vec3d anE1V1; + BVH_Vec3d anE1V2; + BVH_Vec3d anE2V1; + BVH_Vec3d anE2V2; + + thePatch1.GetEdgeVertices (theEdge1, anE1V1, anE1V2); + thePatch2.GetEdgeVertices (theEdge2, anE2V1, anE2V2); + + BVH_Vec3d anEdge1 = anE1V2 - anE1V1; + BVH_Vec3d anEdge2 = anE2V2 - anE2V1; + + Standard_Real aScalarProduct = anEdge1.Dot (anEdge2); + + Standard_Real aDist = BRepMesh_TriangulatedPatch::EdgeSet::EdgeToEdgeDistance (anE1V1, anE1V2, anE2V1, anE2V2); + Standard_Real aSign = -_copysign (1.0, aScalarProduct); + + return aSign * sqrt (std::abs (aScalarProduct)) / (1.0 + aDist); +} + +// ======================================================================= +// function : computeCoherence +// purpose : +// ======================================================================= +Standard_Real BRepMesh_RestoreOrientationTool::computeCoherence (Handle (BRepMesh_TriangulatedPatch)& thePatch0, + Handle (BRepMesh_TriangulatedPatch)& thePatch1) +{ + Standard_Real aTotalCoherence = 0.0; + + // Ensure that BVH already built before parallel section + thePatch0->BoundaryEdgeSet().BVH(); + thePatch1->BoundaryEdgeSet().BVH(); + + #pragma omp parallel for + for (Standard_Integer anIdx = 0; anIdx < thePatch1->BoundaryEdgeSet().Size(); ++anIdx) + { + BRepMesh_OrientedEdge aCurrentEdge = thePatch1->BoundaryEdge (anIdx); + BVH_Vec3d aVertex1; + BVH_Vec3d aVertex2; + thePatch1->GetEdgeVertices (aCurrentEdge, aVertex1, aVertex2); + BRepMesh_OrientedEdge aClosestEdge = thePatch0->BoundaryEdge ( + thePatch0->BoundaryEdgeSet().CoherentEdge (aVertex1, aVertex2)); + + #pragma omp critical + { + aTotalCoherence += edgeCoherence (aCurrentEdge, *thePatch1, aClosestEdge, *thePatch0); + } + } + + return aTotalCoherence; +} + +// ======================================================================= +// function : Perform +// purpose : +// ======================================================================= +void BRepMesh_RestoreOrientationTool::Perform() +{ + Standard_Real aMaxArea = std::numeric_limits::epsilon(); + +#ifdef MEASURE_PERFORMANCE + OSD_Timer aTimer; + aTimer.Start(); +#endif + + std::auto_ptr aTable; + + // Flipped flags for all patches. + std::vector aFlipped (myPatches.size(), 0); + std::vector aFlippedC (myPatches.size(), 0); + + Standard_Real aMaxCoherence = -std::numeric_limits::max(); + + const Standard_Real anEpsilon = std::max (1.0e-4, 1.0e-2 * myTriangulation.Box().Size().Modulus()); + BVH_Vec3d anEpsilonVec (anEpsilon, anEpsilon, anEpsilon); + + if (!myVisibilityOnly) + { + aTable.reset (new CoherencyTable (myPatches.size() * 3)); + aTable->SetSize (myPatches.size()); + + // Compute coherence + for (Standard_Size i = 0; i < myPatches.size() - 1; ++i) + { + Handle (BRepMesh_TriangulatedPatch)& aTriPatch1 = myPatches[i]; + + for (Standard_Size j = i + 1; j < myPatches.size(); ++j) + { + Handle (BRepMesh_TriangulatedPatch)& aTriPatch2 = myPatches[j]; + + if (overlapBoxes (aTriPatch1->Box().CornerMin() - anEpsilonVec, + aTriPatch1->Box().CornerMax() + anEpsilonVec, + aTriPatch2->Box().CornerMin() - anEpsilonVec, + aTriPatch2->Box().CornerMax() + anEpsilonVec)) + { + Standard_Real aCoherence = computeCoherence (aTriPatch1, aTriPatch2) + computeCoherence (aTriPatch2, aTriPatch1); + + aMaxCoherence = std::max (aMaxCoherence, aCoherence); + + aTable->setCoherence (i, j, aCoherence); + aTable->setCoherence (j, i, aCoherence); + } + } + } + + std::vector aMetaPatches; + aMetaPatches.reserve (myPatches.size() * 3); + + // Prepare meta patches + for (Standard_Size i = 0; i < myPatches.size(); ++i) + { + aMetaPatches.push_back (BRepMesh_SuperPatch (static_cast (i))); + } + + myPatchQueue.clear(); + myPatchQueue.reserve (myPatches.size()); + + for (Standard_Size i = 0; i < myPatches.size() - 1; ++i) + { + for (Standard_Size j = i + 1; j < myPatches.size(); ++j) + { + Standard_Real aCoherence = aTable->getCoherence (i, j); + if (aCoherence != 0.0) + { + myPatchQueue.push_back (BRepMesh_SuperPatchPair (&aMetaPatches[i], &aMetaPatches[j], aTable->getCoherence (i, j))); + } + } + } + + std::set aRemovedSuperPatches; + + // Greedy optimization of coherence + while (!myPatchQueue.empty()) + { + std::sort (myPatchQueue.begin(), myPatchQueue.end()); + + const BRepMesh_SuperPatchPair& aBestPair = myPatchQueue.back(); + + if (aBestPair.Coherence < 0) + { + // Do flip of second patch + Standard_Integer anIndex = aBestPair.SecondPatch->Index; + + for (Standard_Size j = 0; j < aBestPair.SecondPatch->PatchIds.size(); ++j) + { + Standard_Integer aPatchId = aBestPair.SecondPatch->PatchIds[j]; + + if (aPatchId == anIndex) + { + aFlippedC[aPatchId] = aFlippedC[aPatchId] == 0 ? 1 : 0; + } + + aFlipped[aPatchId] = aFlipped[aPatchId] == 0 ? 1 : 0; + } + + for (Standard_Size i = 0; i < aTable->Size(); ++i) + { + Standard_Real aCoherence = aTable->getCoherence (i, anIndex); + + aTable->setCoherence (i, anIndex, -aCoherence); + aTable->setCoherence (anIndex, i, -aCoherence); + } + } + + aRemovedSuperPatches.insert (aBestPair.FirstPatch->Index); + aRemovedSuperPatches.insert (aBestPair.SecondPatch->Index); + + aMetaPatches.push_back (BRepMesh_SuperPatch (*aBestPair.FirstPatch, *aBestPair.SecondPatch, (Standard_Integer)aTable->Size())); + + const BRepMesh_SuperPatch& aNewPatch = aMetaPatches.back(); + + // Update coherence + for (Standard_Size i = 0; i < aTable->Size(); ++i) + { + if (aRemovedSuperPatches.find ((Standard_Integer)i) == aRemovedSuperPatches.end()) + { + Standard_Real aNewCoherence = 0.0; + aNewCoherence += aTable->getCoherence (i, aBestPair.FirstPatch->Index); + aNewCoherence += aTable->getCoherence (i, aBestPair.SecondPatch->Index); + + aTable->setCoherence (i, aNewPatch.Index, aNewCoherence); + aTable->setCoherence (aNewPatch.Index, i, aNewCoherence); + } + } + + std::vector aNewPatchQueue; + aNewPatchQueue.reserve (myPatchQueue.size()); + + // Update pairs + for (Standard_Size i = 0; i < myPatchQueue.size(); ++i) + { + if (!myPatchQueue[i].HavePatch (aBestPair.FirstPatch->Index) + && !myPatchQueue[i].HavePatch (aBestPair.SecondPatch->Index)) + { + // Accept a pair + aNewPatchQueue.push_back (myPatchQueue[i]); + } + } + + myPatchQueue.swap (aNewPatchQueue); + + // Add new pairs + for (Standard_Size i = 0; i < aTable->Size(); ++i) + { + if (aRemovedSuperPatches.find ((Standard_Integer)i) == aRemovedSuperPatches.end()) + { + Standard_Real aCoherence = aTable->getCoherence (i, aNewPatch.Index); + + if (aCoherence != 0.0) + { + myPatchQueue.push_back (BRepMesh_SuperPatchPair (&aNewPatch, &aMetaPatches[i], aCoherence)); + } + } + } + + aTable->SetSize (aTable->Size() + 1); + } + +#ifdef MEASURE_PERFORMANCE + aTimer.Stop(); + std::cout << "Coherency time: " << aTimer.ElapsedTime() << std::endl; + aTimer.Reset(); + aTimer.Start(); +#endif + } + + // Compute max area + for (Standard_Size i = 0; i < myPatches.size(); ++i) + { + aMaxArea = std::max (aMaxArea, myPatches[i]->TotalArea()); + } + + const Standard_Real anInvMaxArea = 1.0 / aMaxArea; + + const Standard_Integer aPatchesNb = static_cast (myPatches.size()); + + // Ensure that BVH already built before parallel section + myTriangulation.BVH(); + + // Compute visibility + #pragma omp parallel for + for (Standard_Integer i = 0; i < aPatchesNb; ++i) + { + Handle (BRepMesh_TriangulatedPatch)& aTriPatch = myPatches[i]; + Standard_Integer aRaysNb = Max (MIN_FACE_RAYS, static_cast (aTriPatch->TotalArea() * anInvMaxArea * MAX_FACE_RAYS)); + computeVisibility (myTriangulation, aTriPatch, aRaysNb); + + if (!myVisibilityOnly && aFlipped[i]) + { + aTriPatch->FlipVisibility(); + } + } + +#ifdef MEASURE_PERFORMANCE + aTimer.Stop(); + std::cout << "Visibility time: " << aTimer.ElapsedTime() << std::endl; + aTimer.Reset(); + aTimer.Start(); +#endif + + // Optimization + Energy* aGraph = new Energy(); + + std::vector aVariables (myPatches.size()); + + for (Standard_Size i = 0; i < myPatches.size(); ++i) + { + Handle (BRepMesh_TriangulatedPatch)& aTriPatch = myPatches[i]; + aVariables[i] = aGraph->AddVariable(); + + aGraph->AddTerm (aVariables[i], aTriPatch->BackVisibility() - aTriPatch->FrontVisibility(), + aTriPatch->FrontVisibility() - aTriPatch->BackVisibility()); + } + + if (!myVisibilityOnly) + { + const Standard_Real aInvMaxCoherence = 1.0 / aMaxCoherence; + + for (Standard_Size i = 0; i < myPatches.size() - 1; ++i) + { + Handle (BRepMesh_TriangulatedPatch)& aTriPatch1 = myPatches[i]; + for (Standard_Size j = i + 1; j < myPatches.size(); ++j) + { + Handle (BRepMesh_TriangulatedPatch)& aTriPatch2 = myPatches[j]; + + if (!overlapBoxes (aTriPatch1->Box().CornerMin() - anEpsilonVec, + aTriPatch1->Box().CornerMax() + anEpsilonVec, + aTriPatch2->Box().CornerMin() - anEpsilonVec, + aTriPatch2->Box().CornerMax() + anEpsilonVec)) + { + continue; + } + + bool isCoherenceFlipped = (aFlipped[i] ^ aFlipped[j] ^ aFlippedC[i] ^ aFlippedC[j]) != 0; + Standard_Real aCoherence = aTable->getCoherence (i, j) * aInvMaxCoherence * (isCoherenceFlipped ? -1.0 : 1.0); + + if (aCoherence > 0.0) + { + aCoherence *= COHERENCY_WEIGHT; + + // Prevent setting different labels for already coherent patches. + aGraph->AddPairwise (aVariables[i], aVariables[j], 0, aCoherence, aCoherence, 0); + } + } + } + } + + aGraph->Minimize(); + + for (Standard_Size i = 0; i < myPatches.size(); ++i) + { + if (aGraph->Label (aVariables[i]) == 1) + { + aFlipped[i] = aFlipped[i] == 0 ? 1 : 0; + } + } + +#ifdef MEASURE_PERFORMANCE + aTimer.Stop(); + std::cout << "Optimization time: " << aTimer.ElapsedTime() << std::endl; +#endif + + Handle (BRepMesh_TriangulatedPatch) aMergedPatch = new BRepMesh_TriangulatedPatch(myFaceList[0]); + for (Standard_Size i = 0; i < myPatches.size(); ++i) + { + Handle (BRepMesh_TriangulatedPatch)& aTriPatch = myPatches[i]; + + if (aFlipped[i]) + { + aTriPatch->Flip(); + } + + aMergedPatch->Append (aTriPatch); + } + + BRep_Builder aBuilder; + TopoDS_Compound aCompound; + aBuilder.MakeCompound (aCompound); + for (Standard_Size i = 0; i < myPatches.size(); ++i) + { + TopoDS_Face aFace = myPatches[i]->Face(); + if (aFlipped[i]) + { + aFace.Reverse(); + } + aBuilder.Add (aCompound, aFace); + } + + myShape = aCompound; + myIsDone = true; +} diff --git a/src/BRepMesh/BRepMesh_RestoreOrientationTool.hxx b/src/BRepMesh/BRepMesh_RestoreOrientationTool.hxx new file mode 100644 index 0000000000..1704fc8bb0 --- /dev/null +++ b/src/BRepMesh/BRepMesh_RestoreOrientationTool.hxx @@ -0,0 +1,189 @@ +// Created on: 2015-08-27 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _BRepMesh_RestoreOrientationTool_HeaderFile +#define _BRepMesh_RestoreOrientationTool_HeaderFile + +#include +#include +#include + +#include + +//! Class representing a collection of simple patches. +struct BRepMesh_SuperPatch +{ + //! Create super patch from single patch. + BRepMesh_SuperPatch (Standard_Integer thePatchId) + : Index (thePatchId) + { + PatchIds.push_back (thePatchId); + } + + //! Create super patch from a pair of super patches. + BRepMesh_SuperPatch (const BRepMesh_SuperPatch& theFirst, + const BRepMesh_SuperPatch& theSecond, + Standard_Integer theIndex) + : Index (theIndex) + { + PatchIds.resize (theFirst.PatchIds.size() + theSecond.PatchIds.size()); + + std::vector::iterator aNextPos = + std::copy (theFirst.PatchIds.begin(), theFirst.PatchIds.end(), PatchIds.begin()); + + std::copy (theSecond.PatchIds.begin(), theSecond.PatchIds.end(), aNextPos); + } + + std::vector PatchIds; + Standard_Integer Index; +}; + +//! Class representing a pair of super patches. +struct BRepMesh_SuperPatchPair +{ + BRepMesh_SuperPatchPair (const BRepMesh_SuperPatch* theFirstPatch = NULL, + const BRepMesh_SuperPatch* theSecondPatch = NULL, + Standard_Real theCoherency = 0.0) + : FirstPatch (theFirstPatch), + SecondPatch (theSecondPatch), + Coherence (theCoherency) + { + // + } + + bool operator < (const BRepMesh_SuperPatchPair& theOther) const + { + return (std::abs (Coherence) < std::abs (theOther.Coherence)); + } + + bool HavePatch (Standard_Integer thePatch) + { + return FirstPatch->Index == thePatch || SecondPatch->Index == thePatch; + } + + const BRepMesh_SuperPatch* FirstPatch; + const BRepMesh_SuperPatch* SecondPatch; + Standard_Real Coherence; +}; + +//! List of faces. +typedef std::vector BRepMesh_FaceList; + +//! This tool intended to restore consistent orientation of surfaces. +//! It takes as input OCCT shape and outputs a set of re-oriented faces. +class BRepMesh_RestoreOrientationTool : public Standard_Transient +{ +public: + + DEFINE_STANDARD_ALLOC + + //! Creates uninitialized orientation tool. + Standard_EXPORT BRepMesh_RestoreOrientationTool (const bool theVisibilityOnly = false); + + //! Creates orientation tool from the given shape. + Standard_EXPORT BRepMesh_RestoreOrientationTool (const TopoDS_Shape& theShape, + const bool theVisibilityOnly = false); + +public: + + //! Loads the given shape into orientation tool. + Standard_EXPORT void Init (const TopoDS_Shape& theShape); + + //! Performs restoring of consistent orientation. + Standard_EXPORT void Perform(); + + //! Returns "Visibility only" mode. + bool VisibilityOnly() const + { + return myVisibilityOnly; + } + + //! Set "Visibility only" mode. + void SetVisibilityOnly (bool isEnabled) + { + myVisibilityOnly = isEnabled; + } + + //! True if the algorithm successfully performed. + bool IsDone() const + { + return myIsDone; + } + + //! Returns shape with restored normal orientation. + TopoDS_Shape Shape() const + { + return myShape; + } + +protected: + + //! Calculates the coherence value between two patches. + static Standard_Real computeCoherence (Handle (BRepMesh_TriangulatedPatch)& thePatch0, + Handle (BRepMesh_TriangulatedPatch)& thePatch1); + + //! Calculates the visibility value for the patch. + void computeVisibility (BVH_Triangulation& theTriangulation, + Handle (BRepMesh_TriangulatedPatch)& thePatch, + Standard_Integer theRaysNb); + + //! Returns the closest triangle hit by the ray. + Standard_Integer traceRay (const BVH_Vec3d& theOrigin, + const BVH_Vec3d& theDirect, + BVH_Triangulation& theTriangulation, + Standard_Real& theDistance); + + //! Extracts a triangulation from the face. + static Standard_Integer extractTriangulation (const TopoDS_Face& theShape, + BVH_Triangulation& theTriangulation); + + //! Calculates the coherence value between two edges. + static Standard_Real edgeCoherence (const BRepMesh_OrientedEdge& theEdge1, + const BRepMesh_TriangulatedPatch& thePatch1, + const BRepMesh_OrientedEdge& theEdge2, + const BRepMesh_TriangulatedPatch& thePatch2); + +protected: + + //! Is output ready? + bool myIsDone; + + //! Resulted shape. + TopoDS_Shape myShape; + + //! List of triangulated faces of the shape. + BRepMesh_FaceList myFaceList; + + //! List of patches. + std::vector myPatches; + + //! Vertex and normal data of triangulated shape. + BVH_Triangulation myTriangulation; + + std::vector myPatchQueue; + + //! Maps Id intervals to topological faces. + std::vector myFaceIdIntervals; + + //! Use only visibility metric? + bool myVisibilityOnly; + + DEFINE_STANDARD_RTTI (BRepMesh_RestoreOrientationTool) +}; + +DEFINE_STANDARD_HANDLE(BRepMesh_RestoreOrientationTool, Standard_Transient) + + +#endif diff --git a/src/BRepMesh/BRepMesh_TriangulatedPatch.cxx b/src/BRepMesh/BRepMesh_TriangulatedPatch.cxx new file mode 100644 index 0000000000..75d2eea7f4 --- /dev/null +++ b/src/BRepMesh/BRepMesh_TriangulatedPatch.cxx @@ -0,0 +1,233 @@ +// Created on: 2015-08-28 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include + +#include +#include + +IMPLEMENT_STANDARD_HANDLE (BRepMesh_TriangulatedPatch, Standard_Transient) +IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_TriangulatedPatch, Standard_Transient) + +namespace { + + //! Compare functor for symmetrical pairs of integers. + struct SymmetricalPairComparator + { + bool operator() (const std::pair& theEdge1, + const std::pair& theEdge2) const + { + const int aMin1 = std::min (theEdge1.first, theEdge1.second); + const int aMin2 = std::min (theEdge2.first, theEdge2.second); + + return (aMin1 < aMin2) || + (aMin1 == aMin2 && std::max (theEdge1.first, theEdge1.second) < std::max (theEdge2.first, theEdge2.second)); + } + + bool equal (const std::pair& theEdge1, + const std::pair& theEdge2) const + { + return std::min (theEdge1.first, theEdge1.second) == std::min (theEdge2.first, theEdge2.second) + && std::max (theEdge1.first, theEdge1.second) == std::max (theEdge2.first, theEdge2.second); + } + }; + + // Returns the area of triangle + Standard_Real triangleArea (const BVH_Vec3d& theV1, + const BVH_Vec3d& theV2, + const BVH_Vec3d& theV3) + { + BVH_Vec3d aE1 = theV2 - theV1; + BVH_Vec3d aE2 = theV3 - theV1; + + return BVH_Vec3d::Cross (aE1, aE2).Modulus() * 0.5; + } + + //! Set of integers. + typedef std::set SetOfInteger; + + //! Type of edge-triangles connectivity map. + typedef std::map, SetOfInteger, SymmetricalPairComparator> EdgeFaceMap; +} + +// ======================================================================= +// function : BRepMesh_TriangulatedPatch +// purpose : +// ======================================================================= +BRepMesh_TriangulatedPatch::BRepMesh_TriangulatedPatch (TopoDS_Face theFace) + : BVH_Triangulation(), + myBoundaryEdgeSet (Vertices), + myTotalArea (0.0), + myFrontVisibility (0.0), + myBackVisibility (0.0), + myFrontDistance (0.0), + myBackDistance (0.0), + myFace (theFace) +{ + myBuilder = new BVH_LinearBuilder (5 /* leaf size */, 32 /* max height */); +} + +// ======================================================================= +// function : Append +// purpose : +// ======================================================================= +void BRepMesh_TriangulatedPatch::Append (const Handle (BRepMesh_TriangulatedPatch)& thePatch) +{ + // Vertices + Standard_Size aVertOffset = Vertices.size(); + Vertices.resize (Vertices.size() + thePatch->Vertices.size()); + + std::copy (thePatch->Vertices.begin(), thePatch->Vertices.end(), Vertices.begin() + aVertOffset); + + // Triangles + struct + { + Standard_Integer Offset; + BVH_Vec4i operator() (const BVH_Vec4i& theTriangle) { return theTriangle + BVH_Vec4i (Offset, Offset, Offset, 0); } + } AddOffsetHelper; + + AddOffsetHelper.Offset = static_cast (aVertOffset); + + Standard_Size aFaceOffset = Elements.size(); + Elements.resize (Elements.size() + thePatch->Elements.size()); + + std::transform (thePatch->Elements.begin(), thePatch->Elements.end(), Elements.begin() + aFaceOffset, AddOffsetHelper); + + // Area + myTotalArea += thePatch->TotalArea(); +} + +// ======================================================================= +// function : computeBoundareEdges +// purpose : +// ======================================================================= +void BRepMesh_TriangulatedPatch::computeBoundaryEdges() +{ + myBoundaryEdgeSet.MarkDirty(); + + if (!myBoundaryEdgeSet.Elements.empty()) + { + return; + } + + EdgeFaceMap anEdgeFaceMap; + + // use triangle indices for connectivity extraction + for (Standard_Size aTrgIdx = 0; aTrgIdx < Elements.size(); ++aTrgIdx) + { + const BVH_Vec4i& aTriangle = Elements[aTrgIdx]; + + for (Standard_Integer anEdgeStart = 0, edgeEnd = 1; anEdgeStart < 3; ++anEdgeStart, edgeEnd = (anEdgeStart + 1) % 3) + { + std::pair anEdge (aTriangle[anEdgeStart], aTriangle[edgeEnd]); + + anEdgeFaceMap[anEdge].insert (static_cast (aTrgIdx)); + } + } + + // extract boundary and non-manifold edges + myBoundaryEdgeSet.Elements.clear(); + + for (EdgeFaceMap::iterator anIter = anEdgeFaceMap.begin(); anIter != anEdgeFaceMap.end(); ++anIter) + { + const SetOfInteger& anAdjacentFaces = anIter->second; + + const std::pair& anEdge = anIter->first; + + if (anAdjacentFaces.size() == 1) // boundary + { + myBoundaryEdgeSet.Elements.push_back (BVH_Vec4i (anEdge.first, anEdge.second, 0, 1)); + } + } + + myBoundaryEdgeSet.BVH(); +} + +// ======================================================================= +// function : PreprocessData +// purpose : +// ======================================================================= +void BRepMesh_TriangulatedPatch::PreprocessData() +{ + + computeBoundaryEdges(); + updateArea(); +} + +// ======================================================================= +// function : Flip +// purpose : +// ======================================================================= +void BRepMesh_TriangulatedPatch::Flip() +{ + // mark that triangles was flipped + struct + { + void operator() (BVH_Vec4i& theTriangle) { theTriangle.w() = -theTriangle.w(); } + } TriangleFlipHelper; + + std::for_each (Elements.begin(), Elements.end(), TriangleFlipHelper); + + // mark that edges was flipped + struct + { + void operator() (BVH_Vec4i& theEdge) { theEdge.w() = -theEdge.w(); } + } EdgeFlipHelper; + + std::for_each (myBoundaryEdgeSet.Elements.begin(), myBoundaryEdgeSet.Elements.end(), EdgeFlipHelper); + + FlipVisibility(); +} + +// ======================================================================= +// function : FlipVisibility +// purpose : +// ======================================================================= +void BRepMesh_TriangulatedPatch::FlipVisibility() +{ + std::swap (myFrontVisibility, myBackVisibility); + std::swap (myFrontDistance, myBackDistance); +} + +// ======================================================================= +// function : updateArea +// purpose : +// ======================================================================= +void BRepMesh_TriangulatedPatch::updateArea() +{ + myTotalArea = 0.0; + + myTriangleCDF.resize (Elements.size()); + + for (Standard_Size aTrgIdx = 0; aTrgIdx < Elements.size(); ++aTrgIdx) + { + const BVH_Vec4i& aTriangle = Elements[aTrgIdx]; + + const BVH_Vec3d& aV1 = Vertices[aTriangle.x()]; + const BVH_Vec3d& aV2 = Vertices[aTriangle.y()]; + const BVH_Vec3d& aV3 = Vertices[aTriangle.z()]; + + myTriangleCDF[aTrgIdx] = triangleArea (aV1, aV2, aV3) + (aTrgIdx == 0 ? 0.0 : myTriangleCDF[aTrgIdx - 1]); + } + + myTotalArea = myTriangleCDF.back(); + Standard_Real aInvArea = 1.0 / myTotalArea; + + for (Standard_Size aTrgIdx = 0; aTrgIdx < myTriangleCDF.size(); ++aTrgIdx) + { + myTriangleCDF[aTrgIdx] *= aInvArea; + } +} diff --git a/src/BRepMesh/BRepMesh_TriangulatedPatch.hxx b/src/BRepMesh/BRepMesh_TriangulatedPatch.hxx new file mode 100644 index 0000000000..d864ebee48 --- /dev/null +++ b/src/BRepMesh/BRepMesh_TriangulatedPatch.hxx @@ -0,0 +1,201 @@ +// Created on: 2015-08-28 +// Created by: Danila ULYANOV +// Copyright (c) 2015 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _BRepMesh_TriangulatedPatch_HeaderFile +#define _BRepMesh_TriangulatedPatch_HeaderFile + +#include + +#include +#include + +#include + +class Handle (BRepMesh_TriangulatedPatch); + +//! Class representing a patch composed from triangles. +class BRepMesh_TriangulatedPatch : public BVH_Triangulation, public Standard_Transient +{ + +public: + + DEFINE_STANDARD_ALLOC + + typedef BRepMesh_EdgeSet EdgeSet; + +public: + + //! Creates new triangulated patch. + Standard_EXPORT BRepMesh_TriangulatedPatch(TopoDS_Face theFace); + + //! Updates edge and triangle data. + Standard_EXPORT void PreprocessData(); + + //! Flips the normals orientation of the patch. + Standard_EXPORT void Flip(); + + //! Flips only visibility coefficients of the patch. + Standard_EXPORT void FlipVisibility(); + + //! Returns triangle Id stored in 4th component of triangle record. + Standard_Integer TriangleId (Standard_Integer theIndex) const + { + return Elements[theIndex].w(); + } + + //! Returns set of patch boundary edges. + BRepMesh_EdgeSet& BoundaryEdgeSet() + { + return myBoundaryEdgeSet; + } + + //! Returns set of patch boundary edges. + const BRepMesh_EdgeSet& BoundaryEdgeSet() const + { + return myBoundaryEdgeSet; + } + + //! Returns boundary edge by index. + BRepMesh_OrientedEdge BoundaryEdge (Standard_Integer theIndex) const + { + BVH_Vec4i anEdge = myBoundaryEdgeSet.Elements[theIndex]; + return BRepMesh_OrientedEdge (anEdge.w() > 0.0 ? anEdge.x() : anEdge.y(), + anEdge.w() > 0.0 ? anEdge.y() : anEdge.x()); + } + + //! Returns edge vertices. + void GetEdgeVertices (const BRepMesh_OrientedEdge& theEdge, + BVH_Vec3d& theVertex1, + BVH_Vec3d& theVertex2) const + { + theVertex1 = Vertices[theEdge.FirstNode()]; + theVertex2 = Vertices[theEdge.LastNode()]; + } + + //! Appends another patch geometry. + Standard_EXPORT void Append (const Handle (BRepMesh_TriangulatedPatch)& thePatch); + + //! Returns precomputed CDF for triangles. + const std::vector& TriangleCDF() const + { + return myTriangleCDF; + } + + //! Returns total area of the patch. + Standard_Real TotalArea() const + { + return myTotalArea; + } + + //! Sets total area of the patch. + void SetTotalArea (Standard_Real theVal) + { + myTotalArea = theVal; + } + + //! Returns front visibility of the patch. + Standard_Real FrontVisibility() const + { + return myFrontVisibility; + } + + //! Sets front visibility of the patch. + void SetFrontVisibility (Standard_Real theVal) + { + myFrontVisibility = theVal; + } + + //! Returns back visibility of the patch. + Standard_Real BackVisibility() const + { + return myBackVisibility; + } + + //! Sets back visibility of the patch. + void SetBackVisibility (Standard_Real theVal) + { + myBackVisibility = theVal; + } + + //! Returns front distance for the patch. + Standard_Real FrontDistance() const + { + return myFrontDistance; + } + + //! Sets front distance for the patch. + void SetFrontDistance (Standard_Real theVal) + { + myFrontDistance = theVal; + } + + //! Returns back distance for the patch. + Standard_Real BackDistance() const + { + return myBackDistance; + } + + //! Sets back distance for the patch. + void SetBackDistance (Standard_Real theVal) + { + myBackDistance = theVal; + } + + //! Returns original topological face. + TopoDS_Face Face() + { + return myFace; + } + +protected: + + //! Extracts edges from triangulation. Fills edge -> faces map. Extracts boundary edges. Fills edge set. + void computeBoundaryEdges(); + + //! Computes areas of triangles. Computes total patch area. Computes sampling probabilities. + void updateArea(); + +private: + + //! Set of boundary edges. + EdgeSet myBoundaryEdgeSet; + + //! Total triangle area. + Standard_Real myTotalArea; + + //! Total front visibility. + Standard_Real myFrontVisibility; + + //! Total back visibility. + Standard_Real myBackVisibility; + + //! Accumulated front distance of the patch. + Standard_Real myFrontDistance; + + //! Accumulated back distance of the patch. + Standard_Real myBackDistance; + + //! Cumulative distribution function. + std::vector myTriangleCDF; + + //! Reference to original topological face. + TopoDS_Face myFace; + + DEFINE_STANDARD_RTTI (BRepMesh_TriangulatedPatch) +}; + +DEFINE_STANDARD_HANDLE(BRepMesh_TriangulatedPatch, Standard_Transient) + +#endif diff --git a/src/BRepMesh/FILES b/src/BRepMesh/FILES index a0476779a0..786051bb26 100755 --- a/src/BRepMesh/FILES +++ b/src/BRepMesh/FILES @@ -54,4 +54,14 @@ BRepMesh_IncrementalMesh.cxx BRepMesh_GeomTool.hxx BRepMesh_GeomTool.cxx BRepMesh_PairOfPolygon.hxx +BRepMesh_RestoreOrientationTool.cxx +BRepMesh_RestoreOrientationTool.hxx +BRepMesh_TriangulatedPatch.cxx +BRepMesh_TriangulatedPatch.hxx +BRepMesh_EdgeSet.hxx +BRepMesh_EdgeSet.lxx +BRepMesh_Block.hxx +BRepMesh_MinStCut.cxx +BRepMesh_MinStCut.hxx +BRepMesh_LineBoxDistance.hxx EXTERNLIB diff --git a/src/MeshTest/MeshTest.cxx b/src/MeshTest/MeshTest.cxx index 103ffc6988..b54a8f795d 100644 --- a/src/MeshTest/MeshTest.cxx +++ b/src/MeshTest/MeshTest.cxx @@ -14,87 +14,77 @@ // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. -#include - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include #include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include #include #include -#include -#include #include +#include #include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include -#include -#include #include +#include #include -#include - +#include +#include #include -#include #include -#include - -#include +#include +#include +#include +#include +#include #include -#include +#include +#include +#include +#include +#include +#include #include +#include +#include +#include +#include #include +#include #include +#include +#include +#include #include -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - - -//epa Memory leaks test -#include +#include +#include +#include +#include +#include +#include #include -#include -#include - -//OAN: for triepoints -#include -#include +#include #include +#include +#include +//epa Memory leaks test +//OAN: for triepoints #ifdef WNT Standard_IMPORT Draw_Viewer dout; #endif @@ -1615,6 +1605,45 @@ Standard_Integer correctnormals (Draw_Interpretor& theDI, return 0; } +//======================================================================= +//function : RestoreOrientation +//purpose : Makes normals of a shape consistently oriented +//======================================================================= +Standard_Integer RestoreOrientation (Draw_Interpretor& /*theDI*/, + Standard_Integer theNArg, + const char** theArgVal) +{ + bool aVisibilityOnly = false; + + if (theNArg < 2) + { + std::cout << "restoreorient shape [-visibility]\n"; + std::cout << "Please specify the shape\n"; + } + + if (theNArg > 2) + { + if (!strcmp (theArgVal[2], "-visibility")) + { + aVisibilityOnly = true; + } + } + + TopoDS_Shape aShape = DBRep::Get (theArgVal[1]); + + Handle(BRepMesh_RestoreOrientationTool) aTool = new BRepMesh_RestoreOrientationTool (aShape, aVisibilityOnly); + aTool->Perform(); + + if (aTool->IsDone()) + { + std::string aShapeName = std::string (theArgVal[1]) + "_restored"; + DBRep::Set (aShapeName.c_str(), aTool->Shape()); + std::cout << "Resulting shape: " << aShapeName << std::endl; + } + + return 0; +} + //======================================================================= void MeshTest::Commands(Draw_Interpretor& theCommands) //======================================================================= @@ -1651,6 +1680,8 @@ void MeshTest::Commands(Draw_Interpretor& theCommands) theCommands.Add("correctnormals", "correctnormals shape",__FILE__, correctnormals, g); + theCommands.Add("restoreorient", "restoreorient shape [-visibility]",__FILE__, RestoreOrientation, g); + #if 0 theCommands.Add("extrema","extrema ",__FILE__, extrema, g); theCommands.Add("vb","vb ",__FILE__, vb, g);