1 // Copyright (c) 1997-2001 ETH Zurich (Switzerland). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org). 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Bounding_volumes/include/CGAL/Approximate_min_ellipsoid_d/Approximate_min_ellipsoid_d_debug.h $ 7 // $Id: Approximate_min_ellipsoid_d_debug.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Kaspar Fischer <fischerk@inf.ethz.ch> 12 13 #ifndef CGAL_APPROX_MIN_ELL_D_DEBUG_H 14 #define CGAL_APPROX_MIN_ELL_D_DEBUG_H 15 16 #include <CGAL/license/Bounding_volumes.h> 17 18 19 #include <cmath> 20 #include <map> 21 #include <iostream> 22 #include <fstream> 23 #include <sstream> 24 #include <boost/random/linear_congruential.hpp> 25 26 namespace CGAL { 27 28 // We hide the following debugging utilities in a "private" namespace: 29 namespace Approximate_min_ellipsoid_d_impl { 30 31 // For debugging only: 32 template<typename Iterator> print_matrix(int d,const char * name,Iterator A)33 void print_matrix(int d, const char *name, Iterator A) 34 { 35 std::cout << name << ":= Matrix([\n"; 36 for (int i=0; i<d; ++i) { 37 std::cout << " [ "; 38 for (int j=0; j<d; ++j) { 39 std::cout << std::setprecision(30) << A[i+j*d]; 40 if (j<d-1) 41 std::cout << ", "; 42 } 43 std::cout << " ]"; 44 if (i<d-1) 45 std::cout << ","; 46 std::cout << "\n"; 47 } 48 std::cout << "]);\n"; 49 } 50 51 // For debugging only: 52 template<typename Iterator> print_vector(int d,const char * name,Iterator v)53 void print_vector(int d, const char *name, Iterator v) 54 { 55 std::cout << name << ":= Vector([\n"; 56 for (int j=0; j<d; ++j) { 57 std::cout << std::setprecision(30) << v[j]; 58 if (j<d-1) 59 std::cout << ", "; 60 } 61 std::cout << "]);\n"; 62 } 63 64 #ifdef CGAL_APPEL_LOG_MODE 65 class Logger 66 // A singleton class which sends debugging comments to log files. 67 // (A class is called a "singleton" if at most one instance of it 68 // may ever exist.) By calling void log(channel,msg) you can send 69 // the string msg to the file with name channel. 70 { 71 private: // private members: 72 typedef std::map<std::string,std::ofstream*> Streams; 73 Streams channels; // a collection of pairs (k,v) where 74 // k is the file-name and v is the 75 // (open) stream associated with k 76 77 private: // (Construction and destruction are private to prevent 78 // more than one instantiation.) 79 Logger()80 Logger() {} 81 Logger(const Logger&); 82 Logger& operator=(const Logger&); 83 ~Logger()84 ~Logger() 85 { 86 // we walk through the list of all opened files and close them: 87 for (Streams::iterator it = channels.begin(); 88 it != channels.end(); ++it) { 89 (*it).second->close(); 90 delete (*it).second; 91 } 92 } 93 94 public: // access and routines: 95 instance()96 static Logger& instance() 97 // Returns a reference to the only existing instance of this class: 98 { 99 // Here's where we maintain the only instance: (Notice that it 100 // gets constructed automatically the first time instance() is 101 // called, and that it gets disposed of (if ever contructed) at 102 // program termination.) 103 static Logger instance; 104 return instance; 105 } 106 log(const char * channel,const std::string & msg)107 void log(const char *channel,const std::string& msg) 108 // If this is the first call to log with string channel as the 109 // first parameter, then the file with name channel.log is 110 // opened (at the beginning) for writing and msg is written to 111 // it. Otherwise, the string msg is opened to the already open 112 // file channel.log. 113 { 114 const std::string name(channel); 115 Streams::iterator it = channels.find(name); 116 117 // have we already opened this file? 118 if (it != channels.end()) { 119 // If so, then just append the message: 120 *(*it).second << msg; 121 (*it).second->flush(); 122 123 } else { 124 // If we haven't seen 'name' before, we create a new file: 125 using std::ofstream; 126 ofstream *o = new ofstream((name+".log").c_str(), 127 ofstream::out|ofstream::trunc); 128 channels[name] = o; 129 *o << msg; 130 } 131 } 132 }; 133 #endif // CGAL_APPEL_LOG_MODE 134 135 } // end of namespace Approximate_min_ellipsoid_d_impl 136 } // end of namespace CGAL 137 138 #ifdef CGAL_APPEL_TIMER_MODE 139 #include <sys/time.h> 140 #include <sys/resource.h> 141 142 namespace CGAL { 143 namespace Approximate_min_ellipsoid_d_impl { 144 145 // The following routine is taken from file mptimeval.h from 146 // "Matpack Library Release 1.7.1" which is copyright (C) 1991-2002 147 // by Berndt M. Gammel. It works on the timeval struct defined in 148 // sys/time.h: 149 // 150 // struct timeval { 151 // long tv_sec; /* seconds */ 152 // long tv_usec; /* microseconds */ 153 // }; 154 // 155 inline timeval& operator-=(timeval &t1,const timeval &t2) 156 { 157 t1.tv_sec -= t2.tv_sec; 158 if ( (t1.tv_usec -= t2.tv_usec) < 0 ) { 159 --t1.tv_sec; 160 t1.tv_usec += 1000000; 161 } 162 return t1; 163 } 164 165 // (Adapted from the above.) 166 inline timeval& operator+=(timeval &t1,const timeval &t2) 167 { 168 t1.tv_sec += t2.tv_sec; 169 if ( (t1.tv_usec += t2.tv_usec) > 1000000) { 170 ++t1.tv_sec; 171 t1.tv_usec -= 1000000; 172 } 173 return t1; 174 } 175 176 class Timer 177 // A singleton class which maintains a collection of named timers. 178 // (A class is called a "singleton" if at most one instance of it 179 // may ever exist.) The following routines are provided: 180 // 181 // - start(name): If this is the first time start() has been 182 // called with name as the first parameter, then a new timer is 183 // created and started. Otherwise, the timer with name name is 184 // restarted. 185 // 186 // - lapse(name): Retuns the number of seconds which have elapsed 187 // since start(name) was called last. 188 // Precondition: start(name) has been called once. 189 { 190 private: // (Construction and destruction are private to prevent 191 // more than one instantiation.) 192 Timer()193 Timer() {} 194 Timer(const Timer&); 195 Timer& operator=(const Timer&); 196 197 public: // access and routines: 198 instance()199 static Timer& instance() 200 // Returns a reference to the only existing instance of this class: 201 { 202 // Here's where we maintain the only instance: (Notice that it 203 // gets constructed automatically the first time instance() is 204 // called, and that it gets disposed of (if ever contructed) at 205 // program termination.) 206 static Timer instance; 207 return instance; 208 } 209 start(const char * timer_name)210 void start(const char *timer_name) 211 { 212 // fetch current usage: 213 rusage now; 214 int status = getrusage(RUSAGE_SELF,&now); 215 CGAL_assertion(status == 0); 216 217 // save it: 218 timers[std::string(timer_name)] = now.ru_utime; 219 } 220 lapse(const char * name)221 float lapse(const char *name) 222 { 223 // assert that start(name) has been called before: 224 CGAL_assertion(timers.find(std::string(name)) != timers.end()); 225 226 // get current usage: 227 rusage now; 228 int status = getrusage(RUSAGE_SELF,&now); 229 CGAL_assertion(status == 0); 230 231 // compute elapsed usage: 232 now.ru_utime -= (*timers.find(std::string(name))).second; 233 return now.ru_utime.tv_sec + now.ru_utime.tv_usec * 1e-6; 234 } 235 236 private: // private members: 237 typedef std::map<std::string,timeval> Timers; 238 Timers timers; // a collection of pairs (k,v) where 239 // k is the timer name and v is the 240 // (started) timer associated with k 241 }; 242 } // end of namespace Approximate_min_ellipsoid_d_impl 243 } // end of namespace CGAL 244 #endif // CGAL_APPEL_TIMER_MODE 245 246 namespace CGAL { 247 namespace Approximate_min_ellipsoid_d_impl { 248 249 template <typename T> tostr(const T & t)250 std::string tostr(const T& t) { 251 std::stringstream strm; 252 strm << t; 253 return strm.str(); 254 } 255 256 class Eps_export_2 { 257 // An instance of the following class accepts circles and ellipses 258 // and procudes an Enhanced-PostScript figure. 259 public: 260 enum Stroke_mode { Solid=0, Solid_filled=1, Dashed=2 }; 261 enum Label_mode { None, Angle, Random_angle }; 262 263 private: 264 std::vector<double> bb; // bounding box 265 Stroke_mode bm; // current stroke mode 266 Label_mode lm; // current label mode 267 double next_angle; // next angle to use (for Angle mode) 268 int count; // counts the objects (starting at 1) 269 std::string next_label; // next label to use 270 std::ostringstream body; // buffered output 271 bool adjust_bb; 272 double zoom; 273 std::string filename; 274 boost::rand48 rng; 275 276 public: // construction and destruction: 277 278 // Begins exporting to file filename. Sets the current stroke mode 279 // to Solid and the current label-mode to Random_angle. 280 Eps_export_2(std::string filename, double zoom, int /* seed */ = 0) 281 : bb(4,0.0), bm(Solid), lm(Random_angle), count(1), 282 next_label(default_label(1)), adjust_bb(true), 283 zoom(zoom), filename(filename) 284 {} 285 286 // Ends exporting and closes the file. ~Eps_export_2()287 ~Eps_export_2() 288 { 289 // open output file: 290 using std::endl; 291 std::ofstream f(filename.c_str()); 292 if (!f) 293 std::cerr << "Couldn't open file " << filename << "." << endl; 294 295 // write header: 296 const double margin = 10.0; 297 f << "%!PS-Adobe-2.0 EPSF-2.0" << endl 298 << "%%Title: Some balls" << endl 299 << "%%Creator: a program by kf." << endl 300 << "%%CreationDate: now" << endl 301 << "%%For: you, the unknown" << endl 302 << "%%Pages: 1" << endl 303 << "%%DocumentFonts: Helvetica" << endl 304 << "%%BoundingBox: " 305 << bb[0]*zoom-margin << " " 306 << bb[1]*zoom-margin << " " 307 << bb[2]*zoom+margin << " " 308 << bb[3]*zoom+margin << endl 309 << "%%Magnification: 1.0000" << endl 310 << "%%EndComments" << endl 311 << "%%BeginProlog" << endl 312 << "%" << endl 313 << "/zoom " << zoom << " def" << endl 314 << "zoom zoom scale" << endl 315 << "/Helvetica findfont 7 zoom div scalefont setfont" << endl 316 << "0.06299 7.500 mul zoom div setlinewidth" << endl 317 << "/dashlen 4.5 zoom div def" << endl 318 << "/dashlen2 1.5 zoom div def" << endl 319 << "/ptlen 3.5 zoom div def" << endl << "%" << endl 320 << "% center-x center-y radius bp" << endl 321 << "% creates a path representing a circle" << endl 322 << "/bp { /radius exch def /cy exch def /cx exch def" << endl 323 << "radius 0 eq" << endl 324 << "{ newpath cx ptlen sub cy moveto ptlen ptlen add 0" << endl 325 << "rlineto cx cy ptlen sub moveto 0 ptlen ptlen add" << endl 326 << "rlineto }" << endl 327 << "{ newpath cx radius add cy moveto cx cy radius" << endl 328 << "0 360 arc closepath } ifelse } def" << "%" << endl 329 << "%center-x center-y radius ball" << endl 330 << "% draws a circle, unfilled, no dash" << endl 331 << "/ball { bp stroke } def" << endl 332 << "% center-x center-y radius bball" << endl 333 << "% draws a circle, filled, no dash" << endl 334 << "/bball { /radius exch def /cy exch def /cx exch def" << endl 335 << "cx cy radius radius 0 eq" << endl 336 << "{ gsave bp 0.6 setgray stroke" << endl 337 << "newpath cx cy moveto closepath 1 setlinecap" << endl 338 << "0.0 setgray stroke grestore }" << endl 339 << "{ bp gsave 0.6 setgray fill grestore stroke }" << endl 340 << "ifelse } def" << endl 341 << "% center-x center-y radius dball" << endl 342 << "% draws a circle, unfilled, dashed" << endl 343 << "/dball { gsave bp [dashlen] 0 setdash stroke grestore" << endl 344 << "} def" << endl 345 << "/eball { gsave bp 0.6 setgray stroke grestore } def" << endl 346 << "%" << endl 347 << "% center-x center-y radius label dist angle drawl" << endl 348 << "% Draws the label label." << endl 349 << "/drawl { /angle exch def /dist exch def /label exch def" << endl 350 << "/radius exch def /cy exch def /cx exch def" << endl 351 << "cx radius dist zoom div add angle cos mul add" << endl 352 << "cy radius dist zoom div add angle sin mul add" << endl 353 << "moveto label show } def" << endl 354 << "%" << endl 355 << "% Usage a b c find-roots x1 x2" << endl 356 << "% Finds the two roots of a x^2 + b x + c = 0, assuming" << endl 357 << "% that a is nonzero." << endl 358 << "/find-roots {" << endl 359 << " % compute discriminant:" << endl 360 << " 3 copy 3 2 roll" << endl 361 << " mul 4.0 mul neg exch dup mul add" << endl 362 << " % stack: a b c disc" << endl 363 << " sqrt" << endl 364 << " 2 index 0 ge { neg } { } ifelse" << endl 365 << " % stack: a b c sqrt(disc)" << endl 366 << " % compute first solution (sqrt(disc)-b) / (2*a):" << endl 367 << " dup 3 index sub 4 index 2.0 mul div" << endl 368 << " 5 1 roll" << endl 369 << " % stack: x1 a b c sqrt(disc)" << endl 370 << " % compute second solution 2*c / (sqrt(disc)-b):" << endl 371 << " 3 2 roll sub exch 2.0 mul exch div" << endl 372 << " exch pop" << endl 373 << "} def" << endl 374 << "%" << endl 375 << "% Usage: u v normalize u' v'" << endl 376 << "% Takes the vector [u,v] and normalizes it to length 1." << endl 377 << "/normalize {" << endl 378 << " dup dup mul" << endl 379 << " % stack: u v v^2" << endl 380 << " 2 index dup mul add" << endl 381 << " % stack: u v v^2+u^2" << endl 382 << " sqrt 1.0 exch div dup 4 1 roll mul 3 1 roll mul exch" << endl 383 << "} def" << endl 384 << "%" << endl 385 << "% Usage: a b h cellipse" << endl 386 << "% Draws the ellipse x^T M x <= 1 with M = [[a,h],[h,b]]." << endl 387 << "%" << endl 388 << "/cellipse {" << endl 389 << "% compute Eigen decomposition of M:" << endl 390 << "%" << endl 391 << "% stack: a b h" << endl 392 << "dup 0.0 eq {" << endl 393 << " % ellipse is a sphere:" << endl 394 << " 2 index sqrt" << endl 395 << " 2 index sqrt" << endl 396 << " [ 1.0 0.0 0.0 1.0 0.0 0.0 ]" << endl 397 << "}" << endl 398 << "{" << endl 399 << " 1.0" << endl 400 << " 4 copy pop pop add neg" << endl 401 << " 5 copy pop pop dup mul neg 3 1 roll mul add" << endl 402 << " find-roots" << endl 403 << " % stack: a b h x1 x2" << endl 404 << " %" << endl 405 << " % build first Eigenvector:" << endl 406 << " 2 index 2 index 6 index sub normalize" << endl 407 << " % build second Eigenvector:" << endl 408 << " % stack: a b h x1 x2 ev1x ev1y" << endl 409 << " 4 index 3 index 8 index sub normalize" << endl 410 << " % build matrix containing Eigenvectors:" << endl 411 << " % stack: a b h x1 x2 ev1x ev1y ev2x ev2y" << endl 412 << " 6 array" << endl 413 << " dup 0 6 index put" << endl 414 << " dup 1 4 index put" << endl 415 << " dup 2 5 index put" << endl 416 << " dup 3 3 index put" << endl 417 << " dup 4 0 put" << endl 418 << " dup 5 0 put" << endl 419 << " 5 1 roll pop pop pop pop" << endl 420 << "} ifelse" << endl 421 << " % stack: a b h x1 x2 T" << endl 422 << " %" << endl 423 << " % We now draw an circle scaled by x1 (along the " << endl 424 << " % x-axis) and x2 (along the y-axis):" << endl 425 << " matrix currentmatrix % remember CTM" << endl 426 << " % stack: a b h x1 x2 T old-ctm" << endl 427 << " exch matrix invertmatrix concat" << endl 428 << " 2 index sqrt 1.0 exch div 2 index sqrt 1.0" << endl 429 << " exch div scale" << endl 430 << " newpath" << endl 431 << " 0 0 1 0 360 arc closepath" << endl 432 << " setmatrix % restore CTM" << endl 433 << " pop pop pop pop pop" << endl 434 << "} def" << endl 435 << "%" << endl 436 << "% Usage: a b d u v mu ellipse" << endl 437 << "% Finds the path (i.e., border) of the ellipse" << endl 438 << "%" << endl 439 << "% E = { x | x^T M x + x^T m + mu <= 0 }," << endl 440 << "%" << endl 441 << "% where" << endl 442 << "%" << endl 443 << "% [ a b ] [ u ]" << endl 444 << "% M = [ ], m = [ ]" << endl 445 << "% [ b d ] [ v ]" << endl 446 << "%" << endl 447 << "% and det(M) > 0." << endl 448 << "/ellipse {" << endl 449 << " /emu exch def" << endl 450 << " /ev exch def" << endl 451 << " /eu exch def" << endl 452 << " /ed exch def" << endl 453 << " /eb exch def" << endl 454 << " /ea exch def" << endl 455 << " /ematrix [ ea eb eb ed 0 0 ] matrix invertmatrix def" << endl 456 << " % compute z = M^{-1} m:" << endl 457 << " eu ev ematrix transform" << endl 458 << " /ezy exch def" << endl 459 << " /ezx exch def" << endl 460 << " % translate to center:" << endl 461 << " matrix currentmatrix % remember CTM" << endl 462 << " ezx neg 2 div ezy neg 2 div translate" << endl 463 << " % compute matrix of (now centrally symmetric) ellipse:" << endl 464 << " /efactor ezx eu mul ezy ev mul add 4.0 div emu sub def" << endl 465 << " ea efactor div ed efactor div eb efactor div cellipse" << endl 466 << " setmatrix % restore CTM" << endl 467 << "} def" << endl 468 << "/setcol { 193 255 div 152 255 div 159 255 div" << endl 469 << " setrgbcolor" << endl 470 << "} def" << endl 471 << "%" << endl 472 << "%" << endl 473 << "%%EndProlog" << endl 474 << "%%Page: 1 1" << endl; 475 476 // write content: 477 f << "gsave" << endl 478 << body.str() 479 << "grestore" << endl 480 << "showpage" << endl 481 << "%%Trailer" << endl; 482 483 // write footer: 484 f.close(); 485 } 486 487 public: // modifiers: 488 // Sets the stroke mode for the next objects. set_stroke_mode(Stroke_mode m)489 void set_stroke_mode(Stroke_mode m) 490 { 491 bm = m; 492 } 493 494 // Sets the labelmode for the next objects. set_label_mode(Label_mode m)495 void set_label_mode(Label_mode m) 496 { 497 lm = m; 498 } 499 500 // Sets the label for the next object. (This only makes sense if 501 // the label mode for the next object will be different from None.) set_label(const std::string & l)502 void set_label(const std::string& l) 503 { 504 next_label = l; 505 } 506 507 // Sets the angle for the next object drawn in mode Angle or 508 // Random_angle. (So you can temporarily overwrite the random 509 // angle.) set_angle(double angle)510 void set_angle(double angle) 511 { 512 next_angle = angle; 513 } 514 515 // Enable or disable automatic adjusting of the bounding box. enlarge_bounding_box(bool active)516 void enlarge_bounding_box(bool active) 517 { 518 adjust_bb = active; 519 } 520 521 // Renders the line from (x,y) to (u,v) in the current stroke 522 // and label mode. 523 // 524 // Implementation note: stroke and label mode are not yet implemented. 525 void write_line(double x, double y, double u, double v, 526 bool colored = false) 527 { 528 // set color: 529 if (colored) 530 body << "gsave setcol" << std::endl; 531 532 // output ball in appropriate mode: 533 adjust_bounding_box(x,y); 534 adjust_bounding_box(u,v); 535 body << "newpath\n" 536 << x << ' ' << y << " moveto\n" 537 << u << ' ' << v << " lineto\n" 538 << "stroke\n"; 539 540 // restore color: 541 if (colored) 542 body << "grestore" << std::endl; 543 } 544 545 // Renders the circle with center (x,y) and radius r in 546 // the current stroke and label mode. All coordinates must be of 547 // type double. 548 void write_circle(double x, double y, double r, bool colored = false) 549 { 550 // set color: 551 if (colored) 552 body << "gsave setcol" << std::endl; 553 554 // output ball in appropriate mode: 555 adjust_bounding_box(x-r,y-r); 556 adjust_bounding_box(x+r,y+r); 557 body << x << ' ' << y << ' ' << r << ' '; 558 const char *mode[] = { "ball", "bball", "dball" }; 559 body << mode[static_cast<int>(bm)] << std::endl; 560 561 // output label: 562 if (lm != None) { 563 const double dist = 5.0; 564 const double lx = x + (r+dist/zoom)*std::cos(next_angle); 565 const double ly = y + (r+dist/zoom)*std::sin(next_angle); 566 adjust_bounding_box(lx,ly); 567 body << x << ' ' << y << ' ' << r << " (" << next_label << ") " 568 << dist << ' ' << next_angle << " drawl" << std::endl; 569 } 570 571 // restore color: 572 if (colored) 573 body << "grestore" << std::endl; 574 575 // generate next angle (if needed): 576 if (lm == Random_angle) 577 next_angle = static_cast<double>(rng()); 578 579 // update counter and label: 580 next_label = default_label(++count); 581 } 582 583 // Renders the centrally symmetric ellipse x^T M x <= 1 with M = 584 // [[a,h],[h,b]] in the current stroke and label mode. write_cs_ellipse(double a,double b,double h)585 void write_cs_ellipse(double a,double b,double h) 586 { 587 using std::endl; 588 589 // Adjust the bounding box: For this, we use the fact that the 590 // ellipse has semi-axes of lengths (1/sqrt(c),1/sqrt(d)) where 591 // 592 // (c,d) = find_roots(1.0,-(a+b),a b - h^2) 593 // 594 // (See Sven's thesis, code on page 87.) So it it enclosed in 595 // some (rotated) rectangle of side lengths 2*c,2*d, centered 596 // at the origin. Consequently, the circle of radius 597 // r=sqrt(c^2+d^2) encloses the ellipse, and we use this to 598 // find the bounding box: 599 const std::pair<double,double> semi = find_roots(1.0,-(a+b),a*b-h*h); 600 const double r = std::sqrt(1/semi.first+1/semi.second); 601 adjust_bounding_box(-r,-r); 602 adjust_bounding_box(r,r); 603 604 // draw the ellipse: 605 body << "gsave" << endl 606 << " " << a << ' ' << b << ' ' << h << " cellipse" << endl; 607 if (bm == Solid_filled) 608 body << " gsave 0.6 setgray eofill grestore" << endl; 609 if (bm == Dashed) 610 body << " [dashlen] 0 setdash" << endl; 611 body << " stroke" << endl << "grestore" << endl; 612 613 // output label: 614 // Todo: Not implemented yet. 615 if (false && lm != None) { 616 const double distance = 5.0; 617 body << a << ' ' << b << ' ' << distance << ' ' << next_angle 618 << " cellipse-pt moveto (" << next_label << ") show" << endl; 619 } 620 621 // generate next angle (if needed): 622 if (lm == Random_angle) 623 next_angle = static_cast<double>(rng()); 624 625 // update counter and label: 626 next_label = default_label(++count); 627 } 628 629 // Renders the ellipse 630 // 631 // E = { x | x^T M x + x^T m + mu <= 0 } 632 // 633 // where 634 // 635 // [ a h ] [ u ] 636 // M = [ ], m = [ ] 637 // [ h b ] [ v ] 638 // 639 // and det(M) > 0 in the current stroke and label mode. write_ellipse(double a,double b,double h,double u,double v,double mu)640 void write_ellipse(double a,double b,double h, 641 double u,double v,double mu) 642 { 643 using std::endl; 644 645 // adjust bounding box (see file boundingbox.mw and my file note): 646 const double tmp1 = a*b-h*h; 647 const double tmp2 = ((u*(u*b-2.0*v*h)+v*v*a)/tmp1-4.0*mu)/tmp1; 648 const double hw = 0.5*std::sqrt(b*tmp2), vw = 0.5*std::sqrt(a*tmp2); 649 const double hoff = -0.5*(u*b-v*h)/tmp1, voff = 0.5*(u*h-v*a)/tmp1; 650 adjust_bounding_box((std::min)(hw+hoff,-hw+hoff), 651 (std::min)(vw+voff,-vw+voff)); 652 adjust_bounding_box((std::max)(hw+hoff,-hw+hoff), 653 (std::max)(vw+voff,-vw+voff)); 654 655 #if 0 656 // draw bounding box: 657 body << "newpath " << (std::min)(hw+hoff,-hw+hoff) << " " 658 << (std::min)(vw+voff,-vw+voff) << " moveto " 659 << (std::max)(hw+hoff,-hw+hoff) << " " 660 << (std::min)(vw+voff,-vw+voff) << " lineto " 661 << (std::max)(hw+hoff,-hw+hoff) << " " 662 << (std::max)(vw+voff,-vw+voff) << " lineto " 663 << (std::min)(hw+hoff,-hw+hoff) << " " 664 << (std::max)(vw+voff,-vw+voff) << " lineto closepath stroke\n"; 665 #endif 666 667 // begin drawing: 668 body << "gsave" << endl; 669 670 // draw the ellipse: 671 body << " " << a << ' ' << h << ' ' << b 672 << ' ' << u << ' ' << v << ' ' << mu 673 << " ellipse" << endl; 674 if (bm == Solid_filled) 675 body << " gsave 0.6 setgray eofill grestore" << endl; 676 if (bm == Dashed) 677 body << " [dashlen] 0 setdash" << endl; 678 body << " stroke" << endl; 679 680 // output label: 681 // Todo: Not implemented yet. 682 683 // end drawing: 684 body << "grestore" << endl; 685 686 // generate next angle (if needed): 687 if (lm == Random_angle) 688 next_angle = static_cast<double>(rng()); 689 690 // update counter and label: 691 next_label = default_label(++count); 692 } 693 adjust_bounding_box(double x,double y)694 void adjust_bounding_box(double x,double y) 695 // Make sure the bounding box is large enough to contain the point (x,y). 696 { 697 if (adjust_bb) { 698 bb[0] = (std::min)(x,bb[0]); 699 bb[2] = (std::max)(x,bb[2]); 700 bb[1] = (std::min)(y,bb[1]); 701 bb[3] = (std::max)(y,bb[3]); 702 } 703 } 704 705 private: // utilities: 706 find_roots(double a,double b,double c)707 std::pair<double,double> find_roots(double a,double b,double c) 708 // Assuming a is nonzero, finds the roots of a x^2 + b x + c = 0. 709 { 710 double sd = std::sqrt(b*b-4.0*a*c); 711 if (b >= 0.0) 712 sd = -sd; 713 return std::pair<double,double>( (sd-b)/(2*a), 2*c/(sd-b) ); 714 } 715 default_label(int count)716 static std::string default_label(int count) 717 { 718 return tostr(count); 719 } 720 }; 721 722 } // namespace Approximate_min_ellipsoid_d_impl 723 724 } // namespace CGAL 725 726 #endif // CGAL_APPROX_MIN_ELL_D_DEBUG_H 727