NanoShaper
0.7.2
NanoShaper is a tool able to triangulate and inspect an arbitray triangulated surface or several types of molecular surfaces
|
00001 //--------------------------------------------------------- 00005 //--------------------------------------------------------- 00006 00007 #ifndef ConnollySurface_h 00008 #define ConnollySurface_h 00009 00010 #include "Surface.h" 00011 #include "SurfaceFactory.h" 00012 00013 #ifdef DBGMEM_CRT 00014 #define _CRTDBG_MAP_ALLOC 00015 #define _CRTDBG_MAP_ALLOC_NEW 00016 #endif 00017 00018 00019 //#define DEBUG_CONNOLLY 00020 00021 #define MAX_INCIDENT_PROBES 50 00022 00023 #define DEFAULT_PROBE_RADIUS 1.4 // default probe radius 00024 00025 00026 class ConnollyCell{ 00027 public: 00028 int patch_type; 00029 00030 virtual ~ConnollyCell() 00031 {} 00032 }; 00033 00034 class FacetCell: public ConnollyCell 00035 { 00036 public: 00037 double center[3]; 00038 int id[3]; // atoms triplet indices 00039 bool isSelfIntersecting; 00040 double planes[4][4]; // trimming tethraedron 00041 int plane_indices[3][2]; // every plane depends on a pair of atoms 00042 vector<double*> self_intersection_planes; // additional self intersection planes 00043 FacetCell* mirrorCell; 00044 00045 FacetCell() 00046 { 00047 mirrorCell=NULL; 00048 } 00049 00050 virtual ~FacetCell() 00051 { 00052 for (unsigned int i=0;i<self_intersection_planes.size();i++) 00053 deleteVector<double>(self_intersection_planes[i]); 00054 } 00055 }; 00056 00057 class EdgeCell: public ConnollyCell 00058 { 00059 public: 00060 double center[3]; 00061 int id[2]; // atom pair indices 00062 bool isSelfIntersecting; 00063 double clipping_center[3]; // center of the clipping sphere of the torus 00064 double clipping_radius; 00065 double major_radius; 00066 double u[3],v[3]; 00067 double self_intersection_radius; // radius of the clipping sphere in case of spindle torus 00068 double Rot[3][3]; // rotation matrix of the torus 00069 double invrot[3][3]; // inverse rotation matrix 00070 double cutting_planes[2][4]; // two cutting planes of the torus 00071 vector<double*> additional_planes; // additional clipping planes due to tori clipping by singular facets 00072 vector<bool> flags; // acute / non acute flags for additional planes 00073 bool acute; // if angle between planes is acute perform "and" side test between planes, otherwise "or" side test is sufficient 00074 double rcoi; 00075 00076 virtual ~EdgeCell() 00077 { 00078 for (unsigned int i=0;i<additional_planes.size();i++) 00079 deleteVector<double>(additional_planes[i]); 00080 } 00081 }; 00082 00083 00084 class PointCell: public ConnollyCell 00085 { 00086 public: 00087 int id; // atom index 00088 vector<EdgeCell*> neighbours; 00089 vector<FacetCell*> incidentProbes; 00090 vector<EdgeCell*> buried_neighbours; 00091 00092 virtual ~PointCell() 00093 { 00094 for (unsigned int i=0;i<buried_neighbours.size();i++) 00095 delete buried_neighbours[i]; 00096 } 00097 00098 }; 00099 00100 // row major 00101 // 3d map 00102 #define GRID_CONNOLLY_CELL_MAP(i,j,k,l,NX,NY,NZ) gridConnollyCellMap[(l)+(MAX_CONNOLLY_CELLS-1)*((k)+(NZ)*((j)+(NY)*(i)))] 00103 // 2d map 00104 #define GRID_CONNOLLY_CELL_MAP_2D(i,j,k,NA,NB) gridConnollyCellMap2D[(k)+(MAX_CONNOLLY_CELLS_2D-1)*((j)+(NB)*(i))] 00105 00106 #define SELF_MAP(i,j,k,l,NX,NY,NZ) gridProbesMap[(l)+(MAX_PROBES-1)*((k)+(NZ)*((j)+(NY)*(i)))] 00107 00108 #ifdef ENABLE_CGAL 00109 00110 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> 00111 #include <CGAL/Regular_triangulation_3.h> 00112 #include <CGAL/Regular_triangulation_filtered_traits_3.h> 00113 #include <CGAL/Triangulation_cell_base_with_info_3.h> 00114 #include <CGAL/Regular_triangulation_cell_base_3.h> 00115 #include <CGAL/Triangulation_vertex_base_with_info_3.h> 00116 #include <CGAL/Triangulation_data_structure_3.h> 00117 #include <CGAL/Fixed_alpha_shape_3.h> 00118 #include <CGAL/Fixed_alpha_shape_vertex_base_3.h> 00119 #include <CGAL/Fixed_alpha_shape_cell_base_3.h> 00120 #include <CGAL/Polyhedron_3.h> 00121 #include <CGAL/convex_hull_3.h> 00122 #include <CGAL/Vector_3.h> 00123 #include <CGAL/Plane_3.h> 00124 #include <CGAL/Ray_3.h> 00125 #include <CGAL/Line_3.h> 00126 #include <CGAL/intersections.h> 00127 #include <cassert> 00128 #include <vector> 00129 #include <fstream> 00131 #endif 00132 00133 // here singular/regular is in the alpha shape nomenclature 00134 00135 #define POINT_CELL 0 00136 #define REGULAR_EDGE_CELL 1 00137 #define SINGULAR_EDGE_CELL 2 00138 #define REGULAR_FACE_CELL 3 00139 #define SINGULAR_FACE_CELL 4 00140 00149 class ConnollySurface: public Surface 00150 { 00151 00152 #ifdef ENABLE_CGAL 00153 private: 00154 typedef CGAL::Exact_predicates_inexact_constructions_kernel K; 00155 typedef CGAL::Regular_triangulation_filtered_traits_3<K> Traits; 00156 typedef CGAL::Polyhedron_3<K> Polyhedron; 00157 typedef Polyhedron::Halfedge_around_facet_circulator HF_circulator; 00158 typedef CGAL::Vector_3<K> Vector3; 00159 typedef CGAL::Plane_3<K> Plane3; 00160 typedef CGAL::Ray_3<K> Ray3; 00161 typedef CGAL::Line_3<K> Line3; 00162 typedef CGAL::Triangulation_vertex_base_with_info_3<PointCell*,Traits> Vertex_with_info; 00163 typedef CGAL::Fixed_alpha_shape_vertex_base_3<Traits,Vertex_with_info> Vb; 00164 typedef CGAL::Triangulation_cell_base_with_info_3< FacetCell*[4],Traits> Cell_with_info; 00165 typedef CGAL::Regular_triangulation_cell_base_3<Traits,Cell_with_info> RT_Cell_with_info; 00166 typedef CGAL::Fixed_alpha_shape_cell_base_3<Traits,RT_Cell_with_info> Fb; 00167 typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds; 00168 typedef Tds::Cell_circulator Cell_circulator; 00169 typedef CGAL::Regular_triangulation_3<Traits,Tds> Rt; 00170 typedef CGAL::Fixed_alpha_shape_3<Rt> Fixed_alpha_shape_3; 00171 typedef Fixed_alpha_shape_3::Cell_handle Alpha_Cell_handle; 00172 typedef Fixed_alpha_shape_3::Vertex_handle Alpha_Vertex_handle; 00173 typedef Fixed_alpha_shape_3::Facet Alpha_Facet; 00174 typedef Fixed_alpha_shape_3::Edge Alpha_Edge; 00175 typedef Traits::Bare_point Point3; 00176 typedef Traits::Weighted_point Weighted_point; 00177 00178 typedef Rt::Vertex_iterator Vertex_iterator; 00179 typedef Rt::Finite_vertices_iterator Finite_Vertex_Iterator; 00180 typedef Rt::Finite_cells_iterator Finite_Cells_Iterator; 00181 typedef Rt::Finite_edges_iterator Finite_Edges_Iterator; 00182 typedef Rt::Finite_facets_iterator Finite_Facets_Iterator; 00183 typedef Rt::Vertex_handle Vertex_handle; 00184 typedef Rt::Facet Facet; 00185 00186 00187 #endif 00188 00189 private: 00191 int type[5]; 00193 int* gridConnollyCellMap2D; 00194 int* gridConnollyCellMap; 00196 unsigned int nx,ny,nz; 00197 double scale; 00198 double side; 00200 unsigned int nx_2d,ny_2d,nz_2d; 00201 double scale_2d; 00202 double side_2d; 00203 double xmin_2d,xmax_2d,ymin_2d,ymax_2d,zmin_2d,zmax_2d; 00204 unsigned int** ind_2d; 00205 00206 double*x,*y,*z; 00207 double s; 00208 double xmin,xmax,ymin,ymax,zmin,zmax; 00209 unsigned short*** ind; 00210 00211 unsigned int AUX_GRID_DIM_CONNOLLY; 00212 unsigned int MAX_CONNOLLY_CELLS; 00213 unsigned int AUX_GRID_DIM_CONNOLLY_2D; 00214 unsigned int MAX_CONNOLLY_CELLS_2D; 00215 00217 vector<ConnollyCell*> sesComplex; 00219 PointCell** atomPatches; 00222 bool buildConnolly(); 00224 bool buildAuxiliaryGrid(); 00225 00226 #ifdef ENABLE_CGAL 00227 00228 bool buildConnollyCGAL(); 00229 #endif 00230 00231 bool savePovRay; 00232 int MAX_PROBES; 00234 double si_perfil; 00235 00236 public: 00238 ConnollySurface(); 00240 ConnollySurface(DelPhiShared* ds); 00242 ConnollySurface(ConfigFile* cf,DelPhiShared* ds); 00243 00245 00246 virtual bool build(); 00248 virtual bool save(char* fileName); 00250 virtual bool load(char* fileName); 00252 virtual void printSummary(); 00254 virtual bool getProjection(double p[3],double* proj1,double* proj2, 00255 double* proj3,double* normal1,double* normal2,double* normal3); 00261 virtual void getRayIntersection(double p1[3],double p2[3],vector<pair<double,double*> >& intersections,int thdID,bool computeNormals); 00263 virtual void init(); 00265 virtual void init(ConfigFile* cf); 00267 virtual void clear(); 00269 virtual void preProcessPanel(); 00270 virtual void postRayCasting(); 00271 virtual bool preBoundaryProjection(); 00273 00274 void setProbeRadius(double pr) 00275 { 00276 double e = 1e-2; 00277 if (pr>e) 00278 { 00279 probe_radius=pr; 00280 } 00281 else 00282 { 00283 cout << endl << WARN << "Cannot set " << pr << "<=" << e << ". Setting "<< DEFAULT_PROBE_RADIUS; 00284 probe_radius = DEFAULT_PROBE_RADIUS; 00285 } 00286 } 00287 00288 double getProbeRadius() 00289 { 00290 return probe_radius; 00291 } 00292 00293 void setSavePovRay(bool ss) 00294 { 00295 savePovRay = ss; 00296 } 00297 00298 bool getSavePovRay() 00299 { 00300 return savePovRay; 00301 } 00302 00304 void setAuxGrid(unsigned int dim,unsigned int max) 00305 { 00306 AUX_GRID_DIM_CONNOLLY = dim; 00307 MAX_CONNOLLY_CELLS = max; 00308 } 00309 00313 void setAuxGrid2D(unsigned int dim,unsigned int max) 00314 { 00315 AUX_GRID_DIM_CONNOLLY_2D = dim; 00316 MAX_CONNOLLY_CELLS_2D = (max*dim); 00317 } 00318 00321 void setMaxProbes(int m) 00322 { 00323 if (m<=0) 00324 { 00325 cout << endl << WARN << "Cannot set max probes <0"; 00326 return; 00327 } 00328 MAX_PROBES = m; 00329 } 00330 00331 int getMaxProbes() 00332 { 00333 return MAX_PROBES; 00334 } 00335 00336 void setSIPerfil(double si) 00337 { 00338 si_perfil = si; 00339 } 00340 00341 double getSIPerfil() 00342 { 00343 return si_perfil; 00344 } 00345 00346 virtual ~ConnollySurface(); 00347 00348 private: 00349 00350 bool rayConnollyCellIntersection(double*,double*,ConnollyCell*,double*,double*,double*,double*,int& numInt,int); 00352 bool isFeasible(ConnollyCell* cc,double* point); 00354 void projectToTorus(double* y,EdgeCell* ec,double* proj,double* norm,double& dist); 00357 void getNormalToTorus(double* y,EdgeCell* ec,double* normal); 00361 void getNormal(double* y,ConnollyCell* cc,double* normal); 00363 void saveConcaveSpherePatch(ofstream& of,FacetCell* fc,int i); 00364 void saveSphere(ostream& of,double* center,double radius); 00365 void saveAtomPatch(ofstream& of,PointCell* pc); 00366 void saveEdgePatch(ofstream& of,EdgeCell* ec,int size,double bigR,double* u,double* v,double* w,bool isComplex); 00368 bool orientation(double* pb_center1,double* pb_center2,double* w1,double* w2); 00370 void sortProbes(EdgeCell* ec,FacetCell** fcv,int np,int* sorted); 00371 void getCoi(double* torus_center,double rcoi,double** sampledPoints,int numPoints,double* u,double* v); 00375 void projectToCircle(double* point,double radius,double* center,double* plane,double* proj,double& dist); 00376 }; 00377 00378 00379 // expand it explicitly because Swig is not able to expand it 00380 static class ConnollySurfaceRegister{ 00381 static Surface* createSurface(ConfigFile* conf,DelPhiShared* ds) 00382 { 00383 return new ConnollySurface(conf,ds); 00384 } 00385 public: 00386 ConnollySurfaceRegister() 00387 { 00388 surfaceFactory().add("ses",createSurface); 00389 } 00390 } ConnollySurfaceRegisterObject; 00391 00392 00393 //static SurfaceRecorder<ConnollySurface> sesRecorder("ses"); 00394 00395 #endif