BALL  1.4.1
reducedSurface.h
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines