--- /dev/null
+// 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 <stdlib.h>
+
+template <class Type> 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 Type> 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; anItem<myFirstFree+myBLockSize-1; anItem++)
+ { anItem->NextFree = 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
+
--- /dev/null
+// 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 <BVH_PrimitiveSet.hxx>
+#include <BVH_Types.hxx>
+
+//! Class representing set of boundary edges.
+template<class T, int N>
+class BRepMesh_EdgeSet : public BVH_PrimitiveSet<T, N>
+{
+public:
+
+ typedef typename BVH::VectorType<T, N>::Type BVH_VecNt;
+ typedef typename BVH::ArrayType<T, N>::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<T, N>::Box;
+
+ //! Returns AABB of the given edge.
+ virtual BVH_Box<T, N> 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 <BRepMesh_EdgeSet.lxx>
+#endif // _BRepMesh_EdgeSet_Header
--- /dev/null
+// 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 <BRepMesh_LineBoxDistance.hxx>
+
+// =======================================================================
+// function : BRepMesh_EdgeSet
+// purpose :
+// =======================================================================
+template<class T, int N>
+BRepMesh_EdgeSet<T, N>::BRepMesh_EdgeSet (const BVH_ArrayNt& theTriangulationVertices)
+ : Vertices (theTriangulationVertices)
+{
+ //
+}
+
+// =======================================================================
+// function : ~BRepMesh_EdgeSet
+// purpose :
+// =======================================================================
+template<class T, int N>
+BRepMesh_EdgeSet<T, N>::~BRepMesh_EdgeSet()
+{
+ //
+}
+
+// =======================================================================
+// function : Size
+// purpose :
+// =======================================================================
+template<class T, int N>
+Standard_Integer BRepMesh_EdgeSet<T, N>::Size() const
+{
+ return BVH::Array<Standard_Integer, 4>::Size (Elements);
+}
+
+// =======================================================================
+// function : Box
+// purpose :
+// =======================================================================
+template<class T, int N>
+BVH_Box<T, N> BRepMesh_EdgeSet<T, N>::Box (const Standard_Integer theIndex) const
+{
+ const BVH_Vec4i& anIndex = BVH::Array<Standard_Integer, 4>::Value (Elements, theIndex);
+
+ const BVH_VecNt& aPoint0 = BVH::Array<T, N>::Value (Vertices, anIndex.x());
+ const BVH_VecNt& aPoint1 = BVH::Array<T, N>::Value (Vertices, anIndex.y());
+
+ BVH_VecNt aMinPoint = aPoint0;
+ BVH_VecNt aMaxPoint = aPoint0;
+
+ BVH::BoxMinMax<T, N>::CwiseMin (aMinPoint, aPoint1);
+ BVH::BoxMinMax<T, N>::CwiseMax (aMaxPoint, aPoint1);
+
+ return BVH_Box<T, N> (aMinPoint, aMaxPoint);
+}
+
+// =======================================================================
+// function : Center
+// purpose :
+// =======================================================================
+template<class T, int N>
+T BRepMesh_EdgeSet<T, N>::Center (const Standard_Integer theIndex,
+ const Standard_Integer theAxis) const
+{
+ const BVH_Vec4i& anIndex = BVH::Array<Standard_Integer, 4>::Value (Elements, theIndex);
+
+ const BVH_VecNt& aPoint0 = BVH::Array<T, N>::Value (Vertices, anIndex.x());
+ const BVH_VecNt& aPoint1 = BVH::Array<T, N>::Value (Vertices, anIndex.y());
+
+ return ( BVH::VecComp<T, N>::Get (aPoint0, theAxis) +
+ BVH::VecComp<T, N>::Get (aPoint1, theAxis) ) * static_cast<T> (1.0 / 2.0);
+}
+
+// =======================================================================
+// function : Swap
+// purpose :
+// =======================================================================
+template<class T, int N>
+void BRepMesh_EdgeSet<T, N>::Swap (const Standard_Integer theIndex1,
+ const Standard_Integer theIndex2)
+{
+ BVH_Vec4i& anIndices1 = BVH::Array<Standard_Integer, 4>::ChangeValue (Elements, theIndex1);
+ BVH_Vec4i& anIndices2 = BVH::Array<Standard_Integer, 4>::ChangeValue (Elements, theIndex2);
+
+ std::swap (anIndices1, anIndices2);
+}
+
+// =======================================================================
+// function : ClosestEdge
+// purpose :
+// =======================================================================
+template<class T, int N>
+Standard_Integer BRepMesh_EdgeSet<T, N>::ClosestEdge (const BVH_VecNt& theEdgeV1,
+ const BVH_VecNt& theEdgeV2,
+ T* theDist)
+{
+ Standard_Integer aStack[32];
+
+ const NCollection_Handle<BVH_Tree<T, N> >& aBVH = BVH();
+
+ if (aBVH.IsNull())
+ {
+ return -1;
+ }
+
+ Standard_Integer aNode = 0;
+
+ T aMinDist = std::numeric_limits<T>::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<class T, int N>
+Standard_Integer BRepMesh_EdgeSet<T, N>::CoherentEdge (const BVH_VecNt& theEdgeV1,
+ const BVH_VecNt& theEdgeV2,
+ T* theDist)
+{
+ Standard_Integer aStack[32];
+
+ const NCollection_Handle<BVH_Tree<T, N> >& aBVH = BVH();
+
+ if (aBVH.IsNull())
+ {
+ return -1;
+ }
+
+ Standard_Integer aNode = 0;
+
+ T aMinDist = std::numeric_limits<T>::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<class T, int N>
+T BRepMesh_EdgeSet<T, N>::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<class T, int N>
+T BRepMesh_EdgeSet<T, N>::EdgeToBoxDistance (const BVH_VecNt& theEV1,
+ const BVH_VecNt& theEV2,
+ const BVH_VecNt& theMinPoint,
+ const BVH_VecNt& theMaxPoint)
+{
+ T anSqrDist = std::numeric_limits<T>::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<T, N> aDist;
+ LineToBoxDistance<Standard_Real, 3>::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);
+}
--- /dev/null
+// 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 <BVH_Types.hxx>
+
+template <class T, int N>
+class LineToBoxDistance
+{
+public:
+
+ typedef typename BVH::VectorType<T, N>::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<class T, int N>
+typename LineToBoxDistance<T, N>::Result
+LineToBoxDistance<T, N>::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<class T, int N>
+void LineToBoxDistance<T, N>::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<class T, int N>
+void LineToBoxDistance<T, N>::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<class T, int N>
+void LineToBoxDistance<T, N>::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<class T, int N>
+void LineToBoxDistance<T, N>::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<class T, int N>
+void LineToBoxDistance<T, N>::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<class T, int N>
+void LineToBoxDistance<T, N>::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
--- /dev/null
+// 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 <stdio.h>
+#include <BRepMesh_MinStCut.hxx>
+
+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]; aNode<aNodeBlock->Current; 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]; aNode<aNodeBlock->Current; 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; anArcFor<aForwardLAst; anArcFor++)
+ {
+ Node* to = NEIGHBOR_NODE (aNode, anArcFor->Shift);
+ 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]; aNode<aNodeBlock->Current; 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; anArc0For<anArc0ForLast; anArc0For++)
+ if (anArc0For->ReverseResidualCap)
+ {
+ 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 (aDist<INFINITE_D)
+ {
+ if (aDist<aDistMin)
+ {
+ anArc0Min = anArc0For;
+ aDistMin = aDist;
+ }
+
+ for (aNode=NEIGHBOR_NODE (theNode, anArc0For->Shift); 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<anArc0RevLast; anArc0Rev++)
+ {
+ anArc0For = 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 (aDist<INFINITE_D)
+ {
+ if (aDist<aDistMin)
+ {
+ anArc0Min = MAKE_ODD (anArc0For);
+ aDistMin = aDist;
+ }
+
+ for (aNode=NEIGHBOR_NODE_REV (theNode,anArc0For->Shift); 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<anArc0ForLast; anArc0For++)
+ {
+ aNode = NEIGHBOR_NODE (theNode, 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<anArc0RevLast; anArc0Rev++)
+ {
+ anArc0For = 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; anArc0For<anArc0ForLast; anArc0For++)
+ if (anArc0For->ResidualCap)
+ {
+ 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 (aDist<INFINITE_D)
+ {
+ if (aDist<aDistMin)
+ {
+ anArc0Min = anArc0For;
+ aDistMin = aDist;
+ }
+
+ for (aNode=NEIGHBOR_NODE (theNode, anArc0For->Shift); 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<anArc0RevLast; anArc0Rev++)
+ {
+ anArc0For = 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 (aDist<INFINITE_D)
+ {
+ if (aDist<aDistMin)
+ {
+ anArc0Min = MAKE_ODD (anArc0For);
+ aDistMin = aDist;
+ }
+
+ for (aNode=NEIGHBOR_NODE_REV (theNode,anArc0For->Shift); 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<anArc0ForLast; anArc0For++)
+ {
+ aNode = NEIGHBOR_NODE (theNode, 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<anArc0RevLast; anArc0Rev++)
+ {
+ anArc0For = 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> (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; anArcFor<anArcForLast; anArcFor++)
+ if (anArcFor->ResidualCap)
+ {
+ 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<anArcRevLAst; anArcRev++)
+ {
+ anArcFor = 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; anArcFor<anArcForLast; anArcFor++)
+ if (anArcFor->ReverseResidualCap)
+ {
+ 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<anArcRevLAst; anArcRev++)
+ {
+ anArcFor = 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;
+}
--- /dev/null
+// 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 <assert.h>
+#include <BRepMesh_Block.hxx>
+
+#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<NodePtr>* 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
--- /dev/null
+// 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 <BRepMesh_RestoreOrientationTool.hxx>
+
+#include <OSD_Timer.hxx>
+#include <TCollection_AsciiString.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopoDS_Compound.hxx>
+
+#include <BRep_Builder.hxx>
+#include <BRep_Tool.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <Poly_Triangulation.hxx>
+
+#include <NCollection_Map.hxx>
+#include <NCollection_DataMap.hxx>
+
+#include <Standard_Integer.hxx>
+
+#include <Precision.hxx>
+
+#include <math_BullardGenerator.hxx>
+
+#include <set>
+#include <map>
+
+#include <BRepMesh_MinStCut.hxx>
+
+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<Standard_Real>::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<Standard_Real>& theCDF)
+{
+ return static_cast<Standard_Integer> (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<Standard_ShortReal> (theValue);
+ }
+
+ //! Returns coherence value.
+ Standard_Real getCoherence (size_t theI, size_t theJ)
+ {
+ return static_cast<Standard_Real> (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<Standard_ShortReal> 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<const TopoDS_Face&> (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<Standard_Real, 3>& 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<bool> (isMirrored ^ (theFace.Orientation() == TopAbs_REVERSED));
+
+ BRepAdaptor_Surface aFaceAdaptor (theFace, Standard_False);
+
+ const Standard_Integer aVertOffset =
+ static_cast<Standard_Integer> (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<Standard_Real, 3>& 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<Standard_Real, 3>& theTriangulation,
+ Standard_Real& theDistance)
+{
+ Standard_Integer aStack[32];
+
+ const NCollection_Handle<BVH_Tree<Standard_Real, 3> >& aBVH = theTriangulation.BVH();
+
+ if (aBVH.IsNull())
+ {
+ return INVALID_HIT;
+ }
+
+ Standard_Integer aNode = 0;
+
+ Standard_Real aMinDist = std::numeric_limits<Standard_Real>::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<Standard_Real>::epsilon();
+
+#ifdef MEASURE_PERFORMANCE
+ OSD_Timer aTimer;
+ aTimer.Start();
+#endif
+
+ std::auto_ptr<CoherencyTable> aTable;
+
+ // Flipped flags for all patches.
+ std::vector<char> aFlipped (myPatches.size(), 0);
+ std::vector<char> aFlippedC (myPatches.size(), 0);
+
+ Standard_Real aMaxCoherence = -std::numeric_limits<Standard_Real>::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<BRepMesh_SuperPatch> 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<Standard_Integer> (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<Standard_Integer> 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<BRepMesh_SuperPatchPair> 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<Standard_Integer> (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<Standard_Integer> (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<Energy::Var> 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;
+}
--- /dev/null
+// 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 <Standard.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS_Face.hxx>
+
+#include <BRepMesh_TriangulatedPatch.hxx>
+
+//! 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<Standard_Integer>::iterator aNextPos =
+ std::copy (theFirst.PatchIds.begin(), theFirst.PatchIds.end(), PatchIds.begin());
+
+ std::copy (theSecond.PatchIds.begin(), theSecond.PatchIds.end(), aNextPos);
+ }
+
+ std::vector<Standard_Integer> 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<TopoDS_Face> 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<Standard_Real, 3>& 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<Standard_Real, 3>& theTriangulation,
+ Standard_Real& theDistance);
+
+ //! Extracts a triangulation from the face.
+ static Standard_Integer extractTriangulation (const TopoDS_Face& theShape,
+ BVH_Triangulation<Standard_Real, 3>& 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<Handle (BRepMesh_TriangulatedPatch)> myPatches;
+
+ //! Vertex and normal data of triangulated shape.
+ BVH_Triangulation<Standard_Real, 3> myTriangulation;
+
+ std::vector<BRepMesh_SuperPatchPair> myPatchQueue;
+
+ //! Maps Id intervals to topological faces.
+ std::vector<Standard_Integer> myFaceIdIntervals;
+
+ //! Use only visibility metric?
+ bool myVisibilityOnly;
+
+ DEFINE_STANDARD_RTTI (BRepMesh_RestoreOrientationTool)
+};
+
+DEFINE_STANDARD_HANDLE(BRepMesh_RestoreOrientationTool, Standard_Transient)
+
+
+#endif
--- /dev/null
+// 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 <BRepMesh_TriangulatedPatch.hxx>
+#include <BVH_LinearBuilder.hxx>
+
+#include <set>
+#include <map>
+
+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<int, int>& theEdge1,
+ const std::pair<int, int>& 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<int, int>& theEdge1,
+ const std::pair<int, int>& 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<int> SetOfInteger;
+
+ //! Type of edge-triangles connectivity map.
+ typedef std::map<std::pair<int, int>, SetOfInteger, SymmetricalPairComparator> EdgeFaceMap;
+}
+
+// =======================================================================
+// function : BRepMesh_TriangulatedPatch
+// purpose :
+// =======================================================================
+BRepMesh_TriangulatedPatch::BRepMesh_TriangulatedPatch (TopoDS_Face theFace)
+ : BVH_Triangulation<Standard_Real, 3>(),
+ myBoundaryEdgeSet (Vertices),
+ myTotalArea (0.0),
+ myFrontVisibility (0.0),
+ myBackVisibility (0.0),
+ myFrontDistance (0.0),
+ myBackDistance (0.0),
+ myFace (theFace)
+{
+ myBuilder = new BVH_LinearBuilder<Standard_Real, 3> (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<Standard_Integer> (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<Standard_Integer, Standard_Integer> anEdge (aTriangle[anEdgeStart], aTriangle[edgeEnd]);
+
+ anEdgeFaceMap[anEdge].insert (static_cast<Standard_Integer> (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<Standard_Integer, Standard_Integer>& 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;
+ }
+}
--- /dev/null
+// 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 <BVH_Triangulation.hxx>
+
+#include <BRepMesh_OrientedEdge.hxx>
+#include <BRepMesh_EdgeSet.hxx>
+
+#include <TopoDS_Face.hxx>
+
+class Handle (BRepMesh_TriangulatedPatch);
+
+//! Class representing a patch composed from triangles.
+class BRepMesh_TriangulatedPatch : public BVH_Triangulation<Standard_Real, 3>, public Standard_Transient
+{
+
+public:
+
+ DEFINE_STANDARD_ALLOC
+
+ typedef BRepMesh_EdgeSet<Standard_Real, 3> 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<Standard_Real, 3>& BoundaryEdgeSet()
+ {
+ return myBoundaryEdgeSet;
+ }
+
+ //! Returns set of patch boundary edges.
+ const BRepMesh_EdgeSet<Standard_Real, 3>& 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<Standard_Real>& 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<Standard_Real> myTriangleCDF;
+
+ //! Reference to original topological face.
+ TopoDS_Face myFace;
+
+ DEFINE_STANDARD_RTTI (BRepMesh_TriangulatedPatch)
+};
+
+DEFINE_STANDARD_HANDLE(BRepMesh_TriangulatedPatch, Standard_Transient)
+
+#endif
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
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.
-#include <Standard_Stream.hxx>
-
-#include <stdio.h>
-#include <MeshTest.ixx>
-
-#include <MeshTest_DrawableMesh.hxx>
-#include <TopAbs_ShapeEnum.hxx>
-#include <TopoDS.hxx>
-#include <TopoDS_Edge.hxx>
-#include <TopoDS_Face.hxx>
-#include <TopoDS_Shape.hxx>
-#include <TopoDS_Compound.hxx>
-#include <TopExp_Explorer.hxx>
-#include <TopTools_ListIteratorOfListOfShape.hxx>
-#include <DBRep.hxx>
-#include <BRepTest.hxx>
-#include <GeometryTest.hxx>
-#include <BRep_Tool.hxx>
+#include <AppCont_ContMatrices.hxx>
+#include <Bnd_Box.hxx>
#include <BRep_Builder.hxx>
-#include <Draw_MarkerShape.hxx>
-#include <Draw_Appli.hxx>
-#include <Draw.hxx>
-#include <DrawTrSurf.hxx>
-#include <BRepMesh_Triangle.hxx>
+#include <BRep_Tool.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepBndLib.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakePolygon.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
+#include <BRepLib.hxx>
#include <BRepMesh_DataStructureOfDelaun.hxx>
#include <BRepMesh_Delaun.hxx>
-#include <BRepMesh_FastDiscret.hxx>
-#include <BRepMesh_Vertex.hxx>
#include <BRepMesh_Edge.hxx>
+#include <BRepMesh_FastDiscret.hxx>
#include <BRepMesh_IncrementalMesh.hxx>
-#include <TColStd_ListIteratorOfListOfInteger.hxx>
-#include <TColStd_MapIteratorOfMapOfInteger.hxx>
-#include <Bnd_Box.hxx>
-#include <Precision.hxx>
+#include <BRepMesh_Triangle.hxx>
+#include <BRepMesh_Vertex.hxx>
+#include <BRepTest.hxx>
+#include <BRepTools.hxx>
+#include <CSLib.hxx>
+#include <CSLib_DerivativeStatus.hxx>
+#include <DBRep.hxx>
+#include <Draw.hxx>
+#include <Draw_Appli.hxx>
#include <Draw_Interpretor.hxx>
-#include <Geom_Plane.hxx>
-#include <Geom_Surface.hxx>
#include <Draw_Marker3D.hxx>
+#include <Draw_MarkerShape.hxx>
#include <Draw_Segment2D.hxx>
-#include <TCollection_AsciiString.hxx>
-
+#include <DrawTrSurf.hxx>
+#include <Extrema_LocateExtPC.hxx>
#include <GCPnts_UniformAbscissa.hxx>
-#include <GeomAdaptor_Curve.hxx>
#include <Geom_Curve.hxx>
-#include <Extrema_LocateExtPC.hxx>
-
-#include <TopLoc_Location.hxx>
+#include <Geom_Plane.hxx>
+#include <Geom_Surface.hxx>
+#include <GeomAdaptor_Curve.hxx>
+#include <GeometryTest.hxx>
+#include <gp_Pln.hxx>
#include <gp_Trsf.hxx>
-#include <Poly_Triangulation.hxx>
+#include <math.hxx>
+#include <math_Matrix.hxx>
+#include <math_Vector.hxx>
+#include <MeshTest.hxx>
+#include <MeshTest_DrawableMesh.hxx>
+#include <PLib.hxx>
#include <Poly_Connect.hxx>
+#include <Poly_PolygonOnTriangulation.hxx>
+#include <Poly_Triangulation.hxx>
+#include <Precision.hxx>
+#include <Standard_Stream.hxx>
#include <TColgp_Array1OfPnt2d.hxx>
+#include <TCollection_AsciiString.hxx>
#include <TColStd_HArray1OfInteger.hxx>
+#include <TColStd_ListIteratorOfListOfInteger.hxx>
+#include <TColStd_MapIteratorOfMapOfInteger.hxx>
+#include <TopAbs_ShapeEnum.hxx>
#include <TopExp_Explorer.hxx>
-#include <gp_Pln.hxx>
-
-#include <PLib.hxx>
-#include <AppCont_ContMatrices.hxx>
-#include <math_Vector.hxx>
-#include <math_Matrix.hxx>
-#include <math.hxx>
-
-#include <CSLib_DerivativeStatus.hxx>
-#include <CSLib.hxx>
-#include <BRepAdaptor_Surface.hxx>
-#include <Bnd_Box.hxx>
-#include <BRepBndLib.hxx>
-#include <BRepLib.hxx>
-
-
-//epa Memory leaks test
-#include <BRepBuilderAPI_MakePolygon.hxx>
+#include <TopLoc_Location.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Compound.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Shape.hxx>
#include <TopoDS_Wire.hxx>
-#include <BRepBuilderAPI_MakeFace.hxx>
-#include <BRepTools.hxx>
-
-//OAN: for triepoints
-#include <BRepBuilderAPI_MakeVertex.hxx>
-#include <Poly_PolygonOnTriangulation.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <TopTools_MapIteratorOfMapOfShape.hxx>
+#include <BRepMesh_RestoreOrientationTool.hxx>
+#include <stdio.h>
+//epa Memory leaks test
+//OAN: for triepoints
#ifdef WNT
Standard_IMPORT Draw_Viewer dout;
#endif
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)
//=======================================================================
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);
projects / occt-copy.git / commitdiff