cgma
FBClassify.cpp
Go to the documentation of this file.
00001 /*
00002  *
00003  *
00004  * Copyright (C) 2004 Sandia Corporation.  Under the terms of Contract DE-AC04-94AL85000
00005  * with Sandia Corporation, the U.S. Government retains certain rights in this software.
00006  *
00007  * This file is part of facetbool--contact via cubit@sandia.gov
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  * This library is distributed in the hope that it will be useful,
00015  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00016  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017  * Lesser General Public License for more details.
00018  *
00019  * You should have received a copy of the GNU Lesser General Public
00020  * License along with this library; if not, write to the Free Software
00021  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022  *
00023  *
00024  *
00025  */
00026 
00027 #include "CubitMessage.hpp"
00028 #include "FBClassify.hpp"
00029 #include "IntegerHash.hpp"
00030 #include "GfxDebug.hpp"
00031 #include <stack>
00032 //make this the same CUBIT_RESABS???
00033 //const double EPSILON_CLASSIFY = 1.e-12;
00034 const double EPSILON_CLASSIFY = 1.e-10;
00035 FBClassify::FBClassify()
00036 {
00037   number_of_groups = 0;
00038   polya = polyb = 0;
00039 }
00040 
00041 FBClassify::~FBClassify()
00042 {
00043 
00044 }
00045 
00046 void FBClassify::SetPoly(FBPolyhedron *poly1, FBPolyhedron *poly2)
00047 {
00048   polya = poly1;
00049   polyb = poly2;
00050 }
00051 
00052 CubitStatus FBClassify::Group(int which)
00053 {
00054 unsigned int i;
00055 int hashvalue, v0, v1, v2, tv0, tv1, tv2;
00056 //const int numhashbins = 101;
00057 int *hasharrayptr, hasharraysize, j, itri;
00058 FBPolyhedron *poly;
00059 const int primes[] = { 307, 601, 1009, 3001, 6007, 10007, 30011, 60013, 100003 };
00060 int numhashbins;
00061 
00062   if ( which == 1 ) poly = polya;
00063   else poly = polyb;
00064   
00065   if ( poly == 0 ) {
00066     PRINT_ERROR("ERROR:  Polyhedral object undefined in FBClassify::Group\n");
00067     return CUBIT_FAILURE;
00068   }
00069 
00070   i = poly->verts.size();
00071   if ( i < 3000 ) numhashbins = primes[0];
00072   else if ( i < 6000 ) numhashbins = primes[1];
00073   else if ( i < 10000 ) numhashbins = primes[2];
00074   else if ( i < 30000 ) numhashbins = primes[3];
00075   else if ( i < 60000 ) numhashbins = primes[4];
00076   else if ( i < 100000 ) numhashbins = primes[5];
00077   else if ( i < 300000 ) numhashbins = primes[6];
00078 //  else if ( i < 600000 ) numhashbins = primes[7];
00079   else numhashbins = primes[7]; 
00080 
00081   IntegerHash *hash = new IntegerHash(numhashbins,50);
00082   e0 = new int[poly->tris.size()];
00083   e1 = new int[poly->tris.size()];
00084   e2 = new int[poly->tris.size()];
00085 
00086   for ( i = 0; i < poly->tris.size(); i++ ) {
00087     e0[i] = e1[i] = e2[i]= NO_EDGE_NBR; // initialize 
00088     group.push_back(UNSET);
00089     v0 = poly->tris[i]->v0;
00090     v1 = poly->tris[i]->v1;
00091     v2 = poly->tris[i]->v2;
00092     //  Put the triangle sequence, i, in the hash list for each of the 3 edges.
00093     hash->addtoHashList((v0+v1)%numhashbins,i);
00094     hash->addtoHashList((v1+v2)%numhashbins,i);
00095     hash->addtoHashList((v2+v0)%numhashbins,i);
00096   }
00097 /*  int jmin, jmax, jave;
00098   jmin = 100000000;
00099   jmax = -100000000;
00100   jave = 0;
00101   for ( i = 0; i < poly->tris.size(); i++ ) {
00102     hasharrayptr = hash->getHashBin(i,&hasharraysize);
00103     if ( hasharraysize < jmin ) jmin = hasharraysize;
00104     if ( hasharraysize > jmax ) jmax = hasharraysize;
00105     jave += hasharraysize;
00106     
00107   } 
00108   jave /= poly->tris.size();
00109   printf("jmin jmax jave %d %d %d\n",jmin, jmax, jave);
00110   */
00111   
00112   for ( i = 0; i < poly->tris.size(); i++ ) {
00113     v0 = poly->tris[i]->v0;
00114     v1 = poly->tris[i]->v1;
00115     v2 = poly->tris[i]->v2;
00116     hashvalue = (v0+v1)%numhashbins;  // get the hash value for edge 0
00117     hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize);
00118     for ( j = 0; j < hasharraysize; j++ ) {  
00119       //  Go through the list and find the other triangle, itri, for each edge.
00120       //  Then assign its sequence to the edge array.
00121       itri = hasharrayptr[j];
00122       if ( (unsigned int)itri == i ) continue;
00123       tv0 = poly->tris[itri]->v0;
00124       tv1 = poly->tris[itri]->v1;
00125       tv2 = poly->tris[itri]->v2;
00126       if ( ((v1 == tv0) && (v0 == tv1)) || ((v0 == tv0) && (v1 == tv1)) ) {
00127         e0[i] = itri;
00128       } else if ( ((v1 == tv1) && (v0 == tv2)) || ((v0 == tv1) && (v1 == tv2)) ) {
00129         e0[i] = itri;
00130       } else if ( ((v1 == tv2) && (v0 == tv0)) || ((v0 == tv2) && (v1 == tv0)) ) {
00131         e0[i] = itri;
00132       } 
00133     } 
00134     hashvalue = (v1+v2)%numhashbins;
00135     hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize);
00136     for ( j = 0; j < hasharraysize; j++ ) {
00137       itri = hasharrayptr[j];
00138       if ( (unsigned int)itri == i ) continue;
00139       tv0 = poly->tris[itri]->v0;
00140       tv1 = poly->tris[itri]->v1;
00141       tv2 = poly->tris[itri]->v2;
00142       if ( ((v1 == tv0) && (v2 == tv1)) || ((v2 == tv0) && (v1 == tv1)) ) {
00143         e1[i] = itri;
00144       } else if ( ((v1 == tv1) && (v2 == tv2)) || ((v2 == tv1) && (v1 == tv2)) ) {
00145         e1[i] = itri;
00146       } else if ( ((v1 == tv2) && (v2 == tv0)) || ((v2 == tv2) && (v1 == tv0)) ) {
00147         e1[i] = itri;
00148       } 
00149     } 
00150     hashvalue = (v2+v0)%numhashbins;
00151     hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize);
00152     for ( j = 0; j < hasharraysize; j++ ) {
00153       itri = hasharrayptr[j];
00154       if ( (unsigned int)itri == i ) continue;
00155       tv0 = poly->tris[itri]->v0;
00156       tv1 = poly->tris[itri]->v1;
00157       tv2 = poly->tris[itri]->v2;
00158       if ( ((v0 == tv0) && (v2 == tv1)) || ((v2 == tv0) && (v0 == tv1)) ) {
00159         e2[i] = itri;
00160       } else if ( ((v0 == tv1) && (v2 == tv2)) || ((v2 == tv1) && (v0 == tv2)) ) {
00161         e2[i] = itri;
00162       } else if ( ((v0 == tv2) && (v2 == tv0)) || ((v2 == tv2) && (v0 == tv0)) ) {
00163         e2[i] = itri;
00164       }
00165     }       
00166   }
00167   //  Now we have to remove the other-side triangles where there was an intersection
00168   //  edge. 
00169   for ( i = 0; i < poly->intersection_edges.size(); i++ ) {
00170     v0 = poly->intersection_edges[i]->v0;
00171     v1 = poly->intersection_edges[i]->v1;
00172     hashvalue = (v0+v1)%numhashbins;
00173     hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize);
00174     for ( j = 0; j < hasharraysize; j++ ) {
00175       itri = hasharrayptr[j];
00176       tv0 = poly->tris[itri]->v0;
00177       tv1 = poly->tris[itri]->v1;
00178       tv2 = poly->tris[itri]->v2;
00179       if ( ((v0 == tv0) && (v1 == tv1)) || ((v0 == tv1) && (v1 == tv0)) ) {
00180         e0[itri] = NO_EDGE_NBR;
00181       };
00182       if ( ((v0 == tv0) && (v1 == tv2)) || ((v0 == tv2) && (v1 == tv0)) ) {
00183         e2[itri] = NO_EDGE_NBR;
00184       };
00185       if ( ((v0 == tv1) && (v1 == tv2)) || ((v0 == tv2) && (v1 == tv1)) ) {
00186         e1[itri] = NO_EDGE_NBR;
00187       }
00188     }    
00189   }
00190     
00191   //  Group the triangles that are neighbors.
00192   for ( i = 0; i < poly->tris.size(); i++ ) {
00193     if ( group[i] == UNSET ) {
00194       fill_group(i,number_of_groups++);
00195     }
00196   }
00197 
00198   delete[] e0; delete[] e1; delete[] e2; 
00199   delete hash;
00200     
00201   return CUBIT_SUCCESS;
00202 }
00203 
00204 void FBClassify::fill_group(int itri, int ngroup)
00205 {
00206 std::stack<int> vstack;
00207 int ktri;
00208 
00209   group[itri] = ngroup;
00210   
00211   if ( (e0[itri] != NO_EDGE_NBR) && (group[e0[itri]] == UNSET) ) 
00212     vstack.push(e0[itri]); 
00213   if ( (e1[itri] != NO_EDGE_NBR) && (group[e1[itri]] == UNSET) ) 
00214     vstack.push(e1[itri]); 
00215   if ( (e2[itri] != NO_EDGE_NBR) && (group[e2[itri]] == UNSET) ) 
00216     vstack.push(e2[itri]); 
00217   while (vstack.size() > 0 ) {
00218     ktri = vstack.top();
00219     vstack.pop();
00220     group[ktri] = ngroup;
00221     if ( (e0[ktri] != NO_EDGE_NBR) && (group[e0[ktri]] == UNSET) ) 
00222       vstack.push(e0[ktri]); 
00223     if ( (e1[ktri] != NO_EDGE_NBR) && (group[e1[ktri]] == UNSET) ) 
00224       vstack.push(e1[ktri]); 
00225     if ( (e2[ktri] != NO_EDGE_NBR) && (group[e2[ktri]] == UNSET) ) 
00226       vstack.push(e2[ktri]); 
00227   }
00228 }
00229 
00230 CubitStatus FBClassify::CharacterizeGroups(int which, bool other_is_planar)
00231 {
00232   int numberdone;
00233   unsigned int i;
00234   bool ifoundit;
00235   int itri = -1, type;
00236   FBPolyhedron *poly;
00237 
00238   if ( which == 1 ) poly = polya;
00239   else poly = polyb;
00240   
00241   if ( poly == 0 ) {
00242     PRINT_ERROR("ERROR:  Polyhedral object undefined in FBClassify::Group\n");
00243     return CUBIT_FAILURE;
00244   }
00245   
00246   numberdone = 0;
00247   i = 0; 
00248   while ( numberdone < number_of_groups ) {
00249     ifoundit = false;
00250     while ( i < poly->tris.size() ) {
00251       if ( group[i] == numberdone ) {
00252         ifoundit = true;
00253         itri = (int)i;
00254         break;      
00255       }
00256       i++;
00257     }
00258     if ( ifoundit == true ) {
00259       if ( other_is_planar == false ) 
00260         type = classify(itri, which);
00261       else 
00262         type = classify_against_plane(itri, which);      
00263       group_characterization.push_back(type);
00264       numberdone++;
00265     } else {
00266       PRINT_ERROR("Error in FBClassify::CharacterizeGroups\n");
00267       return CUBIT_FAILURE;
00268     }
00269   }
00270   return CUBIT_SUCCESS;
00271 }
00272 
00273 int FBClassify::classify_against_plane(int itri, int which)
00274 {
00275 FBPolyhedron *polyref, *polyobj; 
00276 int type, v0, v1, v2;
00277 double xbary, ybary, zbary, a, b, c;
00278 ;
00279 
00280   type = FB_ORIENTATION_UNDEFINED;
00281   if ( which == 1 ) {
00282     polyobj = polyb;
00283     polyref = polya;
00284   } else if ( which == 2 ) {
00285     polyobj = polya;
00286     polyref = polya;
00287   } else {
00288     PRINT_ERROR("ERROR in FBClassify::classify\n");
00289     return type;
00290   }
00291   v0 = polyref->tris[itri]->v0;
00292   v1 = polyref->tris[itri]->v1;
00293   v2 = polyref->tris[itri]->v2;
00294 
00295   xbary = ( polyref->verts[v0]->coord[0] + 
00296             polyref->verts[v1]->coord[0] + 
00297             polyref->verts[v2]->coord[0] )/3.;
00298   ybary = ( polyref->verts[v0]->coord[1] + 
00299             polyref->verts[v1]->coord[1] + 
00300             polyref->verts[v2]->coord[1] )/3.;
00301   zbary = ( polyref->verts[v0]->coord[2] + 
00302             polyref->verts[v1]->coord[2] + 
00303             polyref->verts[v2]->coord[2] )/3.;
00304   a = polyref->tris[itri]->a;
00305   b = polyref->tris[itri]->b;
00306   c = polyref->tris[itri]->c;
00307 
00308   //  Figure out which side of the plane we are on.  Since all
00309   //  of the plane's triangles have the same plane equation
00310   //  coefficients, might as well use the first one.
00311   
00312 double obj_tri_a, obj_tri_b, obj_tri_c, obj_tri_d, dotprod, disttoplane;
00313 
00314   obj_tri_a = polyobj->tris[0]->a;
00315   obj_tri_b = polyobj->tris[0]->b;
00316   obj_tri_c = polyobj->tris[0]->c;
00317   obj_tri_d = polyobj->tris[0]->d;
00318       
00319   disttoplane = obj_tri_a*xbary + obj_tri_b*ybary + obj_tri_c*zbary + obj_tri_d;    
00320   if ( disttoplane > EPSILON ) return FB_ORIENTATION_OUTSIDE;
00321   else if ( disttoplane < -EPSILON ) return FB_ORIENTATION_INSIDE; 
00322   
00323   dotprod = obj_tri_a*a + obj_tri_b*b + obj_tri_c*c;
00324   if ( dotprod > 0. ) return FB_ORIENTATION_SAME;
00325   else return FB_ORIENTATION_OPPOSITE;
00326   
00327   return type;
00328 }
00329 
00330 //returns an orientation for the triangle relative to the other body.
00331 //This triangle can be inside or outside the other body.  Or, it can
00332 //be opposite or same.  Opposite means this triangle is very close to
00333 // a triangle in the other body and it has an opposite pointing normal.
00334 // Same means it is very close to a trianglein the other body and it
00335 // their normals point in the same direction.
00336 int FBClassify::classify(int itri, int which)
00337 {
00338     //  "which" is the object, either 1 or 2, that the triangle itri
00339     // belongs to.
00340     //  Th e object to test for inside or outside is the other object.
00341   FBPolyhedron *polyref, *polyobj; 
00342     //  polyref = itri's polyhedron; polyobj = the other object
00343   int type, v0, v1, v2;
00344   double xbary, ybary, zbary, a, b, c;
00345 
00346   type = FB_ORIENTATION_UNDEFINED;
00347   if ( which == 1 ) {
00348     polyobj = polyb;
00349     polyref = polya;
00350   } else if ( which == 2 ) {
00351     polyobj = polya;
00352     polyref = polyb;
00353   } else {
00354     PRINT_ERROR("ERROR in FBClassify::classify\n");
00355     return type;
00356   }
00357   int mydebug = 0;
00358   if(mydebug)
00359     polyref->debug_draw_fb_triangle(polyref->tris[itri]);
00360 
00361   v0 = polyref->tris[itri]->v0;
00362   v1 = polyref->tris[itri]->v1;
00363   v2 = polyref->tris[itri]->v2;
00364 
00365   xbary = ( polyref->verts[v0]->coord[0] + 
00366             polyref->verts[v1]->coord[0] + 
00367             polyref->verts[v2]->coord[0] )/3.;
00368   ybary = ( polyref->verts[v0]->coord[1] + 
00369             polyref->verts[v1]->coord[1] + 
00370             polyref->verts[v2]->coord[1] )/3.;
00371   zbary = ( polyref->verts[v0]->coord[2] + 
00372             polyref->verts[v1]->coord[2] + 
00373             polyref->verts[v2]->coord[2] )/3.;
00374   a = polyref->tris[itri]->a;
00375   b = polyref->tris[itri]->b;
00376   c = polyref->tris[itri]->c;
00377 
00378   unsigned int i, num_perturb;
00379   double obj_tri_a, obj_tri_b, obj_tri_c, obj_tri_d, dotprod;
00380   double distance_to_other_sqr, closest_distance_to_plane, t;
00381   double closest_distance_to_other_sqr;
00382   double xint, yint, zint, distance_to_plane, closest_dotproduct;
00383   bool perturb, done, foundone;
00384   double other_xbar, other_ybar, other_zbar;
00385   int other_tri_0, other_tri_1, other_tri_2;
00386 
00387   perturb = false;
00388   num_perturb = 0;
00389   done = false;
00390   
00391   while ( (done == false) && (num_perturb < 20) ) {
00392     closest_dotproduct = -CUBIT_DBL_MAX + 1.;
00393     closest_distance_to_plane = CUBIT_DBL_MAX;
00394     closest_distance_to_other_sqr = CUBIT_DBL_MAX;
00395     foundone = false;
00396     for ( i = 0; i < polyobj->tris.size(); i++ ) {
00397       obj_tri_a = polyobj->tris[i]->a;
00398       obj_tri_b = polyobj->tris[i]->b;
00399       obj_tri_c = polyobj->tris[i]->c;
00400       obj_tri_d = polyobj->tris[i]->d;
00401   
00402       dotprod = obj_tri_a*a + obj_tri_b*b + obj_tri_c*c;
00403         //calculate the distance to the other triangles plane
00404       distance_to_plane = (obj_tri_a*xbary + obj_tri_b*ybary +
00405                            obj_tri_c*zbary + obj_tri_d);
00406        
00407       
00408       if ( fabs(dotprod) < EPSILON_CLASSIFY ) {
00409           //  Is the point in the plane?
00410         if ( fabs(distance_to_plane) < EPSILON_CLASSIFY ) {
00411             //  Perturb the ray and recast.
00412           perturb = true;
00413           num_perturb += 1;
00414           break;
00415         }
00416         continue;
00417       }
00418       
00419       t =-(distance_to_plane)/dotprod;
00420       if ( t < -EPSILON_CLASSIFY ) continue;
00421       xint = xbary + a*t;
00422       yint = ybary + b*t;
00423       zint = zbary + c*t;
00424 
00425         //  Check whether the intersection point lies in or on
00426         //  the object triangle's
00427         //  bounding box.
00428       if ( (polyobj->tris[i]->boundingbox.xmin - EPSILON > xint) || 
00429            (polyobj->tris[i]->boundingbox.xmax + EPSILON < xint) ||
00430            (polyobj->tris[i]->boundingbox.ymin - EPSILON > yint) || 
00431            (polyobj->tris[i]->boundingbox.ymax + EPSILON < yint) ||
00432            (polyobj->tris[i]->boundingbox.zmin - EPSILON > zint) || 
00433            (polyobj->tris[i]->boundingbox.zmax + EPSILON < zint) ) 
00434         continue;
00435       
00436         //  Is the point (xint, yint, zint) inside or on the triangle?
00437         //  Get a principal projection to make this a 2D problem.
00438       double xp1, yp1, xp2, yp2, xp3, yp3, ptx, pty;
00439       int retval;
00440 
00441       if ( (fabs(obj_tri_b) >= fabs(obj_tri_a)) && 
00442            (fabs(obj_tri_b) >= fabs(obj_tri_c)) ) {
00443         xp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[0];
00444         yp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[2];
00445         xp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[0];
00446         yp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[2];
00447         xp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[0];
00448         yp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[2];
00449         ptx = xint;
00450         pty = zint;        
00451       } else if ( fabs(obj_tri_a) >= fabs(obj_tri_c) ) {
00452         xp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[1];
00453         yp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[2];
00454         xp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[1];
00455         yp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[2];
00456         xp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[1];
00457         yp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[2];
00458         ptx = yint;
00459         pty = zint;   
00460       } else {
00461         xp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[0];
00462         yp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[1];
00463         xp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[0];
00464         yp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[1];
00465         xp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[0];
00466         yp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[1];
00467         ptx = xint;
00468         pty = yint;  
00469       }
00470       retval = pt_in_tri_2d(ptx,pty,xp1,yp1,xp2,yp2,xp3,yp3);
00471       if ( (retval == FB_ORIENTATION_INSIDE) ||
00472            (retval == FB_ORIENTATION_ON) ) {
00473           //calculate the distance to the other triangle's centroid
00474         other_tri_0 = polyobj->tris[i]->v0;
00475         other_tri_1 = polyobj->tris[i]->v1;
00476         other_tri_2 = polyobj->tris[i]->v2;
00477         other_xbar = ( polyobj->verts[other_tri_0]->coord[0] + 
00478                        polyobj->verts[other_tri_1]->coord[0] + 
00479                        polyobj->verts[other_tri_2]->coord[0] )/3.;
00480         other_ybar = ( polyobj->verts[other_tri_0]->coord[1] + 
00481                        polyobj->verts[other_tri_1]->coord[1] + 
00482                        polyobj->verts[other_tri_2]->coord[1] )/3.;
00483         other_zbar = ( polyobj->verts[other_tri_0]->coord[2] + 
00484                        polyobj->verts[other_tri_1]->coord[2] + 
00485                        polyobj->verts[other_tri_2]->coord[2] )/3.;
00486         
00487           //calculate the distance (squared) to the other triangle's centroid
00488         distance_to_other_sqr = ( (xbary-other_xbar)*(xbary-other_xbar) +
00489                                   (ybary-other_ybar)*(ybary-other_ybar) +
00490                                   (zbary-other_zbar)*(zbary-other_zbar) );
00491           //if this is the closest other triangle so far...
00492         if(closest_distance_to_other_sqr > distance_to_other_sqr){
00493             //then we found one, and update the closest distance, dot prod,
00494             // and distance to other plane.
00495           foundone = true;
00496           closest_distance_to_other_sqr = distance_to_other_sqr;
00497           closest_dotproduct = dotprod;
00498           if(mydebug){
00499             polyobj->debug_draw_fb_triangle(polyobj->tris[i]);
00500             GfxDebug::mouse_xforms();
00501           }
00502           closest_distance_to_plane = distance_to_plane;
00503             //  This is the closest triangle.
00504           if ( fabs(closest_distance_to_plane) < EPSILON_CLASSIFY ) 
00505             break;   
00506         }
00507       }       
00508     }
00509     if ( perturb == false ) done = true;
00510     else {
00511         //  perturb the ray and try again.
00512       perturb_the_ray(xbary, ybary, zbary);      
00513       perturb = false;
00514     }
00515   }
00516     //if we are very close to the plane of the closest other triangle
00517     // then we are either going to classify as opposite or same depending
00518     // on the relationship of the normals
00519   if ( (fabs(closest_distance_to_plane) < EPSILON_CLASSIFY ) ){
00520     if ( closest_dotproduct > 0.0 )
00521       type = FB_ORIENTATION_SAME;
00522     else
00523       type = FB_ORIENTATION_OPPOSITE;
00524   }
00525     //otherwise we are going to classify as inside or outside.  If we
00526     // didn't find any triangles that projected into our triangle,
00527     // we must be outside.  Otherwise we compare the normals to determine
00528     // whether we are inside or outside.
00529   else {
00530     if  ( foundone == false )
00531       type = FB_ORIENTATION_OUTSIDE;
00532     else if ( closest_dotproduct > 0.0 )
00533       type = FB_ORIENTATION_INSIDE;
00534     else
00535       type = FB_ORIENTATION_OUTSIDE;  
00536   }
00537   if(mydebug)
00538     GfxDebug::display_all();
00539   return type;
00540 }
00541 
00542 void FBClassify::perturb_the_ray(double &xbary, double &ybary, double &zbary)
00543 {
00544   xbary += 1.e-4*(double(rand())/(RAND_MAX+1.0)-0.5);
00545   ybary += 1.e-4*(double(rand())/(RAND_MAX+1.0)-0.5);
00546   zbary += 1.e-4*(double(rand())/(RAND_MAX+1.0)-0.5);
00547 }
00548 
00549 int FBClassify::pt_in_tri_2d(double xpt, double ypt,
00550                               double x0, double y0,
00551                       double x1, double y1,
00552                       double x2, double y2)
00553 {
00554 //  From Schneider & Eberly, "Geometric Tools for COmputer Graphics",
00555 //  Chap. 13.3.1.  If triangle is needle-thin, CUBIT_FAILURE might be
00556 //  returned, in wich case is_point_in is undefined.
00557 
00558 double c0, c1, c2;
00559 double e0x, e1x, e2x, e0y, e1y, e2y;
00560 double n0x, n1x, n2x, n0y, n1y, n2y;
00561 double denom0, denom1, denom2;
00562 int result;
00563 
00564   e0x = x1 - x0; e0y = y1 - y0;  
00565   e1x = x2 - x1; e1y = y2 - y1;  
00566   e2x = x0 - x2; e2y = y0 - y2;  
00567   n0x = e0y; n0y = -e0x;
00568   n1x = e1y; n1y = -e1x;
00569   n2x = e2y; n2y = -e2x;
00570   denom0 = n1x*e0x + n1y*e0y;
00571   if ( fabs(denom0) < EPSILON_CLASSIFY ) {
00572     PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n");
00573     return FB_ORIENTATION_UNDEFINED;
00574   }
00575   denom1 = n2x*e1x + n2y*e1y;
00576   if ( fabs(denom1) < EPSILON_CLASSIFY ) {
00577     PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n");
00578     return FB_ORIENTATION_UNDEFINED;
00579   }
00580   denom2 = n0x*e2x + n0y*e2y;
00581   if ( fabs(denom2) < EPSILON_CLASSIFY ) {
00582     PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n");
00583     return FB_ORIENTATION_UNDEFINED;
00584   }
00585   
00586   c0 = -( n1x*(xpt-x1) + n1y*(ypt-y1) )/denom0;
00587   c1 = -( n2x*(xpt-x2) + n2y*(ypt-y2) )/denom1;
00588   c2 = -( n0x*(xpt-x0) + n0y*(ypt-y0) )/denom2;
00589 
00590   if ( (c0 > 0.0) && (c1 > 0.0) && (c2 > 0.0) ) result = FB_ORIENTATION_INSIDE;
00591   else if ( (c0 < 0.0) || (c1 < 0.0) || (c2 < 0.0) ) result = FB_ORIENTATION_OUTSIDE;
00592   else result = FB_ORIENTATION_ON;
00593 
00594   return result;
00595 
00596 }
00597 
00598 void FBClassify::get_group(std::vector<int> **this_group,
00599                            std::vector<int> **this_group_characterization)
00600 {
00601   *this_group = &group;
00602   *this_group_characterization = &group_characterization;
00603 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines