BALL
1.4.1
|
00001 // -*- Mode: C++; tab-width: 2; -*- 00002 // vi: set ts=2: 00003 // 00004 00005 #ifndef BALL_STRUCTURE_REDUCEDSURFACE_H 00006 #define BALL_STRUCTURE_REDUCEDSURFACE_H 00007 00008 #ifndef BALL_MATHC_COMMON_H 00009 # include <BALL/MATHS/common.h> 00010 #endif 00011 00012 #ifndef BALL_MATHS_SIMPLEBOX3_H 00013 # include <BALL/MATHS/simpleBox3.h> 00014 #endif 00015 00016 #ifndef BALL_MATHS_CIRCLE3_H 00017 # include <BALL/MATHS/circle3.h> 00018 #endif 00019 00020 #ifndef BALL_MATHS_SPHERE_H 00021 # include <BALL/MATHS/sphere3.h> 00022 #endif 00023 00024 #ifndef BALL_MATHS_VECTOR3_H 00025 # include <BALL/MATHS/vector3.h> 00026 #endif 00027 00028 #ifndef BALL_DATATYPE_HASHSET_H 00029 # include <BALL/DATATYPE/hashMap.h> 00030 #endif 00031 00032 #ifndef BALL_DATATYPE_HASHSET_H 00033 # include <BALL/DATATYPE/hashSet.h> 00034 #endif 00035 00036 #ifndef BALL_COMMON_EXCEPTION_H 00037 # include <BALL/COMMON/exception.h> 00038 #endif 00039 00040 #ifndef BALL_STRUCTURE_RSEDGE_H 00041 # include <BALL/STRUCTURE/RSEdge.h> 00042 #endif 00043 00044 #ifndef BALL_STRUCTURE_RSFACE_H 00045 # include <BALL/STRUCTURE/RSFace.h> 00046 #endif 00047 00048 #ifndef BALL_STRUCTURE_RSVERTEX_H 00049 # include <BALL/STRUCTURE/RSVertex.h> 00050 #endif 00051 00052 #include <set> 00053 #include <list> 00054 #include <deque> 00055 #include <vector> 00056 00057 namespace BALL 00058 { 00059 struct SortedPosition2 00060 { 00061 SortedPosition2(Position a1, Position a2) 00062 : a(a1), b(a2) 00063 { 00064 if (a > b) std::swap(a, b); 00065 } 00066 00067 bool operator==(const SortedPosition2& pos) const 00068 { 00069 return (a == pos.a) && (b == pos.b); 00070 } 00071 00072 Position a; 00073 Position b; 00074 }; 00075 00076 struct SortedPosition3 00077 { 00078 SortedPosition3(Position a1, Position a2, Position a3) 00079 : a(a1), b(a2), c(a3) 00080 { 00081 if (a > b) std::swap(a, b); 00082 if (a > c) std::swap(a, c); 00083 if (b > c) std::swap(b, c); 00084 } 00085 00086 bool operator==(const SortedPosition3& pos) const 00087 { 00088 return (a == pos.a) && (b == pos.b) && (c == pos.c); 00089 } 00090 00091 Position a; 00092 Position b; 00093 Position c; 00094 }; 00095 } 00096 00097 #ifdef BALL_HAS_BOOST_UNORDERED_MAP 00098 namespace boost 00099 { 00100 #endif 00101 00102 #ifdef BALL_EXTEND_HASH_IN_STD_NS 00103 namespace std 00104 { 00105 #endif // BALL_EXTEND_HASH_IN_STD_NS 00106 00107 #ifdef BALL_HAS_TR1_UNORDERED_MAP 00108 namespace tr1 00109 { 00110 #endif // BALL_HAS_TR1_UNORDERED_MAP 00111 template<> 00112 struct hash<BALL::SortedPosition2> : public std::unary_function<BALL::SortedPosition2, size_t> 00113 { 00114 inline size_t operator()(const BALL::SortedPosition2& p) const 00115 { 00116 return 5323 * p.a + 1847 * p.b; 00117 } 00118 }; 00119 00120 template<> 00121 struct hash<BALL::SortedPosition3> : public std::unary_function<BALL::SortedPosition3, size_t> 00122 { 00123 inline size_t operator()(const BALL::SortedPosition3& p) const 00124 { 00125 return 5323 * p.a + 1847 * p.b + 2347 * p.c; 00126 } 00127 }; 00128 #ifdef BALL_HAS_TR1_UNORDERED_MAP 00129 } 00130 #endif // BALL_HAS_TR1_UNORDERED_MAP 00131 00132 #ifdef BALL_EXTEND_HASH_IN_STD_NS 00133 } 00134 #endif // BALL_EXTEND_HASH_IN_STD_NS 00135 00136 #ifdef BALL_HAS_BOOST_UNORDERED_MAP 00137 } 00138 #endif 00139 00140 namespace BALL 00141 { 00142 class RSComputer; 00143 class SolventExcludedSurface; 00144 class SESComputer; 00145 class SESSingularityCleaner; 00146 class TriangulatedSES; 00147 class SolventAccessibleSurface; 00148 class TriangulatedSAS; 00149 class SESTriangulator; 00150 00154 class BALL_EXPORT ReducedSurface 00155 { 00156 public: 00157 00170 friend class RSComputer; 00171 friend class SolventExcludedSurface; 00172 friend class SESComputer; 00173 friend class SESSingularityCleaner; 00174 friend class SolventAccessibleSurface; 00175 friend class TriangulatedSES; 00176 friend class TriangulatedSAS; 00177 friend class SESTriangulator; 00178 00179 BALL_CREATE(ReducedSurface) 00180 00181 00184 00189 ReducedSurface(); 00190 00195 ReducedSurface(const ReducedSurface& reduced_surface, bool = true); 00196 00200 ReducedSurface(const std::vector< TSphere3<double> >& spheres, 00201 const double& probe_radius); 00202 00205 virtual ~ReducedSurface(); 00206 00208 00211 00215 void operator = (const ReducedSurface& reduced_surface); 00216 00220 void set(const ReducedSurface& reduced_surface); 00221 00224 void clear(); 00225 00228 void clean(); 00229 00231 00234 00238 Size numberOfAtoms() const; 00239 00243 Size numberOfVertices() const; 00244 00248 Size numberOfEdges() const; 00249 00253 Size numberOfFaces() const; 00254 00258 double getProbeRadius() const; 00259 00264 TSphere3<double> getSphere(Position i) const 00265 throw(Exception::IndexOverflow); 00266 00271 RSVertex* getVertex(Position i) const 00272 throw(Exception::IndexOverflow); 00273 00278 RSEdge* getEdge(Position i) const 00279 throw(Exception::IndexOverflow); 00280 00285 RSFace* getFace(Position i) const 00286 throw(Exception::IndexOverflow); 00287 00291 void insert(RSVertex* rsvertex); 00292 00296 void insert(RSEdge* rsedge); 00297 00301 void insert(RSFace* rsface); 00302 00306 double getMaximalRadius() const; 00307 00311 TSimpleBox3<double> getBoundingBox() const; 00312 00317 void deleteSimilarFaces(RSFace* face1, RSFace* face2); 00318 00327 bool getAngle(RSFace* face1, RSFace* face2, 00328 RSVertex* vertex1, RSVertex* vertex2, 00329 TAngle<double>& angle, bool check = false) const; 00330 00333 void compute() 00334 throw(Exception::GeneralException, 00335 Exception::DivisionByZero, 00336 Exception::IndexOverflow); 00337 00339 00340 private: 00341 00342 /*_ Test whether a ReducedSurface object can be copied. 00343 */ 00344 bool canBeCopied(const ReducedSurface& reduced_surface); 00345 00346 /*_ Copy a ReducedSurface object. 00347 */ 00348 void copy(const ReducedSurface& reduced_surface); 00349 00350 /*_ 00351 */ 00352 void correctEdges(RSFace* face1, RSFace* face2, 00353 RSEdge* edge1, RSEdge* edge2); 00354 00355 /*_ 00356 */ 00357 void joinVertices(RSFace* face1, RSFace* face2, 00358 RSVertex* vertex1, RSVertex* vertex2); 00359 00360 /*_ 00361 */ 00362 void findSimilarVertices(RSFace* face1, RSFace* face2, 00363 std::vector<RSVertex*>& rsvertex1, 00364 std::vector<RSVertex*>& rsvertex2); 00365 00366 /*_ 00367 */ 00368 void findSimilarEdges(RSFace* face1, RSFace* face2, 00369 std::vector<RSEdge*>& rsedge1, 00370 std::vector<RSEdge*>& rsedge2); 00371 00372 protected: 00373 00374 /*_ the number of atoms of the reduced surface 00375 */ 00376 Size number_of_atoms_; 00377 00378 /*_ the atoms of the molecule 00379 */ 00380 00381 std::vector< TSphere3<double> > atom_; 00382 00383 /*_ probe radius 00384 */ 00385 double probe_radius_; 00386 00387 /*_ the number of vertices of the reduced surface 00388 */ 00389 Size number_of_vertices_; 00390 00391 /*_ the vertices of the reduced surface 00392 */ 00393 std::vector< RSVertex* > vertices_; 00394 00395 /*_ the number of edges of the reduced surface 00396 */ 00397 Size number_of_edges_; 00398 00399 /*_ the edges of the reduced surface 00400 */ 00401 std::vector< RSEdge* > edges_; 00402 00403 /*_ the number of faces of the reduced surface 00404 */ 00405 Size number_of_faces_; 00406 00407 /*_ the faces of the reduced surface 00408 */ 00409 std::vector< RSFace* > faces_; 00410 00411 /*_ maximal radius of all atoms 00412 */ 00413 double r_max_; 00414 00415 /*_ bounding SimpleBox of the atom centers of the molecule 00416 */ 00417 TSimpleBox3<double> bounding_box_; 00418 }; 00419 00423 00427 BALL_EXPORT std::ostream& operator << (std::ostream& s, const ReducedSurface& rs); 00428 00430 00434 class BALL_EXPORT RSComputer 00435 { 00436 public: 00437 00438 BALL_CREATE(RSComputer) 00439 00440 00443 00449 enum ProbeStatus 00450 { 00451 STATUS_OK = 0, 00452 STATUS_NOT_OK, 00453 STATUS_NOT_TESTED 00454 }; 00455 00461 enum AtomStatus 00462 { 00463 STATUS_ON_SURFACE = 0, 00464 STATUS_INSIDE, 00465 STATUS_UNKNOWN 00466 }; 00468 00469 struct ProbePosition 00470 { 00471 ProbeStatus status[2]; 00472 TVector3<double> point[2]; 00473 }; 00474 00478 00483 RSComputer(); 00484 00487 RSComputer(ReducedSurface* rs); 00488 00491 virtual ~RSComputer(); 00492 00494 00497 00500 void run() 00501 throw(Exception::GeneralException, 00502 Exception::DivisionByZero, 00503 Exception::IndexOverflow); 00504 00506 00507 00508 private: 00509 00510 /*_ @name Computing reduced surface 00511 */ 00513 00514 /*_ 00515 */ 00516 void preProcessing(); 00517 00518 /*_ Compute a RSComponent. 00519 */ 00520 void getRSComponent() 00521 throw(Exception::GeneralException, 00522 Exception::DivisionByZero, 00523 Exception::IndexOverflow); 00524 00525 /*_ Treat all edges of a face. 00526 @param face the RSFace to be treated 00527 */ 00528 bool treatFace(RSFace* face) 00529 throw(Exception::GeneralException, 00530 Exception::DivisionByZero, 00531 Exception::IndexOverflow); 00532 00533 /*_ Roll over an edge that belongs to only one face and find the other one. 00534 @param edge the RSEdge to be treated 00535 */ 00536 bool treatEdge(RSEdge* edge) 00537 throw(Exception::GeneralException, 00538 Exception::DivisionByZero, 00539 Exception::IndexOverflow); 00540 00541 /*_ Treat an ambiguous situation. 00542 All vertices on an ambiguous atom are deleted with all its edges and 00543 faces. The radius of the atom is decreased by 10 EPSILON. 00544 @param atom the index of the atom 00545 */ 00546 void correct(Index atom); 00547 00548 /*_ Check all new created vertices for extensions 00549 */ 00550 void extendComponent() 00551 throw(Exception::GeneralException, 00552 Exception::DivisionByZero, 00553 Exception::IndexOverflow); 00554 00555 /*_ Find a third atom rolling over two vertices starting on a face. 00556 From all atoms which can be touched by the probe sphere when it 00557 touches the given two vertices we choose the one with smallest 00558 rotation angle. 00559 If the rotation angle equals zero, the probe sphere can touch four 00560 atoms and an exception is thrown. 00561 If no atom can be found an exception is thrown. 00562 @param vertex1 the first vertex 00563 @param vertex2 the second vertex 00564 @param face the starting face 00565 @param probe the new probe sphere 00566 @param phi the rotation angle 00567 @return Index index of the found atom 00568 */ 00569 Index thirdAtom(RSVertex* vertex1, RSVertex* vertex2, 00570 RSFace* face, TSphere3<double>& probe, TAngle<double>& phi) 00571 throw(Exception::GeneralException, 00572 Exception::DivisionByZero, 00573 Exception::IndexOverflow); 00574 00576 /*_ @name Finding a start position 00577 */ 00579 00580 /*_ Find a start position 00581 @param vertex a pointer to the found vertex, if only a vertex 00582 can be found 00583 @param edge a pointer to the found edge, if only an edge can be 00584 found 00585 @param face a pointer to the found face, if a face can be found 00586 @return Position 0, if no start position is found, 00587 1, if a single vertex is found, 00588 2, if an edge is found, 00589 3, if a face is found 00590 */ 00591 Position getStartPosition() 00592 throw(Exception::DivisionByZero); 00593 00595 /*_ @name Finding a first face 00596 */ 00598 00599 /*_ Try to find a starting face 00600 @return RSFace* a pointer to the found face, if a face can be found, 00601 NULL otherwise 00602 */ 00603 RSFace* findFirstFace() 00604 throw(Exception::DivisionByZero); 00605 00606 /*_ Try to find a starting face in a given direction 00607 @param direction search in x-direction, if direction is 0, 00608 search in y-direction, if direction is 1, 00609 search in z-direction, if direction is 2 00610 @param extrem search in min direction, if extrem is 0, 00611 search in max direction, if extrem is 1 00612 @return RSFace* a pointer to the found face, if a face can be found, 00613 NULL otherwise 00614 */ 00615 RSFace* findFace(Position direction, Position extrem) 00616 throw(Exception::DivisionByZero); 00617 00619 /*_ @name Finding a first edge 00620 */ 00622 00623 /*_ Try to find a starting edge 00624 @return RSEdge* a pointer to the found edge, if a face can be found, 00625 NULL otherwise 00626 */ 00627 RSEdge* findFirstEdge(); 00628 00629 /*_ Try to find a starting edge in a given direction 00630 @param direction search in x-direction, if direction is 0, 00631 search in y-direction, if direction is 1, 00632 search in z-direction, if direction is 2 00633 @param extrem search in min direction, if extrem is 0, 00634 search in max direction, if extrem is 1 00635 @return RSEdge* a pointer to the found edge, if a face can be found, 00636 NULL otherwise 00637 */ 00638 RSEdge* findEdge(Position direction, Position extrem); 00639 00641 /*_ @name Finding a first vertex 00642 */ 00644 00645 /*_ Try to find a single atom 00646 @return RSVertex* a pointer to the found vertex, if a vertex can be 00647 found, NULL otherwise 00648 */ 00649 RSVertex* findFirstVertex(); 00650 00651 /*_ Find a single atom in a given direction 00652 @param direction search in x-direction, if direction is 0, 00653 search in y-direction, if direction is 1, 00654 search in z-direction, if direction is 2 00655 @param extrem search in min direction, if extrem is 0, 00656 search in max direction, if extrem is 1 00657 @return Index the index of the found atom 00658 */ 00659 Index findFirstAtom(Position direction, Position extrem); 00660 00661 /*_ Find a second atom close enougth to the first atom in a given direction 00662 @param atom1 the index of the first atom 00663 @param direction search in x-direction, if direction is 0, 00664 search in y-direction, if direction is 1, 00665 search in z-direction, if direction is 2 00666 @param extrem search in min direction, if extrem is 0, 00667 search in max direction, if extrem is 1 00668 @return Index the index of the found atom 00669 */ 00670 Index findSecondAtom(Index atom, Position direction, Position extrem); 00671 00672 /*_ Find a second atom close enougth to the first two atoms 00673 @param atom1 the index of the first atom 00674 @param atom2 the index of the second atom 00675 @param atom_list a HashSet of the indices of all candidate atoms 00676 @return ::std::list< ::std::pair< Index,TSphere3<double> > > 00677 a list of all candidates with their probe spheres 00678 */ 00679 void findThirdAtom(Index atom1, Index atom2, const std::deque<Index>& third, 00680 std::deque< std::pair< Index,TSphere3<double> > >& atoms); 00681 00683 /*_ @name Some utilities 00684 */ 00686 00687 /*_ Find all atoms close enougth to two given atoms. 00688 The indices of all atoms which can be touched by the probe sphere when 00689 it touches the given atoms are computed. 00690 @param atom1 the index of the first given atom 00691 @param atom2 the index of the second given atom 00692 @param output_list list of all atoms close enougth to the given atoms 00693 */ 00694 const std::deque<Index>& neighboursOfTwoAtoms(const SortedPosition2& pos); 00695 00696 /*_ Find all atoms close enougth to three given atoms. 00697 The indices of all atoms which can be touched by the probe sphere when 00698 it touches the given atoms are computed. 00699 @param atom1 the index of the first given atom 00700 @param atom2 the index of the second given atom 00701 @param atom3 the index of the third given atom 00702 @param output_list list of all atoms close enougth to the given atoms 00703 */ 00704 void neighboursOfThreeAtoms(Index atom1, Index atom2, Index atom3, 00705 std::deque<Index>& output_list); 00706 00707 /*_ Get the extrem coordinate of a circle in a given direction 00708 @param circle the circle 00709 @param direction search in x-direction, if direction is 0, 00710 search in y-direction, if direction is 1, 00711 search in z-direction, if direction is 2 00712 @param extrem search in min direction, if extrem is 0, 00713 search in max direction, if extrem is 1 00714 @return double the extrem coordinate 00715 */ 00716 double getCircleExtremum(const TCircle3<double>& circle, 00717 Position direction, Position extrem); 00718 00720 /*_ @name Creating / updating edges / faces 00721 */ 00723 00724 /*_ Create a free edge from two vertices if it is a free edge 00725 @param vertex1 a pointer to the first vertex 00726 @param vertex2 a pointer to the second vertex 00727 @return RSEdge* a pointer to the created free edge, if there is one, 00728 NULL otherwise 00729 */ 00730 RSEdge* createFreeEdge(RSVertex* vertex1, RSVertex* vertex2); 00731 00732 /*_ Get the circle described by the center of the probe sphere and the two 00733 contact circles with the atoms when the probe sphere rolls over two 00734 atoms 00735 @param atom1 the index of the first atom 00736 @param atom2 the index of the second atom 00737 @param circle1 the circle described by the center of the probe sphere 00738 @param circle2 the contact circle with atom1 00739 @param circle3 the contact circle with atom2 00740 @return bool true, if the probe sphere can touch both atoms, 00741 false, otherwise 00742 */ 00743 bool getCircles(Index atom1, Index atom2, TCircle3<double>& circle1, 00744 TCircle3<double>& circle2, TCircle3<double>& circle3); 00745 00746 /*_ Get the normal vector of the face described by three atoms and a probe 00747 @param atom1 the index of the first atom 00748 @param atom2 the index of the second atom 00749 @param atom3 the index of the third atom 00750 @param probe the probe sphere lying on atom1, atom2, atom3 00751 @return TVector3<double> the normal vector 00752 */ 00753 TVector3<double> getFaceNormal(const TSphere3<double>& atom1, const TSphere3<double>& atom2, 00754 const TSphere3<double>& atom3, const TSphere3<double>& probe); 00755 00756 /*_ Update a face and it's edges 00757 @param v1 the first vertex of the face 00758 @param v2 the second vertex of the face 00759 @param v3 the third vertex of the face 00760 @param e1 the first edge 00761 @param e2 the second edge 00762 @param e3 the third edge 00763 @param f the face 00764 @param probe the probe sphere of the face 00765 */ 00766 void updateFaceAndEdges(RSVertex* v1, RSVertex* v2, RSVertex* v3, 00767 RSEdge* e1, RSEdge* e2, RSEdge* e3, 00768 RSFace* f, const TSphere3<double>& probe); 00769 00770 /*_ Test, whether a face exists or not 00771 @param face a pointer to the face to be tested 00772 @return RSFace* a pointer to the face, if it exists, otherwise NULL 00773 */ 00774 RSFace* faceExists(RSFace* face, const std::list< RSVertex* >& vertices); 00775 00777 /*_ @name Finding a probe sphere 00778 */ 00780 00781 /*_ Get the centers of the probe sphere when it lies on three atoms 00782 @param pos the three atom's indices 00783 @param c1 the first center 00784 @param c2 the second center 00785 @return bool true, if the probe sphere can touch the three atoms, 00786 false, otherwise 00787 */ 00788 bool centerOfProbe(const SortedPosition3& pos, TVector3<double>& c1, TVector3<double>& c2); 00789 00790 /*_ Check,weather a probe sphere is inside an atom 00791 @param probe the probe sphere to be tested 00792 @return bool true, if the probe sphere is not intersecting any atom 00793 false, otherwise 00794 */ 00795 bool checkProbe(const TSphere3<double>& probe, const SortedPosition3& pos); 00796 00797 /*_ 00798 */ 00799 void correctProbePosition(Position atom); 00800 00801 /*_ 00802 */ 00803 void correctProbePosition(const SortedPosition3& pos); 00804 00805 /*_ 00806 */ 00807 void insert(RSVertex* vertex); 00808 00809 /*_ 00810 */ 00811 void insert(RSEdge* edge); 00812 00813 /*_ 00814 */ 00815 void insert(RSFace* face); 00816 00818 00819 protected: 00820 00821 00822 /*_ a pointer to the reduced surface to compute 00823 */ 00824 ReducedSurface* rs_; 00825 00826 /*_ for each atom a list of its neighbours 00827 */ 00828 std::vector< std::deque<Index> > neighbours_; 00829 00830 /*_ for each atom a status 00831 */ 00832 std::vector< AtomStatus > atom_status_; 00833 00834 /*_ for each pair of atoms a list of its neighbours 00835 */ 00836 HashMap< SortedPosition2, std::deque<Index> > neighbours_of_two_; 00837 00838 /*_ for each triple of atoms its probe positions 00839 */ 00840 HashMap< SortedPosition3, ProbePosition* > probe_positions_; 00841 00842 /*_ all new created vertices which are not yet checked for extensions 00843 */ 00844 HashSet<RSVertex*> new_vertices_; 00845 00846 /*_ all new created faces which are not completely treated yet 00847 */ 00848 HashSet<RSFace*> new_faces_; 00849 00850 /*_ for each atom a list of the rsvertices of the atom 00851 */ 00852 std::vector< std::list<RSVertex*> > vertices_; 00853 }; 00854 } // namespace BALL 00855 00856 #endif // BALL_STRUCTURE_REDUCEDSURFACE_H