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
C:/Documents and Settings/sdecherchi/My Documents/Ricerca/software nostro/NanoShaper 0.3.1/NanoShaper/src/SkinSurface_official.h
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; // these are only pointer, no explicit denstructor
00029 
00030         virtual ~MixedCell()
00031         {               
00032                 if (quadric!=NULL)
00033                         deleteVector<double>(quadric);
00034         }
00035 };
00036 
00037 
00038 // voronoi facet cell or del edge
00039 class Del1Cell : public MixedCell
00040 {
00041 public:
00042         int ids[2]; // the two atoms generating the del edge
00043         vector<double*> planes; // lateral planes
00044 
00045         virtual ~Del1Cell()
00046         {               
00047                 for (unsigned int i=0;i<planes.size();i++)
00048                         free(planes[i]);
00049         }
00050 
00051         // #IMPROVE one could build pointers from here to the planes computed in Del2Cell or viceversa
00052         double upper[4]; // upper plane
00053         double lower[4]; // upper plane
00054 };
00055 
00056 // reduced voronoi cell
00057 class Del0Cell : public MixedCell
00058 {
00059 public:
00060         int id; // atom index
00061         vector<double*> planes; // planes pointer that points to upper/lower of Del1Cell
00062 
00063         virtual ~Del0Cell()
00064         {               
00065                 // do nothing, planes are not managed by this class
00066         }
00067 };
00068 
00069 // del facet cell
00070 class Del2Cell : public MixedCell
00071 {
00072 public:
00073         int ids[3]; // atom indices
00074         double planes[3][4]; // lateral planes
00075         double upper[4]; 
00076         double lower[4]; 
00077 
00078         virtual ~Del2Cell()
00079         {               
00080                 // do nothing, no pointers
00081         }
00082 };
00083 
00084 // reduced thethraedron
00085 class Del3Cell: public MixedCell
00086 {
00087 public:
00088         int ids[4]; // atom indices
00089         // planes of the reduced tethraedron
00090         double planes[4][4];
00091         double vor[3]; // Voronoi center
00092         double reduced[4][3]; // reduced points in the same order of atom indices
00093 
00094         virtual ~Del3Cell()
00095         {               
00096                 // do nothing, no pointers
00097         }
00098 };
00099 
00100 
00101 #define GRID_MIXEDCELLMAP_2D(i,j,k,NA,NB) gridMixedCellMap2D[(k)+(MAX_MIXEDCELLS_2D-1)*((j)+(NB)*(i))]
00102 // fortran like
00103 //#define GRIDMIXEDCELLMAP(i,j,k,l,NX,NY,NZ) gridMixedCellMap[(i)+(NX)*((j)+(NY)*((k)+(NZ)*(l)))]
00104 // c like (NX is ignored)
00105 #define GRIDMIXEDCELLMAP(i,j,k,l,NX,NY,NZ) gridMixedCellMap[(l)+(MAX_MIXEDCELLS-1)*((k)+(NZ)*((j)+(NY)*(i)))]
00106 
00107 #ifdef ENABLE_CGAL 
00108 
00109 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
00110 #include <CGAL/Regular_triangulation_3.h>
00111 #include <CGAL/Regular_triangulation_euclidean_traits_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 <cassert>
00118 #include <vector>
00119 #include <fstream>
00120 
00121 // use to translate in Pov-Ray
00122 #include <CGAL/Polyhedron_3.h>
00123 #include <CGAL/convex_hull_3.h>
00124         
00126 #endif
00127 
00128 #define DELAUNAY_POINT_CELL             0
00129 #define DELAUNAY_EDGE_CELL              1
00130 #define DELAUNAY_FACET_CELL             2
00131 #define DELAUNAY_TETRA_CELL             3
00132 
00141 class SkinSurface: public Surface
00142 {
00143 
00144 #ifdef ENABLE_CGAL 
00145 private:
00146         
00147         typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
00148         typedef CGAL::Regular_triangulation_filtered_traits_3<K>    Traits;
00149         typedef CGAL::Triangulation_cell_base_with_info_3< Del3Cell*,Traits> Cell_with_info;
00150         typedef CGAL::Regular_triangulation_cell_base_3<Traits,Cell_with_info> RT_Cell_with_info;
00151         typedef CGAL::Triangulation_vertex_base_3<Traits> Vertex;
00152         typedef CGAL::Triangulation_data_structure_3<Vertex,RT_Cell_with_info> tds;
00153         typedef Traits::RT                                          Weight;
00154         typedef Traits::Bare_point                                  Point3;
00155         typedef Traits::Weighted_point                              Weighted_point;
00156         typedef CGAL::Regular_triangulation_3<Traits,tds>           Rt;
00157         typedef Rt::Vertex_iterator                                 Vertex_iterator;
00158         typedef Rt::Finite_vertices_iterator                        Finite_Vertex_Iterator;
00159         typedef Rt::Finite_cells_iterator                                                       Finite_Cells_Iterator;
00160         typedef Rt::Finite_edges_iterator                                                       Finite_Edges_Iterator;
00161         typedef Rt::Finite_facets_iterator                                                      Finite_Facets_Iterator;
00162         typedef Rt::Vertex_handle                                   Vertex_handle;
00163         typedef tds::Cell_circulator                                                            Cell_circulator;
00164         typedef tds::Cell_handle                                                                        Cell_handle;
00165         typedef Rt::Facet                                                                                       Facet;
00166 
00167         // used to translate in Pov-Ray
00168         typedef CGAL::Polyhedron_3<K>                                           Polyhedron;
00169         typedef Polyhedron::Halfedge_around_facet_circulator        HF_circulator;
00170 
00171         //a functor computing the plane containing a triangular facet
00172         struct Plane_from_facet {
00173                 Polyhedron::Plane_3 operator()(Polyhedron::Facet& f) {
00174                 Polyhedron::Halfedge_handle h = f.halfedge();
00175                 return Polyhedron::Plane_3( h->vertex()->point(),
00176                                     h->next()->vertex()->point(),
00177                                     h->opposite()->vertex()->point());
00178         }
00179 };
00180 
00181 
00182 #endif
00183 
00184 private:
00186         int numMixedCells;
00192         int type[4];
00194         int* gridMixedCellMap;
00196         int nx,ny,nz;
00197         double scale;
00198         double side;
00199         double*x,*y,*z;
00200         double s;
00201         double xmin,xmax,ymin,ymax,zmin,zmax;
00202         short*** ind;
00203         int* gridMixedCellMap2D;
00204 
00206         int nx_2d,ny_2d,nz_2d;
00207         double scale_2d;
00208         double side_2d;
00209         double xmin_2d,xmax_2d,ymin_2d,ymax_2d,zmin_2d,zmax_2d;
00210         int** ind_2d;
00211 
00213         MixedCell** mixedComplex;       
00216         bool buildSkin();
00218         bool buildAuxiliaryGrid();
00219         #ifdef ENABLE_CGAL
00220 
00221                 bool buildSkinCGAL();
00222         #endif
00223 
00224         bool savePovRay;
00225 
00226 public:
00228         SkinSurface();
00230         SkinSurface(DelPhiShared* ds);                  
00231 
00233 
00234         virtual bool build();
00236         virtual bool save(char* fileName);
00238         virtual bool load(char* fileName);
00240         virtual void printSummary();            
00242         virtual bool getProjection(double p[3],double* proj1,double* proj2,
00243                 double* proj3,double* normal1,double* normal2,double* normal3);
00249         virtual void getRayIntersection(double p1[3],double p2[3],vector<pair<double,double*> >& intersections,int thdID);              
00251         virtual void preProcessPanel();
00253 
00254         void setShrinking(double ss) 
00255         { 
00256                 if (ss<=1.0 && ss>0.0) 
00257                 { 
00258                         s=ss; 
00259                 }
00260                 else
00261                 {
00262                         cout << endl << WARN << "Cannot set " << ss << ". s parameter is in (0,1]. Setting "<< DEFAULT_S;
00263                         s = DEFAULT_S;
00264                 }
00265         }
00266         
00267         double getShrinking() 
00268         {
00269                 return s;
00270         }
00271 
00272         void setSavePovRay(bool ss) 
00273         { 
00274                 savePovRay = ss;
00275         }
00276         
00277         bool getSavePovRay() 
00278         {
00279                 return savePovRay;
00280         }
00281         
00282         virtual ~SkinSurface();
00283 
00284 private:
00287         //bool rayQuadricIntersection(double*,double*,double**,double*,double*,int thdID,double* cache);
00288         bool rayQuadricIntersection(double*,double*,double*,double*,double*,int thdID,double* cache);
00290         bool isFeasible(MixedCell* mc,double* point);
00292         void projectToQuadric(double* y,double* Q,int type,double* proj,double* norm,double& dist);
00294         void saveSkinPatch(ofstream& of,MixedCell* mc,int index,vector<Weighted_point>& la);
00295 };
00296 
00297 #endif