Salome HOME
Intersector: optim for SegSeg and ArcCSeg when an extremity is shared.
[tools/medcoupling.git] / src / INTERP_KERNEL / Geometric2D / InterpKernelGeo2DEdge.cxx
index 16625228794d106c2a78d875e0e932208a6d1dee..5d815d9478c70e8d24fc549426112ea3e56eabcd 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -28,6 +28,9 @@
 
 #include <algorithm>
 
+#include <functional>
+
+
 using namespace INTERP_KERNEL;
 
 MergePoints::MergePoints():_ass1Start1(0),_ass1End1(0),_ass1Start2(0),_ass1End2(0),
@@ -288,7 +291,7 @@ void IntersectElement::performMerging(MergePoints& commonNode) const
 }
 
 /*!
- * This methode is const because 'node' is supposed to be equal geomitrically to _node.
+ * This method is const because 'node' is supposed to be equal geometrically to _node.
  */
 void IntersectElement::setNode(Node *node) const
 {
@@ -337,7 +340,7 @@ bool IntersectElement::isIncludedByBoth() const
   return _e1.isIn(_chararct_val_for_e1) && _e2.isIn(_chararct_val_for_e2);
 }
 
-bool EdgeIntersector::intersect(const Bounds *whereToFind, std::vector<Node *>& newNodes, bool& order, MergePoints& commonNode)
+bool EdgeIntersector::intersect(std::vector<Node *>& newNodes, bool& order, MergePoints& commonNode)
 {
   std::list< IntersectElement > listOfIntesc=getIntersectionsCharacteristicVal();
   std::list< IntersectElement >::iterator iter;
@@ -388,6 +391,37 @@ bool EdgeIntersector::intersect(const Bounds *whereToFind, std::vector<Node *>&
   return true;
 }
 
+/*! If the 2 edges share one extremity, we can optimize since we already know where is the intersection.
+ *  In the case of ArcCSegIntersector, this also helps avoid degenerated cases.
+ */
+void EdgeIntersector::identifyEarlyIntersection(bool& i1S2S, bool& i1E2S, bool& i1S2E, bool& i1E2E)
+{
+  i1S2S = _e1.getStartNode() == _e2.getStartNode();
+  i1E2S = _e1.getEndNode() == _e2.getStartNode();
+  i1S2E = _e1.getStartNode() == _e2.getEndNode();
+  i1E2E = _e1.getEndNode() == _e2.getEndNode();
+  if (i1S2S || i1E2S || i1S2E || i1E2E)
+    {
+      Node * node;
+      bool i_1S(false),i_1E(false),i_2S(false),i_2E(false);
+      if (i1S2S || i1E2S)   // Common node is e2 start
+        {
+          node = _e2.getStartNode();
+          i_1S = i1S2S;        i_2S = true;
+          i_1E = i1E2S;        i_2E = false;
+        }
+      else                  // Common node is e2 end
+        {
+          node = _e2.getEndNode();
+          i_1S = i1S2E;        i_2S = false;
+          i_1E = i1E2E;        i_2E = true;
+        }
+      node->incrRef();
+      _earlyInter = new IntersectElement(_e1.getCharactValue(*node), _e2.getCharactValue(*node),
+          i_1S,i_1E,i_2S,i_2E,node,_e1,_e2,keepOrder());
+    }
+}
+
 /*!
  * Locates 'node' regarding edge this->_e1. If node is located close to (with distant lt epsilon) start or end point of _e1,
  * 'node' takes its place. In this case 'obvious' is set to true and 'commonNode' stores information of merge point and finally 'where' is set.
@@ -628,7 +662,7 @@ bool Edge::intersectWith(const Edge *other, MergePoints& commonNode,
   delete merge;
   merge=0;
   EdgeIntersector *intersector=BuildIntersectorWith(this,other);
-  ret=Intersect(this,other,intersector,merge,commonNode,outVal1,outVal2);
+  ret=Intersect(this,other,intersector,commonNode,outVal1,outVal2);
   delete intersector;
   return ret;
 }
@@ -671,7 +705,7 @@ void Edge::Interpolate1DLin(const std::vector<double>& distrib1, const std::vect
                   ComposedEdge *f2=new ComposedEdge;
                   SegSegIntersector inters(*e1,*e2);
                   bool b1,b2;
-                  inters.areOverlappedOrOnlyColinears(0,b1,b2);
+                  inters.areOverlappedOrOnlyColinears(b1,b2);
                   if(IntersectOverlapped(e1,e2,&inters,commonNode,*f1,*f2))
                     {
                       result[i][j]=f1->getCommonLengthWith(*f2)/e1->getCurveLength();
@@ -731,22 +765,22 @@ void Edge::getMiddleOfPointsOriented(const double *p1, const double *p2, double
   return getMiddleOfPoints(p1, p2, mid);
 }
 
-bool Edge::Intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, const Bounds *whereToFind, MergePoints& commonNode,
+bool Edge::Intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, MergePoints& commonNode,
                      ComposedEdge& outValForF1, ComposedEdge& outValForF2)
 {
   bool obviousNoIntersection;
   bool areOverlapped;
-  intersector->areOverlappedOrOnlyColinears(whereToFind,obviousNoIntersection,areOverlapped);
+  intersector->areOverlappedOrOnlyColinears(obviousNoIntersection,areOverlapped);
   if(areOverlapped)
     return IntersectOverlapped(f1,f2,intersector,commonNode,outValForF1,outValForF2);
   if(obviousNoIntersection)
     return false;
   std::vector<Node *> newNodes;
   bool order;
-  if(intersector->intersect(whereToFind,newNodes,order,commonNode))
+  if(intersector->intersect(newNodes,order,commonNode))
     {
       if(newNodes.empty())
-        throw Exception("Internal error occured - error in intersector implementation!");// This case should never happen
+        throw Exception("Internal error occurred - error in intersector implementation!");// This case should never happen
       std::vector<Node *>::iterator iter=newNodes.begin();
       std::vector<Node *>::reverse_iterator iterR=newNodes.rbegin();
       f1->addSubEdgeInVector(f1->getStartNode(),*iter,outValForF1);
@@ -768,7 +802,7 @@ bool Edge::Intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersecto
         }
       return true;
     }
-  else//no intersection inside whereToFind
+  else
     return false;
 }
 
@@ -923,14 +957,25 @@ bool Edge::SplitOverlappedEdges(const Edge *e1, const Edge *e2, Node *nS, Node *
   }
 }
 
-bool Edge::isEqual(const Edge& other) const
+void Edge::dumpToCout(const std::map<INTERP_KERNEL::Node *,int>& mapp, int index) const
 {
-  return _start->isEqual(*other._start) && _end->isEqual(*other._end);
+  auto sI(mapp.find(getStartNode())), eI(mapp.find(getEndNode()));
+  int start = (sI == mapp.end() ? -1 : sI->second), end = (eI == mapp.end() ? -1 : eI->second);
+  std::string locs;
+  switch (getLoc())
+  {
+    case FULL_IN_1: locs="FULL_IN_1"; break;
+    case FULL_ON_1: locs="FULL_ON_1"; break;
+    case FULL_OUT_1: locs="FULL_OUT_1"; break;
+    case FULL_UNKNOWN: locs="FULL_UNKNOWN"; break;
+    default: locs="oh my God! This is so wrong.";
+  }
+  std::cout << "Edge [" << index << "] : ("<<  std::hex << this << std::dec << ") -> (" << start << ", " << end << ")\t" << locs << std::endl;
 }
 
-inline bool eqpair(const std::pair<double,Node *>& p1, const std::pair<double,Node *>& p2)
+bool Edge::isEqual(const Edge& other) const
 {
-  return fabs(p1.first-p2.first)<QUADRATIC_PLANAR::_precision;
+  return _start->isEqual(*other._start) && _end->isEqual(*other._end);
 }
 
 /**
@@ -982,6 +1027,15 @@ bool Edge::sortSubNodesAbs(const double *coo, std::vector<int>& subNodes)
 void Edge::sortIdsAbs(const std::vector<INTERP_KERNEL::Node *>& addNodes, const std::map<INTERP_KERNEL::Node *, int>& mapp1,
                       const std::map<INTERP_KERNEL::Node *, int>& mapp2, std::vector<int>& edgesThis)
 {
+  int startId=(*mapp1.find(_start)).second;
+  int endId=(*mapp1.find(_end)).second;
+  if (! addNodes.size()) // quick way out, no new node to add.
+    {
+      edgesThis.push_back(startId);
+      edgesThis.push_back(endId);
+      return;
+    }
+
   Bounds b;
   b.prepareForAggregation();
   b.aggregate(getBounds());
@@ -998,23 +1052,10 @@ void Edge::sortIdsAbs(const std::vector<INTERP_KERNEL::Node *>& addNodes, const
   for(std::size_t i=0;i<sz;i++)
     an2[i]=std::pair<double,Node *>(getCharactValueBtw0And1(*addNodes[i]),addNodes[i]);
   std::sort(an2.begin(),an2.end());
-  int startId=(*mapp1.find(_start)).second;
-  int endId=(*mapp1.find(_end)).second;
   std::vector<int> tmpp;
-  std::vector< std::pair<double,Node *> >::const_iterator itend=std::unique(an2.begin(),an2.end(),eqpair);
-  for(std::vector< std::pair<double,Node *> >::const_iterator it=an2.begin();it!=itend;it++)
+  for(std::vector< std::pair<double,Node *> >::const_iterator it=an2.begin();it!=an2.end();it++)
     {
       int idd=(*mapp2.find((*it).second)).second;
-      if((*it).first<QUADRATIC_PLANAR::_precision)
-        {
-          startId=idd;
-          continue;
-        }
-      if((*it).first>1-QUADRATIC_PLANAR::_precision)
-        {
-          endId=idd;
-          continue;
-        }
       tmpp.push_back(idd);
     }
   std::vector<int> tmpp2(tmpp.size()+2);
@@ -1030,3 +1071,31 @@ void Edge::sortIdsAbs(const std::vector<INTERP_KERNEL::Node *>& addNodes, const
       edgesThis.push_back(tmpp2[i+1]);
     }
 }
+
+void Edge::fillGlobalInfoAbs(bool direction, const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
+                                std::vector<int>& edgesThis, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int> mapAddCoo) const
+{
+  int tmp[2];
+  _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp);
+  _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1);
+  if(direction)
+    {
+      edgesThis.push_back(tmp[0]);
+      edgesThis.push_back(tmp[1]);
+    }
+  else
+    {
+      edgesThis.push_back(tmp[1]);
+      edgesThis.push_back(tmp[0]);
+    }
+}
+
+void Edge::fillGlobalInfoAbs2(const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
+                              short skipStartOrEnd,
+                              std::vector<int>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int>& mapAddCoo) const
+{
+  if (skipStartOrEnd != -1) // see meaning in splitAbs()
+    _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
+  if (skipStartOrEnd != 1)
+  _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
+}