Tulip  4.6.0
Better Visualization Through Research
library/tulip-core/include/tulip/cxx/Circle.cxx
00001 /*
00002  *
00003  * This file is part of Tulip (www.tulip-software.org)
00004  *
00005  * Authors: David Auber and the Tulip development Team
00006  * from LaBRI, University of Bordeaux
00007  *
00008  * Tulip is free software; you can redistribute it and/or modify
00009  * it under the terms of the GNU Lesser General Public License
00010  * as published by the Free Software Foundation, either version 3
00011  * of the License, or (at your option) any later version.
00012  *
00013  * Tulip is distributed in the hope that it will be useful,
00014  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00015  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00016  * See the GNU General Public License for more details.
00017  *
00018  */
00019 #include <cstdlib>
00020 
00021 template<typename Obj, typename OTYPE>
00022 tlp::Circle<Obj, OTYPE> & tlp::Circle<Obj, OTYPE>::merge(const tlp::Circle<Obj, OTYPE> &c) {
00023   Vector<Obj,2, OTYPE> p1(*this);
00024   Vector<Obj,2, OTYPE> p2(c);
00025   Vector<Obj,2, OTYPE> c12=p2-p1;
00026   double norm=c12.norm();
00027 
00028   if (norm < 0.0000001) {
00029     if (radius<c.radius)
00030       return (*this)=c;
00031     else
00032       return *this;
00033   }
00034 
00035   double dx=c12[0]/norm;
00036   double dy=c12[1]/norm;
00037   p1[0]-=dx*radius;
00038   p1[1]-=dy*radius;
00039   p2[0]+=dx*c.radius;
00040   p2[1]+=dy*c.radius;
00041   Obj tmpDist=p1.dist(p2)/2.;
00042 
00043   if ( (tmpDist<radius) || (tmpDist<c.radius)) {
00044     if (radius>c.radius)
00045       return *this;
00046     else
00047       return (*this)=c;
00048   }
00049   else {
00050     p1+=p2;
00051     p1/=2.;
00052     (*this)[0]=p1[0];
00053     (*this)[1]=p1[1];
00054     radius=tmpDist;
00055     return *this;
00056   }
00057 }
00058 
00059 template<typename Obj, typename OTYPE>
00060 bool tlp::Circle<Obj, OTYPE>::isIncludeIn(const tlp::Circle<Obj, OTYPE> &c) const {
00061   Vector<Obj,2, OTYPE> dir=c-*this;
00062   return (dir.norm()+radius)<=c.radius;
00063 }
00064 
00065 template<typename Obj, typename OTYPE>
00066 tlp::Circle<Obj, OTYPE> tlp::enclosingCircle(const tlp::Circle<Obj, OTYPE> &c1,const tlp::Circle<Obj, OTYPE> &c2) {
00067   Vector<Obj,2, OTYPE> dir=c2-c1;
00068   Obj n=dir.norm();
00069 
00070   if (n==0)
00071     return Circle<Obj, OTYPE>(c1, std::max(c1.radius,  c2.radius));
00072 
00073   dir/=n;
00074   Vector<Obj, 2, OTYPE> ext1=c1-dir*c1.radius;
00075   Vector<Obj, 2, OTYPE> ext2=c2+dir*c2.radius;
00076   return Circle<Obj, OTYPE>((ext1+ext2)/Obj(2),(ext2-ext1).norm()/Obj(2));
00077 }
00078 
00079 template<typename Obj, typename OTYPE>
00080 tlp::Circle<Obj, OTYPE> tlp::lazyEnclosingCircle(const std::vector<tlp::Circle<Obj, OTYPE> > & circles) {
00081   //compute bounding box of a
00082   tlp::Vector<Obj,4, OTYPE> boundingBox;
00083   //  for (int i=0;i<4;++i) boundingBox[i]=0;
00084   typename std::vector< tlp::Circle<Obj, OTYPE> >::const_iterator it=circles.begin();
00085   boundingBox[0]=(*it)[0]-(*it).radius;
00086   boundingBox[1]=(*it)[1]-(*it).radius;
00087   boundingBox[2]=(*it)[0]+(*it).radius;
00088   boundingBox[3]=(*it)[1]+(*it).radius;
00089   ++it;
00090 
00091   for (; it!=circles.end(); ++it) {
00092     boundingBox[0] = std::min(boundingBox[0] , ((*it)[0]-(*it).radius));
00093     boundingBox[1] = std::min(boundingBox[1] , ((*it)[1]-(*it).radius));
00094     boundingBox[2] = std::max(boundingBox[2] , ((*it)[0]+(*it).radius));
00095     boundingBox[3] = std::max(boundingBox[3] , ((*it)[1]+(*it).radius));
00096   }
00097 
00098   tlp::Vector<Obj,2, OTYPE> center;
00099   center[0]=(boundingBox[0]+boundingBox[2])/2.;
00100   center[1]=(boundingBox[1]+boundingBox[3])/2.;
00101   Obj radius = std::max( (boundingBox[2]-boundingBox[0])/2., (boundingBox[3]-boundingBox[1])/2. );
00102   tlp::Circle<Obj, OTYPE> result(center,radius);
00103 
00104   //compute circle hull
00105   for (typename std::vector< tlp::Circle<Obj, OTYPE> >::const_iterator it=circles.begin(); it!=circles.end(); ++it)
00106     result.merge(*it);
00107 
00108   return result;
00109 }
00110 
00111 template<typename Obj, typename OTYPE>
00112 tlp::Circle<Obj, OTYPE> tlp::enclosingCircle(const std::vector<tlp::Circle<Obj, OTYPE> > & circles) {
00113   class OptimumCircleHull {
00114     const std::vector<tlp::Circle<Obj, OTYPE> > *circles;
00115     std::vector<unsigned> enclosedCircles;
00116     unsigned first,last;
00117     unsigned b1,b2;
00118     tlp::Circle<Obj, OTYPE> result;
00119 
00120     static tlp::Circle<Obj, OTYPE> enclosingCircle(const tlp::Circle<Obj, OTYPE> &c1,const tlp::Circle<Obj, OTYPE> &c2,const tlp::Circle<Obj, OTYPE> &c3) {
00121       Obj a1 = c1[0];
00122       Obj b1 = c1[1];
00123       Obj r1 = c1.radius;
00124       Obj a2 = c2[0];
00125       Obj b2 = c2[1];
00126       Obj r2 = c2.radius;
00127       Obj a3 = c3[0];
00128       Obj b3 = c3[1];
00129       Obj r3 = c3.radius;
00130       Obj b = -b1*r2*a1*a1*b3+a3*a3*r2*r2*r2+b2*b2*r3*r3*r3+b1*b1*r3*r3*r3+b1*b1*r2*r2*r2+r1*r1*r1*b3*b3+r1*r1*r1*b2*b2+r2*r2*r2*b3*b3
00131               +a1*a1*r3*r3*r3+a2*a2*r1*r1*r1+a2*a2*r3*r3*r3+a1*a1*r2*r2*r2+a3*a3*r1*r1*r1-b2*b2*r3*a3*a3-b2*b2*r3*b3*b3-2*b2*r3*r3*r3*b1
00132               -b1*b1*r3*a3*a3-b1*b1*r3*b3*b3-b1*b1*r2*a2*a2+2*b1*b1*r2*b3*b3-b1*b1*r2*r3*r3-r1*b3*b3*b3*b2+r1*b3*b3*b3*b1
00133               +2*r1*b2*b2*b3*b3-r1*b2*b2*r3*r3+r2*b3*b3*b3*b2-r2*b3*b3*b3*b1-b2*b2*b2*r3*b1+2*b2*b2*r3*b1*b1+b2*b2*b2*r3*b3
00134               -b2*b2*r3*r1*r1+b1*b1*b1*r3*b3-b1*b1*b1*r3*b2-b1*b1*r3*r2*r2-b1*b1*r2*b2*b2-b1*b1*b1*r2*b3+b1*b1*b1*r2*b2
00135               -2*b1*r2*r2*r2*b3-r1*b3*b3*b1*b1-2*r1*r1*r1*b3*b2-r1*b3*b3*r2*r2-r1*b3*b3*a1*a1+r1*b2*b2*b2*b1
00136               -r1*b2*b2*b1*b1-r1*b2*b2*b2*b3-r1*b2*b2*a1*a1-r2*b3*b3*a2*a2-r2*b3*b3*b2*b2-r2*b3*b3*r1*r1
00137               -b2*r3*b1*b1*b3+b2*r3*a2*a2*b3+b2*r3*r1*r1*b3-b2*r3*r2*r2*b3+b2*r3*a1*a1*b3-2*b2*r3*a1*a2*b3+a1*a3*b2*b2*r3
00138               -2*a1*a3*b2*b1*r3-2*a1*a3*b2*b1*r2-2*a1*a3*b2*r1*b3+a1*a3*b2*b2*r1-2*a1*a3*b2*r2*b3-b2*r3*b1*a2*a2
00139               +2*b2*r3*b1*a3*a3+2*b2*r3*b1*b3*b3+b1*r2*b2*a3*a3-b1*r2*b2*b3*b3+b1*r2*b2*r3*r3+r1*b3*b1*a2*a2
00140               -r1*b3*b2*a3*a3+r1*b3*b2*r3*r3+r1*b3*b1*a3*a3-r1*b3*b1*r3*r3+r1*b2*b1*a2*a2+r1*b2*b1*a3*a3
00141               -r1*b2*b1*b3*b3+r1*b2*b1*r3*r3+2*r2*b3*b1*a2*a2+r2*b3*b2*a3*a3-r2*b3*b2*r3*r3-r2*b3*b1*a3*a3
00142               +r2*b3*b1*r3*r3+b2*r3*b1*r2*r2+4*b2*r3*a1*a2*b1+b1*r3*a2*a2*b3-b1*r3*b2*b2*b3
00143               -b1*r3*r1*r1*b3+b1*r3*r1*r1*b2+b1*r3*r2*r2*b3+b1*r3*a1*a1*b3-b1*r3*a1*a1*b2-2*b1*r3*a1*a2*b3-b1*b1*r3*a1*a2
00144               +b1*b1*r3*a1*a3+2*b1*r2*b2*b2*b3+b1*r2*r1*r1*b3-b1*r2*r1*r1*b2+b1*r2*a1*a1*b2
00145               -2*b1*r2*a1*a2*b3+b1*b1*r2*a1*a2-b1*b1*r2*a1*a3-r1*b3*b1*b2*b2+2*r1*b3*b1*b1*b2+2*r1*b3*a1*a1*b2+r1*b3*b1*r2*r2
00146               +r1*b3*b3*a1*a2-r1*b2*a2*a2*b3+r1*b2*r2*r2*b3-r1*b2*b1*r2*r2-2*r1*b2*a1*a2*b3-r2*b3*b1*b1*b2
00147               +r2*b3*r1*r1*b2+r2*b3*a1*a1*b2+r2*b3*b3*a1*a2+4*r2*b3*a1*a3*b1-a1*a1*a3*a3*r3+2*a1*a1*a3*a3*r2+a1*a3*a3*a3*r1
00148               -a1*a1*b3*b3*r3-a1*b3*b3*a3*r2+a1*r3*r3*a3*r2-a2*a1*a1*a3*r2+a2*b1*b1*a3*r2+a2*b3*b3*a3*r2
00149               -a1*a3*a3*a2*r1+2*a1*a3*a3*a2*r3+2*a1*b3*b3*a2*r3+a1*b3*b3*a3*r1-2*a1*r3*r3*r3*a2-a1*a1*r3*r3*r2
00150               -a2*a1*a1*a1*r3-a2*a2*a1*a1*r1+2*a2*a2*a1*a1*r3+a2*a1*a1*a1*r2-a2*a2*b1*b1*r1+2*a2*a2*a3*a3*r1-a2*a2*a3*a3*r3
00151               -a2*a3*a3*a3*r1-a2*a2*b3*b3*r3-a2*a2*r1*r1*r3-2*a2*r1*r1*r1*a3+a1*r3*r3*a2*r1-a1*r3*r3*a3*r1
00152               +2*a2*a1*a1*a3*r1+2*a2*b1*b1*a3*r1-a2*a3*a3*a1*r2-a2*b3*b3*a3*r1+a2*r1*r1*a1*r3-a2*r1*r1*a1*r2-a2*a2*r3*r3*r1
00153               -a1*a3*a3*a3*r2+a2*a3*a3*a3*r2-a1*a1*r3*r2*r2+a1*a1*a1*r3*a3+a2*a2*a2*r1*a1-a2*a2*a2*r1*a3-a2*a2*a2*r3*a1+a2*a2*a2*r3*a3
00154               -a1*a1*r2*a2*a2-a1*a1*r2*b2*b2-a1*a1*a1*r2*a3-2*a1*r2*r2*r2*a3-a3*a3*r1*a1*a1-a3*a3*r1*b1*b1-a3*a3*r1*r2*r2
00155               -a3*a3*r2*a2*a2-a3*a3*r2*b2*b2-a3*a3*r2*r1*r1+a2*r3*r3*a1*r2+a2*r3*r3*a3*r1-2*a3*r2*b1*a2*b3
00156               +a2*r1*r1*a3*r2-a2*r3*r3*a3*r2-a1*r3*a3*a2*a2-a1*r3*a3*r1*r1+a1*r3*a3*r2*r2+a2*r1*a1*b2*b2
00157               -a2*r1*a1*r2*r2-a2*r1*a3*b2*b2+a2*r1*a3*r2*r2-2*a2*r1*b1*a3*b2-a2*r3*a1*b2*b2+a2*r3*a1*r2*r2
00158               -a2*r3*a3*a1*a1+a2*r3*a3*b1*b1+a2*r3*a3*b2*b2+a2*r3*a3*r1*r1-a2*r3*a3*r2*r2-2*a2*r3*b1*a3*b2+2*a1*r2*a3*a2*a2
00159               +2*a1*r2*a3*b2*b2+a1*r2*a3*r1*r1-a3*r1*a1*a2*a2+a3*r1*a1*r2*r2-2*a3*r1*b1*a2*b3
00160               +4*r1*a3*b2*a2*b3;
00161       Obj tmp=a2*b3-a3*b2-a2*b1-a1*b3+a1*b2+a3*b1;
00162       Obj d = sqrt((a2*a2+b2*b2-r2*r2-2*a2*a3+a3*a3-r3*r3-2*b3*b2+b3*b3+2*r3*r2)*
00163                    (b3*b3+b1*b1-r1*r1-r3*r3-2*b1*b3+2*r3*r1+a3*a3-2*a1*a3+a1*a1)*
00164                    (a1*a1-2*a1*a2-r1*r1+b1*b1+b2*b2-r2*r2-2*b2*b1+2*r2*r1+a2*a2)*
00165                    tmp*tmp);
00166       Obj v = d-b;
00167 
00168       if (v<0) return tlp::Circle<Obj, OTYPE>(0,0,0);
00169 
00170       Obj aa = -2*a3*b2*b2*a1-2*a2*b3*b3*a1
00171                -2*a1*a1*b3*b2+a3*a3*b2*b2+a2*a2*b3*b3-r1*r1*b3*b3-r1*r1*b2*b2+a2*a2*b1*b1+a3*a3*b1*b1-a2*a2*r1*r1
00172                -a2*a2*r3*r3-a3*a3*r2*r2-2*a3*b2*a2*b3+2*a3*b2*a2*b1+2*a2*b3*a3*b1-2*b1*a2*a2*b3-a3*a3*r1*r1
00173                -2*a3*a3*b2*b1+2*b1*r3*r3*b2+2*r1*r1*b3*b2+2*b3*b1*r2*r2-2*a2*b1*b1*a3+2*a2*a3*r1*r1-b2*b2*r3*r3-b1*b1*r3*r3
00174                -b1*b1*r2*r2-r2*r2*b3*b3-a1*a1*r3*r3-a1*a1*r2*r2+a1*a1*b3*b3+a1*a1*b2*b2+2*b2*b2*r3*r1+2*b1*b1*r3*r2
00175                +2*r1*b3*b3*r2-2*b2*r3*b1*r2-2*b2*r3*r1*b3+2*b2*r3*r2*b3+2*b1*r3*r1*b3-2*b1*r3*r1*b2
00176                -2*b1*r3*r2*b3-2*b1*r2*r1*b3+2*b1*r2*r1*b2-2*r1*b2*r2*b3+2*a1*r3*r3*a2+2*a1*a1*r3*r2+2*a2*a2*r1*r3
00177                +2*a1*r2*r2*a3+2*a3*a3*r1*r2-2*a1*r3*a2*r1+2*a1*r3*a3*r1+2*a2*r1*a1*r2-2*a2*r3*a1*r2
00178                -2*a2*r3*a3*r1-2*a1*r2*a3*r1-2*a1*r3*a3*r2-2*a2*r1*a3*r2+2*a2*r3*a3*r2+2*a2*b3*a1*b2
00179                +2*a3*b2*a1*b3+2*a2*b1*a1*b3-2*a2*b1*a1*b2-2*a1*b3*a3*b1+2*a1*b2*a3*b1;
00180       Obj R = 0.5*v/aa;
00181       Obj y = -0.5*(-2*a1*R*r2-2*a3*R*r1+2*a2*R*r1+2*a3*R*r2-2*a2*R*r3+a1*a3*a3+a1*b3*b3-a1*r3*r3
00182                     +a2*a1*a1+a2*b1*b1+a1*r2*r2+2*a1*R*r3-a3*b1*b1+a3*a2*a2+a3*b2*b2+a3*r1*r1-a3*r2*r2-a2*a3*a3-a2*b3*b3
00183                     -a2*r1*r1+a2*r3*r3-a1*a2*a2-a1*b2*b2-a3*a1*a1)/(a2*b3+a3*b1-a2*b1-a1*b3+a1*b2-a3*b2);
00184       Obj x = 0.5*(-a1*a1*b3+a1*a1*b2+2*R*r2*b3+b1*a3*a3+b1*b3*b3+2*b1*R*r3-2*R*r1*b3+2*R*r1*b2+b1*r2*r2
00185                    -b2*a3*a3-b2*b3*b3+b2*r3*r3-2*b2*R*r3-b1*r3*r3-r2*r2*b3+a2*a2*b3-r1*r1*b2-b1*b2*b2+b1*b1*b2-b1*b1*b3
00186                    -b1*a2*a2+b2*b2*b3+r1*r1*b3-2*b1*R*r2)/(a2*b3+a3*b1-a2*b1-a1*b3+a1*b2-a3*b2);
00187       return  tlp::Circle<Obj, OTYPE>(x,y,R);
00188     }
00189 
00190     void process2() {
00191       if (isEmpty()) {
00192         result=tlp::enclosingCircle((*circles)[b1],(*circles)[b2]);
00193       }
00194       else {
00195         unsigned selectedCircle=popBack();
00196         process2();
00197 
00198         if (!(*circles)[selectedCircle].isIncludeIn(result)) {
00199           result=enclosingCircle((*circles)[b1],(*circles)[b2],(*circles)[selectedCircle]);
00200           pushFront(selectedCircle);
00201         }
00202         else {
00203           pushBack(selectedCircle);
00204         }
00205       }
00206     }
00207 
00208     void process1() {
00209       if (isEmpty()) {
00210         result=(*circles)[b1];
00211       }
00212       else {
00213         unsigned selectedCircle=popBack();
00214         process1();
00215 
00216         if (!(*circles)[selectedCircle].isIncludeIn(result)) {
00217           b2=selectedCircle;
00218           process2();
00219           pushFront(selectedCircle);
00220         }
00221         else {
00222           pushBack(selectedCircle);
00223         }
00224       }
00225     }
00226     void process0() {
00227       if (isEmpty()) {
00228         result=Circle<Obj, OTYPE>(0,0,0);
00229       }
00230       else {
00231         unsigned selectedCircle=popBack();
00232         process0();
00233 
00234         if (!(*circles)[selectedCircle].isIncludeIn(result)) {
00235           b1=selectedCircle;
00236           process1();
00237           pushFront(selectedCircle);
00238         }
00239         else {
00240           pushBack(selectedCircle);
00241         }
00242       }
00243     }
00244     bool isEmpty() const {
00245       return first==(last+1)%enclosedCircles.size();
00246     }
00247     void pushFront(unsigned c) {
00248       first=(first+enclosedCircles.size()-1)%enclosedCircles.size();
00249       enclosedCircles[first]=c;
00250     }
00251     unsigned popBack() {
00252       unsigned result=enclosedCircles[last];
00253       last=(last+enclosedCircles.size()-1)%enclosedCircles.size();
00254       return result;
00255     }
00256     void pushBack(unsigned c) {
00257       last=(last+1)%enclosedCircles.size();
00258       enclosedCircles[last]=c;
00259     }
00260   public:
00261     OptimumCircleHull():circles(NULL),first(0),last(0),b1(0),b2(0) {}
00262 
00263     tlp::Circle<Obj, OTYPE> operator()(const std::vector<tlp::Circle<Obj, OTYPE> > &circlesSet) {
00264       circles=&circlesSet;
00265       enclosedCircles.resize(circlesSet.size()+1);
00266       first=0;
00267       last=circlesSet.size()-1;
00268 
00269       for (unsigned i=0; i<circlesSet.size(); ++i)
00270         enclosedCircles[i]=i;
00271 
00272       for (unsigned i=circlesSet.size(); i>0;) {
00273         unsigned idx=(unsigned)(1.0*rand()*i/RAND_MAX);
00274         --i;
00275         std::swap(enclosedCircles[idx],enclosedCircles[i]);
00276       }
00277 
00278       process0();
00279       return result;
00280     }
00281   };
00282   return OptimumCircleHull()(circles);
00283 }
00284 
00285 template <typename Obj, typename OTYPE>
00286 std::ostream& tlp::operator<<(std::ostream &os,const tlp::Circle<Obj, OTYPE> &a) {
00287   os << "((" << a[0] << "," <<  a[1] << ")," << a.radius << ")";
00288   return os;
00289 }
 All Classes Files Functions Variables Enumerations Enumerator Properties