![]() |
Tulip
4.6.0
Better Visualization Through Research
|
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 }