]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx
Salome HOME
Ensure that ShapeRecognition links against cblas
[tools/medcoupling.git] / src / INTERP_KERNEL / Geometric2D / InterpKernelGeo2DEdge.cxx
1 // Copyright (C) 2007-2024  CEA, EDF
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // Author : Anthony Geay (CEA/DEN)
20
21 #include "InterpKernelGeo2DEdge.hxx"
22 #include "InterpKernelGeo2DEdgeLin.hxx"
23 #include "InterpKernelGeo2DEdgeInfLin.hxx"
24 //#include "EdgeParabol.hxx"
25 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
26 #include "InterpKernelGeo2DComposedEdge.hxx"
27 #include "InterpKernelException.hxx"
28
29 #include <algorithm>
30
31 #include <functional>
32
33
34 using namespace INTERP_KERNEL;
35
36 MergePoints::MergePoints():_ass1Start1(0),_ass1End1(0),_ass1Start2(0),_ass1End2(0),
37     _ass2Start1(0),_ass2End1(0),_ass2Start2(0),_ass2End2(0)
38 {
39 }
40
41 void MergePoints::start1Replaced()
42 {
43   unsigned nbOfAsso=getNumberOfAssociations();
44   if(nbOfAsso==0)
45     _ass1Start1=1;
46   else
47     _ass2Start1=1;
48 }
49
50 void MergePoints::end1Replaced()
51 {
52   unsigned nbOfAsso=getNumberOfAssociations();
53   if(nbOfAsso==0)
54     _ass1End1=1;
55   else
56     _ass2End1=1;
57 }
58
59 void MergePoints::start1OnStart2()
60 {
61   unsigned nbOfAsso=getNumberOfAssociations();
62   if(nbOfAsso==0)
63     {
64       _ass1Start1=1;
65       _ass1Start2=1;
66     }
67   else
68     {
69       _ass2Start1=1;
70       _ass2Start2=1;
71     }
72 }
73
74 void MergePoints::start1OnEnd2()
75 {
76   unsigned nbOfAsso=getNumberOfAssociations();
77   if(nbOfAsso==0)
78     {
79       _ass1Start1=1;
80       _ass1End2=1;
81     }
82   else
83     {
84       _ass2Start1=1;
85       _ass2End2=1;
86     }
87 }
88
89 void MergePoints::end1OnStart2()
90 {
91   unsigned nbOfAsso=getNumberOfAssociations();
92   if(nbOfAsso==0)
93     {
94       _ass1End1=1;
95       _ass1Start2=1;
96     }
97   else
98     {
99       _ass2End1=1;
100       _ass2Start2=1;
101     }
102 }
103
104 void MergePoints::end1OnEnd2()
105 {
106   unsigned nbOfAsso=getNumberOfAssociations();
107   if(nbOfAsso==0)
108     {
109       _ass1End1=1;
110       _ass1End2=1;
111     }
112   else
113     {
114       _ass2End1=1;
115       _ass2End2=1;
116     }
117 }
118
119 bool MergePoints::isStart1(unsigned rk) const
120 {
121   if(rk==0)
122     return _ass1Start1;
123   else
124     return _ass2Start1;
125 }
126
127 bool MergePoints::isEnd1(unsigned rk) const
128 {
129   if(rk==0)
130     return _ass1End1;
131   else
132     return _ass2End1;
133 }
134
135 bool MergePoints::isStart2(unsigned rk) const
136 {
137   if(rk==0)
138     return _ass1Start2;
139   else
140     return _ass2Start2;
141 }
142
143 bool MergePoints::isEnd2(unsigned rk) const
144 {
145   if(rk==0)
146     return _ass1End2;
147   else
148     return _ass2End2;
149 }
150
151 void MergePoints::clear()
152 {
153   _ass1Start1=0;_ass1End1=0;_ass1Start2=0;_ass1End2=0;
154   _ass2Start1=0;_ass2End1=0;_ass2Start2=0;_ass2End2=0;
155 }
156
157 unsigned MergePoints::getNumberOfAssociations() const
158 {
159   unsigned ret=0;
160   unsigned subTot=_ass1Start1+_ass1End1+_ass1Start2+_ass1End2;
161   if(subTot!=0)
162     ret++;
163   subTot=_ass2Start1+_ass2End1+_ass2Start2+_ass2End2;
164   if(subTot!=0)
165     ret++;
166   return ret;
167 }
168
169 void MergePoints::PushInMap(mcIdType key, mcIdType value, std::map<mcIdType,mcIdType>& mergedNodes)
170 {
171   if(key!=-1 && value!=-1)
172     mergedNodes[key]=value;
173 }
174
175 void MergePoints::updateMergedNodes(mcIdType e1Start, mcIdType e1End, mcIdType e2Start, mcIdType e2End, std::map<mcIdType,mcIdType>& mergedNodes)
176 {
177   unsigned subTot(_ass1Start1+_ass1End1+_ass1Start2+_ass1End2);
178   if(subTot!=0)
179     {
180       if(_ass1Start1 && _ass1Start2)
181         PushInMap(e2Start,e1Start,mergedNodes);
182       if(_ass1Start1 && _ass1End2)
183         PushInMap(e2End,e1Start,mergedNodes);
184       if(_ass1End1 && _ass1Start2)
185         PushInMap(e2Start,e1End,mergedNodes);
186       if(_ass1End1 && _ass1End2)
187         PushInMap(e2End,e1End,mergedNodes);
188     }
189   subTot=_ass2Start1+_ass2End1+_ass2Start2+_ass2End2;
190   if(subTot!=0)
191     {
192       if(_ass2Start1 && _ass2Start2)
193         PushInMap(e2Start,e1Start,mergedNodes);
194       if(_ass2Start1 && _ass2End2)
195         PushInMap(e2End,e1Start,mergedNodes);
196       if(_ass2End1 && _ass2Start2)
197         PushInMap(e2Start,e1End,mergedNodes);
198       if(_ass2End1 && _ass2End2)
199         PushInMap(e2End,e1End,mergedNodes);
200     }
201 }
202
203 IntersectElement::IntersectElement(double val1, double val2, bool start1, bool end1, bool start2, bool end2, Node *node
204                                    , const Edge& e1, const Edge& e2, bool keepOrder):_1S(keepOrder?start1:start2),
205                                        _1E(keepOrder?end1:end2),
206                                        _2S(keepOrder?start2:start1),
207                                        _2E(keepOrder?end2:end1),
208                                        _chararct_val_for_e1(keepOrder?val1:val2),
209                                        _chararct_val_for_e2(keepOrder?val2:val1),
210                                        _node(node),_loc_of_node(node->getLoc()),_e1(keepOrder?e1:e2),
211                                        _e2(keepOrder?e2:e1)
212 {
213 }
214
215 IntersectElement::IntersectElement(const IntersectElement& other):_1S(other._1S),_1E(other._1E),_2S(other._2S),_2E(other._2E),
216     _chararct_val_for_e1(other._chararct_val_for_e1),
217     _chararct_val_for_e2(other._chararct_val_for_e2),_node(other._node),
218     _loc_of_node(other._loc_of_node),_e1(other._e1), _e2(other._e2)
219 {
220   if(_node)
221     _node->incrRef();
222 }
223
224 IntersectElement& IntersectElement::operator=(const IntersectElement& other)
225 {
226   _1S=other._1S;_1E=other._1E; _2S=other._2S; _2E=other._2E;
227   _chararct_val_for_e1=other._chararct_val_for_e1;
228   _chararct_val_for_e2=other._chararct_val_for_e2;
229   setNode(other._node);
230   return *this;
231 }
232
233 bool IntersectElement::operator<(const IntersectElement& other) const
234 {
235   return _e1.isLower(_chararct_val_for_e1,other._chararct_val_for_e1);
236 }
237
238 IntersectElement::~IntersectElement()
239 {
240   if(_node)
241     _node->decrRef();
242 }
243
244 /*!
245  * Returns 0 or 1.
246  */
247 bool IntersectElement::isOnMergedExtremity() const
248 {
249   if( (_1S && _2S) || (_1S && _2E) || (_1E && _2S) || (_1E && _2E) )
250     return true;
251   return false;
252 }
253
254 /*!
255  * To call if isOnMergedExtremity returned true.
256  */
257 void IntersectElement::performMerging(MergePoints& commonNode) const
258 {
259   if(_1S && _2S)
260     {
261       if(_e1.changeStartNodeWith(_e2.getStartNode()))
262         {
263           _e2.getStartNode()->declareOnLim();
264           commonNode.start1OnStart2();
265         }
266     }
267   else if(_1S && _2E)
268     {
269       if(_e1.changeStartNodeWith(_e2.getEndNode()))
270         {
271           _e2.getEndNode()->declareOnLim();
272           commonNode.start1OnEnd2();
273         }
274     }
275   else if(_1E && _2S)
276     {
277       if(_e1.changeEndNodeWith(_e2.getStartNode()))
278         {
279           _e2.getStartNode()->declareOnLim();
280           commonNode.end1OnStart2();
281         }
282     }
283   else if(_1E && _2E)
284     {
285       if(_e1.changeEndNodeWith(_e2.getEndNode()))
286         {
287           _e2.getEndNode()->declareOnLim();
288           commonNode.end1OnEnd2();
289         }
290     }
291 }
292
293 /*!
294  * This method is const because 'node' is supposed to be equal geometrically to _node.
295  */
296 void IntersectElement::setNode(Node *node) const
297 {
298   if(node!=_node)
299     {
300       if(_node)
301         ((Node *)_node)->decrRef();
302       (const_cast<IntersectElement *>(this))->_node=node;
303       if(_node)
304         _node->incrRef();
305     }
306 }
307
308 bool IntersectElement::isLowerOnOther(const IntersectElement& other) const
309 {
310   return _e2.isLower(_chararct_val_for_e2,other._chararct_val_for_e2);
311 }
312
313 unsigned IntersectElement::isOnExtrForAnEdgeAndInForOtherEdge() const
314 {
315   if(( _1S && !(_2S || _2E) ) || ( _1E && !(_2S || _2E) ))
316     {
317       if(_1S && !(_2S || _2E))
318         setNode(_e1.getStartNode());
319       else
320         setNode(_e1.getEndNode());
321       if(_e2.isIn(_chararct_val_for_e2))
322         return LIMIT_ON;
323       return LIMIT_ALONE;
324     }
325   if(( _2S && !(_1S || _1E) ) || ( _2E && !(_1S || _1E)))
326     {
327       if(_2S && !(_1S || _1E))
328         setNode(_e2.getStartNode());
329       else
330         setNode(_e2.getEndNode());
331       if(_e1.isIn(_chararct_val_for_e1))
332         return LIMIT_ON;
333       return LIMIT_ALONE;
334     }
335   return NO_LIMIT;
336 }
337
338 bool IntersectElement::isIncludedByBoth() const
339 {
340   return _e1.isIn(_chararct_val_for_e1) && _e2.isIn(_chararct_val_for_e2);
341 }
342
343 bool EdgeIntersector::intersect(std::vector<Node *>& newNodes, bool& order, MergePoints& commonNode)
344 {
345   std::list< IntersectElement > listOfIntesc=getIntersectionsCharacteristicVal();
346   std::list< IntersectElement >::iterator iter;
347   for(iter=listOfIntesc.begin();iter!=listOfIntesc.end();)
348     {
349       if((*iter).isOnMergedExtremity())
350         {
351           (*iter).performMerging(commonNode);
352           iter=listOfIntesc.erase(iter);
353           continue;
354         }
355       unsigned tmp=(*iter).isOnExtrForAnEdgeAndInForOtherEdge();
356       if(tmp==IntersectElement::LIMIT_ALONE)
357         {
358           iter=listOfIntesc.erase(iter);
359           continue;
360         }
361       else if(tmp==IntersectElement::LIMIT_ON)
362         {
363           (*iter).attachLoc();
364           iter++;
365           continue;
366         }
367       if(!(*iter).isIncludedByBoth())
368         {
369           iter=listOfIntesc.erase(iter);
370           continue;
371         }
372       (*iter).attachLoc();
373       iter++;
374     }
375   if(listOfIntesc.size()==0)
376     return false;
377   if(listOfIntesc.size()==1)
378     {
379       order=true;//useless
380       newNodes.push_back(listOfIntesc.front().getNodeAndReleaseIt());
381     }
382   else
383     {
384       std::vector<IntersectElement> vecOfIntesc(listOfIntesc.begin(),listOfIntesc.end());
385       listOfIntesc.clear();
386       sort(vecOfIntesc.begin(),vecOfIntesc.end());
387       for(std::vector<IntersectElement>::iterator iterV=vecOfIntesc.begin();iterV!=vecOfIntesc.end();iterV++)
388         newNodes.push_back((*iterV).getNodeAndReleaseIt());
389       order=vecOfIntesc.front().isLowerOnOther(vecOfIntesc.back());
390     }
391   return true;
392 }
393
394 /*! If the 2 edges share one extremity, we can optimize since we already know where is the intersection.
395  *  In the case of ArcCSegIntersector, this also helps avoid degenerated cases.
396  */
397 void EdgeIntersector::identifyEarlyIntersection(bool& i1S2S, bool& i1E2S, bool& i1S2E, bool& i1E2E)
398 {
399   i1S2S = _e1.getStartNode() == _e2.getStartNode();
400   i1E2S = _e1.getEndNode() == _e2.getStartNode();
401   i1S2E = _e1.getStartNode() == _e2.getEndNode();
402   i1E2E = _e1.getEndNode() == _e2.getEndNode();
403   if (i1S2S || i1E2S || i1S2E || i1E2E)
404     {
405       Node * node;
406       bool i_1S(false),i_1E(false),i_2S(false),i_2E(false);
407       if (i1S2S || i1E2S)   // Common node is e2 start
408         {
409           node = _e2.getStartNode();
410           i_1S = i1S2S;        i_2S = true;
411           i_1E = i1E2S;        i_2E = false;
412         }
413       else                  // Common node is e2 end
414         {
415           node = _e2.getEndNode();
416           i_1S = i1S2E;        i_2S = false;
417           i_1E = i1E2E;        i_2E = true;
418         }
419       node->incrRef();
420       _earlyInter = new IntersectElement(_e1.getCharactValue(*node), _e2.getCharactValue(*node),
421           i_1S,i_1E,i_2S,i_2E,node,_e1,_e2,keepOrder());
422     }
423 }
424
425 /*!
426  * Locates 'node' regarding edge this->_e1. If node is located close to (with distant lt epsilon) start or end point of _e1,
427  * 'node' takes its place. In this case 'obvious' is set to true and 'commonNode' stores information of merge point and finally 'where' is set.
428  * Furthermore 'node' is declared as ON LIMIT to indicate in locating process that an absolute location computation will have to be done.
429  * If 'node' is not close to start or end point of _e1, 'obvious' is set to false and 'commonNode' and 'where' are let unchanged.
430  */
431 void EdgeIntersector::obviousCaseForCurvAbscisse(Node *node, TypeOfLocInEdge& where, MergePoints& commonNode, bool& obvious) const
432 {
433   obvious=true;
434   if(node->isEqual(*_e1.getStartNode()))
435     {
436       where=START;
437       if(_e1.changeStartNodeWith(node))
438         {
439           commonNode.start1Replaced();
440           node->declareOnLim();
441         }
442       return ;
443     }
444   if(node->isEqual(*_e1.getEndNode()))
445     {
446       where=END;
447       if(_e1.changeEndNodeWith(node))
448         {
449           commonNode.end1Replaced();
450           node->declareOnLim();
451         }
452       return ;
453     }
454   obvious=false;
455 }
456
457 Edge::Edge(double sX, double sY, double eX, double eY):_cnt(1),_loc(FULL_UNKNOWN),_start(new Node(sX,sY)),_end(new Node(eX,eY))
458 {
459 }
460
461 Edge::~Edge()
462 {
463   _start->decrRef();
464   if(_end)
465     _end->decrRef();
466 }
467
468 bool Edge::decrRef()
469 {
470   bool ret=(--_cnt==0);
471   if(ret)
472     delete this;
473   return ret;
474 }
475
476 void Edge::declareOn() const
477 {
478   if(_loc==FULL_UNKNOWN)
479     {
480       _loc=FULL_ON_1;
481       _start->declareOn();
482       _end->declareOn();
483     }
484 }
485
486 void Edge::declareIn() const
487 {
488   if(_loc==FULL_UNKNOWN)
489     {
490       _loc=FULL_IN_1;
491       _start->declareIn();
492       _end->declareIn();
493     }
494 }
495
496 void Edge::declareOut() const
497 {
498   if(_loc==FULL_UNKNOWN)
499     {
500       _loc=FULL_OUT_1;
501       _start->declareOut();
502       _end->declareOut();
503     }
504 }
505
506 void Edge::fillXfigStreamForLoc(std::ostream& stream) const
507 {
508   switch(_loc)
509   {
510     case FULL_IN_1:
511       stream << '2';//Green
512       break;
513     case FULL_OUT_1:
514       stream << '1';//Bleue
515       break;
516     case FULL_ON_1:
517       stream << '4';//Red
518       break;
519     default:
520       stream << '0';
521   }
522 }
523
524 bool Edge::changeStartNodeWith(Node *otherStartNode) const
525 {
526   if(_start==otherStartNode)
527     return true;
528   if(_start->isEqual(*otherStartNode))
529     {
530       ((const_cast<Edge *>(this))->_start)->decrRef();//un-const cast Ok thanks to 2 lines above.
531       ((const_cast<Edge *>(this))->_start)=otherStartNode;
532       _start->incrRef();
533       return true;
534     }
535   return false;
536 }
537
538 bool Edge::changeStartNodeWithAndKeepTrack(Node *otherStartNode, std::vector<Node *>& track) const
539 {
540   if(_start==otherStartNode)
541     return true;
542   if(_start->isEqualAndKeepTrack(*otherStartNode,track))
543     {
544       ((const_cast<Edge *>(this))->_start)->decrRef();//un-const cast Ok thanks to 2 lines above.
545       ((const_cast<Edge *>(this))->_start)=otherStartNode;
546       otherStartNode->incrRef();
547       return true;
548     }
549   return false;
550 }
551
552 bool Edge::changeEndNodeWith(Node *otherEndNode) const
553 {
554   if(_end==otherEndNode)
555     return true;
556   if(_end->isEqual(*otherEndNode))
557     {
558       ((const_cast<Edge *>(this))->_end)->decrRef();
559       ((const_cast<Edge *>(this))->_end)=otherEndNode;
560       _end->incrRef();
561       return true;
562     }
563   return false;
564 }
565
566 bool Edge::changeEndNodeWithAndKeepTrack(Node *otherEndNode, std::vector<Node *>& track) const
567 {
568   if(_end==otherEndNode)
569     return true;
570   if(_end->isEqualAndKeepTrack(*otherEndNode,track))
571     {
572       ((const_cast<Edge *>(this))->_end)->decrRef();
573       ((const_cast<Edge *>(this))->_end)=otherEndNode;
574       otherEndNode->incrRef();
575       return true;
576     }
577   return false;
578 }
579
580 /*!
581  * Precondition : 'start' and 'end' are lying on the same curve than 'this'.
582  * Add in vec the sub edge lying on this.
583  * If 'start' is equal (by pointer) to '_end' and 'end' is equal to '_end' too nothing is added.
584  * If 'start' is equal (by pointer) to '_start' and 'end' is equal to '_start' too nothing is added.
585  * If 'start' is equal (by pointer) to '_start' and 'end' is equal to '_end' this is added in vec.
586  */
587 void Edge::addSubEdgeInVector(Node *start, Node *end, ComposedEdge& vec) const
588 {
589   if((start==_start && end==_start) || (start==_end && end==_end))
590     return ;
591   if(start==_start && end==_end)
592     {
593       incrRef();
594       vec.pushBack(const_cast<Edge *>(this));
595       return ;
596     }
597   vec.pushBack(buildEdgeLyingOnMe(start,end,true));
598 }
599
600 /*!
601  * Retrieves a vector 'vectOutput' that is normal to 'this'. 'vectOutput' is normalized.
602  */
603 void Edge::getNormalVector(double *vectOutput) const
604 {
605   std::copy((const double *)(*_end),(const double *)(*_end)+2,vectOutput);
606   std::transform(vectOutput,vectOutput+2,(const double *)(*_start),vectOutput,std::minus<double>());
607   double norm=1./Node::norm(vectOutput);
608   std::transform(vectOutput,vectOutput+2,vectOutput,bind(std::multiplies<double>(),std::placeholders::_1,norm));
609   double tmp=vectOutput[0];
610   vectOutput[0]=vectOutput[1];
611   vectOutput[1]=-tmp;
612 }
613
614 Edge *Edge::BuildEdgeFrom3Points(const double *start, const double *middle, const double *end)
615 {
616   Node *b(new Node(start[0],start[1])),*m(new Node(middle[0],middle[1])),*e(new Node(end[0],end[1]));
617   EdgeLin *e1(new EdgeLin(b,m)),*e2(new EdgeLin(m,e));
618   SegSegIntersector inters(*e1,*e2); bool colinearity=inters.areColinears(); delete e1; delete e2;
619   Edge *ret=0;
620   if(colinearity)
621     ret=new EdgeLin(b,e);
622   else
623     ret=new EdgeArcCircle(b,m,e);
624   b->decrRef(); m->decrRef(); e->decrRef();
625   return ret;
626 }
627
628 Edge *Edge::BuildEdgeFrom(Node *start, Node *end)
629 {
630   return new EdgeLin(start,end);
631 }
632
633 Edge *Edge::BuildFromXfigLine(std::istream& str)
634 {
635   unsigned char type;
636   str >> type;
637   if(type=='2')
638     return new EdgeLin(str);
639   else if(type=='5')
640     return new EdgeArcCircle(str);
641   else
642     {
643       std::cerr << "Unknown line found...";
644       return 0;
645     }
646 }
647
648 /*!
649  * \param other The Edge with which we are going to intersect.
650  * \param commonNode Output. The common nodes found during operation of intersecting.
651  * \param outVal1 Output filled in case true is returned. It specifies the new or not new edges by which 'this' is replaced after intersecting op.
652  * \param outVal2 Output filled in case true is returned. It specifies the new or not new edges by which 'other' is replaced after intersecting op.
653  * return true if the intersection between this.
654  */
655 bool Edge::intersectWith(const Edge *other, MergePoints& commonNode,
656                          ComposedEdge& outVal1, ComposedEdge& outVal2) const
657 {
658   bool ret=true;
659   Bounds *merge=_bounds.nearlyAmIIntersectingWith(other->getBounds());
660   if(!merge)
661     return false;
662   delete merge;
663   merge=0;
664   EdgeIntersector *intersector=BuildIntersectorWith(this,other);
665   ret=Intersect(this,other,intersector,commonNode,outVal1,outVal2);
666   delete intersector;
667   return ret;
668 }
669
670 bool Edge::IntersectOverlapped(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, MergePoints& commonNode,
671                                ComposedEdge& outValForF1, ComposedEdge& outValForF2)
672 {
673   bool rev=intersector->haveTheySameDirection();
674   Node *f2Start=f2->getNode(rev?START:END);
675   Node *f2End=f2->getNode(rev?END:START);
676   TypeOfLocInEdge place1, place2;
677   intersector->getPlacements(f2Start,f2End,place1,place2,commonNode);
678   int codeForIntersectionCase=CombineCodes(place1,place2);
679   return SplitOverlappedEdges(f1,f2,f2Start,f2End,rev,codeForIntersectionCase,outValForF1,outValForF2);
680 }
681
682 /*!
683  * Perform 1D linear interpolation. Warning distrib1 and distrib2 are expected to be in ascending mode.
684  */
685 void Edge::Interpolate1DLin(const std::vector<double>& distrib1, const std::vector<double>& distrib2, std::map<int, std::map<int,double> >& result)
686 {
687   std::size_t nbOfV1=distrib1.size()-1;
688   std::size_t nbOfV2=distrib2.size()-1;
689   Node *n1=new Node(0.,0.); Node *n3=new Node(0.,0.);
690   Node *n2=new Node(0.,0.); Node *n4=new Node(0.,0.);
691   MergePoints commonNode;
692   for(unsigned int i=0;i<nbOfV1;i++)
693     {
694       std::vector<double>::const_iterator iter=find_if(distrib2.begin()+1,distrib2.end(),bind(std::greater_equal<double>(),std::placeholders::_1,distrib1[i]));
695       if(iter!=distrib2.end())
696         {
697           for(unsigned int j=(unsigned)((iter-1)-distrib2.begin());j<nbOfV2;j++)
698             {
699               if(distrib2[j]<=distrib1[i+1])
700                 {
701                   EdgeLin *e1=new EdgeLin(n1,n2); EdgeLin *e2=new EdgeLin(n3,n4);
702                   n1->setNewCoords(distrib1[i],0.); n2->setNewCoords(distrib1[i+1],0.);
703                   n3->setNewCoords(distrib2[j],0.); n4->setNewCoords(distrib2[j+1],0.);
704                   ComposedEdge *f1=new ComposedEdge;
705                   ComposedEdge *f2=new ComposedEdge;
706                   SegSegIntersector inters(*e1,*e2);
707                   bool b1,b2;
708                   inters.areOverlappedOrOnlyColinears(b1,b2);
709                   if(IntersectOverlapped(e1,e2,&inters,commonNode,*f1,*f2))
710                     {
711                       result[i][j]=f1->getCommonLengthWith(*f2)/e1->getCurveLength();
712                     }
713                   ComposedEdge::Delete(f1); ComposedEdge::Delete(f2);
714                   e1->decrRef(); e2->decrRef();
715                 }
716             }
717         }
718     }
719   n1->decrRef(); n2->decrRef(); n3->decrRef(); n4->decrRef();
720 }
721
722 EdgeIntersector *Edge::BuildIntersectorWith(const Edge *e1, const Edge *e2)
723 {
724   EdgeIntersector *ret=0;
725   const EdgeLin *tmp1=0;
726   const EdgeArcCircle *tmp2=0;
727   unsigned char type1=e1->getTypeOfFunc();
728   e1->dynCastFunction(tmp1,tmp2);
729   unsigned char type2=e2->getTypeOfFunc();
730   e2->dynCastFunction(tmp1,tmp2);
731   type1|=type2;
732   switch(type1)
733   {
734     case 1:// Intersection seg/seg
735       ret=new SegSegIntersector((const EdgeLin &)(*e1),(const EdgeLin &)(*e2));
736       break;
737     case 5:// Intersection seg/arc of circle
738       ret=new ArcCSegIntersector(*tmp2,*tmp1,tmp2==e1);
739       break;
740     case 4:// Intersection arc/arc of circle
741       ret=new ArcCArcCIntersector((const EdgeArcCircle &)(*e1),(const EdgeArcCircle &)(*e2));
742       break;
743     default:
744       //Should never happen
745       throw Exception("A non managed association of edge has been detected. Go work for intersection computation implementation.");
746   }
747   return ret;
748 }
749
750 /*!
751  * See Node::applySimilarity to see signification of params.
752  */
753 void Edge::applySimilarity(double xBary, double yBary, double dimChar)
754 {
755   _bounds.applySimilarity(xBary,yBary,dimChar);
756 }
757
758 void Edge::unApplySimilarity(double xBary, double yBary, double dimChar)
759 {
760   _bounds.unApplySimilarity(xBary,yBary,dimChar);
761 }
762
763 void Edge::getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const
764 {
765   return getMiddleOfPoints(p1, p2, mid);
766 }
767
768 bool Edge::Intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, MergePoints& commonNode,
769                      ComposedEdge& outValForF1, ComposedEdge& outValForF2)
770 {
771   bool obviousNoIntersection;
772   bool areOverlapped;
773   intersector->areOverlappedOrOnlyColinears(obviousNoIntersection,areOverlapped);
774   if(areOverlapped)
775     return IntersectOverlapped(f1,f2,intersector,commonNode,outValForF1,outValForF2);
776   if(obviousNoIntersection)
777     return false;
778   std::vector<Node *> newNodes;
779   bool order;
780   if(intersector->intersect(newNodes,order,commonNode))
781     {
782       if(newNodes.empty())
783         throw Exception("Internal error occurred - error in intersector implementation!");// This case should never happen
784       std::vector<Node *>::iterator iter=newNodes.begin();
785       std::vector<Node *>::reverse_iterator iterR=newNodes.rbegin();
786       f1->addSubEdgeInVector(f1->getStartNode(),*iter,outValForF1);
787       f2->addSubEdgeInVector(f2->getStartNode(),order?*iter:*iterR,outValForF2);
788       for(std::vector<Node *>::iterator iter2=newNodes.begin();iter2!=newNodes.end();iter2++,iterR++)
789         {
790           if((iter2+1)==newNodes.end())
791             {
792               f1->addSubEdgeInVector(*iter2,f1->getEndNode(),outValForF1);
793               (*iter2)->decrRef();
794               f2->addSubEdgeInVector(order?*iter2:*iterR,f2->getEndNode(),outValForF2);
795             }
796           else
797             {
798               f1->addSubEdgeInVector(*iter2,*(iter2+1),outValForF1);
799               (*iter2)->decrRef();
800               f2->addSubEdgeInVector(order?*iter2:*iterR,order?*(iter2+1):*(iterR+1),outValForF2);
801             }
802         }
803       return true;
804     }
805   else
806     return false;
807 }
808
809 int Edge::CombineCodes(TypeOfLocInEdge code1, TypeOfLocInEdge code2)
810 {
811   int ret=(int)code1;
812   ret*=OFFSET_FOR_TYPEOFLOCINEDGE;
813   ret+=(int)code2;
814   return ret;
815 }
816
817 /*!
818  * This method splits e1 and e2 into pieces as much sharable as possible. The precondition to the call of this method
819  * is that e1 and e2 have been declared as overlapped by corresponding intersector built from e1 and e2 type.
820  *
821  * @param nS start node of e2 with the SAME DIRECTION as e1. The pointer nS should be equal to start node of e2 or to its end node.
822  * @param nE end node of e2 with the SAME DIRECTION as e1. The pointer nE should be equal to start node of e2 or to its end node.
823  * @param direction is param that specifies if e2 and e1 have same directions (true) or opposed (false).
824  * @param code is the code returned by method Edge::combineCodes.
825  */
826 bool Edge::SplitOverlappedEdges(const Edge *e1, const Edge *e2, Node *nS, Node *nE, bool direction, int code,
827                                 ComposedEdge& outVal1, ComposedEdge& outVal2)
828 {
829   Edge *tmp;
830   switch(code)
831   {
832     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+START:      // OUT_BEFORE - START
833     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_BEFORE: // OUT_BEFORE - OUT_BEFORE
834     case OUT_AFTER*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:   // OUT_AFTER - OUT_AFTER
835     case END*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:         // END - OUT_AFTER
836     case END*OFFSET_FOR_TYPEOFLOCINEDGE+START:             // END - START
837     return false;
838     case INSIDE*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:      // INSIDE - OUT_AFTER
839     outVal1.pushBack(e1->buildEdgeLyingOnMe(e1->getStartNode(),nS,true));
840     tmp=e1->buildEdgeLyingOnMe(nS,e1->getEndNode()); tmp->incrRef();
841     outVal1.pushBack(tmp);
842     outVal2.resize(2);
843     outVal2.setValueAt(direction?0:1,tmp,direction); tmp->declareOn();
844     outVal2.setValueAt(direction?1:0,e1->buildEdgeLyingOnMe(e1->getEndNode(),nE,direction));
845     return true;
846     case INSIDE*OFFSET_FOR_TYPEOFLOCINEDGE+INSIDE:         // INSIDE - INSIDE
847     {
848       if(!e2->isIn(e2->getCharactValue(*(e1->getStartNode()))))
849         {
850           e2->incrRef(); e2->incrRef();
851           outVal1.resize(3);
852           outVal1.setValueAt(0,e1->buildEdgeLyingOnMe(e1->getStartNode(),nS));
853           outVal1.setValueAt(1,const_cast<Edge*>(e2),direction);
854           outVal1.setValueAt(2,e1->buildEdgeLyingOnMe(nE,e1->getEndNode()));
855           outVal2.pushBack(const_cast<Edge*>(e2)); e2->declareOn();
856           return true;
857         }
858       else
859         {
860           outVal1.resize(3);
861           outVal2.resize(3);
862           tmp=e1->buildEdgeLyingOnMe(e1->getStartNode(),nE); tmp->incrRef(); tmp->declareOn();
863           outVal1.setValueAt(0,tmp,true); outVal2.setValueAt(direction?2:0,tmp,direction);
864           outVal1.setValueAt(1,e1->buildEdgeLyingOnMe(nE,nS));
865           tmp=e1->buildEdgeLyingOnMe(nS,e1->getEndNode()); tmp->incrRef(); tmp->declareOn();
866           outVal1.setValueAt(2,tmp,true); outVal2.setValueAt(direction?0:2,tmp,direction);
867           tmp=e1->buildEdgeLyingOnMe(e1->getEndNode(),e1->getStartNode());
868           outVal2.setValueAt(1,tmp,direction);
869           return true;
870         }
871     }
872     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+INSIDE:     // OUT_BEFORE - INSIDE
873     {
874       tmp=e1->buildEdgeLyingOnMe(e1->getStartNode(),nE); tmp->incrRef();
875       outVal1.pushBack(tmp);
876       outVal1.pushBack(e1->buildEdgeLyingOnMe(nE,e1->getEndNode()));
877       outVal2.resize(2);
878       outVal2.setValueAt(direction?0:1,e1->buildEdgeLyingOnMe(nS,e1->getStartNode(),direction));
879       outVal2.setValueAt(direction?1:0,tmp,direction); tmp->declareOn();
880       return true;
881     }
882     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:  // OUT_BEFORE - OUT_AFTER
883     {
884       e1->incrRef(); e1->incrRef();
885       outVal1.pushBack(const_cast<Edge*>(e1));
886       outVal2.resize(3);
887       outVal2.setValueAt(direction?0:2,e1->buildEdgeLyingOnMe(nS,e1->getStartNode(),direction));
888       outVal2.setValueAt(1,const_cast<Edge*>(e1),direction); e1->declareOn();
889       outVal2.setValueAt(direction?2:0,e1->buildEdgeLyingOnMe(e1->getEndNode(),nE,direction));
890       return true;
891     }
892     case START*OFFSET_FOR_TYPEOFLOCINEDGE+END:             // START - END
893     {
894       e1->incrRef(); e1->incrRef();
895       outVal1.pushBack(const_cast<Edge*>(e1));
896       outVal2.pushBack(const_cast<Edge*>(e1),direction); e1->declareOn();
897       return true;
898     }
899     case START*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:       // START - OUT_AFTER
900     {
901       e1->incrRef(); e1->incrRef();
902       outVal1.pushBack(const_cast<Edge*>(e1));
903       outVal2.resize(2);
904       outVal2.setValueAt(direction?0:1,const_cast<Edge*>(e1),direction); e1->declareOn();
905       outVal2.setValueAt(direction?1:0,e1->buildEdgeLyingOnMe(e1->getEndNode(),nE,direction));
906       return true;
907     }
908     case INSIDE*OFFSET_FOR_TYPEOFLOCINEDGE+END:            // INSIDE - END
909     {
910       e2->incrRef(); e2->incrRef();
911       outVal1.pushBack(e1->buildEdgeLyingOnMe(e1->getStartNode(),nS,true));
912       outVal1.pushBack(const_cast<Edge*>(e2),direction);
913       outVal2.pushBack(const_cast<Edge*>(e2)); e2->declareOn();
914       return true;
915     }
916     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+END:        // OUT_BEFORE - END
917     {
918       e1->incrRef(); e1->incrRef();
919       outVal1.pushBack(const_cast<Edge*>(e1));
920       outVal2.resize(2);
921       outVal2.setValueAt(direction?0:1,e1->buildEdgeLyingOnMe(nS,e1->getStartNode(),direction));
922       outVal2.setValueAt(direction?1:0,const_cast<Edge*>(e1),direction); e1->declareOn();
923       return true;
924     }
925     case START*OFFSET_FOR_TYPEOFLOCINEDGE+INSIDE:          // START - INSIDE
926     {
927       e2->incrRef(); e2->incrRef();
928       outVal1.pushBack(const_cast<Edge*>(e2),direction);
929       outVal1.pushBack(e1->buildEdgeLyingOnMe(nE,e1->getEndNode()));
930       outVal2.pushBack(const_cast<Edge*>(e2)); e2->declareOn();
931       return true;
932     }
933     case INSIDE*OFFSET_FOR_TYPEOFLOCINEDGE+START:          // INSIDE - START
934     {
935       outVal1.resize(2);
936       outVal2.resize(2);
937       tmp=e1->buildEdgeLyingOnMe(nS,e1->getEndNode()); tmp->incrRef(); tmp->declareOn();
938       outVal1.setValueAt(0,e1->buildEdgeLyingOnMe(e1->getStartNode(),nS));
939       outVal1.setValueAt(1,tmp);
940       outVal2.setValueAt(direction?0:1,tmp,direction);
941       outVal2.setValueAt(direction?1:0,e1->buildEdgeLyingOnMe(e1->getEndNode(),nE,direction));
942       return true;
943     }
944     case END*OFFSET_FOR_TYPEOFLOCINEDGE+INSIDE:            // END - INSIDE
945     {
946       outVal1.resize(2);
947       outVal2.resize(2);
948       tmp=e1->buildEdgeLyingOnMe(e1->getStartNode(),nE); tmp->incrRef(); tmp->declareOn();
949       outVal1.setValueAt(0,tmp);
950       outVal1.setValueAt(1,e1->buildEdgeLyingOnMe(nE,e1->getEndNode()));
951       outVal2.setValueAt(direction?0:1,e1->buildEdgeLyingOnMe(e1->getEndNode(),e1->getStartNode(),direction));
952       outVal2.setValueAt(direction?1:0,tmp,direction);
953       return true;
954     }
955     default:
956       throw Exception("Unexpected situation of overlapping edges : internal error occurs ! ");
957   }
958 }
959
960 void Edge::dumpToCout(const std::map<INTERP_KERNEL::Node *,int>& mapp, int index) const
961 {
962   auto sI(mapp.find(getStartNode())), eI(mapp.find(getEndNode()));
963   int start = (sI == mapp.end() ? -1 : sI->second), end = (eI == mapp.end() ? -1 : eI->second);
964   std::string locs;
965   switch (getLoc())
966   {
967     case FULL_IN_1: locs="FULL_IN_1"; break;
968     case FULL_ON_1: locs="FULL_ON_1"; break;
969     case FULL_OUT_1: locs="FULL_OUT_1"; break;
970     case FULL_UNKNOWN: locs="FULL_UNKNOWN"; break;
971     default: locs="oh my God! This is so wrong.";
972   }
973   std::cout << "Edge [" << index << "] : ("<<  std::hex << this << std::dec << ") -> (" << start << ", " << end << ")\t" << locs << std::endl;
974 }
975
976 bool Edge::isEqual(const Edge& other) const
977 {
978   return _start->isEqual(*other._start) && _end->isEqual(*other._end);
979 }
980
981 /**
982  * This method takes in input nodes in \a subNodes (using \a coo)
983  *
984  * \param [in,out] subNodes to be sorted
985  * \return true if a reordering was necessary false if not.
986  */
987 bool Edge::sortSubNodesAbs(const double *coo, std::vector<mcIdType>& subNodes)
988 {
989   Bounds b;
990   b.prepareForAggregation();
991   b.aggregate(getBounds());
992   double xBary,yBary;
993   double dimChar(b.getCaracteristicDim());
994   b.getBarycenter(xBary,yBary);
995   applySimilarity(xBary,yBary,dimChar);
996   _start->applySimilarity(xBary,yBary,dimChar);
997   _end->applySimilarity(xBary,yBary,dimChar);
998   //
999   std::size_t sz(subNodes.size()),i(0);
1000   std::vector< std::pair<double,Node *> > an2(sz);
1001   std::map<Node *, mcIdType> m;
1002   for(std::vector<mcIdType>::const_iterator it=subNodes.begin();it!=subNodes.end();it++,i++)
1003     {
1004       Node *n(new Node(coo[2*(*it)],coo[2*(*it)+1]));
1005       n->applySimilarity(xBary,yBary,dimChar);
1006       m[n]=*it;
1007       an2[i]=std::pair<double,Node *>(getCharactValueBtw0And1(*n),n);
1008     }
1009   std::sort(an2.begin(),an2.end());
1010   //
1011   bool ret(false);
1012   for(i=0;i<sz;i++)
1013     {
1014       mcIdType id(m[an2[i].second]);
1015       if(id!=subNodes[i])
1016         { subNodes[i]=id; ret=true; }
1017     }
1018   //
1019   for(std::map<INTERP_KERNEL::Node *,mcIdType>::const_iterator it2=m.begin();it2!=m.end();it2++)
1020     (*it2).first->decrRef();
1021   return ret;
1022 }
1023
1024 /**
1025  * Sort nodes so that they all lie consecutively on the edge that has been cut.
1026  */
1027 void Edge::sortIdsAbs(const std::vector<INTERP_KERNEL::Node *>& addNodes, const std::map<INTERP_KERNEL::Node *, mcIdType>& mapp1,
1028                       const std::map<INTERP_KERNEL::Node *, mcIdType>& mapp2, std::vector<mcIdType>& edgesThis)
1029 {
1030   mcIdType startId=(*mapp1.find(_start)).second;
1031   mcIdType endId=(*mapp1.find(_end)).second;
1032   if (! addNodes.size()) // quick way out, no new node to add.
1033     {
1034       edgesThis.push_back(startId);
1035       edgesThis.push_back(endId);
1036       return;
1037     }
1038
1039   Bounds b;
1040   b.prepareForAggregation();
1041   b.aggregate(getBounds());
1042   double xBary,yBary;
1043   double dimChar=b.getCaracteristicDim();
1044   b.getBarycenter(xBary,yBary);
1045   for(std::vector<Node *>::const_iterator iter=addNodes.begin();iter!=addNodes.end();iter++)
1046     (*iter)->applySimilarity(xBary,yBary,dimChar);
1047   applySimilarity(xBary,yBary,dimChar);
1048   _start->applySimilarity(xBary,yBary,dimChar);
1049   _end->applySimilarity(xBary,yBary,dimChar);
1050   std::size_t sz=addNodes.size();
1051   std::vector< std::pair<double,Node *> > an2(sz);
1052   for(std::size_t i=0;i<sz;i++)
1053     an2[i]=std::pair<double,Node *>(getCharactValueBtw0And1(*addNodes[i]),addNodes[i]);
1054   std::sort(an2.begin(),an2.end());
1055   std::vector<mcIdType> tmpp;
1056   for(std::vector< std::pair<double,Node *> >::const_iterator it=an2.begin();it!=an2.end();it++)
1057     {
1058       mcIdType idd=(*mapp2.find((*it).second)).second;
1059       tmpp.push_back(idd);
1060     }
1061   std::vector<mcIdType> tmpp2(tmpp.size()+2);
1062   tmpp2[0]=startId;
1063   std::copy(tmpp.begin(),tmpp.end(),tmpp2.begin()+1);
1064   tmpp2[tmpp.size()+1]=endId;
1065   std::vector<mcIdType>::iterator itt=std::unique(tmpp2.begin(),tmpp2.end());
1066   tmpp2.resize(std::distance(tmpp2.begin(),itt));
1067   std::size_t nbOfEdges=tmpp2.size()-1;
1068   for(std::size_t i=0;i<nbOfEdges;i++)
1069     {
1070       edgesThis.push_back(tmpp2[i]);
1071       edgesThis.push_back(tmpp2[i+1]);
1072     }
1073 }
1074
1075 void Edge::fillGlobalInfoAbs(bool direction, const std::map<INTERP_KERNEL::Node *,mcIdType>& mapThis, const std::map<INTERP_KERNEL::Node *,mcIdType>& mapOther, mcIdType offset1, mcIdType offset2, double fact, double baryX, double baryY,
1076                                 std::vector<mcIdType>& edgesThis, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,mcIdType> mapAddCoo) const
1077 {
1078   mcIdType tmp[2];
1079   _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp);
1080   _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1);
1081   if(direction)
1082     {
1083       edgesThis.push_back(tmp[0]);
1084       edgesThis.push_back(tmp[1]);
1085     }
1086   else
1087     {
1088       edgesThis.push_back(tmp[1]);
1089       edgesThis.push_back(tmp[0]);
1090     }
1091 }
1092
1093 void Edge::fillGlobalInfoAbs2(const std::map<INTERP_KERNEL::Node *,mcIdType>& mapThis, const std::map<INTERP_KERNEL::Node *,mcIdType>& mapOther, mcIdType offset1, mcIdType offset2, double fact, double baryX, double baryY,
1094                               short skipStartOrEnd,
1095                               std::vector<mcIdType>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,mcIdType>& mapAddCoo) const
1096 {
1097   if (skipStartOrEnd != -1) // see meaning in splitAbs()
1098     _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
1099   if (skipStartOrEnd != 1)
1100   _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
1101 }