NanoShaper  0.7.2
NanoShaper is a tool able to triangulate and inspect an arbitray triangulated surface or several types of molecular surfaces
C:/Documents and Settings/sdecherchi/My Documents/Ricerca/software nostro/NanoShaper0.7/src/Surface.h
Go to the documentation of this file.
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