NanoShaper  0.3.1
NanoShaper is a tool able to triangulate and inspect an arbitray triangulated surface or several types of molecular surfaces such as the Gaussian, Skin and the SES
C:/Documents and Settings/sdecherchi/My Documents/Ricerca/software nostro/NanoShaper 0.3.1/NanoShaper/src/ConnollySurface.h
Go to the documentation of this file.
00001 //---------------------------------------------------------
00005 //---------------------------------------------------------
00006 
00007 #ifndef ConnollySurface_h
00008 #define ConnollySurface_h
00009 
00010 #include "Surface.h"
00011 
00012 //#define DEBUG_CONNOLLY
00013 
00014 #define MAX_INCIDENT_PROBES 50
00015 #define MAX_PROBES 50
00016 
00017 #define DEFAULT_PROBE_RADIUS 1.4 // default probe radius
00018 
00019 
00020 class ConnollyCell{
00021 public: 
00022         int patch_type;
00023 
00024         virtual ~ConnollyCell()
00025         {}
00026 };
00027 
00028 class FacetCell: public ConnollyCell
00029 {
00030 public:
00031         double center[3];
00032         int id[3]; // atoms triplet indices
00033         bool isSelfIntersecting;        
00034         double planes[4][4]; // trimming tethraedron
00035         int plane_indices[3][2]; // every plane depends on a pair of atoms
00036         vector<double*> self_intersection_planes; // additional self intersection planes
00037         FacetCell* mirrorCell;  
00038 
00039         FacetCell()
00040         {
00041                 mirrorCell=NULL;
00042         }
00043 
00044         virtual ~FacetCell()
00045         {
00046                 for (unsigned int i=0;i<self_intersection_planes.size();i++)
00047                         deleteVector<double>(self_intersection_planes[i]);
00048         }
00049 };
00050 
00051 class EdgeCell: public ConnollyCell
00052 {
00053 public:
00054         double center[3];
00055         int id[2]; // atom pair indices
00056         bool isSelfIntersecting;
00057         double clipping_center[3]; // center of the clipping sphere of the torus
00058         double clipping_radius;
00059         double major_radius;
00060         double u[3],v[3];
00061         double self_intersection_radius; // radius of the clipping sphere in case of spindle torus
00062         double Rot[3][3]; // rotation matrix of the torus
00063         double invrot[3][3]; // inverse rotation matrix
00064         double cutting_planes[2][4]; // two cutting planes of the torus
00065         vector<double*> additional_planes; // additional clipping planes due to tori clipping by singular facets
00066         vector<bool> flags; // acute / non acute flags for additional planes    
00067         bool acute; // if angle between planes is acute perform "and" side test between planes, otherwise "or" side test is sufficient
00068         double rcoi;
00069 
00070         virtual ~EdgeCell()
00071         {
00072                 for (unsigned int i=0;i<additional_planes.size();i++)
00073                         deleteVector<double>(additional_planes[i]);
00074         }
00075 };
00076 
00077 
00078 class PointCell: public ConnollyCell
00079 {
00080 public:
00081         int id; // atom index
00082         vector<EdgeCell*> neighbours;
00083         vector<FacetCell*> incidentProbes;
00084         vector<EdgeCell*> buried_neighbours;
00085 
00086         virtual ~PointCell()
00087         {
00088                 for (unsigned int i=0;i<buried_neighbours.size();i++)
00089                         delete buried_neighbours[i];
00090         }
00091 
00092 };
00093 
00094 // row major
00095 // 3d map
00096 #define GRID_CONNOLLY_CELL_MAP(i,j,k,l,NX,NY,NZ) gridConnollyCellMap[(l)+(MAX_CONNOLLY_CELLS-1)*((k)+(NZ)*((j)+(NY)*(i)))]
00097 // 2d map
00098 #define GRID_CONNOLLY_CELL_MAP_2D(i,j,k,NA,NB) gridConnollyCellMap2D[(k)+(MAX_CONNOLLY_CELLS_2D-1)*((j)+(NB)*(i))]
00099 
00100 #define SELF_MAP(i,j,k,l,NX,NY,NZ) gridProbesMap[(l)+(MAX_PROBES-1)*((k)+(NZ)*((j)+(NY)*(i)))]
00101 
00102 #ifdef ENABLE_CGAL 
00103 
00104 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
00105 #include <CGAL/Regular_triangulation_3.h>
00106 #include <CGAL/Regular_triangulation_filtered_traits_3.h>
00107 #include <CGAL/Triangulation_cell_base_with_info_3.h>
00108 #include <CGAL/Regular_triangulation_cell_base_3.h>
00109 #include <CGAL/Triangulation_vertex_base_with_info_3.h>
00110 #include <CGAL/Triangulation_data_structure_3.h>
00111 #include <CGAL/Fixed_alpha_shape_3.h>
00112 #include <CGAL/Fixed_alpha_shape_vertex_base_3.h>
00113 #include <CGAL/Fixed_alpha_shape_cell_base_3.h>
00114 #include <CGAL/Polyhedron_3.h>
00115 #include <CGAL/convex_hull_3.h>
00116 #include <CGAL/Vector_3.h>
00117 #include <CGAL/Plane_3.h>
00118 #include <CGAL/Ray_3.h>
00119 #include <CGAL/Line_3.h>
00120 #include <CGAL/intersections.h>
00121 #include <cassert>
00122 #include <vector>
00123 #include <fstream>
00124 
00126 #endif
00127 
00128 // here singular/regular is in the alpha shape nomenclature
00129 
00130 #define POINT_CELL                              0
00131 #define REGULAR_EDGE_CELL               1
00132 #define SINGULAR_EDGE_CELL              2
00133 #define REGULAR_FACE_CELL               3
00134 #define SINGULAR_FACE_CELL              4
00135 
00143 class ConnollySurface: public Surface
00144 {
00145 
00146 #ifdef ENABLE_CGAL 
00147 private:
00148         typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
00149         typedef CGAL::Regular_triangulation_filtered_traits_3<K>        Traits;
00150         typedef CGAL::Polyhedron_3<K>                                           Polyhedron;     
00151         typedef Polyhedron::Halfedge_around_facet_circulator        HF_circulator;
00152         typedef CGAL::Vector_3<K>                                                                       Vector3;
00153         typedef CGAL::Plane_3<K>                                                                        Plane3;
00154         typedef CGAL::Ray_3<K>                                                                      Ray3;
00155         typedef CGAL::Line_3<K>                                                                     Line3;      
00156         typedef CGAL::Triangulation_vertex_base_with_info_3<PointCell*,Traits> Vertex_with_info;                
00157         typedef CGAL::Fixed_alpha_shape_vertex_base_3<Traits,Vertex_with_info>          Vb;                     
00158         typedef CGAL::Triangulation_cell_base_with_info_3< FacetCell*[4],Traits> Cell_with_info;
00159         typedef CGAL::Regular_triangulation_cell_base_3<Traits,Cell_with_info> RT_Cell_with_info;
00160         typedef CGAL::Fixed_alpha_shape_cell_base_3<Traits,RT_Cell_with_info>         Fb;
00161         typedef CGAL::Triangulation_data_structure_3<Vb,Fb>                     Tds;
00162         typedef Tds::Cell_circulator                                                            Cell_circulator;
00163         typedef Tds::Cell_handle                                                                        Cell_handle;
00164         typedef CGAL::Regular_triangulation_3<Traits,Tds>                       Rt;
00165         typedef CGAL::Fixed_alpha_shape_3<Rt>                                   Fixed_alpha_shape_3;
00166         typedef Fixed_alpha_shape_3::Cell_handle                                        Alpha_Cell_handle;
00167         typedef Fixed_alpha_shape_3::Vertex_handle                                      Alpha_Vertex_handle;
00168         typedef Fixed_alpha_shape_3::Facet                                                      Alpha_Facet;
00169         typedef Fixed_alpha_shape_3::Edge                                                       Alpha_Edge;
00170         typedef Traits::Bare_point                                  Point3;
00171         typedef Traits::Weighted_point                              Weighted_point;
00172 
00173         typedef Rt::Vertex_iterator                                 Vertex_iterator;
00174         typedef Rt::Finite_vertices_iterator                        Finite_Vertex_Iterator;
00175         typedef Rt::Finite_cells_iterator                                                       Finite_Cells_Iterator;
00176         typedef Rt::Finite_edges_iterator                                                       Finite_Edges_Iterator;
00177         typedef Rt::Finite_facets_iterator                                                      Finite_Facets_Iterator;
00178         typedef Rt::Vertex_handle                                   Vertex_handle;
00179         typedef Rt::Facet                                                                                       Facet;
00180         
00181 
00182 #endif
00183 
00184 private:
00186         int type[5];
00188         int* gridConnollyCellMap2D;
00189         int* gridConnollyCellMap;       
00191         unsigned int nx,ny,nz;
00192         double scale;
00193         double side;
00195         unsigned int nx_2d,ny_2d,nz_2d;
00196         double scale_2d;
00197         double side_2d;
00198         double xmin_2d,xmax_2d,ymin_2d,ymax_2d,zmin_2d,zmax_2d;
00199         unsigned int** ind_2d;
00200 
00201         double*x,*y,*z;
00202         double s;
00203         double xmin,xmax,ymin,ymax,zmin,zmax;
00204         unsigned short*** ind;
00205 
00206         unsigned int AUX_GRID_DIM_CONNOLLY;
00207         unsigned int MAX_CONNOLLY_CELLS;
00208         unsigned int AUX_GRID_DIM_CONNOLLY_2D;
00209         unsigned int MAX_CONNOLLY_CELLS_2D;
00210 
00212         vector<ConnollyCell*> sesComplex;       
00214         PointCell** atomPatches;
00217         bool buildConnolly();
00219         bool buildAuxiliaryGrid();
00220 
00221         #ifdef ENABLE_CGAL
00222 
00223                 bool buildConnollyCGAL();
00224         #endif
00225 
00226         bool savePovRay;
00227 
00228 public:
00230         ConnollySurface();
00232         ConnollySurface(DelPhiShared* ds);                      
00233 
00235 
00236         virtual bool build();
00238         virtual bool save(char* fileName);
00240         virtual bool load(char* fileName);
00242         virtual void printSummary();            
00244         virtual bool getProjection(double p[3],double* proj1,double* proj2,
00245                 double* proj3,double* normal1,double* normal2,double* normal3);
00251         virtual void getRayIntersection(double p1[3],double p2[3],vector<pair<double,double*> >& intersections,int thdID);              
00253         virtual void preProcessPanel();
00254         virtual void postRayCasting();
00255         virtual bool preBoundaryProjection();
00257 
00258         void setProbeRadius(double pr) 
00259         { 
00260                 double e = 1e-2;
00261                 if (pr>e) 
00262                 { 
00263                         probe_radius=pr; 
00264                 }
00265                 else
00266                 {
00267                         cout << endl << WARN << "Cannot set " << pr << "<=" << e << ". Setting "<< DEFAULT_PROBE_RADIUS;
00268                         probe_radius = DEFAULT_PROBE_RADIUS;
00269                 }
00270         }
00271         
00272         double getProbeRadius() 
00273         {
00274                 return probe_radius;
00275         }
00276 
00277         void setSavePovRay(bool ss) 
00278         { 
00279                 savePovRay = ss;
00280         }
00281         
00282         bool getSavePovRay() 
00283         {
00284                 return savePovRay;
00285         }
00286 
00288         void setAuxGrid(unsigned int dim,unsigned int max)
00289         {               
00290                 AUX_GRID_DIM_CONNOLLY = dim;
00291                 MAX_CONNOLLY_CELLS = max;
00292         }
00293 
00297         void setAuxGrid2D(unsigned int dim,unsigned int max)
00298         {               
00299                 AUX_GRID_DIM_CONNOLLY_2D = dim;
00300                 MAX_CONNOLLY_CELLS_2D = (max*dim);
00301         }
00302         
00303         virtual ~ConnollySurface();
00304 
00305 private:
00306 
00307         bool rayConnollyCellIntersection(double*,double*,ConnollyCell*,double*,double*,double*,double*,int& numInt,int);
00309         bool isFeasible(ConnollyCell* cc,double* point);
00311         void projectToTorus(double* y,EdgeCell* ec,double* proj,double* norm,double& dist);
00313         void projectToSphere(double* y,double* center,double radius,double* proj,double* norm,double& dist);
00315         void saveConcaveSpherePatch(ofstream& of,FacetCell* fc,int i);
00316         void saveSphere(ostream& of,double* center,double radius);
00317         void saveAtomPatch(ofstream& of,PointCell* pc);
00318         void saveEdgePatch(ofstream& of,EdgeCell* ec,int size,double bigR,double* u,double* v,double* w,bool isComplex);
00320         bool orientation(double* pb_center1,double* pb_center2,double* w1,double* w2);
00322         void sortProbes(EdgeCell* ec,FacetCell** fcv,int np,int* sorted);       
00323         bool raySphere(double* orig,double* dir,double* sphere_center,double sphere_radius,double* t1,double* t2);
00324         void getCoi(double* torus_center,double rcoi,double** sampledPoints,int numPoints,double* u,double* v);
00328         void projectToCircle(double* point,double radius,double* center,double* plane,double* proj,double& dist);
00329 };
00330 
00331 #endif