face.cc

Go to the documentation of this file.
00001 #include "face.h"
00002 
00003 #include <cmath>
00004 #include <iostream>
00005 
00006 #include "box.h"
00007 #include "misc.h"
00008 #include "object.h"
00009 #include "vertex.h"
00010 
00011 using std::cout;
00012 using std::endl;
00013 
00014 Face::Face (const Face& rhs)
00015   :index(rhs.index),bb(rhs.bb),e(rhs.e),v(rhs.v)
00016 {
00017 }
00018 
00019 Face& Face::operator = (const Face& rhs)
00020 {
00021   cout << "Copy assignment operator prohibited on instances of Face class.\n";
00022   cout << "Face " << rhs.index << endl;
00023   exit(0);
00024 }
00025 
00026 Face::Face (char *triplet,vec_vp &vp)
00027   :index(0),bb(6,0.0),e(3,static_cast<Edge*>(NULL)),v(3,static_cast<Vertex*>(NULL))
00028 {
00029   char val[80];
00030   char *eptr;
00031   int i;
00032 
00033   // get past 'Face'
00034   while (strchr("Face",*triplet)!=NULL) {triplet++;}
00035 
00036   // grab Face index
00037   while (strchr(" \t,",*triplet)!=NULL) {triplet++;}
00038   i=0;
00039   while (strchr("0123456789+-eE.",*triplet)!=NULL)
00040   {
00041     val[i++] = *triplet++;
00042   }
00043   val[i]=0;
00044   index = (int) strtod(val,&eptr);
00045   if (val==eptr)
00046   {
00047     index=0;
00048     v[0]=v[1]=v[2]=0;
00049     printf("Error in reading face index\n");
00050     return;
00051   }
00052 
00053   // grab first vertex index
00054   while (strchr(" \t,",*triplet)!=NULL) {triplet++;}
00055   i=0;
00056   while (strchr("0123456789+-eE.",*triplet)!=NULL)
00057   {
00058     val[i++] = *triplet++;
00059   }
00060   val[i]=0;
00061   v[0] = vp[(int)strtod(val,&eptr)-1];
00062   if (val==eptr)
00063   {
00064     v[0]=v[1]=v[2]=0;
00065     printf("Error in reading vertex index\n");
00066     return;
00067   }
00068 
00069   // grab second vertex index
00070   while (strchr(" \t,",*triplet)!=NULL) {triplet++;}
00071   i=0;
00072   while (strchr("0123456789+-eE.",*triplet))
00073   {
00074     val[i++] = *triplet++;
00075   }
00076   val[i]=0;
00077   v[1] = vp[(int)strtod(val,&eptr)-1];
00078   if (val==eptr)
00079   {
00080     v[0]=v[1]=v[2]=0;
00081     printf("Error in reading vertex index\n");
00082     return;
00083   }
00084 
00085   // grab third vertex index
00086   while (strchr(" \t,",*triplet)!=NULL) {triplet++;}
00087   i=0;
00088   while (strchr("0123456789+-eE.",*triplet))
00089   {
00090     val[i++] = *triplet++;
00091   }
00092   val[i]=0;
00093   v[2] = vp[(int)strtod(val,&eptr)-1];
00094   if (val==eptr)
00095   {
00096     v[0]=v[1]=v[2]=0;
00097     printf("Error in reading vertex index\n");
00098     return;
00099   }
00100 
00101   // update face bounding box
00102   setBoundingBox();
00103 }
00104 
00105 bool Face::isMatch (int i, str line) const
00106 {
00107   return i==index && line==v[0]->getObject()->getName();
00108 }
00109 
00110 void Face::print (std::ostream &target) const
00111 {
00112   char str[1024];
00113   sprintf(str,"Face <obj> %s <ind> %d\n",v[0]->getObject()->getName().c_str(),index);
00114   target << str;
00115   v[0]->print(target);
00116   v[1]->print(target);
00117   v[2]->print(target);
00118 }
00119 
00120 void Face::printCP (std::ostream &target) const
00121 {
00122   v[0]->printCP(target);
00123   v[1]->printCP(target);
00124   v[2]->printCP(target);
00125 }
00126 
00127 void Face::getVertexCoord (vec_d const * & v0,vec_d const * & v1,vec_d const * & v2) const
00128 {
00129   v0 = v[0]->getPos();
00130   v1 = v[1]->getPos();
00131   v2 = v[2]->getPos();
00132 }
00133 
00134 double Face::getAngle (Vertex const * const vv) const
00135 {
00136   // identify face vertices
00137   Vertex *vA=const_cast<Vertex*>(vv),*vB=NULL,*vC=NULL;
00138   if (v[0]!=vv){vB=v[0];}
00139   else if (v[1]!=vv){vB=v[1];}
00140   else if (v[2]!=vv){vB=v[2];}
00141   if (v[0]!=vv && v[0]!=vB){vC=v[0];}
00142   else if (v[1]!=vv && v[1]!=vB){vC=v[1];}
00143   else if (v[2]!=vv && v[2]!=vB){vC=v[2];}
00144   // AB,AC
00145   double AB[3] = {vB->getCoord(0)-vA->getCoord(0),
00146                   vB->getCoord(1)-vA->getCoord(1),
00147                   vB->getCoord(2)-vA->getCoord(2)};
00148   double AC[3] = {vC->getCoord(0)-vA->getCoord(0),
00149                   vC->getCoord(1)-vA->getCoord(1),
00150                   vC->getCoord(2)-vA->getCoord(2)};
00151   // lengths
00152   double acL=sqrt( dot(AC,AC) );
00153   double abL=sqrt( dot(AB,AB) );
00154   double costheta=(dot(AB,AC))/abL/acL;
00155   if      (costheta >  1) costheta=1;
00156   else if (costheta < -1) costheta=-1;
00157   return acos(costheta);
00158 }
00159 
00160 double Face::getAngleProxy (Vertex const * const vv) const
00161 {
00162   // identify face vertices
00163   Vertex *vA=const_cast<Vertex*>(vv),*vB=NULL,*vC=NULL;
00164   if (v[0]!=vv){vB=v[0];}
00165   else if (v[1]!=vv){vB=v[1];}
00166   else if (v[2]!=vv){vB=v[2];}
00167   if (v[0]!=vv && v[0]!=vB){vC=v[0];}
00168   else if (v[1]!=vv && v[1]!=vB){vC=v[1];}
00169   else if (v[2]!=vv && v[2]!=vB){vC=v[2];}
00170   // AB,AC
00171   double AB[3] = {vB->getCoord(0)-vA->getCoord(0),
00172                   vB->getCoord(1)-vA->getCoord(1),
00173                   vB->getCoord(2)-vA->getCoord(2)};
00174   double AC[3] = {vC->getCoord(0)-vA->getCoord(0),
00175                   vC->getCoord(1)-vA->getCoord(1),
00176                   vC->getCoord(2)-vA->getCoord(2)};
00177   double BC[3] = {vC->getCoord(0)-vB->getCoord(0),
00178                   vC->getCoord(1)-vB->getCoord(1),
00179                   vC->getCoord(2)-vB->getCoord(2)};
00180   // lengths
00181   double acL= dot(AC,AC);
00182   double abL= dot(AB,AB);
00183   double cbL= dot(BC,BC);
00184   return cbL/(acL+abL);
00185 }
00186 
00187 double Face::getAspectRatio (void) const
00188 {
00189   // Make triangle edge vectors
00190   double va[3]={v[1]->getCoord(0)-v[0]->getCoord(0),
00191                 v[1]->getCoord(1)-v[0]->getCoord(1),
00192                 v[1]->getCoord(2)-v[0]->getCoord(2)};
00193   double vb[3]={v[2]->getCoord(0)-v[1]->getCoord(0),
00194                 v[2]->getCoord(1)-v[1]->getCoord(1),
00195                 v[2]->getCoord(2)-v[1]->getCoord(2)};
00196   double vc[3]={v[0]->getCoord(0)-v[2]->getCoord(0),
00197                 v[0]->getCoord(1)-v[2]->getCoord(1),
00198                 v[0]->getCoord(2)-v[2]->getCoord(2)};
00199   double vbase[3]={0,0,0};
00200   double vopp[3]={0,0,0};
00201 
00202   // Find length of longest edge
00203   double lmax=-1e30;
00204   double la=sqrt(dot(va,va));
00205   double lb=sqrt(dot(vb,vb));
00206   double lc=sqrt(dot(vc,vc));
00207   if (la>lmax)
00208   {
00209     lmax=la;
00210     vbase[0]=va[0];
00211     vbase[1]=va[1];
00212     vbase[2]=va[2];
00213     vc[0]=v[2]->getCoord(0)-v[0]->getCoord(0);
00214     vc[1]=v[2]->getCoord(1)-v[0]->getCoord(1);
00215     vc[2]=v[2]->getCoord(2)-v[0]->getCoord(2);
00216     vopp[0]=vc[0];
00217     vopp[1]=vc[1];
00218     vopp[2]=vc[2];
00219   }
00220   if (lb>lmax)
00221   {
00222     lmax=lb;
00223     vbase[0]=vb[0];
00224     vbase[1]=vb[1];
00225     vbase[2]=vb[2];
00226     va[0]=v[0]->getCoord(0)-v[1]->getCoord(0);
00227     va[1]=v[0]->getCoord(1)-v[1]->getCoord(1);
00228     va[2]=v[0]->getCoord(2)-v[1]->getCoord(2);
00229     vopp[0]=va[0];
00230     vopp[1]=va[1];
00231     vopp[2]=va[2];
00232   }
00233   if (lc>lmax)
00234   {
00235     lmax=lc;
00236     vbase[0]=vc[0];
00237     vbase[1]=vc[1];
00238     vbase[2]=vc[2];
00239     vb[0]=v[1]->getCoord(0)-v[2]->getCoord(0);
00240     vb[1]=v[1]->getCoord(1)-v[2]->getCoord(1);
00241     vb[2]=v[1]->getCoord(2)-v[2]->getCoord(2);
00242     vopp[0]=vb[0];
00243     vopp[1]=vb[1];
00244     vopp[2]=vb[2];
00245   }
00246 
00247   // Find shortest altitude
00248   double ll = 1.0/sqrt(dot(vbase,vbase));
00249   vbase[0]=vbase[0]*ll;
00250   vbase[1]=vbase[1]*ll;
00251   vbase[2]=vbase[2]*ll;
00252   double dot_prod = dot(vbase,vopp);
00253   double alt[3]={vopp[0]-(dot_prod*vbase[0]),
00254                  vopp[1]-(dot_prod*vbase[1]),
00255                  vopp[2]-(dot_prod*vbase[2])};
00256   double amin=sqrt(dot(alt,alt));
00257   return lmax/amin;
00258 }
00259 
00260 void Face::addEdge (Edge * ptr)
00261 {
00262   if     (e[0]==NULL){e[0]=ptr;}
00263   else if(e[1]==NULL){e[1]=ptr;}
00264   else if(e[2]==NULL){e[2]=ptr;}
00265   else
00266   {
00267     cout << "Error. Tried to add fourth edge to face.\n"
00268           << "Face " << index 
00269           << " " << (int)v[0]->getIndex()
00270           << " " << (int)v[1]->getIndex()
00271           << " " << (int)v[2]->getIndex()
00272           << endl;
00273     exit(1); 
00274   }
00275 }
00276 
00277 vec_d Face::getNormal (void) const
00278 {
00279   double uX, uY, uZ, vX, vY, vZ;
00280   // compute vectors 01 and 12
00281   uX = v[1]->getCoord(0)-v[0]->getCoord(0);
00282   uY = v[1]->getCoord(1)-v[0]->getCoord(1);
00283   uZ = v[1]->getCoord(2)-v[0]->getCoord(2);
00284   vX = v[2]->getCoord(0)-v[0]->getCoord(0);
00285   vY = v[2]->getCoord(1)-v[0]->getCoord(1);
00286   vZ = v[2]->getCoord(2)-v[0]->getCoord(2);
00287   // compute cross product (u x v)
00288   vec_d n(3,0);
00289   n[0] = uY*vZ-uZ*vY;
00290   n[1] = uZ*vX-uX*vZ;
00291   n[2] = uX*vY-uY*vX;
00292   return n;
00293 }
00294 
00295 void Face::updateBoxes (mmap_fp_bp &ob,mmap_fp_bp &nb)
00296 {
00297   // this = adjacent face of current vertex
00298   // ob = multimap of current vertex adjacent faces
00299   //      and the boxes in which they perviously lay
00300   // nb = multimap of current vertex adjacent faces
00301   //      and the boxes in which they now lie
00302 
00304   bp_list ovec;
00305   //ovec.reserve(VECTOR_RESERVE);
00306   std::pair<tf_it,tf_it> j=ob.equal_range(this);
00307   // for each box* in old list (ob)
00308   for(tf_it i=j.first;i!=j.second;i++)
00309   {
00310     ovec.push_back((*i).second);
00311   }
00312 
00314   j=nb.equal_range(this);
00315   // for each box* in new list (nb)
00316   tf_it i=j.first;
00317   while(i!=j.second)
00318   {
00319     // if new box* found in old list (ovec)
00320     bpl_it q = find(ovec.begin(),ovec.end(),(*i).second);
00321     // if found, then remove box* from old and new list
00322     if(q!=ovec.end())
00323     {
00324       q=ovec.erase(q);
00325       tf_it k = i;k++;
00326       nb.erase(i);
00327       i=k;
00328     }
00329     else {i++;}
00330   }
00331 
00333   if (!ovec.empty())
00334   {
00335     // for each remaining Box* in ovec (the old list)
00336     for (bpl_it k=ovec.begin();k!=ovec.end();k++)
00337     {
00338       // sort Face* vector
00339       (*k)->sortFaces();
00340       // look for this face in Box* face list
00341       // and remove if found
00342       (*k)->removeFace(this);
00343 //      //sort((*j)->f.begin(),(*j)->f.end());
00344 //      // look for this face in Box* face list
00345 //      std::pair<fp_it,fp_it> q;
00346 //      q=equal_range((*j)->f.begin(),(*j)->f.end(),this);
00347 //      // if found, then remove
00348 //      if(q.first!=q.second)
00349 //      {
00350 //        (*j)->f.erase(q.first);
00351 //      }
00352     }
00353   }
00354 
00356   j=nb.equal_range(this);
00357   // for each box* in new list (nb)
00358   for(i=j.first;i!=j.second;i++)
00359   {
00360     // add this face to Box* face vector
00361     (*i).second->addFace(this);
00362     //(*i).second->f.push_back(this);
00363   }
00364 }
00365 
00366 void Face::setBoundingBox (void)
00367 {
00368   // assume bb=bounding box = [xmn xmax ymin ymax zmin zmax]
00369   // for each vertex of face
00370   for(int j=0;j<3;j++)
00371   {
00372     if ( v[0]->getCoord(j) < v[1]->getCoord(j) )
00373     {
00374       if ( v[0]->getCoord(j) < v[2]->getCoord(j) )
00375       {
00376         if (v[1]->getCoord(j) < v[2]->getCoord(j))
00377         {
00378           bb[2*j]=v[0]->getCoord(j);
00379           bb[2*j+1]=v[2]->getCoord(j);
00380         }
00381         else
00382         {
00383           bb[2*j]=v[0]->getCoord(j);
00384           bb[2*j+1]=v[1]->getCoord(j);
00385         }
00386       }
00387       else
00388       {
00389         bb[2*j]=v[2]->getCoord(j);
00390         bb[2*j+1]=v[1]->getCoord(j);
00391       }
00392     }
00393     else
00394     {
00395       if ( v[0]->getCoord(j) < v[2]->getCoord(j) )
00396       {
00397         bb[2*j]=v[1]->getCoord(j);
00398         bb[2*j+1]=v[2]->getCoord(j);
00399       }
00400       else
00401       {
00402         if (v[1]->getCoord(j) < v[2]->getCoord(j))
00403         {
00404           bb[2*j]=v[1]->getCoord(j);
00405           bb[2*j+1]=v[0]->getCoord(j);
00406         }
00407         else
00408         {
00409           bb[2*j]=v[2]->getCoord(j);
00410           bb[2*j+1]=v[0]->getCoord(j);
00411         }
00412       }
00413     }
00414   }
00415 }
00416 
00417 double Face::getBoundingBox (int i) const
00418 {
00419   return bb[i];
00420 }
00421 
00422 Edge * Face::getEdge (int i)
00423 {
00424   return e[i];
00425 }
00426 
00427 Vertex * Face::getVertex (int i) const
00428 {
00429   return v[i];
00430 }
00431 
00432 int Face::getIndex (void) const
00433 {
00434   return index;
00435 }

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