1 #include <U2Core/disable-warnings.h> 2 U2_DISABLE_WARNINGS 3 4 // -*- Mode: C++; tab-width: 2; -*- 5 // vi: set ts=2: 6 // 7 8 #include <BALL/STRUCTURE/solventExcludedSurface.h> 9 #include <BALL/STRUCTURE/SESEdge.h> 10 #include <BALL/STRUCTURE/SESFace.h> 11 #include <BALL/STRUCTURE/SESVertex.h> 12 #include <BALL/STRUCTURE/reducedSurface.h> 13 #include <BALL/MATHS/analyticalGeometry.h> 14 #include <BALL/MATHS/angle.h> 15 #include <BALL/MATHS/circle3.h> 16 #include <BALL/MATHS/matrix44.h> 17 #include <BALL/MATHS/quaternion.h> 18 #include <BALL/MATHS/sphere3.h> 19 #include <BALL/MATHS/vector3.h> 20 #include <BALL/MATHS/vector4.h> 21 #include <BALL/DATATYPE/hashGrid.h> 22 #include <BALL/DATATYPE/hashMap.h> 23 #include <vector> 24 #include <list> 25 #include <iostream> 26 27 28 namespace BALL 29 { 30 SolventExcludedSurface()31 SolventExcludedSurface::SolventExcludedSurface() 32 : number_of_vertices_(0), 33 vertices_(), 34 number_of_edges_(0), 35 edges_(), 36 number_of_singular_edges_(0), 37 singular_edges_(), 38 number_of_contact_faces_(0), 39 contact_faces_(), 40 number_of_toric_faces_(0), 41 toric_faces_(), 42 number_of_spheric_faces_(0), 43 spheric_faces_(), 44 reduced_surface_(NULL) 45 { 46 } 47 48 SolventExcludedSurface(const SolventExcludedSurface & ses,bool)49 SolventExcludedSurface::SolventExcludedSurface 50 (const SolventExcludedSurface& ses, bool) 51 : number_of_vertices_(ses.vertices_.size()), 52 vertices_(number_of_vertices_), 53 number_of_edges_(ses.edges_.size()), 54 edges_(number_of_edges_), 55 number_of_singular_edges_(0), 56 singular_edges_(), 57 number_of_contact_faces_(ses.contact_faces_.size()), 58 contact_faces_(number_of_contact_faces_), 59 number_of_toric_faces_(ses.toric_faces_.size()), 60 toric_faces_(number_of_toric_faces_), 61 number_of_spheric_faces_(ses.spheric_faces_.size()), 62 spheric_faces_(number_of_spheric_faces_), 63 reduced_surface_(NULL) 64 { 65 //for (Position i = 0; i < number_of_vertices_; i++) 66 //{ 67 // vertices_[i] = new SESVertex(*ses.vertices_[i],false); 68 //} 69 //for (Position i = 0; i < number_of_edges_; i++) 70 //{ 71 // edges_[i] = new SESVertex(*ses.edges_[i],false); 72 //} 73 //::std::list<SESEdge*>::const_iterator se; 74 //for (se = ses.singular_edges_.begin(); se != ses.singular_edges_.end(); se++) 75 //{ 76 // singular_edges_.push_back(edges_[(*se)->index_]); 77 //} 78 //for (Position i = 0; i < number_of_contact_faces_; i++) 79 //{ 80 // contact_faces_[i] = new SESFace(*ses.contact_faces_[i],false); 81 //} 82 //for (Position i = 0; i < number_of_toric_faces_; i++) 83 //{ 84 // toric_faces_[i] = new SESFace(*ses.toric_faces_[i],false); 85 //} 86 //for (Position i = 0; i < number_of_spheric_faces_; i++) 87 //{ 88 // spheric_faces_[i] = new SESFace(*ses.spheric_faces_[i],false); 89 //} 90 } 91 92 SolventExcludedSurface(ReducedSurface * reduced_surface)93 SolventExcludedSurface::SolventExcludedSurface 94 (ReducedSurface* reduced_surface) 95 : number_of_vertices_(0), 96 vertices_(), 97 number_of_edges_(0), 98 edges_(), 99 number_of_singular_edges_(0), 100 singular_edges_(), 101 number_of_contact_faces_(0), 102 contact_faces_(), 103 number_of_toric_faces_(0), 104 toric_faces_(), 105 number_of_spheric_faces_(0), 106 spheric_faces_(), 107 reduced_surface_(reduced_surface) 108 { 109 } 110 111 ~SolventExcludedSurface()112 SolventExcludedSurface::~SolventExcludedSurface() 113 { 114 clear(); 115 } 116 117 clear()118 void SolventExcludedSurface::clear() 119 { 120 Position i; 121 for (i = 0; i < number_of_vertices_; i++) 122 { 123 delete vertices_[i]; 124 } 125 for (i = 0; i < number_of_edges_; i++) 126 { 127 delete edges_[i]; 128 } 129 for (i = 0; i < number_of_contact_faces_; i++) 130 { 131 delete contact_faces_[i]; 132 } 133 for (i = 0; i < number_of_toric_faces_; i++) 134 { 135 delete toric_faces_[i]; 136 } 137 for (i = 0; i < number_of_spheric_faces_; i++) 138 { 139 delete spheric_faces_[i]; 140 } 141 vertices_.clear(); 142 edges_.clear(); 143 singular_edges_.clear(); 144 contact_faces_.clear(); 145 toric_faces_.clear(); 146 spheric_faces_.clear(); 147 number_of_vertices_ = 0; 148 number_of_edges_ = 0; 149 number_of_singular_edges_ = 0; 150 number_of_contact_faces_ = 0; 151 number_of_toric_faces_ = 0; 152 number_of_spheric_faces_ = 0; 153 } 154 155 clean(const double & density)156 void SolventExcludedSurface::clean(const double& density) 157 { 158 SESFace* face(0); 159 bool done = false; 160 double sqrt_density = sqrt(density); 161 while (!done) 162 { 163 done = true; 164 for (Position i = 0; i < toric_faces_.size(); i++) 165 { 166 face = toric_faces_[i]; 167 if (face != NULL) 168 { 169 if (!face->isFree()) 170 { 171 if (face->type_ == SESFace::TYPE_TORIC_SINGULAR) 172 { 173 if (cleanSingularToricFace(face,sqrt_density) == false) 174 { 175 done = false; 176 } 177 } 178 else 179 { 180 if (cleanToricFace(face,sqrt_density) == false) 181 { 182 done = false; 183 } 184 } 185 } 186 } 187 } 188 } 189 190 cleanVertices(); 191 cleanEdges(); 192 cleanContactFaces(); 193 cleanToricFaces(); 194 cleanSphericFaces(); 195 } 196 197 cleanSingularToricFace(SESFace * face,const double & sqrt_density)198 bool SolventExcludedSurface::cleanSingularToricFace(SESFace* face, const double& sqrt_density) 199 { 200 face->normalize(true); 201 std::list<SESEdge*>::iterator e = face->edge_.begin(); 202 SESEdge* edge0 = *e; 203 e++; 204 e++; 205 e++; 206 SESEdge* edge3 = *e; 207 std::list<SESVertex*>::iterator v = face->vertex_.begin(); 208 SESVertex* v0 = *v; 209 v++; 210 v++; 211 SESVertex* v2 = *v; 212 v++; 213 SESVertex* v3 = *v; 214 v++; 215 v++; 216 SESVertex* v5 = *v; 217 bool del = false; 218 bool set = false; 219 SESEdge* edge = NULL; 220 if (v0 == v2) 221 { 222 del = (edge0->rsedge_->angle_.value < Constants::PI); 223 set = !del; 224 edge = edge0; 225 } 226 else 227 { 228 if (v3 == v5) 229 { 230 del = (edge3->rsedge_->angle_.value < Constants::PI); 231 set = !del; 232 edge = edge3; 233 } 234 else 235 { 236 double exact_number_of_segments(face->rsedge_->angle_.value* 237 edge3->circle_.radius* 238 sqrt_density); 239 del = (Maths::isLess(exact_number_of_segments,0.1)); 240 } 241 } 242 if (del) 243 { 244 deleteSmallSingularToricFace(face); 245 } 246 if (set) 247 { 248 edge->rsedge_->angle_.value = 2*Constants::PI; 249 } 250 return !del; 251 } 252 253 cleanToricFace(SESFace * face,const double & sqrt_density)254 bool SolventExcludedSurface::cleanToricFace 255 (SESFace* face, 256 const double& sqrt_density) 257 258 { 259 face->normalize(false); 260 std::list<SESEdge*>::iterator e = face->edge_.begin(); 261 e++; 262 SESEdge* edge1 = *e; 263 e++; 264 e++; 265 SESEdge* edge3 = *e; 266 std::list<SESVertex*>::iterator v = face->vertex_.begin(); 267 SESVertex* v0 = *v; 268 v++; 269 SESVertex* v1 = *v; 270 v++; 271 SESVertex* v2 = *v; 272 v++; 273 SESVertex* v3 = *v; 274 bool del = false; 275 bool set = false; 276 SESEdge* edge = NULL; 277 if (v0 == v3) 278 { 279 del = (edge3->rsedge_->angle_.value < Constants::PI); 280 set = !del; 281 edge = edge3; 282 } 283 else 284 { 285 if (v1 == v2) 286 { 287 del = (edge1->rsedge_->angle_.value < Constants::PI); 288 set = !del; 289 edge = edge1; 290 } 291 else 292 { 293 double exact_number_of_segments(face->rsedge_->angle_.value* 294 edge3->circle_.radius* 295 sqrt_density); 296 del = (Maths::isLess(exact_number_of_segments,0.1)); 297 } 298 } 299 if (del) 300 { 301 deleteSmallToricFace(face); 302 } 303 if (set) 304 { 305 edge->rsedge_->angle_.value = 2*Constants::PI; 306 } 307 return !del; 308 } 309 310 cleanVertices()311 void SolventExcludedSurface::cleanVertices() 312 313 { 314 while ((number_of_vertices_ > 0) && 315 (vertices_[number_of_vertices_-1] == NULL)) 316 { 317 vertices_.pop_back(); 318 number_of_vertices_--; 319 } 320 for (Position i = 0; i < number_of_vertices_; i++) 321 { 322 if (vertices_[i] == NULL) 323 { 324 vertices_[i] = vertices_[number_of_vertices_-1]; 325 vertices_[i]->index_ = i; 326 vertices_.pop_back(); 327 number_of_vertices_--; 328 while (vertices_[number_of_vertices_-1] == NULL) 329 { 330 vertices_.pop_back(); 331 number_of_vertices_--; 332 } 333 } 334 } 335 } 336 337 cleanEdges()338 void SolventExcludedSurface::cleanEdges() 339 340 { 341 while ((number_of_edges_ > 0) && 342 (edges_[number_of_edges_-1] == NULL)) 343 { 344 edges_.pop_back(); 345 number_of_edges_--; 346 } 347 for (Position i = 0; i < number_of_edges_; i++) 348 { 349 if (edges_[i] == NULL) 350 { 351 edges_[i] = edges_[number_of_edges_-1]; 352 edges_[i]->index_ = i; 353 edges_.pop_back(); 354 number_of_edges_--; 355 while (edges_[number_of_edges_-1] == NULL) 356 { 357 edges_.pop_back(); 358 number_of_edges_--; 359 } 360 } 361 } 362 } 363 364 cleanContactFaces()365 void SolventExcludedSurface::cleanContactFaces() 366 367 { 368 while ((number_of_contact_faces_ > 0) && 369 (contact_faces_[number_of_contact_faces_-1] == NULL)) 370 { 371 contact_faces_.pop_back(); 372 number_of_contact_faces_--; 373 } 374 for (Position i = 0; i < number_of_contact_faces_; i++) 375 { 376 if (contact_faces_[i] == NULL) 377 { 378 contact_faces_[i] = contact_faces_[number_of_contact_faces_-1]; 379 contact_faces_[i]->index_ = i; 380 contact_faces_.pop_back(); 381 number_of_contact_faces_--; 382 while (contact_faces_[number_of_contact_faces_-1] == NULL) 383 { 384 contact_faces_.pop_back(); 385 number_of_contact_faces_--; 386 } 387 } 388 } 389 } 390 391 cleanToricFaces()392 void SolventExcludedSurface::cleanToricFaces() 393 394 { 395 while ((number_of_toric_faces_ > 0) && 396 (toric_faces_[number_of_toric_faces_-1] == NULL)) 397 { 398 toric_faces_.pop_back(); 399 number_of_toric_faces_--; 400 } 401 for (Position i = 0; i < number_of_toric_faces_; i++) 402 { 403 if (toric_faces_[i] == NULL) 404 { 405 toric_faces_[i] = toric_faces_[number_of_toric_faces_-1]; 406 toric_faces_[i]->index_ = i; 407 toric_faces_.pop_back(); 408 number_of_toric_faces_--; 409 while (toric_faces_[number_of_toric_faces_-1] == NULL) 410 { 411 toric_faces_.pop_back(); 412 number_of_toric_faces_--; 413 } 414 } 415 } 416 } 417 418 cleanSphericFaces()419 void SolventExcludedSurface::cleanSphericFaces() 420 421 { 422 while ((number_of_spheric_faces_ > 0) && 423 (spheric_faces_[number_of_spheric_faces_-1] == NULL)) 424 { 425 spheric_faces_.pop_back(); 426 number_of_spheric_faces_--; 427 } 428 for (Position i = 0; i < number_of_spheric_faces_; i++) 429 { 430 if (spheric_faces_[i] == NULL) 431 { 432 spheric_faces_[i] = spheric_faces_[number_of_spheric_faces_-1]; 433 spheric_faces_[i]->index_ = i; 434 spheric_faces_.pop_back(); 435 number_of_spheric_faces_--; 436 while (spheric_faces_[number_of_spheric_faces_-1] == NULL) 437 { 438 spheric_faces_.pop_back(); 439 number_of_spheric_faces_--; 440 } 441 } 442 } 443 } 444 445 deleteSmallToricFace(SESFace * face)446 void SolventExcludedSurface::deleteSmallToricFace(SESFace* face) 447 { 448 SESEdge* edge[4]; 449 std::list<SESEdge*>::iterator e = face->edge_.begin(); 450 for (Position i = 0; i < 4; i++) 451 { 452 edge[i] = *e; 453 e++; 454 } 455 SESVertex* p[4]; 456 std::list<SESVertex*>::iterator v = face->vertex_.begin(); 457 for (Position i = 0; i < 4; i++) 458 { 459 p[i] = *v; 460 v++; 461 } 462 SESFace* neighbour1 = edge[1]->other(face); 463 SESFace* neighbour2 = edge[2]->other(face); 464 SESFace* neighbour3 = edge[3]->other(face); 465 if (p[0] != p[3]) 466 { 467 p[0]->join(*p[3]); 468 neighbour3->vertex_.remove(p[3]); 469 p[3]->substitute(p[0]); 470 } 471 if (p[1] != p[2]) 472 { 473 p[1]->join(*p[2]); 474 neighbour1->vertex_.remove(p[2]); 475 p[2]->substitute(p[1]); 476 } 477 p[0]->edges_.erase(edge[2]); 478 p[0]->edges_.erase(edge[3]); 479 p[1]->edges_.erase(edge[2]); 480 p[1]->edges_.erase(edge[1]); 481 p[0]->faces_.erase(face); 482 p[1]->faces_.erase(face); 483 edge[0]->substitute(face,neighbour2); 484 neighbour2->substitute(edge[2],edge[0]); 485 if (p[2] != p[1]) 486 { 487 vertices_[p[2]->index_] = NULL; 488 neighbour1->edge_.remove(edge[1]); 489 delete p[2]; 490 } 491 else 492 { 493 p[1]->faces_.erase(neighbour1); 494 contact_faces_[neighbour1->index_] = NULL; 495 delete neighbour1; 496 } 497 if (p[3] != p[0]) 498 { 499 vertices_[p[3]->index_] = NULL; 500 neighbour3->edge_.remove(edge[3]); 501 delete p[3]; 502 } 503 else 504 { 505 p[0]->faces_.erase(neighbour3); 506 contact_faces_[neighbour3->index_] = NULL; 507 delete neighbour3; 508 } 509 edges_[edge[1]->index_] = NULL; 510 delete edge[1]; 511 edges_[edge[2]->index_] = NULL; 512 delete edge[2]; 513 edges_[edge[3]->index_] = NULL; 514 delete edge[3]; 515 toric_faces_[face->index_] = NULL; 516 delete face; 517 edge[0]->type_ = SESEdge::TYPE_SINGULAR; 518 TAngle<double> phi 519 (getOrientedAngle(edge[0]->vertex_[0]->point_-edge[0]->circle_.p, 520 edge[0]->vertex_[1]->point_-edge[0]->circle_.p, 521 edge[0]->circle_.n)); 522 if (phi.value > Constants::PI) 523 { 524 edge[0]->circle_.n.negate(); 525 } 526 singular_edges_.push_back(edge[0]); 527 } 528 529 deleteSmallSingularToricFace(SESFace * face)530 void SolventExcludedSurface::deleteSmallSingularToricFace 531 (SESFace* face) 532 533 { 534 SESEdge* edge[6]; 535 std::list<SESEdge*>::iterator e = face->edge_.begin(); 536 for (Position i = 0; i < 6; i++) 537 { 538 edge[i] = *e; 539 e++; 540 } 541 SESVertex* p[6]; 542 std::list<SESVertex*>::iterator v = face->vertex_.begin(); 543 for (Position i = 0; i < 6; i++) 544 { 545 p[i] = *v; 546 v++; 547 } 548 SESFace* neighbour0 = edge[0]->other(face); 549 SESFace* neighbour2 = edge[2]->other(face); 550 SESFace* neighbour3 = edge[3]->other(face); 551 SESFace* neighbour5 = edge[5]->other(face); 552 if (p[0] != p[2]) 553 { 554 p[0]->join(*p[2]); 555 neighbour0->vertex_.remove(p[2]); 556 p[2]->substitute(p[0]); 557 } 558 if (p[3] != p[5]) 559 { 560 p[3]->join(*p[5]); 561 neighbour3->vertex_.remove(p[5]); 562 p[5]->substitute(p[3]); 563 } 564 p[0]->edges_.erase(edge[0]); 565 p[0]->edges_.erase(edge[2]); 566 p[1]->edges_.erase(edge[2]); 567 p[3]->edges_.erase(edge[3]); 568 p[3]->edges_.erase(edge[5]); 569 p[4]->edges_.erase(edge[5]); 570 p[0]->faces_.erase(face); 571 p[1]->faces_.erase(face); 572 p[3]->faces_.erase(face); 573 p[4]->faces_.erase(face); 574 edge[1]->substitute(face,neighbour2); 575 edge[4]->substitute(face,neighbour5); 576 neighbour2->substitute(edge[2],edge[1]); 577 neighbour5->substitute(edge[5],edge[4]); 578 if (p[2] != p[0]) 579 { 580 vertices_[p[2]->index_] = NULL; 581 neighbour0->edge_.remove(edge[0]); 582 delete p[2]; 583 } 584 else 585 { 586 p[0]->faces_.erase(neighbour0); 587 contact_faces_[neighbour0->index_] = NULL; 588 delete neighbour0; 589 } 590 if (p[3] != p[5]) 591 { 592 vertices_[p[5]->index_] = NULL; 593 neighbour3->edge_.remove(edge[3]); 594 delete p[5]; 595 } 596 else 597 { 598 p[3]->faces_.erase(neighbour3); 599 contact_faces_[neighbour3->index_] = NULL; 600 delete neighbour3; 601 } 602 edges_[edge[0]->index_] = NULL; 603 delete edge[0]; 604 edges_[edge[2]->index_] = NULL; 605 delete edge[2]; 606 edges_[edge[3]->index_] = NULL; 607 delete edge[3]; 608 edges_[edge[5]->index_] = NULL; 609 delete edge[5]; 610 toric_faces_[face->index_] = NULL; 611 delete face; 612 edge[1]->type_ = SESEdge::TYPE_SINGULAR; 613 TAngle<double> phi 614 (getOrientedAngle(edge[1]->vertex_[0]->point_-edge[1]->circle_.p, 615 edge[1]->vertex_[1]->point_-edge[1]->circle_.p, 616 edge[1]->circle_.n)); 617 if (phi.value > Constants::PI) 618 { 619 edge[1]->circle_.n.negate(); 620 } 621 edge[4]->type_ = SESEdge::TYPE_SINGULAR; 622 phi = getOrientedAngle(edge[4]->vertex_[0]->point_-edge[4]->circle_.p, 623 edge[4]->vertex_[1]->point_-edge[4]->circle_.p, 624 edge[4]->circle_.n); 625 if (phi.value > Constants::PI) 626 { 627 edge[4]->circle_.n.negate(); 628 } 629 singular_edges_.push_back(edge[1]); 630 singular_edges_.push_back(edge[4]); 631 } 632 633 compute()634 void SolventExcludedSurface::compute() 635 throw(Exception::GeneralException) 636 { 637 SESComputer sesc(this); 638 sesc.run(); 639 } 640 641 splitSphericFaces()642 void SolventExcludedSurface::splitSphericFaces() 643 { 644 for (Position i = 0; i < number_of_spheric_faces_; i++) 645 { 646 splitSphericFace(i); 647 } 648 } 649 650 splitSphericFace(Position i)651 void SolventExcludedSurface::splitSphericFace(Position i) 652 653 { 654 SESFace* face = spheric_faces_[i]; 655 std::list<SESEdge*>::iterator e = face->edge_.begin(); 656 while (e != face->edge_.end()) 657 { 658 if ((*e)->vertex_[0] == NULL) 659 { 660 return; 661 } 662 e++; 663 } 664 std::list<SESEdge*> contour; 665 std::list<SESVertex*> vertices; 666 e = face->edge_.begin(); 667 SESVertex* start = (*e)->vertex_[0]; 668 SESVertex* next = (*e)->vertex_[1]; 669 SESEdge* current = *e; 670 contour.push_back(current); 671 vertices.push_back(next); 672 while (next != start) 673 { 674 for (e = face->edge_.begin(); e != face->edge_.end(); e++) 675 { 676 if (*e != current) 677 { 678 if ((*e)->vertex_[0] == next) 679 { 680 contour.push_back(*e); 681 next = (*e)->vertex_[1]; 682 vertices.push_back(next); 683 current = *e; 684 } 685 else 686 { 687 if ((*e)->vertex_[1] == next) 688 { 689 contour.push_back(*e); 690 next = (*e)->vertex_[0]; 691 vertices.push_back(next); 692 current = *e; 693 } 694 } 695 } 696 } 697 } 698 if (contour.size() != face->edge_.size()) 699 { 700 SESFace* new_face = new SESFace(*face,true); 701 for (e = contour.begin(); e != contour.end(); e++) 702 { 703 new_face->edge_.remove(*e); 704 } 705 std::list<SESVertex*>::iterator v; 706 for (v = vertices.begin(); v != vertices.end(); v++) 707 { 708 new_face->vertex_.remove(*v); 709 } 710 new_face->index_ = number_of_spheric_faces_; 711 spheric_faces_.push_back(new_face); 712 number_of_spheric_faces_++; 713 face->edge_ = contour; 714 face->vertex_ = vertices; 715 } 716 } 717 718 check()719 bool SolventExcludedSurface::check() 720 721 { 722 for (Position i = 0; i < number_of_vertices_; i++) 723 { 724 if (vertices_[i]->edges_.size() != vertices_[i]->faces_.size()) 725 { 726 return false; 727 } 728 } 729 for (Position i = 0; i < number_of_spheric_faces_; i++) 730 { 731 if (spheric_faces_[i]->edge_.size() != spheric_faces_[i]->vertex_.size()) 732 { 733 Index test = spheric_faces_[i]->edge_.size()- 734 spheric_faces_[i]->vertex_.size(); 735 std::list<SESEdge*>::iterator e = spheric_faces_[i]->edge_.begin(); 736 while (e != spheric_faces_[i]->edge_.end()) 737 { 738 if ((*e)->vertex_[0] == NULL) 739 { 740 test--; 741 } 742 e++; 743 } 744 if (test != 0) 745 { 746 return false; 747 } 748 } 749 } 750 return true; 751 } 752 753 754 SolventExcludedSurface::VertexIterator beginVertex()755 SolventExcludedSurface::beginVertex() 756 757 { 758 return vertices_.begin(); 759 } 760 761 762 SolventExcludedSurface::ConstVertexIterator beginVertex() const763 SolventExcludedSurface::beginVertex() const 764 765 { 766 return vertices_.begin(); 767 } 768 769 770 SolventExcludedSurface::VertexIterator endVertex()771 SolventExcludedSurface::endVertex() 772 773 { 774 return vertices_.end(); 775 } 776 777 778 SolventExcludedSurface::ConstVertexIterator endVertex() const779 SolventExcludedSurface::endVertex() const 780 781 { 782 return vertices_.end(); 783 } 784 785 786 SolventExcludedSurface::EdgeIterator beginEdge()787 SolventExcludedSurface::beginEdge() 788 789 { 790 return edges_.begin(); 791 } 792 793 794 SolventExcludedSurface::ConstEdgeIterator beginEdge() const795 SolventExcludedSurface::beginEdge() const 796 797 { 798 return edges_.begin(); 799 } 800 801 802 SolventExcludedSurface::EdgeIterator endEdge()803 SolventExcludedSurface::endEdge() 804 805 { 806 return edges_.end(); 807 } 808 809 810 SolventExcludedSurface::ConstEdgeIterator endEdge() const811 SolventExcludedSurface::endEdge() const 812 813 { 814 return edges_.end(); 815 } 816 817 818 SolventExcludedSurface::SingularEdgeIterator beginSingularEdge()819 SolventExcludedSurface::beginSingularEdge() 820 821 { 822 return singular_edges_.begin(); 823 } 824 825 826 SolventExcludedSurface::ConstSingularEdgeIterator beginSingularEdge() const827 SolventExcludedSurface::beginSingularEdge() const 828 829 { 830 return singular_edges_.begin(); 831 } 832 833 834 SolventExcludedSurface::SingularEdgeIterator endSingularEdge()835 SolventExcludedSurface::endSingularEdge() 836 837 { 838 return singular_edges_.end(); 839 } 840 841 842 SolventExcludedSurface::ConstSingularEdgeIterator endSingularEdge() const843 SolventExcludedSurface::endSingularEdge() const 844 845 { 846 return singular_edges_.end(); 847 } 848 849 850 SolventExcludedSurface::ContactFaceIterator beginContactFace()851 SolventExcludedSurface::beginContactFace() 852 853 { 854 return contact_faces_.begin(); 855 } 856 857 858 SolventExcludedSurface::ConstContactFaceIterator beginContactFace() const859 SolventExcludedSurface::beginContactFace() const 860 861 { 862 return contact_faces_.begin(); 863 } 864 865 866 SolventExcludedSurface::ContactFaceIterator endContactFace()867 SolventExcludedSurface::endContactFace() 868 869 { 870 return contact_faces_.end(); 871 } 872 873 874 SolventExcludedSurface::ConstContactFaceIterator endContactFace() const875 SolventExcludedSurface::endContactFace() const 876 877 { 878 return contact_faces_.end(); 879 } 880 881 882 SolventExcludedSurface::SphericFaceIterator beginSphericFace()883 SolventExcludedSurface::beginSphericFace() 884 885 { 886 return spheric_faces_.begin(); 887 } 888 889 890 SolventExcludedSurface::ConstSphericFaceIterator beginSphericFace() const891 SolventExcludedSurface::beginSphericFace() const 892 893 { 894 return spheric_faces_.begin(); 895 } 896 897 898 SolventExcludedSurface::SphericFaceIterator endSphericFace()899 SolventExcludedSurface::endSphericFace() 900 901 { 902 return spheric_faces_.end(); 903 } 904 905 906 SolventExcludedSurface::ConstSphericFaceIterator endSphericFace() const907 SolventExcludedSurface::endSphericFace() const 908 909 { 910 return spheric_faces_.end(); 911 } 912 913 914 SolventExcludedSurface::ToricFaceIterator beginToricFace()915 SolventExcludedSurface::beginToricFace() 916 917 { 918 return toric_faces_.begin(); 919 } 920 921 922 SolventExcludedSurface::ConstToricFaceIterator beginToricFace() const923 SolventExcludedSurface::beginToricFace() const 924 925 { 926 return toric_faces_.begin(); 927 } 928 929 930 SolventExcludedSurface::ToricFaceIterator endToricFace()931 SolventExcludedSurface::endToricFace() 932 933 { 934 return toric_faces_.end(); 935 } 936 937 938 SolventExcludedSurface::ConstToricFaceIterator endToricFace() const939 SolventExcludedSurface::endToricFace() const 940 941 { 942 return toric_faces_.end(); 943 } 944 945 operator <<(std::ostream & s,const SolventExcludedSurface & ses)946 std::ostream& operator << (std::ostream& s, const SolventExcludedSurface& ses) 947 { 948 s << "Vertices:\n"; 949 SolventExcludedSurface::ConstVertexIterator v; 950 for (v = ses.beginVertex(); v != ses.endVertex(); v++) 951 { 952 if (*v != NULL) 953 { 954 s << " " << **v << "\n"; 955 } 956 else 957 { 958 s << " --\n"; 959 } 960 } 961 s << "Edges:\n"; 962 SolventExcludedSurface::ConstEdgeIterator e; 963 for (e = ses.beginEdge(); e != ses.endEdge(); e++) 964 { 965 if (*e != NULL) 966 { 967 s << " " << **e << "\n"; 968 } 969 else 970 { 971 s << " --\n"; 972 } 973 } 974 s << "singular Edges:\n"; 975 SolventExcludedSurface::ConstSingularEdgeIterator se; 976 for (se = ses.beginSingularEdge(); se != ses.endSingularEdge(); se++) 977 { 978 if (*se != NULL) 979 { 980 s << " " << **se << "\n"; 981 } 982 else 983 { 984 s << " --\n"; 985 } 986 } 987 s << "contact Faces:\n"; 988 SolventExcludedSurface::ConstContactFaceIterator cf; 989 for (cf = ses.beginContactFace(); cf != ses.endContactFace(); cf++) 990 { 991 if (*cf != NULL) 992 { 993 s << " " << **cf << "\n"; 994 } 995 else 996 { 997 s << " --\n"; 998 } 999 } 1000 s << "toric Faces:\n"; 1001 SolventExcludedSurface::ConstToricFaceIterator tf; 1002 for (tf = ses.beginToricFace(); tf != ses.endToricFace(); tf++) 1003 { 1004 if (*tf != NULL) 1005 { 1006 s << " " << **tf << "\n"; 1007 } 1008 else 1009 { 1010 s << " --\n"; 1011 } 1012 } 1013 s << "spheric Faces:\n"; 1014 SolventExcludedSurface::ConstSphericFaceIterator sf; 1015 for (sf = ses.beginSphericFace(); sf != ses.endSphericFace(); sf++) 1016 { 1017 if (*sf != NULL) 1018 { 1019 s << " " << **sf << "\n"; 1020 } 1021 else 1022 { 1023 s << " --\n"; 1024 } 1025 } 1026 return s; 1027 } 1028 1029 /////////////////////////////// 1030 1031 SESComputer()1032 SESComputer::SESComputer() 1033 : ses_(), 1034 vertex_grid_() 1035 { 1036 } 1037 1038 SESComputer(SolventExcludedSurface * ses)1039 SESComputer::SESComputer(SolventExcludedSurface* ses) 1040 : ses_(ses), 1041 vertex_grid_() 1042 { 1043 } 1044 1045 ~SESComputer()1046 SESComputer::~SESComputer() 1047 { 1048 } 1049 1050 run()1051 void SESComputer::run() 1052 throw(Exception::GeneralException) 1053 { 1054 preProcessing(); 1055 get(); 1056 SESSingularityCleaner sessc(ses_,&vertex_grid_); 1057 while (!sessc.run()) 1058 { 1059 ses_->clear(); 1060 vertex_grid_.clear(); 1061 preProcessing(); 1062 get(); 1063 sessc.vertex_grid_ = &vertex_grid_; 1064 } 1065 } 1066 1067 preProcessing()1068 void SESComputer::preProcessing() 1069 1070 { 1071 ses_->clear(); 1072 ses_->number_of_contact_faces_ 1073 = ses_->reduced_surface_->number_of_vertices_; 1074 ses_->number_of_toric_faces_ 1075 = ses_->reduced_surface_->number_of_edges_; 1076 ses_->number_of_spheric_faces_ 1077 = ses_->reduced_surface_->number_of_faces_; 1078 SESFace* face; 1079 Position i; 1080 for (i = 0; i < ses_->number_of_contact_faces_; i++) 1081 { 1082 face = new SESFace; 1083 face->type_ = SESFace::TYPE_CONTACT; 1084 face->rsvertex_ = ses_->reduced_surface_->vertices_[i]; 1085 face->rsedge_ = NULL; 1086 face->rsface_ = NULL; 1087 face->index_ = i; 1088 ses_->contact_faces_.push_back(face); 1089 } 1090 for (i = 0; i < ses_->number_of_toric_faces_; i++) 1091 { 1092 face = new SESFace; 1093 face->type_ = SESFace::TYPE_TORIC; 1094 face->rsvertex_ = NULL; 1095 face->rsedge_ = ses_->reduced_surface_->edges_[i]; 1096 face->rsface_ = NULL; 1097 face->index_ = i; 1098 ses_->toric_faces_.push_back(face); 1099 } 1100 for (i = 0; i < ses_->number_of_spheric_faces_; i++) 1101 { 1102 face = new SESFace; 1103 face->type_ = SESFace::TYPE_SPHERIC; 1104 face->rsvertex_ = NULL; 1105 face->rsedge_ = NULL; 1106 face->rsface_ = ses_->reduced_surface_->faces_[i]; 1107 face->index_ = i; 1108 ses_->spheric_faces_.push_back(face); 1109 } 1110 double x_min = ses_->reduced_surface_->bounding_box_.a.x; 1111 double y_min = ses_->reduced_surface_->bounding_box_.a.y; 1112 double z_min = ses_->reduced_surface_->bounding_box_.a.z; 1113 double x_max = ses_->reduced_surface_->bounding_box_.b.x; 1114 double y_max = ses_->reduced_surface_->bounding_box_.b.y; 1115 double z_max = ses_->reduced_surface_->bounding_box_.b.z; 1116 double dist = ses_->reduced_surface_->r_max_; 1117 Position nx = (Position)((x_max-x_min)/dist+5); 1118 Position ny = (Position)((y_max-y_min)/dist+5); 1119 Position nz = (Position)((z_max-z_min)/dist+5); 1120 Vector3 origin(x_min-2*dist,y_min-2*dist,z_min-2*dist); 1121 vertex_grid_.set(origin, Vector3(dist), nx, ny, nz); 1122 } 1123 1124 get()1125 void SESComputer::get() 1126 { 1127 for (Position i = 0; i < ses_->number_of_spheric_faces_; i++) 1128 { 1129 createSphericFace(i); 1130 } 1131 for (Position i = 0; i < ses_->number_of_toric_faces_; i++) 1132 { 1133 createToricFace(i); 1134 } 1135 } 1136 1137 createSphericFace(Position j)1138 void SESComputer::createSphericFace(Position j) 1139 { 1140 SESFace* face = ses_->spheric_faces_[j]; 1141 RSFace* rsface = face->rsface_; 1142 double probe_radius = ses_->reduced_surface_->probe_radius_; 1143 TSphere3<double> probe(rsface->center_,probe_radius); 1144 // create three vertices and push them to their faces 1145 // and in the list of vertices 1146 pushVertex(face,probe,rsface->vertex_[0]); 1147 pushVertex(face,probe,rsface->vertex_[1]); 1148 pushVertex(face,probe,rsface->vertex_[2]); 1149 // create three concace edges and push them to their faces 1150 // and in the list of edges 1151 pushConcaveEdge(face,0,1,probe_radius); 1152 pushConcaveEdge(face,1,2,probe_radius); 1153 pushConcaveEdge(face,2,0,probe_radius); 1154 } 1155 1156 pushVertex(SESFace * face,const TSphere3<double> & probe,RSVertex * rsvertex)1157 void SESComputer::pushVertex 1158 (SESFace* face, 1159 const TSphere3<double>& probe, 1160 RSVertex* rsvertex) 1161 { 1162 // Create a new vertex on the correct position ... 1163 SESVertex* vertex(createVertex(probe.p,rsvertex->atom_)); 1164 // ... and push it to the face's vertices. 1165 face->vertex_.push_back(vertex); 1166 vertex->faces_.insert(face); 1167 // Get the RSEdges of the corresponding RSVertex ... 1168 RSEdge* tf1(0); 1169 RSEdge* tf2(0); 1170 face->rsface_->getEdges(rsvertex,tf1,tf2); 1171 // ... and push the vertex to these toric faces. 1172 ses_->toric_faces_[tf1->index_]->vertex_.push_back(vertex); 1173 vertex->faces_.insert(ses_->toric_faces_[tf1->index_]); 1174 ses_->toric_faces_[tf2->index_]->vertex_.push_back(vertex); 1175 vertex->faces_.insert(ses_->toric_faces_[tf2->index_]); 1176 // Push the vertex to the contact face of the corresponding RSVertex ... 1177 ses_->contact_faces_[rsvertex->index_]->vertex_.push_back(vertex); 1178 vertex->faces_.insert(ses_->contact_faces_[rsvertex->index_]); 1179 // ... and to the vertices of the SES. 1180 ses_->vertices_.push_back(vertex); 1181 Vector3 pos(vertex->point_.x, 1182 vertex->point_.y, 1183 vertex->point_.z); 1184 vertex_grid_.insert(pos,vertex->index_); 1185 ses_->number_of_vertices_++; 1186 } 1187 1188 pushConcaveEdge(SESFace * face,Position p1,Position p2,const double & radius_of_probe)1189 void SESComputer::pushConcaveEdge 1190 (SESFace* face, 1191 Position p1, 1192 Position p2, 1193 const double& radius_of_probe) 1194 1195 { 1196 RSFace* rsface = face->rsface_; 1197 RSEdge* rsedge(0); 1198 // get the corresponding RSEdge 1199 rsface->getEdge(rsface->getVertex(p1),rsface->getVertex(p2),rsedge); 1200 Index index = rsedge->index_; 1201 // create a new SESEdge 1202 SESEdge* edge(createConcaveEdge(face,p1,p2,index,radius_of_probe)); 1203 // and push it to all it's places 1204 face->edge_.push_back(edge); 1205 ses_->toric_faces_[index]->edge_.push_back(edge); 1206 edge->vertex_[0]->edges_.insert(edge); 1207 edge->vertex_[1]->edges_.insert(edge); 1208 ses_->edges_.push_back(edge); 1209 ses_->number_of_edges_++; 1210 } 1211 1212 createVertex(const TVector3<double> & probe_center,Index index)1213 SESVertex* SESComputer::createVertex 1214 (const TVector3<double>& probe_center, 1215 Index index) 1216 1217 { 1218 SESVertex* vertex = new SESVertex; 1219 TSphere3<double>* atom = &(ses_->reduced_surface_->atom_[index]); 1220 // get the position of the new vertex 1221 getPoint(atom->p,probe_center,atom->radius,vertex->point_); 1222 vertex->normal_.set(probe_center-vertex->point_); 1223 vertex->atom_ = index; 1224 vertex->index_ = ses_->number_of_vertices_; 1225 return vertex; 1226 } 1227 1228 createConcaveEdge(SESFace * face,Position p1,Position p2,Index index,const double & radius_of_probe)1229 SESEdge* SESComputer::createConcaveEdge 1230 (SESFace* face, 1231 Position p1, 1232 Position p2, 1233 Index index, 1234 const double& radius_of_probe) 1235 1236 { 1237 SESEdge* edge = new SESEdge; 1238 // set the vertices of the edge 1239 std::list<SESVertex*>::iterator v = face->vertex_.begin(); 1240 for (Position i = 0; i < p1; i++) 1241 { 1242 v++; 1243 } 1244 edge->vertex_[0] = *v; 1245 v = face->vertex_.begin(); 1246 for (Position i = 0; i < p2; i++) 1247 { 1248 v++; 1249 } 1250 edge->vertex_[1] = *v; 1251 // set the faces of the edge 1252 edge->face_[0] = face; 1253 edge->face_[1] = ses_->toric_faces_[index]; 1254 // 1255 edge->rsedge_ = NULL; 1256 edge->type_ = SESEdge::TYPE_CONCAVE; 1257 edge->index_ = ses_->number_of_edges_; 1258 // compute the circle on which the edge lies 1259 RSFace* rsface = face->rsface_; 1260 TVector3<double> normal = (edge->vertex_[0]->point_-rsface->center_)% 1261 (edge->vertex_[1]->point_-rsface->center_); 1262 edge->circle_.set(rsface->center_,normal,radius_of_probe); 1263 return edge; 1264 } 1265 1266 createToricFace(Position i)1267 void SESComputer::createToricFace(Position i) 1268 1269 { 1270 SESFace* face = ses_->toric_faces_[i]; 1271 if (face->isFree()) 1272 { 1273 createFreeToricFace(i); 1274 } 1275 else 1276 { 1277 SESEdge* edge1 = createConvexEdge(face,face->rsedge_->vertex_[0]); 1278 SESEdge* edge2 = createConvexEdge(face,face->rsedge_->vertex_[1]); 1279 if (Maths::isEqual(face->rsedge_->angle_.value,Constants::PI)) 1280 { 1281 RSFace* rsface1 = face->rsedge_->face_[0]; 1282 RSFace* rsface2 = face->rsedge_->face_[1]; 1283 RSVertex* rsvertex1 = face->rsedge_->vertex_[0]; 1284 RSVertex* rsvertex2 = face->rsedge_->vertex_[1]; 1285 RSVertex* rsvertex3 = rsface1->third(rsvertex1,rsvertex2); 1286 RSVertex* rsvertex4 = rsface2->third(rsvertex1,rsvertex2); 1287 TVector3<double> point1(ses_->reduced_surface_->atom_[rsvertex1->atom_].p); 1288 TVector3<double> point2(ses_->reduced_surface_->atom_[rsvertex2->atom_].p); 1289 TVector3<double> point3(ses_->reduced_surface_->atom_[rsvertex3->atom_].p); 1290 TVector3<double> point4(ses_->reduced_surface_->atom_[rsvertex4->atom_].p); 1291 TVector3<double> middle(edge1->circle_.n% 1292 (edge1->vertex_[0]->point_-edge1->circle_.p)); 1293 middle.normalize(); 1294 middle *= edge1->circle_.radius; 1295 middle += edge1->circle_.p; 1296 TPlane3<double> plane(point1,point2,rsface1->center_); 1297 double test = plane.n*(point3-plane.p); 1298 if (test*(plane.n*(middle-plane.p)) > 0) 1299 { 1300 edge1->revert(); 1301 } 1302 middle.set(edge2->circle_.n% 1303 (edge2->vertex_[0]->point_-edge2->circle_.p)); 1304 middle.normalize(); 1305 middle *= edge2->circle_.radius; 1306 middle += edge2->circle_.p; 1307 if (test*(plane.n*(middle-plane.p)) > 0) 1308 { 1309 edge2->revert(); 1310 } 1311 } 1312 if (face->rsedge_->singular_) 1313 { 1314 treatSingularToricFace(i); 1315 } 1316 } 1317 } 1318 1319 createConvexEdge(SESFace * face,RSVertex * rsvertex)1320 SESEdge* SESComputer::createConvexEdge(SESFace* face, RSVertex* rsvertex) 1321 { 1322 SESEdge* edge = new SESEdge; 1323 Index atom = rsvertex->atom_; 1324 Index index = rsvertex->index_; 1325 // find the first vertex of the toric face 1326 // which lies on the surface of the rsvertex 1327 std::list<SESVertex*>::iterator v = face->vertex_.begin(); 1328 while ((*v)->atom_ != atom) 1329 { 1330 v++; 1331 } 1332 edge->vertex_[0] = *v; 1333 // find the second vertex of the toric face 1334 // which lies on the surface of the rsvertex 1335 v++; 1336 while ((*v)->atom_ != atom) 1337 { 1338 v++; 1339 } 1340 edge->vertex_[1] = *v; 1341 // set the faces of the edge 1342 edge->face_[0] = face; 1343 edge->face_[1] = ses_->contact_faces_[index]; 1344 // set the rsedge, type and index of the edge 1345 edge->rsedge_ = face->rsedge_; 1346 edge->type_ = SESEdge::TYPE_CONVEX; 1347 edge->index_ = ses_->number_of_edges_; 1348 // compute the circle on which the edge lies 1349 TCircle3<double> circle0(edge->rsedge_->circle0_); 1350 TCircle3<double> circle1(edge->rsedge_->circle1_); 1351 if (edge->rsedge_->vertex_[0]->index_ == index) 1352 { 1353 edge->circle_.set(circle0.p,circle0.p-circle1.p,circle0.radius); 1354 } 1355 else 1356 { 1357 edge->circle_.set(circle1.p,circle1.p-circle0.p,circle1.radius); 1358 } 1359 TVector3<double> v1(edge->vertex_[0]->point_-edge->circle_.p); 1360 TVector3<double> v2(edge->vertex_[1]->point_-edge->circle_.p); 1361 TVector3<double> n(edge->circle_.n); 1362 TAngle<double> test_phi(getOrientedAngle(v1,v2,n)); 1363 if ((test_phi.value-Constants::PI)* 1364 (edge->rsedge_->angle_.value-Constants::PI) < 0) 1365 { 1366 // test_phi smaller than PI, but expected to be greater or 1367 // test_phi greater than PI, but expected to be smaller 1368 edge->revert(); 1369 } 1370 face->edge_.push_back(edge); 1371 ses_->contact_faces_[index]->edge_.push_back(edge); 1372 edge->vertex_[0]->edges_.insert(edge); 1373 edge->vertex_[1]->edges_.insert(edge); 1374 ses_->edges_.push_back(edge); 1375 ses_->number_of_edges_++; 1376 return edge; 1377 } 1378 1379 treatSingularToricFace(Position i)1380 void SESComputer::treatSingularToricFace(Position i) 1381 1382 { 1383 SESFace* face = ses_->toric_faces_[i]; 1384 face->normalize(false); 1385 SESEdge* edge[4]; 1386 std::list<SESEdge*>::iterator e = face->edge_.begin(); 1387 for (Position i = 0; i < 4; i++) 1388 { 1389 edge[i] = *e; 1390 e++; 1391 } 1392 SESVertex* p[4]; 1393 std::list<SESVertex*>::iterator v = face->vertex_.begin(); 1394 for (Position i = 0; i < 4; i++) 1395 { 1396 p[i] = *v; 1397 v++; 1398 } 1399 // compute the circle on which the singular edge lies 1400 SESFace* neighbour0 = edge[0]->other(face); 1401 SESFace* neighbour2 = edge[2]->other(face); 1402 double probe_radius = ses_->reduced_surface_->probe_radius_; 1403 TSphere3<double> probe1(neighbour0->rsface_->center_,probe_radius); 1404 TSphere3<double> probe2(neighbour2->rsface_->center_,probe_radius); 1405 TCircle3<double> intersection_circle; 1406 GetIntersection(probe1,probe2,intersection_circle); 1407 // create the new edges 1408 SESEdge* new_edge0 = new SESEdge(*edge[0],true); 1409 SESEdge* new_edge2 = new SESEdge(*edge[2],true); 1410 SESEdge* new_edge = new SESEdge(NULL,NULL,neighbour0,neighbour2, 1411 intersection_circle,face->rsedge_, 1412 SESEdge::TYPE_SINGULAR,-1); 1413 // create the new points 1414 Position ip1 = ((p[1]->atom_ == face->rsedge_->vertex_[0]->atom_) ? 0 : 1); 1415 Position ip3 = 1-ip1; 1416 SESVertex* new_point1 1417 = createSingularVertex(ip3,intersection_circle.p, 1418 face,neighbour0,neighbour2, 1419 edge[0],edge[2],new_edge); 1420 SESVertex* new_point3 1421 = createSingularVertex(ip1,intersection_circle.p, 1422 face,neighbour0,neighbour2, 1423 new_edge0,new_edge2,new_edge); 1424 // update the new edges 1425 updateEdge(edge[0],p[0],new_point1,false); 1426 updateEdge(edge[2],p[3],new_point1,false); 1427 updateEdge(new_edge0,p[1],new_point3,true); 1428 updateEdge(new_edge2,p[2],new_point3,true); 1429 // update the singular edge 1430 updateEdge(new_edge,new_point3,new_point1,true); 1431 ses_->singular_edges_.push_back(new_edge); 1432 ses_->number_of_singular_edges_++; 1433 // swap normal of new edge if necessary 1434 TAngle<double> phi(getOrientedAngle(new_point1->point_-intersection_circle.p, 1435 new_point3->point_-intersection_circle.p, 1436 intersection_circle.n)); 1437 if ((face->rsedge_->angle_.value-Constants::PI)* 1438 (phi.value -Constants::PI) < 0) 1439 { 1440 new_edge->circle_.n.negate(); 1441 } 1442 // update the neighbour faces 1443 neighbour0->edge_.push_back(new_edge0); 1444 neighbour0->edge_.push_back(new_edge); 1445 neighbour0->vertex_.push_back(new_point1); 1446 neighbour0->vertex_.push_back(new_point3); 1447 neighbour2->edge_.push_back(new_edge2); 1448 neighbour2->edge_.push_back(new_edge); 1449 neighbour2->vertex_.push_back(new_point1); 1450 neighbour2->vertex_.push_back(new_point3); 1451 // update the toric face 1452 face->type_ = SESFace::TYPE_TORIC_SINGULAR; 1453 face->vertex_.push_back(new_point1); 1454 face->vertex_.push_back(new_point3); 1455 face->edge_.push_back(new_edge0); 1456 face->edge_.push_back(new_edge2); 1457 // update the vertices 1458 p[1]->edges_.erase(edge[0]); 1459 p[1]->edges_.insert(new_edge0); 1460 p[2]->edges_.erase(edge[2]); 1461 p[2]->edges_.insert(new_edge2); 1462 } 1463 1464 createSingularVertex(Position ip,const TVector3<double> & dir,SESFace * face0,SESFace * face1,SESFace * face2,SESEdge * edge0,SESEdge * edge1,SESEdge * edge2)1465 SESVertex* SESComputer::createSingularVertex 1466 (Position ip, 1467 const TVector3<double>& dir, 1468 SESFace* face0, 1469 SESFace* face1, 1470 SESFace* face2, 1471 SESEdge* edge0, 1472 SESEdge* edge1, 1473 SESEdge* edge2) 1474 1475 { 1476 SESVertex* vertex(0); 1477 TVector3<double> intersection_point(face0->rsedge_->getIntersectionPoint(ip)); 1478 Index test = vertexExists(intersection_point); 1479 if (test == -1) 1480 { 1481 vertex = new SESVertex(intersection_point, 1482 dir-intersection_point, 1483 face0->rsedge_->getVertex(ip)->atom_, 1484 ses_->number_of_vertices_); 1485 ses_->vertices_.push_back(vertex); 1486 Vector3 pos(vertex->point_.x, 1487 vertex->point_.y, 1488 vertex->point_.z); 1489 vertex_grid_.insert(pos,vertex->index_); 1490 ses_->number_of_vertices_++; 1491 } 1492 else 1493 { 1494 vertex = ses_->vertices_[test]; 1495 } 1496 vertex->edges_.insert(edge0); 1497 vertex->edges_.insert(edge1); 1498 vertex->edges_.insert(edge2); 1499 vertex->faces_.insert(face0); 1500 vertex->faces_.insert(face1); 1501 vertex->faces_.insert(face2); 1502 return vertex; 1503 } 1504 1505 updateEdge(SESEdge * edge,SESVertex * vertex1,SESVertex * vertex2,bool is_new)1506 void SESComputer::updateEdge 1507 (SESEdge* edge, 1508 SESVertex* vertex1, 1509 SESVertex* vertex2, 1510 bool is_new) 1511 1512 { 1513 if (edge->vertex_[0] == vertex1) 1514 { 1515 edge->vertex_[0] = vertex1; 1516 edge->vertex_[1] = vertex2; 1517 } 1518 else 1519 { 1520 edge->vertex_[0] = vertex2; 1521 edge->vertex_[1] = vertex1; 1522 } 1523 if (is_new) 1524 { 1525 edge->index_ = ses_->number_of_edges_; 1526 ses_->edges_.push_back(edge); 1527 ses_->number_of_edges_++; 1528 } 1529 } 1530 1531 createFreeToricFace(Position i)1532 void SESComputer::createFreeToricFace(Position i) 1533 1534 { 1535 SESFace* face(ses_->toric_faces_[i]); 1536 TCircle3<double> circle1(face->rsedge_->circle0_); 1537 TCircle3<double> circle2(face->rsedge_->circle1_); 1538 Index index1(face->rsedge_->vertex_[0]->index_); 1539 Index index2(face->rsedge_->vertex_[1]->index_); 1540 SESEdge* edge = new SESEdge; 1541 edge->type_ = SESEdge::TYPE_CONVEX; 1542 edge->vertex_[0] = NULL; 1543 edge->vertex_[1] = NULL; 1544 edge->rsedge_ = face->rsedge_; 1545 edge->face_[0] = face; 1546 edge->face_[1] = ses_->contact_faces_[index1]; 1547 edge->circle_.set(circle1.p,circle1.p-circle2.p,circle1.radius); 1548 edge->index_ = ses_->number_of_edges_; 1549 face->edge_.push_back(edge); 1550 ses_->contact_faces_[index1]->edge_.push_back(edge); 1551 ses_->edges_.push_back(edge); 1552 ses_->number_of_edges_++; 1553 edge = new SESEdge; 1554 edge->type_ = SESEdge::TYPE_CONVEX; 1555 edge->vertex_[0] = NULL; 1556 edge->vertex_[1] = NULL; 1557 edge->rsedge_ = face->rsedge_; 1558 edge->face_[0] = face; 1559 edge->face_[1] = ses_->contact_faces_[index2]; 1560 edge->circle_.set(circle2.p,circle2.p-circle1.p,circle2.radius); 1561 edge->index_ = ses_->number_of_edges_; 1562 face->edge_.push_back(edge); 1563 ses_->contact_faces_[index2]->edge_.push_back(edge); 1564 ses_->edges_.push_back(edge); 1565 ses_->number_of_edges_++; 1566 } 1567 1568 getPoint(const TVector3<double> & p1,const TVector3<double> & p2,const double & dist,TVector3<double> & result)1569 void SESComputer::getPoint 1570 (const TVector3<double>& p1, 1571 const TVector3<double>& p2, 1572 const double& dist, 1573 TVector3<double>& result) 1574 1575 { 1576 result.set(p2-p1); 1577 result.normalize(); 1578 result *= dist; 1579 result += p1; 1580 } 1581 1582 vertexExists(const TVector3<double> & point)1583 Index SESComputer::vertexExists(const TVector3<double>& point) 1584 1585 { 1586 double epsilon = Constants::EPSILON; 1587 Constants::EPSILON = 0.001; 1588 Vector3 p(point.x,point.y,point.z); 1589 HashGridBox3<Index>* box = vertex_grid_.getBox(p); 1590 HashGridBox3<Index>::ConstBoxIterator b; 1591 HashGridBox3<Index>::ConstDataIterator d; 1592 if (box != NULL) 1593 { 1594 for (b = box->beginBox(); b != box->endBox(); b++) 1595 { 1596 for (d = b->beginData(); d != b->endData(); d++) 1597 { 1598 if (ses_->vertices_[*d]->point_ == point) 1599 { 1600 Constants::EPSILON = epsilon; 1601 return *d; 1602 } 1603 } 1604 } 1605 } 1606 Constants::EPSILON = epsilon; 1607 return -1; 1608 } 1609 1610 1611 1612 /////////////////////////////// 1613 1614 SESSingularityCleaner()1615 SESSingularityCleaner::SESSingularityCleaner() 1616 : ses_(), 1617 vertex_grid_(), 1618 probe_intersections_() 1619 { 1620 } 1621 1622 SESSingularityCleaner(SolventExcludedSurface * ses,HashGrid3<Index> * vertex_grid)1623 SESSingularityCleaner::SESSingularityCleaner(SolventExcludedSurface* ses, HashGrid3<Index>* vertex_grid) 1624 : ses_(ses), 1625 vertex_grid_(vertex_grid), 1626 probe_intersections_() 1627 { 1628 } 1629 1630 ~SESSingularityCleaner()1631 SESSingularityCleaner::~SESSingularityCleaner() 1632 { 1633 // delete probe_intersections 1634 HashMap< Position, 1635 HashMap< Position, 1636 HashMap< Position, 1637 ProbeIntersection* > > >::Iterator pi1; 1638 HashMap< Position, 1639 HashMap< Position, 1640 ProbeIntersection* > >::Iterator pi2; 1641 HashMap< Position,ProbeIntersection* >::Iterator pi3; 1642 for (pi1 = probe_intersections_.begin(); 1643 pi1 != probe_intersections_.end(); 1644 pi1++) 1645 { 1646 for (pi2 = pi1->second.begin(); pi2 != pi1->second.end(); pi2++) 1647 { 1648 for (pi3 = pi2->second.begin(); pi3 != pi2->second.end(); pi3++) 1649 { 1650 delete pi3->second; 1651 } 1652 } 1653 } 1654 } 1655 1656 run()1657 bool SESSingularityCleaner::run() 1658 throw(Exception::GeneralException) 1659 { 1660 if (!treatFirstCategory()) 1661 { 1662 return false; 1663 } 1664 if (ses_->number_of_singular_edges_ > 0) 1665 { 1666 treatSecondCategory(); 1667 } 1668 return true; 1669 } 1670 1671 treatFirstCategory()1672 bool SESSingularityCleaner::treatFirstCategory() 1673 { 1674 std::list<SESFace*> first_category_faces; 1675 getFirstCategoryFaces(first_category_faces); 1676 1677 SESFace* face1(0); 1678 SESFace* face2(0); 1679 bool modified = false; 1680 std::list<SESFace*>::iterator f 1681 = first_category_faces.begin(); 1682 while (f != first_category_faces.end()) 1683 { 1684 face1 = *f; 1685 f++; 1686 face2 = *f; 1687 f++; 1688 switch (face1->edge_.size()) 1689 { 1690 case 3 : noCut(face1,face2); 1691 break; 1692 case 5 : break; 1693 case 7 : twoCuts(face1,face2); 1694 break; 1695 case 9 : ses_->reduced_surface_->deleteSimilarFaces(face1->rsface_, 1696 face2->rsface_); 1697 modified = true; 1698 break; 1699 } 1700 } 1701 if (modified) 1702 { 1703 ses_->reduced_surface_->clean(); 1704 return false; 1705 } 1706 else 1707 { 1708 return true; 1709 } 1710 } 1711 1712 getFirstCategoryFaces(std::list<SESFace * > & first_category_faces)1713 void SESSingularityCleaner::getFirstCategoryFaces(std::list<SESFace*>& first_category_faces) 1714 { 1715 std::list<SESFace*> singular_faces; 1716 getSingularFaces(singular_faces); 1717 1718 while (!singular_faces.empty()) 1719 { 1720 SESFace* current = singular_faces.front(); 1721 singular_faces.pop_front(); 1722 std::list<SESFace*>::iterator i = singular_faces.begin(); 1723 while (i != singular_faces.end()) 1724 { 1725 if (*current->rsface_ *= *((*i)->rsface_)) 1726 { 1727 first_category_faces.push_back(current); 1728 first_category_faces.push_back(*i); 1729 singular_faces.erase(i); 1730 1731 break; 1732 } 1733 else 1734 { 1735 i++; 1736 } 1737 } 1738 } 1739 } 1740 1741 getSingularFaces(std::list<SESFace * > & faces)1742 void SESSingularityCleaner::getSingularFaces(std::list<SESFace*>& faces) 1743 { 1744 for (Position i = 0; i < ses_->number_of_spheric_faces_; i++) 1745 { 1746 if (ses_->spheric_faces_[i]->rsface_->singular_ == true) 1747 { 1748 faces.push_back(ses_->spheric_faces_[i]); 1749 } 1750 } 1751 } 1752 1753 noCut(SESFace * face1,SESFace * face2)1754 void SESSingularityCleaner::noCut(SESFace* face1, SESFace* face2) 1755 { 1756 TCircle3<double> circle; 1757 double probe_radius = ses_->reduced_surface_->probe_radius_; 1758 TSphere3<double> s1(face1->rsface_->center_,probe_radius); 1759 TSphere3<double> s2(face2->rsface_->center_,probe_radius); 1760 GetIntersection(s1,s2,circle); 1761 // test whether the circle is really an edge 1762 TVector3<double> normal(face1->rsface_->normal_); 1763 TVector3<double> point1(ses_->reduced_surface_ 1764 ->atom_[face1->rsface_->vertex_[0]->atom_].p); 1765 TVector3<double> point2(ses_->reduced_surface_ 1766 ->atom_[face1->rsface_->vertex_[1]->atom_].p); 1767 TVector3<double> point3(ses_->reduced_surface_ 1768 ->atom_[face1->rsface_->vertex_[2]->atom_].p); 1769 TVector3<double> u(normal%(point1-point2)); 1770 TVector3<double> v(normal%(point2-point3)); 1771 TVector3<double> w(normal%(point3-point1)); 1772 TVector3<double> diff1(point1-circle.p); 1773 TVector3<double> diff2(point2-circle.p); 1774 double test1 = u*diff1; 1775 double test2 = v*diff2; 1776 double test3 = w*diff1; 1777 if ((Maths::isLess(test1,0.0) && 1778 Maths::isLess(test2,0.0) && 1779 Maths::isLess(test3,0.0) ) || 1780 (Maths::isGreater(test1,0.0) && 1781 Maths::isGreater(test2,0.0) && 1782 Maths::isGreater(test3,0.0) ) ) 1783 { 1784 SESEdge* edge 1785 = new SESEdge(NULL,NULL,face1,face2,circle,NULL, 1786 SESEdge::TYPE_SINGULAR, 1787 ses_->number_of_edges_); 1788 ses_->edges_.push_back(edge); 1789 ses_->singular_edges_.push_back(edge); 1790 ses_->number_of_edges_++; 1791 face1->edge_.push_back(edge); 1792 face2->edge_.push_back(edge); 1793 } 1794 } 1795 1796 twoCuts(SESFace * face1,SESFace * face2)1797 void SESSingularityCleaner::twoCuts 1798 (SESFace* face1, SESFace* face2) 1799 1800 { 1801 std::vector<SESEdge*> sesedge1(7); 1802 std::vector<SESEdge*> sesedge2(7); 1803 std::vector<SESVertex*> sesvertex1(7); 1804 std::vector<SESVertex*> sesvertex2(7); 1805 sort(face1,face2,sesedge1,sesedge2,sesvertex1,sesvertex2); 1806 TCircle3<double> circle; 1807 TSphere3<double> sphere1(face1->rsface_->center_, 1808 ses_->reduced_surface_->probe_radius_); 1809 TSphere3<double> sphere2(face2->rsface_->center_, 1810 ses_->reduced_surface_->probe_radius_); 1811 GetIntersection(sphere1,sphere2,circle); 1812 TAngle<double> phi(getOrientedAngle(sesvertex1[0]->point_-circle.p, 1813 sesvertex1[2]->point_-circle.p, 1814 circle.n)); 1815 if (phi.value > Constants::PI) 1816 { 1817 circle.n.negate(); 1818 } 1819 SESEdge* new_edge1 1820 = new SESEdge(sesvertex1[0],sesvertex1[2],face1,face2, 1821 circle,NULL,SESEdge::TYPE_SINGULAR, 1822 ses_->number_of_edges_); 1823 ses_->edges_.push_back(new_edge1); 1824 ses_->singular_edges_.push_back(new_edge1); 1825 face1->edge_.push_back(new_edge1); 1826 face2->edge_.push_back(new_edge1); 1827 sesvertex1[0]->edges_.insert(new_edge1); 1828 sesvertex1[2]->edges_.insert(new_edge1); 1829 ses_->number_of_edges_++; 1830 ses_->number_of_singular_edges_++; 1831 SESEdge* new_edge2 1832 = new SESEdge(sesvertex1[3],sesvertex1[6],face1,face2, 1833 circle,NULL,SESEdge::TYPE_SINGULAR, 1834 ses_->number_of_edges_); 1835 ses_->edges_.push_back(new_edge2); 1836 ses_->singular_edges_.push_back(new_edge2); 1837 face1->edge_.push_back(new_edge2); 1838 face2->edge_.push_back(new_edge2); 1839 sesvertex1[3]->edges_.insert(new_edge2); 1840 sesvertex1[6]->edges_.insert(new_edge2); 1841 ses_->number_of_edges_++; 1842 ses_->number_of_singular_edges_++; 1843 if (sesedge1[2] == sesedge2[2]) 1844 { 1845 ses_->edges_[sesedge1[2]->index_] = NULL; 1846 ses_->singular_edges_.remove(sesedge1[2]); 1847 sesvertex1[2]->edges_.erase(sesedge1[2]); 1848 sesvertex1[3]->edges_.erase(sesedge1[2]); 1849 face1->edge_.remove(sesedge1[2]); 1850 face2->edge_.remove(sesedge1[2]); 1851 delete sesedge1[2]; 1852 } 1853 if (sesedge1[6] == sesedge2[6]) 1854 { 1855 ses_->edges_[sesedge1[6]->index_] = NULL; 1856 ses_->singular_edges_.remove(sesedge1[6]); 1857 sesvertex1[6]->edges_.erase(sesedge1[6]); 1858 sesvertex1[0]->edges_.erase(sesedge1[6]); 1859 face1->edge_.remove(sesedge1[6]); 1860 face2->edge_.remove(sesedge1[6]); 1861 delete sesedge1[6]; 1862 } 1863 } 1864 1865 sort(SESFace * face1,SESFace * face2,std::vector<SESEdge * > & sesedge1,std::vector<SESEdge * > & sesedge2,std::vector<SESVertex * > & sesvertex1,std::vector<SESVertex * > & sesvertex2)1866 void SESSingularityCleaner::sort 1867 (SESFace* face1, 1868 SESFace* face2, 1869 std::vector<SESEdge*>& sesedge1, 1870 std::vector<SESEdge*>& sesedge2, 1871 std::vector<SESVertex*>& sesvertex1, 1872 std::vector<SESVertex*>& sesvertex2) 1873 1874 { 1875 // find two equal vertices 1876 std::list<SESVertex*>::iterator v1 = face1->vertex_.begin(); 1877 std::list<SESVertex*>::iterator v2; 1878 bool found = false; 1879 while (!found) 1880 { 1881 v2 = face2->vertex_.begin(); 1882 while (!found && (v2 != face2->vertex_.end())) 1883 { 1884 if (*v2 == *v1) 1885 { 1886 sesvertex1[0] = *v1; 1887 sesvertex2[0] = *v2; 1888 found = true; 1889 } 1890 v2++; 1891 } 1892 v1++; 1893 } 1894 // find first corresponding edges 1895 face1->getEdges(sesvertex1[0],sesedge1[0],sesedge1[1]); 1896 face2->getEdges(sesvertex2[0],sesedge2[0],sesedge2[1]); 1897 if (*sesedge1[0] == *sesedge2[1]) 1898 { 1899 sesedge1[0] = sesedge1[1]; 1900 } 1901 else 1902 { 1903 if (*sesedge1[1] == *sesedge2[0]) 1904 { 1905 sesedge2[0] = sesedge2[1]; 1906 } 1907 else 1908 { 1909 if (*sesedge1[0] == *sesedge2[0]) 1910 { 1911 sesedge1[0] = sesedge1[1]; 1912 sesedge2[0] = sesedge2[1]; 1913 } 1914 } 1915 } 1916 // find remaining edges and vertices 1917 SESEdge* sesedge(0); 1918 sesvertex1[1] = sesedge1[0]->other(sesvertex1[0]); 1919 sesvertex2[1] = sesedge2[0]->other(sesvertex2[0]); 1920 for (Position i = 1; i < 6; i++) 1921 { 1922 face1->getEdges(sesvertex1[i],sesedge1[i],sesedge); 1923 if (sesedge != sesedge1[i-1]) 1924 { 1925 sesedge1[i] = sesedge; 1926 } 1927 face2->getEdges(sesvertex2[i],sesedge2[i],sesedge); 1928 if (sesedge != sesedge2[i-1]) 1929 { 1930 sesedge2[i] = sesedge; 1931 } 1932 sesvertex1[i+1] = sesedge1[i]->other(sesvertex1[i]); 1933 sesvertex2[i+1] = sesedge2[i]->other(sesvertex2[i]); 1934 } 1935 face1->getEdge(sesvertex1[0],sesvertex1[6],sesedge1[6]); 1936 face2->getEdge(sesvertex2[0],sesvertex2[6],sesedge2[6]); 1937 // 1938 SESVertex* sesvertex(0); 1939 if (sesvertex1[2] != sesvertex2[2]) 1940 { 1941 for (Position i = 0; i < 3; i++) 1942 { 1943 sesvertex = sesvertex1[i]; 1944 sesvertex1[i] = sesvertex1[6-i]; 1945 sesvertex1[6-i] = sesvertex; 1946 sesvertex = sesvertex2[i]; 1947 sesvertex2[i] = sesvertex2[6-i]; 1948 sesvertex2[6-i] = sesvertex; 1949 sesedge = sesedge1[i]; 1950 sesedge1[i] = sesedge1[5-i]; 1951 sesedge1[5-i] = sesedge; 1952 sesedge = sesedge2[i]; 1953 sesedge2[i] = sesedge2[5-i]; 1954 sesedge2[5-i] = sesedge; 1955 } 1956 } 1957 } 1958 1959 treatSecondCategory()1960 void SESSingularityCleaner::treatSecondCategory() 1961 { 1962 double x_min = ses_->spheric_faces_[0]->rsface_->center_.x; 1963 double y_min = ses_->spheric_faces_[0]->rsface_->center_.y; 1964 double z_min = ses_->spheric_faces_[0]->rsface_->center_.z; 1965 double x_max = ses_->spheric_faces_[0]->rsface_->center_.x; 1966 double y_max = ses_->spheric_faces_[0]->rsface_->center_.y; 1967 double z_max = ses_->spheric_faces_[0]->rsface_->center_.z; 1968 1969 for (Position i = 1; i != ses_->number_of_spheric_faces_; i++) 1970 { 1971 x_min = std::min(x_min, ses_->spheric_faces_[i]->rsface_->center_.x); 1972 y_min = std::min(y_min, ses_->spheric_faces_[i]->rsface_->center_.y); 1973 z_min = std::min(z_min, ses_->spheric_faces_[i]->rsface_->center_.z); 1974 1975 x_max = std::max(x_max, ses_->spheric_faces_[i]->rsface_->center_.x); 1976 y_max = std::max(y_max, ses_->spheric_faces_[i]->rsface_->center_.y); 1977 z_max = std::max(z_max, ses_->spheric_faces_[i]->rsface_->center_.z); 1978 } 1979 1980 double dist = 2*ses_->reduced_surface_->probe_radius_; 1981 1982 Position nx = (Position)((x_max-x_min)/dist+5); 1983 Position ny = (Position)((y_max-y_min)/dist+5); 1984 Position nz = (Position)((z_max-z_min)/dist+5); 1985 1986 Vector3 origin(x_min-2*dist,y_min-2*dist,z_min-2*dist); 1987 HashGrid3<Position> grid(origin,nx,ny,nz,dist); 1988 Vector3 pos; 1989 for (Position i = 0; i != ses_->number_of_spheric_faces_; i++) 1990 { 1991 pos.set(ses_->spheric_faces_[i]->rsface_->center_.x, 1992 ses_->spheric_faces_[i]->rsface_->center_.y, 1993 ses_->spheric_faces_[i]->rsface_->center_.z); 1994 grid.insert(pos,i); 1995 } 1996 1997 std::list<SESEdge*>::iterator edge; 1998 std::list<SESEdge*> deletable_edges; 1999 for (edge = ses_->singular_edges_.begin(); 2000 edge != ses_->singular_edges_.end(); 2001 edge++) 2002 { 2003 treatSingularEdge(*edge, grid, deletable_edges); 2004 } 2005 2006 for (edge = deletable_edges.begin(); edge != deletable_edges.end(); edge++) 2007 { 2008 (*edge)->face_[0]->edge_.remove(*edge); 2009 (*edge)->face_[1]->edge_.remove(*edge); 2010 (*edge)->vertex_[0]->edges_.erase(*edge); 2011 (*edge)->vertex_[1]->edges_.erase(*edge); 2012 ses_->edges_[(*edge)->index_] = NULL; 2013 ses_->singular_edges_.remove(*edge); 2014 delete *edge; 2015 } 2016 } 2017 treatSingularEdge(SESEdge * edge,HashGrid3<Position> & grid,std::list<SESEdge * > & deletable_edges)2018 void SESSingularityCleaner::treatSingularEdge(SESEdge* edge, HashGrid3<Position>& grid, 2019 std::list<SESEdge*>& deletable_edges) 2020 { 2021 if (edge->vertex_[0] == NULL) 2022 { 2023 return; 2024 } 2025 TAngle<double> phi = getOrientedAngle(edge->vertex_[0]->point_-edge->circle_.p, 2026 edge->vertex_[1]->point_-edge->circle_.p, 2027 edge->circle_.n); 2028 std::list<Intersection> intersections; 2029 getIntersectionsOfSingularEdge(edge,phi,grid,intersections); 2030 if (!intersections.empty()) 2031 { 2032 std::list<Intersection> min; 2033 std::list<Intersection> max; 2034 getExtrema(intersections,min,max); 2035 HashSet<Index> indices; 2036 std::list<Intersection>::const_iterator i; 2037 for (i = min.begin(); i != min.end(); i++) 2038 { 2039 indices.insert(i->first.second); 2040 } 2041 for (i = max.begin(); i != max.end(); i++) 2042 { 2043 indices.insert(i->first.second); 2044 } 2045 Index face1 = edge->face_[0]->index_; 2046 Index face2 = edge->face_[1]->index_; 2047 indices.insert(face1); 2048 indices.insert(face2); 2049 SESVertex* end_vertex1(0); 2050 SESVertex* end_vertex2(0); 2051 Index actual_min(0); 2052 Index actual_max(0); 2053 buildEndEdges(edge,min,max,end_vertex1,end_vertex2, 2054 actual_min,actual_max); 2055 Index next_face = actual_min; 2056 SESVertex* vertex = end_vertex1; 2057 while ((next_face != face2) && (vertex != NULL)) 2058 { 2059 buildEdge(edge,face1,next_face,face2,vertex,indices,true); 2060 } 2061 if (next_face != face2) 2062 { 2063 next_face = actual_max; 2064 vertex = end_vertex2; 2065 while ((next_face != face2) && (vertex != NULL)) 2066 { 2067 buildEdge(edge,face1,next_face,face2,vertex,indices,false); 2068 } 2069 } 2070 face2 = face1; 2071 face1 = edge->face_[1]->index_; 2072 next_face = actual_min; 2073 vertex = end_vertex1; 2074 while ((next_face != face2) && (vertex != NULL)) 2075 { 2076 buildEdge(edge,face1,next_face,face2,vertex,indices,true); 2077 } 2078 if (next_face != face2) 2079 { 2080 next_face = actual_max; 2081 vertex = end_vertex2; 2082 while ((next_face != face2) && (vertex != NULL)) 2083 { 2084 buildEdge(edge,face1,next_face,face2,vertex,indices,false); 2085 } 2086 } 2087 deletable_edges.push_back(edge); 2088 } 2089 } 2090 2091 getExtrema(const std::list<Intersection> & intersections,std::list<Intersection> & min,std::list<Intersection> & max)2092 void SESSingularityCleaner::getExtrema 2093 (const std::list<Intersection>& intersections, 2094 std::list<Intersection>& min, 2095 std::list<Intersection>& max) 2096 2097 { 2098 TAngle<double> min_phi(2*Constants::PI); 2099 TAngle<double> max_phi(0,true); 2100 std::list<Intersection>::const_iterator i; 2101 double epsilon = Constants::EPSILON; 2102 Constants::EPSILON = 0.0001; 2103 for (i = intersections.begin(); i != intersections.end(); i++) 2104 { 2105 if (i->first.first <= min_phi) 2106 { 2107 if (i->first.first < min_phi) 2108 { 2109 min.clear(); 2110 min_phi = i->first.first; 2111 } 2112 min.push_back(*i); 2113 } 2114 if (i->first.first >= max_phi) 2115 { 2116 if (i->first.first > max_phi) 2117 { 2118 max.clear(); 2119 max_phi = i->first.first; 2120 } 2121 max.push_back(*i); 2122 } 2123 } 2124 Constants::EPSILON = epsilon; 2125 } 2126 2127 getIntersectionsOfSingularEdge(SESEdge * edge,const TAngle<double> & phi,HashGrid3<Position> & grid,std::list<Intersection> & intersections)2128 void SESSingularityCleaner::getIntersectionsOfSingularEdge 2129 (SESEdge* edge, 2130 const TAngle<double>& phi, 2131 HashGrid3<Position>& grid, 2132 std::list<Intersection>& intersections) 2133 2134 { 2135 TQuaternion<double> rotate(edge->circle_.n,phi/2); 2136 TMatrix4x4<double> rotation; 2137 rotate.getRotationMatrix(rotation); 2138 TVector4<double> middle_(edge->vertex_[0]->point_.x-edge->circle_.p.x, 2139 edge->vertex_[0]->point_.y-edge->circle_.p.y, 2140 edge->vertex_[0]->point_.z-edge->circle_.p.z, 2141 0); 2142 middle_ = rotation*middle_; 2143 TVector3<double> middle(middle_.x+edge->circle_.p.x, 2144 middle_.y+edge->circle_.p.y, 2145 middle_.z+edge->circle_.p.z); 2146 Index face1 = edge->face_[0]->index_; 2147 Index face2 = edge->face_[1]->index_; 2148 TAngle<double> phi1; 2149 TAngle<double> phi2; 2150 TVector3<double> point1; 2151 TVector3<double> point2; 2152 TSphere3<double> probe; 2153 probe.radius = ses_->reduced_surface_->probe_radius_; 2154 Intersection intersection; 2155 HashGridBox3<Position>* box(0); 2156 HashGridBox3<Position>::ConstBoxIterator b; 2157 HashGridBox3<Position>::ConstDataIterator d; 2158 Vector3 pos(edge->circle_.p.x, 2159 edge->circle_.p.y, 2160 edge->circle_.p.z); 2161 box = grid.getBox(pos); 2162 if (box != NULL) 2163 { 2164 for (b = box->beginBox(); b != box->endBox(); b++) 2165 { 2166 for (d = b->beginData(); d != b->endData(); d++) 2167 { 2168 if (((Index)*d != face1) && ((Index)*d != face2)) 2169 { 2170 if (getIntersectionPointsAndAngles 2171 (edge->circle_,edge->vertex_[0]->point_, 2172 edge->face_[0]->index_,edge->face_[1]->index_, 2173 ses_->spheric_faces_[*d]->index_, 2174 phi1,point1,phi2,point2)) 2175 { 2176 probe.p = ses_->spheric_faces_[*d]->rsface_->center_; 2177 if (isIntersection(phi1,phi2,phi,middle,probe)) 2178 { 2179 intersection.first.second = *d; 2180 intersection.first.first = phi1; 2181 intersection.second = point1; 2182 intersections.push_back(intersection); 2183 intersection.first.first = phi2; 2184 intersection.second = point2; 2185 intersections.push_back(intersection); 2186 } 2187 } 2188 } 2189 } 2190 } 2191 } 2192 } 2193 2194 getIntersectionPointsAndAngles(const TCircle3<double> & circle,const TVector3<double> & point,Position index1,Position index2,Position probe_index,TAngle<double> & phi1,TVector3<double> & point1,TAngle<double> & phi2,TVector3<double> & point2)2195 bool SESSingularityCleaner::getIntersectionPointsAndAngles 2196 (const TCircle3<double>& circle, 2197 const TVector3<double>& point, 2198 Position index1, 2199 Position index2, 2200 Position probe_index, 2201 TAngle<double>& phi1, 2202 TVector3<double>& point1, 2203 TAngle<double>& phi2, 2204 TVector3<double>& point2) 2205 2206 { 2207 if (probeIntersection(index1, 2208 index2, 2209 probe_index, 2210 point1, 2211 point2)) 2212 { 2213 phi1 = getOrientedAngle(point-circle.p, 2214 point1-circle.p, 2215 circle.n); 2216 phi2 = getOrientedAngle(point-circle.p, 2217 point2-circle.p, 2218 circle.n); 2219 double epsilon = Constants::EPSILON; 2220 Constants::EPSILON = 0.001; 2221 if (Maths::isEqual(phi1.value,2*Constants::PI)) 2222 { 2223 phi1.value = 0.0; 2224 } 2225 if (Maths::isEqual(phi2.value,2*Constants::PI)) 2226 { 2227 phi2.value = 0.0; 2228 } 2229 Constants::EPSILON = epsilon; 2230 if (phi2 < phi1) 2231 { 2232 phi1.swap(phi2); 2233 point1.swap(point2); 2234 } 2235 return true; 2236 } 2237 else 2238 { 2239 return false; 2240 } 2241 } 2242 2243 isIntersection(const TAngle<double> & min_phi,const TAngle<double> & max_phi,const TAngle<double> & phi,const TVector3<double> & middle,const TSphere3<double> & probe)2244 bool SESSingularityCleaner::isIntersection 2245 (const TAngle<double>& min_phi, 2246 const TAngle<double>& max_phi, 2247 const TAngle<double>& phi, 2248 const TVector3<double>& middle, 2249 const TSphere3<double>& probe) 2250 2251 { 2252 double epsilon = Constants::EPSILON; 2253 Constants::EPSILON = 0.001; 2254 bool back; 2255 if (max_phi > phi) 2256 { 2257 back = false; 2258 } 2259 else 2260 { 2261 if (Maths::isNotZero(min_phi.value) || (max_phi < phi)) 2262 { 2263 back = true; 2264 } 2265 else 2266 { 2267 double epsilon = Constants::EPSILON; 2268 Constants::EPSILON = 1e-6; 2269 if (Maths::isGreater(probe.p.getSquareDistance(middle), 2270 probe.radius*probe.radius)) 2271 { 2272 back = false; 2273 } 2274 else 2275 { 2276 back = true; 2277 } 2278 Constants::EPSILON = epsilon; 2279 } 2280 } 2281 Constants::EPSILON = epsilon; 2282 return back; 2283 } 2284 2285 buildEndEdges(SESEdge * edge,const std::list<Intersection> & min,const std::list<Intersection> & max,SESVertex * & vertex1,SESVertex * & vertex2,Index & actual_min,Index & actual_max)2286 void SESSingularityCleaner::buildEndEdges 2287 (SESEdge* edge, 2288 const std::list<Intersection>& min, 2289 const std::list<Intersection>& max, 2290 SESVertex*& vertex1, 2291 SESVertex*& vertex2, 2292 Index& actual_min, 2293 Index& actual_max) 2294 2295 { 2296 buildEndEdge(edge,min,vertex1,actual_min,true); 2297 buildEndEdge(edge,max,vertex2,actual_max,false); 2298 } 2299 2300 buildEndEdge(SESEdge * edge,const std::list<Intersection> & extrema,SESVertex * & vertex,Index & actual_extremum,bool min)2301 void SESSingularityCleaner::buildEndEdge 2302 (SESEdge* edge, 2303 const std::list<Intersection>& extrema, 2304 SESVertex*& vertex, 2305 Index& actual_extremum, 2306 bool min) 2307 2308 { 2309 vertex = NULL; 2310 std::list<Intersection>::const_iterator m; 2311 for (m = extrema.begin(); m != extrema.end(); m++) 2312 { 2313 Index test = vertexExists(m->second); 2314 if (test != -1) 2315 { 2316 vertex = ses_->vertices_[test]; 2317 actual_extremum = m->first.second; 2318 } 2319 } 2320 if (vertex == NULL) 2321 { 2322 Intersection absolute_extremum = extrema.front(); 2323 if (min) 2324 { 2325 for (m = extrema.begin(); m != extrema.end(); m++) 2326 { 2327 if (m->first.first.value < absolute_extremum.first.first.value) 2328 { 2329 absolute_extremum = *m; 2330 } 2331 } 2332 } 2333 else 2334 { 2335 for (m = extrema.begin(); m != extrema.end(); m++) 2336 { 2337 if (m->first.first.value > absolute_extremum.first.first.value) 2338 { 2339 absolute_extremum = *m; 2340 } 2341 } 2342 } 2343 actual_extremum = absolute_extremum.first.second; 2344 vertex = new SESVertex(absolute_extremum.second, 2345 edge->circle_.p-absolute_extremum.second, 2346 -2,ses_->number_of_vertices_); 2347 ses_->vertices_.push_back(vertex); 2348 Vector3 pos(vertex->point_.x, 2349 vertex->point_.y, 2350 vertex->point_.z); 2351 vertex_grid_->insert(pos,vertex->index_); 2352 ses_->number_of_vertices_++; 2353 } 2354 Position v = (min ? 0 : 1); 2355 if (vertex != edge->vertex_[v]) 2356 { 2357 SESEdge* new_edge(0); 2358 new_edge = new SESEdge(*edge,true); 2359 new_edge->vertex_[1-v] = vertex; 2360 new_edge->rsedge_ = NULL; 2361 new_edge->index_ = ses_->number_of_edges_; 2362 ses_->edges_.push_back(new_edge); 2363 ses_->number_of_edges_++; 2364 ses_->singular_edges_.push_front(new_edge); 2365 ses_->number_of_singular_edges_++; 2366 new_edge->vertex_[0]->edges_.insert(new_edge); 2367 new_edge->vertex_[1]->edges_.insert(new_edge); 2368 new_edge->face_[0]->edge_.push_back(new_edge); 2369 new_edge->face_[1]->edge_.push_back(new_edge); 2370 new_edge->face_[0]->insert(new_edge->vertex_[1-v]); 2371 new_edge->face_[1]->insert(new_edge->vertex_[1-v]); 2372 vertex->faces_.insert(new_edge->face_[0]); 2373 vertex->faces_.insert(new_edge->face_[1]); 2374 } 2375 } 2376 2377 buildEdge(SESEdge * edge,Index face1,Index & face2,Index end,SESVertex * & vertex,const HashSet<Index> & indices,bool minimum)2378 void SESSingularityCleaner::buildEdge 2379 (SESEdge* edge, 2380 Index face1, 2381 Index& face2, 2382 Index end, 2383 SESVertex*& vertex, 2384 const HashSet<Index>& indices, 2385 bool minimum) 2386 2387 { 2388 SESFace* spheric_face1 = ses_->spheric_faces_[face1]; 2389 SESFace* spheric_face2 = ses_->spheric_faces_[face2]; 2390 TSphere3<double> probe1(spheric_face1->rsface_->center_, 2391 ses_->reduced_surface_->probe_radius_); 2392 TSphere3<double> probe2(spheric_face2->rsface_->center_, 2393 ses_->reduced_surface_->probe_radius_); 2394 TCircle3<double> circle; 2395 GetIntersection(probe1,probe2,circle); 2396 Index sign = (minimum ? -1 : 1); 2397 if (((probe1.p-edge->circle_.p)*edge->circle_.n)* 2398 ((probe1.p-circle.p)*circle.n)*sign > 0) 2399 { 2400 circle.n.negate(); 2401 } 2402 TVector3<double> point(vertex->point_); 2403 TVector3<double> point1; 2404 TVector3<double> point2; 2405 TAngle<double> phi1; 2406 TAngle<double> phi2; 2407 std::list< std::pair<TVector3<double>,Index> > min; 2408 TAngle<double> min_phi(2*Constants::PI,true); 2409 std::pair<TVector3<double>,Index> new_min; 2410 double epsilon = Constants::EPSILON; 2411 Constants::EPSILON = 0.001; 2412 HashSet<Index>::ConstIterator i; 2413 for (i = indices.begin(); i != indices.end(); i++) 2414 { 2415 if ((*i != face1) && (*i != face2)) 2416 { 2417 if (getIntersectionPointsAndAngles(circle,point,face1,face2,*i, 2418 phi1,point1,phi2,point2)) 2419 { 2420 if ((phi1 <= min_phi) && Maths::isGreater(phi1.value,0.0)) 2421 { 2422 if (phi1 < min_phi) 2423 { 2424 min.clear(); 2425 } 2426 min_phi = phi1; 2427 new_min.first = point1; 2428 new_min.second = *i; 2429 min.push_back(new_min); 2430 } 2431 if ((phi2 <= min_phi) && Maths::isGreater(phi2.value,0.0)) 2432 { 2433 if (phi2 < min_phi) 2434 { 2435 min.clear(); 2436 } 2437 min_phi = phi2; 2438 new_min.first = point2; 2439 new_min.second = *i; 2440 min.push_back(new_min); 2441 } 2442 } 2443 } 2444 } 2445 Constants::EPSILON = epsilon; 2446 SESVertex* new_vertex = NULL; 2447 std::list< std::pair<TVector3<double>,Index> >::iterator m 2448 = min.begin(); 2449 bool not_found = true; 2450 while (not_found && (m != min.end())) 2451 { 2452 if (m->second == end) 2453 { 2454 point2 = m->first; 2455 face2 = end; 2456 Index index = vertexExists(point2); 2457 if (index != -1) 2458 { 2459 new_vertex = ses_->vertices_[index]; 2460 } 2461 not_found = false; 2462 } 2463 m++; 2464 } 2465 if (not_found) 2466 { 2467 Index index = -1; 2468 m = min.begin(); 2469 while ((index == -1) && (m != min.end())) 2470 { 2471 point2 = m->first; 2472 face2 = m->second; 2473 index = vertexExists(point2); 2474 if (index != -1) 2475 { 2476 new_vertex = ses_->vertices_[index]; 2477 } 2478 m++; 2479 } 2480 } 2481 if (spheric_face1->isNeighbouredTo(spheric_face2) == false) 2482 { 2483 point1.set(vertex->point_); 2484 if (new_vertex == NULL) 2485 { 2486 new_vertex = new SESVertex(point2,circle.p-point2,-2, 2487 ses_->number_of_vertices_); 2488 ses_->vertices_.push_back(new_vertex); 2489 Vector3 pos(new_vertex->point_.x, 2490 new_vertex->point_.y, 2491 new_vertex->point_.z); 2492 vertex_grid_->insert(pos,new_vertex->index_); 2493 ses_->number_of_vertices_++; 2494 } 2495 SESEdge* new_edge = new SESEdge; 2496 new_edge->vertex_[0] = vertex; 2497 new_edge->vertex_[1] = new_vertex; 2498 new_edge->face_[0] = spheric_face1; 2499 new_edge->face_[1] = spheric_face2; 2500 new_edge->type_ = SESEdge::TYPE_SINGULAR; 2501 new_edge->circle_ = circle; 2502 new_edge->rsedge_ = NULL; 2503 new_edge->index_ = ses_->number_of_edges_; 2504 ses_->edges_.push_back(new_edge); 2505 ses_->number_of_edges_++; 2506 ses_->singular_edges_.push_back(new_edge); 2507 ses_->number_of_singular_edges_++; 2508 spheric_face1->edge_.push_back(new_edge); 2509 spheric_face2->edge_.push_back(new_edge); 2510 vertex->edges_.insert(new_edge); 2511 new_vertex->edges_.insert(new_edge); 2512 spheric_face1->insert(vertex); 2513 spheric_face2->insert(vertex); 2514 spheric_face1->insert(new_vertex); 2515 spheric_face2->insert(new_vertex); 2516 vertex->faces_.insert(spheric_face1); 2517 vertex->faces_.insert(spheric_face2); 2518 new_vertex->faces_.insert(spheric_face1); 2519 new_vertex->faces_.insert(spheric_face2); 2520 vertex = new_vertex; 2521 } 2522 else 2523 { 2524 vertex = NULL; 2525 } 2526 } 2527 2528 probeIntersection(Index face1,Index face2,Index face3,TVector3<double> & point1,TVector3<double> & point2)2529 bool SESSingularityCleaner::probeIntersection 2530 (Index face1, 2531 Index face2, 2532 Index face3, 2533 TVector3<double>& point1, 2534 TVector3<double>& point2) 2535 2536 { 2537 // sort the indices of the spheric faces 2538 sort(face1,face2,face3,face1,face2,face3); 2539 // try to find the intersection points 2540 HashMap< Position, 2541 HashMap< Position, 2542 HashMap< Position, 2543 ProbeIntersection* > > >::Iterator i1; 2544 HashMap< Position, 2545 HashMap< Position, 2546 ProbeIntersection* > >::Iterator i2; 2547 HashMap< Position,ProbeIntersection* >::Iterator i3; 2548 bool back = false; 2549 bool found = false; 2550 i1 = probe_intersections_.find(face1); 2551 if (i1 != probe_intersections_.end()) 2552 { 2553 i2 = i1->second.find(face2); 2554 if (i2 != i1->second.end()) 2555 { 2556 i3 = i2->second.find(face3); 2557 if (i3 != i2->second.end()) 2558 { 2559 found = true; 2560 if (i3->second == NULL) 2561 { 2562 back = false; 2563 } 2564 else 2565 { 2566 point1 = i3->second->point[0]; 2567 point2 = i3->second->point[1]; 2568 back = true; 2569 } 2570 } 2571 } 2572 } 2573 // if the intersection points are not computed yet, compute them now 2574 if (found == false) 2575 { 2576 TSphere3<double> s1(ses_->spheric_faces_[face1]->rsface_->center_, 2577 ses_->reduced_surface_->probe_radius_); 2578 TSphere3<double> s2(ses_->spheric_faces_[face2]->rsface_->center_, 2579 ses_->reduced_surface_->probe_radius_); 2580 TSphere3<double> s3(ses_->spheric_faces_[face3]->rsface_->center_, 2581 ses_->reduced_surface_->probe_radius_); 2582 if (GetIntersection(s1,s2,s3,point1,point2,false)) 2583 { 2584 ProbeIntersection* intersection = new ProbeIntersection; 2585 intersection->point[0] = point1; 2586 intersection->point[1] = point2; 2587 probe_intersections_[face1][face2][face3] = intersection; 2588 back = true; 2589 } 2590 else 2591 { 2592 probe_intersections_[face1][face2][face3] = NULL; 2593 back = false; 2594 } 2595 } 2596 return back; 2597 } 2598 2599 sort(Index u1,Index u2,Index u3,Index & s1,Index & s2,Index & s3)2600 void SESSingularityCleaner::sort 2601 (Index u1, 2602 Index u2, 2603 Index u3, 2604 Index& s1, 2605 Index& s2, 2606 Index& s3) 2607 2608 { 2609 s1 = u1; 2610 s2 = u2; 2611 s3 = u3; 2612 Index tmp; 2613 if (s1 > s2) 2614 { 2615 tmp = s1; 2616 s1 = s2; 2617 s2 = tmp; 2618 } 2619 if (s1 > s3) 2620 { 2621 tmp = s1; 2622 s1 = s3; 2623 s3 = tmp; 2624 } 2625 if (s2 > s3) 2626 { 2627 tmp = s2; 2628 s2 = s3; 2629 s3 = tmp; 2630 } 2631 } 2632 2633 vertexExists(TVector3<double> point)2634 Index SESSingularityCleaner::vertexExists(TVector3<double> point) 2635 2636 { 2637 double epsilon = Constants::EPSILON; 2638 Constants::EPSILON = 0.001; 2639 Vector3 p(point.x,point.y,point.z); 2640 HashGridBox3<Index>* box = vertex_grid_->getBox(p); 2641 HashGridBox3<Index>::ConstBoxIterator b; 2642 HashGridBox3<Index>::ConstDataIterator d; 2643 if (box != NULL) 2644 { 2645 for (b = box->beginBox(); b != box->endBox(); b++) 2646 { 2647 for (d = b->beginData(); d != b->endData(); d++) 2648 { 2649 if (ses_->vertices_[*d]->point_ == point) 2650 { 2651 Constants::EPSILON = epsilon; 2652 return *d; 2653 } 2654 } 2655 } 2656 } 2657 Constants::EPSILON = epsilon; 2658 return -1; 2659 } 2660 2661 2662 2663 } // namespace BALL 2664