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