Tulip  4.1.0
Better Visualization Through Research
 All Classes Files Functions Variables Enumerations Enumerator Properties Groups Pages
Circle.cxx
1 /*
2  *
3  * This file is part of Tulip (www.tulip-software.org)
4  *
5  * Authors: David Auber and the Tulip development Team
6  * from LaBRI, University of Bordeaux 1 and Inria Bordeaux - Sud Ouest
7  *
8  * Tulip is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU Lesser General Public License
10  * as published by the Free Software Foundation, either version 3
11  * of the License, or (at your option) any later version.
12  *
13  * Tulip is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16  * See the GNU General Public License for more details.
17  *
18  */
19 #include <cstdlib>
20 
21 template<typename Obj, typename OTYPE>
22 tlp::Circle<Obj, OTYPE> & tlp::Circle<Obj, OTYPE>::merge(const tlp::Circle<Obj, OTYPE> &c) {
23  Vector<Obj,2, OTYPE> p1(*this);
24  Vector<Obj,2, OTYPE> p2(c);
25  Vector<Obj,2, OTYPE> c12=p2-p1;
26  double norm=c12.norm();
27 
28  if (norm < 0.0000001) {
29  if (radius<c.radius)
30  return (*this)=c;
31  else
32  return *this;
33  }
34 
35  double dx=c12[0]/norm;
36  double dy=c12[1]/norm;
37  p1[0]-=dx*radius;
38  p1[1]-=dy*radius;
39  p2[0]+=dx*c.radius;
40  p2[1]+=dy*c.radius;
41  Obj tmpDist=p1.dist(p2)/2.;
42 
43  if ( (tmpDist<radius) || (tmpDist<c.radius)) {
44  if (radius>c.radius)
45  return *this;
46  else
47  return (*this)=c;
48  }
49  else {
50  p1+=p2;
51  p1/=2.;
52  (*this)[0]=p1[0];
53  (*this)[1]=p1[1];
54  radius=tmpDist;
55  return *this;
56  }
57 }
58 
59 template<typename Obj, typename OTYPE>
60 bool tlp::Circle<Obj, OTYPE>::isIncludeIn(const tlp::Circle<Obj, OTYPE> &c) const {
61  Vector<Obj,2, OTYPE> dir=c-*this;
62  return (dir.norm()+radius)<=c.radius;
63 }
64 
65 template<typename Obj, typename OTYPE>
66 tlp::Circle<Obj, OTYPE> tlp::enclosingCircle(const tlp::Circle<Obj, OTYPE> &c1,const tlp::Circle<Obj, OTYPE> &c2) {
67  Vector<Obj,2, OTYPE> dir=c2-c1;
68  Obj n=dir.norm();
69 
70  if (n==0)
71  return Circle<Obj, OTYPE>(c1, std::max(c1.radius, c2.radius));
72 
73  dir/=n;
74  Vector<Obj, 2, OTYPE> ext1=c1-dir*c1.radius;
75  Vector<Obj, 2, OTYPE> ext2=c2+dir*c2.radius;
76  return Circle<Obj, OTYPE>((ext1+ext2)/Obj(2),(ext2-ext1).norm()/Obj(2));
77 }
78 
79 template<typename Obj, typename OTYPE>
80 tlp::Circle<Obj, OTYPE> tlp::lazyEnclosingCircle(const std::vector<tlp::Circle<Obj, OTYPE> > & circles) {
81  //compute bounding box of a
82  tlp::Vector<Obj,4, OTYPE> boundingBox;
83  // for (int i=0;i<4;++i) boundingBox[i]=0;
84  typename std::vector< tlp::Circle<Obj, OTYPE> >::const_iterator it=circles.begin();
85  boundingBox[0]=(*it)[0]-(*it).radius;
86  boundingBox[1]=(*it)[1]-(*it).radius;
87  boundingBox[2]=(*it)[0]+(*it).radius;
88  boundingBox[3]=(*it)[1]+(*it).radius;
89  ++it;
90 
91  for (; it!=circles.end(); ++it) {
92  boundingBox[0] = std::min(boundingBox[0] , ((*it)[0]-(*it).radius));
93  boundingBox[1] = std::min(boundingBox[1] , ((*it)[1]-(*it).radius));
94  boundingBox[2] = std::max(boundingBox[2] , ((*it)[0]+(*it).radius));
95  boundingBox[3] = std::max(boundingBox[3] , ((*it)[1]+(*it).radius));
96  }
97 
98  tlp::Vector<Obj,2, OTYPE> center;
99  center[0]=(boundingBox[0]+boundingBox[2])/2.;
100  center[1]=(boundingBox[1]+boundingBox[3])/2.;
101  Obj radius = std::max( (boundingBox[2]-boundingBox[0])/2., (boundingBox[3]-boundingBox[1])/2. );
102  tlp::Circle<Obj, OTYPE> result(center,radius);
103 
104  //compute circle hull
105  for (typename std::vector< tlp::Circle<Obj, OTYPE> >::const_iterator it=circles.begin(); it!=circles.end(); ++it)
106  result.merge(*it);
107 
108  return result;
109 }
110 
111 template<typename Obj, typename OTYPE>
112 tlp::Circle<Obj, OTYPE> tlp::enclosingCircle(const std::vector<tlp::Circle<Obj, OTYPE> > & circles) {
113  class OptimumCircleHull {
114  const std::vector<tlp::Circle<Obj, OTYPE> > *circles;
115  std::vector<unsigned> enclosedCircles;
116  unsigned first,last;
117  unsigned b1,b2;
118  tlp::Circle<Obj, OTYPE> result;
119 
120  static tlp::Circle<Obj, OTYPE> enclosingCircle(const tlp::Circle<Obj, OTYPE> &c1,const tlp::Circle<Obj, OTYPE> &c2,const tlp::Circle<Obj, OTYPE> &c3) {
121  Obj a1 = c1[0];
122  Obj b1 = c1[1];
123  Obj r1 = c1.radius;
124  Obj a2 = c2[0];
125  Obj b2 = c2[1];
126  Obj r2 = c2.radius;
127  Obj a3 = c3[0];
128  Obj b3 = c3[1];
129  Obj r3 = c3.radius;
130  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
131  +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
132  -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
133  +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
134  -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
135  -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
136  -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
137  -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
138  -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
139  +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
140  -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
141  -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
142  +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
143  -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
144  +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
145  -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
146  +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
147  +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
148  -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
149  -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
150  -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
151  -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
152  +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
153  -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
154  -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
155  -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
156  +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
157  -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
158  -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
159  +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
160  +4*r1*a3*b2*a2*b3;
161  Obj tmp=a2*b3-a3*b2-a2*b1-a1*b3+a1*b2+a3*b1;
162  Obj d = sqrt((a2*a2+b2*b2-r2*r2-2*a2*a3+a3*a3-r3*r3-2*b3*b2+b3*b3+2*r3*r2)*
163  (b3*b3+b1*b1-r1*r1-r3*r3-2*b1*b3+2*r3*r1+a3*a3-2*a1*a3+a1*a1)*
164  (a1*a1-2*a1*a2-r1*r1+b1*b1+b2*b2-r2*r2-2*b2*b1+2*r2*r1+a2*a2)*
165  tmp*tmp);
166  Obj v = d-b;
167 
168  if (v<0) return tlp::Circle<Obj, OTYPE>(0,0,0);
169 
170  Obj aa = -2*a3*b2*b2*a1-2*a2*b3*b3*a1
171  -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
172  -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
173  -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
174  -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
175  +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
176  -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
177  +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
178  -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
179  +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;
180  Obj R = 0.5*v/aa;
181  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
182  +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
183  -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);
184  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
185  -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
186  -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);
187  return tlp::Circle<Obj, OTYPE>(x,y,R);
188  }
189 
190  void process2() {
191  if (isEmpty()) {
192  result=tlp::enclosingCircle((*circles)[b1],(*circles)[b2]);
193  }
194  else {
195  unsigned selectedCircle=popBack();
196  process2();
197 
198  if (!(*circles)[selectedCircle].isIncludeIn(result)) {
199  result=enclosingCircle((*circles)[b1],(*circles)[b2],(*circles)[selectedCircle]);
200  pushFront(selectedCircle);
201  }
202  else {
203  pushBack(selectedCircle);
204  }
205  }
206  }
207 
208  void process1() {
209  if (isEmpty()) {
210  result=(*circles)[b1];
211  }
212  else {
213  unsigned selectedCircle=popBack();
214  process1();
215 
216  if (!(*circles)[selectedCircle].isIncludeIn(result)) {
217  b2=selectedCircle;
218  process2();
219  pushFront(selectedCircle);
220  }
221  else {
222  pushBack(selectedCircle);
223  }
224  }
225  }
226  void process0() {
227  if (isEmpty()) {
228  result=Circle<Obj, OTYPE>(0,0,0);
229  }
230  else {
231  unsigned selectedCircle=popBack();
232  process0();
233 
234  if (!(*circles)[selectedCircle].isIncludeIn(result)) {
235  b1=selectedCircle;
236  process1();
237  pushFront(selectedCircle);
238  }
239  else {
240  pushBack(selectedCircle);
241  }
242  }
243  }
244  bool isEmpty() const {
245  return first==(last+1)%enclosedCircles.size();
246  }
247  void pushFront(unsigned c) {
248  first=(first+enclosedCircles.size()-1)%enclosedCircles.size();
249  enclosedCircles[first]=c;
250  }
251  unsigned popBack() {
252  unsigned result=enclosedCircles[last];
253  last=(last+enclosedCircles.size()-1)%enclosedCircles.size();
254  return result;
255  }
256  void pushBack(unsigned c) {
257  last=(last+1)%enclosedCircles.size();
258  enclosedCircles[last]=c;
259  }
260  public:
261  tlp::Circle<Obj, OTYPE> operator()(const std::vector<tlp::Circle<Obj, OTYPE> > &circlesSet) {
262  circles=&circlesSet;
263  enclosedCircles.resize(circlesSet.size()+1);
264  first=0;
265  last=circlesSet.size()-1;
266 
267  for (unsigned i=0; i<circlesSet.size(); ++i)
268  enclosedCircles[i]=i;
269 
270  for (unsigned i=circlesSet.size(); i>0;) {
271  unsigned idx=(unsigned)(1.0*rand()*i/RAND_MAX);
272  --i;
273  std::swap(enclosedCircles[idx],enclosedCircles[i]);
274  }
275 
276  process0();
277  return result;
278  }
279  };
280  return OptimumCircleHull()(circles);
281 }
282 
283 template <typename Obj, typename OTYPE>
284 std::ostream& tlp::operator<<(std::ostream &os,const tlp::Circle<Obj, OTYPE> &a) {
285  os << "((" << a[0] << "," << a[1] << ")," << a.radius << ")";
286  return os;
287 }