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
|
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