cgma
|
#include <FBClassify.hpp>
Public Member Functions | |
FBClassify () | |
~FBClassify () | |
void | SetPoly (FBPolyhedron *poly1, FBPolyhedron *poly2) |
CubitStatus | Group (int which) |
CubitStatus | CharacterizeGroups (int which, bool other_is_planar) |
void | get_group (std::vector< int > **this_group, std::vector< int > **this_group_characterization) |
Private Member Functions | |
void | fill_group (int itri, int ngroup) |
void | perturb_the_ray (double &xbary, double &ybary, double &zbary) |
int | pt_in_tri_2d (double xpt, double ypt, double x0, double y0, double x1, double y1, double x2, double y2) |
int | classify (int itri, int which) |
int | classify_against_plane (int itri, int which) |
Private Attributes | |
FBPolyhedron * | polya |
FBPolyhedron * | polyb |
std::vector< int > | group |
std::vector< int > | group_characterization |
int * | e0 |
int * | e1 |
int * | e2 |
int | number_of_groups |
Definition at line 38 of file FBClassify.hpp.
Definition at line 35 of file FBClassify.cpp.
{ number_of_groups = 0; polya = polyb = 0; }
Definition at line 41 of file FBClassify.cpp.
{ }
CubitStatus FBClassify::CharacterizeGroups | ( | int | which, |
bool | other_is_planar | ||
) |
Definition at line 230 of file FBClassify.cpp.
{ int numberdone; unsigned int i; bool ifoundit; int itri = -1, type; FBPolyhedron *poly; if ( which == 1 ) poly = polya; else poly = polyb; if ( poly == 0 ) { PRINT_ERROR("ERROR: Polyhedral object undefined in FBClassify::Group\n"); return CUBIT_FAILURE; } numberdone = 0; i = 0; while ( numberdone < number_of_groups ) { ifoundit = false; while ( i < poly->tris.size() ) { if ( group[i] == numberdone ) { ifoundit = true; itri = (int)i; break; } i++; } if ( ifoundit == true ) { if ( other_is_planar == false ) type = classify(itri, which); else type = classify_against_plane(itri, which); group_characterization.push_back(type); numberdone++; } else { PRINT_ERROR("Error in FBClassify::CharacterizeGroups\n"); return CUBIT_FAILURE; } } return CUBIT_SUCCESS; }
int FBClassify::classify | ( | int | itri, |
int | which | ||
) | [private] |
Definition at line 336 of file FBClassify.cpp.
{ // "which" is the object, either 1 or 2, that the triangle itri // belongs to. // Th e object to test for inside or outside is the other object. FBPolyhedron *polyref, *polyobj; // polyref = itri's polyhedron; polyobj = the other object int type, v0, v1, v2; double xbary, ybary, zbary, a, b, c; type = FB_ORIENTATION_UNDEFINED; if ( which == 1 ) { polyobj = polyb; polyref = polya; } else if ( which == 2 ) { polyobj = polya; polyref = polyb; } else { PRINT_ERROR("ERROR in FBClassify::classify\n"); return type; } int mydebug = 0; if(mydebug) polyref->debug_draw_fb_triangle(polyref->tris[itri]); v0 = polyref->tris[itri]->v0; v1 = polyref->tris[itri]->v1; v2 = polyref->tris[itri]->v2; xbary = ( polyref->verts[v0]->coord[0] + polyref->verts[v1]->coord[0] + polyref->verts[v2]->coord[0] )/3.; ybary = ( polyref->verts[v0]->coord[1] + polyref->verts[v1]->coord[1] + polyref->verts[v2]->coord[1] )/3.; zbary = ( polyref->verts[v0]->coord[2] + polyref->verts[v1]->coord[2] + polyref->verts[v2]->coord[2] )/3.; a = polyref->tris[itri]->a; b = polyref->tris[itri]->b; c = polyref->tris[itri]->c; unsigned int i, num_perturb; double obj_tri_a, obj_tri_b, obj_tri_c, obj_tri_d, dotprod; double distance_to_other_sqr, closest_distance_to_plane, t; double closest_distance_to_other_sqr; double xint, yint, zint, distance_to_plane, closest_dotproduct; bool perturb, done, foundone; double other_xbar, other_ybar, other_zbar; int other_tri_0, other_tri_1, other_tri_2; perturb = false; num_perturb = 0; done = false; while ( (done == false) && (num_perturb < 20) ) { closest_dotproduct = -CUBIT_DBL_MAX + 1.; closest_distance_to_plane = CUBIT_DBL_MAX; closest_distance_to_other_sqr = CUBIT_DBL_MAX; foundone = false; for ( i = 0; i < polyobj->tris.size(); i++ ) { obj_tri_a = polyobj->tris[i]->a; obj_tri_b = polyobj->tris[i]->b; obj_tri_c = polyobj->tris[i]->c; obj_tri_d = polyobj->tris[i]->d; dotprod = obj_tri_a*a + obj_tri_b*b + obj_tri_c*c; //calculate the distance to the other triangles plane distance_to_plane = (obj_tri_a*xbary + obj_tri_b*ybary + obj_tri_c*zbary + obj_tri_d); if ( fabs(dotprod) < EPSILON_CLASSIFY ) { // Is the point in the plane? if ( fabs(distance_to_plane) < EPSILON_CLASSIFY ) { // Perturb the ray and recast. perturb = true; num_perturb += 1; break; } continue; } t =-(distance_to_plane)/dotprod; if ( t < -EPSILON_CLASSIFY ) continue; xint = xbary + a*t; yint = ybary + b*t; zint = zbary + c*t; // Check whether the intersection point lies in or on // the object triangle's // bounding box. if ( (polyobj->tris[i]->boundingbox.xmin - EPSILON > xint) || (polyobj->tris[i]->boundingbox.xmax + EPSILON < xint) || (polyobj->tris[i]->boundingbox.ymin - EPSILON > yint) || (polyobj->tris[i]->boundingbox.ymax + EPSILON < yint) || (polyobj->tris[i]->boundingbox.zmin - EPSILON > zint) || (polyobj->tris[i]->boundingbox.zmax + EPSILON < zint) ) continue; // Is the point (xint, yint, zint) inside or on the triangle? // Get a principal projection to make this a 2D problem. double xp1, yp1, xp2, yp2, xp3, yp3, ptx, pty; int retval; if ( (fabs(obj_tri_b) >= fabs(obj_tri_a)) && (fabs(obj_tri_b) >= fabs(obj_tri_c)) ) { xp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[0]; yp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[2]; xp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[0]; yp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[2]; xp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[0]; yp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[2]; ptx = xint; pty = zint; } else if ( fabs(obj_tri_a) >= fabs(obj_tri_c) ) { xp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[1]; yp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[2]; xp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[1]; yp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[2]; xp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[1]; yp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[2]; ptx = yint; pty = zint; } else { xp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[0]; yp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[1]; xp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[0]; yp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[1]; xp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[0]; yp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[1]; ptx = xint; pty = yint; } retval = pt_in_tri_2d(ptx,pty,xp1,yp1,xp2,yp2,xp3,yp3); if ( (retval == FB_ORIENTATION_INSIDE) || (retval == FB_ORIENTATION_ON) ) { //calculate the distance to the other triangle's centroid other_tri_0 = polyobj->tris[i]->v0; other_tri_1 = polyobj->tris[i]->v1; other_tri_2 = polyobj->tris[i]->v2; other_xbar = ( polyobj->verts[other_tri_0]->coord[0] + polyobj->verts[other_tri_1]->coord[0] + polyobj->verts[other_tri_2]->coord[0] )/3.; other_ybar = ( polyobj->verts[other_tri_0]->coord[1] + polyobj->verts[other_tri_1]->coord[1] + polyobj->verts[other_tri_2]->coord[1] )/3.; other_zbar = ( polyobj->verts[other_tri_0]->coord[2] + polyobj->verts[other_tri_1]->coord[2] + polyobj->verts[other_tri_2]->coord[2] )/3.; //calculate the distance (squared) to the other triangle's centroid distance_to_other_sqr = ( (xbary-other_xbar)*(xbary-other_xbar) + (ybary-other_ybar)*(ybary-other_ybar) + (zbary-other_zbar)*(zbary-other_zbar) ); //if this is the closest other triangle so far... if(closest_distance_to_other_sqr > distance_to_other_sqr){ //then we found one, and update the closest distance, dot prod, // and distance to other plane. foundone = true; closest_distance_to_other_sqr = distance_to_other_sqr; closest_dotproduct = dotprod; if(mydebug){ polyobj->debug_draw_fb_triangle(polyobj->tris[i]); GfxDebug::mouse_xforms(); } closest_distance_to_plane = distance_to_plane; // This is the closest triangle. if ( fabs(closest_distance_to_plane) < EPSILON_CLASSIFY ) break; } } } if ( perturb == false ) done = true; else { // perturb the ray and try again. perturb_the_ray(xbary, ybary, zbary); perturb = false; } } //if we are very close to the plane of the closest other triangle // then we are either going to classify as opposite or same depending // on the relationship of the normals if ( (fabs(closest_distance_to_plane) < EPSILON_CLASSIFY ) ){ if ( closest_dotproduct > 0.0 ) type = FB_ORIENTATION_SAME; else type = FB_ORIENTATION_OPPOSITE; } //otherwise we are going to classify as inside or outside. If we // didn't find any triangles that projected into our triangle, // we must be outside. Otherwise we compare the normals to determine // whether we are inside or outside. else { if ( foundone == false ) type = FB_ORIENTATION_OUTSIDE; else if ( closest_dotproduct > 0.0 ) type = FB_ORIENTATION_INSIDE; else type = FB_ORIENTATION_OUTSIDE; } if(mydebug) GfxDebug::display_all(); return type; }
int FBClassify::classify_against_plane | ( | int | itri, |
int | which | ||
) | [private] |
Definition at line 273 of file FBClassify.cpp.
{ FBPolyhedron *polyref, *polyobj; int type, v0, v1, v2; double xbary, ybary, zbary, a, b, c; ; type = FB_ORIENTATION_UNDEFINED; if ( which == 1 ) { polyobj = polyb; polyref = polya; } else if ( which == 2 ) { polyobj = polya; polyref = polya; } else { PRINT_ERROR("ERROR in FBClassify::classify\n"); return type; } v0 = polyref->tris[itri]->v0; v1 = polyref->tris[itri]->v1; v2 = polyref->tris[itri]->v2; xbary = ( polyref->verts[v0]->coord[0] + polyref->verts[v1]->coord[0] + polyref->verts[v2]->coord[0] )/3.; ybary = ( polyref->verts[v0]->coord[1] + polyref->verts[v1]->coord[1] + polyref->verts[v2]->coord[1] )/3.; zbary = ( polyref->verts[v0]->coord[2] + polyref->verts[v1]->coord[2] + polyref->verts[v2]->coord[2] )/3.; a = polyref->tris[itri]->a; b = polyref->tris[itri]->b; c = polyref->tris[itri]->c; // Figure out which side of the plane we are on. Since all // of the plane's triangles have the same plane equation // coefficients, might as well use the first one. double obj_tri_a, obj_tri_b, obj_tri_c, obj_tri_d, dotprod, disttoplane; obj_tri_a = polyobj->tris[0]->a; obj_tri_b = polyobj->tris[0]->b; obj_tri_c = polyobj->tris[0]->c; obj_tri_d = polyobj->tris[0]->d; disttoplane = obj_tri_a*xbary + obj_tri_b*ybary + obj_tri_c*zbary + obj_tri_d; if ( disttoplane > EPSILON ) return FB_ORIENTATION_OUTSIDE; else if ( disttoplane < -EPSILON ) return FB_ORIENTATION_INSIDE; dotprod = obj_tri_a*a + obj_tri_b*b + obj_tri_c*c; if ( dotprod > 0. ) return FB_ORIENTATION_SAME; else return FB_ORIENTATION_OPPOSITE; return type; }
void FBClassify::fill_group | ( | int | itri, |
int | ngroup | ||
) | [private] |
Definition at line 204 of file FBClassify.cpp.
{ std::stack<int> vstack; int ktri; group[itri] = ngroup; if ( (e0[itri] != NO_EDGE_NBR) && (group[e0[itri]] == UNSET) ) vstack.push(e0[itri]); if ( (e1[itri] != NO_EDGE_NBR) && (group[e1[itri]] == UNSET) ) vstack.push(e1[itri]); if ( (e2[itri] != NO_EDGE_NBR) && (group[e2[itri]] == UNSET) ) vstack.push(e2[itri]); while (vstack.size() > 0 ) { ktri = vstack.top(); vstack.pop(); group[ktri] = ngroup; if ( (e0[ktri] != NO_EDGE_NBR) && (group[e0[ktri]] == UNSET) ) vstack.push(e0[ktri]); if ( (e1[ktri] != NO_EDGE_NBR) && (group[e1[ktri]] == UNSET) ) vstack.push(e1[ktri]); if ( (e2[ktri] != NO_EDGE_NBR) && (group[e2[ktri]] == UNSET) ) vstack.push(e2[ktri]); } }
void FBClassify::get_group | ( | std::vector< int > ** | this_group, |
std::vector< int > ** | this_group_characterization | ||
) |
Definition at line 598 of file FBClassify.cpp.
{ *this_group = &group; *this_group_characterization = &group_characterization; }
CubitStatus FBClassify::Group | ( | int | which | ) |
Definition at line 52 of file FBClassify.cpp.
{ unsigned int i; int hashvalue, v0, v1, v2, tv0, tv1, tv2; //const int numhashbins = 101; int *hasharrayptr, hasharraysize, j, itri; FBPolyhedron *poly; const int primes[] = { 307, 601, 1009, 3001, 6007, 10007, 30011, 60013, 100003 }; int numhashbins; if ( which == 1 ) poly = polya; else poly = polyb; if ( poly == 0 ) { PRINT_ERROR("ERROR: Polyhedral object undefined in FBClassify::Group\n"); return CUBIT_FAILURE; } i = poly->verts.size(); if ( i < 3000 ) numhashbins = primes[0]; else if ( i < 6000 ) numhashbins = primes[1]; else if ( i < 10000 ) numhashbins = primes[2]; else if ( i < 30000 ) numhashbins = primes[3]; else if ( i < 60000 ) numhashbins = primes[4]; else if ( i < 100000 ) numhashbins = primes[5]; else if ( i < 300000 ) numhashbins = primes[6]; // else if ( i < 600000 ) numhashbins = primes[7]; else numhashbins = primes[7]; IntegerHash *hash = new IntegerHash(numhashbins,50); e0 = new int[poly->tris.size()]; e1 = new int[poly->tris.size()]; e2 = new int[poly->tris.size()]; for ( i = 0; i < poly->tris.size(); i++ ) { e0[i] = e1[i] = e2[i]= NO_EDGE_NBR; // initialize group.push_back(UNSET); v0 = poly->tris[i]->v0; v1 = poly->tris[i]->v1; v2 = poly->tris[i]->v2; // Put the triangle sequence, i, in the hash list for each of the 3 edges. hash->addtoHashList((v0+v1)%numhashbins,i); hash->addtoHashList((v1+v2)%numhashbins,i); hash->addtoHashList((v2+v0)%numhashbins,i); } /* int jmin, jmax, jave; jmin = 100000000; jmax = -100000000; jave = 0; for ( i = 0; i < poly->tris.size(); i++ ) { hasharrayptr = hash->getHashBin(i,&hasharraysize); if ( hasharraysize < jmin ) jmin = hasharraysize; if ( hasharraysize > jmax ) jmax = hasharraysize; jave += hasharraysize; } jave /= poly->tris.size(); printf("jmin jmax jave %d %d %d\n",jmin, jmax, jave); */ for ( i = 0; i < poly->tris.size(); i++ ) { v0 = poly->tris[i]->v0; v1 = poly->tris[i]->v1; v2 = poly->tris[i]->v2; hashvalue = (v0+v1)%numhashbins; // get the hash value for edge 0 hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize); for ( j = 0; j < hasharraysize; j++ ) { // Go through the list and find the other triangle, itri, for each edge. // Then assign its sequence to the edge array. itri = hasharrayptr[j]; if ( (unsigned int)itri == i ) continue; tv0 = poly->tris[itri]->v0; tv1 = poly->tris[itri]->v1; tv2 = poly->tris[itri]->v2; if ( ((v1 == tv0) && (v0 == tv1)) || ((v0 == tv0) && (v1 == tv1)) ) { e0[i] = itri; } else if ( ((v1 == tv1) && (v0 == tv2)) || ((v0 == tv1) && (v1 == tv2)) ) { e0[i] = itri; } else if ( ((v1 == tv2) && (v0 == tv0)) || ((v0 == tv2) && (v1 == tv0)) ) { e0[i] = itri; } } hashvalue = (v1+v2)%numhashbins; hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize); for ( j = 0; j < hasharraysize; j++ ) { itri = hasharrayptr[j]; if ( (unsigned int)itri == i ) continue; tv0 = poly->tris[itri]->v0; tv1 = poly->tris[itri]->v1; tv2 = poly->tris[itri]->v2; if ( ((v1 == tv0) && (v2 == tv1)) || ((v2 == tv0) && (v1 == tv1)) ) { e1[i] = itri; } else if ( ((v1 == tv1) && (v2 == tv2)) || ((v2 == tv1) && (v1 == tv2)) ) { e1[i] = itri; } else if ( ((v1 == tv2) && (v2 == tv0)) || ((v2 == tv2) && (v1 == tv0)) ) { e1[i] = itri; } } hashvalue = (v2+v0)%numhashbins; hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize); for ( j = 0; j < hasharraysize; j++ ) { itri = hasharrayptr[j]; if ( (unsigned int)itri == i ) continue; tv0 = poly->tris[itri]->v0; tv1 = poly->tris[itri]->v1; tv2 = poly->tris[itri]->v2; if ( ((v0 == tv0) && (v2 == tv1)) || ((v2 == tv0) && (v0 == tv1)) ) { e2[i] = itri; } else if ( ((v0 == tv1) && (v2 == tv2)) || ((v2 == tv1) && (v0 == tv2)) ) { e2[i] = itri; } else if ( ((v0 == tv2) && (v2 == tv0)) || ((v2 == tv2) && (v0 == tv0)) ) { e2[i] = itri; } } } // Now we have to remove the other-side triangles where there was an intersection // edge. for ( i = 0; i < poly->intersection_edges.size(); i++ ) { v0 = poly->intersection_edges[i]->v0; v1 = poly->intersection_edges[i]->v1; hashvalue = (v0+v1)%numhashbins; hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize); for ( j = 0; j < hasharraysize; j++ ) { itri = hasharrayptr[j]; tv0 = poly->tris[itri]->v0; tv1 = poly->tris[itri]->v1; tv2 = poly->tris[itri]->v2; if ( ((v0 == tv0) && (v1 == tv1)) || ((v0 == tv1) && (v1 == tv0)) ) { e0[itri] = NO_EDGE_NBR; }; if ( ((v0 == tv0) && (v1 == tv2)) || ((v0 == tv2) && (v1 == tv0)) ) { e2[itri] = NO_EDGE_NBR; }; if ( ((v0 == tv1) && (v1 == tv2)) || ((v0 == tv2) && (v1 == tv1)) ) { e1[itri] = NO_EDGE_NBR; } } } // Group the triangles that are neighbors. for ( i = 0; i < poly->tris.size(); i++ ) { if ( group[i] == UNSET ) { fill_group(i,number_of_groups++); } } delete[] e0; delete[] e1; delete[] e2; delete hash; return CUBIT_SUCCESS; }
void FBClassify::perturb_the_ray | ( | double & | xbary, |
double & | ybary, | ||
double & | zbary | ||
) | [private] |
Definition at line 542 of file FBClassify.cpp.
{ xbary += 1.e-4*(double(rand())/(RAND_MAX+1.0)-0.5); ybary += 1.e-4*(double(rand())/(RAND_MAX+1.0)-0.5); zbary += 1.e-4*(double(rand())/(RAND_MAX+1.0)-0.5); }
int FBClassify::pt_in_tri_2d | ( | double | xpt, |
double | ypt, | ||
double | x0, | ||
double | y0, | ||
double | x1, | ||
double | y1, | ||
double | x2, | ||
double | y2 | ||
) | [private] |
Definition at line 549 of file FBClassify.cpp.
{ // From Schneider & Eberly, "Geometric Tools for COmputer Graphics", // Chap. 13.3.1. If triangle is needle-thin, CUBIT_FAILURE might be // returned, in wich case is_point_in is undefined. double c0, c1, c2; double e0x, e1x, e2x, e0y, e1y, e2y; double n0x, n1x, n2x, n0y, n1y, n2y; double denom0, denom1, denom2; int result; e0x = x1 - x0; e0y = y1 - y0; e1x = x2 - x1; e1y = y2 - y1; e2x = x0 - x2; e2y = y0 - y2; n0x = e0y; n0y = -e0x; n1x = e1y; n1y = -e1x; n2x = e2y; n2y = -e2x; denom0 = n1x*e0x + n1y*e0y; if ( fabs(denom0) < EPSILON_CLASSIFY ) { PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n"); return FB_ORIENTATION_UNDEFINED; } denom1 = n2x*e1x + n2y*e1y; if ( fabs(denom1) < EPSILON_CLASSIFY ) { PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n"); return FB_ORIENTATION_UNDEFINED; } denom2 = n0x*e2x + n0y*e2y; if ( fabs(denom2) < EPSILON_CLASSIFY ) { PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n"); return FB_ORIENTATION_UNDEFINED; } c0 = -( n1x*(xpt-x1) + n1y*(ypt-y1) )/denom0; c1 = -( n2x*(xpt-x2) + n2y*(ypt-y2) )/denom1; c2 = -( n0x*(xpt-x0) + n0y*(ypt-y0) )/denom2; if ( (c0 > 0.0) && (c1 > 0.0) && (c2 > 0.0) ) result = FB_ORIENTATION_INSIDE; else if ( (c0 < 0.0) || (c1 < 0.0) || (c2 < 0.0) ) result = FB_ORIENTATION_OUTSIDE; else result = FB_ORIENTATION_ON; return result; }
void FBClassify::SetPoly | ( | FBPolyhedron * | poly1, |
FBPolyhedron * | poly2 | ||
) |
Definition at line 46 of file FBClassify.cpp.
int* FBClassify::e0 [private] |
Definition at line 59 of file FBClassify.hpp.
int * FBClassify::e1 [private] |
Definition at line 59 of file FBClassify.hpp.
int * FBClassify::e2 [private] |
Definition at line 59 of file FBClassify.hpp.
std::vector<int> FBClassify::group [private] |
Definition at line 52 of file FBClassify.hpp.
std::vector<int> FBClassify::group_characterization [private] |
Definition at line 52 of file FBClassify.hpp.
int FBClassify::number_of_groups [private] |
Definition at line 60 of file FBClassify.hpp.
FBPolyhedron* FBClassify::polya [private] |
Definition at line 51 of file FBClassify.hpp.
FBPolyhedron * FBClassify::polyb [private] |
Definition at line 51 of file FBClassify.hpp.