Tulip  5.1.0
Large graphs analysis and drawing
Circle.cxx
1 /*
2  *
3  * This file is part of Tulip (http://tulip.labri.fr)
4  *
5  * Authors: David Auber and the Tulip development Team
6  * from LaBRI, University of Bordeaux
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 #include <tulip/TlpTools.h>
21 
22 template<typename Obj, typename OTYPE>
23 tlp::Circle<Obj, OTYPE> & tlp::Circle<Obj, OTYPE>::merge(const tlp::Circle<Obj, OTYPE> &c) {
24  Vector<Obj,2, OTYPE> p1(*this);
25  Vector<Obj,2, OTYPE> p2(c);
26  Vector<Obj,2, OTYPE> c12=p2-p1;
27  double norm=c12.norm();
28 
29  if (norm < 0.0000001) {
30  if (radius<c.radius)
31  return (*this)=c;
32  else
33  return *this;
34  }
35 
36  double dx=c12[0]/norm;
37  double dy=c12[1]/norm;
38  p1[0]-=dx*radius;
39  p1[1]-=dy*radius;
40  p2[0]+=dx*c.radius;
41  p2[1]+=dy*c.radius;
42  Obj tmpDist=p1.dist(p2)/2.;
43 
44  if ( (tmpDist<radius) || (tmpDist<c.radius)) {
45  if (radius>c.radius)
46  return *this;
47  else
48  return (*this)=c;
49  }
50  else {
51  p1+=p2;
52  p1/=2.;
53  (*this)[0]=p1[0];
54  (*this)[1]=p1[1];
55  radius=tmpDist;
56  return *this;
57  }
58 }
59 
60 template<typename Obj, typename OTYPE>
61 bool tlp::Circle<Obj, OTYPE>::isIncludeIn(const tlp::Circle<Obj, OTYPE> &c) const {
62  Vector<Obj,2, OTYPE> dir=c-*this;
63  return (dir.norm()+radius)<=c.radius;
64 }
65 
66 template<typename Obj, typename OTYPE>
67 tlp::Circle<Obj, OTYPE> tlp::enclosingCircle(const tlp::Circle<Obj, OTYPE> &c1,const tlp::Circle<Obj, OTYPE> &c2) {
68  Vector<Obj,2, OTYPE> dir=c2-c1;
69  Obj n=dir.norm();
70 
71  if (n==0)
72  return Circle<Obj, OTYPE>(c1, std::max(c1.radius, c2.radius));
73 
74  dir/=n;
75  Vector<Obj, 2, OTYPE> ext1=c1-dir*c1.radius;
76  Vector<Obj, 2, OTYPE> ext2=c2+dir*c2.radius;
77  return Circle<Obj, OTYPE>((ext1+ext2)/Obj(2),(ext2-ext1).norm()/Obj(2));
78 }
79 
80 template<typename Obj, typename OTYPE>
81 tlp::Circle<Obj, OTYPE> tlp::lazyEnclosingCircle(const std::vector<tlp::Circle<Obj, OTYPE> > & circles) {
82  //compute bounding box of a
83  tlp::Vector<Obj,4, OTYPE> boundingBox;
84  // for (int i=0;i<4;++i) boundingBox[i]=0;
85  typename std::vector< tlp::Circle<Obj, OTYPE> >::const_iterator it=circles.begin();
86  boundingBox[0]=(*it)[0]-(*it).radius;
87  boundingBox[1]=(*it)[1]-(*it).radius;
88  boundingBox[2]=(*it)[0]+(*it).radius;
89  boundingBox[3]=(*it)[1]+(*it).radius;
90  ++it;
91 
92  for (; it!=circles.end(); ++it) {
93  boundingBox[0] = std::min(boundingBox[0] , ((*it)[0]-(*it).radius));
94  boundingBox[1] = std::min(boundingBox[1] , ((*it)[1]-(*it).radius));
95  boundingBox[2] = std::max(boundingBox[2] , ((*it)[0]+(*it).radius));
96  boundingBox[3] = std::max(boundingBox[3] , ((*it)[1]+(*it).radius));
97  }
98 
99  tlp::Vector<Obj,2, OTYPE> center;
100  center[0]=(boundingBox[0]+boundingBox[2])/2.;
101  center[1]=(boundingBox[1]+boundingBox[3])/2.;
102  Obj radius = std::max( (boundingBox[2]-boundingBox[0])/2., (boundingBox[3]-boundingBox[1])/2. );
103  tlp::Circle<Obj, OTYPE> result(center,radius);
104 
105  //compute circle hull
106  for (typename std::vector< tlp::Circle<Obj, OTYPE> >::const_iterator it=circles.begin(); it!=circles.end(); ++it)
107  result.merge(*it);
108 
109  return result;
110 }
111 
112 template<typename Obj, typename OTYPE>
113 tlp::Circle<Obj, OTYPE> tlp::enclosingCircle(const std::vector<tlp::Circle<Obj, OTYPE> > & circles) {
114  class OptimumCircleHull {
115  const std::vector<tlp::Circle<Obj, OTYPE> > *circles;
116  std::vector<unsigned> enclosedCircles;
117  unsigned first,last;
118  unsigned b1,b2;
119  tlp::Circle<Obj, OTYPE> result;
120 
121  static tlp::Circle<Obj, OTYPE> enclosingCircle(const tlp::Circle<Obj, OTYPE> &c1,const tlp::Circle<Obj, OTYPE> &c2,const tlp::Circle<Obj, OTYPE> &c3) {
122  Obj a1 = c1[0];
123  Obj b1 = c1[1];
124  Obj r1 = c1.radius;
125  Obj a2 = c2[0];
126  Obj b2 = c2[1];
127  Obj r2 = c2.radius;
128  Obj a3 = c3[0];
129  Obj b3 = c3[1];
130  Obj r3 = c3.radius;
131  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
132  +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
133  -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
134  +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
135  -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
136  -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
137  -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
138  -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
139  -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
140  +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
141  -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
142  -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
143  +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
144  -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
145  +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
146  -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
147  +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
148  +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
149  -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
150  -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
151  -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
152  -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
153  +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
154  -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
155  -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
156  -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
157  +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
158  -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
159  -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
160  +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
161  +4*r1*a3*b2*a2*b3;
162  Obj tmp=a2*b3-a3*b2-a2*b1-a1*b3+a1*b2+a3*b1;
163  Obj d = sqrt((a2*a2+b2*b2-r2*r2-2*a2*a3+a3*a3-r3*r3-2*b3*b2+b3*b3+2*r3*r2)*
164  (b3*b3+b1*b1-r1*r1-r3*r3-2*b1*b3+2*r3*r1+a3*a3-2*a1*a3+a1*a1)*
165  (a1*a1-2*a1*a2-r1*r1+b1*b1+b2*b2-r2*r2-2*b2*b1+2*r2*r1+a2*a2)*
166  tmp*tmp);
167  Obj v = d-b;
168 
169  if (v<0) return tlp::Circle<Obj, OTYPE>(0,0,0);
170 
171  Obj aa = -2*a3*b2*b2*a1-2*a2*b3*b3*a1
172  -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
173  -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
174  -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
175  -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
176  +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
177  -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
178  +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
179  -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
180  +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;
181  Obj R = 0.5*v/aa;
182  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
183  +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
184  -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);
185  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
186  -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
187  -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);
188  return tlp::Circle<Obj, OTYPE>(x,y,R);
189  }
190 
191  void process2() {
192  if (isEmpty()) {
193  result=tlp::enclosingCircle((*circles)[b1],(*circles)[b2]);
194  }
195  else {
196  unsigned selectedCircle=popBack();
197  process2();
198 
199  if (!(*circles)[selectedCircle].isIncludeIn(result)) {
200  result=enclosingCircle((*circles)[b1],(*circles)[b2],(*circles)[selectedCircle]);
201  pushFront(selectedCircle);
202  }
203  else {
204  pushBack(selectedCircle);
205  }
206  }
207  }
208 
209  void process1() {
210  if (isEmpty()) {
211  result=(*circles)[b1];
212  }
213  else {
214  unsigned selectedCircle=popBack();
215  process1();
216 
217  if (!(*circles)[selectedCircle].isIncludeIn(result)) {
218  b2=selectedCircle;
219  process2();
220  pushFront(selectedCircle);
221  }
222  else {
223  pushBack(selectedCircle);
224  }
225  }
226  }
227  void process0() {
228  if (isEmpty()) {
229  result=Circle<Obj, OTYPE>(0,0,0);
230  }
231  else {
232  unsigned selectedCircle=popBack();
233  process0();
234 
235  if (!(*circles)[selectedCircle].isIncludeIn(result)) {
236  b1=selectedCircle;
237  process1();
238  pushFront(selectedCircle);
239  }
240  else {
241  pushBack(selectedCircle);
242  }
243  }
244  }
245  bool isEmpty() const {
246  return first==(last+1)%enclosedCircles.size();
247  }
248  void pushFront(unsigned c) {
249  first=(first+enclosedCircles.size()-1)%enclosedCircles.size();
250  enclosedCircles[first]=c;
251  }
252  unsigned popBack() {
253  unsigned result=enclosedCircles[last];
254  last=(last+enclosedCircles.size()-1)%enclosedCircles.size();
255  return result;
256  }
257  void pushBack(unsigned c) {
258  last=(last+1)%enclosedCircles.size();
259  enclosedCircles[last]=c;
260  }
261  public:
262  OptimumCircleHull():circles(NULL),first(0),last(0),b1(0),b2(0) {}
263 
264  tlp::Circle<Obj, OTYPE> operator()(const std::vector<tlp::Circle<Obj, OTYPE> > &circlesSet) {
265  circles=&circlesSet;
266  enclosedCircles.resize(circlesSet.size()+1);
267  first=0;
268  last=circlesSet.size()-1;
269 
270  for (unsigned i=0; i<circlesSet.size(); ++i)
271  enclosedCircles[i]=i;
272 
273  for (unsigned i=circlesSet.size(); i>0;) {
274  unsigned idx = tlp::randomUnsignedInteger(i-1);
275  --i;
276  std::swap(enclosedCircles[idx],enclosedCircles[i]);
277  }
278 
279  process0();
280  return result;
281  }
282  };
283  return OptimumCircleHull()(circles);
284 }
285 
286 template <typename Obj, typename OTYPE>
287 std::ostream& tlp::operator<<(std::ostream &os,const tlp::Circle<Obj, OTYPE> &a) {
288  os << "((" << a[0] << "," << a[1] << ")," << a.radius << ")";
289  return os;
290 }
unsigned int randomUnsignedInteger(unsigned int max)
Returns a random unsigned integer in the range [0, max].