Branch data Line data Source code
1 : : /*
2 : : *
3 : : *
4 : : * Copyright (C) 2004 Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000
5 : : * with Sandia Corporation, the U.S. Government retains certain rights in this software.
6 : : *
7 : : * This file is part of facetbool--contact via [email protected]
8 : : *
9 : : * This library is free software; you can redistribute it and/or
10 : : * modify it under the terms of the GNU Lesser General Public
11 : : * License as published by the Free Software Foundation; either
12 : : * version 2.1 of the License, or (at your option) any later version.
13 : : *
14 : : * This library is distributed in the hope that it will be useful,
15 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : * Lesser General Public License for more details.
18 : : *
19 : : * You should have received a copy of the GNU Lesser General Public
20 : : * License along with this library; if not, write to the Free Software
21 : : * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : : *
23 : : *
24 : : *
25 : : */
26 : :
27 : : #include "CubitMessage.hpp"
28 : : #include "FBClassify.hpp"
29 : : #include "IntegerHash.hpp"
30 : : #include "GfxDebug.hpp"
31 : : #include <stack>
32 : : //make this the same CUBIT_RESABS???
33 : : //const double EPSILON_CLASSIFY = 1.e-12;
34 : : const double EPSILON_CLASSIFY = 1.e-10;
35 [ # # ]: 0 : FBClassify::FBClassify()
36 : : {
37 : 0 : number_of_groups = 0;
38 : 0 : polya = polyb = 0;
39 : 0 : }
40 : :
41 [ # # ]: 0 : FBClassify::~FBClassify()
42 : : {
43 : :
44 : 0 : }
45 : :
46 : 0 : void FBClassify::SetPoly(FBPolyhedron *poly1, FBPolyhedron *poly2)
47 : : {
48 : 0 : polya = poly1;
49 : 0 : polyb = poly2;
50 : 0 : }
51 : :
52 : 0 : CubitStatus FBClassify::Group(int which)
53 : : {
54 : : unsigned int i;
55 : : int hashvalue, v0, v1, v2, tv0, tv1, tv2;
56 : : //const int numhashbins = 101;
57 : : int *hasharrayptr, hasharraysize, j, itri;
58 : : FBPolyhedron *poly;
59 : : const int primes[] = { 307, 601, 1009, 3001, 6007, 10007, 30011, 60013, 100003 };
60 : : int numhashbins;
61 : :
62 [ # # ]: 0 : if ( which == 1 ) poly = polya;
63 : 0 : else poly = polyb;
64 : :
65 [ # # ]: 0 : if ( poly == 0 ) {
66 [ # # ][ # # ]: 0 : PRINT_ERROR("ERROR: Polyhedral object undefined in FBClassify::Group\n");
[ # # ][ # # ]
67 : 0 : return CUBIT_FAILURE;
68 : : }
69 : :
70 [ # # ]: 0 : i = poly->verts.size();
71 [ # # ]: 0 : if ( i < 3000 ) numhashbins = primes[0];
72 [ # # ]: 0 : else if ( i < 6000 ) numhashbins = primes[1];
73 [ # # ]: 0 : else if ( i < 10000 ) numhashbins = primes[2];
74 [ # # ]: 0 : else if ( i < 30000 ) numhashbins = primes[3];
75 [ # # ]: 0 : else if ( i < 60000 ) numhashbins = primes[4];
76 [ # # ]: 0 : else if ( i < 100000 ) numhashbins = primes[5];
77 [ # # ]: 0 : else if ( i < 300000 ) numhashbins = primes[6];
78 : : // else if ( i < 600000 ) numhashbins = primes[7];
79 : 0 : else numhashbins = primes[7];
80 : :
81 [ # # ][ # # ]: 0 : IntegerHash *hash = new IntegerHash(numhashbins,50);
82 [ # # ][ # # ]: 0 : e0 = new int[poly->tris.size()];
[ # # ]
83 [ # # ][ # # ]: 0 : e1 = new int[poly->tris.size()];
[ # # ]
84 [ # # ][ # # ]: 0 : e2 = new int[poly->tris.size()];
[ # # ]
85 : :
86 [ # # ][ # # ]: 0 : for ( i = 0; i < poly->tris.size(); i++ ) {
87 : 0 : e0[i] = e1[i] = e2[i]= NO_EDGE_NBR; // initialize
88 [ # # ]: 0 : group.push_back(UNSET);
89 [ # # ]: 0 : v0 = poly->tris[i]->v0;
90 [ # # ]: 0 : v1 = poly->tris[i]->v1;
91 [ # # ]: 0 : v2 = poly->tris[i]->v2;
92 : : // Put the triangle sequence, i, in the hash list for each of the 3 edges.
93 [ # # ]: 0 : hash->addtoHashList((v0+v1)%numhashbins,i);
94 [ # # ]: 0 : hash->addtoHashList((v1+v2)%numhashbins,i);
95 [ # # ]: 0 : hash->addtoHashList((v2+v0)%numhashbins,i);
96 : : }
97 : : /* int jmin, jmax, jave;
98 : : jmin = 100000000;
99 : : jmax = -100000000;
100 : : jave = 0;
101 : : for ( i = 0; i < poly->tris.size(); i++ ) {
102 : : hasharrayptr = hash->getHashBin(i,&hasharraysize);
103 : : if ( hasharraysize < jmin ) jmin = hasharraysize;
104 : : if ( hasharraysize > jmax ) jmax = hasharraysize;
105 : : jave += hasharraysize;
106 : :
107 : : }
108 : : jave /= poly->tris.size();
109 : : printf("jmin jmax jave %d %d %d\n",jmin, jmax, jave);
110 : : */
111 : :
112 [ # # ][ # # ]: 0 : for ( i = 0; i < poly->tris.size(); i++ ) {
113 [ # # ]: 0 : v0 = poly->tris[i]->v0;
114 [ # # ]: 0 : v1 = poly->tris[i]->v1;
115 [ # # ]: 0 : v2 = poly->tris[i]->v2;
116 : 0 : hashvalue = (v0+v1)%numhashbins; // get the hash value for edge 0
117 [ # # ]: 0 : hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize);
118 [ # # ]: 0 : for ( j = 0; j < hasharraysize; j++ ) {
119 : : // Go through the list and find the other triangle, itri, for each edge.
120 : : // Then assign its sequence to the edge array.
121 : 0 : itri = hasharrayptr[j];
122 [ # # ]: 0 : if ( (unsigned int)itri == i ) continue;
123 [ # # ]: 0 : tv0 = poly->tris[itri]->v0;
124 [ # # ]: 0 : tv1 = poly->tris[itri]->v1;
125 [ # # ]: 0 : tv2 = poly->tris[itri]->v2;
126 [ # # ][ # # ]: 0 : if ( ((v1 == tv0) && (v0 == tv1)) || ((v0 == tv0) && (v1 == tv1)) ) {
[ # # ][ # # ]
127 : 0 : e0[i] = itri;
128 [ # # ][ # # ]: 0 : } else if ( ((v1 == tv1) && (v0 == tv2)) || ((v0 == tv1) && (v1 == tv2)) ) {
[ # # ][ # # ]
129 : 0 : e0[i] = itri;
130 [ # # ][ # # ]: 0 : } else if ( ((v1 == tv2) && (v0 == tv0)) || ((v0 == tv2) && (v1 == tv0)) ) {
[ # # ][ # # ]
131 : 0 : e0[i] = itri;
132 : : }
133 : : }
134 : 0 : hashvalue = (v1+v2)%numhashbins;
135 [ # # ]: 0 : hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize);
136 [ # # ]: 0 : for ( j = 0; j < hasharraysize; j++ ) {
137 : 0 : itri = hasharrayptr[j];
138 [ # # ]: 0 : if ( (unsigned int)itri == i ) continue;
139 [ # # ]: 0 : tv0 = poly->tris[itri]->v0;
140 [ # # ]: 0 : tv1 = poly->tris[itri]->v1;
141 [ # # ]: 0 : tv2 = poly->tris[itri]->v2;
142 [ # # ][ # # ]: 0 : if ( ((v1 == tv0) && (v2 == tv1)) || ((v2 == tv0) && (v1 == tv1)) ) {
[ # # ][ # # ]
143 : 0 : e1[i] = itri;
144 [ # # ][ # # ]: 0 : } else if ( ((v1 == tv1) && (v2 == tv2)) || ((v2 == tv1) && (v1 == tv2)) ) {
[ # # ][ # # ]
145 : 0 : e1[i] = itri;
146 [ # # ][ # # ]: 0 : } else if ( ((v1 == tv2) && (v2 == tv0)) || ((v2 == tv2) && (v1 == tv0)) ) {
[ # # ][ # # ]
147 : 0 : e1[i] = itri;
148 : : }
149 : : }
150 : 0 : hashvalue = (v2+v0)%numhashbins;
151 [ # # ]: 0 : hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize);
152 [ # # ]: 0 : for ( j = 0; j < hasharraysize; j++ ) {
153 : 0 : itri = hasharrayptr[j];
154 [ # # ]: 0 : if ( (unsigned int)itri == i ) continue;
155 [ # # ]: 0 : tv0 = poly->tris[itri]->v0;
156 [ # # ]: 0 : tv1 = poly->tris[itri]->v1;
157 [ # # ]: 0 : tv2 = poly->tris[itri]->v2;
158 [ # # ][ # # ]: 0 : if ( ((v0 == tv0) && (v2 == tv1)) || ((v2 == tv0) && (v0 == tv1)) ) {
[ # # ][ # # ]
159 : 0 : e2[i] = itri;
160 [ # # ][ # # ]: 0 : } else if ( ((v0 == tv1) && (v2 == tv2)) || ((v2 == tv1) && (v0 == tv2)) ) {
[ # # ][ # # ]
161 : 0 : e2[i] = itri;
162 [ # # ][ # # ]: 0 : } else if ( ((v0 == tv2) && (v2 == tv0)) || ((v2 == tv2) && (v0 == tv0)) ) {
[ # # ][ # # ]
163 : 0 : e2[i] = itri;
164 : : }
165 : : }
166 : : }
167 : : // Now we have to remove the other-side triangles where there was an intersection
168 : : // edge.
169 [ # # ][ # # ]: 0 : for ( i = 0; i < poly->intersection_edges.size(); i++ ) {
170 [ # # ]: 0 : v0 = poly->intersection_edges[i]->v0;
171 [ # # ]: 0 : v1 = poly->intersection_edges[i]->v1;
172 : 0 : hashvalue = (v0+v1)%numhashbins;
173 [ # # ]: 0 : hasharrayptr = hash->getHashBin(hashvalue,&hasharraysize);
174 [ # # ]: 0 : for ( j = 0; j < hasharraysize; j++ ) {
175 : 0 : itri = hasharrayptr[j];
176 [ # # ]: 0 : tv0 = poly->tris[itri]->v0;
177 [ # # ]: 0 : tv1 = poly->tris[itri]->v1;
178 [ # # ]: 0 : tv2 = poly->tris[itri]->v2;
179 [ # # ][ # # ]: 0 : if ( ((v0 == tv0) && (v1 == tv1)) || ((v0 == tv1) && (v1 == tv0)) ) {
[ # # ][ # # ]
180 : 0 : e0[itri] = NO_EDGE_NBR;
181 : : };
182 [ # # ][ # # ]: 0 : if ( ((v0 == tv0) && (v1 == tv2)) || ((v0 == tv2) && (v1 == tv0)) ) {
[ # # ][ # # ]
183 : 0 : e2[itri] = NO_EDGE_NBR;
184 : : };
185 [ # # ][ # # ]: 0 : if ( ((v0 == tv1) && (v1 == tv2)) || ((v0 == tv2) && (v1 == tv1)) ) {
[ # # ][ # # ]
186 : 0 : e1[itri] = NO_EDGE_NBR;
187 : : }
188 : : }
189 : : }
190 : :
191 : : // Group the triangles that are neighbors.
192 [ # # ][ # # ]: 0 : for ( i = 0; i < poly->tris.size(); i++ ) {
193 [ # # ][ # # ]: 0 : if ( group[i] == UNSET ) {
194 [ # # ]: 0 : fill_group(i,number_of_groups++);
195 : : }
196 : : }
197 : :
198 [ # # ][ # # ]: 0 : delete[] e0; delete[] e1; delete[] e2;
[ # # ]
199 [ # # ][ # # ]: 0 : delete hash;
200 : :
201 : 0 : return CUBIT_SUCCESS;
202 : : }
203 : :
204 : 0 : void FBClassify::fill_group(int itri, int ngroup)
205 : : {
206 [ # # ][ # # ]: 0 : std::stack<int> vstack;
[ # # ]
207 : : int ktri;
208 : :
209 [ # # ]: 0 : group[itri] = ngroup;
210 : :
211 [ # # ][ # # ]: 0 : if ( (e0[itri] != NO_EDGE_NBR) && (group[e0[itri]] == UNSET) )
[ # # ][ # # ]
212 [ # # ]: 0 : vstack.push(e0[itri]);
213 [ # # ][ # # ]: 0 : if ( (e1[itri] != NO_EDGE_NBR) && (group[e1[itri]] == UNSET) )
[ # # ][ # # ]
214 [ # # ]: 0 : vstack.push(e1[itri]);
215 [ # # ][ # # ]: 0 : if ( (e2[itri] != NO_EDGE_NBR) && (group[e2[itri]] == UNSET) )
[ # # ][ # # ]
216 [ # # ]: 0 : vstack.push(e2[itri]);
217 [ # # ][ # # ]: 0 : while (vstack.size() > 0 ) {
218 [ # # ]: 0 : ktri = vstack.top();
219 [ # # ]: 0 : vstack.pop();
220 [ # # ]: 0 : group[ktri] = ngroup;
221 [ # # ][ # # ]: 0 : if ( (e0[ktri] != NO_EDGE_NBR) && (group[e0[ktri]] == UNSET) )
[ # # ][ # # ]
222 [ # # ]: 0 : vstack.push(e0[ktri]);
223 [ # # ][ # # ]: 0 : if ( (e1[ktri] != NO_EDGE_NBR) && (group[e1[ktri]] == UNSET) )
[ # # ][ # # ]
224 [ # # ]: 0 : vstack.push(e1[ktri]);
225 [ # # ][ # # ]: 0 : if ( (e2[ktri] != NO_EDGE_NBR) && (group[e2[ktri]] == UNSET) )
[ # # ][ # # ]
226 [ # # ]: 0 : vstack.push(e2[ktri]);
227 [ # # ]: 0 : }
228 : 0 : }
229 : :
230 : 0 : CubitStatus FBClassify::CharacterizeGroups(int which, bool other_is_planar)
231 : : {
232 : : int numberdone;
233 : : unsigned int i;
234 : : bool ifoundit;
235 : 0 : int itri = -1, type;
236 : : FBPolyhedron *poly;
237 : :
238 [ # # ]: 0 : if ( which == 1 ) poly = polya;
239 : 0 : else poly = polyb;
240 : :
241 [ # # ]: 0 : if ( poly == 0 ) {
242 [ # # ][ # # ]: 0 : PRINT_ERROR("ERROR: Polyhedral object undefined in FBClassify::Group\n");
[ # # ][ # # ]
243 : 0 : return CUBIT_FAILURE;
244 : : }
245 : :
246 : 0 : numberdone = 0;
247 : 0 : i = 0;
248 [ # # ]: 0 : while ( numberdone < number_of_groups ) {
249 : 0 : ifoundit = false;
250 [ # # ][ # # ]: 0 : while ( i < poly->tris.size() ) {
251 [ # # ][ # # ]: 0 : if ( group[i] == numberdone ) {
252 : 0 : ifoundit = true;
253 : 0 : itri = (int)i;
254 : 0 : break;
255 : : }
256 : 0 : i++;
257 : : }
258 [ # # ]: 0 : if ( ifoundit == true ) {
259 [ # # ]: 0 : if ( other_is_planar == false )
260 [ # # ]: 0 : type = classify(itri, which);
261 : : else
262 [ # # ]: 0 : type = classify_against_plane(itri, which);
263 [ # # ]: 0 : group_characterization.push_back(type);
264 : 0 : numberdone++;
265 : : } else {
266 [ # # ][ # # ]: 0 : PRINT_ERROR("Error in FBClassify::CharacterizeGroups\n");
[ # # ][ # # ]
267 : 0 : return CUBIT_FAILURE;
268 : : }
269 : : }
270 : 0 : return CUBIT_SUCCESS;
271 : : }
272 : :
273 : 0 : int FBClassify::classify_against_plane(int itri, int which)
274 : : {
275 : : FBPolyhedron *polyref, *polyobj;
276 : : int type, v0, v1, v2;
277 : : double xbary, ybary, zbary, a, b, c;
278 : : ;
279 : :
280 : 0 : type = FB_ORIENTATION_UNDEFINED;
281 [ # # ]: 0 : if ( which == 1 ) {
282 : 0 : polyobj = polyb;
283 : 0 : polyref = polya;
284 [ # # ]: 0 : } else if ( which == 2 ) {
285 : 0 : polyobj = polya;
286 : 0 : polyref = polya;
287 : : } else {
288 [ # # ][ # # ]: 0 : PRINT_ERROR("ERROR in FBClassify::classify\n");
289 : 0 : return type;
290 : : }
291 : 0 : v0 = polyref->tris[itri]->v0;
292 : 0 : v1 = polyref->tris[itri]->v1;
293 : 0 : v2 = polyref->tris[itri]->v2;
294 : :
295 : 0 : xbary = ( polyref->verts[v0]->coord[0] +
296 : 0 : polyref->verts[v1]->coord[0] +
297 : 0 : polyref->verts[v2]->coord[0] )/3.;
298 : 0 : ybary = ( polyref->verts[v0]->coord[1] +
299 : 0 : polyref->verts[v1]->coord[1] +
300 : 0 : polyref->verts[v2]->coord[1] )/3.;
301 : 0 : zbary = ( polyref->verts[v0]->coord[2] +
302 : 0 : polyref->verts[v1]->coord[2] +
303 : 0 : polyref->verts[v2]->coord[2] )/3.;
304 : 0 : a = polyref->tris[itri]->a;
305 : 0 : b = polyref->tris[itri]->b;
306 : 0 : c = polyref->tris[itri]->c;
307 : :
308 : : // Figure out which side of the plane we are on. Since all
309 : : // of the plane's triangles have the same plane equation
310 : : // coefficients, might as well use the first one.
311 : :
312 : : double obj_tri_a, obj_tri_b, obj_tri_c, obj_tri_d, dotprod, disttoplane;
313 : :
314 : 0 : obj_tri_a = polyobj->tris[0]->a;
315 : 0 : obj_tri_b = polyobj->tris[0]->b;
316 : 0 : obj_tri_c = polyobj->tris[0]->c;
317 : 0 : obj_tri_d = polyobj->tris[0]->d;
318 : :
319 : 0 : disttoplane = obj_tri_a*xbary + obj_tri_b*ybary + obj_tri_c*zbary + obj_tri_d;
320 [ # # ]: 0 : if ( disttoplane > EPSILON ) return FB_ORIENTATION_OUTSIDE;
321 [ # # ]: 0 : else if ( disttoplane < -EPSILON ) return FB_ORIENTATION_INSIDE;
322 : :
323 : 0 : dotprod = obj_tri_a*a + obj_tri_b*b + obj_tri_c*c;
324 [ # # ]: 0 : if ( dotprod > 0. ) return FB_ORIENTATION_SAME;
325 : 0 : else return FB_ORIENTATION_OPPOSITE;
326 : :
327 : : return type;
328 : : }
329 : :
330 : : //returns an orientation for the triangle relative to the other body.
331 : : //This triangle can be inside or outside the other body. Or, it can
332 : : //be opposite or same. Opposite means this triangle is very close to
333 : : // a triangle in the other body and it has an opposite pointing normal.
334 : : // Same means it is very close to a trianglein the other body and it
335 : : // their normals point in the same direction.
336 : 0 : int FBClassify::classify(int itri, int which)
337 : : {
338 : : // "which" is the object, either 1 or 2, that the triangle itri
339 : : // belongs to.
340 : : // Th e object to test for inside or outside is the other object.
341 : : FBPolyhedron *polyref, *polyobj;
342 : : // polyref = itri's polyhedron; polyobj = the other object
343 : : int type, v0, v1, v2;
344 : : double xbary, ybary, zbary, a, b, c;
345 : :
346 : 0 : type = FB_ORIENTATION_UNDEFINED;
347 [ # # ]: 0 : if ( which == 1 ) {
348 : 0 : polyobj = polyb;
349 : 0 : polyref = polya;
350 [ # # ]: 0 : } else if ( which == 2 ) {
351 : 0 : polyobj = polya;
352 : 0 : polyref = polyb;
353 : : } else {
354 [ # # ][ # # ]: 0 : PRINT_ERROR("ERROR in FBClassify::classify\n");
[ # # ][ # # ]
355 : 0 : return type;
356 : : }
357 : 0 : int mydebug = 0;
358 [ # # ]: 0 : if(mydebug)
359 [ # # ][ # # ]: 0 : polyref->debug_draw_fb_triangle(polyref->tris[itri]);
360 : :
361 [ # # ]: 0 : v0 = polyref->tris[itri]->v0;
362 [ # # ]: 0 : v1 = polyref->tris[itri]->v1;
363 [ # # ]: 0 : v2 = polyref->tris[itri]->v2;
364 : :
365 [ # # ]: 0 : xbary = ( polyref->verts[v0]->coord[0] +
366 [ # # ]: 0 : polyref->verts[v1]->coord[0] +
367 [ # # ]: 0 : polyref->verts[v2]->coord[0] )/3.;
368 [ # # ]: 0 : ybary = ( polyref->verts[v0]->coord[1] +
369 [ # # ]: 0 : polyref->verts[v1]->coord[1] +
370 [ # # ]: 0 : polyref->verts[v2]->coord[1] )/3.;
371 [ # # ]: 0 : zbary = ( polyref->verts[v0]->coord[2] +
372 [ # # ]: 0 : polyref->verts[v1]->coord[2] +
373 [ # # ]: 0 : polyref->verts[v2]->coord[2] )/3.;
374 [ # # ]: 0 : a = polyref->tris[itri]->a;
375 [ # # ]: 0 : b = polyref->tris[itri]->b;
376 [ # # ]: 0 : c = polyref->tris[itri]->c;
377 : :
378 : : unsigned int i, num_perturb;
379 : : double obj_tri_a, obj_tri_b, obj_tri_c, obj_tri_d, dotprod;
380 : : double distance_to_other_sqr, closest_distance_to_plane, t;
381 : : double closest_distance_to_other_sqr;
382 : : double xint, yint, zint, distance_to_plane, closest_dotproduct;
383 : : bool perturb, done, foundone;
384 : : double other_xbar, other_ybar, other_zbar;
385 : : int other_tri_0, other_tri_1, other_tri_2;
386 : :
387 : 0 : perturb = false;
388 : 0 : num_perturb = 0;
389 : 0 : done = false;
390 : :
391 [ # # ][ # # ]: 0 : while ( (done == false) && (num_perturb < 20) ) {
392 : 0 : closest_dotproduct = -CUBIT_DBL_MAX + 1.;
393 : 0 : closest_distance_to_plane = CUBIT_DBL_MAX;
394 : 0 : closest_distance_to_other_sqr = CUBIT_DBL_MAX;
395 : 0 : foundone = false;
396 [ # # ][ # # ]: 0 : for ( i = 0; i < polyobj->tris.size(); i++ ) {
397 [ # # ]: 0 : obj_tri_a = polyobj->tris[i]->a;
398 [ # # ]: 0 : obj_tri_b = polyobj->tris[i]->b;
399 [ # # ]: 0 : obj_tri_c = polyobj->tris[i]->c;
400 [ # # ]: 0 : obj_tri_d = polyobj->tris[i]->d;
401 : :
402 : 0 : dotprod = obj_tri_a*a + obj_tri_b*b + obj_tri_c*c;
403 : : //calculate the distance to the other triangles plane
404 : 0 : distance_to_plane = (obj_tri_a*xbary + obj_tri_b*ybary +
405 : 0 : obj_tri_c*zbary + obj_tri_d);
406 : :
407 : :
408 [ # # ]: 0 : if ( fabs(dotprod) < EPSILON_CLASSIFY ) {
409 : : // Is the point in the plane?
410 [ # # ]: 0 : if ( fabs(distance_to_plane) < EPSILON_CLASSIFY ) {
411 : : // Perturb the ray and recast.
412 : 0 : perturb = true;
413 : 0 : num_perturb += 1;
414 : 0 : break;
415 : : }
416 : 0 : continue;
417 : : }
418 : :
419 : 0 : t =-(distance_to_plane)/dotprod;
420 [ # # ]: 0 : if ( t < -EPSILON_CLASSIFY ) continue;
421 : 0 : xint = xbary + a*t;
422 : 0 : yint = ybary + b*t;
423 : 0 : zint = zbary + c*t;
424 : :
425 : : // Check whether the intersection point lies in or on
426 : : // the object triangle's
427 : : // bounding box.
428 [ # # ][ # # ]: 0 : if ( (polyobj->tris[i]->boundingbox.xmin - EPSILON > xint) ||
[ # # ]
429 [ # # ][ # # ]: 0 : (polyobj->tris[i]->boundingbox.xmax + EPSILON < xint) ||
430 [ # # ][ # # ]: 0 : (polyobj->tris[i]->boundingbox.ymin - EPSILON > yint) ||
431 [ # # ][ # # ]: 0 : (polyobj->tris[i]->boundingbox.ymax + EPSILON < yint) ||
432 [ # # ][ # # ]: 0 : (polyobj->tris[i]->boundingbox.zmin - EPSILON > zint) ||
[ # # ]
433 [ # # ]: 0 : (polyobj->tris[i]->boundingbox.zmax + EPSILON < zint) )
434 : 0 : continue;
435 : :
436 : : // Is the point (xint, yint, zint) inside or on the triangle?
437 : : // Get a principal projection to make this a 2D problem.
438 : : double xp1, yp1, xp2, yp2, xp3, yp3, ptx, pty;
439 : : int retval;
440 : :
441 [ # # ][ # # ]: 0 : if ( (fabs(obj_tri_b) >= fabs(obj_tri_a)) &&
442 : 0 : (fabs(obj_tri_b) >= fabs(obj_tri_c)) ) {
443 [ # # ][ # # ]: 0 : xp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[0];
444 [ # # ][ # # ]: 0 : yp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[2];
445 [ # # ][ # # ]: 0 : xp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[0];
446 [ # # ][ # # ]: 0 : yp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[2];
447 [ # # ][ # # ]: 0 : xp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[0];
448 [ # # ][ # # ]: 0 : yp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[2];
449 : 0 : ptx = xint;
450 : 0 : pty = zint;
451 [ # # ]: 0 : } else if ( fabs(obj_tri_a) >= fabs(obj_tri_c) ) {
452 [ # # ][ # # ]: 0 : xp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[1];
453 [ # # ][ # # ]: 0 : yp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[2];
454 [ # # ][ # # ]: 0 : xp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[1];
455 [ # # ][ # # ]: 0 : yp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[2];
456 [ # # ][ # # ]: 0 : xp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[1];
457 [ # # ][ # # ]: 0 : yp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[2];
458 : 0 : ptx = yint;
459 : 0 : pty = zint;
460 : : } else {
461 [ # # ][ # # ]: 0 : xp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[0];
462 [ # # ][ # # ]: 0 : yp1 = polyobj->verts[polyobj->tris[i]->v0]->coord[1];
463 [ # # ][ # # ]: 0 : xp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[0];
464 [ # # ][ # # ]: 0 : yp2 = polyobj->verts[polyobj->tris[i]->v1]->coord[1];
465 [ # # ][ # # ]: 0 : xp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[0];
466 [ # # ][ # # ]: 0 : yp3 = polyobj->verts[polyobj->tris[i]->v2]->coord[1];
467 : 0 : ptx = xint;
468 : 0 : pty = yint;
469 : : }
470 [ # # ]: 0 : retval = pt_in_tri_2d(ptx,pty,xp1,yp1,xp2,yp2,xp3,yp3);
471 [ # # ][ # # ]: 0 : if ( (retval == FB_ORIENTATION_INSIDE) ||
472 : : (retval == FB_ORIENTATION_ON) ) {
473 : : //calculate the distance to the other triangle's centroid
474 [ # # ]: 0 : other_tri_0 = polyobj->tris[i]->v0;
475 [ # # ]: 0 : other_tri_1 = polyobj->tris[i]->v1;
476 [ # # ]: 0 : other_tri_2 = polyobj->tris[i]->v2;
477 [ # # ]: 0 : other_xbar = ( polyobj->verts[other_tri_0]->coord[0] +
478 [ # # ]: 0 : polyobj->verts[other_tri_1]->coord[0] +
479 [ # # ]: 0 : polyobj->verts[other_tri_2]->coord[0] )/3.;
480 [ # # ]: 0 : other_ybar = ( polyobj->verts[other_tri_0]->coord[1] +
481 [ # # ]: 0 : polyobj->verts[other_tri_1]->coord[1] +
482 [ # # ]: 0 : polyobj->verts[other_tri_2]->coord[1] )/3.;
483 [ # # ]: 0 : other_zbar = ( polyobj->verts[other_tri_0]->coord[2] +
484 [ # # ]: 0 : polyobj->verts[other_tri_1]->coord[2] +
485 [ # # ]: 0 : polyobj->verts[other_tri_2]->coord[2] )/3.;
486 : :
487 : : //calculate the distance (squared) to the other triangle's centroid
488 : 0 : distance_to_other_sqr = ( (xbary-other_xbar)*(xbary-other_xbar) +
489 : 0 : (ybary-other_ybar)*(ybary-other_ybar) +
490 : 0 : (zbary-other_zbar)*(zbary-other_zbar) );
491 : : //if this is the closest other triangle so far...
492 [ # # ]: 0 : if(closest_distance_to_other_sqr > distance_to_other_sqr){
493 : : //then we found one, and update the closest distance, dot prod,
494 : : // and distance to other plane.
495 : 0 : foundone = true;
496 : 0 : closest_distance_to_other_sqr = distance_to_other_sqr;
497 : 0 : closest_dotproduct = dotprod;
498 [ # # ]: 0 : if(mydebug){
499 [ # # ][ # # ]: 0 : polyobj->debug_draw_fb_triangle(polyobj->tris[i]);
500 [ # # ]: 0 : GfxDebug::mouse_xforms();
501 : : }
502 : 0 : closest_distance_to_plane = distance_to_plane;
503 : : // This is the closest triangle.
504 [ # # ]: 0 : if ( fabs(closest_distance_to_plane) < EPSILON_CLASSIFY )
505 : 0 : break;
506 : : }
507 : : }
508 : : }
509 [ # # ]: 0 : if ( perturb == false ) done = true;
510 : : else {
511 : : // perturb the ray and try again.
512 [ # # ]: 0 : perturb_the_ray(xbary, ybary, zbary);
513 : 0 : perturb = false;
514 : : }
515 : : }
516 : : //if we are very close to the plane of the closest other triangle
517 : : // then we are either going to classify as opposite or same depending
518 : : // on the relationship of the normals
519 [ # # ]: 0 : if ( (fabs(closest_distance_to_plane) < EPSILON_CLASSIFY ) ){
520 [ # # ]: 0 : if ( closest_dotproduct > 0.0 )
521 : 0 : type = FB_ORIENTATION_SAME;
522 : : else
523 : 0 : type = FB_ORIENTATION_OPPOSITE;
524 : : }
525 : : //otherwise we are going to classify as inside or outside. If we
526 : : // didn't find any triangles that projected into our triangle,
527 : : // we must be outside. Otherwise we compare the normals to determine
528 : : // whether we are inside or outside.
529 : : else {
530 [ # # ]: 0 : if ( foundone == false )
531 : 0 : type = FB_ORIENTATION_OUTSIDE;
532 [ # # ]: 0 : else if ( closest_dotproduct > 0.0 )
533 : 0 : type = FB_ORIENTATION_INSIDE;
534 : : else
535 : 0 : type = FB_ORIENTATION_OUTSIDE;
536 : : }
537 [ # # ]: 0 : if(mydebug)
538 [ # # ]: 0 : GfxDebug::display_all();
539 : 0 : return type;
540 : : }
541 : :
542 : 0 : void FBClassify::perturb_the_ray(double &xbary, double &ybary, double &zbary)
543 : : {
544 : 0 : xbary += 1.e-4*(double(rand())/(RAND_MAX+1.0)-0.5);
545 : 0 : ybary += 1.e-4*(double(rand())/(RAND_MAX+1.0)-0.5);
546 : 0 : zbary += 1.e-4*(double(rand())/(RAND_MAX+1.0)-0.5);
547 : 0 : }
548 : :
549 : 0 : int FBClassify::pt_in_tri_2d(double xpt, double ypt,
550 : : double x0, double y0,
551 : : double x1, double y1,
552 : : double x2, double y2)
553 : : {
554 : : // From Schneider & Eberly, "Geometric Tools for COmputer Graphics",
555 : : // Chap. 13.3.1. If triangle is needle-thin, CUBIT_FAILURE might be
556 : : // returned, in wich case is_point_in is undefined.
557 : :
558 : : double c0, c1, c2;
559 : : double e0x, e1x, e2x, e0y, e1y, e2y;
560 : : double n0x, n1x, n2x, n0y, n1y, n2y;
561 : : double denom0, denom1, denom2;
562 : : int result;
563 : :
564 : 0 : e0x = x1 - x0; e0y = y1 - y0;
565 : 0 : e1x = x2 - x1; e1y = y2 - y1;
566 : 0 : e2x = x0 - x2; e2y = y0 - y2;
567 : 0 : n0x = e0y; n0y = -e0x;
568 : 0 : n1x = e1y; n1y = -e1x;
569 : 0 : n2x = e2y; n2y = -e2x;
570 : 0 : denom0 = n1x*e0x + n1y*e0y;
571 [ # # ]: 0 : if ( fabs(denom0) < EPSILON_CLASSIFY ) {
572 [ # # ][ # # ]: 0 : PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n");
573 : 0 : return FB_ORIENTATION_UNDEFINED;
574 : : }
575 : 0 : denom1 = n2x*e1x + n2y*e1y;
576 [ # # ]: 0 : if ( fabs(denom1) < EPSILON_CLASSIFY ) {
577 [ # # ][ # # ]: 0 : PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n");
578 : 0 : return FB_ORIENTATION_UNDEFINED;
579 : : }
580 : 0 : denom2 = n0x*e2x + n0y*e2y;
581 [ # # ]: 0 : if ( fabs(denom2) < EPSILON_CLASSIFY ) {
582 [ # # ][ # # ]: 0 : PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n");
583 : 0 : return FB_ORIENTATION_UNDEFINED;
584 : : }
585 : :
586 : 0 : c0 = -( n1x*(xpt-x1) + n1y*(ypt-y1) )/denom0;
587 : 0 : c1 = -( n2x*(xpt-x2) + n2y*(ypt-y2) )/denom1;
588 : 0 : c2 = -( n0x*(xpt-x0) + n0y*(ypt-y0) )/denom2;
589 : :
590 [ # # ][ # # ]: 0 : if ( (c0 > 0.0) && (c1 > 0.0) && (c2 > 0.0) ) result = FB_ORIENTATION_INSIDE;
[ # # ]
591 [ # # ][ # # ]: 0 : else if ( (c0 < 0.0) || (c1 < 0.0) || (c2 < 0.0) ) result = FB_ORIENTATION_OUTSIDE;
[ # # ]
592 : 0 : else result = FB_ORIENTATION_ON;
593 : :
594 : 0 : return result;
595 : :
596 : : }
597 : :
598 : 0 : void FBClassify::get_group(std::vector<int> **this_group,
599 : : std::vector<int> **this_group_characterization)
600 : : {
601 : 0 : *this_group = &group;
602 : 0 : *this_group_characterization = &group_characterization;
603 [ + - ][ + - ]: 6540 : }
|
