//=========================================================================== // GoTools Core - SINTEF Geometry Tools Core library, version 2.0.1 // // Copyright (C) 2000-2007, 2010 SINTEF ICT, Applied Mathematics, Norway. // // This program is free software; you can redistribute it and/or // modify it under the terms of the GNU General Public License // as published by the Free Software Foundation version 2 of the License. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program; if not, write to the Free Software // Foundation, Inc., // 59 Temple Place - Suite 330, // Boston, MA 02111-1307, USA. // // Contact information: E-mail: tor.dokken@sintef.no // SINTEF ICT, Department of Applied Mathematics, // P.O. Box 124 Blindern, // 0314 Oslo, Norway. // // Other licenses are also available for this software, notably licenses // for: // - Building commercial software. // - Building software whose source code you wish to keep private. //=========================================================================== 00106 #ifndef DBG 00107 const bool dbg = false; 00108 #endif 00109 00110 if (!degenerate_triangle(vert_p[c1], vert_p[c2], vert_p[c3] DBG_FLAG)) 00111 { 00112 mesh.push_back(c1); 00113 mesh.push_back(c2); 00114 mesh.push_back(c3); 00115 } 00116 else 00117 if (dbg) 00118 { 00119 printf(" add_triangle(): Not adding triangle. Degenerate.\n"); 00120 printf("\n"); 00121 printf(" hold on; plot(%f, %f, 'ro', 'markersize', 10, 'linewidth', 2);\n", vert_p[c1][0], vert_p[c1][1]); 00122 printf(" plot(%f, %f, 'go', 'markersize', 10, 'linewidth', 2);\n", vert_p[c2][0], vert_p[c2][1]); 00123 printf(" plot(%f, %f, 'bo', 'markersize', 10, 'linewidth', 2); hold off\n", vert_p[c3][0], vert_p[c3][1]); 00124 printf("\n"); 00125 // const bool tmp = degenerate_triangle(vert_p[c1], vert_p[c2], vert_p[c3] DBG_FLAG); 00126 } 00127 } 00128 00129 00130 00131 00132 00133 00134 #ifdef DBG 00135 00136 //============================================================================================================== 00137 // 00138 // 090216: For debugging. 00139 // 00140 //============================================================================================================== 00141 00142 bool same_point(const Vector2D &c1, const Vector2D &d1, const double tol) 00143 { 00144 return (fabs(c1[0]-d1[0])<tol) && (fabs(c1[1]-d1[1])<tol); 00145 } 00146 00147 bool same_triangle(const Vector2D &c1, const Vector2D &c2, const Vector2D &c3, 00148 const Vector2D &d1, const Vector2D &d2, const Vector2D &d3, 00149 const double tol) 00150 { 00151 return ( ((same_point(c1, d1, tol)) && (same_point(c2, d2, tol)) && (same_point(c3, d3, tol))) || 00152 ((same_point(c1, d1, tol)) && (same_point(c2, d3, tol)) && (same_point(c3, d2, tol))) || 00153 ((same_point(c1, d2, tol)) && (same_point(c2, d1, tol)) && (same_point(c3, d3, tol))) || 00154 ((same_point(c1, d2, tol)) && (same_point(c2, d3, tol)) && (same_point(c3, d1, tol))) || 00155 ((same_point(c1, d3, tol)) && (same_point(c2, d1, tol)) && (same_point(c3, d2, tol))) || 00156 ((same_point(c1, d3, tol)) && (same_point(c2, d2, tol)) && (same_point(c3, d1, tol))) ); 00157 } 00158 00159 #endif 00160 00161 00162 00163 00164 00165 00166 //-------------------------------------------------------------------------------------------------------------- 00167 // 00168 // This subroutine assumes that either one, two or three points are inside, and the rest outside, so that two 00169 // intersections (between "grid lines" and trim curve) are to be found, and either one, two or three triangles 00170 // are to be added, respectively. 00171 // 00172 // 081208: Note that 'srf' is not evaluated to create new vertices, the old ones are simply interpolated! This 00173 // could be improved! 00174 // 00175 // 081208: Hmm... Can this interpolation result in non-normalized normals?! Think so! But is it specified that 00176 // they should be normalized?! Safest thing would probably be to normalize them... Not doing it for 00177 // split_quad, but for split_triangle... 00178 // 00179 // 081208: Adding 'vert_p' for use by later 'split_triangle'-calls. This will hold parameter pairs. 00180 // 00181 // 090117: Adding 'with_boundary' so that the caller may decide whether or not the boundary (with an 00182 // "epsilon-tube" around it) should be defined as "inside". 00183 // 00184 //-------------------------------------------------------------------------------------------------------------- 00185 00186 void split_quad(const int i0, const int j0, 00187 const int i1, const int j1, // May not be used, if so, -1 is used. 00188 const int i2, const int j2, // May not be used, if so, -1 is used. 00189 const vector<int> &si0, const vector<int> &sj0, 00190 const vector<int> &si1, const vector<int> &sj1, 00191 00192 // 090115: All the vectors above have two elements each, s?0[i] and s?1[i] are the two end 00193 // points of two of the sides of some quad, that the contour is assumed to intersect. 00194 // (The contour is given in 'contour', and 'vert', together. I think.) 00195 // (i?, j?) (actually, it makes more sense to say (j?, i?)...) specifies corners of 00196 // the quad. Which corners, is dependent on the configuration in question... (sigh.) 00197 00198 boost::shared_ptr<SplineSurface> srf, 00199 vector< Vector3D > &vert, vector<Vector2D> &vert_p, 00200 vector< int > &bd, vector< Vector3D > &norm, 00201 //vector< Vector3D > &col, 00202 vector<int> &mesh, const int dn, const int dm, 00203 //vector< Vector3D > &extra_v, 00204 const vector< Vector3D > &trim_curve_p, 00205 const vector<int> &contour, 00206 const bool with_boundary, const bool dbg=false) 00207 { 00208 // 00209 // This one always get two existing points and two segments between two pairs of points in, and intersect 00210 // the trimming curve twice against these segments for the two remaining points. The triangles are not 00211 // always oriented the same way, but that shouldn't matter too much. [Could easily be fixed by a switch 00212 // among the arguments...] 00213 // 00214 00215 int i; 00216 00217 // 00218 // Do the intersections. 00219 // 00220 00221 const double u0=srf->startparam_u(); 00222 const double u1=srf->endparam_u(); 00223 const double v0=srf->startparam_v(); 00224 const double v1=srf->endparam_v(); 00225 for (i=0; i<2; i++) 00226 { 00227 double x, y, s; 00228 if (dbg) 00229 printf("%f %f %f %f\n", u0*(1.0-sj0[i]/double(dn)) + u1*sj0[i]/double(dn), 00230 v0*(1.0-si0[i]/double(dm)) + v1*si0[i]/double(dm), 00231 u0*(1.0-sj1[i]/double(dn)) + u1*sj1[i]/double(dn), 00232 v0*(1.0-si1[i]/double(dm)) + v1*si1[i]/double(dm)); 00233 bool inters_found = segment_contour_intersection_for_s2m(u0*(1.0-sj0[i]/double(dn)) + u1*sj0[i]/double(dn), 00234 v0*(1.0-si0[i]/double(dm)) + v1*si0[i]/double(dm), 00235 u0*(1.0-sj1[i]/double(dn)) + u1*sj1[i]/double(dn), 00236 v0*(1.0-si1[i]/double(dm)) + v1*si1[i]/double(dm), 00237 &trim_curve_p[0][0], contour, 00238 x, y, s, 00239 false // 100222: New 'snap_ends' flag. 00240 DBG_FLAG); 00241 if (dbg) 00242 printf("=1===> i=%d, inters=%d\n", i, inters_found?1:0); 00243 00244 if ((!inters_found) && (with_boundary)) 00245 { 00246 // 090117: Didn't find intersection, even though one should be here. Checking if the segment's ends 00247 // are on the contour. 00248 if (is_on_contour(trim_curve_p, contour, 00249 u0*(1.0-sj0[i]/double(dn)) + u1*sj0[i]/double(dn), 00250 v0*(1.0-si0[i]/double(dm)) + v1*si0[i]/double(dm))) 00251 { 00252 x=u0*(1.0-sj0[i]/double(dn)) + u1*sj0[i]/double(dn); 00253 y=v0*(1.0-si0[i]/double(dm)) + v1*si0[i]/double(dm); 00254 s=0.0; 00255 inters_found=true; 00256 if (dbg) 00257 printf("=2===> i=%d, inters=%d\n", i, inters_found?1:0); 00258 } 00259 else 00260 if (is_on_contour(trim_curve_p, contour, 00261 u0*(1.0-sj1[i]/double(dn)) + u1*sj1[i]/double(dn), 00262 v0*(1.0-si1[i]/double(dm)) + v1*si1[i]/double(dm))) 00263 { 00264 x=u0*(1.0-sj1[i]/double(dn)) + u1*sj1[i]/double(dn); 00265 y=v0*(1.0-si1[i]/double(dm)) + v1*si1[i]/double(dm); 00266 s=1.0; 00267 inters_found=true; 00268 if (dbg) 00269 printf("=3===> i=%d, inters=%d\n", i, inters_found?1:0); 00270 } 00271 } 00272 00273 if (!inters_found) 00274 { 00275 MESSAGE("Failed finding intersection. Expected one."); 00276 fflush(stderr); 00277 } 00278 else 00279 { 00280 const double eps=1e-13; // 090115: Used for zero-tests for distances in the parameter domain. 00281 // Hmm... these should *really*, *really* be taken from some global 00282 // variable or something 00283 // 090117: Anyway... s will never be out of range here... 00284 ASSERT2((s>=-eps) && (s<=1.0+eps), printf("Huh?! s=%f\n", s)); 00285 } 00286 // eval surf i x og y for aa faa vert og norm, evt. returnere param 00287 // for skjaering mellom punktene s.a. vi igjen kan danne lin komb. 00288 if (inters_found) 00289 { 00290 vert.push_back((1.0-s)*vert[si0[i]*(dn+1)+sj0[i]] + 00291 s*vert[si1[i]*(dn+1)+sj1[i]] ); 00292 vert_p.push_back((1.0-s)*vert_p[si0[i]*(dn+1)+sj0[i]] + 00293 s*vert_p[si1[i]*(dn+1)+sj1[i]] ); 00294 if (dbg) 00295 printf("=4===> inters: %f %f\n", vert_p[vert_p.size()-1][0], vert_p[vert_p.size()-1][1]); 00296 bd.push_back(1); 00297 Vector3D new_norm = (1.0-s)*norm[si0[i]*(dn+1)+sj0[i]] + 00298 s*norm[si1[i]*(dn+1)+sj1[i]]; 00299 new_norm.normalize(); 00300 norm.push_back(new_norm); 00301 } 00302 //extra_v.push_back(vert[vert.size()-1]); // for debugging 00303 } 00304 00305 // 00306 // Add triangles to the mesh. 00307 // 00308 // Adding three triangles: (A, B, 0), (B, 1, 0), (0, 2, A). 00309 // 00310 // (i1, j1) B 00311 // o------x----o 00312 // | \ | 00313 // | \ x A 00314 // | | 00315 // o-----------o 00316 // (i0, j0) (i2, j2) 00317 // 00318 // 00319 // Adding two triangles: 00320 // 00321 // "-1" -----o "-2" Note that we want the triangles 00322 // o-----/ | oriented ccw., i.e., 00323 // | | (-2, -1, 0), and (-2, 0, 1). 00324 // | | 00325 // o-----------o 00326 // (i0, j0) (i1, j1) 00327 // 00328 00329 // This triangle can be used in all three modes. 00330 add_triangle(vert_p, mesh, vert.size()-2, vert.size()-1, i0*(dn+1)+j0); 00331 //col.push_back(Vector3D(1.0, 0.6, 0.6)); 00332 00333 if ((i1==-1) || (j1==-1)) 00334 // Only one triangle to add. 00335 return; 00336 00337 if ((i2==-1) || (j2==-1)) 00338 { 00339 // Two triangles to add, i.e., one more. 00340 add_triangle(vert_p, mesh, vert.size()-2, i0*(dn+1)+j0, i1*(dn+1)+j1); 00341 //col.push_back(Vector3D(0.0, 0.0, 1.0)); 00342 } 00343 else 00344 { 00345 // Three triangles to add, i.e., two more. 00346 add_triangle(vert_p, mesh, i0*(dn+1)+j0, vert.size()-1, i1*(dn+1)+j1); 00347 //col.push_back(Vector3D(0.6, 1.0, 0.6)); 00348 00349 add_triangle(vert_p, mesh, vert.size()-2, i0*(dn+1)+j0, i2*(dn+1)+j2); 00350 //col.push_back(Vector3D(0.6, 0.6, 1.0)); 00351 } 00352 } 00353 00354 00355 00356 00357 00358 00359 //============================================================================================================== 00360 // 00361 // 081208: Splitting triangles. The original split_quad should probably have been implemented by two calls to 00362 // such a routine as this, it would have made for smaller and more compact code. 00363 // 00364 // This subroutine assumes that either one or two points are inside, and the rest outside, so that 00365 // exactly two intersections (between triangle edges and trim curve) are to be found, and either one 00366 // or two triangles are to be added. 00367 // 00368 // (If there are 2n, for n>1, edge-curve intersections, the triangle should be split recursively, or 00369 // the a priori refinement was too coarse...) 00370 // 00371 // 100218: The first corner, c1, is assumed to be inside, and we may then assume that one of three 00372 // possible situations has occured: c2 inside and c3 outside, or c2 outside and c3 inside, 00373 // or both c2 and c3 outside. (If all are inside, we shouldn't be here in the first place.) 00374 // 00375 // Since at least one corner is inside and one is outside, per assumption, we may also 00376 // assume now (this is made sure by the caller, who may re-order corners) that c1 is inside 00377 // and c2 is outside. This leaves just two cases: 00378 // 00379 // 1) c1 inside, c2 outside and c3 inside, or 00380 // 2) c1 inside, c2 outside and c3 outside. 00381 // 00382 // There should then be exactly two intersections, one on 1-2, and one on either 2-3 or on 00383 // 1-3. 00384 // 00385 // 081208: Note that 'srf' is not evaluated to create new vertices, the old ones are simply interpolated! This 00386 // could be improved! 00387 // 00388 // 081208: Hmm... Can this interpolation result in non-normalized normals?! Think so! But is it specified that 00389 // they should be normalized?! Safest thing would probably be to normalize them... Not doing it for 00390 // split_quad, but for split_tirangle... 00391 // 00392 // 090129: Fixed bug; the routine did not update 'vert_p', hence buggy later splitting of the newly produced 00393 // triangles. 00394 // 00395 // 090203: Hmm... There is a problem here, when the trimming curve passes through one or more of the corners 00396 // of the triangle. How can this be solved? 00397 // First idea: Simply disallow the second intersection to be the same as the first. 00398 // One advantage of this, if possible, is that 'segment_contour_intersection_for_s2m' will 00399 // not have to be modified. Such a modification could affect everything else in a somewhat 00400 // unpredictable manner... 00401 // 00402 // 090217: Adding 'forced_skipping_of_second_edge' so that we get more control. This makes sure some triangles 00403 // with more than two intersections are correctly (i.e., better) handled. Unfortunately, it doesn't 00404 // solve all problems... 00405 // 00406 // 090218: Return value introduced. Will be set to true if further processing of produced triangles is needed. 00407 // 00408 // 100218: There is a problem with cases where the contour intersects (or touches) the triangles in more 00409 // than two places. This may for instance happen when the contour has a "thin corner" like this: 00410 // 00411 // \ | 00412 // c2 +---\---------o c3 (c3=) i2 (2-3: s2=1) 00413 // \ \ /| 00414 // \ \ / | 00415 // \ \ / | 00416 // \\ / | 00417 // X | 00418 // c1 \ | 00419 // (c1=) i1 (1-2: s1=0) \ | 00420 // \ | 00421 // \| 00422 // 00423 // Here, c2 is outside, c1 and c3 inside (really, "on", but that info is not passed on to 00424 // 'split_triangle' as of today) and we get intersections i1 and i2 on the edges 1-2 and 2-3, 00425 // respectively. This will result in the triangles (c1, i1, i2) and (c1, i2, c3), but 00426 // 00427 // 100222: Another problem: Previously to getting here, a point may have been classified as "inside" due to 00428 // in fact being "on" the contour, and when we later get here for splitting, the point was suddenly 00429 // not being classified as "on" anymore because the contour segment which it is "on" is parallel to 00430 // the ray being used in the "crossings loop" in the inside-test. To fix this, we must (can) snap 00431 // points to ends of (such) segments. This is now done by sending the appropriate flag ('snap_ends') 00432 // to 'segment_contour_intersection_for_s2m'. Currently, this routine is the only one setting that 00433 // flag to true. 00434 // 00435 // 100223: Adding a flag for "innerness" of contours, because this is needed in order to take the 00436 // appropriate action when the first intersection ('first_inters') is not found, see below. 00437 // 00438 // 100223: New situation not anticipated arises for 'bootm_face.g2'... c1 inside, c2 and c3 outside, none 00439 // "on", but there are four intersections on the edges, and nothing really goes wrong, we just get 00440 // an incorrect result due to this unforeseen and complex situation: 00441 // 00442 // 00443 // | 1 /contour 00444 // | /| / 00445 // | / | / 00446 // | / | / 00447 // | / T |/ 00448 // i4 | /....... / i1 00449 // | /| 00450 // / | Q / | 00451 // / | / | 00452 // 3 ----|------/-- 2 00453 // i3 | / i2 00454 // | / 00455 // | / 00456 // | / 00457 // | / 00458 // |/ 00459 // 00460 // The result will be the triangle T, and we miss the quad Q, which could have been composed of two 00461 // triangles. This is just one possible hard configuration. An implicit assumption of the current 00462 // approach is that the contours cannot intersect triangles in more than at most two places. Of 00463 // course, a real life contour might intersect a triangle any number of times... 00464 // 00465 // I have a feeling that 'forced_skipping_of_second_edge' possibly should be removed, and maybe 00466 // replaced with some recursive refinement variation for cases where more than the expected maximal 00467 // two intersections can be detected... Adding an experimental switch for this. 00468 // 00469 // 100223: In order to try not to break too much of the old code, and to improve readability, we'll try to 00470 // resolve the situation with another processing stage: Before entering the main processing of the 00471 // "old" 'split_triangle', we will check if we can find three distinct intersections of the edges of 00472 // the triangle. If so, we split the triangle first, then call the "old" routine for each of these 00473 // triangles. The simplest way to do this is to just split and return 'true' (meaning that further 00474 // processing must be done) and then let 'trim_a_triangle' take care of the rest! 00475 // 00476 // 100224: It is really, really time to get rid of 'forced_skipping_of_second_edge'. In a very few cases, it 00477 // may still produce better results, even though most of the time it is now (with the new 00478 // three-intersecting-edges-preprocessing stage) obsolete... If those few cases could be solved by 00479 // other means, it would make for a nice cleanup in the removal of 00480 // 'forced_skipping_of_second_edge'... 00481 // 00482 // 100224: As an attempt to fix the above mentioned situation, and possibly many similar ones, a new 00483 // fallback for triangles of case 's==0' in 'trim_a_triangle' will be added. It is also an advantage 00484 // in not cluttering up the current 'split_triangle', which is complex enough as it is. 00485 // 00486 //============================================================================================================== 00487 00488 00489 00490 00491 //============================================================================================================== 00492 // 00493 // 100223: For increased readability... 00494 // Not using references, in fear of what push_backing on the vertex lists will do... 00495 // (Hmm, should be safe, I think...) 00496 // 00497 // 100224: Hmm... It is probably not necessary to do this linear interpolation with 's', as 00498 // 'segment_contour_intersection_for_s2m' seems to return the same result in its 'x' and 'y' 00499 // parameters... Remember to fix this after the stuff works as it should... 00500 // 00501 // 100225: Now evaluating the surface instead of interpolating triangle corners in 3D. 00502 // 00503 //============================================================================================================== 00504 00505 inline void push_an_intersection(const boost::shared_ptr<SplineSurface> srf, 00506 const Vector3D &c1, const Vector3D &c2, 00507 const Vector2D &c1_p, const Vector2D &c2_p, 00508 const Vector3D &n1, const Vector3D &n2, 00509 const double &s, 00510 vector<Vector3D> &vert, vector<Vector2D> &vert_p, 00511 vector<Vector3D> &norm, 00512 vector<int> &bd) 00513 { 00514 vert_p.push_back( (1.0-s)*c1_p + s*c2_p ); 00515 00516 #ifndef EVAL_SRF_NOT_INTERP 00517 vert .push_back( (1.0-s)*c1 + s*c2 ); 00518 Vector3D tmp = (1.0-s)*n1 + s*n2; 00519 tmp.normalize(); 00520 norm.push_back(tmp); 00521 bd.push_back(1); 00522 #else 00523 vector<Point> res(3); 00524 srf->point(res, vert_p.back()[0], vert_p.back()[1], 1); // Could we have used 0 here? 00525 00526 // The point on the surface: 00527 vert.push_back( Vector3D(res[0].begin()) ); 00528 00529 // The normal in the same point: 00530 const Point nrm = res[1].cross(res[2]); 00531 Vector3D tmp = Vector3D(nrm); 00532 tmp.normalize(); 00533 norm.push_back(tmp); 00534 00535 // This is Vibeke's, don't know what it's for... 00536 bd.push_back(1); 00537 #endif 00538 00539 } 00540 00541 00542 00543 00544 bool split_triangle(const int c1_indx, const int c2_indx, const int c3_indx, 00545 boost::shared_ptr<SplineSurface> srf, 00546 vector< Vector3D > &vert, vector<Vector2D> &vert_p, 00547 vector< int > &bd, vector< Vector3D > &norm, 00548 vector<int> &newmesh, 00549 const vector< Vector3D > &trim_curve_p, 00550 const vector<int> &contour, 00551 const bool dbg = false, 00552 /* const */ bool forced_skipping_of_second_edge = false, 00553 const bool inner_trimming_curve = false) // 100223: See above. 00554 { 00555 const double abs_eps = 1e-12; // 100218: Used for distance between two intersections in the param domain. 00556 00557 // 100223: Testing a new approach. Want to get rid of 'forced_skipping_of_second_edge', but trying to do 00558 // it one step at a time, through these temporary switches... 00559 00560 const bool try_without_forced_stuff = true; 00561 00562 if (try_without_forced_stuff) 00563 forced_skipping_of_second_edge = false; 00564 00565 // 090129: Not using references in fear of what 'push_back' will do... 00566 const Vector2D c1_p=vert_p[c1_indx], c2_p=vert_p[c2_indx], c3_p=vert_p[c3_indx]; 00567 const Vector3D c1=vert[c1_indx], c2=vert[c2_indx], c3=vert[c3_indx]; 00568 const Vector3D n1=norm[c1_indx], n2=norm[c2_indx], n3=norm[c3_indx]; 00569 00570 if (dbg) 00571 { 00572 puts("\n\n ----- split_triangle starting -----------------------------------------------------------------"); 00573 printf("\n The corners:\n\n"); 00574 printf("hold on; plot(%f, %f, 'rd', 'markersize', %d, 'linewidth', 2); hold off\n", c1_p[0], c1_p[1], 14); 00575 printf("hold on; plot(%f, %f, 'gd', 'markersize', %d, 'linewidth', 2); hold off\n", c2_p[0], c2_p[1], 14); 00576 printf("hold on; plot(%f, %f, 'bd', 'markersize', %d, 'linewidth', 2); hold off\n\n", c3_p[0], c3_p[1], 14); 00577 } 00578 00579 double x1, y1, s1, x2, y2, s2; 00580 00581 if (dbg) printf(" Looking for intersection on edge 1-2...\n"); 00582 const bool first_inters = 00583 segment_contour_intersection_for_s2m(c1_p[0], c1_p[1], c2_p[0], c2_p[1], 00584 &trim_curve_p[0][0], // Pointer to vertices is enough, size not needed 00585 contour, x1, y1, s1, 00586 true // 100222: New 'snap_ends' flag. 00587 DBG_FLAG); 00588 if (!first_inters) 00589 { 00590 if (!inner_trimming_curve) 00591 { 00592 // 100223: Disabling this warning, for with the new branch depending on "innerness" of contour, 00593 // the situation will probably (hopefully) be correctly treated in all cases now. 00594 if ( (dbg) && (0) ) 00595 puts(" Warning:\n It is likely that the 'inside'-test gave a false positive for a corner of a\n" 00596 " triangle that is really fully outside a trimming curve. We continue with this\n" 00597 " assumption, and discard the triangle in question, without attempting to clip it.\n"); 00598 } 00599 else 00600 // 100223: The new and alternate branch, which should be more likely to be correct for inner curves. 00601 add_triangle(vert_p, newmesh, c1_indx, c2_indx, c3_indx); 00602 if (dbg) puts(" ----- split_triangle done, no new triangles produced ---------------------------------\n\n"); 00603 return false; 00604 } 00605 00606 // 00607 // The first intersection is now between corner c1 and c2. The second may be between c2 and c3, or between 00608 // c3 and c1. Per assumption, c1 is *inside* the contour. (And c2 outside, c3 unknown.) 00609 // 00610 // 100223: But remember, now, that there may in fact be a lot more than the anticipated/assumed two 00611 // intersections... trying to cope better with such situations from now on. 00612 // 00613 00614 // Pushing inters1. 00615 push_an_intersection(srf, c1, c2, c1_p, c2_p, n1, n2, s1, vert, vert_p, norm, bd); 00616 if (dbg) 00617 printf(" First intersection found: s1=%f\n\n hold on; plot(%f, %f, 'm*', 'markersize', " 00618 "7, 'linewidth', 2); hold off\n\n", s1, vert_p.back()[0], vert_p.back()[1]); 00619 00620 // 100218: Note that we only have to *test* for a 2-3 intersection, since if there is none, we should have a 3-1. 00621 if (dbg) printf(" Looking for intersection on edge 2-3...\n"); 00622 const bool second_inters = 00623 segment_contour_intersection_for_s2m(c2_p[0], c2_p[1], c3_p[0], c3_p[1], &trim_curve_p[0][0], contour, x2, y2, s2, 00624 true // 100222: New 'snap_ends'-flag 00625 DBG_FLAG); 00626 00627 const double i1_i2_dist_squared = second_inters ? (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) : 0.0; 00628 if (dbg) printf(" second_inters (2-3) =%d, s2=%f, dist=%g\n", second_inters, s2, sqrt(i1_i2_dist_squared)); 00629 const bool inters1_and_inters2_distinct = i1_i2_dist_squared > abs_eps*abs_eps; 00630 00631 00632 //-------------------------------------------------------------------------------------------------------------- 00633 // 00634 // 100223: When three intersections between triangle edges and the contour are detected, the triangle is 00635 // split in four, and flagged for further trimming in a new iteration. This should take care of a 00636 // lot of previously unsolved cases with minimal new code! 00637 // 00638 //-------------------------------------------------------------------------------------------------------------- 00639 00640 bool third_inters; 00641 const bool look_for_a_third_inters = first_inters && second_inters && inters1_and_inters2_distinct; 00642 if (dbg) printf(" look_for_a_third_inters = %d\n", look_for_a_third_inters); 00643 double x3, y3, s3; 00644 if ( look_for_a_third_inters ) 00645 { 00646 if (dbg) printf(" PRE_SPLITTING_3_INTERS: Looking for intersection on edge 3-1...\n"); 00647 third_inters = segment_contour_intersection_for_s2m(c3_p[0], c3_p[1], c1_p[0], c1_p[1], 00648 &trim_curve_p[0][0], contour, x3, y3, s3, 00649 true // 100222: New 'snap_ends'-flag 00650 DBG_FLAG); 00651 if (third_inters) 00652 { 00653 const double i1_i3_dist_squared = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1); 00654 const bool inters1_and_inters3_distinct = i1_i3_dist_squared > abs_eps*abs_eps; 00655 const double i2_i3_dist_squared = (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2); 00656 const bool inters2_and_inters3_distinct = i2_i3_dist_squared > abs_eps*abs_eps; 00657 00658 if (inters1_and_inters3_distinct && inters2_and_inters3_distinct) 00659 { 00660 // puts("JEPP, FANT ET KJIPT LITE TRIANGEL..."); 00661 00662 // Pushing inters2. 00663 push_an_intersection(srf, c2, c3, c2_p, c3_p, n2, n3, s2, vert, vert_p, norm, bd); 00664 // Pushing inters3. 00665 push_an_intersection(srf, c3, c1, c3_p, c1_p, n3, n1, s3, vert, vert_p, norm, bd); 00666 00667 // Triangles... 00668 add_triangle(vert_p, newmesh, c1_indx, vert.size()-3, vert.size()-1 DBG_FLAG); 00669 add_triangle(vert_p, newmesh, c2_indx, vert.size()-2, vert.size()-3 DBG_FLAG); 00670 add_triangle(vert_p, newmesh, c3_indx, vert.size()-1, vert.size()-2 DBG_FLAG); 00671 add_triangle(vert_p, newmesh, vert.size()-3, vert.size()-2, vert.size()-1 DBG_FLAG); 00672 00673 return true; 00674 } 00675 } 00676 } 00677 00678 if ( second_inters && inters1_and_inters2_distinct && (!forced_skipping_of_second_edge) ) 00679 { 00680 //---------------------------------------------------------------------------------------------------- 00681 // 00682 // The second intersection is between corners c2 and c3, and it is distinct from the first 00683 // intersection. We can try to make the triangles as round as possible... 00684 // 00685 // 100223: What happens when 'forced_skipping_of_second_edge=true', is that this branch is not taken, 00686 // although a perfectly fine intersection distinct from inters1 was found on 2-3. That will 00687 // cause an intersection on 3-1 to be searched for, and used instead. (And who knows what will 00688 // happen if that 3-1 intersection does not exist...) If we were to enter this branch, it 00689 // seems a lot of triangles will be wrongly discarded, for some reasone. Cannot actually see 00690 // how that will happen in this branch, right now, though... Could be that one of the four 00691 // 'add_triangle' calls below try to add a degenerate triangle. Hmm... seems not so unlikely. 00692 // 00693 //---------------------------------------------------------------------------------------------------- 00694 00695 // Pushing inters2. 00696 push_an_intersection(srf, c2, c3, c2_p, c3_p, n2, n3, s2, vert, vert_p, norm, bd); 00697 if (dbg) 00698 printf(" branch 1:\n\n hold on; plot(%f, %f, 'g*', 'markersize', 7, 'linewidth', 2); hold off\n\n", 00699 vert_p.back()[0], vert_p.back()[1]); 00700 00701 if ( (x2-c1_p[0])*(x2-c1_p[0])+(y2-c1_p[1])*(y2-c1_p[1]) < (x1-c3_p[0])*(x1-c3_p[0])+(y1-c3_p[1])*(y1-c3_p[1]) ) 00702 // Splitting between c1 and inters2. 00703 { 00704 if (dbg) printf(" branch 2: Adding triangles (c1, i1, i2) and (c1, i2, c3).\n"); 00705 if (dbg) printf(" adding the first:\n"); 00706 add_triangle(vert_p, newmesh, c1_indx, vert.size()-2, vert.size()-1 DBG_FLAG); 00707 if (dbg) printf(" adding the second:\n"); 00708 add_triangle(vert_p, newmesh, c1_indx, vert.size()-1, c3_indx DBG_FLAG); 00709 } 00710 else 00711 // Splitting between c3 and inters1. 00712 { 00713 if (dbg) printf(" branch 3: Adding triangles (c1, i1, c3) and (i1, i2, c3).\n"); 00714 add_triangle(vert_p, newmesh, c1_indx, vert.size()-2, c3_indx DBG_FLAG); 00715 add_triangle(vert_p, newmesh, vert.size()-2, vert.size()-1, c3_indx DBG_FLAG); 00716 } 00717 } 00718 else 00719 { 00720 //---------------------------------------------------------------------------------------------------- 00721 // 00722 // Ok, no second intersection 2-3, then it should be on 3-1. If not, that's an error... (There should 00723 // always be exactly two intersections!) (100223: Yeah... well... not entirely true, really...) 00724 // 00725 //---------------------------------------------------------------------------------------------------- 00726 00727 if (dbg) printf(" branch 4:\n Looking for intersection on edge 3-1...\n"); 00728 bool tmp; 00729 if (look_for_a_third_inters) // If so, we have already called 'segment_contour_intersection_for_s2m'... 00730 tmp = third_inters, x2 = x3, y2 = y3, s2 = s3; 00731 else 00732 tmp = segment_contour_intersection_for_s2m(c3_p[0], c3_p[1], c1_p[0], c1_p[1], 00733 &trim_curve_p[0][0], contour, x2, y2, s2, 00734 true); // 100222: New 'snap_ends'-flag 00735 00736 // 100224: Hmm... tmp only used for debugging output? (Hence the name 'tmp'?!) 00737 00738 if (dbg) 00739 printf(" 2nd try: 3-1-inters=%d, dist=%g\n forced_skipping_of_second_edge=%d\n\n" 00740 " hold on; plot(%f, %f, 'g*', 'markersize', 8, 'linewidth', 2); hold off\n\n", 00741 tmp?1:0, sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)), forced_skipping_of_second_edge, x2, y2); 00742 00743 // 00744 // 100224: Ok, before, we had 'if (!forced...) { ... }' here, because it would be detected as a 00745 // problem if we got here, i.e., one 1-2-intersection and no 2-3-intersection *without* the 00746 // caller expecting this and therefore using 'forced...'. 00747 // 00748 // If we now want to get rid of the 'forced...' flag, we cannot trap this as a potential 00749 // problem, since we have no way of knowing if the caller expected it. If we get here, we must 00750 // simply assume that it is correct, and continue with the splitting... Therefore: We wrap the 00751 // test in the new 'do_ignore_forced_skipping'-test: 00752 // 00753 00754 if (!try_without_forced_stuff) 00755 if (!forced_skipping_of_second_edge) 00756 { 00757 MESSAGE(" Warning:\nShould not happen 2 ! Where did the second intersection go?!" 00758 //"\nJust keeping the full triangle...\n"); 00759 "\n Discarding the full triangle...\n"); 00760 fflush(stdout); 00761 fflush(stderr); 00762 //add_triangle(vert_p, newmesh, c1_indx, c2_indx, c3_indx); 00763 00764 if (dbg) 00765 { 00766 if (degenerate_triangle(vert_p[c1_indx], vert_p[c2_indx], vert_p[c3_indx] DBG_FLAG)) 00767 puts("Hmm... triangle is degenerate, so maybe appropriate to discard it then..."); 00768 int m=5*(c1_indx%4+1); 00769 printf("\n\nwwwwwwwwwwww\n\n\n" 00770 "hold on; plot(%f, %f, 'rs', 'markersize', %d, 'linewidth', 2); hold off\n" 00771 "hold on; plot(%f, %f, 'gs', 'markersize', %d, 'linewidth', 2); hold off\n" 00772 "hold on; plot(%f, %f, 'bs', 'markersize', %d, 'linewidth', 2); hold off\n\n", 00773 c1_p[0], c1_p[1], m, c2_p[0], c2_p[1], m, c3_p[0], c3_p[1], m); 00774 puts(" ----- split_triangle done ------------------------------------------------------------\n\n"); 00775 } 00776 return false; 00777 } 00778 00779 //============================================================================================================== 00780 // 090217: Now, if the first intersection is through the first corner, and there is another "non-corner" 00781 // intersection in the "middle" edge (and 'forced_skipping_of_second_edge' is true) we do some 00782 // special treatment: We split the triangle in two new ones covering the old. Note that these 00783 // must then be recursively re-processed by the caller! To make this simpler, we introduce a 00784 // boolean return value, which will be true when further processing may be required. 00785 00786 if ( ((fabs(s1)<1e-12) || (fabs(s2-1.0)<1e-12)) && forced_skipping_of_second_edge ) 00787 { 00788 double x3, y3, s3; 00789 if (dbg) printf(" Looking for intersection on edge 2-3 again...\n"); 00790 const bool tmp2= 00791 segment_contour_intersection_for_s2m(c2_p[0], c2_p[1], c3_p[0], c3_p[1], 00792 &trim_curve_p[0][0], contour, x3, y3, s3, 00793 true); // 100222: New 'snap_ends' flag. 00794 if ( tmp2 && (fabs(s3)>1e-12) && (fabs(s3-1.0)<1-1e-12) ) 00795 { 00796 // Pushing inters3. 00797 push_an_intersection(srf, c2, c3, c2_p, c3_p, n2, n3, s3, vert, vert_p, norm, bd); 00798 // Hmm... now we have added a totally unneeded point, for intersection 1... 00799 add_triangle(vert_p, newmesh, c1_indx, vert.size()-1, c3_indx); 00800 add_triangle(vert_p, newmesh, c1_indx, c2_indx, vert.size()-1); 00801 00802 if (dbg) puts(" ----- split_triangle done -----------------------------------------------------\n\n"); 00803 return true; 00804 } 00805 } 00806 00807 // Pushing inters2. 00808 push_an_intersection(srf, c3, c1, c3_p, c1_p, n3, n1, s2, vert, vert_p, norm, bd); 00809 if (dbg) 00810 printf(" s2=%g\n hold on; plot(%f, %f, 'b*', 'markersize', 7, 'linewidth', 2); hold off\n", 00811 s2, vert_p[vert_p.size()-1][0], vert_p[vert_p.size()-1][1]); 00812 00813 // There is only one triangle to push on the list in this case. 00814 if (dbg) printf(" adding the only triangle...\n"); 00815 add_triangle(vert_p, newmesh, c1_indx, vert.size()-2, vert.size()-1); 00816 00817 } // end of the block for the second intersection being of 3-1 kind instead of 2-3... 00818 00819 if (dbg) puts(" ----- split_triangle done ----------- xxx ------------------------------------------------\n\n"); 00820 return false; 00821 } 00822 00823 00824 00825 00826 00827 00828 #ifdef SPLIT_LARGE_TRIANGLES 00829 00830 //============================================================================================================== 00831 // 00832 // 100224: Making a version for splitting a triangle which would otherwise be discarded, due to all corners 00833 // being classified as outside the contour. Such a triangle might still contain intersections. 00834 // 00835 // Hmm... there is a danger of inifinite recursion here, how do we handle this? 00836 // 00837 // 100225: Testing for the need of splitting in this function too. (Note that in the 00838 // 'trim_a_triangle'/'split_triangle'-case this is divided between the two functions.) 00839 // 00840 // New name: 'split_triangle_with_all_corners_outside', better reflecting the functionality. 00841 // 00842 // Assuming all corners are outside. Not assuming anything about their "on-ness". (100225: Returning 00843 // of none are "on".) 00844 // 00845 // Adding functionality for new cases as they appear. First case: 00846 // 00847 // A corner is "on", and there is an *inner* intersection on the opposite edge. Then we split the 00848 // triangle in two here. Hopefully, this approach is so conservative that nothing will be broken 00849 // of those things already working... 00850 // 00851 // 00852 // Return value: False if no further processing of produced triangles are needed. Else true. 00853 // 00854 //============================================================================================================== 00855 00856 bool split_triangle_with_all_corners_outside(const int c1_indx, const int c2_indx, const int c3_indx, 00857 const bool c1_inside, const bool c2_inside, const bool c3_inside, 00858 const bool c1_on, const bool c2_on, const bool c3_on, 00859 boost::shared_ptr<SplineSurface> srf, 00860 vector< Vector3D > &vert, vector<Vector2D> &vert_p, 00861 vector< int > &bd, vector< Vector3D > &norm, 00862 vector<int> &newmesh, 00863 const vector< Vector3D > &trim_curve_p, 00864 const vector<int> &contour, 00865 const bool dbg = false, 00866 const bool inner_trimming_curve = false) 00867 { 00868 const double abs_eps = 1e-12; // 100218: Used for distance between two intersections in the param domain. 00869 const double snap_eps = 1e-5; 00870 00871 // 090129: Not using references in fear of what 'push_back' will do... 00872 const Vector2D c1_p=vert_p[c1_indx], c2_p=vert_p[c2_indx], c3_p=vert_p[c3_indx]; 00873 const Vector3D c1=vert[c1_indx], c2=vert[c2_indx], c3=vert[c3_indx]; 00874 const Vector3D n1=norm[c1_indx], n2=norm[c2_indx], n3=norm[c3_indx]; 00875 00876 if (dbg) 00877 { 00878 puts("\n\n ----- split_triangle_with_all_corners_outside ------------------------------------------------"); 00879 printf("\n The corners:\n\n"); 00880 printf("hold on; plot(%f, %f, 'rd', 'markersize', %d, 'linewidth', 2); hold off\n", c1_p[0], c1_p[1], 14); 00881 printf("hold on; plot(%f, %f, 'gd', 'markersize', %d, 'linewidth', 2); hold off\n", c2_p[0], c2_p[1], 14); 00882 printf("hold on; plot(%f, %f, 'bd', 'markersize', %d, 'linewidth', 2); hold off\n\n", c3_p[0], c3_p[1], 14); 00883 } 00884 00885 if (!(c1_on || c2_on || c3_on)) 00886 // 100225: At this time, we don't do anything with such a triangle. 00887 return false; 00888 00889 // 100225: Ok, at least one corner is "on", so let's make sure c1 is "on" 00890 if (!c1_on) 00891 { 00892 if (c2_on) 00893 return split_triangle_with_all_corners_outside(c2_indx, c3_indx, c1_indx, 00894 c2_inside, c3_inside, c1_inside, 00895 c2_on, c3_on, c1_on, 00896 srf, vert, vert_p, bd, norm, newmesh, trim_curve_p, 00897 contour, dbg, inner_trimming_curve); 00898 else 00899 return split_triangle_with_all_corners_outside(c3_indx, c1_indx, c2_indx, 00900 c3_inside, c1_inside, c2_inside, 00901 c3_on, c1_on, c2_on, 00902 srf, vert, vert_p, bd, norm, newmesh, trim_curve_p, 00903 contour, dbg, inner_trimming_curve); 00904 } 00905 00906 // 100225: Now c1 is "on". Checking for an intersection (inner) on edge 2-3. 00907 00908 double x1, y1, s1; 00909 if (dbg) printf(" Looking for intersection on edge 2-3...\n"); 00910 const bool inters = 00911 segment_contour_intersection_for_s2m(c2_p[0], c2_p[1], c3_p[0], c3_p[1], 00912 &trim_curve_p[0][0], // Pointer to vertices is enough, size not needed 00913 contour, x1, y1, s1, 00914 true // 100222: New 'snap_ends' flag. 00915 DBG_FLAG); 00916 if (!inters) 00917 return false; 00918 00919 const double dist2_squared = (x1-c2_p[0])*(x1-c2_p[0]) + (y1-c2_p[1])*(y1-c2_p[1]); 00920 const double dist3_squared = (x1-c3_p[0])*(x1-c3_p[0]) + (y1-c3_p[1])*(y1-c3_p[1]); 00921 const bool interior = (dist2_squared>snap_eps*snap_eps) && (dist3_squared>snap_eps*snap_eps); 00922 if (dbg) printf(" dist2=%f, dist3=%f, interior=%d\n", sqrt(dist2_squared), sqrt(dist3_squared), interior); 00923 if (!interior) 00924 return false; 00925 00926 // 100225: Ok, all conditions are met, we split. 00927 00928 push_an_intersection(srf, c2, c3, c2_p, c3_p, n2, n3, s1, vert, vert_p, norm, bd); 00929 add_triangle(vert_p, newmesh, c1_indx, c2_indx, vert.size()-1 DBG_FLAG); 00930 add_triangle(vert_p, newmesh, c1_indx, vert.size()-1, c2_indx DBG_FLAG); 00931 if (dbg) puts(" ----- split_triangle_with_all_corners_outside, produced two new triangles ----------------\n\n"); 00932 return true; 00933 } 00934 00935 #endif 00936 00937 00938 00939 00940 00941 00942 //============================================================================================================== 00943 // 00944 // 090202: Hmm... Should maybe introduce some tolerances here too... 00945 // 00946 // 090202: Note that this function is supposed to handle 3D points, even though we in this case have 2D points 00947 // in the parameter domain. (In the 3D case, the point tested will be the projection onto the plane 00948 // containing the triangle, if I remember correctly...) I think just replacing Vector3D with Vector2D 00949 // actually will do the trick... 00950 // 00951 //============================================================================================================== 00952 00953 inline bool pt_inside_tri(const Vector3D &p, const Vector3D &c1, const Vector3D &c2, const Vector3D &c3) 00954 { 00955 const Vector3D e=c2-c1, f=c3-c1, pc=p-c1; 00956 const double ee=e*e, ff=f*f, ef=e*f; 00957 #if 0 00958 // This is faster on x86 than allocating new variables for pc*e and pc*f, it seems... 00959 const double u = ff*(pc*e) - ef*(pc*f); 00960 const double v = ee*(pc*f) - ef*(pc*e); 00961 #else 00962 // Not on core2 is seems... 00963 const double pce=pc*e, pcf=pc*f; 00964 const double u = ff*pce - ef*pcf; 00965 const double v = ee*pcf - ef*pce; 00966 #endif 00967 00968 return ((u>=0.0) && (v>=0.0) && (u+v<=ee*ff-ef*ef)); 00969 } 00970 00971 inline bool pt_inside_tri(const Vector2D &p, const Vector2D &c1, const Vector2D &c2, const Vector2D &c3) 00972 { 00973 const Vector2D e=c2-c1, f=c3-c1, pc=p-c1; 00974 const double ee=e*e, ff=f*f, ef=e*f; 00975 #if 0 00976 // This is faster on x86 than allocating new variables for pc*e and pc*f, it seems... 00977 const double u = ff*(pc*e) - ef*(pc*f); 00978 const double v = ee*(pc*f) - ef*(pc*e); 00979 #else 00980 // Not on core2 is seems... 00981 const double pce=pc*e, pcf=pc*f; 00982 const double u = ff*pce - ef*pcf; 00983 const double v = ee*pcf - ef*pce; 00984 #endif 00985 00986 return ((u>=0.0) && (v>=0.0) && (u+v<=ee*ff-ef*ef)); 00987 } 00988 00989 //============================================================================================================== 00990 // 00991 // 090203: To avoid splitting triangles into almost degenerate ones. 00992 // 00993 // Assuming 'pt_inside_tri' is true, so that we can measure distance to infinite extensions of edges 00994 // for simplicity. 00995 // 00996 //============================================================================================================== 00997 00998 inline double pt_dist_from_line(const Vector2D &p, const Vector2D &A, const Vector2D &B) 00999 { 01000 const double t = ((p-A)*(B-A))/((B-A)*(B-A)); 01001 const Vector2D proj = A+t*(B-A); 01002 return sqrt((proj-p)*(proj-p)); 01003 } 01004 01005 inline double pt_dist_from_tri_edge(const Vector2D &p, const Vector2D &c1, const Vector2D &c2, const Vector2D &c3) 01006 { 01007 return std::min(pt_dist_from_line(p, c1, c2), 01008 std::min(pt_dist_from_line(p, c2, c3), pt_dist_from_line(p, c3, c1))); 01009 } 01010 01011 01012 01013 01014 01015 01016 //============================================================================================================== 01017 // 01018 // 090217: Encapsulating this for clarity. 01019 // 01020 //============================================================================================================== 01021 01022 bool trim_a_triangle(boost::shared_ptr<SplineSurface> srf, 01023 vector<Vector3D> &vert, 01024 vector<Vector2D> &vert_p, 01025 vector< int > &bd, vector<Vector3D> &norm, 01026 const vector<Vector3D> &trim_curve_p, 01027 const vector<int> &contour, 01028 const int c1_indx, const int c2_indx, const int c3_indx, 01029 vector<int> &mesh, 01030 const bool inner_trimming_curve = false 01031 #ifdef DBG 01032 , const bool dbg = false 01033 #endif 01034 ) 01035 { 01036 #ifndef DBG 01037 const bool dbg=false; 01038 #endif 01039 if (dbg) 01040 printf("\n\n= trim_a_triangle starting =====================================================================\n"); 01041 01042 vector<int> new_mesh; 01043 01044 const Vector3D c1=vert[c1_indx], c2=vert[c2_indx], c3=vert[c3_indx]; 01045 const Vector3D n1=norm[c1_indx], n2=norm[c2_indx], n3=norm[c3_indx]; 01046 const Vector2D c1_p=vert_p[c1_indx], c2_p=vert_p[c2_indx], c3_p=vert_p[c3_indx]; 01047 int s0, s, s3; 01048 01049 if (dbg) 01050 { 01051 printf("\n The corners:\n\n"); 01052 printf("hold on; plot(%f, %f, 'rx', 'markersize', %d, 'linewidth', 2); hold off\n", c1_p[0], c1_p[1], 14); 01053 printf("hold on; plot(%f, %f, 'gx', 'markersize', %d, 'linewidth', 2); hold off\n", c2_p[0], c2_p[1], 14); 01054 printf("hold on; plot(%f, %f, 'bx', 'markersize', %d, 'linewidth', 2); hold off\n\n", c3_p[0], c3_p[1], 14); 01055 } 01056 01057 if (dbg) printf(" = Calling is_inside (point_inside_contour) for corner c1... ============================\n"); 01058 const int c1_inside = is_inside(trim_curve_p, contour, c1_p[0], c1_p[1] DBG_FLAG); 01059 if (dbg) printf(" = Calling is_inside (point_inside_contour) for corner c2... ============================\n"); 01060 const int c2_inside = is_inside(trim_curve_p, contour, c2_p[0], c2_p[1] DBG_FLAG); 01061 if (dbg) printf(" = Calling is_inside (point_inside_contour) for corner c3... ============================\n"); 01062 const int c3_inside = is_inside(trim_curve_p, contour, c3_p[0], c3_p[1] DBG_FLAG); 01063 01064 if (dbg) printf(" = Calling is_on_contour (point_on_contour) for corner c1... ============================\n"); 01065 const int c1_on = is_on_contour(trim_curve_p, contour, c1_p[0], c1_p[1] DBG_FLAG); 01066 if (dbg) printf(" = Calling is_on_contour (point_on_contour) for corner c2... ============================\n"); 01067 const int c2_on = is_on_contour(trim_curve_p, contour, c2_p[0], c2_p[1] DBG_FLAG); 01068 if (dbg) printf(" = Calling is_on_contour (point_on_contour) for corner c3... ============================\n"); 01069 const int c3_on = is_on_contour(trim_curve_p, contour, c3_p[0], c3_p[1] DBG_FLAG); 01070 01071 if (dbg) 01072 printf("----------------------------------------------------------------------------------------------------\n" 01073 "IN/ON-tests before considering whether curve is outer or inner:\n" 01074 "c1_inside=%d, c2_inside=%d, c3_inside=%d, c1_on=%d, c2_on=%d, c3_on=%d\n" 01075 "----------------------------------------------------------------------------------------------------\n", 01076 c1_inside, c2_inside, c3_inside, c1_on, c2_on, c3_on); 01077 01078 // 100211: s2==0 <=> no corner on curve itself. 01079 const int s2 = (c1_on << 0) + (c2_on << 1) + (c3_on << 2); 01080 01081 const Vector2D centroid = (1.0/3.0)*(c1_p+c2_p+c3_p); 01082 const int c_inside = is_inside(trim_curve_p, contour, centroid[0], centroid[1]); 01083 01084 if (inner_trimming_curve) 01085 { 01086 // 100211: s0==0 <=> all corners are not inside. (Is "on edge" included in "inside"?) 01087 s0 = ((!c1_inside) << 0) + ((!c2_inside) << 1) + ((!c3_inside) << 2); 01088 01089 // 100211: 01090 s = (((!c1_inside) | c1_on) << 0) + (((!c2_inside) | c2_on) << 1) + (((!c3_inside) | c3_on) << 2); 01091 01092 // 100211: s3==true <=> "centroid is not inside". 01093 s3 = !c_inside; 01094 01095 if (dbg) 01096 // if ((s0!=7) || (s2!=0) || (s!=7)) 01097 printf(" s0=%d, s2=%d, s=%d, s3=%d\n", s0, s2, s, s3); 01098 01099 if ( (s0==0) && ((s2==1) || (s2==2) || (s2==4)) ) 01100 { 01101 if (dbg) 01102 printf("No points on the inside, exactly one *on* the contour. Setting s=0 and skipping triangle. 1\n"); 01103 s=0; 01104 } 01105 if ( ((s0==1) && (s2==1)) || ((s0==2) && (s2==2)) || ((s0==4) && (s2==4)) ) 01106 { 01107 if (dbg) 01108 printf("Two points not on the inside, the remaining both inside and *on*. Skipping triangle. 2\n"); 01109 s=0; 01110 } 01111 01112 if (dbg) 01113 if ( (s==7) && (s3==0) ) 01114 { 01115 printf("All corners on the inside, mean not, skipping triangle.\n"); 01116 // 090220: This is not so straightforward... Not all these triangles are either fully inside or 01117 // outside. They should really be split, but then we should split neighbours also, if we 01118 // ever are to get a valid triangulation. (Which we are currently not getting, since we 01119 // produce lots of duplicate nodes and so on...) 01120 //s=0; 01121 } 01122 01123 if (dbg) 01124 if ( ((s2==3) && (s0==4)) || 01125 ((s2==5) && (s0==2)) || 01126 ((s2==6) && (s0==1)) ) 01127 { 01128 printf("Two points on the curve, one inside. Skipping triangle.\n"); 01129 // s=0; 01130 } 01131 01132 if (dbg) 01133 if (s2==7) 01134 { 01135 printf("All points on contour, guessing that the triangle should be excluded. \n"); 01136 printf(" s0=%d, s2=%d, s=%d, s3=%d\n", s0, s2, s, s3); 01137 // s=0; 01138 } 01139 01140 } 01141 else 01142 { 01143 s0 = (c1_inside << 0) + (c2_inside << 1) + (c3_inside << 2); 01144 s = ((c1_inside | c1_on) << 0) + ((c2_inside | c2_on) << 1) + ((c3_inside | c3_on) << 2); 01145 s3 = c_inside; 01146 01147 if (dbg) 01148 printf("---> s0=%d s2=%d s=%d s3=%d\n", s0, s2, s, s3), fflush(stdout); 01149 01150 if ( (s0==0) && ((s2==1) || (s2==2) || (s2==4)) ) 01151 { 01152 if (dbg) 01153 printf("No points on the inside, exactly one *on* the contour. Setting s=0 and skipping triangle. 4\n"); 01154 // s=0; 01155 } 01156 01157 if ( ((s0==1) && (s2==1)) || ((s0==2) && (s2==2)) || ((s0==4) && (s2==4)) ) 01158 { 01159 if (dbg) 01160 printf("Two points not inside, the remaining both inside and *on*. Skipping triangle. s=%d 5\n", s); 01161 //s=0; 01162 //s=7; 01163 } 01164 01165 if (dbg) 01166 if ((s==7) && (!s3)) 01167 { 01168 printf("====================> HUH?! s0=%d s2=%d s=%d s3=%d\n", s0, s2, s, s3); 01169 // s=0; 01170 } 01171 01172 if (dbg) 01173 printf("---> s0=%d s2=%d s=%d s3=%d\n", s0, s2, s, s3), fflush(stdout); 01174 01175 } 01176 01177 if (degenerate_triangle(c1_p, c2_p, c3_p)) 01178 { 01179 printf("Degenerate triangle! setting s=-1 and skipping triangle. 6\n"); 01180 s=-1; 01181 } 01182 01183 01184 // 01185 // 100223: Hmm... What is really the pattern with respect to the usage of 01186 // 'forced_skipping_of_second_edge'-values used below? 01187 // 01188 // Note that if s==0, the triangle will silently be forgotten... 01189 // 01190 // 100225: No, adding a stage for this situation too. 01191 // 01192 // Hmm... Maybe it's time to stop letting "on" imply "inside"...?! 01193 // 01194 01195 #ifdef SPLIT_LARGE_TRIANGLES 01196 if ( (!(c1_inside || c2_inside || c3_inside)) && (!inner_trimming_curve) && (c1_on || c2_on || c3_on) ) 01197 s=0; 01198 #endif 01199 01200 01201 01202 01203 bool redo=false; 01204 if (dbg) 01205 printf(" Going into the switch, s=%d\n", s); 01206 switch (s) 01207 { 01208 #ifdef SPLIT_LARGE_TRIANGLES 01209 case 0: 01210 if (dbg) printf("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n"); 01211 if (!inner_trimming_curve) 01212 redo = split_triangle_with_all_corners_outside(c1_indx, c2_indx, c3_indx, 01213 c1_inside, c2_inside, c3_inside, 01214 c1_on, c2_on, c3_on, 01215 srf, vert, vert_p, bd, norm, new_mesh, trim_curve_p, contour, 01216 dbg, inner_trimming_curve); 01217 if (dbg) printf("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n"); 01218 break; 01219 #endif 01220 case 5: // Corners 1 and 3 are inside, intersections must be on edges 1-2 and 2-3. Compatible with 1-case. 01221 // 090217: Not now, don't want the 'forced skip'. 01222 redo=split_triangle(c1_indx, c2_indx, c3_indx, srf, vert, vert_p, bd, norm, new_mesh, 01223 trim_curve_p, contour, dbg, false, inner_trimming_curve); 01224 break; 01225 case 1: // Only corner 1 is inside, so intersections must be on edges 1-2 and 3-1. 01226 redo=split_triangle(c1_indx, c2_indx, c3_indx, srf, vert, vert_p, bd, norm, new_mesh, 01227 trim_curve_p, contour, dbg, true, inner_trimming_curve); 01228 break; 01229 case 2: // Corner 2 is inside, intersections must be on edges 2-3 and 3-1. Can go straight into the 3-case. 01230 01231 // 090217: Huh?! Should be 2-3 and 1-2, shouldn't it?! Think maybe the comment was wrong, but 01232 // conclusion and code correct. Sadly, this will only work when the assumption about no 01233 // "spurious" intersections is met. I.e., that there are no "superfluous" 01234 // intersections. In many cases this does not hold. Trying to fix this by not combining 01235 // cases here, and at the same time, not letting split_triangle choose to continue on the 01236 // third edge if it does not find any intersection on the second edge. The problem seems 01237 // to be that such a spurious intersection can be on edge 2, causing split_triangle not to 01238 // look for the (correct and sought) intersection on edge 3... 01239 01240 redo=split_triangle(c2_indx, c3_indx, c1_indx, srf, vert, vert_p, bd, norm, new_mesh, 01241 trim_curve_p, contour, dbg, true, inner_trimming_curve); 01242 break; 01243 case 3: // Corners 1 and 2 are inside, so intersections must be on edges 2-3 and 3-1. 01244 redo=split_triangle(c2_indx, c3_indx, c1_indx, srf, vert, vert_p, bd, norm, new_mesh, 01245 trim_curve_p, contour, dbg, false, inner_trimming_curve); 01246 break; 01247 case 6: // Corners 2 and 3 are inside, intersections must be on edges 3-1 and 1-2. Compatible with 4-case. 01248 // 090217: Not now, don't want the 'forced skip'. 01249 redo=split_triangle(c3_indx, c1_indx, c2_indx, srf, vert, vert_p, bd, norm, new_mesh, 01250 trim_curve_p, contour, dbg, false, inner_trimming_curve); 01251 break; 01252 case 4: // Corner 3 is inside, so intersections must be on edges 3-1 and 2-3. 01253 redo=split_triangle(c3_indx, c1_indx, c2_indx, srf, vert, vert_p, bd, norm, new_mesh, 01254 trim_curve_p, contour, dbg, true, inner_trimming_curve); 01255 break; 01256 case 7: // All corners are inside, not splitting, just keeping it. 01257 add_triangle(vert_p, new_mesh, c1_indx, c2_indx, c3_indx); 01258 redo=false; 01259 break; 01260 } 01261 mesh=new_mesh; 01262 01263 if (dbg) 01264 printf(" redo=%d\n= trim_a_triangle done =======================================================" 01265 "==================\n\n\n", redo); 01266 01267 return redo; 01268 } 01269 01270 01271 01272 01273 01274 01275 //============================================================================================================== 01276 // 01277 // More or less just a loop over 'triam_a_triangle' for a list of triangles. 01278 // 01279 // 100211: For debugging, a matlab script is written for each triangle in 'mesh', rendering the triangle 01280 // plus those produced. Note that (currently) this is only done during the first iteration. 01281 // 01282 // Magenta: This is being processed. One for each tats-file. 01283 // Red: These were produced from the red one. 01284 // 01285 //============================================================================================================== 01286 01287 void trim_a_triangle_soup(boost::shared_ptr<SplineSurface> srf, 01288 vector<Vector3D> &vert, 01289 vector<Vector2D> &vert_p, 01290 vector< int > &bd, vector<Vector3D> &norm, 01291 const vector<Vector3D> &trim_curve_p, 01292 const vector<int> &contour, 01293 vector<int> &mesh, 01294 #ifdef DBG 01295 vector<char> &col, 01296 #endif 01297 const bool inner_trimming_curve = false) 01298 { 01299 vector<int> new_mesh, mesh_to_reiterate, mesh_ok; 01300 const int maxiter=10; 01301 int iter=0; 01302 01303 #ifdef DBG 01304 system("rm tats*.m"); 01305 #endif 01306 01307 do 01308 { 01309 01310 for (int ilim=mesh.size()/3, i=0; i<ilim; i++) 01311 //for (int ilim=mesh.size()/3, i=0; i<1; i++) 01312 { 01313 const int c1_indx=mesh[3*i], c2_indx=mesh[3*i+1], c3_indx=mesh[3*i+2]; 01314 bool redo=trim_a_triangle(srf, vert, vert_p, bd, norm, trim_curve_p, contour, 01315 c1_indx, c2_indx, c3_indx, new_mesh, inner_trimming_curve 01316 #ifdef DBG 01317 // 100223: Remember that when having more trimming curves, this may be 01318 // triggered once for each of them... 01319 ,false // , i==423 // , false // , i==1069 // , false // , i==94 01320 #endif 01321 ); 01322 if (redo) 01323 mesh_to_reiterate.insert(mesh_to_reiterate.end(), new_mesh.begin(), new_mesh.end()); 01324 else 01325 mesh_ok.insert(mesh_ok.end(), new_mesh.begin(), new_mesh.end()); 01326 #ifdef DBG 01327 // 100210: Showing both the triangle being processed (magenta) and the ones being produced (red.) 01328 if ( (iter==0) || (iter==1) ) 01329 #if 0 01330 if ( 01331 01332 ( ((fabs(vert_p[mesh[3*i ]][0]-(-22.56))<1e-5) && (fabs(vert_p[mesh[3*i ]][1]-(30.888))<1e-5)) || 01333 ((fabs(vert_p[mesh[3*i+1]][0]-(-22.56))<1e-5) && (fabs(vert_p[mesh[3*i+1]][1]-(30.888))<1e-5)) || 01334 ((fabs(vert_p[mesh[3*i+2]][0]-(-22.56))<1e-5) && (fabs(vert_p[mesh[3*i+2]][1]-(30.888))<1e-5)) ) && 01335 01336 (vert_p[mesh[3*i]][0] > -26) && 01337 (vert_p[mesh[3*i]][0] < -17) && 01338 (vert_p[mesh[3*i]][1] > 28) && 01339 (vert_p[mesh[3*i]][1] < 36) ) 01340 #endif 01341 { 01342 char fname[1000]; 01343 if (iter==0) 01344 sprintf(fname, "tats%04d.m", i); 01345 else 01346 sprintf(fname, "xtats%04d.m", i); 01347 FILE *f = fopen(fname, "w"); 01348 fprintf(f, "hold on\ngrid on\n"); 01349 const int u=0, v=1; 01350 const char col1 = iter==0 ? 'm' : 'b'; 01351 const char col2 = iter==0 ? 'r' : 'g'; 01352 fprintf(f, "patch([%f; %f; %f], [%f; %f; %f], '%c');\n", 01353 vert_p[mesh[3*i]][u], vert_p[mesh[3*i+1]][u], vert_p[mesh[3*i+2]][u], 01354 vert_p[mesh[3*i]][v], vert_p[mesh[3*i+1]][v], vert_p[mesh[3*i+2]][v], col1); 01355 fprintf(f, "line([%f; %f; %f; %f], [%f; %f; %f; %f], 'color', 'k');\n", 01356 vert_p[mesh[3*i]][u], vert_p[mesh[3*i+1]][u], vert_p[mesh[3*i+2]][u], vert_p[mesh[3*i]][u], 01357 vert_p[mesh[3*i]][v], vert_p[mesh[3*i+1]][v], vert_p[mesh[3*i+2]][v], vert_p[mesh[3*i]][v]); 01358 for (int j=0; j<int(new_mesh.size())/3; j++) 01359 { 01360 fprintf(f, "patch([%f; %f; %f], [%f; %f; %f], '%c');\n", 01361 vert_p[new_mesh[3*j]][u], vert_p[new_mesh[3*j+1]][u], vert_p[new_mesh[3*j+2]][u], 01362 vert_p[new_mesh[3*j]][v], vert_p[new_mesh[3*j+1]][v], vert_p[new_mesh[3*j+2]][v], col2); 01363 fprintf(f, "line([%f; %f; %f; %f], [%f; %f; %f; %f], 'color', 'k');\n", 01364 vert_p[new_mesh[3*j ]][u], vert_p[new_mesh[3*j+1]][u], 01365 vert_p[new_mesh[3*j+2]][u], vert_p[new_mesh[3*j ]][u], 01366 vert_p[new_mesh[3*j ]][v], vert_p[new_mesh[3*j+1]][v], 01367 vert_p[new_mesh[3*j+2]][v], vert_p[new_mesh[3*j ]][v]); 01368 } 01369 fprintf(f, "xlabel('magenta/blue=old triangle, red/green=new ones');\n"); 01370 fprintf(f, "hold off\n"); 01371 fclose(f); 01372 } 01373 #endif 01374 } 01375 01376 #ifdef DBG 01377 if (iter==0) 01378 if (inner_trimming_curve) 01379 col=vector<char>(mesh_ok.size()/3, 'm'); 01380 else 01381 col=vector<char>(mesh_ok.size()/3, 'c'); 01382 #endif 01383 01384 mesh=mesh_to_reiterate; 01385 iter++; 01386 01387 } 01388 while ( (mesh.size()>0) && (iter<maxiter) ); 01389 01390 mesh=mesh_ok; 01391 #ifdef DBG 01392 while (col.size()<mesh.size()/3) 01393 if (inner_trimming_curve) 01394 col.push_back('y'); 01395 else 01396 col.push_back('b'); 01397 #endif 01398 } 01399 01400 01401 01402 01403 01404 01405 //============================================================================================================== 01406 // 01407 // 090218: Pulling this out in a routine of it's own. 01408 // 090219: Making sure duplicate points do not survive. 01409 // 01410 //============================================================================================================== 01411 01412 void construct_corner_lists(boost::shared_ptr<SplineSurface> srf, 01413 vector<boost::shared_ptr<SplineCurve> >& crv_set, 01414 const vector< Vector3D > &vert, 01415 const vector< Vector2D > &vert_p, 01416 const vector< Vector3D > &norm, 01417 const vector<int> &mesh, 01418 const vector< vector< Vector3D > > &trim_curve_all, 01419 const vector< vector< Vector3D > > &trim_curve_p_all, 01420 const int dn, const int dm, // Will generate an (dn+1)x(dm+1)-mesh... 01421 const double bd_res_ratio, 01422 const double eps, 01423 const double cosangle_limit, 01424 const double pt_edge_degen_limit, 01425 const double pt_mult_def, 01426 const double min_corner_dist, const double max_corner_dist, 01427 01428 vector<int> &skip_quad, 01429 vector< vector<Vector3D> > &quad_corner_trim_curve, 01430 vector< vector<Vector2D> > &quad_corner_trim_curve_p, 01431 int &quad_corners) 01432 { 01433 // 090130: "Corner-handling", 1st step: Identify quads needing special treatment, i.e. extra refinement. 01434 // 090203: Note that we really want to avoid making almost-degenerate triangles... 01435 // Hmm... This should really be something relative to the extent of the quad... 01436 //const double delta=1e-8; // See above. 01437 skip_quad=vector<int>(dn*dm, 0); 01438 quad_corner_trim_curve.resize(dn*dm); 01439 quad_corner_trim_curve_p.resize(dn*dm); 01440 01441 const double u0=srf->startparam_u(); 01442 const double u1=srf->endparam_u(); 01443 const double v0=srf->startparam_v(); 01444 const double v1=srf->endparam_v(); 01445 const double dv=(v1-v0)/dm, du=(u1-u0)/dn; 01446 01447 #ifdef DBG 01448 // 100222: Producing matlab-scripts for debugging. 01449 { 01450 const double u02=vert_p[0][0], v02=vert_p[0][1]; 01451 const double u12=vert_p[(dn+1)*(dm+1)-1][0], v12=vert_p[(dn+1)*(dn+1)-1][1]; 01452 if (fabs(u02-u0)>eps) THROW("Huh?!"); 01453 if (fabs(u12-u1)>eps) THROW("Huh?!"); 01454 if (fabs(v02-v0)>eps) THROW("Huh?!"); 01455 if (fabs(v12-v1)>eps) THROW("Huh?!"); 01456 printf("Parameter domain: [%g, %g] x [%g, %g]\n", u02, u12, v02, v12); 01457 } 01458 FILE *f1=fopen("c1.m", "w"); // corners in quads 01459 FILE *f2=fopen("c2.m", "w"); // corner skipped due to small angle 01460 FILE *f3=fopen("c3.m", "w"); // corner skipped, is a duplicate 01461 FILE *f4=fopen("c4.m", "w"); // all points 01462 fprintf(f1, "hold on\ngrid on\n"); 01463 fprintf(f2, "hold on\ngrid on\n"); 01464 fprintf(f3, "hold on\ngrid on\n"); 01465 fprintf(f4, "hold on\ngrid on\n"); 01466 #endif 01467 01468 for (int curve=0; curve<int(crv_set.size()); curve++) 01469 { 01470 #ifdef CLC_DBG 01471 printf("curve %d\n", curve); 01472 #endif 01473 const vector<Vector3D> &trim_curve_p=trim_curve_p_all[curve]; 01474 const vector<Vector3D> &trim_curve=trim_curve_all[curve]; 01475 Vector2D last_dir(1.0, 0.0); // Will be kept normalized. The initial value will not be used, 01476 // ref. 'a_first_point_added' 01477 bool a_first_point_added=false; 01478 Vector2D last_used_p(0.0, 0.0); 01479 const int contour_points=trim_curve_p.size(); 01480 01481 for (int i=0; i<contour_points; i++) 01482 { 01483 double u=trim_curve_p[i][0], v=trim_curve_p[i][1]; 01484 // const int /* ip1=(i+1)%contour_points, */ im1=(i+contour_points-1)%contour_points; 01485 #ifdef CLC_DBG 01486 printf(" i=%4d u=%7.3f, v=%7.3f", i, u, v); 01487 #endif 01488 if ((u>=u0) && (u<=u1) && (v>=v0) && (v<=v1)) 01489 { 01490 const int p=std::min(int(floor((v-v0)/dv)), dm-1), q=std::min(int(floor((u-u0)/du)), dn-1); 01491 const int vert_indx=p*(dn+1)+q, quad_indx=p*dn+q; 01492 const double &u_left=vert_p[vert_indx][0], &v_left=vert_p[vert_indx][1]; 01493 const double &u_right=vert_p[vert_indx+dn+2][0], &v_right=vert_p[vert_indx+dn+2][1]; 01494 // if (u<u_left) 01495 // { 01496 // printf("\nu=%g, u_left=%g\n", u, u_left); 01497 // THROW("Huh?!"); 01498 // } 01499 // if (v<v_left) THROW("Huh?!"); 01500 // if ((u>u_right) && (q<n-1)) THROW("Huh?!"); 01501 // if ((v>v_right) && (p<n-1)) THROW("Huh?!"); 01502 // Note that due to rounding errors, discrete arithmetic, or floor's behaviour, u or v might 01503 // actually be outside the quads parameter (sub)domain. This should only happen on the right and 01504 // upper rim... 01505 // 090219: Doing a small adjustment for this... 01506 u=std::min(std::max(u_left, u), u_right); 01507 v=std::min(std::max(v_left, v), v_right); 01508 01509 // Need two things: 1) angle!=pi, 2) not duplicate point. And 3) not on the edges. 01510 // 090218: New criterium: Accumulating angle changes. 01511 // 090220: Taking care of weeding out duplicates in the calling function. (?) 01512 01513 if ( ((u-u_left<pt_edge_degen_limit) || (u_right-u<pt_edge_degen_limit) || 01514 (v-v_left<pt_edge_degen_limit) || (v_right-v<pt_edge_degen_limit) ) && (0) ) 01515 { 01516 //printf(" quad-corner too close to edge to split.\n"); 01517 // 090202: Hmmm... On second thought... We should still split here, it's just that we 01518 // shouldn't split into four new triangles. 01519 // 090203: Hmm... On third thought... It's ok *not* to "corner split" here, because 01520 // everything will be ok when we do the "ordinary" curve-splitting later. The 01521 // point about the corners is to force splitting in the corner, and when the 01522 // corner happens to be on a mesh-edge, this will happen automagically! 01523 // 090203: Note that even if we end up here, the angle, distance to previous point 01524 // etc. may still not warrant a splitting. This is likely often the case when the 01525 // trimming curve follows an edge of the surface! 01526 // 090203: No! If we don't add a corner even though on an edge, splitting may not occur 01527 // later either. (But why not? Was this not solved with the "boundary-fix" some 01528 // time earlier?) Another solution is to split, and then remove (or not add) 01529 // degenerate triangles... 01530 } 01531 else 01532 { 01533 const Vector2D current_p(trim_curve_p[i][0], trim_curve_p[i][1]); 01534 const Vector2D previous_p = 01535 i>0 ? Vector2D(trim_curve_p[i-1][0], trim_curve_p[i-1][1]) : Vector2D(1e99, 1e99); 01536 const double dist_since_last_used_p=sqrt((current_p-last_used_p)*(current_p-last_used_p)); 01537 const double dist_since_prev_p=sqrt((current_p-previous_p)*(current_p-previous_p)); 01538 01539 #ifdef DBG 01540 fprintf(f4, "plot(%f, %f, 'kd', 'markersize', 5);\n", current_p[0], current_p[1]); 01541 #endif 01542 01543 if ( (dist_since_prev_p>pt_mult_def) // a new distinct point 01544 || (!a_first_point_added) ) 01545 { 01546 const double cosangle = (last_dir*(current_p-last_used_p))/dist_since_last_used_p; 01547 // (Note: We have made sure 'dist_since_last_used_p' is non-zero!) 01548 #ifdef CLC_DBG 01549 printf(" [%7.3f %7.3f %7.3f %7.3f]", last_dir[0], last_dir[1], 01550 current_p[0]-last_used_p[0], current_p[1]-last_used_p[1]); 01551 printf(" angle=%7.3f", acos(std::min(cosangle, 1.0))/M_PI*180.0); 01552 printf(" cosangle=%6.3f", cosangle); 01553 #endif 01554 // if ( ((cosangle<cosangle_limit) || (b*b>max_corner_dist*max_corner_dist)) && 01555 // (b*b>min_corner_dist*max_corner_dist) ) 01556 if ((cosangle<cosangle_limit) || (!a_first_point_added)) 01557 { 01558 //const Vector2D previous_p(trim_curve_p[im1][0], trim_curve_p[im1][1]); 01559 // if ( (current_p-previous_p)*(current_p-previous_p) < 1e-15*1e-15 ) 01560 // { 01561 // #ifdef DBG 01562 // printf(" skipping point, same as previous. "); 01563 // #endif 01564 // } 01565 // else 01566 // { 01567 // Now we flag the quad as being split, and append the corner to the quad's list. 01568 //printf("dist from previous: %g ", sqrt((current_p-previous_p)*(current_p-previous_p))); 01569 skip_quad[quad_indx]=1; 01570 quad_corner_trim_curve[quad_indx].push_back(trim_curve[i]); 01571 quad_corner_trim_curve_p[quad_indx].push_back(current_p); 01572 quad_corners++; 01573 last_dir=current_p-last_used_p; 01574 //printf("\n last_dir 1 =%g %g\n", last_dir[0], last_dir[1]); 01575 last_dir.normalize(); 01576 //printf("\n last_dir 2 =%g %g\n", last_dir[0], last_dir[1]); 01577 last_used_p=current_p; 01578 a_first_point_added=true; 01579 #ifdef DBG 01580 fprintf(f1, "plot(%f, %f, 'm*');\n", current_p[0], current_p[1]); 01581 char *tmpstr=new char[100]; 01582 sprintf(tmpstr, "%d", i); 01583 fprintf(f1, "text(%f, %f, '%s');\n", current_p[0], current_p[1]+0.00005, tmpstr); 01584 # ifdef CLC_DBG 01585 printf(" (c1) USED"); 01586 # endif 01587 #endif 01588 // } 01589 } 01590 #ifdef DBG 01591 else 01592 { 01593 fprintf(f2, "plot(%f, %f, 'b+');\n", current_p[0], current_p[1]); 01594 char *tmpstr=new char[100]; 01595 sprintf(tmpstr, "%d", i); 01596 fprintf(f2, "text(%f, %f, '%s');\n", current_p[0], current_p[1]+0.00010, tmpstr); 01597 # ifdef CLC_DBG 01598 printf(" (c2) TOO SMALL ANGLE"); 01599 # endif 01600 } 01601 #endif 01602 } 01603 else 01604 { 01605 // Use it, but not if already used!!!! 01606 if (dist_since_last_used_p>1e-14) 01607 { 01608 01609 #ifdef DBG 01610 fprintf(f3, "plot(%f, %f, 'go', 'markersize', 4);\n", current_p[0], current_p[1]); 01611 char *tmpstr=new char[100]; 01612 sprintf(tmpstr, "%d", i); 01613 fprintf(f3, "text(%f, %f, '%s');\n", current_p[0], current_p[1]+0.00015, tmpstr); 01614 # ifdef CLC_DBG 01615 printf(" (c3) DUPLICATE: %f", dist_since_last_used_p); 01616 # endif 01617 #endif 01618 01619 // Multiple point, definitely use this one, since we assume it is due to a kink in 01620 // the curve. 01621 skip_quad[quad_indx]=1; 01622 quad_corner_trim_curve[quad_indx].push_back(trim_curve[i]); 01623 quad_corner_trim_curve_p[quad_indx].push_back(current_p); 01624 quad_corners++; 01625 last_dir=current_p-last_used_p; 01626 //printf("\n last_dir 3 =%g %g\n", last_dir[0], last_dir[1]); 01627 last_dir.normalize(); 01628 //printf("\n last_dir 4 =%g %g\n", last_dir[0], last_dir[1]); 01629 last_used_p=current_p; 01630 a_first_point_added=true; 01631 } 01632 // else 01633 // { 01634 // printf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ %f\n", dist_since_last_used_p); 01635 // } 01636 01637 } 01638 } 01639 } // end of if-test ensuring that [u, v] is in the parameter domain. 01640 #ifdef CLC_DBG 01641 printf("\n"); 01642 #endif 01643 } // end of i-loop over points on curve 'curve'. 01644 } // end of 'curve'-loop 01645 01646 #ifdef DBG 01647 fprintf(f1, "hold off\n"); 01648 fprintf(f2, "hold off\n"); 01649 fprintf(f3, "hold off\n"); 01650 fprintf(f4, "hold off\n"); 01651 fclose(f1); 01652 fclose(f2); 01653 fclose(f3); 01654 fclose(f4); 01655 #endif 01656 } 01657 01658 01659 01660 01661 01662 01663 //============================================================================================================== 01664 // 01665 // Main routine returning a trimmed mesh in form of a list of triangles. 01666 // 01667 // 081206: Adding support for more than one curve. 01668 // 01669 //============================================================================================================== 01670 01671 void make_trimmed_mesh(boost::shared_ptr<SplineSurface> srf, 01672 vector<boost::shared_ptr<SplineCurve> >& crv_set, 01673 vector< Vector3D > &vert, 01674 vector< Vector2D > &vert_p, 01675 vector< int > &bd, 01676 vector< Vector3D > &norm, 01677 //vector< Vector3D > &col, // colour // 090130: Not in use 01678 vector<int> &mesh, 01679 vector< Vector3D > &trim_curve_0, // Output 01680 vector< Vector3D > &trim_curve_p_0, // Output 01681 const int dn, // Initially dn+1 nodes in the u-direction 01682 const int dm, // Initially dm+1 nodes in the v-direction 01683 //vector< Vector3D > &extra_v, // 090130: Not in use 01684 double bd_res_ratio) 01685 { 01686 const double duplicate_tolerance = 1e-8; // 100210: Absolute number. Was 1e-12. 01687 01688 vert.resize((dn+1)*(dm+1)); 01689 vert_p.resize((dn+1)*(dm+1)); 01690 bd.resize((dn+1)*(dm+1)); 01691 norm.resize((dn+1)*(dm+1)); 01692 mesh.resize(0); 01693 01694 const double u0=srf->startparam_u(); 01695 const double u1=srf->endparam_u(); 01696 const double v0=srf->startparam_v(); 01697 const double v1=srf->endparam_v(); 01698 int i, j; 01699 const double ustep = (u1 - u0)/(dn-1); 01700 const double vstep = (v1 - v0)/(dm-1); 01701 01702 //-------------------------------------------------------------------------------------------------------------- 01703 // 01704 // Discretizing the trim curve... 01705 // 040621: We use the 'trim_curve_p' to store parameter pairs with z=0.0 in the form of triples for the 01706 // points in 'trim_curve'. 01707 // 081208: Now discretizing a set of curves. 01708 // 01709 //-------------------------------------------------------------------------------------------------------------- 01710 01711 vector< vector<Vector3D> > trim_curve_all(crv_set.size()), trim_curve_p_all(crv_set.size()); 01712 for (int c=0; c<int(crv_set.size()); c++) 01713 { 01714 vector<Vector3D> &trim_curve = trim_curve_all[c]; 01715 vector<Vector3D> &trim_curve_p = trim_curve_p_all[c]; 01716 01717 boost::shared_ptr<SplineCurve> crv = crv_set[c]; 01718 01719 vector<double> corner_pars; 01720 int kk = crv->order(); 01721 int kn = crv->numCoefs(); 01722 std::vector<double>::iterator st = crv->knotsBegin(); 01723 int kstat = 0; 01724 corner_pars.push_back(st[kk-1]); 01725 int knot_ind = kk; 01726 while (knot_ind < kn) 01727 { 01728 int knot_mult = 1; 01729 while (st[knot_ind] == st[knot_ind+knot_mult]) 01730 ++knot_mult; 01731 if (kstat < 0) 01732 { 01733 THROW("Unexpected incident."); 01734 } 01735 if (knot_mult > kk - 2) // A kink 01736 corner_pars.push_back(st[knot_ind]); 01737 01738 knot_ind += knot_mult; 01739 } 01740 01741 corner_pars.push_back(st[kn]); 01742 01743 int n2 = 200/(corner_pars.size() - 1); //200; //std::max(200, 4*n); @@sbr Should be const. 01744 01745 // 100212: For debugging/testing: 01746 // int n2 = 30/(corner_pars.size() - 1); //200; //std::max(200, 4*n); @@sbr Should be const. 01747 01748 for (i = 0; i < int(corner_pars.size()) - 1; ++i) 01749 { 01750 double t0=corner_pars[i]; 01751 double t1=corner_pars[i+1]; 01752 double tstep = (t1 - t0)/(n2-1); 01753 // for (i=0; i<n2; i++) 01754 double t = t0; 01755 double ref_step = tstep; 01756 01757 while (t <= t1) 01758 { 01759 // double t=i/double(n2); // Assuming the curve is closed, so we don't 01760 // need t==1.0... 01761 // t = ( crv->et[crv->ik-1] * (1.0-t) + 01762 // crv->et[crv->in] * t ); 01763 01764 Point result; 01765 Point result2; 01766 crv->point(result, t); 01767 01768 if ((bd_res_ratio > 0.0) && (trim_curve_p.size() > 0) && 01769 ((fabs(trim_curve_p[trim_curve_p.size()-1][0] - result[0]) > bd_res_ratio*ustep) || 01770 (fabs(trim_curve_p[trim_curve_p.size()-1][1] - result[1]) > bd_res_ratio*vstep))) 01771 { 01772 t -= ref_step; // Return to t of last approved parameter pt. 01773 ref_step *= 0.5; 01774 t += ref_step; 01775 } 01776 else 01777 // 090219: We only allow a point to be added if it is distinctly different from the last one 01778 // added! 01779 { 01780 // if ( (trim_curve_p.size()==0) || 01781 // ( ((result[0]-trim_curve_p[trim_curve_p.size()-1][0])* 01782 // (result[1]-trim_curve_p[trim_curve_p.size()-1][1]) ) > 1e-14*1e-14 ) ) 01783 { 01784 trim_curve_p.push_back(Vector3D(result[0], result[1], 0.0)); 01785 // printf("added %g %g\n", result[0], result[1]); 01786 srf->point(result2, result[0], result[1]); 01787 trim_curve.push_back(Vector3D(result2.begin())); 01788 } 01789 if (t == t1) 01790 break; 01791 t += tstep; 01792 if (t > t1) { 01793 tstep = t1 - (t - tstep); 01794 t = t1; 01795 } 01796 ref_step = tstep; 01797 } 01798 // printf("==================== t=%g\n", t); 01799 } 01800 } 01801 } // closing of loop over the curves 01802 01803 01804 01805 #ifdef DBG 01806 // 100222: Just to separate different runs from one another... 01807 printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"); 01808 printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"); 01809 printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"); 01810 #endif 01811 01812 01813 01814 01815 #ifdef DBG 01816 puts("===================================================================================================="); 01817 printf("Weeding out duplicates from contours, tolerance=%e\n", duplicate_tolerance); 01818 #endif 01819 01820 //-------------------------------------------------------------------------------------------------------------- 01821 // 01822 // 090220: Weeding out duplicates. Note that we want the duplicates taken out for the lists used for the 01823 // insideness-testing routines, but we want them in for the corner construction! Hence... making two 01824 // sets of lists... 01825 // 100212: There was a bug here, when more than two points were equal, duplicates would survive. 01826 // 01827 //-------------------------------------------------------------------------------------------------------------- 01828 01829 vector< vector<Vector3D> > trim_curve_all_orig=trim_curve_all, trim_curve_p_all_orig=trim_curve_p_all; 01830 01831 for (int c=0; c<int(crv_set.size()); c++) 01832 { 01833 // printf("Length before: %d %d\n", int(trim_curve_all[c].size()), int(trim_curve_p_all[c].size())); 01834 01835 ASSERT2(trim_curve_all[c].size()>0, printf("Should not be here with degenerate curve.\n")); 01836 ASSERT2(trim_curve_p_all[c].size()>0, printf("Should not be here with degenerate curve.\n")); 01837 01838 vector<Vector3D> trim_curve, trim_curve_p; 01839 trim_curve.push_back(trim_curve_all[c][0]); 01840 trim_curve_p.push_back(trim_curve_p_all[c][0]); 01841 01842 for (int i=1; i<int(trim_curve_all[c].size()); i++) 01843 { 01844 const double last_u = trim_curve_p[trim_curve.size()-1][0]; 01845 const double last_v = trim_curve_p[trim_curve.size()-1][1]; 01846 const double dist_squared = ( (trim_curve_p_all[c][i][0]-last_u) * (trim_curve_p_all[c][i][0]-last_u) + 01847 (trim_curve_p_all[c][i][1]-last_v) * (trim_curve_p_all[c][i][1]-last_v) ); 01848 if ( dist_squared > duplicate_tolerance*duplicate_tolerance ) 01849 { 01850 // printf("%3d:%5d: prev point was (%f, %f), new is (%f, %f), dist=%e, pushing.\n", 01851 // c, i, last_u, last_v, trim_curve_p_all[c][i][0], trim_curve_p_all[c][i][1], sqrt(dist_squared)); 01852 trim_curve.push_back(trim_curve_all[c][i]); 01853 trim_curve_p.push_back(trim_curve_p_all[c][i]); 01854 } 01855 // else 01856 // printf("%3d:%5d: dist = %e, point skipped\n", c, i, sqrt(dist_squared)); 01857 } 01858 01859 // printf("Length after: %d %d\n", int(trim_curve.size()), int(trim_curve_p.size())); 01860 01861 // 100212: Must now handle the first point! It can still be coincident with the last one. 01862 const double last_u = trim_curve_p[trim_curve.size()-1][0]; 01863 const double last_v = trim_curve_p[trim_curve.size()-1][1]; 01864 const double dist_squared = ( (trim_curve_p[0][0]-last_u) * (trim_curve_p[0][0]-last_u) + 01865 (trim_curve_p[0][1]-last_v) * (trim_curve_p[0][1]-last_v) ); 01866 if ( dist_squared <= duplicate_tolerance*duplicate_tolerance ) 01867 { 01868 trim_curve.pop_back(); 01869 trim_curve_p.pop_back(); 01870 } 01871 01872 // Finally, we replace the original curves with the new ones. 01873 trim_curve_p_all[c] = trim_curve_p; 01874 trim_curve_all[c] = trim_curve; 01875 01876 // printf("Length after: %d %d\n", int(trim_curve.size()), int(trim_curve_p.size())); 01877 } 01878 01879 01880 01881 01882 //-------------------------------------------------------------------------------------------------------------- 01883 // 01884 // 100212: A check to make sure we have no duplicates. 01885 // 01886 //-------------------------------------------------------------------------------------------------------------- 01887 01888 { 01889 bool duplicates_found = false; 01890 for (int c=0; c<int(crv_set.size()); c++) 01891 { 01892 vector<Vector3D> &trim_curve = trim_curve_all[c]; 01893 vector<Vector3D> &trim_curve_p = trim_curve_p_all[c]; 01894 01895 for (int i=0; i<int(trim_curve.size()); i++) 01896 { 01897 const int im1 = (i-1+trim_curve.size())%trim_curve.size(); 01898 const double dist_to_previous = 01899 sqrt( (trim_curve_p[i][0]-trim_curve_p[im1][0])*(trim_curve_p[i][0]-trim_curve_p[im1][0]) + 01900 (trim_curve_p[i][1]-trim_curve_p[im1][1])*(trim_curve_p[i][1]-trim_curve_p[im1][1]) ); 01901 01902 // printf("%2d:%5d: dist = %e\n", c, i, dist_to_previous); 01903 if ( dist_to_previous <= duplicate_tolerance ) 01904 { 01905 #ifdef DBG 01906 printf("%2d:%5d: DUPLICATE POINT! dist_to_previous = %e\n", c, i, dist_to_previous); 01907 #endif 01908 duplicates_found = true; 01909 } 01910 } 01911 } 01912 if (duplicates_found) 01913 MESSAGE("\n\n WARNING: Something is fishy. At this point all duplicates should have been eliminated.\n" 01914 " This may lead to a funny mesh...\n\n"); 01915 } 01916 01917 01918 01919 01920 //-------------------------------------------------------------------------------------------------------------- 01921 // 01922 // Setting up the structure for fast inside-testing, and the 'contour' cursor array. 01923 // 081208: Again, extending to a set of contours... 01924 // 100212: The "fast inside-testing" stuff is long gone. 01925 // 01926 //-------------------------------------------------------------------------------------------------------------- 01927 01928 vector< vector<int> > contour_all(crv_set.size()); 01929 for (int c=0; c<int(crv_set.size()); c++) 01930 { 01931 vector<int> &contour = contour_all[c]; 01932 const vector<Vector3D> &trim_curve_p = trim_curve_p_all[c]; 01933 01934 contour.resize(trim_curve_p.size()); 01935 for (i=0; i<(int)trim_curve_p.size(); i++) 01936 contour[i]=3*i; 01937 } 01938 01939 #ifdef DBG 01940 // Just printing out matlab-code for rendering of all contours. 01941 if (1) 01942 for (int c=0; c<int(crv_set.size()); c++) 01943 { 01944 const vector<int> &contour = contour_all[c]; 01945 const vector<Vector3D> &trim_curve_p = trim_curve_p_all[c]; 01946 01947 char *tmp=new char[1000]; 01948 sprintf(tmp, "d%d.m", c); 01949 FILE *f=fopen(tmp, "w"); 01950 for (int i=0; i<int(trim_curve_p.size())-1; i++) 01951 { 01952 const double &x0 = trim_curve_p[contour[i]/3 ][0], &y0=trim_curve_p[contour[i]/3 ][1]; 01953 const double &x1 = trim_curve_p[contour[i]/3+1][0], &y1=trim_curve_p[contour[i]/3+1][1]; 01954 fprintf(f, "line([%f; %f], [%f; %f], 'color', 'b', 'linewidth', 3);\n", x0, x1, y0, y1); 01955 } 01956 fclose(f); 01957 } 01958 #endif 01959 01960 01961 01962 //-------------------------------------------------------------------------------------------------------------- 01963 // 01964 // Cache'ing all 'is_inside' results... 01965 // 081208: Again, extending to a set of contours... 01966 // 01967 // 100210: It seems like a (harmless) "bug" that the 'vert-arrays are filled once for every curve. 01968 // 01969 //-------------------------------------------------------------------------------------------------------------- 01970 01971 vector< vector<int> > inside_all(crv_set.size()); 01972 for (int c=0; c<int(crv_set.size()); c++) 01973 { 01974 const vector<int> &contour = contour_all[c]; 01975 const vector<Vector3D> &trim_curve_p = trim_curve_p_all[c]; 01976 vector<int> &inside = inside_all[c]; 01977 inside = vector<int>((dn+1)*(dm+1)); 01978 double uv[2], s; 01979 vector<Point> res(3); 01980 Point nrm; 01981 01982 ASSERT2(srf->dimension()==3, printf("Huh?! dim=%d\n", srf->dimension())); 01983 01984 for (i=0; i<=dm; i++) 01985 { 01986 double t=i/double(dm); 01987 uv[1]=v0*(1.0-t) + v1*t; 01988 for (j=0; j<=dn; j++) 01989 { 01990 s=j/double(dn); 01991 uv[0]=u0*(1.0-s) + u1*s; 01992 01993 inside[i*(dn+1)+j] = is_inside(trim_curve_p, contour, uv[0], uv[1]); 01994 01995 // 090115: 01996 if (s2m_with_boundary) 01997 inside[i*(dn+1)+j] |= is_on_contour(trim_curve_p, contour, uv[0], uv[1]); 01998 01999 // 090203: Enable this to see the trimming curve inside the untrimmed surface, for debugging purposes. 02000 // 100210: This does not work, or this enabling is not enough. 02001 // inside[i*(dn+1)+j] = 1; // !!! 02002 02003 srf->point(res, uv[0], uv[1], 1); 02004 nrm = res[1].cross(res[2]); 02005 02006 // 100210: Why on earth is this done for every curve?! 02007 vert[i*(dn+1)+j] = Vector3D(res[0].begin()); 02008 vert_p[i*(dn+1)+j] = Vector2D(uv); 02009 02010 bd[i*(dn+1)+j]= (i==0 || i==dm || j==0 || j==dn) ? 1 : 0; 02011 if (nrm.length() < 1.0e-12) 02012 nrm.setValue(0.0, 0.0, 1.0); 02013 norm[i*(dn+1)+j] = Vector3D(nrm.begin()); 02014 norm[i*(dn+1)+j].normalize(); 02015 } 02016 } 02017 } 02018 02019 02020 #ifdef DBG 02021 { 02022 FILE *f5=fopen("r1.m", "w"); // regular mesh 02023 fprintf(f5, "hold on\ngrid on\n"); 02024 for (int i=0; i<dm; i++) 02025 for (int j=0; j<=dn; j++) 02026 fprintf(f5, "line([%f; %f], [%f; %f], 'color', 'k');\n", 02027 vert_p[i*(dn+1)+j][0], vert_p[(i+1)*(dn+1)+j][0], 02028 vert_p[i*(dn+1)+j][1], vert_p[(i+1)*(dn+1)+j][1]); 02029 for (int i=0; i<=dm; i++) 02030 for (int j=0; j<dn; j++) 02031 fprintf(f5, "line([%f; %f], [%f; %f], 'color', 'k');\n", 02032 vert_p[i*(dn+1)+j][0], vert_p[i*(dn+1)+j+1][0], 02033 vert_p[i*(dn+1)+j][1], vert_p[i*(dn+1)+j+1][1]); 02034 fprintf(f5, "hold off\n"); 02035 fclose(f5); 02036 } 02037 #endif 02038 02039 02040 02041 //-------------------------------------------------------------------------------------------------------------- 02042 // 02043 // 090130: We want to preserve corner points on the contour. The idea: Identify all quads containing such 02044 // corners, then treat these quads separately: Split in the first corner, producing four new 02045 // triangles. These must then be processed further, until they contain no corner points. Then these 02046 // triangles must still be processed the "old way". 02047 // 02048 // 090203: Note that by checking that we do not make degenerate primitives, it does not become as important 02049 // to check for duplicate corner points any more. 02050 // 02051 //-------------------------------------------------------------------------------------------------------------- 02052 02053 02054 const double pt_edge_degen_limit=1e-6; // 090203: Should really be relative to the extent of the 02055 // geometry. Used to avoid constructing 02056 // almost-degenerate polygons. Applied in the 02057 // parameter domain. 02058 02059 // const double cosangle_limit= // 090203: The angle must be greater than this for a "corner" 02060 // cos(M_PI/180.0); 02061 const double cosangle_limit= // 090203: The angle must be greater than this for a "corner" 02062 cos(M_PI/180.0*0.1); // to be defined as a corner, i.e., for splitting to 02063 // be done. 02064 02065 const double pt_mult_def=1e-14; // 090203: In the parameter domain, two points closer than 02066 // this are identified. 02067 02068 const double eps=1e-13; // 090203: Used for zero-tests for distances in the parameter domain. 02069 02070 //const double tau=1e-12; // 090203: Used for zero-tests for distances in the geometric space. 02071 02072 const double max_corner_dist= // 090218: When corners are this far away, they are used no matter what 02073 100.0*std::max((u1-u0)/dn, (v1-v0)/dm); // the angle is. Must be larger than 'min_corner_dist'. 02074 02075 const double min_corner_dist= // 090219: When corners are this close, they are not used no 02076 0.5*std::min((u1-u0)/dn, (v1-v0)/dm); // matter what the angle is. 02077 02078 02079 02080 02081 02082 #ifdef DBG 02083 puts("===================================================================================================="); 02084 puts("Generating list of corners to force splitting in."); 02085 #endif 02086 02087 02088 02089 // 090218: Each of these will have dim nxn. 02090 vector<int> skip_quad; // Contains a 1 if quad is to be skipped, i.e., it is 02091 // replaced by custom list of triangles, since it did 02092 // contain corners or such. 02093 02094 vector< vector<Vector3D> > quad_corner_trim_curve; // The 3D corner point, fetched from the proper 02095 // trimming curve. Is this really used afterward? 02096 // Yes. Need it for updating 'vert'. 02097 02098 vector< vector<Vector2D> > quad_corner_trim_curve_p; // The 2D parameter corner point, constructed from a 02099 // 3D point by ignoring the third coordinate. Fetched 02100 // from the proper trimming curve. Not smart to call 02101 // it '...trim_curve...', should have been just 02102 // 'quad_corner_p' or something... 02103 02104 int quad_corners=0; 02105 construct_corner_lists(srf, crv_set, vert, vert_p, norm, mesh, 02106 trim_curve_all_orig, trim_curve_p_all_orig, 02107 dn, dm, bd_res_ratio, eps, cosangle_limit, pt_edge_degen_limit, pt_mult_def, min_corner_dist, 02108 max_corner_dist, skip_quad, quad_corner_trim_curve, quad_corner_trim_curve_p, quad_corners); 02109 02110 02111 02112 02113 02114 //-------------------------------------------------------------------------------------------------------------- 02115 // 02116 // 100210: Ok, in this stage, quads containing corners are split. They are initially split into four 02117 // triangles, which may again be further split in the innermost loop. This is to take care of cases 02118 // where quads contain more than just one "corner points". 02119 // 02120 //-------------------------------------------------------------------------------------------------------------- 02121 02122 #ifdef CORNER_SPLITTING 02123 02124 #ifdef DBG 02125 puts("===================================================================================================="); 02126 printf("\"Pre-splitting\" mesh, %d corners was found.\n", quad_corners); 02127 #endif 02128 02129 // 02130 // 090206: We build structures like 'quad_corner' for needed data. The problem with the "merged list" 02131 // approach is that information for each curve, like where it ends, must be preserved. This is 02132 // because the curves are closed. 02133 // 02134 02135 // Now splitting first quads, then triangles. 02136 02137 #ifdef DBG 02138 printf("dn=%d, dm=%d\n", dn, dm); 02139 #endif 02140 for (int i=0; i<dn*dm; i++) 02141 { 02142 // printf("i=%3d: corners: %d\n", i, int(quad_corner[i].size())); 02143 if (quad_corner_trim_curve_p[i].size()>0) 02144 { 02145 // First, split the quad. 02146 02147 const int p=i/dn, q=i%dn; 02148 02149 vert.push_back(quad_corner_trim_curve[i][0]); 02150 vert_p.push_back(quad_corner_trim_curve_p[i][0]); 02151 bd.push_back(1); 02152 // Unfortunately, we don't have any normals readily available. What we can do, is to just take the 02153 // average of the quad. What we *should* do, is to evaluate the surface... 02154 { 02155 Vector3D new_norm = norm[p*(dn+1)+q] + norm[p*(dn+1)+q+1] + norm[(p+1)*(dn+1)+q] + norm[(p+1)*(dn+1)+q+1]; 02156 new_norm.normalize(); 02157 norm.push_back(new_norm); 02158 } 02159 02160 int first_new_triangle=mesh.size()/3; 02161 02162 // printf("i=%d, p=%d, q=%d, ndx=%d: splitting this quad in %g, %g.\n", i, p, q, p*(n+1)+q, 02163 // vert_p[vert_p.size()-1][0], vert_p[vert_p.size()-1][1]); 02164 02165 add_triangle(vert_p, mesh, vert.size()-1, (p+1)*(dn+1)+q , p *(dn+1)+q ); 02166 add_triangle(vert_p, mesh, vert.size()-1, (p+1)*(dn+1)+q+1, (p+1)*(dn+1)+q ); 02167 add_triangle(vert_p, mesh, vert.size()-1, p *(dn+1)+q+1, (p+1)*(dn+1)+q+1); 02168 add_triangle(vert_p, mesh, vert.size()-1, p *(dn+1)+q , p *(dn+1)+q+1); 02169 02170 // Now, if there are more corners, we have to go through the triangles. And we have to do this again 02171 // for each new corner, and apply the split-testing and eventual splitting to all newly generated 02172 // triangles also. We assume that only one triangle can contain a given corner, so when a corner is 02173 // used, it does not need to be tested for again. 02174 02175 if (int(quad_corner_trim_curve_p[i].size())>1) 02176 { 02177 02178 //printf(" Corners in total: %d, continuing with triangle-splitting.\n", int(quad_corner[i].size())); 02179 02180 for (int j=1; j<int(quad_corner_trim_curve_p[i].size()); j++) 02181 { 02182 bool done_splitting=false; 02183 // 100210: Seems to be enough to enable this to skip the triangle-splitting stage. 02184 // done_splitting=true; 02185 02186 Vector2D &corner = quad_corner_trim_curve_p[i][j]; 02187 //printf(" j=%d, checking corner %f, %f.\n", j, corner[0], corner[1]); 02188 02189 // Ok, must now check, and possibly split triangles from 'first_new_triangle'. 02190 // Note that we quit after we split, since we want to start over again asap. 02191 for (int k=first_new_triangle; (k<int(mesh.size())/3) && (!done_splitting); k++) 02192 { 02193 //printf(" k=%d\n", k); 02194 if ( 02195 pt_inside_tri(corner, vert_p[mesh[3*k]], vert_p[mesh[3*k+1]], vert_p[mesh[3*k+2]]) 02196 // && (pt_dist_from_tri_edge(corner, vert_p[mesh[3*k]], 02197 // vert_p[mesh[3*k+1]], vert_p[mesh[3*k+2]]) 02198 // > pt_edge_degen_limit) 02199 // 090204: This test caused corners on the outer boundary not to trigger 02200 // splitting. Enabling this test is probably not a good idea, unless these 02201 // splitting causes other problems... 02202 ) 02203 // Splitting the triangle. 02204 { 02205 //printf(" splitting, i=%d, j=%d, k=%d.\n", i, j, k); 02206 vert.push_back(quad_corner_trim_curve[i][j]); 02207 vert_p.push_back(corner); 02208 bd.push_back(1); 02209 const int a=mesh[3*k], b=mesh[3*k+1], c=mesh[3*k+2]; 02210 // Unfortunately, we don't have any normals readily available. What we can do, is 02211 // to just take the average of the quad. What we should do, is to evaluate the 02212 // surface... 02213 { 02214 Vector3D tmp = norm[a] + norm[b] + norm[c]; 02215 tmp.normalize(); 02216 norm.push_back(tmp); 02217 } 02218 mesh.erase(mesh.begin()+3*k, mesh.begin()+3*k+3); 02219 add_triangle(vert_p, mesh, vert.size()-1, a, b); 02220 add_triangle(vert_p, mesh, vert.size()-1, b, c); 02221 add_triangle(vert_p, mesh, vert.size()-1, c, a); 02222 done_splitting=true; // No need to check more triangles for a corner already used. 02223 // 090204: Hmm... maybe not entirely true. If splitting was done on an edge, there may 02224 // be a triangle on the other side of that edge that also needs splitting...?! 02225 } 02226 } 02227 } 02228 } 02229 } 02230 } 02231 02232 #endif // end of #ifdef CORNER_SPLITTING 02233 02234 #ifdef DBG 02235 vector<char> col(mesh.size()/3, 'g'); 02236 #endif 02237 02238 02239 02240 02241 #ifdef SPLIT_NON_CORNER_QUADS 02242 02243 #ifdef DBG 02244 puts("===================================================================================================="); 02245 puts("Splitting remaining quads according to the first trimming curve."); 02246 #endif 02247 02248 for (i=0; i<dm; i++) 02249 for (j=0; j<dn; j++) 02250 if (!skip_quad[i*dn+j]) 02251 { 02252 // if ( (inside1[i*(n+1)+j]!=0) && (inside1[i*(n+1)+j]!=1) ) 02253 // MESSAGE("noe er muffens, inside1 er hverken 0 eller 1?!"); 02254 // if ( (inside[i*(n+1)+j]!=0) && (inside[i*(n+1)+j]!=1) ) 02255 // MESSAGE("noe er muffens, inside er hverken 0 eller 1?!"); 02256 02257 int s=((inside_all[0][ i *(dn+1)+j ] << 3) + // 8 02258 (inside_all[0][ i *(dn+1)+j+1] << 2) + // 4 02259 (inside_all[0][(i+1)*(dn+1)+j ] << 1) + // 2 02260 (inside_all[0][(i+1)*(dn+1)+j+1] << 0)); // 1 02261 // if (s!=0) // ((i==6) && (j==17)) 02262 // printf("s=%d, i=%d, j=%d\n", s, i, j); 02263 02264 // 090216: keeping all to test if the missing triangle was made at all. Ok, it was made, but then wrongly split here. 02265 //s=15; 02266 02267 // 090216: checking the lower right corner for debugging of problematic quad... 02268 //if ((fabs(vert_p[i*(n+1)+j+1][0]-.01)<1e-5) && (fabs(vert_p[i*(n+1)+j+1][1]-0.65)<1e-5)) 02269 // if ((i==12) && (j==3)) 02270 // printf("oooooooooooooooooooooooooooooooooooooooooooooooooo> s=%d\n", s), s=15; 02271 02272 // 100223: Shit! Didn't realize this override was here... This means that all the quad-splitting 02273 // code below is actually not used... Wonder when this was done... 02274 s=15; 02275 02276 switch (s) 02277 { 02278 case 15: 02279 // The whole rectangle is inside. 02280 { 02281 add_triangle(vert_p, mesh, (i+1)*(dn+1)+j, i*(dn+1)+j , i *(dn+1)+j+1); 02282 add_triangle(vert_p, mesh, (i+1)*(dn+1)+j, i*(dn+1)+j+1, (i+1)*(dn+1)+j+1); 02283 } 02284 break; 02285 //------------------------------------------------------------------- 02286 case 10: 02287 // The two left (smaller u, smaller j) points are inside. 02288 { 02289 // first segm this col 02290 vector<int> si0; si0.push_back(i ); si0.push_back(i+1); 02291 vector<int> sj0; sj0.push_back(j ); sj0.push_back(j ); 02292 vector<int> si1; si1.push_back(i ); si1.push_back(i+1); 02293 vector<int> sj1; sj1.push_back(j+1); sj1.push_back(j+1); 02294 // 02295 // s0[1] B s1[1] (Husk: i+1 er her!) 02296 // o----x------o 02297 // | | 02298 // | | 02299 // | | 02300 // o----x------o 02301 // s0[0] A s1[0] (Husk: i er her!) 02302 // 02303 split_quad(i+1, j, // s0[1] 02304 i, j, // s0[0] 02305 -1, -1, 02306 // Then the two intersection points are 02307 si0, sj0, si1, // somewhere on s0[0]-s1[0] (A) and 02308 sj1, // somewhere on s0[1]-s1[1] (B). 02309 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02310 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02311 } 02312 break; 02313 //------------------------------------------------------------------- 02314 case 5: 02315 // The two right (larger u, larger j) points are inside. 02316 { 02317 // first segm this col 02318 vector<int> si0; si0.push_back(i+1); si0.push_back(i ); 02319 vector<int> sj0; sj0.push_back(j ); sj0.push_back(j ); 02320 vector<int> si1; si1.push_back(i+1); si1.push_back(i ); 02321 vector<int> sj1; sj1.push_back(j+1); sj1.push_back(j+1); 02322 // 02323 // s0[0] A s1[0] (Husk: i+1 er her!) 02324 // o----x------o 02325 // | | 02326 // | | 02327 // | | 02328 // o----x------o 02329 // s0[1] B s1[1] (Husk: i er her!) 02330 // 02331 split_quad(i, j+1, // s1[1] 02332 i+1, j+1, // s1[0] 02333 -1, -1, 02334 // Then the two intersection points are 02335 si0, sj0, si1, // somewhere on s0[0]-s1[0] (A) and 02336 sj1, // somewhere on s0[1]-s1[1] (B). 02337 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02338 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02339 } 02340 break; 02341 //------------------------------------------------------------------- 02342 case 12: 02343 // The two lower (smaller v, smaller i) points are inside. 02344 { 02345 // (si0[0], sj0[0]) - (si1[0], sj1[0]) is the first segment, etc. 02346 // first segm this col 02347 vector<int> si0; si0.push_back(i ); si0.push_back(i+1); 02348 vector<int> sj0; sj0.push_back(j+1); sj0.push_back(j ); 02349 vector<int> si1; si1.push_back(i+1); si1.push_back(i ); 02350 vector<int> sj1; sj1.push_back(j+1); sj1.push_back(j ); 02351 // 02352 // s0[1] s1[0] (Husk: i+1 er her!) 02353 // o-----------o 02354 // | | 02355 // B x x A 02356 // | | 02357 // o-----------o 02358 // s1[1] s0[0] (Husk: i er her!) 02359 // 02360 split_quad(i, j, // s1[1] 02361 i, j+1, // s0[0] 02362 -1, -1, 02363 // Then the two intersection points are 02364 si0, sj0, si1, // somewhere on s0[0]-s1[0] (A) and 02365 sj1, // somewhere on s0[1]-s1[1] (B). 02366 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02367 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02368 } 02369 break; 02370 //------------------------------------------------------------------- 02371 case 1: 02372 // The upper, right (larger i, larger j) point is inside. 02373 { 02374 // first segm this col 02375 vector<int> si0; si0.push_back(i+1); si0.push_back(i+1); 02376 vector<int> sj0; sj0.push_back(j+1); sj0.push_back(j+1); 02377 vector<int> si1; si1.push_back(i+1); si1.push_back(i ); 02378 vector<int> sj1; sj1.push_back(j ); sj1.push_back(j+1); 02379 // 02380 // s1[0] A s0[0]=s0[1] 02381 // o-------x---o 02382 // | \ | 02383 // | \ x B 02384 // | | 02385 // o-----------o 02386 // s1[1] 02387 // 02388 split_quad(i+1, j+1, // s0[0] 02389 -1, -1, -1, -1, 02390 // Then the two intersection points are 02391 si0, sj0, si1, // somewhere on s0[0]-s1[0] (A) and 02392 sj1, // somewhere on s0[1]-s1[1] (B). 02393 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02394 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02395 } 02396 break; 02397 //------------------------------------------------------------------- 02398 case 4: 02399 // The lower, right (smaller i, larger j) point is inside. 02400 { 02401 // first segm this col 02402 vector<int> si0; si0.push_back(i+1); si0.push_back(i ); 02403 vector<int> sj0; sj0.push_back(j+1); sj0.push_back(j ); 02404 vector<int> si1; si1.push_back(i ); si1.push_back(i ); 02405 vector<int> sj1; sj1.push_back(j+1); sj1.push_back(j+1); 02406 split_quad(i, j+1, 02407 -1, -1, -1, -1, 02408 // Then the two intersection points are 02409 si0, sj0, si1, // somewhere on s0[0]-s1[0] (A) and 02410 sj1, // somewhere on s0[1]-s1[1] (B). 02411 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02412 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02413 } 02414 break; 02415 //------------------------------------------------------------------- 02416 case 8: 02417 // The lower, left (smaller i, smaller j) point is inside. 02418 { 02419 // first segm this col 02420 vector<int> si0; si0.push_back(i ); si0.push_back(i ); 02421 vector<int> sj0; sj0.push_back(j ); sj0.push_back(j ); 02422 vector<int> si1; si1.push_back(i ); si1.push_back(i+1); 02423 vector<int> sj1; sj1.push_back(j+1); sj1.push_back(j ); 02424 split_quad(i, j, 02425 -1, -1, -1, -1, 02426 // Then the two intersection points are 02427 si0, sj0, si1, // somewhere on s0[0]-s1[0] (A) and 02428 sj1, // somewhere on s0[1]-s1[1] (B). 02429 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02430 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02431 } 02432 break; 02433 //------------------------------------------------------------------- 02434 case 2: 02435 // The upper, left (larger i, smaller j) point is inside. 02436 { 02437 // first segm this col 02438 vector<int> si0; si0.push_back(i+1); si0.push_back(i+1); 02439 vector<int> sj0; sj0.push_back(j ); sj0.push_back(j ); 02440 vector<int> si1; si1.push_back(i ); si1.push_back(i+1); 02441 vector<int> sj1; sj1.push_back(j ); sj1.push_back(j+1); 02442 split_quad(i+1, j, 02443 -1, -1, -1, -1, 02444 // Then the two intersection points are 02445 si0, sj0, si1, // somewhere on s0[0]-s1[0] (A) and 02446 sj1, // somewhere on s0[1]-s1[1] (B). 02447 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02448 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02449 } 02450 break; 02451 //------------------------------------------------------------------- 02452 case 3: 02453 // The two upper points are inside. 02454 { 02455 // Same as for 12, but swapped the two segments. 02456 // Also using (i+1, j+1) and (i+1, j) as the inside points. 02457 vector<int> si0; si0.push_back(i+1); si0.push_back(i ); 02458 vector<int> sj0; sj0.push_back(j ); sj0.push_back(j+1); 02459 vector<int> si1; si1.push_back(i ); si1.push_back(i+1); 02460 vector<int> sj1; sj1.push_back(j ); sj1.push_back(j+1); 02461 split_quad(i+1, j+1, // s1[1] 02462 i+1, j, // s0[0] 02463 -1, -1, 02464 // Then the two intersection points are 02465 si0, sj0, si1, // somewhere on s0[0]-s1[0] (A) and 02466 sj1, // somewhere on s0[1]-s1[1] (B). 02467 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02468 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02469 } 02470 break; 02471 //------------------------------------------------------------------- 02472 case 14: 02473 // All except the upper right point are inside. 02474 { 02475 // first segm this col 02476 vector<int> si0; si0.push_back(i+1); si0.push_back(i+1); 02477 vector<int> sj0; sj0.push_back(j+1); sj0.push_back(j+1); 02478 vector<int> si1; si1.push_back(i ); si1.push_back(i+1); 02479 vector<int> sj1; sj1.push_back(j+1); sj1.push_back(j ); 02480 split_quad(i, j, // 02481 i+1, j, // 02482 i, j+1, // 02483 // Then the two intersection points are 02484 si0, sj0, // somewhere on s0[0]-s1[0] (A) and 02485 si1, sj1, // somewhere on s0[1]-s1[1] (B). 02486 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02487 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02488 } 02489 break; 02490 //------------------------------------------------------------------- 02491 case 11: 02492 // All except the lower right point are inside. 02493 { 02494 // first segm this col 02495 vector<int> si0; si0.push_back(i ); si0.push_back(i ); 02496 vector<int> sj0; sj0.push_back(j+1); sj0.push_back(j+1); 02497 vector<int> si1; si1.push_back(i ); si1.push_back(i+1); 02498 vector<int> sj1; sj1.push_back(j ); sj1.push_back(j+1); 02499 split_quad(i+1, j, // 02500 i+1, j+1, // 02501 i, j, // 02502 // Then the two intersection points are 02503 si0, sj0, // somewhere on s0[0]-s1[0] (A) and 02504 si1, sj1, // somewhere on s0[1]-s1[1] (B). 02505 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02506 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02507 } 02508 break; 02509 //------------------------------------------------------------------- 02510 case 7: 02511 // All except the lower left point are inside. 02512 { 02513 // first segm this col 02514 vector<int> si0; si0.push_back(i ); si0.push_back(i ); 02515 vector<int> sj0; sj0.push_back(j ); sj0.push_back(j ); 02516 vector<int> si1; si1.push_back(i+1); si1.push_back(i ); 02517 vector<int> sj1; sj1.push_back(j ); sj1.push_back(j+1); 02518 split_quad(i+1, j+1, // 02519 i, j+1, // 02520 i+1, j, // 02521 // Then the two intersection points are 02522 si0, sj0, // somewhere on s0[0]-s1[0] (A) and 02523 si1, sj1, // somewhere on s0[1]-s1[1] (B). 02524 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02525 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02526 } 02527 break; 02528 //------------------------------------------------------------------- 02529 case 13: 02530 // All except the upper left point are inside. 02531 { 02532 // first segm this col 02533 vector<int> si0; si0.push_back(i+1); si0.push_back(i+1); 02534 vector<int> sj0; sj0.push_back(j ); sj0.push_back(j ); 02535 vector<int> si1; si1.push_back(i+1); si1.push_back(i ); 02536 vector<int> sj1; sj1.push_back(j+1); sj1.push_back(j ); 02537 split_quad(i, j+1, // 02538 i, j, // 02539 i+1, j+1, // 02540 // Then the two intersection points are 02541 si0, sj0, // somewhere on s0[0]-s1[0] (A) and 02542 si1, sj1, // somewhere on s0[1]-s1[1] (B). 02543 srf, vert, vert_p, bd, norm, mesh, dn, dm, 02544 trim_curve_p_all[0], contour_all[0], s2m_with_boundary); 02545 } 02546 break; 02547 //------------------------------------------------------------------- 02548 } 02549 } 02550 02551 #ifdef DBG 02552 while (col.size()<mesh.size()/3) 02553 col.push_back('r'); 02554 #endif 02555 02556 #endif // end of #ifdef SPLIT_NON_CORNER_QUADS 02557 02558 02559 02560 02561 02562 #ifdef SPLIT_TRIANGLES_AT_OUTER_CURVE 02563 02564 # ifdef DBG 02565 puts("===================================================================================================="); 02566 puts("Splitting all triangles just generated, according to the first trimming curve."); 02567 # endif 02568 02569 // 02570 // 090203: Here we need a new pass to split triangles produced by the new "corner-stuff", and now we should 02571 // split according to the first trimming curve, only. (Note, this all is based on the assumption 02572 // that we should keep whatever is *inside* curve 0, and *outside* curves 1, 2, ...) 02573 // 02574 // 090203: Hmm... maybe this code can somehow be combined with the stage below, for more compact code?! 02575 // 02576 // 090203: Have a suspicion that h.g2 produces a denegerate triangle in this part of the code. 02577 // 02578 02579 for (int c=0; c<std::min(1, int(crv_set.size())); c++) 02580 trim_a_triangle_soup(srf, vert, vert_p, bd, norm, trim_curve_p_all[c], contour_all[c], mesh 02581 # ifdef DBG 02582 , col 02583 # endif 02584 ); 02585 02586 #endif // end of #ifdef SPLIT_TRIANGLES_AT_OUTER_CURVE 02587 02588 02589 02590 02591 02592 #ifdef SPLIT_TRIANGLES_AT_INNER_CURVES 02593 02594 #ifdef DBG 02595 puts("===================================================================================================="); 02596 puts("Splitting triangles according to the second and further trimming curves."); 02597 #endif 02598 02599 //printf("Triangles added in 'make_trimmed_mesh': %d\n", int(mesh.size())/3); 02600 02601 // 02602 // 081208: Adding a post-processing step in which we take into consideration further trimming curves. Now we 02603 // have a mesh consisting of a number of triangles. We loop through all of them, and split them 02604 // similarly to what was done with the first trimming curve above. 02605 // 02606 // (For max code reuse, clarity and orthogonality, we should really replace the above, older code, 02607 // with a pass similar to the one following below.) 02608 // 02609 // First a version without the cached inside-results, since these are kind of messed up now, when we 02610 // must free us from the regular grid... can be fixed later... 02611 // 02612 // 081210: What the h... did I mean by this? Shouldn't using the cached values be 02613 // straightforward?! 02614 // 02615 // 081210: Note that the number of vertices are not reduced, even when it could be. Should be relatively 02616 // easy to add a pass for doing such reduction later, if needed... 02617 // 02618 02619 #ifdef DBG 02620 printf("\n\n================================================== crv_set.size=%d\n\n\n", int(crv_set.size())); 02621 #endif 02622 02623 for (int c=1; c<int(crv_set.size()); c++) 02624 //for (int c=3; c<5; c++) 02625 { 02626 #ifdef DBG 02627 printf("----------==========########## trimming against inner curve %d (range=[%d, %d]) ##########==========----------\n", c, 1, int(crv_set.size())-1); 02628 #endif 02629 trim_a_triangle_soup(srf, vert, vert_p, bd, norm, trim_curve_p_all[c], contour_all[c], mesh, 02630 #ifdef DBG 02631 col, 02632 #endif 02633 true); 02634 } 02635 02636 #endif // end of #ifdef SPLIT_TRIANGLES_AT_INNER_CURVES 02637 02638 02639 02640 02641 #if 0 02642 02643 # if 1 02644 02645 // 100210: What was this?!?! probably some debugging... Seems to have been used for some hardcoded inner 02646 // trimming curves... (number 2 and 3..) 02647 02648 trim_a_triangle_soup(srf, vert, vert_p, bd, norm, trim_curve_p_all[2], contour_all[2], mesh, col, true); 02649 { 02650 FILE *f=fopen("e1.m", "w"); 02651 for (int i=0; i<int(trim_curve_p_all[2].size())-1; i++) 02652 { 02653 const double &x0=trim_curve_p_all[2][contour_all[2][i]/3 ][0]; 02654 const double &y0=trim_curve_p_all[2][contour_all[2][i]/3 ][1]; 02655 const double &x1=trim_curve_p_all[2][contour_all[2][i]/3+1][0]; 02656 const double &y1=trim_curve_p_all[2][contour_all[2][i]/3+1][1]; 02657 fprintf(f, "line([%f; %f], [%f; %f], 'color', 'r', 'linewidth', 3);\n", x0, x1, y0, y1); 02658 } 02659 fclose(f); 02660 } 02661 # endif 02662 02663 puts("####################################################################################################"); 02664 puts("####################################################################################################"); 02665 puts("####################################################################################################"); 02666 02667 02668 # if 0 02669 trim_a_triangle_soup(srf, vert, vert_p, bd, norm, trim_curve_p_all[3], contour_all[3], mesh, col, true); 02670 { 02671 FILE *f=fopen("e2.m", "w"); 02672 for (int i=0; i<int(trim_curve_p_all[3].size())-1; i++) 02673 { 02674 const double &x0=trim_curve_p_all[3][contour_all[3][i]/3 ][0]; 02675 const double &y0=trim_curve_p_all[3][contour_all[3][i]/3 ][1]; 02676 const double &x1=trim_curve_p_all[3][contour_all[3][i]/3+1][0]; 02677 const double &y1=trim_curve_p_all[3][contour_all[3][i]/3+1][1]; 02678 fprintf(f, "line([%f; %f], [%f; %f], 'color', 'r', 'linewidth', 3);\n", x0, x1, y0, y1); 02679 } 02680 fclose(f); 02681 } 02682 # endif 02683 02684 #endif 02685 02686 02687 02688 02689 02690 02691 #ifdef DBG 02692 02693 //-------------------------------------------------------------------------------------------------------------- 02694 // 02695 // 100210: Showing all remaining triangles. Colour coding: 02696 // 02697 // yellow: 02698 // green: 02699 // red: 02700 // cyan: 02701 // magenta: 02702 // 02703 //-------------------------------------------------------------------------------------------------------------- 02704 02705 FILE *f=fopen("vis.m", "w"); 02706 fprintf(f, "clf\n"); 02707 for (int i=0; i<int(mesh.size())/3; i++) 02708 { 02709 const int u=0, v=1; 02710 #if 0 02711 if ( (vert_p[mesh[3*i]][u] > -40) && 02712 (vert_p[mesh[3*i]][u] < -10) && 02713 (vert_p[mesh[3*i]][v] > 25) && 02714 (vert_p[mesh[3*i]][v] < 60) ) 02715 #endif 02716 { 02717 fprintf(f, "patch([%f; %f; %f], [%f; %f; %f], '%c');\n", 02718 vert_p[mesh[3*i]][u], vert_p[mesh[3*i+1]][u], vert_p[mesh[3*i+2]][u], 02719 vert_p[mesh[3*i]][v], vert_p[mesh[3*i+1]][v], vert_p[mesh[3*i+2]][v], col[i]); 02720 fprintf(f, "line([%f; %f; %f; %f], [%f; %f; %f; %f], 'color', 'k');\n", 02721 vert_p[mesh[3*i]][u], vert_p[mesh[3*i+1]][u], vert_p[mesh[3*i+2]][u], vert_p[mesh[3*i]][u], 02722 vert_p[mesh[3*i]][v], vert_p[mesh[3*i+1]][v], vert_p[mesh[3*i+2]][v], vert_p[mesh[3*i]][v]); 02723 } 02724 } 02725 fprintf(f, "zoom on\n"); 02726 fprintf(f, "grid on\n"); 02727 // fprintf(f, "axis('equal');\n"); 02728 fprintf(f, "xlabel('u');\n"); 02729 fprintf(f, "ylabel('v');\n"); 02730 fprintf(f, "title('" 02731 "green = post-corner-splitting\\newline" 02732 "red=post-first-curve-quad-splitting\\newline" 02733 "cyan = post-outer-curve-triangle-splitting, first iteration\\newline" 02734 "magenta = post-inner-curve-triangle-splitting, first iteration\\newline" 02735 "yellow = post-inner-curve-triangle-splitting\\newline" 02736 "blue = post-outer-curve-triangle-splitting');\n"); 02737 // fprintf(f, "set(gca, 'ydir', 'reverse');\n"); // 100210: To get the same view as default 'goview'. 02738 fclose(f); 02739 02740 #endif 02741 02742 02743 02744 02745 // 02746 // 081210: Note that it does not make as much sense as before to return these now, if we have more than one 02747 // trimming curve. But the caller is probably not using them... As it is, we return only for the 02748 // outer trimming curve. 02749 // 02750 trim_curve_0 = trim_curve_all[0]; 02751 trim_curve_p_0 = trim_curve_p_all[0]; 02752 02753 return; 02754 } 02755 02756 02757 02758 02759 02760 02761 } // namespace Go
Generated on Tue Sep 21 15:44:24 2010 for GoTools Core by  doxygen 1.6.3