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_fast2.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                 deleteVector<double>(quadric);
00033         }
00034 };
00035 
00036 
00037 // voronoi facet cell or del edge
00038 class Del1Cell : public MixedCell
00039 {
00040 public:
00041         int ids[2]; // the two atoms generating the del edge
00042         vector<double*> planes; // lateral planes
00043 
00044         virtual ~Del1Cell()
00045         {               
00046                 for (unsigned int i=0;i<planes.size();i++)
00047                         free(planes[i]);
00048         }
00049 
00050         // #IMPROVE one could build pointers from here to the planes computed in Del2Cell or viceversa
00051         double upper[4]; // upper plane
00052         double lower[4]; // upper plane
00053 };
00054 
00055 // reduced voronoi cell
00056 class Del0Cell : public MixedCell
00057 {
00058 public:
00059         int id; // atom index
00060         vector<double*> planes; // planes pointer that points to upper/lower of Del1Cell
00061 
00062         virtual ~Del0Cell()
00063         {               
00064                 // do nothing, planes are not managed by this class
00065         }
00066 };
00067 
00068 // del facet cell
00069 class Del2Cell : public MixedCell
00070 {
00071 public:
00072         int ids[3]; // atom indices
00073         double planes[3][4]; // lateral planes
00074         double upper[4]; 
00075         double lower[4]; 
00076 
00077         virtual ~Del2Cell()
00078         {               
00079                 // do nothing, no pointers
00080         }
00081 };
00082 
00083 // reduced thethraedron
00084 class Del3Cell: public MixedCell
00085 {
00086 public:
00087         int ids[4]; // atom indices
00088         // planes of the reduced tethraedron
00089         double planes[4][4];
00090         double vor[3]; // Voronoi center
00091         double reduced[4][3]; // reduced points in the same order of atom indices
00092 
00093         virtual ~Del3Cell()
00094         {               
00095                 // do nothing, no pointers
00096         }
00097 };
00098 
00099 
00100 #define GRID_MIXEDCELLMAP_2D(i,j,k,NA,NB) gridMixedCellMap2D[(k)+(MAX_MIXEDCELLS_2D-1)*((j)+(NB)*(i))]
00101 // fortran like
00102 //#define GRIDMIXEDCELLMAP(i,j,k,l,NX,NY,NZ) gridMixedCellMap[(i)+(NX)*((j)+(NY)*((k)+(NZ)*(l)))]
00103 // c like (NX is ignored)
00104 #define GRIDMIXEDCELLMAP(i,j,k,l,NX,NY,NZ) gridMixedCellMap[(l)+(MAX_MIXEDCELLS-1)*((k)+(NZ)*((j)+(NY)*(i)))]
00105 
00106 #ifdef ENABLE_CGAL 
00107 
00108 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
00109 #include <CGAL/Regular_triangulation_3.h>
00110 #include <CGAL/Regular_triangulation_euclidean_traits_3.h>
00111 #include <CGAL/Regular_triangulation_filtered_traits_3.h>
00112 #include <CGAL/Triangulation_cell_base_with_info_3.h>
00113 #include <CGAL/Regular_triangulation_cell_base_3.h>
00114 #include <CGAL/Triangulation_vertex_base_with_info_3.h>
00115 #include <CGAL/Triangulation_data_structure_3.h>
00116 #include <CGAL/Polyhedron_3.h>
00117 #include <CGAL/convex_hull_3.h>
00118 #include <CGAL/Vector_3.h>
00119 #include <CGAL/Plane_3.h>
00120 #include <CGAL/Ray_3.h>
00121 #include <CGAL/Line_3.h>
00122 #include <CGAL/intersections.h>
00123 #include <CGAL/Linear_algebraCd.h>
00124 #include <cassert>
00125 #include <vector>
00126 #include <fstream>
00127 
00129 #endif
00130 
00131 #define DELAUNAY_POINT_CELL             0
00132 #define DELAUNAY_EDGE_CELL              1
00133 #define DELAUNAY_FACET_CELL             2
00134 #define DELAUNAY_TETRA_CELL             3
00135 
00144 class SkinSurface: public Surface
00145 {
00146 
00147 #ifdef ENABLE_CGAL 
00148 private:
00149         
00150         typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
00151         typedef CGAL::Regular_triangulation_filtered_traits_3<K>    Traits;
00152 
00153         typedef CGAL::Triangulation_cell_base_with_info_3< Del3Cell*,Traits> Cell_with_info;
00154         typedef CGAL::Regular_triangulation_cell_base_3<Traits,Cell_with_info> RT_Cell_with_info;
00155         typedef CGAL::Triangulation_vertex_base_3<Traits> Vertex;
00156         typedef CGAL::Triangulation_data_structure_3<Vertex,RT_Cell_with_info> tds;
00157         typedef Traits::RT                                          Weight;
00158         typedef Traits::Bare_point                                  Point3;
00159         typedef Traits::Weighted_point                              Weighted_point;
00160         typedef CGAL::Polyhedron_3<K>                                           Polyhedron;
00161         typedef CGAL::Regular_triangulation_3<Traits,tds>           Rt;
00162         typedef Rt::Vertex_iterator                                 Vertex_iterator;
00163         typedef Rt::Finite_vertices_iterator                        Finite_Vertex_Iterator;
00164         typedef Rt::Finite_cells_iterator                                                       Finite_Cells_Iterator;
00165         typedef Rt::Finite_edges_iterator                                                       Finite_Edges_Iterator;
00166         typedef Rt::Finite_facets_iterator                                                      Finite_Facets_Iterator;
00167         typedef Rt::Vertex_handle                                   Vertex_handle;
00168         typedef tds::Cell_circulator                                                            Cell_circulator;
00169         typedef tds::Cell_handle                                                                        Cell_handle;
00170         typedef Polyhedron::Halfedge_around_facet_circulator        HF_circulator;
00171         typedef Rt::Facet                                                                                       Facet;
00172         typedef CGAL::Linear_algebraCd<double>                                          LA;
00173 
00174 #endif
00175 
00176 private:
00178         int numMixedCells;
00184         int type[4];
00186         int* gridMixedCellMap;
00188         int nx,ny,nz;
00189         double scale;
00190         double side;
00191         double*x,*y,*z;
00192         double s;
00193         double xmin,xmax,ymin,ymax,zmin,zmax;
00194         short*** ind;
00195         int* gridMixedCellMap2D;
00196 
00198         int nx_2d,ny_2d,nz_2d;
00199         double scale_2d;
00200         double side_2d;
00201         double xmin_2d,xmax_2d,ymin_2d,ymax_2d,zmin_2d,zmax_2d;
00202         int** ind_2d;
00203 
00205         MixedCell** mixedComplex;       
00208         bool buildSkin();
00210         bool buildAuxiliaryGrid();
00211         #ifdef ENABLE_CGAL
00212 
00213                 bool buildSkinCGAL();
00214         #endif
00215 
00216         bool savePovRay;
00217 
00218 public:
00220         SkinSurface();
00222         SkinSurface(DelPhiShared* ds);                  
00223 
00225 
00226         virtual bool build();
00228         virtual bool save(char* fileName);
00230         virtual bool load(char* fileName);
00232         virtual void printSummary();            
00234         virtual bool getProjection(double p[3],double* proj1,double* proj2,
00235                 double* proj3,double* normal1,double* normal2,double* normal3);
00241         virtual void getRayIntersection(double p1[3],double p2[3],vector<pair<double,double*> >& intersections,int thdID);              
00243         virtual void preProcessPanel();
00245 
00246         void setShrinking(double ss) 
00247         { 
00248                 if (ss<=1.0 && ss>0.0) 
00249                 { 
00250                         s=ss; 
00251                 }
00252                 else
00253                 {
00254                         cout << endl << WARN << "Cannot set " << ss << ". s parameter is in (0,1]. Setting "<< DEFAULT_S;
00255                         s = DEFAULT_S;
00256                 }
00257         }
00258         
00259         double getShrinking() 
00260         {
00261                 return s;
00262         }
00263 
00264         void setSavePovRay(bool ss) 
00265         { 
00266                 savePovRay = ss;
00267         }
00268         
00269         bool getSavePovRay() 
00270         {
00271                 return savePovRay;
00272         }
00273         
00274         virtual ~SkinSurface();
00275 
00276 private:
00279         //bool rayQuadricIntersection(double*,double*,double**,double*,double*,int thdID,double* cache);
00280         bool rayQuadricIntersection(double*,double*,double*,double*,double*,int thdID,double* cache);
00282         bool isFeasible(MixedCell* mc,double* point);
00284         void projectToQuadric(double* y,double* Q,int type,double* proj,double* norm,double& dist);
00285 };
00286 
00287 #endif