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/reducedSurface.h> 9 10 #include <BALL/MATHS/analyticalGeometry.h> 11 #include <BALL/DATATYPE/hashGrid.h> 12 13 #include <iterator> 14 15 namespace BALL 16 { 17 ReducedSurface()18 ReducedSurface::ReducedSurface() 19 : number_of_atoms_(0), 20 atom_(), 21 probe_radius_(0.0), 22 number_of_vertices_(0), 23 vertices_(), 24 number_of_edges_(0), 25 edges_(), 26 number_of_faces_(0), 27 faces_(), 28 r_max_(0), 29 bounding_box_() 30 { 31 } 32 ReducedSurface(const ReducedSurface & reduced_surface,bool)33 ReducedSurface::ReducedSurface(const ReducedSurface& reduced_surface, bool) 34 : number_of_atoms_(0), 35 atom_(), 36 probe_radius_(0.0), 37 number_of_vertices_(0), 38 vertices_(), 39 number_of_edges_(0), 40 edges_(), 41 number_of_faces_(0), 42 faces_(), 43 r_max_(0), 44 bounding_box_() 45 { 46 copy(reduced_surface); 47 } 48 ReducedSurface(const std::vector<TSphere3<double>> & spheres,const double & probe_radius)49 ReducedSurface::ReducedSurface(const std::vector< TSphere3<double> >& spheres, 50 const double& probe_radius) 51 : number_of_atoms_(spheres.size()), 52 atom_(spheres), 53 probe_radius_(probe_radius), 54 number_of_vertices_(0), 55 vertices_(), 56 number_of_edges_(0), 57 edges_(), 58 number_of_faces_(0), 59 faces_(), 60 r_max_(0), 61 bounding_box_() 62 { 63 } 64 ~ReducedSurface()65 ReducedSurface::~ReducedSurface() 66 { 67 clear(); 68 } 69 set(const ReducedSurface & reduced_surface)70 void ReducedSurface::set(const ReducedSurface& reduced_surface) 71 { 72 copy(reduced_surface); 73 } 74 operator =(const ReducedSurface & reduced_surface)75 void ReducedSurface::operator = (const ReducedSurface& reduced_surface) 76 { 77 copy(reduced_surface); 78 } 79 clear()80 void ReducedSurface::clear() 81 { 82 for (Position i = 0; i < number_of_vertices_; i++) 83 { 84 delete vertices_[i]; 85 } 86 for (Position i = 0; i < number_of_edges_; i++) 87 { 88 delete edges_[i]; 89 } 90 for (Position i = 0; i < number_of_faces_; i++) 91 { 92 delete faces_[i]; 93 } 94 vertices_.clear(); 95 edges_.clear(); 96 faces_.clear(); 97 number_of_vertices_ = 0; 98 number_of_edges_ = 0; 99 number_of_faces_ = 0; 100 } 101 clean()102 void ReducedSurface::clean() 103 { 104 while ((number_of_vertices_ > 0) && 105 (vertices_[number_of_vertices_-1] == NULL)) 106 { 107 vertices_.pop_back(); 108 number_of_vertices_--; 109 } 110 111 for (Position i = 0; i < number_of_vertices_; i++) 112 { 113 if (vertices_[i] == NULL) 114 { 115 vertices_[i] = vertices_[number_of_vertices_-1]; 116 vertices_[i]->index_ = i; 117 vertices_.pop_back(); 118 number_of_vertices_--; 119 while (vertices_[number_of_vertices_-1] == NULL) 120 { 121 vertices_.pop_back(); 122 number_of_vertices_--; 123 } 124 } 125 } 126 while ((number_of_edges_ > 0) && 127 (edges_[number_of_edges_-1] == NULL)) 128 { 129 edges_.pop_back(); 130 number_of_edges_--; 131 } 132 for (Position i = 0; i < number_of_edges_; i++) 133 { 134 if (edges_[i] == NULL) 135 { 136 edges_[i] = edges_[number_of_edges_-1]; 137 edges_[i]->index_ = i; 138 edges_.pop_back(); 139 number_of_edges_--; 140 while (edges_[number_of_edges_-1] == NULL) 141 { 142 edges_.pop_back(); 143 number_of_edges_--; 144 } 145 } 146 } 147 while ((number_of_faces_ > 0) && 148 (faces_[number_of_faces_-1] == NULL)) 149 { 150 faces_.pop_back(); 151 number_of_faces_--; 152 } 153 for (Position i = 0; i < number_of_faces_; i++) 154 { 155 if (faces_[i] == NULL) 156 { 157 faces_[i] = faces_[number_of_faces_-1]; 158 faces_[i]->index_ = i; 159 faces_.pop_back(); 160 number_of_faces_--; 161 while (faces_[number_of_faces_-1] == NULL) 162 { 163 faces_.pop_back(); 164 number_of_faces_--; 165 } 166 } 167 } 168 } 169 170 numberOfAtoms() const171 Size ReducedSurface::numberOfAtoms() const 172 { 173 return number_of_atoms_; 174 } 175 numberOfVertices() const176 Size ReducedSurface::numberOfVertices() const 177 { 178 return number_of_vertices_; 179 } 180 numberOfEdges() const181 Size ReducedSurface::numberOfEdges() const 182 { 183 return number_of_edges_; 184 } 185 numberOfFaces() const186 Size ReducedSurface::numberOfFaces() const 187 { 188 return number_of_faces_; 189 } 190 getProbeRadius() const191 double ReducedSurface::getProbeRadius() const 192 { 193 return probe_radius_; 194 } 195 getSphere(Position i) const196 TSphere3<double> ReducedSurface::getSphere(Position i) const 197 throw(Exception::IndexOverflow) 198 { 199 if (i < number_of_atoms_) 200 { 201 return atom_[i]; 202 } 203 else 204 { 205 throw Exception::IndexOverflow(__FILE__, __LINE__, i, number_of_atoms_-1); 206 } 207 } 208 getVertex(Position i) const209 RSVertex* ReducedSurface::getVertex(Position i) const 210 throw(Exception::IndexOverflow) 211 { 212 if (i < number_of_vertices_) 213 { 214 return vertices_[i]; 215 } 216 else 217 { 218 throw Exception::IndexOverflow(__FILE__, __LINE__, i, number_of_vertices_-1); 219 } 220 } 221 getEdge(Position i) const222 RSEdge* ReducedSurface::getEdge(Position i) const 223 throw(Exception::IndexOverflow) 224 { 225 if (i < number_of_edges_) 226 { 227 return edges_[i]; 228 } 229 else 230 { 231 throw Exception::IndexOverflow(__FILE__, __LINE__, i, number_of_edges_-1); 232 } 233 } 234 getFace(Position i) const235 RSFace* ReducedSurface::getFace(Position i) const 236 throw(Exception::IndexOverflow) 237 { 238 if (i < number_of_faces_) 239 { 240 return faces_[i]; 241 } 242 else 243 { 244 throw Exception::IndexOverflow(__FILE__, __LINE__, i, number_of_faces_-1); 245 } 246 } 247 insert(RSVertex * rsvertex)248 void ReducedSurface::insert(RSVertex* rsvertex) 249 { 250 rsvertex->index_ = number_of_vertices_; 251 vertices_.push_back(rsvertex); 252 number_of_vertices_++; 253 } 254 insert(RSEdge * rsedge)255 void ReducedSurface::insert(RSEdge* rsedge) 256 { 257 rsedge->index_ = number_of_edges_; 258 edges_.push_back(rsedge); 259 number_of_edges_++; 260 } 261 insert(RSFace * rsface)262 void ReducedSurface::insert(RSFace* rsface) 263 { 264 rsface->index_ = number_of_faces_; 265 faces_.push_back(rsface); 266 number_of_faces_++; 267 } 268 getMaximalRadius() const269 double ReducedSurface::getMaximalRadius() const 270 { 271 return r_max_; 272 } 273 getBoundingBox() const274 TSimpleBox3<double> ReducedSurface::getBoundingBox() const 275 { 276 return bounding_box_; 277 } 278 deleteSimilarFaces(RSFace * face1,RSFace * face2)279 void ReducedSurface::deleteSimilarFaces(RSFace* face1, RSFace* face2) 280 { 281 if ((*face1) *= (*face2)) 282 { 283 // find the similar edges 284 std::vector<RSEdge*> rsedge1(3); 285 std::vector<RSEdge*> rsedge2(3); 286 287 findSimilarEdges(face1, face2, rsedge1, rsedge2); 288 289 // find the similar vertices 290 std::vector<RSVertex*> rsvertex1(3); 291 std::vector<RSVertex*> rsvertex2(3); 292 293 findSimilarVertices(face1, face2, rsvertex1, rsvertex2); 294 295 // join the similar vertices and delete the faces in their face lists 296 for (Position i = 0; i < 3; i++) 297 { 298 joinVertices(face1, face2, rsvertex1[i], rsvertex2[i]); 299 } 300 301 // correct the edges 302 for (Position i = 0; i < 3; i++) 303 { 304 correctEdges(face1, face2, rsedge1[i], rsedge2[i]); 305 } 306 307 faces_[face1->index_] = NULL; 308 faces_[face2->index_] = NULL; 309 310 delete face1; 311 delete face2; 312 } 313 } 314 findSimilarEdges(RSFace * face1,RSFace * face2,std::vector<RSEdge * > & rsedge1,std::vector<RSEdge * > & rsedge2)315 void ReducedSurface::findSimilarEdges(RSFace* face1, RSFace* face2, 316 std::vector<RSEdge*>& rsedge1, std::vector<RSEdge*>& rsedge2) 317 { 318 rsedge1[0] = face1->edge_[0]; 319 rsedge1[1] = face1->edge_[1]; 320 rsedge1[2] = face1->edge_[2]; 321 322 RSEdge* edge; 323 for (Position j = 0; j < 3; j++) 324 { 325 for (Position i = 0; i < 3; i++) 326 { 327 edge = face2->getEdge(i); 328 if (*edge *= *rsedge1[j]) 329 { 330 rsedge2[j] = edge; 331 } 332 } 333 } 334 } 335 findSimilarVertices(RSFace * face1,RSFace * face2,std::vector<RSVertex * > & rsvertex1,std::vector<RSVertex * > & rsvertex2)336 void ReducedSurface::findSimilarVertices(RSFace* face1, RSFace* face2, 337 std::vector<RSVertex*>& rsvertex1, 338 std::vector<RSVertex*>& rsvertex2) 339 { 340 rsvertex1[0] = face1->vertex_[0]; 341 rsvertex1[1] = face1->vertex_[1]; 342 rsvertex1[2] = face1->vertex_[2]; 343 344 RSVertex* vertex; 345 for (Position j = 0; j < 3; j++) 346 { 347 for (Position i = 0; i < 3; i++) 348 { 349 vertex = face2->getVertex(i); 350 if (vertex->atom_ == rsvertex1[j]->atom_) 351 { 352 rsvertex2[j] = vertex; 353 } 354 } 355 } 356 } 357 joinVertices(RSFace * face1,RSFace * face2,RSVertex * vertex1,RSVertex * vertex2)358 void ReducedSurface::joinVertices(RSFace* face1, RSFace* face2, 359 RSVertex* vertex1, RSVertex* vertex2) 360 { 361 if (vertex1 != vertex2) 362 { 363 vertex1->join(*vertex2); 364 vertex2->substitute(vertex1); 365 vertices_[vertex2->index_] = NULL; 366 367 delete vertex2; 368 } 369 370 vertex1->faces_.erase(face1); 371 vertex1->faces_.erase(face2); 372 } 373 correctEdges(RSFace * face1,RSFace * face2,RSEdge * edge1,RSEdge * edge2)374 void ReducedSurface::correctEdges(RSFace* face1, RSFace* face2, 375 RSEdge* edge1, RSEdge* edge2) 376 { 377 if (edge1 == edge2) 378 { 379 if (edge1->singular_) 380 { 381 edge1->vertex_[0]->edges_.erase(edge1); 382 edge1->vertex_[1]->edges_.erase(edge1); 383 edges_[edge1->index_] = NULL; 384 385 delete edge1; 386 } 387 else 388 { 389 edge1->face_[0] = NULL; 390 edge1->face_[1] = NULL; 391 edge1->angle_.value = 2*Constants::PI; 392 } 393 } 394 else 395 { 396 RSFace* neighbour2 = edge2->other(face2); 397 if (edge1->face_[0] == face1) 398 { 399 edge1->face_[0] = neighbour2; 400 } 401 else 402 { 403 edge1->face_[1] = neighbour2; 404 } 405 406 for (Position j = 0; j < 3; j++) 407 { 408 if (neighbour2->getEdge(j) == edge2) 409 { 410 neighbour2->setEdge(j,edge1); 411 } 412 } 413 414 edge2->vertex_[0]->edges_.erase(edge2); 415 edge2->vertex_[1]->edges_.erase(edge2); 416 edges_[edge2->index_] = NULL; 417 418 delete edge2; 419 420 // recomputation of rsedge1[i]->angle_ 421 RSFace* neighbour1 = edge1->face_[0]; 422 neighbour2 = edge1->face_[1]; 423 getAngle(neighbour1, neighbour2, 424 edge1->vertex_[0], edge1->vertex_[1], edge1->angle_, false); 425 } 426 } 427 getAngle(RSFace * face1,RSFace * face2,RSVertex * vertex1,RSVertex * vertex2,TAngle<double> & angle,bool check) const428 bool ReducedSurface::getAngle(RSFace* face1, RSFace* face2, 429 RSVertex* vertex1, RSVertex* vertex2, 430 TAngle<double>& angle, bool check) const 431 { 432 if (check) 433 { 434 if (! (face1->has(vertex1) && face1->has(vertex2) 435 && face2->has(vertex1) && face2->has(vertex2))) 436 { 437 return false; 438 } 439 } 440 441 RSVertex* vertex3(face1->third(vertex1, vertex2)); 442 443 TSphere3<double> atom1(atom_[vertex1->atom_]); 444 TSphere3<double> atom2(atom_[vertex2->atom_]); 445 TSphere3<double> atom3(atom_[vertex3->atom_]); 446 447 TVector3<double> axis(atom1.p - atom2.p); 448 TVector3<double> test(axis % face1->normal_); 449 450 if (Maths::isLess(test*(atom1.p - atom3.p), 0.0)) 451 { 452 axis.negate(); 453 } 454 455 atom1.radius += probe_radius_; 456 atom2.radius += probe_radius_; 457 458 TCircle3<double> circle; 459 GetIntersection(atom1, atom2, circle); 460 461 TVector3<double> v1 = face1->center_ - circle.p; 462 TVector3<double> v2 = face2->center_ - circle.p; 463 angle = getOrientedAngle(v1, v2, axis); 464 465 return true; 466 } 467 copy(const ReducedSurface & reduced_surface)468 void ReducedSurface::copy(const ReducedSurface& reduced_surface) 469 { 470 if (canBeCopied(reduced_surface)) 471 { 472 number_of_atoms_ = reduced_surface.number_of_atoms_; 473 atom_ = reduced_surface.atom_; 474 probe_radius_ = reduced_surface.probe_radius_; 475 number_of_vertices_ = reduced_surface.number_of_vertices_; 476 number_of_edges_ = reduced_surface.number_of_edges_; 477 number_of_faces_ = reduced_surface.number_of_faces_; 478 479 RSVertex* vertex; 480 for (Position i = 0; i < number_of_vertices_; i++) 481 { 482 vertex = new RSVertex(*reduced_surface.vertices_[i],false); 483 vertices_.push_back(vertex); 484 } 485 486 RSEdge* edge; 487 for (Position i = 0; i < number_of_edges_; i++) 488 { 489 edge = new RSEdge(*reduced_surface.edges_[i],false); 490 edges_.push_back(edge); 491 } 492 493 RSFace* face; 494 for (Position i = 0; i < number_of_faces_; i++) 495 { 496 face = new RSFace(*reduced_surface.faces_[i],false); 497 faces_.push_back(face); 498 } 499 500 HashSet<RSEdge*>::ConstIterator e; 501 HashSet<RSFace*>::ConstIterator f; 502 for (Position i = 0; i < number_of_vertices_; i++) 503 { 504 for (e = reduced_surface.vertices_[i]->edges_.begin(); 505 e != reduced_surface.vertices_[i]->edges_.end(); 506 e++) 507 { 508 vertices_[i]->edges_.insert(edges_[(*e)->index_]); 509 } 510 511 for (f = reduced_surface.vertices_[i]->faces_.begin(); 512 f != reduced_surface.vertices_[i]->faces_.end(); 513 f++) 514 { 515 vertices_[i]->faces_.insert(faces_[(*f)->index_]); 516 } 517 } 518 519 for (Position i = 0; i < number_of_edges_; i++) 520 { 521 edge = reduced_surface.edges_[i]; 522 edges_[i]->vertex_[0] = vertices_[edge->vertex_[0]->index_]; 523 edges_[i]->vertex_[1] = vertices_[edge->vertex_[1]->index_]; 524 edges_[i]->face_[0] = faces_[edge->face_[0]->index_]; 525 edges_[i]->face_[1] = faces_[edge->face_[1]->index_]; 526 } 527 528 for (Position i = 0; i < number_of_faces_; i++) 529 { 530 face = reduced_surface.faces_[i]; 531 faces_[i]->vertex_[0] = vertices_[face->vertex_[0]->index_]; 532 faces_[i]->vertex_[1] = vertices_[face->vertex_[1]->index_]; 533 faces_[i]->vertex_[2] = vertices_[face->vertex_[2]->index_]; 534 faces_[i]->edge_[0] = edges_[face->edge_[0]->index_]; 535 faces_[i]->edge_[1] = edges_[face->edge_[1]->index_]; 536 faces_[i]->edge_[2] = edges_[face->edge_[2]->index_]; 537 } 538 } 539 } 540 canBeCopied(const ReducedSurface & reduced_surface)541 bool ReducedSurface::canBeCopied(const ReducedSurface& reduced_surface) 542 { 543 for (Position i = 0; i < number_of_vertices_; i++) 544 { 545 if (reduced_surface.vertices_[i] == NULL) 546 { 547 return false; 548 } 549 550 if (reduced_surface.vertices_[i]->index_ < 0) 551 { 552 return false; 553 } 554 } 555 556 for (Position i = 0; i < number_of_edges_; i++) 557 { 558 if (reduced_surface.edges_[i] == NULL) 559 { 560 return false; 561 } 562 563 if (reduced_surface.edges_[i]->index_ < 0) 564 { 565 return false; 566 } 567 } 568 569 for (Position i = 0; i < number_of_faces_; i++) 570 { 571 if (reduced_surface.faces_[i] == NULL) 572 { 573 return false; 574 } 575 576 if (reduced_surface.faces_[i]->index_ < 0) 577 { 578 return false; 579 } 580 } 581 582 return true; 583 } 584 compute()585 void ReducedSurface::compute() 586 throw(Exception::GeneralException, 587 Exception::DivisionByZero, 588 Exception::IndexOverflow) 589 { 590 RSComputer rsc(this); 591 rsc.run(); 592 } 593 operator <<(std::ostream & s,const ReducedSurface & rs)594 std::ostream& operator << (std::ostream& s, const ReducedSurface& rs) 595 { 596 s << "Spheres:\n"; 597 for (Position i = 0; i < rs.numberOfAtoms(); i++) 598 { 599 s << " " << rs.getSphere(i) << "\n"; 600 } 601 s << "RSVertices:\n"; 602 for (Position i = 0; i < rs.numberOfVertices(); i++) 603 { 604 if (rs.getVertex(i) == NULL) 605 { 606 s << " --\n"; 607 } 608 else 609 { 610 s << " " << rs.getVertex(i) << " " << *(rs.getVertex(i)) << "\n"; 611 } 612 } 613 s << "RSEdges:\n"; 614 for (Position i = 0; i < rs.numberOfEdges(); i++) 615 { 616 if (rs.getEdge(i) == NULL) 617 { 618 s << " --\n"; 619 } 620 else 621 { 622 s << " " << rs.getEdge(i) << " " << *(rs.getEdge(i)) << "\n"; 623 } 624 } 625 s << "RSFaces:\n"; 626 for (Position i = 0; i < rs.numberOfFaces(); i++) 627 { 628 if (rs.getFace(i) == NULL) 629 { 630 s << " --\n"; 631 } 632 else 633 { 634 s << " " << rs.getFace(i) << " " << *(rs.getFace(i)) << "\n"; 635 } 636 } 637 return s; 638 } 639 640 //////////////////////////////////// 641 RSComputer()642 RSComputer::RSComputer() 643 : rs_(NULL), 644 neighbours_(), 645 atom_status_(), 646 neighbours_of_two_(), 647 probe_positions_(), 648 new_vertices_(), 649 new_faces_(), 650 vertices_() 651 { 652 } 653 654 RSComputer(ReducedSurface * rs)655 RSComputer::RSComputer(ReducedSurface* rs) 656 : rs_(rs), 657 neighbours_(rs->number_of_atoms_), 658 atom_status_(rs->number_of_atoms_,STATUS_UNKNOWN), 659 neighbours_of_two_(), 660 probe_positions_(), 661 new_vertices_(), 662 new_faces_(), 663 vertices_(rs->number_of_atoms_) 664 { 665 } 666 667 ~RSComputer()668 RSComputer::~RSComputer() 669 { 670 // delete probe_positions 671 HashMap<SortedPosition3, ProbePosition*>::Iterator pp1; 672 for (pp1 = probe_positions_.begin(); pp1 != probe_positions_.end(); ++pp1) 673 { 674 delete pp1->second; 675 } 676 } 677 run()678 void RSComputer::run() 679 throw(Exception::GeneralException, 680 Exception::DivisionByZero, 681 Exception::IndexOverflow) 682 { 683 // we need to use a larger value for epsilon for our computation 684 // so store the old value and restore it when we are done. 685 // TODO: this is bad in several ways; it is not thread safe, and 686 // I do not see why a smaller value shouldn't work also 687 double epsilon = Constants::EPSILON; 688 Constants::EPSILON = 1e-4; 689 690 // for each atom, find its neighbours 691 preProcessing(); 692 693 // start the computation 694 Position start = 1; 695 while (start != 0) 696 { 697 start = getStartPosition(); 698 699 switch (start) 700 { 701 case 2 : extendComponent(); 702 break; 703 case 3 : getRSComponent(); 704 break; 705 default : break; 706 } 707 } 708 709 rs_->clean(); 710 Constants::EPSILON = epsilon; 711 } 712 713 getRSComponent()714 void RSComputer::getRSComponent() 715 throw(Exception::GeneralException, 716 Exception::DivisionByZero, 717 Exception::IndexOverflow) 718 { 719 Position i = 0; 720 while (i < rs_->number_of_faces_) 721 { 722 if (rs_->faces_[i] != NULL) 723 { 724 if (treatFace(rs_->faces_[i]) == false) 725 { 726 i = 0; 727 } 728 else 729 { 730 i++; 731 } 732 } 733 else 734 { 735 i++; 736 } 737 } 738 extendComponent(); 739 //while (new_faces_.size() > 0) 740 //{ 741 // treatFace(*new_faces_.begin()); 742 //} 743 //extendComponent(); 744 } 745 treatFace(RSFace * face)746 bool RSComputer::treatFace(RSFace* face) 747 throw(Exception::GeneralException, 748 Exception::DivisionByZero, 749 Exception::IndexOverflow) 750 { 751 if (face->edge_[0]->face_[1] == NULL) 752 { 753 if (!treatEdge(face->edge_[0])) 754 { 755 return false; 756 } 757 } 758 if (face->edge_[1]->face_[1] == NULL) 759 { 760 if (!treatEdge(face->edge_[1])) 761 { 762 return false; 763 } 764 } 765 if (face->edge_[2]->face_[1] == NULL) 766 { 767 if (!treatEdge(face->edge_[2])) 768 { 769 return false; 770 } 771 } 772 new_faces_.erase(face); 773 return true; 774 } 775 treatEdge(RSEdge * edge)776 bool RSComputer::treatEdge(RSEdge* edge) 777 throw(Exception::GeneralException, 778 Exception::DivisionByZero, 779 Exception::IndexOverflow) 780 { 781 // This function rolls the probe sphere over a RSEdge. 782 // From all atoms that can be touced by the probe sphere when it touches 783 // the two atoms of the edge is this one selected for which the rotation 784 // angle is the smallest. A new face is found. 785 // If this face already exists the edge exists twice, too. These two 786 // edges and their vertices are joined. 787 // If the face does not exist yet, it will be created. A new vertex and 788 // two new edges will be created, too. 789 // In both cases the treated edge will be updated. It has not to be 790 // considerd again. 791 792 // find third atom 793 TAngle<double> phi; 794 TSphere3<double> probe; 795 RSFace* start_face(edge->face_[0]); // the edge already knows the 796 RSVertex* vertex1(edge->vertex_[0]); // starting face and their 797 RSVertex* vertex2(edge->vertex_[1]); // two vertices 798 RSVertex* vertex3(NULL); 799 Index atom1(vertex1->atom_); 800 Index atom2(vertex2->atom_); 801 Index atom3; 802 try 803 { 804 atom3 = thirdAtom(vertex1,vertex2,start_face,probe,phi); 805 } 806 catch (Exception::GeneralException& e) 807 { 808 String message = e.getMessage(); 809 String test_message = "PROBE SPHERE TOUCHES FOUR ATOMS"; 810 if (message == test_message) 811 { 812 return false; 813 } 814 else 815 { 816 throw; 817 } 818 } 819 TSphere3<double> sphere1(rs_->atom_[atom1]); 820 TSphere3<double> sphere2(rs_->atom_[atom2]); 821 TSphere3<double> sphere3(rs_->atom_[atom3]); 822 // build a new face and two new edges 823 vertex3 = new RSVertex(atom3); 824 RSEdge* edge1; 825 RSEdge* edge2; 826 RSFace* new_face // provisorial new face 827 = new RSFace(vertex1,vertex2,vertex3,NULL,NULL,NULL, 828 probe.p,getFaceNormal(sphere1,sphere2,sphere3,probe), 829 false,-1); 830 RSFace* test = faceExists(new_face,vertices_[vertex1->atom_]); 831 if (test == NULL) 832 { 833 // built face doesn't exist yet 834 // The new vertex has to be created since we don't know at this time 835 // whether it is a new vertex or not. 836 // Attention: one atom can build more than one vertex! 837 insert(vertex3); 838 edge1 = new RSEdge; 839 edge1->vertex_[0] = vertex2; 840 edge1->vertex_[1] = vertex3; 841 edge1->face_[0] = new_face; 842 edge2 = new RSEdge; 843 edge2->vertex_[0] = vertex3; 844 edge2->vertex_[1] = vertex1; 845 edge2->face_[0] = new_face; 846 new_face->edge_[0] = edge; 847 new_face->edge_[1] = edge1; 848 new_face->edge_[2] = edge2; 849 TPlane3<double> plane(sphere1.p,sphere2.p,sphere3.p); 850 new_face->singular_ = Maths::isLess(GetDistance(probe.p,plane), 851 rs_->probe_radius_); 852 insert(new_face); 853 } 854 else 855 { 856 // built face exitsts already 857 // the corresponding edge in the existing face has to be found 858 RSEdge* test_edge; 859 Index i = test->getSimilarEdge(edge,test_edge); 860 // Now the corresponding vertices of the corresponding edges have to be 861 // joined and one of them has to be deleted (if they are not equal). This 862 // is neccessary since creating a new face always creates a new vertex. 863 RSVertex* test_vertex1 = test_edge->vertex_[0]; 864 RSVertex* test_vertex2 = test_edge->vertex_[1]; 865 if (test_vertex1->atom_ == vertex2->atom_) 866 { 867 RSVertex* tmp = test_vertex1; 868 test_vertex1 = test_vertex2; 869 test_vertex2 = tmp; 870 } 871 // now we know which vertices are corresponding 872 if (*vertex1 != *test_vertex1) 873 { 874 // the vertices only have to be joined if they are not equal 875 vertex1->join(*test_vertex1); 876 test_vertex1->substitute(vertex1); 877 rs_->vertices_[test_vertex1->index_] = NULL; 878 new_vertices_.erase(test_vertex1); 879 vertices_[test_vertex1->atom_].remove(test_vertex1); 880 delete test_vertex1; 881 } 882 if (*vertex2 != *test_vertex2) 883 { 884 // the vertices only have to be joined if they are not equal 885 vertex2->join(*test_vertex2); 886 test_vertex2->substitute(vertex2); 887 rs_->vertices_[test_vertex2->index_] = NULL; 888 new_vertices_.erase(test_vertex2); 889 vertices_[test_vertex2->atom_].remove(test_vertex2); 890 delete test_vertex2; 891 } 892 // The vertices should have only one of the two corresponding edges. 893 // The other will be deleted later. 894 vertex1->edges_.erase(test_edge); 895 vertex2->edges_.erase(test_edge); 896 // The face should have only one of the two corresponding edges, too. 897 test->setEdge(i,edge); 898 // Now can we delete the build face and vertex and the double edge. 899 delete new_face; 900 if (test_edge->index_ != -1) // this can happens after a correct step 901 { 902 rs_->edges_[test_edge->index_] = NULL; 903 } 904 delete test_edge; 905 delete vertex3; 906 new_face = test; 907 } // face exitsts test 908 // update edge 909 TCircle3<double> circle1; 910 TCircle3<double> circle2; 911 TCircle3<double> circle3; 912 getCircles(atom1,atom2,circle1,circle2,circle3); 913 TVector3<double> ip1; // intersection points between 914 TVector3<double> ip2; // the edge and the probe sphere 915 TLine3<double> line(sphere1.p,sphere2.p,TLine3<double>::FORM__TWO_POINTS); 916 bool singular(GetIntersection(probe,line,ip1,ip2)); 917 if (singular && 918 Maths::isLess(ip1.getSquareDistance(sphere2.p), 919 ip2.getSquareDistance(sphere2.p))) 920 { // ip1 is the intersection point next to the first 921 ip1.swap(ip2); // vertex of the edge 922 } 923 edge->face_[1] = new_face; 924 edge->center_of_torus_ = circle1.p; 925 edge->radius_of_torus_ = circle1.radius; 926 edge->angle_ = phi; 927 edge->circle0_ = circle2; 928 edge->circle1_ = circle3; 929 edge->intersection_point0_ = ip1; 930 edge->intersection_point1_ = ip2; 931 edge->singular_ = singular; 932 if (edge->index_ == -1) 933 { 934 rs_->insert(edge); 935 } 936 return true; 937 } 938 correct(Index atom)939 void RSComputer::correct(Index atom) 940 { 941 std::list<RSVertex*>::iterator v; 942 RSVertex* vertex; 943 HashSet<RSFace*> faces; 944 HashSet<RSFace*> treat_faces; 945 HashSet<RSFace*>::Iterator f; 946 HashSet<RSVertex*> test_vertices; 947 HashSet<RSEdge*> delete_edges; 948 v = vertices_[atom].begin(); 949 while (v != vertices_[atom].end()) 950 { 951 treat_faces.clear(); 952 test_vertices.clear(); 953 delete_edges.clear(); 954 vertex = *v; 955 v++; 956 faces = vertex->faces_; 957 for (f = faces.begin(); f != faces.end(); f++) 958 { 959 (*f)->remove(delete_edges,test_vertices,treat_faces); 960 } 961 for (f = faces.begin(); f != faces.end(); f++) 962 { 963 treat_faces.erase(*f); 964 new_faces_.erase(*f); 965 rs_->faces_[(*f)->index_] = NULL; 966 delete *f; 967 } 968 for (f = treat_faces.begin(); f != treat_faces.end(); f++) 969 { 970 rs_->faces_[(*f)->index_] = NULL; 971 rs_->faces_.push_back(*f); 972 (*f)->index_ = rs_->number_of_faces_; 973 rs_->number_of_faces_++; 974 } 975 HashSet<RSEdge*>::Iterator edge; 976 for (edge = delete_edges.begin(); edge != delete_edges.end(); edge++) 977 { 978 Index index = (*edge)->index_; 979 if (index != -1) 980 { 981 rs_->edges_[index] = NULL; 982 } 983 delete *edge; 984 } 985 test_vertices.erase(vertex); 986 HashSet<RSVertex*>::Iterator test; 987 for (test = test_vertices.begin(); test != test_vertices.end(); test++) 988 { 989 if ((*test)->hasEdges() == false) 990 { 991 rs_->vertices_[(*test)->index_] = NULL; 992 vertices_[(*test)->atom_].remove(*test); 993 new_vertices_.erase(*test); 994 delete *test; 995 } 996 } 997 rs_->vertices_[vertex->index_] = NULL; 998 vertices_[atom].remove(vertex); 999 new_vertices_.erase(vertex); 1000 delete vertex; 1001 } 1002 rs_->atom_[atom].radius -= 10*Constants::EPSILON; 1003 atom_status_[atom] = STATUS_UNKNOWN; 1004 correctProbePosition(atom); 1005 } 1006 1007 extendComponent()1008 void RSComputer::extendComponent() 1009 throw(Exception::GeneralException, 1010 Exception::DivisionByZero, 1011 Exception::IndexOverflow) 1012 { 1013 std::deque<RSVertex*> new_vertices; 1014 std::copy(new_vertices_.begin(), new_vertices_.end(), std::back_inserter(new_vertices)); 1015 1016 while(!new_vertices.empty()) 1017 { 1018 RSFace* face = NULL; 1019 RSVertex* vertex1 = new_vertices.front(); 1020 new_vertices.pop_front(); 1021 Index atom1(vertex1->atom_); 1022 std::deque<Index>::const_iterator i = neighbours_[atom1].begin(); 1023 bool stop = false; 1024 while (!stop && i != neighbours_[atom1].end()) 1025 { 1026 if (atom_status_[*i] == STATUS_UNKNOWN) 1027 { 1028 Index atom2 = *i; 1029 const std::deque<Index>& s = neighboursOfTwoAtoms(SortedPosition2(atom1,atom2)); 1030 std::deque< std::pair< Index,TSphere3<double> > > candidates; 1031 findThirdAtom(atom1, atom2, s, candidates); 1032 if (candidates.empty()) 1033 { 1034 RSVertex* vertex2 = new RSVertex(atom2); 1035 RSEdge* edge = createFreeEdge(vertex1,vertex2); 1036 if (edge != NULL) 1037 { 1038 insert(edge); 1039 insert(vertex2); 1040 new_vertices.push_back(vertex1); 1041 new_vertices.push_back(vertex2); 1042 // i = neighbours_[atom1].end()--; ??? 1043 break; 1044 } 1045 else 1046 { 1047 delete vertex2; 1048 } 1049 } 1050 else 1051 { 1052 std::deque< std::pair< Index,TSphere3<double> > >::iterator j = candidates.begin(); 1053 while (j != candidates.end()) 1054 { 1055 if (atom_status_[j->first] == STATUS_UNKNOWN) 1056 { 1057 Index atom3 = j->first; 1058 TSphere3<double> probe = j->second; 1059 if (checkProbe(probe,SortedPosition3(atom1,atom2,atom3)) == true) 1060 { 1061 face = new RSFace; 1062 RSEdge* edge1 = new RSEdge; 1063 RSEdge* edge2 = new RSEdge; 1064 RSEdge* edge3 = new RSEdge; 1065 RSVertex* vertex2 = new RSVertex(atom2); 1066 RSVertex* vertex3 = new RSVertex(atom3); 1067 updateFaceAndEdges(vertex1,vertex2,vertex3, 1068 edge1,edge2,edge3, 1069 face,probe); 1070 insert(face); 1071 insert(vertex2); 1072 insert(vertex3); 1073 new_vertices.push_back(vertex1); 1074 new_vertices.push_back(vertex2); 1075 new_vertices.push_back(vertex3); 1076 // i = neighbours_[atom1].end()--; 1077 // j = candidates.end()--; 1078 // ???? 1079 stop = true; 1080 break; 1081 } 1082 } 1083 ++j; 1084 } // while j 1085 } 1086 } 1087 ++i; 1088 } // while i 1089 if (face != NULL) 1090 { 1091 getRSComponent(); 1092 } 1093 } 1094 1095 new_vertices_.clear(); 1096 } 1097 thirdAtom(RSVertex * vertex1,RSVertex * vertex2,RSFace * face,TSphere3<double> & probe,TAngle<double> & phi)1098 Index RSComputer::thirdAtom(RSVertex* vertex1, RSVertex* vertex2, 1099 RSFace* face, TSphere3<double>& probe, TAngle<double>& phi) 1100 throw(Exception::GeneralException, 1101 Exception::DivisionByZero, 1102 Exception::IndexOverflow) 1103 { 1104 // This function chooses from all atoms which can be touced by the probe 1105 // sphere when it touches the given two vertices this one, for which is 1106 // the rotation angle the smalest. 1107 // If the rotation angle equals zero, the probe sphere can touch four or 1108 // more atoms an an exception is thrown. 1109 // If no atom can be found an exception is thrown. 1110 Index atom1(vertex1->atom_); 1111 Index atom2(vertex2->atom_); 1112 const std::deque<Index>& s = neighboursOfTwoAtoms(SortedPosition2(atom1,atom2)); 1113 std::deque<std::pair<Index,TSphere3<double> > > candidates; 1114 findThirdAtom(atom1, atom2, s, candidates); 1115 std::deque<std::pair<Index,TSphere3<double> > >::iterator k; 1116 TAngle<double> old_angle(3*Constants::PI,true); 1117 TAngle<double> new_angle; 1118 TAngle<double> two_pi(2*Constants::PI,true); 1119 TVector3<double> axis = rs_->atom_[atom1].p-rs_->atom_[atom2].p; 1120 TVector3<double> test_vector = face->normal_%axis; 1121 Index third_face_atom = face->third(vertex1,vertex2)->atom_; 1122 if (Maths::isLess(test_vector*rs_->atom_[third_face_atom].p, 1123 test_vector*rs_->atom_[atom1].p) ) 1124 { 1125 axis.negate(); 1126 } 1127 TSphere3<double> sphere1(rs_->atom_[atom1]); 1128 TSphere3<double> sphere2(rs_->atom_[atom2]); 1129 sphere1.radius += rs_->probe_radius_; 1130 sphere2.radius += rs_->probe_radius_; 1131 TCircle3<double> circle; 1132 GetIntersection(sphere1,sphere2,circle); 1133 TVector3<double> start_probe = face->center_; 1134 TVector3<double> v1 = start_probe-circle.p; 1135 TVector3<double> face_normal = face->normal_; 1136 std::deque<std::pair<Index,TSphere3<double> > > third; 1137 for (k = candidates.begin(); k != candidates.end(); k++) 1138 { 1139 if ((k->first != third_face_atom) || (k->second.p != start_probe)) 1140 // not found the starting face 1141 { 1142 TVector3<double> v2(k->second.p-circle.p); 1143 new_angle = getOrientedAngle(v1,v2,axis); 1144 if (Maths::isZero(new_angle.value) || (new_angle == two_pi)) 1145 { 1146 correct(k->first); 1147 throw Exception::GeneralException 1148 (__FILE__,__LINE__,"CAN'T COMPUTE RS", 1149 "PROBE SPHERE TOUCHES FOUR ATOMS"); 1150 } 1151 if (new_angle <= old_angle) 1152 { 1153 if (new_angle < old_angle) 1154 { 1155 old_angle = new_angle; 1156 std::deque<std::pair<Index,TSphere3<double> > >::iterator t; 1157 for (t = third.begin(); t != third.end(); t++) 1158 { 1159 if (atom_status_[t->first] == STATUS_UNKNOWN) 1160 { 1161 atom_status_[t->first] = STATUS_INSIDE; 1162 } 1163 } 1164 third.clear(); 1165 } 1166 third.push_back(*k); 1167 } 1168 else 1169 { 1170 if (atom_status_[k->first] == STATUS_UNKNOWN) 1171 { 1172 atom_status_[k->first] = STATUS_INSIDE; 1173 } 1174 } 1175 } 1176 } 1177 if (third.size() > 1) 1178 { 1179 k = third.begin(); 1180 k++; 1181 while (k != third.end()) 1182 { 1183 correct(k->first); 1184 k++; 1185 } 1186 throw Exception::GeneralException 1187 (__FILE__,__LINE__,"CAN'T COMPUTE RS", 1188 "PROBE SPHERE TOUCHES FOUR ATOMS"); 1189 } 1190 probe = third.front().second; 1191 phi.set(old_angle.value,true); 1192 atom_status_[third.front().first] = STATUS_ON_SURFACE; 1193 return third.front().first; 1194 } 1195 getStartPosition()1196 Position RSComputer::getStartPosition() 1197 throw(Exception::DivisionByZero) 1198 { 1199 if (findFirstFace() != NULL) 1200 { 1201 return 3; 1202 } 1203 if (findFirstEdge() != NULL) 1204 { 1205 return 2; 1206 } 1207 if (findFirstVertex() != NULL) 1208 { 1209 return 1; 1210 } 1211 return 0; 1212 } 1213 findFirstFace()1214 RSFace* RSComputer::findFirstFace() 1215 throw(Exception::DivisionByZero) 1216 { 1217 for (Position direction = 0; direction < 3; direction++) 1218 { 1219 for (Position extreme = 0; extreme < 1; extreme++) 1220 { 1221 RSFace* face = findFace(direction, extreme); 1222 1223 if (face != NULL) 1224 { 1225 return face; 1226 } 1227 } 1228 } 1229 return NULL; 1230 } 1231 findFirstEdge()1232 RSEdge* RSComputer::findFirstEdge() 1233 { 1234 for (Position direction = 0; direction < 3; direction++) 1235 { 1236 for (Position extrem = 0; extrem < 1; extrem++) 1237 { 1238 RSEdge* edge = findEdge(direction,extrem); 1239 1240 if (edge != NULL) 1241 { 1242 return edge; 1243 } 1244 } 1245 } 1246 return NULL; 1247 } 1248 findFirstVertex()1249 RSVertex* RSComputer::findFirstVertex() 1250 { 1251 for (Position i = 0; i < rs_->number_of_atoms_; i++) 1252 { 1253 if (atom_status_[i] == STATUS_UNKNOWN) 1254 { 1255 if (neighbours_[i].size() == 0) 1256 { 1257 RSVertex* vertex = new RSVertex(i); 1258 insert(vertex); 1259 1260 return vertex; 1261 } 1262 } 1263 } 1264 return NULL; 1265 } 1266 findFace(Position direction,Position extrem)1267 RSFace* RSComputer::findFace(Position direction, Position extrem) 1268 throw(Exception::DivisionByZero) 1269 { 1270 Index a1 = findFirstAtom(direction, extrem); 1271 if (a1 == -1) 1272 { 1273 return NULL; 1274 } 1275 1276 Index a2 = findSecondAtom(a1, direction, extrem); 1277 if (a2 == -1) 1278 { 1279 return NULL; 1280 } 1281 1282 const std::deque<Index>& s = neighboursOfTwoAtoms(SortedPosition2(a1, a2)); 1283 std::deque<std::pair<Index,TSphere3<double> > > candidates; 1284 findThirdAtom(a1, a2, s, candidates); 1285 if (candidates.empty()) 1286 { 1287 return NULL; 1288 } 1289 1290 std::deque<std::pair<Index,TSphere3<double> > >::iterator i = candidates.begin(); 1291 Index a3 = -1; 1292 TSphere3<double> probe; 1293 bool found = false; 1294 while ((found == false) && (i != candidates.end())) 1295 { 1296 a3 = i->first; 1297 probe = i->second; 1298 found = (atom_status_[a3] == STATUS_UNKNOWN) && 1299 checkProbe(probe,SortedPosition3(a1,a2,a3)); 1300 i++; 1301 } 1302 if (found) 1303 { 1304 RSVertex* vertex1 = new RSVertex(a1); 1305 RSVertex* vertex2 = new RSVertex(a2); 1306 RSVertex* vertex3 = new RSVertex(a3); 1307 1308 RSEdge* e1 = new RSEdge; 1309 RSEdge* e2 = new RSEdge; 1310 RSEdge* e3 = new RSEdge; 1311 1312 RSFace* face = new RSFace; 1313 1314 updateFaceAndEdges(vertex1,vertex2,vertex3,e1,e2,e3,face,probe); 1315 1316 insert(face); 1317 insert(vertex1); 1318 insert(vertex2); 1319 insert(vertex3); 1320 1321 return face; 1322 } 1323 else 1324 { 1325 atom_status_[a1] = STATUS_INSIDE; 1326 atom_status_[a2] = STATUS_INSIDE; 1327 return NULL; 1328 } 1329 } 1330 findEdge(Position direction,Position extrem)1331 RSEdge* RSComputer::findEdge(Position direction, Position extrem) 1332 { 1333 Index a1 = findFirstAtom(direction,extrem); 1334 if (a1 == -1) 1335 { 1336 return NULL; 1337 } 1338 1339 Index a2 = findSecondAtom(a1,direction,extrem); 1340 if (a2 == -1) 1341 { 1342 return NULL; 1343 } 1344 1345 RSVertex* vertex1 = new RSVertex(a1); 1346 RSVertex* vertex2 = new RSVertex(a2); 1347 neighboursOfTwoAtoms(SortedPosition2(a1,a2)); 1348 1349 RSEdge* edge = createFreeEdge(vertex1,vertex2); 1350 if (edge != NULL) 1351 { 1352 insert(edge); 1353 insert(vertex1); 1354 insert(vertex2); 1355 1356 return edge; 1357 } 1358 else 1359 { 1360 delete vertex1; 1361 delete vertex2; 1362 neighbours_[a1].erase(std::remove(neighbours_[a1].begin(), neighbours_[a1].end(), a2), neighbours_[a1].end()); 1363 neighbours_[a2].erase(std::remove(neighbours_[a2].begin(), neighbours_[a2].end(), a1), neighbours_[a2].end()); 1364 1365 return NULL; 1366 } 1367 } 1368 findFirstAtom(Position direction,Position extrem)1369 Index RSComputer::findFirstAtom(Position direction, Position extrem) 1370 { 1371 Index extrem_atom = -1; 1372 // find the first atom of unknown status 1373 Index i = 0; 1374 bool found = false; 1375 while ((found == false) && (i < (Index)rs_->number_of_atoms_)) 1376 { 1377 if (atom_status_[i] == STATUS_UNKNOWN) 1378 { 1379 found = true; 1380 } 1381 else 1382 { 1383 i++; 1384 } 1385 } 1386 if (found) 1387 { 1388 extrem_atom = i; 1389 TSphere3<double>* next_atom = &(rs_->atom_[i]); 1390 double extrem_value 1391 = ((extrem == 0) ? next_atom->p[direction]-next_atom->radius 1392 : next_atom->p[direction]+next_atom->radius); 1393 i++; 1394 // find the atom of unknown status lying on the extrem position 1395 while (i < (Index)rs_->number_of_atoms_) 1396 { 1397 if (atom_status_[i] == STATUS_UNKNOWN) 1398 { 1399 next_atom = &(rs_->atom_[i]); 1400 double extremum 1401 = ((extrem == 0) ? next_atom->p[direction]-next_atom->radius 1402 : next_atom->p[direction]+next_atom->radius); 1403 if (((extrem == 0) && Maths::isLess(extremum, extrem_value)) || 1404 ((extrem != 0) && Maths::isGreater(extremum, extrem_value))) 1405 { 1406 extrem_value = extremum; 1407 extrem_atom = i; 1408 } 1409 } 1410 i++; 1411 } 1412 } 1413 return extrem_atom; 1414 } 1415 findSecondAtom(Index atom,Position direction,Position extrem)1416 Index RSComputer::findSecondAtom(Index atom, Position direction, Position extrem) 1417 { 1418 Index second_atom = -1; 1419 // find the first neighbour atom of unknown status 1420 std::deque<Index>::const_iterator i = neighbours_[atom].begin(); 1421 bool found = false; 1422 while ((found == false) && (i != neighbours_[atom].end())) 1423 { 1424 if (atom_status_[*i] == STATUS_UNKNOWN) 1425 { 1426 found = true; 1427 } 1428 else 1429 { 1430 i++; 1431 } 1432 } 1433 if (found) 1434 { 1435 second_atom = *i; 1436 TSphere3<double> first_atom(rs_->atom_[atom]); 1437 first_atom.radius += rs_->probe_radius_; 1438 double extrem_value 1439 = ((extrem == 0) ? first_atom.p[direction]+first_atom.radius 1440 : first_atom.p[direction]-first_atom.radius); 1441 TSphere3<double> next_atom; 1442 TCircle3<double> intersection_circle; 1443 // find the neighbour atom of unknown status lying on the extrem position 1444 while (i != neighbours_[atom].end()) 1445 { 1446 if (atom_status_[*i] == STATUS_UNKNOWN) 1447 { 1448 next_atom = rs_->atom_[*i]; 1449 next_atom.radius += rs_->probe_radius_; 1450 if (GetIntersection(first_atom,next_atom,intersection_circle)) 1451 { 1452 double next_extrem 1453 = getCircleExtremum(intersection_circle,direction,extrem); 1454 if (((extrem == 0) && Maths::isLess(next_extrem,extrem_value)) || 1455 ((extrem != 0) && Maths::isGreater(next_extrem,extrem_value))) 1456 { 1457 extrem_value = next_extrem; 1458 second_atom = *i; 1459 } 1460 } 1461 } 1462 i++; 1463 } 1464 } 1465 return second_atom; 1466 } 1467 findThirdAtom(Index atom1,Index atom2,const std::deque<Index> & third,std::deque<std::pair<Index,TSphere3<double>>> & atoms)1468 void RSComputer::findThirdAtom(Index atom1, Index atom2, const std::deque<Index>& third, 1469 std::deque<std::pair<Index,TSphere3<double> > >& atoms) 1470 { 1471 // This function computes a list of all atoms (with its probe positions) 1472 // which can be touched by the probe sphere when it touches the two given 1473 // atoms 1474 std::pair<Index, TSphere3<double> > candidate; 1475 std::deque<Index>::const_iterator i = third.begin(); 1476 TVector3<double> center1, center2; 1477 TSphere3<double> probe; 1478 probe.radius = rs_->probe_radius_; 1479 while (i != third.end()) 1480 { 1481 if (centerOfProbe(SortedPosition3(atom1,atom2,*i),center1,center2)) 1482 { 1483 if (!(Maths::isNan(center1.x) || Maths::isNan(center1.y) || Maths::isNan(center1.z))) 1484 { 1485 probe.p.set(center1); 1486 candidate.first = *i; 1487 candidate.second = probe; 1488 atoms.push_back(candidate); 1489 } 1490 1491 if (!(Maths::isNan(center2.x) || Maths::isNan(center2.y) || Maths::isNan(center2.z))) 1492 { 1493 probe.p.set(center2); 1494 candidate.first = *i; 1495 candidate.second = probe; 1496 atoms.push_back(candidate); 1497 } 1498 } 1499 i++; 1500 } 1501 } 1502 neighboursOfTwoAtoms(const SortedPosition2 & pos)1503 const std::deque<Index>& RSComputer::neighboursOfTwoAtoms(const SortedPosition2& pos) 1504 { 1505 HashMap<SortedPosition2, std::deque<Index> >::Iterator n1 = neighbours_of_two_.find(pos); 1506 1507 if (n1 == neighbours_of_two_.end()) 1508 { 1509 n1 = neighbours_of_two_.insert(std::make_pair(pos, std::deque<Index>())).first; 1510 1511 std::set_intersection(neighbours_[pos.a].begin(), neighbours_[pos.a].end(), 1512 neighbours_[pos.b].begin(), neighbours_[pos.b].end(), 1513 std::back_inserter(n1->second)); 1514 } 1515 1516 return n1->second; 1517 } 1518 neighboursOfThreeAtoms(Index atom1,Index atom2,Index atom3,std::deque<Index> & output_list)1519 void RSComputer::neighboursOfThreeAtoms(Index atom1, Index atom2, Index atom3, 1520 std::deque<Index>& output_list) 1521 { 1522 SortedPosition2 pos12(atom1, atom2); 1523 SortedPosition2 pos13(atom1, atom3); 1524 1525 const std::deque<Index>& s1 = neighboursOfTwoAtoms(pos12); 1526 const std::deque<Index>& s2 = neighboursOfTwoAtoms(pos13); 1527 1528 std::set_intersection(s1.begin(), s1.end(), s2.begin(), s2.end(), std::back_inserter(output_list)); 1529 } 1530 getCircleExtremum(const TCircle3<double> & circle,Position direction,Position extrem)1531 double RSComputer::getCircleExtremum(const TCircle3<double>& circle, 1532 Position direction, Position extrem) 1533 { 1534 double min = 0; 1535 double max = 0; 1536 1537 TVector3<double> norm2(circle.n.x * circle.n.x, 1538 circle.n.y * circle.n.y, 1539 circle.n.z * circle.n.z); 1540 1541 switch (direction) 1542 { 1543 case 0 : 1544 if (Maths::isZero(circle.n.y) && Maths::isZero(circle.n.z)) 1545 { 1546 min = max = circle.p.x; 1547 } 1548 else 1549 { 1550 double x_norm = norm2.y + norm2.z; 1551 x_norm /= norm2.x+x_norm; 1552 x_norm = circle.radius * sqrt(x_norm); 1553 min = (circle.p.x) - x_norm; 1554 max = (circle.p.x) + x_norm; 1555 } 1556 break; 1557 case 1 : 1558 if (Maths::isZero(circle.n.x) && Maths::isZero(circle.n.z)) 1559 { 1560 min = max = circle.p.y; 1561 } 1562 else 1563 { 1564 double y_norm = norm2.x + norm2.z; 1565 y_norm /= norm2.y + y_norm; 1566 y_norm = circle.radius * sqrt(y_norm); 1567 min = (circle.p.y)-y_norm; 1568 max = (circle.p.y)+y_norm; 1569 } 1570 break; 1571 case 2 : 1572 if (Maths::isZero(circle.n.x) && Maths::isZero(circle.n.y)) 1573 { 1574 min = max = circle.p.z; 1575 } 1576 else 1577 { 1578 double z_norm = norm2.x + norm2.y; 1579 z_norm /= norm2.z + z_norm; 1580 z_norm = circle.radius * sqrt(z_norm); 1581 min = circle.p.z - z_norm; 1582 max = circle.p.z + z_norm; 1583 } 1584 break; 1585 } 1586 if (extrem == 0) 1587 { 1588 return min; 1589 } 1590 else 1591 { 1592 return max; 1593 } 1594 } 1595 createFreeEdge(RSVertex * vertex1,RSVertex * vertex2)1596 RSEdge* RSComputer::createFreeEdge(RSVertex* vertex1, RSVertex* vertex2) 1597 { 1598 Index atom1 = vertex1->atom_; 1599 Index atom2 = vertex2->atom_; 1600 1601 TCircle3<double> circle1, circle2, circle3; 1602 1603 // compute the three circles describing the toric face 1604 if (getCircles(atom1, atom2, circle1, circle2, circle3) && // the probe hulls intersect 1605 Maths::isGreater(circle1.radius,rs_->probe_radius_)) // the radius of the toric edge is > 0 1606 { 1607 TPlane3<double> plane(circle1.p, circle1.n); 1608 1609 std::deque<Index>::const_iterator i; 1610 TCircle3<double> test_circle; 1611 TSphere3<double> sphere; 1612 1613 const std::deque<Index>& s = neighboursOfTwoAtoms(SortedPosition2(atom1, atom2)); 1614 1615 // find the mutual neighbours of both atoms 1616 for (i = s.begin(); i != s.end(); i++) 1617 { 1618 // put a sphere into the neighbour 1619 sphere.set(rs_->atom_[*i].p, rs_->atom_[*i].radius+rs_->probe_radius_); 1620 1621 if (GetIntersection(sphere,plane,test_circle)) 1622 { 1623 double radius_dist = test_circle.radius-circle1.radius; 1624 double radius_sum = test_circle.radius+circle1.radius; 1625 double center_dist = test_circle.p.getSquareDistance(circle1.p); 1626 1627 if ( Maths::isLessOrEqual(radius_dist*radius_dist, center_dist) 1628 && Maths::isGreaterOrEqual(radius_sum*radius_sum, center_dist) ) // the circles intersect 1629 { 1630 return NULL; 1631 } 1632 } 1633 } 1634 TVector3<double> vector(0,0,0); 1635 RSEdge* edge = new RSEdge(vertex1, vertex2, NULL, NULL, circle1.p, circle1.radius, 1636 TAngle<double>(2*Constants::PI, true), circle2, circle3, 1637 vector, vector, false, -1); 1638 1639 return edge; 1640 } 1641 1642 return NULL; 1643 } 1644 1645 // circle1 will be the circle of the edge of the probe sphere along a toric edge 1646 // circle2 will be the rim of the cut sphere around atom1 1647 // circle3 will be the rim of the cut sphere around atom2 getCircles(Index atom1,Index atom2,TCircle3<double> & circle1,TCircle3<double> & circle2,TCircle3<double> & circle3)1648 bool RSComputer::getCircles(Index atom1, Index atom2, TCircle3<double>& circle1, 1649 TCircle3<double>& circle2, TCircle3<double>& circle3) 1650 { 1651 TSphere3<double> sphere1(rs_->atom_[atom1]); 1652 TSphere3<double> sphere2(rs_->atom_[atom2]); 1653 1654 sphere1.radius += rs_->probe_radius_; 1655 sphere2.radius += rs_->probe_radius_; 1656 1657 // intersect the spheres surrounding the atoms to yield circle1 1658 if (GetIntersection(sphere1, sphere2, circle1)) 1659 { 1660 // the intercept theorem yields a simple relationship between the radii 1661 double ratio = rs_->atom_[atom1].radius/sphere1.radius; 1662 1663 circle2.radius = circle1.radius*ratio; 1664 circle2.p = sphere1.p+(circle1.p-sphere1.p)*ratio; 1665 1666 ratio = rs_->atom_[atom2].radius/sphere2.radius; 1667 1668 circle3.radius = circle1.radius*ratio; 1669 circle3.p = sphere2.p+(circle1.p-sphere2.p)*ratio; 1670 1671 return true; 1672 } 1673 1674 return false; 1675 } 1676 getFaceNormal(const TSphere3<double> & atom1,const TSphere3<double> & atom2,const TSphere3<double> & atom3,const TSphere3<double> & probe)1677 TVector3<double> RSComputer::getFaceNormal(const TSphere3<double>& atom1, 1678 const TSphere3<double>& atom2, 1679 const TSphere3<double>& atom3, 1680 const TSphere3<double>& probe) 1681 { 1682 TPlane3<double> plane(atom1.p, atom2.p, atom3.p); 1683 TVector3<double> norm(plane.n); 1684 if (Maths::isLess(norm*probe.p, norm*atom1.p)) 1685 { 1686 norm.negate(); 1687 } 1688 1689 return norm; 1690 } 1691 updateFaceAndEdges(RSVertex * v1,RSVertex * v2,RSVertex * v3,RSEdge * e1,RSEdge * e2,RSEdge * e3,RSFace * f,const TSphere3<double> & probe)1692 void RSComputer::updateFaceAndEdges(RSVertex* v1, RSVertex* v2, RSVertex* v3, 1693 RSEdge* e1, RSEdge* e2, RSEdge* e3, 1694 RSFace* f, const TSphere3<double>& probe) 1695 { 1696 e1->vertex_[0] = v1; e1->vertex_[1] = v2; e1->face_[0] = f; 1697 e2->vertex_[0] = v2; e2->vertex_[1] = v3; e2->face_[0] = f; 1698 e3->vertex_[0] = v3; e3->vertex_[1] = v1; e3->face_[0] = f; 1699 1700 f->vertex_[0] = v1; f->vertex_[1] = v2; f->vertex_[2] = v3; 1701 f->edge_[0] = e1; f->edge_[1] = e2; f->edge_[2] = e3; 1702 1703 f->center_ = probe.p; 1704 1705 TPlane3<double> plane(rs_->atom_[v1->atom_].p, 1706 rs_->atom_[v2->atom_].p, 1707 rs_->atom_[v3->atom_].p); 1708 f->normal_ = plane.n; 1709 1710 if (Maths::isLess(f->normal_*probe.p, f->normal_*rs_->atom_[v1->atom_].p)) 1711 { 1712 f->normal_.negate(); 1713 } 1714 1715 f->singular_ = Maths::isLess(GetDistance(probe.p, plane), probe.radius); 1716 } 1717 faceExists(RSFace * face,const std::list<RSVertex * > & vertices)1718 RSFace* RSComputer::faceExists(RSFace* face, const std::list<RSVertex*>& vertices) 1719 { 1720 std::list<RSVertex*>::const_iterator v; 1721 RSFace* f; 1722 for (v = vertices.begin(); v != vertices.end(); v++) 1723 { 1724 f = (*v)->has(face); 1725 if (f != NULL) 1726 { 1727 return f; 1728 } 1729 } 1730 return NULL; 1731 } 1732 centerOfProbe(const SortedPosition3 & pos,TVector3<double> & c1,TVector3<double> & c2)1733 bool RSComputer::centerOfProbe(const SortedPosition3& pos, TVector3<double>& c1, TVector3<double>& c2) 1734 { 1735 1736 bool back = false; 1737 HashMap<SortedPosition3, ProbePosition* >::Iterator pp = probe_positions_.find(pos); 1738 if (pp != probe_positions_.end()) 1739 { 1740 if (pp->second != NULL) 1741 { 1742 c1 = pp->second->point[0]; 1743 c2 = pp->second->point[1]; 1744 back = true; 1745 } 1746 } else { 1747 TSphere3<double> s1(rs_->atom_[pos.a]); 1748 TSphere3<double> s2(rs_->atom_[pos.b]); 1749 TSphere3<double> s3(rs_->atom_[pos.c]); 1750 1751 s1.radius += rs_->probe_radius_; 1752 s2.radius += rs_->probe_radius_; 1753 s3.radius += rs_->probe_radius_; 1754 1755 if (GetIntersection(s1, s2, s3, c1, c2, false)) 1756 { 1757 ProbePosition* position = new ProbePosition; 1758 position->status[0] = STATUS_NOT_TESTED; 1759 position->status[1] = STATUS_NOT_TESTED; 1760 position->point[0] = c1; 1761 position->point[1] = c2; 1762 probe_positions_.insert(std::make_pair(pos, position)); 1763 back = true; 1764 } 1765 else 1766 { 1767 probe_positions_.insert(std::make_pair(pos, (ProbePosition*)NULL)); 1768 } 1769 } 1770 1771 return back; 1772 } 1773 checkProbe(const TSphere3<double> & probe,const SortedPosition3 & pos)1774 bool RSComputer::checkProbe(const TSphere3<double>& probe, const SortedPosition3& pos) 1775 { 1776 Position index; 1777 ProbePosition* position = probe_positions_[pos]; 1778 if (probe.p == position->point[0]) 1779 { 1780 index = 0; 1781 } 1782 else 1783 { 1784 index = 1; 1785 } 1786 1787 if (position->status[index] == STATUS_NOT_TESTED) 1788 { 1789 bool ok = true; 1790 std::deque<Index> atom_list; 1791 neighboursOfThreeAtoms(pos.a, pos.b, pos.c, atom_list); 1792 double dist; 1793 std::deque<Index>::iterator i = atom_list.begin(); 1794 while (ok && (i != atom_list.end())) 1795 { 1796 dist = probe.radius+rs_->atom_[*i].radius; 1797 if (Maths::isLess(probe.p.getSquareDistance(rs_->atom_[*i].p), dist*dist)) 1798 { 1799 position->status[index] = STATUS_NOT_OK; 1800 ok = false; 1801 } 1802 i++; 1803 } 1804 if (ok) 1805 { 1806 position->status[index] = STATUS_OK; 1807 } 1808 } 1809 return (position->status[index] == STATUS_OK); 1810 } 1811 correctProbePosition(Position atom)1812 void RSComputer::correctProbePosition(Position atom) 1813 { 1814 HashMap<SortedPosition3, ProbePosition* >::Iterator pp; 1815 for (pp = probe_positions_.begin(); pp != probe_positions_.end(); ++pp) 1816 { 1817 if ((pp->first.a == atom) || (pp->first.b == atom) || (pp->first.c == atom)) 1818 { 1819 correctProbePosition(pp->first); 1820 } 1821 } 1822 } 1823 correctProbePosition(const SortedPosition3 & pos)1824 void RSComputer::correctProbePosition(const SortedPosition3& pos) 1825 { 1826 TSphere3<double> s1(rs_->atom_[pos.a]); 1827 s1.radius += rs_->probe_radius_; 1828 1829 TSphere3<double> s2(rs_->atom_[pos.b]); 1830 s2.radius += rs_->probe_radius_; 1831 1832 TSphere3<double> s3(rs_->atom_[pos.c]); 1833 s3.radius += rs_->probe_radius_; 1834 1835 TVector3<double> c1, c2; 1836 if (GetIntersection(s1, s2, s3, c1, c2)) 1837 { 1838 ProbePosition* position = probe_positions_[pos]; 1839 if (position == NULL) 1840 { 1841 position = probe_positions_[pos] = new ProbePosition; 1842 } 1843 position->status[0] = STATUS_NOT_TESTED; 1844 position->status[1] = STATUS_NOT_TESTED; 1845 position->point[0] = c1; 1846 position->point[1] = c2; 1847 } 1848 else 1849 { 1850 delete probe_positions_[pos]; 1851 probe_positions_[pos] = NULL; 1852 } 1853 } 1854 preProcessing()1855 void RSComputer::preProcessing() 1856 { 1857 rs_->r_max_ = rs_->atom_[0].radius; 1858 1859 double x_min = rs_->atom_[0].p.x; 1860 double y_min = rs_->atom_[0].p.y; 1861 double z_min = rs_->atom_[0].p.z; 1862 1863 double x_max = x_min; 1864 double y_max = y_min; 1865 double z_max = z_min; 1866 1867 for (Position i = 1; i < rs_->number_of_atoms_; i++) 1868 { 1869 rs_->r_max_ = std::max(rs_->r_max_, rs_->atom_[i].radius); 1870 1871 x_min = std::min(x_min, rs_->atom_[i].p.x); 1872 y_min = std::min(y_min, rs_->atom_[i].p.y); 1873 z_min = std::min(z_min, rs_->atom_[i].p.z); 1874 1875 x_max = std::max(x_max, rs_->atom_[i].p.x); 1876 y_max = std::max(y_max, rs_->atom_[i].p.y); 1877 z_max = std::max(z_max, rs_->atom_[i].p.z); 1878 } 1879 1880 rs_->bounding_box_.set(x_min, y_min, z_min, x_max, y_max, z_max); 1881 1882 double dist = 2*(rs_->r_max_ + rs_->probe_radius_); 1883 1884 Position nx = (Position)((x_max-x_min)/dist+5); 1885 Position ny = (Position)((y_max-y_min)/dist+5); 1886 Position nz = (Position)((z_max-z_min)/dist+5); 1887 1888 Vector3 origin(x_min - 2*dist, y_min - 2*dist, z_min - 2*dist); 1889 HashGrid3<Position> grid(origin, nx, ny, nz, dist); 1890 1891 HashGridBox3<Position>* box; 1892 HashGridBox3<Position>::ConstBoxIterator b; 1893 HashGridBox3<Position>::ConstDataIterator d; 1894 1895 std::list<Position> to_delete; 1896 Size num_deleted = 0; 1897 Vector3 pos; 1898 for (Position i = 0; i < rs_->number_of_atoms_; i++) 1899 { 1900 pos.set(rs_->atom_[i].p.x, rs_->atom_[i].p.y, rs_->atom_[i].p.z); 1901 1902 // remove atoms that are fully contained in another atom 1903 double radius_i = rs_->atom_[i].radius; 1904 bool too_close = false; 1905 box = grid.getBox(pos); 1906 if (box->isEmpty()) { 1907 continue; 1908 } 1909 for (b = box->beginBox(); b != box->endBox() && !too_close; b++) 1910 { 1911 for (d = b->beginData(); d != b->endData() && !too_close; d++) 1912 { 1913 double radius_d = rs_->atom_[*d].radius; 1914 1915 // our algorithm becomes instable somewhere if two atoms are too close... 1916 // TODO: fix it so this safe guard is no longer necessary 1917 if (Maths::isLessOrEqual(rs_->atom_[i].p.getDistance(rs_->atom_[*d].p), 0.05*std::max(radius_d, radius_i))) 1918 { 1919 too_close = true; 1920 to_delete.push_back(i); 1921 num_deleted++; 1922 if (radius_i > radius_d) 1923 { 1924 rs_->atom_[*d].p = rs_->atom_[i].p; 1925 rs_->atom_[*d].radius = rs_->atom_[i].radius; 1926 } 1927 } 1928 } 1929 } 1930 1931 if (!too_close) 1932 grid.insert(pos, i-num_deleted); 1933 } 1934 1935 for (std::list<Position>::reverse_iterator si = to_delete.rbegin(); si != to_delete.rend(); ++si) 1936 { 1937 rs_->atom_.erase(rs_->atom_.begin()+*si); 1938 rs_->number_of_atoms_--; 1939 } 1940 1941 double offset; 1942 1943 for (Position i = 0; i < rs_->number_of_atoms_-1; i++) 1944 { 1945 offset = rs_->atom_[i].radius + 2*rs_->probe_radius_; 1946 pos.set(rs_->atom_[i].p.x, rs_->atom_[i].p.y, rs_->atom_[i].p.z); 1947 1948 box = grid.getBox(pos); 1949 for (b = box->beginBox(); b != box->endBox(); b++) 1950 { 1951 for (d = b->beginData(); d != b->endData(); d++) 1952 { 1953 // we only need to count every pair twice 1954 if (*d > i) 1955 { 1956 TSphere3<double> const& next_atom = rs_->atom_[*d]; 1957 1958 double dist = next_atom.p.getSquareDistance(rs_->atom_[i].p); 1959 double max_dist = next_atom.radius+offset; 1960 1961 max_dist *= max_dist; 1962 if (!Maths::isGreater(dist, max_dist)) 1963 { 1964 neighbours_[i].push_back(*d); 1965 neighbours_[*d].push_back(i); 1966 } 1967 } 1968 } 1969 } 1970 1971 sort(neighbours_[i].begin(), neighbours_[i].end()); 1972 } 1973 } 1974 insert(RSVertex * vertex)1975 void RSComputer::insert(RSVertex* vertex) 1976 { 1977 rs_->insert(vertex); 1978 new_vertices_.insert(vertex); 1979 vertices_[vertex->atom_].push_back(vertex); 1980 atom_status_[vertex->atom_] = STATUS_ON_SURFACE; 1981 } 1982 insert(RSEdge * edge)1983 void RSComputer::insert(RSEdge* edge) 1984 { 1985 rs_->insert(edge); 1986 edge->vertex_[0]->edges_.insert(edge); 1987 edge->vertex_[1]->edges_.insert(edge); 1988 } 1989 insert(RSFace * face)1990 void RSComputer::insert(RSFace* face) 1991 { 1992 rs_->insert(face); 1993 new_faces_.insert(face); 1994 1995 face->vertex_[0]->faces_.insert(face); 1996 face->vertex_[1]->faces_.insert(face); 1997 face->vertex_[2]->faces_.insert(face); 1998 1999 RSEdge* edge = face->edge_[0]; 2000 edge->vertex_[0]->edges_.insert(edge); 2001 edge->vertex_[1]->edges_.insert(edge); 2002 2003 edge = face->edge_[1]; 2004 edge->vertex_[0]->edges_.insert(edge); 2005 edge->vertex_[1]->edges_.insert(edge); 2006 2007 edge = face->edge_[2]; 2008 edge->vertex_[0]->edges_.insert(edge); 2009 edge->vertex_[1]->edges_.insert(edge); 2010 } 2011 2012 } // namespace BALL 2013