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 SkinSurface_h 00008 #define SkinSurface_h 00009 00010 #include "Surface.h" 00011 00012 //#define DEBUG_SKIN 00013 00014 #define AUX_GRID_DIM_SKIN 100 00015 #define MAX_MIXEDCELLS 400 00016 #define DEFAULT_S 0.45 // default shrink factor 00017 00018 // aggressive settings can be used in a 64 bit machine with 6 GB of memory 00019 00020 #define AUX_GRID_DIM_SKIN_2D 50 // aggressive setting is 150 00021 #define MAX_MIXEDCELLS_2D (400*AUX_GRID_DIM_SKIN_2D) // aggressive setting is (200*AUX_GRID_DIM_SKIN_2D) 00022 00024 class MixedCell{ 00025 public: 00026 double* quadric; 00027 int surface_type; 00028 vector<double*> mc_points; 00029 vector<double*> legacy_planes; 00030 }; 00031 00032 00033 // voronoi facet cell or del edge 00034 class Del1Cell : public MixedCell 00035 { 00036 public: 00037 int ids[2]; // the two atoms generating the del edge 00038 vector<double*> planes; // lateral planes 00039 // #IMPROVE one could build pointers from here to the planes computed in Del2Cell 00040 double upper[4]; // upper plane 00041 double lower[4]; // upper plane 00042 }; 00043 00044 // reduced voronoi cell 00045 class Del0Cell : public MixedCell 00046 { 00047 public: 00048 int id; // atom index 00049 vector<double*> planes; // planes pointer that points to upper/lower of Del1Cell 00050 vector<Del1Cell*> neighbours; //Del1 neighbours in the same order of the planes 00051 }; 00052 00053 // del facet cell 00054 class Del2Cell : public MixedCell 00055 { 00056 public: 00057 int ids[3]; // atom indices 00058 double planes[3][4]; // lateral planes 00059 double upper[4]; // pointer to the upper plane stored in a del0 cell 00060 double lower[4]; // pointer to the lower plane stored in a del0 cell 00061 // ids says which are the three atoms involved in the facets 00062 }; 00063 00064 // reduced thethraedron 00065 class Del3Cell: public MixedCell 00066 { 00067 public: 00068 int ids[4]; // atom indices 00069 //Del2Cell* neighbours[4]; // linked del2 cells 00070 // planes of the reduced tethraedron 00071 double planes[4][4]; 00072 double vor[3]; // Voronoi center 00073 double reduced[4][3]; // reduced points in the same order of atom indices 00074 double radius2 ; // if is negative the sphere is immaginary 00075 }; 00076 00077 00078 #define GRID_MIXEDCELLMAP_2D(i,j,k,NA,NB) gridMixedCellMap2D[(k)+(MAX_MIXEDCELLS_2D-1)*((j)+(NB)*(i))] 00079 // fortran like 00080 //#define GRIDMIXEDCELLMAP(i,j,k,l,NX,NY,NZ) gridMixedCellMap[(i)+(NX)*((j)+(NY)*((k)+(NZ)*(l)))] 00081 // c like (NX is ignored) 00082 #define GRIDMIXEDCELLMAP(i,j,k,l,NX,NY,NZ) gridMixedCellMap[(l)+(MAX_MIXEDCELLS-1)*((k)+(NZ)*((j)+(NY)*(i)))] 00083 00084 #ifdef ENABLE_CGAL 00085 00086 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> 00087 #include <CGAL/Regular_triangulation_3.h> 00088 #include <CGAL/Regular_triangulation_euclidean_traits_3.h> 00089 #include <CGAL/Regular_triangulation_filtered_traits_3.h> 00090 #include <CGAL/Triangulation_cell_base_with_info_3.h> 00091 #include <CGAL/Regular_triangulation_cell_base_3.h> 00092 #include <CGAL/Triangulation_vertex_base_with_info_3.h> 00093 #include <CGAL/Triangulation_data_structure_3.h> 00094 #include <CGAL/Polyhedron_3.h> 00095 #include <CGAL/convex_hull_3.h> 00096 #include <CGAL/Vector_3.h> 00097 #include <CGAL/Plane_3.h> 00098 #include <CGAL/Ray_3.h> 00099 #include <CGAL/Line_3.h> 00100 #include <CGAL/intersections.h> 00101 #include <CGAL/Linear_algebraCd.h> 00102 #include <cassert> 00103 #include <vector> 00104 #include <fstream> 00105 00107 #endif 00108 00109 #define DELAUNAY_POINT_CELL 0 00110 #define DELAUNAY_EDGE_CELL 1 00111 #define DELAUNAY_FACET_CELL 2 00112 #define DELAUNAY_TETRA_CELL 3 00113 00122 class SkinSurface: public Surface 00123 { 00124 00125 #ifdef ENABLE_CGAL 00126 private: 00127 00128 typedef CGAL::Exact_predicates_inexact_constructions_kernel K; 00129 typedef CGAL::Regular_triangulation_filtered_traits_3<K> Traits; 00130 00131 typedef CGAL::Triangulation_cell_base_with_info_3< Del3Cell*,Traits> Cell_with_info; 00132 typedef CGAL::Regular_triangulation_cell_base_3<Traits,Cell_with_info> RT_Cell_with_info; 00133 typedef CGAL::Triangulation_vertex_base_3<Traits> Vertex; 00134 typedef CGAL::Triangulation_data_structure_3<Vertex,RT_Cell_with_info> tds; 00135 typedef Traits::RT Weight; 00136 typedef Traits::Bare_point Point3; 00137 typedef Traits::Weighted_point Weighted_point; 00138 typedef CGAL::Polyhedron_3<K> Polyhedron; 00139 typedef CGAL::Regular_triangulation_3<Traits,tds> Rt; 00140 typedef Rt::Vertex_iterator Vertex_iterator; 00141 typedef Rt::Finite_vertices_iterator Finite_Vertex_Iterator; 00142 typedef Rt::Finite_cells_iterator Finite_Cells_Iterator; 00143 typedef Rt::Finite_edges_iterator Finite_Edges_Iterator; 00144 typedef Rt::Finite_facets_iterator Finite_Facets_Iterator; 00145 typedef Rt::Vertex_handle Vertex_handle; 00146 typedef tds::Cell_circulator Cell_circulator; 00147 typedef tds::Cell_handle Cell_handle; 00148 typedef Polyhedron::Halfedge_around_facet_circulator HF_circulator; 00149 typedef Rt::Facet Facet; 00150 typedef CGAL::Linear_algebraCd<double> LA; 00151 00152 #endif 00153 00154 private: 00156 int numMixedCells; 00162 int type[4]; 00164 int* gridMixedCellMap; 00166 int nx,ny,nz; 00167 double scale; 00168 double side; 00169 double*x,*y,*z; 00170 double s; 00171 double xmin,xmax,ymin,ymax,zmin,zmax; 00172 short*** ind; 00173 int* gridMixedCellMap2D; 00174 00176 int nx_2d,ny_2d,nz_2d; 00177 double scale_2d; 00178 double side_2d; 00179 double xmin_2d,xmax_2d,ymin_2d,ymax_2d,zmin_2d,zmax_2d; 00180 int** ind_2d; 00181 00183 MixedCell** mixedComplex; 00186 bool buildSkin(); 00188 bool buildAuxiliaryGrid(); 00189 #ifdef ENABLE_CGAL 00190 00191 bool buildSkinCGAL(); 00192 #endif 00193 00194 bool savePovRay; 00195 00196 public: 00198 SkinSurface(); 00200 SkinSurface(DelPhiShared* ds); 00201 00203 00204 virtual bool build(); 00206 virtual bool save(char* fileName); 00208 virtual bool load(char* fileName); 00210 virtual void printSummary(); 00212 virtual bool getProjection(double p[3],double* proj1,double* proj2, 00213 double* proj3,double* normal1,double* normal2,double* normal3); 00219 virtual void getRayIntersection(double p1[3],double p2[3],vector<pair<double,double*> >& intersections,int thdID); 00221 virtual void preProcessPanel(); 00223 00224 void setShrinking(double ss) 00225 { 00226 if (ss<=1.0 && ss>0.0) 00227 { 00228 s=ss; 00229 } 00230 else 00231 { 00232 cout << endl << WARN << "Cannot set " << ss << ". s parameter is in (0,1]. Setting "<< DEFAULT_S; 00233 s = DEFAULT_S; 00234 } 00235 } 00236 00237 double getShrinking() 00238 { 00239 return s; 00240 } 00241 00242 void setSavePovRay(bool ss) 00243 { 00244 savePovRay = ss; 00245 } 00246 00247 bool getSavePovRay() 00248 { 00249 return savePovRay; 00250 } 00251 00252 virtual ~SkinSurface(); 00253 00254 private: 00257 //bool rayQuadricIntersection(double*,double*,double**,double*,double*,int thdID,double* cache); 00258 bool rayQuadricIntersection(double*,double*,double*,double*,double*,int thdID,double* cache); 00260 //bool isFeasible(int num,double** planes,double* point); 00261 bool SkinSurface::isFeasible(int num,vector<double*> planes,double* point); 00262 00264 void projectToQuadric(double* y,double* Q,int type,double* proj,double* norm,double& dist); 00265 00266 void sortPointsAroundEdge(double* ec,vector<double*> fcv,int np,int* sorted) 00267 { 00268 vector<indexed_double> angles; 00269 double ref[3]; 00270 double tt; 00271 SUB(ref,fcv[0],ec) 00272 NORMALIZE_S(ref,tt); 00273 double ws[3],temp[3],vs[3]; 00274 SUB(temp,fcv[1],ec) 00275 NORMALIZE_S(temp,tt); 00276 CROSS(vs,ref,temp) 00277 CROSS(ws,vs,ref) 00278 00279 for (int i=1;i<np;i++) 00280 { 00281 double temp[3]; 00282 SUB(temp,fcv[i],ec) 00283 NORMALIZE_S(temp,tt) 00284 double a1 = DOT(ws,temp); 00285 double a2 = DOT(ref,temp); 00286 double angle = atan2(a1,a2); 00287 // [0,2pi] 00288 if (angle<0) 00289 angle = TWO_PI+angle; 00290 angles.push_back(indexed_double(angle,i)); 00291 } 00292 // remember indexing and sort 00293 sort(angles.begin(),angles.end(),index_double_comparator); 00294 // save 00295 sorted[0]=0; 00296 for (int i=1;i<np;i++) 00297 sorted[i] = angles[i-1].second; 00298 return; 00299 } 00300 00301 }; 00302 00303 #endif