nice.cc

Go to the documentation of this file.
00001 #include "nice.h"
00002 
00003 #include <cassert>
00004 #include <cmath>
00005 #include <iostream>
00006 
00007 #include "box.h"
00008 #include "container.h"
00009 #include "face.h"
00010 #include "misc.h"
00011 #include "object.h"
00012 #include "space.h"
00013 #include "vertex.h"
00014 
00015 using std::cout;
00016 using std::endl;
00017 
00018 Nice * Nice::only_one = NULL;
00019 
00020 Nice & Nice::instance()
00021 {
00022   // Not thread-safe.
00023   // -- lock mutex
00024   if (only_one == NULL)
00025     only_one = new Nice();
00026   // -- unlock mutex
00027   return *only_one;
00028 }
00029 
00030 Nice::Nice (const Nice & rhs)
00031   :nice()
00032 {
00033   cout << "Copy constructor prohibited on instances of Nice class.\n";
00034   cout << "Nice " << rhs.nice.size() << endl;
00035   exit(0);
00036 }
00037 
00038 Nice::Nice(void)
00039   :nice()
00040 {
00041 }
00042 
00043 Nice& Nice::operator= (const Nice& rhs)
00044 {
00045   cout << "Copy assignment operator prohibited on instances of Nice class.\n";
00046   cout << "Nice " << rhs.nice.size() << endl;
00047   exit(0);
00048 }
00049 
00050 bool Nice::vertexIsNice (Vertex const * const vv) const
00051 {
00052   return nice.find(const_cast<Vertex*>(vv))==nice.end();        
00053 }
00054 
00055 void Nice::setVertexNiceVal (int newval,Vertex * const vv)
00056 {
00057   // newval == 0, nice
00058   // newval == 1, nonnice in different object
00059   // newval == 2, nonnice in same object
00060   if (newval==0)
00061   {
00062     vhm_it i = nice.find(vv);
00063     if(i!=nice.end())
00064     {
00065       nice.erase(i);
00066     }
00067   }
00068   else
00069   {
00070       nice[vv]=newval;
00071   }
00072 }
00073 
00074 int Nice::getVertexNiceVal (Vertex const * const vv) const
00075 {
00076   if (vertexIsNice(const_cast<Vertex*>(vv)))
00077   {
00078     return 0;
00079   }
00080   else
00081   {
00082      vhm_cit i = nice.find(const_cast<Vertex*>(vv));    
00083      return (*i).second;
00084   }
00085 }
00086 
00087 int Nice::findExtraPoint (Vertex const * const v,vec_d & p,int const index) const
00088 {
00089   vec_d origin,end;
00090   getExtraRay(v,origin,end,index); // returns ray
00091   // look for intersected faces along ray, i.e. between face centroid and RAY ORIGIN
00092   vec_fp crossed_faces,unique_faces,edge_faces;
00093   collectNiceFaces(origin,end,unique_faces);
00094   findIntersectedFaces(origin,end,unique_faces,crossed_faces,edge_faces);
00095   sort(crossed_faces.begin(),crossed_faces.end());
00096   // look for adjacent face in crossed_faces
00097   std::pair<fp_it,fp_it> j;
00098   j=equal_range(crossed_faces.begin(),crossed_faces.end(),v->getAdjacentFace(index));
00099   // if found, then remove
00100   if(j.first!=j.second){crossed_faces.erase(j.first);}
00101   if( (edge_faces.size()%2)==1 )
00102   {
00103     cout << "\n\nError: how can # edge_faces be odd?\n"
00104           << "Nice::findExtraPoint: "
00105           << "current vertex:\n";
00106     v->print(std::cout);
00107     cout << "Nice::findExtraPoint: edge_faces:\n";
00108     for(fp_it i=edge_faces.begin();i!=edge_faces.end();i++)
00109     {
00110       (*i)->print(std::cout);
00111     }
00112     cout << endl;
00113     cout << "Nice::findExtraPoint: crossed_faces:\n";
00114     for(fp_it i=crossed_faces.begin();i!=crossed_faces.end();i++)
00115     {
00116       (*i)->print(std::cout);
00117     }
00118     cout << endl;
00119     exit(0);
00120   }
00121   vec_op tmp;
00122   findOddMeshes(crossed_faces,edge_faces,tmp);
00123   p[0]=end[0];
00124   p[1]=end[1];
00125   p[2]=end[2];
00126   if (!tmp.empty())
00127   {
00128     std::pair<op_it,op_it> i;
00129     i=equal_range(tmp.begin(),tmp.end(),v->getObject());
00130     if(i.first!=i.second)
00131     {
00132       return 0;
00133     }
00134   }
00135   return 1;
00136 }
00137 
00138 void Nice::findNonniceVertices (void)
00139 {
00140   cout << "Find nice vertices.............................";
00141   cout.flush();
00142   int k=1;
00143   double goal = 0.2;
00144   printf("0%%..");
00145   fflush(stdout);
00146   double a = 1.0/Container::instance().o.size();
00147   // for each object in container
00148   for (o_it i=Container::instance().o.begin();i!=Container::instance().o.end();i++)
00149   {
00150     // for each vertex in object
00151     for (v_it j=i->v.begin();j!=i->v.end();j++)
00152     {
00153       updateVertexNiceness(&(*j));
00154     }
00155     // track progress
00156     double progress = static_cast<double>(k++)*a;
00157     if(progress>goal)
00158     {
00159       printf("%d%%..",static_cast<int>(goal*100));
00160       fflush(stdout);
00161       goal+=0.2;
00162     }
00163   }
00164   printf("100%%..");
00165   fflush(stdout);
00166   cout << "complete.\n";
00167   cout.flush();
00168 }
00169 
00170 bool Nice::updateVertexNiceness (Vertex * const v)
00171 {
00172   vec_op crossed_objects;
00173   // collect objects inside which vertex lies
00174   collectCrossedObjects(v,crossed_objects);
00175   // update niceness of vertex based on crossed_objects
00176   return setVertexNiceness(v,crossed_objects);
00177 }
00178 
00179 
00180 bool Nice::getPointOutsideObject(Vertex const * const v,vec_d & p,vec_op & crossed_objects) const
00181 {
00182   // find point, p, outside of current object
00183   int index = 0;
00184   bool good = false;
00185   while(good==false)
00186   {
00187     if (index>(v->getNumAdjFaces()-1))
00188     {
00189       cout << "\n\nNice::collectCrossed: "
00190             << "Intersected edge along path from p to ECP on every adjacent face.\n"
00191             << "Warning! Setting undetermined vertex to 'nice'.\n";
00192       v->print(std::cout);
00193       crossed_objects.clear();
00194       return true;
00195     }
00196     while (!findExtraPoint(v,p,index))
00197     {
00198       index++;
00199       if (index>(v->getNumAdjFaces()-1))
00200       {
00201         cout << "\n\nNice::collectCrossed: "
00202               << "Failed to find valid extracellular point on any adjacent face.\n"
00203               << "Warning! Setting undetermined vertex to 'nice'.\n";
00204         v->print(std::cout);
00205         crossed_objects.clear();
00206         return true;
00207       }
00208     }
00209     // grab intersected objects between p and RAY ORIGIN and return as ca
00210     vec_d origin,end;
00211     origin.push_back(v->getCoord(0));
00212     origin.push_back(v->getCoord(1));
00213     origin.push_back(v->getCoord(2));
00214     end.push_back(p[0]);
00215     end.push_back(p[1]);
00216     end.push_back(p[2]);
00217     good=findCrossedObjects1(origin,end,crossed_objects,v);
00218     if(good==false){index++;}
00219   }
00220   return false;
00221 }
00222 
00223 void Nice::collectCrossedObjects (Vertex const * const v,vec_op & crossed_objects) const
00224 {
00225   vec_op ca;
00226   vec_d p(3,0);
00227   // find point, p, outside of current object
00228   // and grab intersected objects between p and RAY ORIGIN
00229   // and return as ca
00230   if(getPointOutsideObject(v,p,ca)==true) return;
00231   // grab intersected objects along RAY as cb
00232   vec_d origin(3,0),end(3,0);
00233   origin[0] = p[0];
00234   origin[1] = p[1];
00235   origin[2] = p[2];
00236   end[0]    = p[0];
00237   end[1]    = p[1];
00238   end[2]    = p[2];
00239   // get all rays sorted from shortest to longest
00240   double rays[6][3];
00241   findClosestAxis(v,origin,rays); // returns ray
00242   // initialize axis selection index
00243   // applies to list of distances to world boundary sorted shortest to longest
00244   // 0 == shortest direction, i.e. zeroth element
00245   // 1 == second shortest direction, i.e. first element
00246   // ...
00247   // 5 == longest direction, i.e. fifth element
00248   int select=0;
00249   bool keep = false;
00250   // while an axis needs to be found and select is less than max
00251   while(keep==false && select<6)
00252   {
00253     origin[0]=p[0];
00254     origin[1]=p[1];
00255     origin[2]=p[2];
00256     end[0]=rays[select][0];
00257     end[1]=rays[select][1];
00258     end[2]=rays[select][2];
00259     keep=findCrossedObjects2(origin,end,crossed_objects);
00260     if(keep==false){select++;}
00261   }
00262   if(select==6)
00263   {
00264     cout << "\n\nNice::collectCrossed: ERROR: All six principal coordinates"
00265           << " resulted in face edge collisions.\n";
00266     v->print(std::cout);
00267     cout << endl;
00268     exit(0);
00269   }
00270   // merge ca elements into crossed_objects
00271   sort(crossed_objects.begin(),crossed_objects.end());
00272   for (op_it i=ca.begin();i!=ca.end();i++)
00273   {
00274     std::pair<op_it,op_it> j;
00275     j=equal_range(crossed_objects.begin(),crossed_objects.end(),*i);
00276     // if object in ca is found in crossed_objects
00277     // then remove from crossed_objects
00278     if(j.first!=j.second){crossed_objects.erase(j.first);}
00279     // else add it
00280     else
00281     {
00282       crossed_objects.push_back(*i);
00283     }
00284   }
00285 }
00286 
00287 bool Nice::setVertexNiceness (Vertex * const v,vec_op &cb)
00288 {
00289   // new_nice==0 means nice
00290   // new_nice==1 means nonnice in different object
00291   // new_nice==2 means nonnice in same object
00292   int old_nice = getVertexNiceVal(v);
00293   int new_nice = 0;
00294   // if cb is not empty
00295   if (cb.empty()==false)
00296   {
00297     // vertex is not nice
00298     new_nice++;
00299     // if vertex is inside self object
00300     std::pair<op_it,op_it> i;
00301     i=equal_range(cb.begin(),cb.end(),v->getObject());
00302     if(i.first!=i.second){new_nice++;}
00303   }
00304   setVertexNiceVal(new_nice,v);
00305   // if vertex niceness changes then return true, else return false
00306   return old_nice!=new_nice;
00307 }
00308 
00309 void Nice::findClosestAxis (Vertex const * const v,vec_d const &origin,double rays[6][3]) const
00310 {
00311   // map: (distance to world bounds)double->(direction code)int
00312   map_di dist;
00313   dist.insert(std::make_pair(fabs(v->getCoord(0)-Space::instance().getWorld(0)),0));
00314   dist.insert(std::make_pair(fabs(v->getCoord(0)-Space::instance().getWorld(1)),1));
00315   dist.insert(std::make_pair(fabs(v->getCoord(1)-Space::instance().getWorld(2)),2));
00316   dist.insert(std::make_pair(fabs(v->getCoord(1)-Space::instance().getWorld(3)),3));
00317   dist.insert(std::make_pair(fabs(v->getCoord(2)-Space::instance().getWorld(4)),4));
00318   dist.insert(std::make_pair(fabs(v->getCoord(2)-Space::instance().getWorld(5)),5));
00319   int j=0;
00320   // for each map element
00321   for(di_it i=dist.begin();i!=dist.end();i++)
00322   {
00323     int index=(*i).second;
00324     if          (index==0)
00325     {
00326       rays[j][0] = origin[0]-2*(Space::instance().getWorld(1)-Space::instance().getWorld(0));   // end x
00327       rays[j][1] = origin[1];                                                   // end y
00328       rays[j][2] = origin[2];                                                   // end z
00329     }
00330     else if (index==1)
00331     {
00332       rays[j][0] = origin[0]+2*(Space::instance().getWorld(1)-Space::instance().getWorld(0));   // end x
00333       rays[j][1] = origin[1];                                                   // end y
00334       rays[j][2] = origin[2];                                                   // end z
00335     }
00336     else if (index==2)
00337     {
00338       rays[j][0] = origin[0];                                                   // end x
00339       rays[j][1] = origin[1]-2*(Space::instance().getWorld(3)-Space::instance().getWorld(2));   // end y
00340       rays[j][2] = origin[2];                                                   // end z
00341     }
00342     else if (index==3)
00343     {
00344       rays[j][0] = origin[0];                                                   // end x
00345       rays[j][1] = origin[1]+2*(Space::instance().getWorld(3)-Space::instance().getWorld(2));   // end y
00346       rays[j][2] = origin[2];                                                   // end z
00347     }
00348     else if (index==4)
00349     {
00350       rays[j][0] = origin[0];                                                   // end x
00351       rays[j][1] = origin[1];                                                   // end y
00352       rays[j][2] = origin[2]-2*(Space::instance().getWorld(5)-Space::instance().getWorld(4));   // end z
00353     }
00354     else if (index==5)
00355     {
00356       rays[j][0] = origin[0];                                                   // end x
00357       rays[j][1] = origin[1];                                                   // end y
00358       rays[j][2] = origin[2]+2*(Space::instance().getWorld(5)-Space::instance().getWorld(4));   // end z
00359     }
00360     j++;
00361   }
00362 }       
00363 
00364 bool Nice::findCrossedObjects1 (vec_d const & origin,vec_d const & end,vec_op &c,Vertex const * const v) const
00365 {
00366   // find and return crossed objects between origin point and end point
00367   vec_fp crossed_faces,unique_faces,edge_faces;
00368   collectNiceFaces(origin,end,unique_faces);
00369   findIntersectedFaces(origin,end,unique_faces,crossed_faces,edge_faces);
00370   // remove current vertex adjacent faces from crossed_faces
00371   if(crossed_faces.empty()==false)
00372   {
00373     sort(crossed_faces.begin(),crossed_faces.end());
00374     // for each adjacent face
00375     for(fp_cit i=v->begin();i!=v->end();i++)
00376     {
00377       std::pair<fp_it,fp_it> j; 
00378       j=equal_range(crossed_faces.begin(),crossed_faces.end(),*i);
00379       // if adjacent face is found in crossed_faces, then remove from crossed_faces
00380       if(j.first!=j.second){crossed_faces.erase(j.first);}
00381     }
00382   }
00383   // remove current vertex adjacent faces from edge_faces
00384   if(edge_faces.empty()==false)
00385   {
00386     sort(edge_faces.begin(),edge_faces.end());
00387     // for each adjacent face
00388     for(fp_cit i=v->begin();i!=v->end();i++)
00389     {
00390       std::pair<fp_it,fp_it> j; 
00391       j=equal_range(edge_faces.begin(),edge_faces.end(),*i);
00392       // if adjacent face is found in edge_faces, then remove from edge_faces
00393       if(j.first!=j.second){edge_faces.erase(j.first);}
00394     }
00395   }
00396 
00397   if(edge_faces.empty()==false){return false;}
00398 
00399   findOddMeshes(crossed_faces,edge_faces,c);
00400   return true;
00401 }
00402 
00403 bool Nice::findCrossedObjects2 (vec_d const & origin,vec_d const & end,vec_op &c) const
00404 {
00405   // find and return crossed objects between pN and extracellular point
00406   vec_fp crossed_faces,unique_faces,edge_faces;
00407   collectNiceFaces(origin,end,unique_faces);
00408   findIntersectedFaces(origin,end,unique_faces,crossed_faces,edge_faces);
00409   if(edge_faces.empty()==false){return false;}
00410   if( (edge_faces.size()%2)==1 )
00411   {     
00412     cout.precision(12);
00413     cout << "\n\nError: how can # edge_faces be odd?\n"
00414           << "Nice::findCrossedObjects2: "
00415           << "origin ["
00416           << origin[0] << " "
00417           << origin[1] << " "
00418           << origin[2] << "], end ["
00419           << end[0] << " "
00420           << end[1] << " "
00421           << end[2] << "]\n";
00422     cout << "Nice::findCrossed2: edge_faces:\n";
00423     for(fp_it i=edge_faces.begin();i!=edge_faces.end();i++)
00424     {
00425       (*i)->print(std::cout);
00426     }
00427     cout << endl;
00428     cout << "Nice::findCrossed2: crossed_faces:\n";
00429     for(fp_it i=crossed_faces.begin();i!=crossed_faces.end();i++)
00430     {
00431       (*i)->print(std::cout);
00432       cout << endl;
00433     }
00434     cout << endl;
00435     exit(0);
00436   }
00437   findOddMeshes(crossed_faces,edge_faces,c);
00438   return true;
00439 }
00440 
00441 void Nice::getExtraRay (Vertex const * const v,vec_d & origin,vec_d & end,const int index) const
00442 {
00443   // get normal info
00444   //double n[3];
00445   vec_d n = v->getAdjacentFace(index)->getNormal();
00446   // compute centroid of first adjacent face
00447   double a = 1.0/3.0;
00448   double cx = (v->getAdjacentFace(index)->getVertex(0)->getCoord(0)+
00449                v->getAdjacentFace(index)->getVertex(1)->getCoord(0)+
00450                v->getAdjacentFace(index)->getVertex(2)->getCoord(0))*a;
00451   double cy = (v->getAdjacentFace(index)->getVertex(0)->getCoord(1)+
00452                v->getAdjacentFace(index)->getVertex(1)->getCoord(1)+
00453                v->getAdjacentFace(index)->getVertex(2)->getCoord(1))*a;
00454   double cz = (v->getAdjacentFace(index)->getVertex(0)->getCoord(2)+
00455                v->getAdjacentFace(index)->getVertex(1)->getCoord(2)+
00456                v->getAdjacentFace(index)->getVertex(2)->getCoord(2))*a;
00457   double L = DOUBLE_EPSILON*SCALE/sqrt( dot(n,n) );
00458   origin.push_back(cx);
00459   origin.push_back(cy);
00460   origin.push_back(cz);
00461   end.push_back(origin[0]+n[0]*L);
00462   end.push_back(origin[1]+n[1]*L);
00463   end.push_back(origin[2]+n[2]*L);
00464 }
00465 
00466 void Nice::collectNiceFaces (vec_d const & origin,vec_d const & end,vec_fp & uf) const
00467 {
00468   // collect boxes to check
00469   vec_bp b;
00470   Space::instance().getBoxesFromPosition(origin,end,b);
00471   // gather faces in boxes
00472   // for each box
00473   for (bp_it i=b.begin();i!=b.end();i++)
00474   {
00475     // for each face in box
00476     for (fp_cit j=(*i)->begin();j!=(*i)->end();j++)
00477     {
00478       uf.push_back(*j);
00479     }
00480   }
00481   // keep unique faces
00482   sort(uf.begin(),uf.end());
00483   uf.assign(uf.begin(),unique(uf.begin(),uf.end()));
00484 }
00485 
00486 void Nice::findIntersectedFaces (vec_d const & origin,vec_d const & end,vec_fp &uf,vec_fp &cf,vec_fp &ef) const
00487 {
00488   // uf = unique_faces
00489   // cf = crossed_faces
00490   // ff = face_flags
00491   fp_it j;
00492   //bool line_flag, poly_flag, poly_edge_flag;
00493   // for each unique polygon
00494   for (j=uf.begin();j!=uf.end();j++)
00495   {
00496     //checkLineFaceInt(*j,origin,end,line_flag,poly_flag,poly_edge_flag,false);
00497     result r = checkLineFaceInt(*j,origin,end,false);
00498     // does point intersect polygon
00499     if (r.line_flag==true)
00500     {
00501       if (r.poly_edge_flag==true)
00502       {
00503         // add polygon_index to edge array
00504         ef.push_back(*j);
00505       }
00506       else if (r.poly_flag==true)
00507       {
00508         // add polygon_index to crossed array
00509         cf.push_back(*j);
00510       }
00511     }
00512   }
00513 }
00514 
00515 void Nice::findOddMeshes (vec_fp &cf,vec_fp &ef,vec_op &tmp) const
00516 {
00517   // find mesh objects crossed an odd number of times by ray
00518 
00519   vec_op ol; // object list
00520 
00521   // for each edge crossed face
00522   if(ef.empty()==false)
00523   {
00524     vec_op ee;
00525     for (fp_it i=ef.begin();i!=ef.end();i++)
00526     {
00527       ee.push_back(const_cast<Object*>((*i)->getVertex(0)->getObject()));
00528     }
00529     findEvenObjects(ee,ol);
00530   }
00531 
00532   // for each crossed face
00533   for (fp_it i=cf.begin();i!=cf.end();i++)
00534   {
00535     ol.push_back(const_cast<Object*>((*i)->getVertex(0)->getObject()));
00536   }
00537   findOddObjects(ol,tmp);
00538 }
00539 
00540 void Nice::findOddObjects (vec_op &ol,vec_op &odd) const
00541 {
00542   // sort object list
00543   sort(ol.begin(),ol.end());
00544 
00545   odd.clear();
00546   int count=ol.size();
00547   int k=0;
00548   // while more Object*s to process
00549   while (k<count)
00550   {
00551     // if not on last Object*
00552     if (k+1!=count)
00553     {
00554       // skip identical pairs of object indices
00555       if (ol[k]==ol[k+1]) {k++;k++;}
00556       // odd object
00557       else {
00558         odd.push_back( ol[k]);
00559         k++;
00560       }
00561     }
00562     else
00563     { // add remaining object to odd object list
00564       odd.push_back( ol[k]);
00565       k++;
00566     }
00567   }
00568 }
00569 
00570 void Nice::findEvenObjects (vec_op &ol,vec_op &even) const
00571 {
00572   // sort object list
00573   sort(ol.begin(),ol.end());
00574 
00575   even.clear();
00576   int count=ol.size();
00577   int k=0;
00578   // while more Object*s to process
00579   while (k<count)
00580   {
00581     // if not on last Object*
00582     if (k+1!=count)
00583     {
00584       // skip identical pairs of object indices
00585       if (ol[k]==ol[k+1])
00586       {
00587         even.push_back( ol[k]);
00588         k++;k++;
00589       }
00590       // odd object
00591       else
00592       {
00593         even.push_back( ol[k]);
00594         k++;
00595         cout << "\n\nNice::findEvenObjects: "
00596               << "Weird. Was expecting only even numbers of edge_faces.\n";
00597       }
00598     }
00599     else
00600     { // add remaining object to even object list
00601       even.push_back( ol[k]);
00602       cout << "\n\nNice::findEvenObjects: "
00603             << "Weird. Was expecting only even numbers of edge_faces.\n";
00604       k++;
00605     }
00606   }
00607 }
00608 
00609 int Nice::getNonniceCount (bool self_switch)
00610 {
00611   // if self_switch == true
00612   // then only count self intersections
00613   // else count all intersections
00614   int count = 0;
00615   // for each element of nice
00616   for (vhm_it j=nice.begin();j!=nice.end();j++)
00617   {
00618     assert( ((*j).second>0) && ((*j).second<3) ); // if zero then should have been removed
00619     if (self_switch==true)
00620     {
00621       if((*j).second==2)
00622       {
00623         count++;
00624       }
00625     }
00626     else
00627     {
00628       count++;
00629     }
00630   }
00631   return count;
00632 }
00633 
00634 vhm_cit Nice::beginNice (void) const
00635 {
00636   return nice.begin();
00637 }
00638 
00639 vhm_cit Nice::endNice (void) const
00640 {
00641   return nice.end();
00642 }
00643 

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