misc.cc

Go to the documentation of this file.
00001 #include "misc.h"
00002 
00003 #include <cmath>
00004 #include <fstream>
00005 #include <iostream>
00006 
00007 #include "face.h"
00008 #include "vertex.h"
00009 
00010 using std::cout;
00011 using std::endl;
00012 
00013 str processDir (str line)
00014 {
00015   if (line.find_last_of('/')==str::npos)
00016   {
00017     line.push_back('/');
00018     return line;
00019   }
00020   size_t pos = line.find_last_of('/');
00021   if ((pos+1)==line.size())
00022   {
00023     return line;
00024   }
00025   else
00026   {
00027     line.push_back('/');
00028     return line;
00029   }
00030 }
00031 
00032 int distinguishable (double a,double b,double epsilon)
00033 {
00034   double c;
00035   c=a-b;
00036   if (c<0) c=-c;
00037   if (a<0) a=-a;
00038   if (a<1) a=1;
00039   if (b<0) b=-b;
00040   if (b<a) return (c>a*epsilon);
00041   else return (c>b*epsilon);
00042 }
00043 
00044 int distinguishable (double a,double b)
00045 {
00046   double c;
00047   c=a-b;
00048   if (c<0) c=-c;
00049   if (a<0) a=-a;
00050   if (a<1) a=1;
00051   if (b<0) b=-b;
00052   if (b<a) return (c>a*DOUBLE_EPSILON);
00053   else return (c>b*DOUBLE_EPSILON);
00054 }
00055 
00056 double dot (vec_d const & a,vec_d const & b)
00057 {
00058   return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
00059 }
00060 
00061 double dot (double const * const a,double const * const b)
00062 {
00063   return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
00064 }
00065 
00066 bool parseCommandLine (int argc,char **argv,str message)
00067 {
00068   bool freeze=false;
00069   int c;
00070   opterr=0;
00071   char *eptr=NULL;
00072   str line;
00073   while((c=getopt(argc,argv,"i:o:t:s:a:f:hx:y:z:g:")) != -1)
00074     switch(c)
00075     {
00076       case 'i':
00077         // specify input directory 
00078         line = optarg;
00079         INPUT_DATA_DIR = processDir(line);
00080         break;
00081       case 'o':
00082         // specify input directory 
00083         line = optarg;
00084         OUTPUT_DATA_DIR = processDir(line);
00085         break;
00086       case 't':
00087         // specify target extracellular width 
00088         TARGET_SEPARATION = strtod(optarg,&eptr);
00089         break;
00090       case 'x':
00091         // specify separation gain 
00092         SEPARATION_GAIN     = strtod(optarg,&eptr);
00093         break;
00094       case 'y':
00095         // specify angle gain 
00096         ANGLE_GAIN     = strtod(optarg,&eptr);
00097         break;
00098       case 'z':
00099         // specify stretch gain 
00100         STRETCH_GAIN     = strtod(optarg,&eptr);
00101         break;
00102       case 'g':
00103         // specify number of groups 
00104         NUM_GROUPS     = static_cast<int>(strtod(optarg,&eptr));
00105         break;
00106       case 's':
00107         // specify separation weight 
00108         SEPARATION_WEIGHT     = strtod(optarg,&eptr)/100.0;
00109         break;
00110       case 'a':
00111         // specify angle weight
00112         ANGLE_STRETCH_WEIGHT  = strtod(optarg,&eptr);
00113         break;
00114       case 'f':
00115         // specify frozen vertices file
00116         FROZEN_VERTICES_FILE = optarg;
00117         freeze = true;
00118         break;
00119       case 'h':
00120         cout << message << endl;
00121         abort ();
00122       case '?':
00123         if (optopt == 'c')
00124           fprintf (stderr,"Option -%c requires an argument.\n", optopt);
00125         else if (isprint (optopt))
00126           fprintf (stderr, "Unknown option '-%c.\n", optopt);
00127         else
00128           fprintf (stderr,
00129                    "Unknown option character '\\x%x'.\n",
00130                    optopt);
00131         exit(0);
00132       default:
00133         cout << message << endl;
00134         abort ();
00135     }
00136 
00137   if(optind<argc)
00138   {
00139     for (int index=optind;index<argc;index++)
00140     {
00141       cout << "\nError: extra argument found = " << argv[index] << endl;
00142     }
00143     cout << message << endl;
00144     exit(0);
00145   }
00146   return freeze;
00147 }
00148 
00149 bool checkIntSize (void)
00150 {
00152   if (32==sizeof(int)*8){ return true; }
00153   else
00154   {
00155     cout << "Error. Int is not 32 bits, sizeof(int) "
00156           << sizeof(int) << endl;
00157     return false;
00158   }
00159 }
00160 
00161 void copyControlFile (void)
00162 {
00163   cout << "Copy control settings..........................";
00164   cout.flush();
00165   str line;
00166   // open input data file
00167   str filein("controls.cc");
00168   std::ifstream inFile(filein.c_str());
00169   if (inFile.fail()) // if stream cannot be opened
00170   { 
00171     cout << "\n\ncopyControlFile: Can not open input file ["
00172           << filein << "]\n\n";
00173     exit(1); 
00174   }
00175   // open output data file
00176   str fileout = OUTPUT_DATA_DIR + CONTROL_FILE;
00177   std::ofstream outFile(fileout.c_str());
00178   if (outFile.fail()) // if stream cannot be opened
00179   {
00180     cout << "\n\ncopyControlFile: Can not open output file ["
00181           << fileout << "]\n\n";
00182     exit(1); 
00183   }
00184   // foreach line in input file
00185   while (inFile.eof()==false)
00186   {
00187     outFile << line << endl;
00188   }
00189   // close files
00190   inFile.close();
00191   outFile.close();
00192   cout << "complete.\n";
00193   cout.flush();
00194 }
00195 
00196 //int biggest (double *pn)
00197 int biggest (vec_d const & pn)
00198 {
00199   const double x[3] = {fabs(pn[2]),fabs(pn[0]),fabs(pn[1])};
00200   if ( x[0] < x[1] )
00201   {
00202     // x[0] < x[1]
00203     if ( x[0] < x[2] )
00204     {
00205       // x[0] < x[1]
00206       // x[0] < x[2]
00207       if (x[1] < x[2])
00208       {
00209         // x[0] < x[1]
00210         // x[1] < x[2]
00211         // x[0] is smallest
00212         // x[2] is biggest
00213         return 2;
00214       }
00215       else
00216       {
00217         // x[0] < x[1]
00218         // x[2] <= x[1]
00219         // x[0] is smallest
00220         // x[1] is biggest or tied with x[2]
00221         return 1;
00222       }
00223     }
00224     else
00225     {
00226       // x[0] < x[1]
00227       // x[2] <= x[0]
00228       // x[2] is smallest or tied with x[0]
00229       // x[1] is biggest
00230       return 1;
00231     }
00232   }
00233   else
00234   {
00235     // x[1] <= x[0]
00236     if ( x[0] < x[2] )
00237     {
00238       // x[1] <= x[0]
00239       // x[0] < x[2]
00240       // x1 is smallest or tied with x[0]
00241       // x2 is biggest
00242       return 2;
00243     }
00244     else
00245     {
00246       return 0;
00247     }
00248   }
00249 }
00250 
00251 //void checkLineFaceInt(Face const * const f,vec_d const & origin,vec_d const & end,bool &line_flag,
00252 //                               bool &poly_flag, bool &poly_edge_flag,bool ignore_line_flag)
00253 //{
00254 //  // Use line points, polygon normal and one polygon vertex
00255 //  // to compute where line intersects plane of polygon (if not parallel).
00256 //  //
00257 //  // pI = p1+u*(p2-p1) where u=dot(N,(p3-p1))/dot(N,(p2-p1))
00258 //  // pI = intersection of line and plane
00259 //  // p3 = any point on plane
00260 //  // p2,p1 = end points of line
00261 //  // N = plane normal
00262 //  // lp  = line_points
00263 //
00264 //  // initialize flags
00265 //  line_flag=poly_flag=poly_edge_flag=false;
00266 //  // get face normal
00267 //  double pn[3];
00268 //  f->getNormal(pn);
00269 //  // denominator of u
00270 //  const double den = pn[0]*(end[0]-origin[0])+pn[1]*(end[1]-origin[1])+pn[2]*(end[2]-origin[2]);
00271 //  // if line and polygon plane are not parallel
00272 //  if(den)
00273 //  {
00274 //    //pvc = polygon_vertex_coordinates
00275 //    Vertex * const v0=f->getVertex(0),* const v1=f->getVertex(1),* const v2=f->getVertex(2); 
00276 //    // point of intersection
00277 //    const double u = (pn[0]*(v0->getCoord(0)-origin[0]) 
00278 //                      + pn[1]*(v0->getCoord(1)-origin[1]) 
00279 //                      + pn[2]*(v0->getCoord(2)-origin[2]));
00280 //    // if polygon cuts through line
00281 //    if(den>0)
00282 //    {
00283 //      if((u>0.0) && (u<den)){ line_flag = true;}
00284 //    }
00285 //    else
00286 //    {
00287 //      if((u<0.0) && (u>den)){ line_flag = true;}
00288 //    }
00289 //    // check early exit
00290 //    if(line_flag==false && ignore_line_flag==false){return;}
00291 //    // compute polygon double-area on each of three principal planes
00292 //    // area = {xy,yz,zx}
00293 //    const int big = biggest(pn);
00294 //    // pI = p1+u*(p2-p1)
00295 //    // where end[] is p2 and origin[] is p1
00296 //    const int t = (big+1)%3;
00297 //    const int tv[2] = {big,t};
00298 //    const double pI[2] = { (den-u)*origin[big] + u*end[big],
00299 //      (den-u)*origin[t]   + u*end[t]};
00300 //    ////////// is point of intersection on other polygon? //////////
00301 //    //pvc = polygon_vertex_coordinates
00302 //    // compute three determinants
00303 //    const double a  = den*v0->getCoord(tv[0])-pI[0];
00304 //    const double b  = den*v0->getCoord(tv[1])-pI[1];
00305 //    const double c  = den*v1->getCoord(tv[0])-pI[0];
00306 //    const double d  = den*v1->getCoord(tv[1])-pI[1];
00307 //    const double e  = den*v2->getCoord(tv[0])-pI[0];
00308 //    const double ff = den*v2->getCoord(tv[1])-pI[1];
00309 //    const double det = a*d-c*b ;
00310 //    const double val1=det*(c*ff-e*d);
00311 //    const double val2=det*(e*b-a*ff);
00312 //    // proceed if determinants are DOUBLE_EPSILON away from zero
00313 //    if (val1>0)
00314 //    {
00315 //      if (val2>0)
00316 //      {
00317 //        // line intersects polygon plane inside polygon
00318 //        poly_flag = true;
00319 //      }
00320 //    }
00321 //    if(fabs(val1)<DOUBLE_EPSILON*den*den*den*den && fabs(val2)<DOUBLE_EPSILON*den*den*den*den)
00322 //    {
00323 //      // does point of intersection lie on polygon edge?
00324 //      poly_edge_flag = true;
00325 //    }
00326 //  }
00327 //}
00328 
00329 result checkLineFaceInt (Face const * const f,vec_d const & origin,
00330                           vec_d const & end,bool ignore_line_flag)
00331 {
00332   // Use line points, polygon normal and one polygon vertex
00333   // to compute where line intersects plane of polygon (if not parallel).
00334   //
00335   // pI = p1+u*(p2-p1) where u=dot(N,(p3-p1))/dot(N,(p2-p1))
00336   // pI = intersection of line and plane
00337   // p3 = any point on plane
00338   // p2,p1 = end points of line
00339   // N = plane normal
00340   // lp  = line_points
00341 
00342   result r;
00343   // initialize flags
00344   r.line_flag=r.poly_flag=r.poly_edge_flag=false;
00345   // get face normal
00346   //double pn[3];
00347   vec_d pn = f->getNormal();
00348   // denominator of u
00349   const double den = pn[0]*(end[0]-origin[0])+pn[1]*(end[1]-origin[1])+pn[2]*(end[2]-origin[2]);
00350   // if line and polygon plane are not parallel
00351   if(den)
00352   {
00353     //pvc = polygon_vertex_coordinates
00354     Vertex * const v0=f->getVertex(0),* const v1=f->getVertex(1),* const v2=f->getVertex(2); 
00355     // point of intersection
00356     const double u = (pn[0]*(v0->getCoord(0)-origin[0]) 
00357                       + pn[1]*(v0->getCoord(1)-origin[1]) 
00358                       + pn[2]*(v0->getCoord(2)-origin[2]));
00359     // if polygon cuts through line
00360     if(den>0)
00361     {
00362       if((u>0.0) && (u<den)){ r.line_flag = true;}
00363     }
00364     else
00365     {
00366       if((u<0.0) && (u>den)){ r.line_flag = true;}
00367     }
00368     // check early exit
00369     if(r.line_flag==false && ignore_line_flag==false){return r;}
00370     // compute polygon double-area on each of three principal planes
00371     // area = {xy,yz,zx}
00372     const int big = biggest(pn);
00373     // pI = p1+u*(p2-p1)
00374     // where end[] is p2 and origin[] is p1
00375     const int t = (big+1)%3;
00376     const int tv[2] = {big,t};
00377     const double pI[2] = { (den-u)*origin[big] + u*end[big],
00378       (den-u)*origin[t]   + u*end[t]};
00380     //pvc = polygon_vertex_coordinates
00381     // compute three determinants
00382     const double a  = den*v0->getCoord(tv[0])-pI[0];
00383     const double b  = den*v0->getCoord(tv[1])-pI[1];
00384     const double c  = den*v1->getCoord(tv[0])-pI[0];
00385     const double d  = den*v1->getCoord(tv[1])-pI[1];
00386     const double e  = den*v2->getCoord(tv[0])-pI[0];
00387     const double ff = den*v2->getCoord(tv[1])-pI[1];
00388     const double det = a*d-c*b ;
00389     const double val1=det*(c*ff-e*d);
00390     const double val2=det*(e*b-a*ff);
00391     // proceed if determinants are DOUBLE_EPSILON away from zero
00392     if (val1>0)
00393     {
00394       if (val2>0)
00395       {
00396         // line intersects polygon plane inside polygon
00397         r.poly_flag = true;
00398       }
00399     }
00400     if(fabs(val1)<DOUBLE_EPSILON*den*den*den*den && fabs(val2)<DOUBLE_EPSILON*den*den*den*den)
00401     {
00402       // does point of intersection lie on polygon edge?
00403       r.poly_edge_flag = true;
00404     }
00405   }
00406   return r;
00407 }
00408 
00409 double l2 (double * const t)
00410 {
00411   // return length of t squared
00412   return t[0]*t[0]+t[1]*t[1]+t[2]*t[2];
00413 }
00414 
00415 void cross (double * const t1,double * const t2,double * const cr)
00416 {
00417   // normals cross product = t1 x t2
00418   cr[0] = t1[1]*t2[2]-t1[2]*t2[1];
00419   cr[1] = t1[2]*t2[0]-t1[0]*t2[2];
00420   cr[2] = t1[0]*t2[1]-t1[1]*t2[0];
00421 }

Generated on Fri Jul 18 19:43:40 2008 for meshmorph by  doxygen 1.5.1