vertex.cc

Go to the documentation of this file.
00001 #include "vertex.h"
00002 
00003 #include <cassert>
00004 #include <cmath>
00005 #include <iostream>
00006 
00007 #include "container.h"
00008 #include "edge.h"
00009 #include "face.h"
00010 #include "intersecting_faces.h"
00011 #include "misc.h"
00012 #include "nice.h"
00013 #include "object.h"
00014 
00015 using std::cout;
00016 using std::endl;
00017 
00018 Vertex::Vertex (const Vertex& rhs)
00019   :index(rhs.index),p(rhs.p),cl(rhs.cl),o(rhs.o),
00020   f(rhs.f),energy(rhs.energy)
00021 {
00022 }
00023 
00024 Vertex& Vertex::operator = (const Vertex& rhs)
00025 {
00026   cout << "Assignment operator prohibited on instances of vertex class.\n";
00027   cout << "Vertex " << rhs.index << endl;
00028   exit(0);
00029 }
00030 
00037 Vertex::Vertex (char const * triplet,Object * const q)
00038   :index(0),p(),cl(NULL),o(q),f(),energy(0)
00039 {
00040   int i;
00041   char val[80];
00042   char * eptr;
00043   char const * cp = triplet;
00044 
00045   // get past 'Vertex'
00046   while (strchr("Vertx",*triplet)!=NULL) {triplet++;}
00047 
00048   // grab vertex index
00049   while (strchr(" \t,",*triplet)!=NULL) {triplet++;}
00050   i=0;
00051   while (strchr("0123456789+-eE.",*triplet)!=NULL)
00052   {
00053     val[i++] = *triplet++;
00054   }
00055   val[i]=0;
00056   index = (int) strtod(val,&eptr);
00057   if (val==eptr)
00058   {
00059     index=0;
00060     printf("Error in reading vertex index\n");
00061     return;
00062   }
00063 
00064   // grab x coord
00065   while (strchr(" \t,",*triplet)!=NULL) {triplet++;}
00066   i=0;
00067   while (strchr("0123456789+-eE.",*triplet)!=NULL)
00068   {
00069     val[i++] = *triplet++;
00070   }
00071   val[i]=0;
00072   p.push_back(strtod(val,&eptr));
00073   if (val==eptr)
00074   {
00075     printf("Error in reading vertex\n");
00076     printf("Error in reading vertex: string %s\n",cp);
00077     return;
00078   }
00079 
00080   // grab y coord
00081   while (strchr(" \t,",*triplet)!=NULL) triplet++;
00082   i=0;
00083   while (strchr("0123456789+-eE.",*triplet))
00084   {
00085     val[i++] = *triplet++;
00086   }
00087   val[i]=0;
00088   p.push_back(strtod(val,&eptr));
00089   if (val==eptr)
00090   {
00091     printf("Error in reading vertex\n");
00092     printf("Error in reading vertex: string %s\n",cp);
00093     return;
00094   }
00095 
00096   // grab z coord
00097   while (strchr(" \t,",*triplet)!=NULL) triplet++;
00098   i=0;
00099   while (strchr("0123456789+-eE.",*triplet))
00100   {
00101     val[i++] = *triplet++;
00102   }
00103   val[i]=0;
00104   p.push_back(strtod(val,&eptr));
00105   if (val==eptr)
00106   {
00107     printf("Error in reading vertex\n");
00108     printf("Error in reading vertex: string %s\n",cp);
00109     return;
00110   }
00111 }
00112 
00118 int Vertex::getIndex (void) const
00119 {
00120   return index;
00121 }
00122 
00128 void Vertex::print (std::ostream & target) const
00129 {
00130   char str[1024];
00131   sprintf(str,"Vertex <obj> %s <ind> %d",o->getName().c_str(),index);
00132   sprintf(str,"%s [%.15g %.15g %.15g]\n",str,p[0],p[1],p[2]);
00133   target << str;
00134 }
00135 
00142 void Vertex::printCP (std::ostream &target) const
00143 {
00144   char str[1024];
00145   sprintf(str,"%.15g %.15g %.15g 1 0 0 1\n",p[0],p[1],p[2]);
00146   target << str;
00147 }
00148 
00157 bool Vertex::isMatch (int const i, str const name) const
00158 {
00159   return i==index && name==o->getName();
00160 }
00168 vec_d Vertex::getNormal (void) const
00169 {
00170   double theta,thetaT=0.0,L;
00171   //n[0]=n[1]=n[2]=0;
00172   vec_d n(3,.00);
00173   // for each adjacent face
00174   for (fp_cit i=f.begin();i!=f.end();i++)
00175   {
00176     // get coordinates of polygon normal
00177     vec_d t = (*i)->getNormal();
00178     L=sqrt(t[0]*t[0]+t[1]*t[1]+t[2]*t[2]);
00179     //theta=(*i)->getAngle(this);
00180     theta=(*i)->getAngleProxy(this);
00181     thetaT+=theta;
00182     // and add to sum
00183     double a = theta/L;
00184     n[0] += t[0]*a;
00185     n[1] += t[1]*a;
00186     n[2] += t[2]*a;
00187   }
00188   double a = 1.0/f.size()/thetaT;
00189   n[0] = n[0]*a;
00190   n[1] = n[1]*a;
00191   n[2] = n[2]*a;
00192   return n;
00193 }
00194 
00204 double Vertex::getSeparationForceEnergy (double * const &force,bool compute_force) const
00205 {
00206   if(cl!=NULL)
00207   {
00208     // get closest point
00209     //double pC[3];
00210     //vec_d pC = Container::instance().getNearPtOnFaceToVertex(cl,this);
00211     vec_fp fp(1,cl);
00212     vec_d dummy_p;
00213     double sqd=0.0;
00214     Face *dummy_ncl = NULL;
00215     Container::instance().findClosestPtToVertexAmongFaces(this,fp,dummy_p,sqd,dummy_ncl);
00216     // DEBUG
00217     //assert(sqd<1E10);
00218     //assert(sqd==sqd);
00219     // DEBUG
00220     // compute separation vector
00221     //double s[3];
00222     vec_d s(3,0);
00223     //for(int i=0;i<3;i++){ s[i]=pC[i]-p[i]; }
00224     // compute separation distance //////////
00225     double sd = sqrt(sqd);
00226     // compute separation error (signed value)
00227     // NOTE THIS ASSUMES VERTEX IS INSIDE SAME OBJECT
00228     // ON WHICH CLOSEST POINT WAS FOUND, OTHERWISE
00229     // OTHERWISE SE!=sd+TARGET_SEPARATION
00230     // also note that se=sd+TARGET_SEPARATION is correct
00231     // for nonnice vertex, since pC-p is oriented so as
00232     // to move vertex out of violated object
00233 
00234     // set target separation
00235     double TS = TARGET_SEPARATION*SCALE;
00236     // if closest face is inside vertex neighborhood
00237     //if(Container::instance().faceInNeighborhood(cl,this)==true){ TS = LOOP_TARGET_SEPARATION*SCALE; }
00238     //else { TS = TARGET_SEPARATION*SCALE; }
00239 
00240     double se;
00241     if(Nice::instance().vertexIsNice(this)==false)
00242     {
00243       se=sd+TS;
00244     }
00245     else
00246     {
00247       se=sd-TS;
00248     }
00249     if(TARGET_SEPARATION==0)
00250     {
00251       cout << "Vertex::getSeparationForceEnergy: we got problems."
00252             << "TARGET_SEPARATION==0.\n";
00253       exit(0);
00254     }
00255     se = se/TARGET_SEPARATION;
00256     // if pC==p, i.e. curretn vertex is on surface of neighbor object
00257     // then use current vertex outward normal as separation vector
00258     // recompute separation distance since it is used as normalization factor
00259     if (!sd)
00260     {
00261       s = getNormal();
00262       sd = sqrt(dot(s,s));
00263     }
00264     if(compute_force==true)
00265     {
00266       // spring_force = spring constant * stretch
00267       // let scaled_spring force = spring_force/separation_distance
00268       // Note SEPARATION_WEIGHT has already been normalized to 0 to 1.
00269       double scaled_spring_force = SEPARATION_GAIN*SEPARATION_WEIGHT*se/sd;
00270       // force cartesian component = spring_force * unit vector
00271       // where unit vector = cartesian separation_distance component / separation_distance
00272       // force cartesian component = scaled_spring_force * cartesian separation_distance component
00273       for (int i=0;i<3;i++) force[i]+=scaled_spring_force*s[i];
00274     }
00275     // return energy
00276     return SEPARATION_WEIGHT/2.0*se*se;
00277   }
00278   else
00279   {
00280     // cl==NULL, i.e. no closest point
00281     // therefore no contribution to force
00282     // peg energy at maximum possible value, i.e. max separation error
00283     // maximum separation error would occur when closest point found
00284     // at SEARCH_RADIUS distance from current vertex AND vertex is not nice
00285     // so max(se) = SEARCH_RADIUS + TARGET_SEPARATION
00286     double max_se = sqrt(SEARCH_RADIUS_SQ*SCALE*SCALE) + LOOP_TARGET_SEPARATION*SCALE;
00287     return SEPARATION_WEIGHT/2.0*max_se*max_se;
00288   }
00289 }
00290 
00291 double Vertex::getEdgeStretchForceEnergy (double * const & force,bool flag) const
00292 {
00293   double en_ergy=0;
00294   vec_ep e;
00295   e.reserve(VECTOR_RESERVE);
00296   getAdjacentEdges(e);
00297   // for each adjacent edge of current vertex
00298   for (ep_it j=e.begin();j!=e.end();j++)
00299   {
00300     // compute force and energy of edge
00301     en_ergy+=(*j)->getStretchForceEnergy(this,force,flag);
00302   }
00303   return en_ergy;
00304 }
00305 
00306 double Vertex::getEdgeAngleForceEnergy (double * const &force,bool flag) const
00307 {
00308   // NOTE : I previously decided NOT to use nonadjacent edges
00309   // so oe is not being used at this time
00310   Vertex *v1=NULL,*v2=NULL,*o1=NULL,*o2=NULL;
00311   vec_ep ae,oe;
00312   ae.reserve(VECTOR_RESERVE);
00313   oe.reserve(VECTOR_RESERVE);
00314 
00315   // for each adjacent face of current vertex
00316   for(fp_cit i=f.begin();i!=f.end();i++)
00317   {
00318     // for each face edge
00319     for(int j=0;j<3;j++)
00320     {
00321       (*i)->getEdge(j)->getVertices(v1,v2,o1,o2);
00322       // if either edge vertex is current vertex
00323       if(v1==this || v2==this)
00324       {
00325         // add edge to adjacent edge vector, ae
00326         ae.push_back((*i)->getEdge(j));
00327       }
00328       else
00329       {
00330         // add edge to non-adjacent edge vector, oe
00331         oe.push_back((*i)->getEdge(j));
00332       }
00333     }
00334   }
00335   // keep unique edges
00336   sort(ae.begin(),ae.end());
00337   ae.assign(ae.begin(),unique(ae.begin(),ae.end()));
00338   sort(oe.begin(),oe.end());
00339   oe.assign(oe.begin(),unique(oe.begin(),oe.end()));
00340   // for each adjacent edge
00341   double en_ergy=0;
00342   for (ep_it i=ae.begin();i!=ae.end();i++)
00343   {
00344     // force and energy contribution
00345     en_ergy+=(*i)->getAngleReForceEnergy(force,flag);
00346   }
00347   return en_ergy;
00348 }
00349 
00350 void Vertex::getAdjFaceIntForce (double * const & force) const
00351 {
00352   // for each adjacent face of current vertex
00353   for(fp_cit i=f.begin();i!=f.end();i++)
00354   {
00355     // if adjacent face is intersected
00356     if(Intersecting_Faces::instance().faceIsIntersectedLHS(*i)==true)
00357     {
00358       // add intersection_force
00359       Intersecting_Faces::instance().getFaceIntersectionForce(*i,force);
00360     }
00361   }
00362 }
00363 
00364 void Vertex::getTotalForceEnergy (double * const &force) const
00365 {
00367   // true -> compute force
00368   getSeparationForceEnergy(force,true);
00369   // DEBUG
00370   if (false) {cout << Container::instance().o.size() << endl;}
00371   bool print_switch = false;
00372   double fmag = 0.0;
00373   double f_old[3];
00374   f_old[0] = force[0];
00375   f_old[1] = force[1];
00376   f_old[2] = force[2];
00377   if(print_switch==true)
00378   {
00379     cout << "Vertex::getTotalForceEnergy: force ["
00380           << force[0] << " "
00381           << force[1] << " "
00382           << force[2] << "\n";
00383   }
00384   // DEBUG
00385   if(print_switch==true)
00386   {
00387     double diff = sqrt(dot(force,force))-fmag;
00388     fmag = fmag+diff;
00389     f_old[0] = force[0];
00390     f_old[1] = force[1];
00391     f_old[2] = force[2];
00392     cout << "Vertex::getTotalForceEnergy: sep force       ["
00393           << force[0] << " "
00394           << force[1] << " "
00395           << force[2] << "]"
00396           << ", diff = " << diff << endl;
00397   }
00398 
00399   // DEBUG
00401   getEdgeStretchForceEnergy(force,true);
00402   // DEBUG
00403   if(print_switch==true)
00404   {
00405     double diff = sqrt(dot(force,force))-fmag;
00406     fmag = fmag+diff;
00407     cout << "Vertex::getTotalForceEnergy: stretch force   ["
00408           << force[0]-f_old[0] << " "
00409           << force[1]-f_old[1] << " "
00410           << force[2]-f_old[2] << "]"
00411           << ", diff = " << diff << endl;
00412     f_old[0] = force[0];
00413     f_old[1] = force[1];
00414     f_old[2] = force[2];
00415   }
00416   // DEBUG
00418   getEdgeAngleForceEnergy(force,true);
00419   // DEBUG
00420   if(print_switch==true)
00421   {
00422     double diff = sqrt(dot(force,force))-fmag;
00423     fmag = fmag+diff;
00424     cout << "Vertex::getTotalForceEnergy: angle force     ["
00425           << force[0]-f_old[0] << " "
00426           << force[1]-f_old[1] << " "
00427           << force[2]-f_old[2] << "]"
00428           << ", diff = " << diff << endl;
00429     cout << "Vertex::getTotalForceEnergy: resultant force ["
00430           << force[0] << " "
00431           << force[1] << " "
00432           << force[2] << "]\n\n";
00433   }
00434   // DEBUG
00436   getAdjFaceIntForce(force);
00437   // DEBUG
00438   if(print_switch==true)
00439   {
00440     double diff = sqrt(dot(force,force))-fmag;
00441     fmag = fmag+diff;
00442     cout << "Vertex::getTotalForceEnergy: face intersection force     ["
00443           << force[0]-f_old[0] << " "
00444           << force[1]-f_old[1] << " "
00445           << force[2]-f_old[2] << "]"
00446           << ", diff = " << diff << endl;
00447     cout << "Vertex::getTotalForceEnergy: resultant force ["
00448           << force[0] << " "
00449           << force[1] << " "
00450           << force[2] << "]\n\n";
00451   }
00452   // DEBUG
00453 }
00454 
00455 vec_d Vertex::getNewPos (double gain)
00456 {
00457   // get force
00458   double force[3]={0,0,0};
00459   getTotalForceEnergy(force);
00460   // compute new vertex coordinates
00461   vec_d pt(3,0);
00462   for (int k=0;k<3;k++) { pt[k] = gain*force[k]+p[k]; }
00463   // update container quantities
00464   Container::instance().addForce(sqrt(force[0]*force[0]+force[1]*force[1]+force[2]*force[2]));
00465   return pt;
00466 }
00467 
00468 void Vertex::getAdjVertices (vec_vp &a) const
00469 {
00470   a.clear();
00471   vec_ep e;
00472   getAdjacentEdges(e);
00473   // for each adjacent edge
00474   for (ep_it i=e.begin();i!=e.end();i++)
00475   {
00476     Vertex *v1=NULL,*v2=NULL,*o1=NULL,*o2=NULL;
00477     (*i)->getVertices(v1,v2,o1,o2);
00478     // find vertex different from self and add different vertex to vector
00479     if (v1!=this){a.push_back(v1);}
00480     else if (v2!=this) {a.push_back(v2);}
00481     else { printf("Error. both vertices of edge are equal to current vertex.\n"); exit(1); }
00482   }
00483 }
00484 
00485 //void Vertex:: getAdjacentFaces(hashset_f &fset) const
00486 //{
00487 //  vec_ep e;
00488 //  getAdjacentEdges(e);
00489 //  // for each adjacent edge
00490 //  for (ep_it i=e.begin();i!=e.end();i++)
00491 //  {
00492 //    // add edge faces to set
00493 //    fset.insert((*i)->f1);
00494 //    fset.insert((*i)->f2);
00495 //  }
00496 //}
00497 
00498 double Vertex::getSqSepDist (void) const
00499 {
00500   if (cl==NULL) return SEARCH_RADIUS_SQ;
00501   // get closest point
00502   //double pC[3];
00503   //vec_d pC = Container::instance().getNearPtOnFaceToVertex(cl,this);
00504   vec_fp fp(1,cl);
00505   vec_d pp;
00506   double sqd=0.0;
00507   Face *ncl = NULL;
00508   Container::instance().findClosestPtToVertexAmongFaces(this,fp,pp,sqd,ncl);
00509   //return (pC[0]-p[0])*(pC[0]-p[0])
00510   //      +(pC[1]-p[1])*(pC[1]-p[1])
00511   //      +(pC[2]-p[2])*(pC[2]-p[2]);
00512   return sqd;
00513 }
00514 
00515 void Vertex::setNewPos (vec_d const & pt)
00516 {
00517   double d=0;
00519   vec_vp v;
00520   getAdjVertices(v);
00521   // for each adjacent vertex
00522   for(vp_it i=v.begin();i!=v.end();i++)
00523   {
00524     // if adjacent vertex is indistinguishable from p, then displace p
00525     if ( !distinguishable((*i)->p[0],pt[0])&&
00526          !distinguishable((*i)->p[1],pt[1])&&
00527          !distinguishable((*i)->p[2],pt[2])){ d=DOUBLE_EPSILON*SCALE;}
00528   }
00529   p[0]=pt[0]+d;
00530   p[1]=pt[1]+d;
00531   p[2]=pt[2]+d;
00532 }
00533 
00534 void Vertex::getBoundingBox (double * const & bb) const
00535 {
00536   // assume bb=bounding box = [xmin xmax ymin ymax zmin zmax]
00537   bb[0]=bb[2]=bb[4]=1E30;
00538   bb[1]=bb[3]=bb[5]=-1E30;
00539   // for each adjacent face
00540   for(fp_cit i=f.begin();i!=f.end();i++)
00541   {
00542     // for each vertex of face
00543     for(int j=0;j<3;j++)
00544     {
00545       Vertex *vv=(*i)->getVertex(j);
00546       if      (vv->p[0]<bb[0]){bb[0]=vv->p[0];}
00547       else if (vv->p[0]>bb[1]){bb[1]=vv->p[0];}
00548       if      (vv->p[1]<bb[2]){bb[2]=vv->p[1];}
00549       else if (vv->p[1]>bb[3]){bb[3]=vv->p[1];}
00550       if      (vv->p[2]<bb[4]){bb[4]=vv->p[2];}
00551       else if (vv->p[2]>bb[5]){bb[5]=vv->p[2];}
00552     }
00553   }
00554 }
00555 
00556 void Vertex::getAdjacentEdges (vec_ep &vec) const
00557 {
00558   vec.clear();
00559   // for each adjacent face
00560   for(fp_cit i=f.begin();i!=f.end();i++)
00561   {
00562     // for each face edge
00563     for(int j=0;j<3;j++)
00564     {
00565       Vertex *v1=NULL,*v2=NULL,*o1=NULL,*o2=NULL;
00566       (*i)->getEdge(j)->getVertices(v1,v2,o1,o2);
00567       // if either edge vertex is current vertex
00568       if(v1==this || v2==this)
00569       {
00570         // add edge to vector
00571         vec.push_back((*i)->getEdge(j));
00572       }
00573     }
00574   }
00575 }
00576 
00577 vec_d const * Vertex::getPos (void) const
00578 {
00579   return &p;
00580 }
00581 
00582 void Vertex::setPos (double a, double b, double c)
00583 {
00584   p[0] = a;
00585   p[1] = b;
00586   p[2] = c;
00587 }
00588 
00589 double Vertex::getCoord (int i) const
00590 {
00591   return p[i];
00592 }
00593 
00594 Face const * Vertex::getFace (void) const
00595 {
00596   return cl;
00597 }
00598 
00599 void Vertex::setFace (Face * const ff)
00600 {
00601   cl=ff;
00602 }
00603 
00604 Object const * Vertex::getObject (void) const
00605 {
00606   return o;
00607 }
00608 
00609 void Vertex::setObject (Object * const oo)
00610 {
00611   o=oo;
00612 }
00613 
00614 fp_cit Vertex::begin (void) const
00615 {
00616   return f.begin();
00617 }
00618 
00619 fp_cit Vertex::end (void) const
00620 {
00621   return f.end();
00622 }
00623 
00624 void Vertex::addAdjacentFace (Face *ff)
00625 {
00626   f.push_back(ff);
00627 }
00628 
00629 Face* Vertex::getAdjacentFace (int i) const
00630 {
00631   return f[i];
00632 }
00633 
00634 int Vertex::getNumAdjFaces (void) const
00635 {
00636   return f.size();
00637 }
00638 
00639 void Vertex::sortAdjacentFaces (void)
00640 {
00641   sort(f.begin(),f.end());
00642 }
00643 
00644 double Vertex::getEnergy (void) const
00645 {
00646   return energy;
00647 }
00648 
00649 void Vertex::setEnergy (double ee)
00650 {
00651   energy = ee;
00652 }
00653 
00654 bool Vertex::faceIsAdjacent (Face const * const ff) const
00655 {
00656   return binary_search(f.begin(),f.end(),ff);
00657 }
00658 
00659 double Vertex::getSqVirtualDisp (double gain)
00660 {
00661   // compute new vertex coords
00662   // q = new holding position coordinates (x,y,z)
00663   vec_d q = getNewPos(gain);
00664   return (q[0]-getCoord(0))*(q[0]-getCoord(0))+
00665          (q[1]-getCoord(1))*(q[1]-getCoord(1))+
00666          (q[2]-getCoord(2))*(q[2]-getCoord(2));
00667 }
00668 
00669 Face* Vertex::getFaceNotAdjToVertex (Vertex const * const cv) const
00670 {
00671   // for each adjacent face of tv
00672   for(fp_cit i=begin();i!=end();i++)
00673   {
00674     // if face is not adjacent to vertex cv
00675     if(binary_search(cv->begin(),cv->end(),*i)==false)
00676     {
00677       return *i;
00678     }
00679   }
00680   cout << "Vertex::getFaceNotAdjToVertex: Error:"
00681         << "no acceptable face was found!\n";
00682   exit(0);
00683 }
00684 

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