1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
21 #include "InterpKernelGeo2DElementaryEdge.hxx"
22 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
23 #include "InterpKernelGeo2DAbstractEdge.hxx"
24 #include "InterpKernelGeo2DEdgeLin.hxx"
25 #include "InterpKernelGeo2DBounds.hxx"
26 #include "InterpKernelGeo2DEdge.txx"
28 #include "NormalizedUnstructuredMesh.hxx"
35 using namespace INTERP_KERNEL;
37 namespace INTERP_KERNEL
39 const unsigned MAX_SIZE_OF_LINE_XFIG_FILE=1024;
42 QuadraticPolygon::QuadraticPolygon(const char *file)
44 char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
45 std::ifstream stream(file);
46 stream.exceptions(std::ios_base::eofbit);
50 stream.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
51 while(strcmp(currentLine,"1200 2")!=0);
54 Edge *newEdge=Edge::BuildFromXfigLine(stream);
56 newEdge->changeStartNodeWith(back()->getEndNode());
61 catch(std::ifstream::failure&)
64 front()->changeStartNodeWith(back()->getEndNode());
67 QuadraticPolygon::~QuadraticPolygon()
71 QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes)
73 QuadraticPolygon *ret=new QuadraticPolygon;
74 std::size_t size=nodes.size();
75 for(std::size_t i=0;i<size;i++)
77 ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%size]));
83 QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes)
85 QuadraticPolygon *ret=new QuadraticPolygon;
86 std::size_t size=nodes.size();
87 for(std::size_t i=0;i<size/2;i++)
90 e1=new EdgeLin(nodes[i],nodes[i+size/2]);
91 e2=new EdgeLin(nodes[i+size/2],nodes[(i+1)%(size/2)]);
92 SegSegIntersector inters(*e1,*e2);
93 bool colinearity=inters.areColinears();
96 ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%(size/2)]));
98 ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+size/2],nodes[(i+1)%(size/2)]));
99 nodes[i]->decrRef(); nodes[i+size/2]->decrRef();
104 void QuadraticPolygon::BuildDbgFile(const std::vector<Node *>& nodes, const char *fileName)
106 std::ofstream file(fileName);
107 file << std::setprecision(16);
108 file << " double coords[]=" << std::endl << " { ";
109 for(std::vector<Node *>::const_iterator iter=nodes.begin();iter!=nodes.end();iter++)
111 if(iter!=nodes.begin())
112 file << "," << std::endl << " ";
113 file << (*(*iter))[0] << ", " << (*(*iter))[1];
115 file << "};" << std::endl;
118 void QuadraticPolygon::closeMe() const
120 if(!front()->changeStartNodeWith(back()->getEndNode()))
121 throw(Exception("big error: not closed polygon..."));
124 void QuadraticPolygon::circularPermute()
126 if(_sub_edges.size()>1)
128 ElementaryEdge *first=_sub_edges.front();
129 _sub_edges.pop_front();
130 _sub_edges.push_back(first);
134 bool QuadraticPolygon::isButterflyAbs()
136 INTERP_KERNEL::Bounds b;
138 b.prepareForAggregation();
140 double dimChar=b.getCaracteristicDim();
141 b.getBarycenter(xBary,yBary);
142 applyGlobalSimilarity(xBary,yBary,dimChar);
144 return isButterfly();
147 bool QuadraticPolygon::isButterfly() const
149 for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
151 Edge *e1=(*it)->getPtr();
152 std::list<ElementaryEdge *>::const_iterator it2=it;
154 for(;it2!=_sub_edges.end();it2++)
156 MergePoints commonNode;
157 ComposedEdge *outVal1=new ComposedEdge;
158 ComposedEdge *outVal2=new ComposedEdge;
159 Edge *e2=(*it2)->getPtr();
160 if(e1->intersectWith(e2,commonNode,*outVal1,*outVal2))
173 void QuadraticPolygon::dumpInXfigFileWithOther(const ComposedEdge& other, const char *fileName) const
175 std::ofstream file(fileName);
176 const int resolution=1200;
178 box.prepareForAggregation();
180 other.fillBounds(box);
181 dumpInXfigFile(file,resolution,box);
182 other.ComposedEdge::dumpInXfigFile(file,resolution,box);
185 void QuadraticPolygon::dumpInXfigFile(const char *fileName) const
187 std::ofstream file(fileName);
188 const int resolution=1200;
190 box.prepareForAggregation();
192 dumpInXfigFile(file,resolution,box);
195 void QuadraticPolygon::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const
197 stream << "#FIG 3.2 Produced by xfig version 3.2.5-alpha5" << std::endl;
198 stream << "Landscape" << std::endl;
199 stream << "Center" << std::endl;
200 stream << "Metric" << std::endl;
201 stream << "Letter" << std::endl;
202 stream << "100.00" << std::endl;
203 stream << "Single" << std::endl;
204 stream << "-2" << std::endl;
205 stream << resolution << " 2" << std::endl;
206 ComposedEdge::dumpInXfigFile(stream,resolution,box);
210 * Warning contrary to intersectWith method this method is \b NOT const. 'this' and 'other' are modified after call of this method.
212 double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other)
214 double ret=0.,xBaryBB,yBaryBB;
215 double fact=normalize(&other,xBaryBB,yBaryBB);
216 std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
217 for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
219 ret+=fabs((*iter)->getArea());
222 return ret*fact*fact;
226 * This method splits 'this' with 'other' into smaller pieces localizable. 'mapThis' is a map that gives the correspondance between nodes contained in 'this' and node ids in a global mesh.
227 * In the same way, 'mapOther' gives the correspondance between nodes contained in 'other' and node ids in a global mesh from wich 'other' is extracted.
228 * This method has 1 out paramater : 'edgesThis', After the call of this method contains nodal connectivity (including type) of 'this' into globlal "this mesh".
229 * This method has 2 in/out parameters : 'subDivOther' and 'addCoo'.'otherEdgeIds' is useful to put values in 'edgesThis', 'subDivOther' and 'addCoo'.
230 * Size of 'otherEdgeIds' has to be equal to number of ElementaryEdges in 'other'. No check of that will be done.
231 * @param offset1 is the number of nodes contained in global mesh from which 'this' is extracted.
232 * @param offset2 is the sum of nodes contained in global mesh from which 'this' is extracted and 'other' is extracted.
233 * @param edgesInOtherColinearWithThis will be appended at the end of the vector with colinear edge ids of other (if any)
234 * @otherEdgeIds is a vector with the same size than other before calling this method. It gives in the same order the cell id in global other mesh.
236 void QuadraticPolygon::splitAbs(QuadraticPolygon& other, const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2 , const std::vector<int>& otherEdgeIds,
237 std::vector<int>& edgesThis, int cellIdThis, std::vector< std::vector<int> >& edgesInOtherColinearWithThis, std::vector< std::vector<int> >& subDivOther, std::vector<double>& addCoo)
239 double xBaryBB, yBaryBB;
240 double fact=normalizeExt(&other, xBaryBB, yBaryBB);
242 IteratorOnComposedEdge it1(this),it3(&other);
244 ComposedEdge *c1=new ComposedEdge;
245 ComposedEdge *c2=new ComposedEdge;
247 std::map<INTERP_KERNEL::Node *,int> mapAddCoo;
248 for(it3.first();!it3.finished();it3.next(),i++)//iteration over 'other' _sub_edges
250 QuadraticPolygon otherTmp;
251 ElementaryEdge* curE3=it3.current();
252 otherTmp.pushBack(new ElementaryEdge(curE3->getPtr(),curE3->getDirection())); curE3->getPtr()->incrRef();
253 IteratorOnComposedEdge it2(&otherTmp);
254 for(it2.first();!it2.finished();it2.next())//iteration on subedges of 'other->_sub_edge'
256 ElementaryEdge* curE2=it2.current();
257 if(!curE2->isThereStartPoint())
260 it1=curE2->getIterator();
261 for(;!it1.finished();)//iteration over 'this' _sub_edges
263 ElementaryEdge* curE1=it1.current();
265 if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
267 if(!curE1->getDirection()) c1->reverse();
268 if(!curE2->getDirection()) c2->reverse();
269 UpdateNeighbours(merge,it1,it2,c1,c2);
270 //Substitution of simple edge by sub-edges.
271 delete curE1; // <-- destroying simple edge coming from pol1
272 delete curE2; // <-- destroying simple edge coming from pol2
273 it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next.
274 it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next.
277 it1.assignMySelfToAllElems(c2);//To avoid that others
285 UpdateNeighbours(merge,it1,it2,curE1,curE2);
290 if(otherTmp.presenceOfOn())
291 edgesInOtherColinearWithThis[otherEdgeIds[i]].push_back(cellIdThis);
292 if(otherTmp._sub_edges.size()>1)
294 for(std::list<ElementaryEdge *>::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++)
295 (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo);
301 for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
302 (*it)->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/edgesThis,addCoo,mapAddCoo);
307 * This method builds from descending conn of a quadratic polygon stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order to give the
308 * orientation of edge.
310 void QuadraticPolygon::buildFromCrudeDataArray(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords,
311 const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
313 std::size_t nbOfSeg=std::distance(descBg,descEnd);
314 for(std::size_t i=0;i<nbOfSeg;i++)
316 appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
320 void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords,
321 const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
325 bool direct=descBg[edgePos]>0;
326 int edgeId=abs(descBg[edgePos])-1;
327 const std::vector<int>& subEdge=intersectEdges[edgeId];
328 std::size_t nbOfSubEdges=subEdge.size()/2;
329 for(std::size_t j=0;j<nbOfSubEdges;j++)
330 appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
334 std::size_t nbOfSeg=std::distance(descBg,descEnd);
335 const double *st=coords+2*(nodalBg[edgePos]);
336 INTERP_KERNEL::Node *st0=new INTERP_KERNEL::Node(st[0],st[1]);
337 const double *endd=coords+2*(nodalBg[(edgePos+1)%nbOfSeg]);
338 INTERP_KERNEL::Node *endd0=new INTERP_KERNEL::Node(endd[0],endd[1]);
339 const double *middle=coords+2*(nodalBg[edgePos+nbOfSeg]);
340 INTERP_KERNEL::Node *middle0=new INTERP_KERNEL::Node(middle[0],middle[1]);
342 e1=new EdgeLin(st0,middle0);
343 e2=new EdgeLin(middle0,endd0);
344 SegSegIntersector inters(*e1,*e2);
345 bool colinearity=inters.areColinears();
346 delete e1; delete e2;
348 bool direct=descBg[edgePos]>0;
349 int edgeId=abs(descBg[edgePos])-1;
350 const std::vector<int>& subEdge=intersectEdges[edgeId];
351 std::size_t nbOfSubEdges=subEdge.size()/2;
354 for(std::size_t j=0;j<nbOfSubEdges;j++)
355 appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
359 Edge *e=new EdgeArcCircle(st0,middle0,endd0,true);
360 for(std::size_t j=0;j<nbOfSubEdges;j++)
361 appendSubEdgeFromCrudeDataArray(e,j,direct,edgeId,subEdge,mapp);
364 st0->decrRef(); endd0->decrRef(); middle0->decrRef();
368 void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, int edgeId, const std::vector<int>& subEdge, const std::map<int,INTERP_KERNEL::Node *>& mapp)
370 std::size_t nbOfSubEdges=subEdge.size()/2;
372 {//it is not a quadratic subedge
373 Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
374 Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
375 ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
379 {//it is a quadratic subedge
380 Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
381 Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
382 Edge *ee=baseEdge->buildEdgeLyingOnMe(start,end);
383 ElementaryEdge *eee=new ElementaryEdge(ee,true);
389 * This method builds from descending conn of a quadratic polygon stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order to give the
390 * orientation of edge.
392 void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords, const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges,
393 const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
394 const std::vector< std::vector<int> >& colinear1)
396 std::size_t nbOfSeg=std::distance(descBg,descEnd);
397 for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
399 bool direct=descBg[i]>0;
400 int edgeId=abs(descBg[i])-1;//current edge id of pol2
401 bool directos=colinear1[edgeId].empty();
402 std::vector<std::pair<int,std::pair<bool,int> > > idIns1;
405 {// if the current edge of pol2 has one or more colinear edges part into pol1
406 const std::vector<int>& c=colinear1[edgeId];
407 std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
408 for(std::size_t j=0;j<nbOfEdgesIn1;j++)
410 int edgeId1=abs(descBg1[j])-1;
411 if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
413 idIns1.push_back(std::pair<int,std::pair<bool,int> >(edgeId1,std::pair<bool,int>(descBg1[j]>0,offset1)));// it exists an edge into pol1 given by tuple (idIn1,direct1) that is colinear at edge 'edgeId' in pol2
414 //std::pair<edgeId1); direct1=descBg1[j]>0;
416 offset1+=intersectEdges1[edgeId1].size()/2;//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2
418 directos=idIns1.empty();
421 {//no subpart of edge 'edgeId' of pol2 is in pol1 so let's operate the same thing that QuadraticPolygon::buildFromCrudeDataArray method
422 appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
425 {//there is subpart of edge 'edgeId' of pol2 inside pol1
426 const std::vector<int>& subEdge=intersectEdges[edgeId];
427 std::size_t nbOfSubEdges=subEdge.size()/2;
428 for(std::size_t j=0;j<nbOfSubEdges;j++)
430 int idBg=direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1];
431 int idEnd=direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2];
432 bool direction11,found=false;
433 bool direct1;//store if needed the direction in 1
435 std::size_t nbOfSubEdges1;
436 for(std::vector<std::pair<int,std::pair<bool,int> > >::const_iterator it=idIns1.begin();it!=idIns1.end();it++)
438 int idIn1=(*it).first;//store if needed the cell id in 1
439 direct1=(*it).second.first;
440 offset1=(*it).second.second;
441 const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
442 nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
444 for(std::size_t k=0;k<nbOfSubEdges1 && !found;k++)
445 {//perform a loop on all subedges of pol1 that includes edge 'edgeId' of pol2. For the moment we iterate only on subedges of ['idIn1']... To improve
446 if(subEdge1PossiblyAlreadyIn1[2*k]==idBg && subEdge1PossiblyAlreadyIn1[2*k+1]==idEnd)
447 { direction11=true; found=true; }
448 else if(subEdge1PossiblyAlreadyIn1[2*k]==idEnd && subEdge1PossiblyAlreadyIn1[2*k+1]==idBg)
449 { direction11=false; found=true; }
455 {//the current subedge of edge 'edgeId' of pol2 is not a part of the colinear edge 'idIn1' of pol1 -> build new Edge instance
456 //appendEdgeFromCrudeDataArray(j,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
457 Node *start=(*mapp.find(idBg)).second;
458 Node *end=(*mapp.find(idEnd)).second;
459 ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
463 {//the current subedge of edge 'edgeId' of pol2 is part of the colinear edge 'idIn1' of pol1 -> reuse Edge instance of pol1
464 ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
465 Edge *ee=e->getPtr();
466 ee->incrRef(); ee->declareOn();
467 pushBack(new ElementaryEdge(ee,!(direct1^direction11)));
474 void QuadraticPolygon::appendCrudeData(const std::map<INTERP_KERNEL::Node *,int>& mapp, double xBary, double yBary, double fact, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI) const
477 bool presenceOfQuadratic=presenceOfQuadraticEdge();
478 conn.push_back(presenceOfQuadratic?NORM_QPOLYG:NORM_POLYGON);
479 for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
482 tmp=(*it)->getStartNode();
483 std::map<INTERP_KERNEL::Node *,int>::const_iterator it1=mapp.find(tmp);
484 conn.push_back((*it1).second);
487 if(presenceOfQuadratic)
490 int off=offset+((int)addCoordsQuadratic.size())/2;
491 for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,j++,nbOfNodesInPg++)
493 INTERP_KERNEL::Node *node=(*it)->getPtr()->buildRepresentantOfMySelf();
494 node->unApplySimilarity(xBary,yBary,fact);
495 addCoordsQuadratic.push_back((*node)[0]);
496 addCoordsQuadratic.push_back((*node)[1]);
497 conn.push_back(off+j);
501 connI.push_back(connI.back()+nbOfNodesInPg+1);
505 * This method make the hypothesis that 'this' and 'other' are splited at the minimum into edges that are fully IN, OUT or ON.
506 * This method returns newly created polygons in 'conn' and 'connI' and the corresponding ids ('idThis','idOther') are stored respectively into 'nbThis' and 'nbOther'.
508 void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, const std::map<INTERP_KERNEL::Node *,int>& mapp, int idThis, int idOther, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI, std::vector<int>& nbThis, std::vector<int>& nbOther)
510 double xBaryBB, yBaryBB;
511 double fact=normalizeExt(&other, xBaryBB, yBaryBB);
512 //Locate 'this' relative to 'other'
513 other.performLocatingOperation(*this);
514 std::vector<QuadraticPolygon *> res=buildIntersectionPolygons(other,*this);
515 for(std::vector<QuadraticPolygon *>::iterator it=res.begin();it!=res.end();it++)
517 (*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,offset,addCoordsQuadratic,conn,connI);
518 nbThis.push_back(idThis);
519 nbOther.push_back(idOther);
522 unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact);
526 * Warning This method is \b NOT const. 'this' and 'other' are modified after call of this method.
527 * 'other' is a QuadraticPolygon of \b non closed edges.
529 double QuadraticPolygon::intersectWithAbs1D(QuadraticPolygon& other, bool& isColinear)
531 double ret = 0., xBaryBB, yBaryBB;
532 double fact = normalize(&other, xBaryBB, yBaryBB);
534 QuadraticPolygon cpyOfThis(*this);
535 QuadraticPolygon cpyOfOther(other);
537 SplitPolygonsEachOther(cpyOfThis, cpyOfOther, nbOfSplits);
538 //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
539 performLocatingOperation(cpyOfOther);
541 for(std::list<ElementaryEdge *>::const_iterator it=cpyOfOther._sub_edges.begin();it!=cpyOfOther._sub_edges.end();it++)
543 switch((*it)->getLoc())
547 ret += fabs((*it)->getPtr()->getCurveLength());
553 ret += fabs((*it)->getPtr()->getCurveLength());
565 * Warning contrary to intersectWith method this method is \b NOT const. 'this' and 'other' are modified after call of this method.
567 double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other, double* barycenter)
569 double ret=0.,bary[2],area,xBaryBB,yBaryBB;
570 barycenter[0] = barycenter[1] = 0.;
571 double fact=normalize(&other,xBaryBB,yBaryBB);
572 std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
573 for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
575 area=fabs((*iter)->getArea());
576 (*iter)->getBarycenter(bary);
579 barycenter[0] += bary[0]*area;
580 barycenter[1] += bary[1]*area;
582 if ( ret > std::numeric_limits<double>::min() )
584 barycenter[0]=barycenter[0]/ret*fact+xBaryBB;
585 barycenter[1]=barycenter[1]/ret*fact+yBaryBB;
588 return ret*fact*fact;
592 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
593 * This is possible because loc attribute in Edge class is mutable.
594 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
596 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other) const
599 std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
600 for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
602 ret+=fabs((*iter)->getArea());
609 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
610 * This is possible because loc attribute in Edge class is mutable.
611 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
613 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other, double* barycenter) const
615 double ret=0., bary[2];
616 barycenter[0] = barycenter[1] = 0.;
617 std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
618 for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
620 double area = fabs((*iter)->getArea());
621 (*iter)->getBarycenter(bary);
624 barycenter[0] += bary[0]*area;
625 barycenter[1] += bary[1]*area;
627 if ( ret > std::numeric_limits<double>::min() )
629 barycenter[0] /= ret;
630 barycenter[1] /= ret;
636 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
637 * This is possible because loc attribute in Edge class is mutable.
638 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
640 void QuadraticPolygon::intersectForPerimeter(const QuadraticPolygon& other, double& perimeterThisPart, double& perimeterOtherPart, double& perimeterCommonPart) const
642 perimeterThisPart=0.; perimeterOtherPart=0.; perimeterCommonPart=0.;
643 QuadraticPolygon cpyOfThis(*this);
644 QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
645 SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
646 performLocatingOperation(cpyOfOther);
647 other.performLocatingOperation(cpyOfThis);
648 cpyOfThis.dispatchPerimeterExcl(perimeterThisPart,perimeterCommonPart);
649 cpyOfOther.dispatchPerimeterExcl(perimeterOtherPart,perimeterCommonPart);
650 perimeterCommonPart/=2.;
654 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
655 * This is possible because loc attribute in Edge class is mutable.
656 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
658 * polThis.size()==this->size() and polOther.size()==other.size().
659 * For each ElementaryEdge of 'this', the corresponding contribution in resulting polygon is in 'polThis'.
660 * For each ElementaryEdge of 'other', the corresponding contribution in resulting polygon is in 'polOther'.
661 * As consequence common part are counted twice (in polThis \b and in polOther).
663 void QuadraticPolygon::intersectForPerimeterAdvanced(const QuadraticPolygon& other, std::vector< double >& polThis, std::vector< double >& polOther) const
665 polThis.resize(size());
666 polOther.resize(other.size());
667 IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
669 for(it1.first();!it1.finished();it1.next(),edgeId++)
671 ElementaryEdge* curE1=it1.current();
672 QuadraticPolygon cpyOfOther(other);
673 QuadraticPolygon tmp;
674 tmp.pushBack(curE1->clone());
676 SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
677 other.performLocatingOperation(tmp);
678 tmp.dispatchPerimeter(polThis[edgeId]);
681 IteratorOnComposedEdge it2(const_cast<QuadraticPolygon *>(&other));
683 for(it2.first();!it2.finished();it2.next(),edgeId++)
685 ElementaryEdge* curE2=it2.current();
686 QuadraticPolygon cpyOfThis(*this);
687 QuadraticPolygon tmp;
688 tmp.pushBack(curE2->clone());
690 SplitPolygonsEachOther(tmp,cpyOfThis,tmp2);
691 performLocatingOperation(tmp);
692 tmp.dispatchPerimeter(polOther[edgeId]);
698 * numberOfCreatedPointsPerEdge is resized to the number of edges of 'this'.
699 * This method returns in ordered maner the number of newly created points per edge.
700 * This method performs a split process between 'this' and 'other' that gives the result PThis.
701 * Then for each edges of 'this' this method counts how many edges in Pthis have the same id.
703 void QuadraticPolygon::intersectForPoint(const QuadraticPolygon& other, std::vector< int >& numberOfCreatedPointsPerEdge) const
705 numberOfCreatedPointsPerEdge.resize(size());
706 IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
708 for(it1.first();!it1.finished();it1.next(),edgeId++)
710 ElementaryEdge* curE1=it1.current();
711 QuadraticPolygon cpyOfOther(other);
712 QuadraticPolygon tmp;
713 tmp.pushBack(curE1->clone());
715 SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
716 numberOfCreatedPointsPerEdge[edgeId]=tmp.recursiveSize()-1;
721 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
722 * This is possible because loc attribute in Edge class is mutable.
723 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
725 std::vector<QuadraticPolygon *> QuadraticPolygon::intersectMySelfWith(const QuadraticPolygon& other) const
727 QuadraticPolygon cpyOfThis(*this);
728 QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
729 SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
730 //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
731 performLocatingOperation(cpyOfOther);
732 return other.buildIntersectionPolygons(cpyOfThis,cpyOfOther);
736 * This method is typically the first step of boolean operations between pol1 and pol2.
737 * This method perform the minimal splitting so that at the end each edges constituting pol1 are fully either IN or OUT or ON.
738 * @param pol1 IN/OUT param that is equal to 'this' when called.
740 void QuadraticPolygon::SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits)
742 IteratorOnComposedEdge it1(&pol1),it2(&pol2);
744 ComposedEdge *c1=new ComposedEdge;
745 ComposedEdge *c2=new ComposedEdge;
746 for(it2.first();!it2.finished();it2.next())
748 ElementaryEdge* curE2=it2.current();
749 if(!curE2->isThereStartPoint())
752 it1=curE2->getIterator();
753 for(;!it1.finished();)
756 ElementaryEdge* curE1=it1.current();
757 merge.clear(); nbOfSplits++;
758 if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
760 if(!curE1->getDirection()) c1->reverse();
761 if(!curE2->getDirection()) c2->reverse();
762 UpdateNeighbours(merge,it1,it2,c1,c2);
763 //Substitution of simple edge by sub-edges.
764 delete curE1; // <-- destroying simple edge coming from pol1
765 delete curE2; // <-- destroying simple edge coming from pol2
766 it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next.
767 it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next.
770 it1.assignMySelfToAllElems(c2);//To avoid that others
778 UpdateNeighbours(merge,it1,it2,curE1,curE2);
787 void QuadraticPolygon::performLocatingOperation(QuadraticPolygon& pol2) const
789 IteratorOnComposedEdge it(&pol2);
790 TypeOfEdgeLocInPolygon loc=FULL_ON_1;
791 for(it.first();!it.finished();it.next())
793 ElementaryEdge *cur=it.current();
794 loc=cur->locateFullyMySelf(*this,loc);
799 * Given 2 polygons 'pol1' and 'pol2' (localized) the resulting polygons are returned.
801 * this : pol2 simplified.
802 * @param pol1 pol1 split.
803 * @param pol2 pol2 split.
805 std::vector<QuadraticPolygon *> QuadraticPolygon::buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const
807 std::vector<QuadraticPolygon *> ret;
808 std::list<QuadraticPolygon *> pol2Zip=pol2.zipConsecutiveInSegments();
810 closePolygons(pol2Zip,pol1,ret);
812 {//borders of pol2 do not cross pol1,and pol2 borders are outside of pol1. That is to say, either pol2 and pol1
813 //do not overlap or pol1 is fully inside pol2. So in the first case no intersection, in the other case
814 //the intersection is pol1.
815 ElementaryEdge *e1FromPol1=pol1[0];
816 TypeOfEdgeLocInPolygon loc=FULL_ON_1;
817 loc=e1FromPol1->locateFullyMySelf(*this,loc);
819 ret.push_back(new QuadraticPolygon(pol1));
825 * Returns parts of potentially non closed-polygons. Each returned polygons are not mergeable.
826 * this : pol2 split and locallized.
828 std::list<QuadraticPolygon *> QuadraticPolygon::zipConsecutiveInSegments() const
830 std::list<QuadraticPolygon *> ret;
831 IteratorOnComposedEdge it((ComposedEdge *)this);
832 int nbOfTurns=recursiveSize();
834 if(!it.goToNextInOn(false,i,nbOfTurns))
840 QuadraticPolygon *tmp1=new QuadraticPolygon;
841 TypeOfEdgeLocInPolygon loc=it.current()->getLoc();
842 while(loc!=FULL_OUT_1 && i<nbOfTurns)
844 ElementaryEdge *tmp3=it.current()->clone();
845 tmp1->pushBack(tmp3);
847 loc=it.current()->getLoc();
855 it.goToNextInOn(true,i,nbOfTurns);
861 * 'this' should be considered as pol2Simplified.
862 * @param pol2zip is a list of set of edges (openned polygon) coming from split polygon 2.
863 * @param pol1 is split pol1.
864 * @param results the resulting \b CLOSED polygons.
866 void QuadraticPolygon::closePolygons(std::list<QuadraticPolygon *>& pol2Zip, const QuadraticPolygon& pol1,
867 std::vector<QuadraticPolygon *>& results) const
869 bool directionKnownInPol1=false;
870 bool directionInPol1;
871 for(std::list<QuadraticPolygon *>::iterator iter=pol2Zip.begin();iter!=pol2Zip.end();)
873 if((*iter)->completed())
875 results.push_back(*iter);
876 directionKnownInPol1=false;
877 iter=pol2Zip.erase(iter);
880 if(!directionKnownInPol1)
882 if(!(*iter)->amIAChanceToBeCompletedBy(pol1,*this,directionInPol1))
883 { delete *iter; iter=pol2Zip.erase(iter); continue; }
885 directionKnownInPol1=true;
887 std::list<QuadraticPolygon *>::iterator iter2=iter; iter2++;
888 std::list<QuadraticPolygon *>::iterator iter3=(*iter)->fillAsMuchAsPossibleWith(pol1,iter2,pol2Zip.end(),directionInPol1);
889 if(iter3!=pol2Zip.end())
891 (*iter)->pushBack(*iter3);
893 pol2Zip.erase(iter3);
899 * 'this' is expected to be set of edges (not closed) of pol2 split.
901 bool QuadraticPolygon::amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted,const QuadraticPolygon& pol2NotSplitted, bool& direction)
903 IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1Splitted));
905 Node *n=getEndNode();
906 ElementaryEdge *cur=it.current();
907 for(it.first();!it.finished() && !found;)
910 found=(cur->getStartNode()==n);
915 throw Exception("Internal error : polygons uncompatible each others. Should never happend");
916 //Ok we found correspondance between this and pol1. Searching for right direction to close polygon.
917 ElementaryEdge *e=_sub_edges.back();
918 if(e->getLoc()==FULL_ON_1)
920 if(e->getPtr()==cur->getPtr())
925 Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
926 bool ret=pol2NotSplitted.isInOrOut(repr);
933 Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
934 bool ret=pol2NotSplitted.isInOrOut(repr);
940 direction=cur->locateFullyMySelfAbsolute(pol2NotSplitted)==FULL_IN_1;
945 * This method fills as much as possible 'this' (part of pol2 split) with edges of 'pol1Splitted'.
947 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted,
948 std::list<QuadraticPolygon *>::iterator iStart,
949 std::list<QuadraticPolygon *>::iterator iEnd,
952 IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1Splitted));
954 Node *n=getEndNode();
956 for(it.first();!it.finished() && !found;)
959 found=(cur->getStartNode()==n);
966 std::list<QuadraticPolygon *>::iterator ret;
970 ElementaryEdge *tmp=cur->clone();
974 nodeToTest=tmp->getEndNode();
975 direction?it.nextLoop():it.previousLoop();
976 ret=CheckInList(nodeToTest,iStart,iEnd);
984 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::CheckInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
985 std::list<QuadraticPolygon *>::iterator iEnd)
987 for(std::list<QuadraticPolygon *>::iterator iter=iStart;iter!=iEnd;iter++)
988 if((*iter)->isNodeIn(n))