NanoShaper
0.7.2
NanoShaper is a tool able to triangulate and inspect an arbitray triangulated surface or several types of molecular surfaces
|
00001 00002 //--------------------------------------------------------- 00006 //--------------------------------------------------------- 00007 00008 #ifndef Surface_h 00009 #define Surface_h 00010 00011 #include "octree.h" 00012 #include "globals.h" 00013 #include "ConfigFile.h" 00014 #include "SurfaceFactory.h" 00015 00016 #ifdef ENABLE_CGAL 00017 00018 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> 00019 #include <CGAL/Regular_triangulation_3.h> 00020 #include <CGAL/Regular_triangulation_euclidean_traits_3.h> 00021 #include <CGAL/Regular_triangulation_filtered_traits_3.h> 00022 #include <CGAL/Triangulation_cell_base_with_info_3.h> 00023 #include <CGAL/Regular_triangulation_cell_base_3.h> 00024 #include <CGAL/Triangulation_vertex_base_with_info_3.h> 00025 #include <CGAL/Triangulation_data_structure_3.h> 00027 00028 // vor point container 00029 class VorPoint 00030 { 00031 public: 00032 double vor[3]; // Voronoi center 00033 bool visited; 00034 00035 VorPoint() 00036 { 00037 visited = false; 00038 } 00039 }; 00040 #endif 00041 00042 00044 class coordVec 00045 { 00046 public: 00047 00048 coordVec(const int x,const int y,const int z,double* const v,const int d) 00049 { 00050 ix = x; 00051 iy = y; 00052 iz = z; 00053 vec = v; 00054 dir = d; 00055 } 00056 00058 coordVec(const coordVec& cv) 00059 { 00060 ix = cv.ix; 00061 iy = cv.iy; 00062 iz = cv.iz; 00063 vec = cv.vec; 00064 dir = cv.dir; 00065 } 00066 00067 int ix; 00068 int iy; 00069 int iz; 00070 int dir; 00071 double* vec; 00072 }; 00073 00074 class packet 00075 { 00076 public: 00077 packet() 00078 {} 00079 packet(const packet& p) 00080 { 00081 first = p.first; 00082 second = p.second; 00083 } 00084 vector<coordVec>* first; 00085 vector<coordVec>* second; 00086 }; 00087 00088 00089 #ifdef DBGMEM_CRT 00090 #define _CRTDBG_MAP_ALLOC 00091 #define _CRTDBG_MAP_ALLOC_NEW 00092 #endif 00093 00094 #include "tools.h" 00095 #include "DelphiShared.h" 00096 00097 // ids for rays directions 00098 #define X_DIR 0 00099 #define Y_DIR 1 00100 #define Z_DIR 2 00101 00102 // ids for panels 00103 #define PANEL_YZ 0 00104 #define PANEL_XY 1 00105 #define PANEL_XZ 2 00106 00107 extern int num_cores; 00108 00109 #define DEFAULT_VOLUME 11.4 00110 // tollerance on two intersections. If less than EPS_INT then the intersection is the same. def 1e-8 00111 #define EPS_INT 1e-8 00112 #define MAX_ATOMS_MULTI_GRID 100 00113 //#define DEBUG_SURFACE 00114 00115 #define DELTA 1e-10 00116 #define RAND_DISPLACEMENT 1e-4 00117 00118 // boundary grid point codes 00119 #define INTERNAL_BGP 0 00120 #define EXTERNAL_BGP 1 00121 00122 // molecular surface 00123 #define MOLECULAR_SURFACE 0 00124 // analytical object 00125 #define OBJECT_SURFACE 3 00126 // unknown surface, e.g., mesh 00127 #define GENERIC_SURFACE 2 00128 // 'sum' of two different surface types 00129 #define HYBRID_SURFACE 4 00130 00131 // file format 00132 #define OFF 0 00133 #define OFF_A 1 00134 #define OFF_N 2 00135 #define OFF_N_A 3 00136 #define MSMS_NO_A 4 00137 #define MSMS 5 00138 00139 // now using inline templates 00140 //#define GRID_MULTI_MAP(i,j,k,l,NX,NY,NZ) gridMultiMap[(l)+(MAX_ATOMS_MULTI_GRID)*((k)+(NZ)*((j)+(NY)*(i)))] 00141 00142 #ifdef ENABLE_BOOST_THREADS 00143 00144 #define THREAD_SAFE_SCOPE(SSS) \ 00145 { \ 00146 boost::mutex::scoped_lock scopedLock(mutex); \ 00147 SSS \ 00148 cout.flush(); \ 00149 } 00150 00151 #else 00152 00153 #define THREAD_SAFE_SCOPE(SSS) \ 00154 { \ 00155 SSS \ 00156 cout.flush(); \ 00157 } 00158 #endif 00159 00160 00161 00162 00196 class Surface 00197 { 00198 #ifdef ENABLE_CGAL 00199 private: 00200 00201 typedef CGAL::Exact_predicates_inexact_constructions_kernel _K; 00202 typedef CGAL::Regular_triangulation_filtered_traits_3<_K> _Traits; 00203 typedef CGAL::Triangulation_cell_base_with_info_3< VorPoint*,_Traits> _Cell_with_info; 00204 typedef CGAL::Regular_triangulation_cell_base_3<_Traits,_Cell_with_info> _RT_Cell_with_info; 00205 typedef CGAL::Triangulation_vertex_base_3<_Traits> _Vertex; 00206 typedef CGAL::Triangulation_data_structure_3<_Vertex,_RT_Cell_with_info> _tds; 00207 typedef _Traits::RT _Weight; 00208 typedef _Traits::Bare_point _Point3; 00209 typedef _Traits::Weighted_point _Weighted_point; 00210 typedef CGAL::Regular_triangulation_3<_Traits,_tds> _Rt; 00211 typedef _Rt::Vertex_iterator _Vertex_iterator; 00212 typedef _Rt::Finite_vertices_iterator _Finite_Vertex_Iterator; 00213 typedef _Rt::Finite_cells_iterator _Finite_Cells_Iterator; 00214 typedef _Rt::Vertex_handle _Vertex_handle; 00215 typedef _tds::Cell_circulator _Cell_circulator; 00216 typedef _tds::Cell_handle _Cell_handle; 00217 #endif 00218 00219 protected: 00220 00222 00223 int surfType; 00227 bool isRCbased; 00230 bool providesAnalyticalNormals; 00232 00234 bool projBGP; 00235 double sternLayer; 00236 bool accurateTriangulation; 00237 bool fillCavitiesFlag; 00238 bool computeNormals; 00240 bool checkDuplicatedVertices; 00243 bool wellShaped; 00245 double probe_radius; 00246 bool useLoadBalancing; 00247 bool vertexAtomsMapFlag; 00248 bool saveMSMS; 00250 00251 // current panel under analysis 00252 int panel; 00253 // delphi environment 00254 DelPhiShared* delphi; 00255 double totalSurfaceArea; 00256 double totalVolume; 00257 00258 // last row and column sizes of ind_2d matrix for ray casting acceleration 00259 int last_rows_ind; 00260 int last_cols_ind; 00262 double randDisplacement; 00263 // last nx,ny,nz dimensions seen by Surface class 00264 int last_nx,last_ny,last_nz; 00266 Octree<int>* intersectionsMatrixAlongX; 00268 Octree<int>* intersectionsMatrixAlongY; 00270 Octree<int>* intersectionsMatrixAlongZ; 00272 Octree<int>* normalsMatrixAlongX; 00274 Octree<int>* normalsMatrixAlongY; 00276 Octree<int>* normalsMatrixAlongZ; 00277 00279 bool* activeCubes; 00280 00281 bool*** verticesInsidenessMap; 00282 //bool *verticesInsidenessMap; 00283 00284 double*** scalarField; 00285 00286 // when a scalar field is available scalarField is used instead of verticesInsidenessMap for vertex interpolation 00287 bool isAvailableScalarField; 00288 00290 00291 vector<int*> triList; 00293 vector<double*> vertList; 00295 vector<double*> normalsList; 00297 00298 double delta_accurate_triangulation; 00299 00301 int* bgp_type; 00302 00304 int* gridMultiMap; 00305 double gxmin,gymin,gzmin,gside,gscale; 00306 unsigned int ggrid; 00307 00309 int* gridLoad; 00310 int totalLoad; 00311 int* vertexAtomsMap; 00312 00314 int inside; 00315 00316 #ifdef ENABLE_BOOST_THREADS 00317 boost::mutex mutex; 00318 #endif 00319 00321 00323 void floodFill(int ix,int iy,int iz,int idold,int idnew); 00324 00326 void floodFill2(int ix,int iy,int iz,int idold,int idnew); 00327 00329 void floodFill4(int ix,int iy,int iz,int idold,int idnew); 00330 00332 void floodFill3( pair<pair<int,int>,int > ind, 00333 pair<int,int> z_limits, 00334 pair<int,int> old_new, 00335 pair< queue< pair<pair<int,int>,int> > *, queue< pair<pair<int,int>,int > > * > queues, 00336 queue< pair<pair<int,int>,int> > * in_queue); 00337 00344 bool innerFloodFill(int* start,int* target,short* status1,short* status2,int& maxMoves,bool debug); 00345 00347 bool vdwAccessible(double* p,int& nearest); 00348 00357 void intersector(double* volPanel,int nb,int start,int end,int jump,int* numIntersections,packet pack); 00358 00361 void projector(int start,int end); 00362 00365 char getInsidness(int i, int j, int k, int vertInd); 00366 00368 void vertexInterp(double isolevel, double* p1, double* p2, double valp1, double valp2,double* p); 00369 00371 //int getTriangles(double* vertexValues,double** vertexPos,double isolevel, int** triangles,int ix,int iy,int iz, 00372 // int NX, int NY,int NZ,int* xedge,int* yedge,int* zedge,int* xedge_down,int* yedge_down); 00373 int getTriangles(double* vertexValues,double** vertexPos,double isolevel, int** triangles,int ix,int iy,int iz, 00374 int NX, int NY,int NZ); 00375 00377 void getVertices(double isolevel,int start_z,int end_z,int jump,vector<coordVec>*,vector<double*>*); 00378 00380 int classifyCube(double* vertexValues,double isolevel); 00381 00383 void approximateNormals(vector<int>& appNormals,bool doOnlyList); 00384 00386 double triangulationKernel(double isolevel,bool revert,int start_z,int end_z,int jump,vector<int*>* localTriList,double* localArea); 00387 00389 void buildAtomsMap(); 00390 00392 void disposeAtomsMap(); 00393 00396 void applyMultidielectric(); 00397 00399 void swap2multi(double gxmin,double gymin,double gzmin,double gside,unsigned int ggrid,int* gridMultiMap,int i,int j,int k,int l); 00400 00402 void buildSternLayer(); 00403 00405 void allocIntersectionsMatrices(); 00406 00408 void allocNormalsMatrices(); 00409 00411 int deduceFormat(); 00412 00413 protected: 00414 00415 Surface(); 00416 Surface(ConfigFile* cf); 00417 00418 public: 00419 00420 // In order to use the set of default surface methods, the following 00421 // interface methods must be provided. 00422 // If a fully custom solution is built the mandatory methods can be fake methods 00423 // and all the work can be made in getSurf(); 00424 00426 00429 virtual bool build() = 0; 00430 00433 virtual bool save(char* fileName) = 0; 00434 00437 virtual bool load(char* fileName) = 0; 00438 00440 virtual void printSummary() = 0; 00441 00443 virtual bool getProjection(double p[3],double* proj1,double* proj2, 00444 double* proj3,double* normal1,double* normal2,double* normal3)=0; 00445 00451 virtual void getRayIntersection(double p1[3],double p2[3],vector<pair<double,double*> >& intersections,int thdID,bool computeNormals)=0; 00452 00454 virtual void init() = 0; 00455 00457 virtual void init(ConfigFile* cf) = 0; 00458 00460 virtual void clear() = 0; 00462 00464 00473 virtual bool getSurf(bool fill=false,double vol=0); 00474 00477 virtual double getVolume(); 00478 00480 virtual double getArea(); 00481 00485 virtual int getCavities(int idStart=STATUS_POINT_TEMPORARY_OUT); 00486 00490 virtual void fillCavities(double vol=0,bool silent=false); 00491 00494 virtual void filterCavities(bool modStatus=false); 00495 00497 virtual void cav2out(); 00498 00504 virtual double triangulateSurface(double iso=0.0,const char* fileName="triangulatedSurf",bool revert=false); 00505 00507 virtual bool saveMesh(int format,bool revert,const char* fileName); 00508 00511 virtual void smoothSurface(const char* fn="triangulatedSurf",bool revert=false); 00512 00515 virtual void preProcessPanel() 00516 {} 00517 00520 virtual void postRayCasting() 00521 {} 00522 00524 virtual bool preBoundaryProjection() 00525 { 00526 return true; 00527 } 00528 00530 void tri2Balls(); 00531 00533 void backupStatus(); 00534 00536 void removeBackupStatus(); 00537 00539 int linkCavities(short* st1,short* st2); 00540 00545 bool difference(Surface*); 00546 00549 Surface& operator-=(Surface& surf2); 00550 00552 void getCavitiesAtoms(); 00553 00555 00558 void setProjBGP(bool flag) 00559 { 00560 projBGP = flag; 00561 } 00562 bool getProjBGP() 00563 { 00564 return projBGP; 00565 } 00566 double getSternLayer() 00567 { 00568 return sternLayer; 00569 } 00570 void setSternLayer(double l) 00571 { 00572 if (l<0) 00573 { 00574 cout << endl << WARN << "Cannot set a negative Stern Layer"; 00575 return; 00576 } 00577 else 00578 sternLayer = l; 00579 } 00580 00581 void setSaveMSMS(bool m) 00582 { 00583 saveMSMS = m; 00584 } 00585 00586 bool getSaveMSMS() 00587 { 00588 return saveMSMS; 00589 } 00590 00591 void setTriangulationFlag(bool flag) 00592 { 00593 accurateTriangulation = flag; 00594 } 00595 bool getTriangulationFlag() 00596 { 00597 return accurateTriangulation; 00598 } 00599 00600 void setVertexAtomsMap(bool f) 00601 { 00602 vertexAtomsMapFlag = f; 00603 } 00604 00605 bool getVertexAtomsMap() 00606 { 00607 return vertexAtomsMapFlag; 00608 } 00609 00610 void setComputeNormals(bool cn) 00611 { 00612 computeNormals = cn; 00613 } 00614 00615 bool getComputeNormals() 00616 { 00617 return computeNormals; 00618 } 00619 00620 void setCheckDuplicatedVertices(bool cd) 00621 { 00622 checkDuplicatedVertices = cd; 00623 } 00624 bool getCheckDuplicatedVertices() 00625 { 00626 return checkDuplicatedVertices; 00627 } 00628 void setKeepWellShapedCavities(bool kwsc) 00629 { 00630 wellShaped = kwsc; 00631 } 00632 bool getKeepWellShapedCavities() 00633 { 00634 return wellShaped; 00635 } 00636 virtual void setProbeRadius(double probeRadius) 00637 { 00638 probe_radius = probeRadius; 00639 } 00640 virtual double getProbeRadius() 00641 { 00642 return probe_radius; 00643 } 00644 00645 virtual double getRandDisplacement() 00646 { 00647 return randDisplacement; 00648 } 00649 00650 virtual void setRandDisplacement(const double r) 00651 { 00652 randDisplacement = r; 00653 } 00654 00655 virtual void setLoadBalancing(bool doLoadBalancing) 00656 { 00657 useLoadBalancing = doLoadBalancing; 00658 } 00659 00660 virtual bool getLoadBalancing() 00661 { 00662 return useLoadBalancing; 00663 } 00664 00665 virtual int getNumTriangles() 00666 { 00667 return triList.size(); 00668 } 00669 00670 virtual int getNumVertices() 00671 { 00672 return vertList.size(); 00673 } 00674 00675 virtual void setInsideCode(int i) 00676 { 00677 inside = i; 00678 } 00679 00680 virtual int getInsideCode() 00681 { 00682 return inside; 00683 } 00684 00686 00688 virtual ~Surface(); 00689 }; 00690 00692 00693 #endif