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