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_fast.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;
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