//=========================================================================== // 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. //=========================================================================== 00282 #ifndef DBG 00283 const bool dbg=false; // This way, all tests on 'dbg' is removed compile-time when DBG is not defined. 00284 #endif 00285 00286 // 00287 // At last we enter the old loop for counting intersections. 00288 // 00289 // 090117: Note that the rest of the code is very similar to 'segment_contour_intersection_for_s2m'. In 00290 // fact, we could have called that one. (Maybe we should have?) The difference is that now we know 00291 // that the segment is horizontal, (the ray,) and hence the code here can be made (slightly) simpler 00292 // and faster. One disadvantage is that the two functions (this and 00293 // 'segment_contour_intersection_for_s2m') must agree on how to handle degenerate cases and so 00294 // on. Otherwise strange problems will occur. 00295 00296 int crossings=0; 00297 00298 // const double eps=1e-8; // 100210: Increasing from 1e-12 to 1e-10. 00299 00300 const double abs_eps = 1e-8; // 100223: Trying to convert to these 00301 // const double snap_eps = 1e-5; 00302 00303 const int n=contour.size(); 00304 00305 if (dbg) 00306 printf("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ n=%d\n", n); 00307 00308 for (int i=0; i<n; i++) 00309 { 00310 const int j=(i+1)%n, pre_i=(i+n-1)%n; 00311 const double &preb=vertices[contour[pre_i]+1]; 00312 const double &prea=vertices[contour[pre_i]]; 00313 const double &a=vertices[contour[i]], &b=vertices[contour[i]+1]; 00314 const double &c=vertices[contour[j]], &d=vertices[contour[j]+1]; 00315 00316 const double len_squared = (a-c)*(a-c) + (b-d)*(b-d); 00317 if ( (len_squared<abs_eps*abs_eps) && (dbg) ) 00318 printf("\n\n HUH?! Segment with length %e, duplicate corner!\n\n\n", sqrt(len_squared)); 00319 00320 // if (dbg) 00321 // printf(" i=%3d: a=%g b=%g c=%g d=%g y0=%g abs_eps=%g !horiz=%d y0-inrange=%d(%d)\n", 00322 // i, a, b, c, d, y0, abs_eps, 00323 // (fabs(b-d)>abs_eps), 00324 // ((b-abs_eps<=y0) && (y0<d-abs_eps)), 00325 // ((d+abs_eps<y0) && (y0<b+abs_eps))); 00326 // if ( dbg && ((i==77) || (i==79)) ) 00327 // { 00328 // printf("\n\n\nhold on; plot(%f, %f, 'md', 'markersize', 10, 'linewidth', 3); hold off;\n", x0, y0); 00329 // printf("line([%f; %f], [%f; %f], 'color', 'm', 'linewidth', 5);\n\n\n", a, c, b, d); 00330 // } 00331 00332 // 00333 // 100223: Remember, crossings are counted when occuring at the *start* of segments, not the ends, 00334 // *except* when the curve has a vertical turn in this "start", i.e., there is a corner 00335 // "pointing" upward or downward. 00336 // 00337 00338 const bool not_horizontal_segment = fabs(b-d) > abs_eps; 00339 const bool y0_in_range_1 = (b-y0<=abs_eps) && (abs_eps<d-y0); // y0 in [b-eps, d-eps) 00340 const bool y0_in_range_2 = (abs_eps<y0-d) && (y0-b<=abs_eps); // y0 in (d+eps, b+eps] 00341 const bool y0_in_range = y0_in_range_1 || y0_in_range_2; 00342 00343 if ( not_horizontal_segment && y0_in_range ) 00344 { 00345 const double t=(y0-b)/(d-b); // Where on the segment does the ray cross? t=0 for b, and t=1 for d. 00346 00347 if (dbg) 00348 printf(" i=%3d: a=%g b=%g c=%g d=%g y0=%g abs_eps=%g !horiz=%d y0-inrange=%d %d\n", 00349 i, a, b, c, d, y0, abs_eps, 00350 not_horizontal_segment, y0_in_range_1, y0_in_range_2); 00351 00352 const bool interior = ((t>=abs_eps) && (t<=(1.0-abs_eps))); 00353 const bool start = (fabs(t)<=abs_eps); 00354 const bool going_up = ((preb<b) && (b<d)); 00355 const bool going_down = ((preb>b) && (b>d)); 00356 const bool not_turning_vertically = going_up || going_down; 00357 const bool actually_crossing = (interior || (start && not_turning_vertically)); 00358 const double intersection_x = a+t*(c-a); 00359 if (dbg) 00360 printf("\n\n t=%g, interior: %d, start: %d, going_up: %d, going_down: %d, " 00361 "actually_crossing: %d, inters_x: %g", 00362 t, interior, start, going_up, going_down, actually_crossing, intersection_x); 00363 if (dbg) 00364 printf("\n preb=%g, b=%g, d=%g x-values: %g %g %g", preb, b, d, prea, a, c); 00365 if ( actually_crossing && (intersection_x > x0+abs_eps) ) 00366 { 00367 crossings++; 00368 if (dbg) printf(", CROSSES "); 00369 if (dbg) 00370 { 00371 printf("\n\n\nhold on; plot(%f, %f, 'md', 'markersize', 10, 'linewidth', 3); hold off;\n", x0, y0); 00372 printf("line([%f; %f], [%f; %f], 'color', 'm', 'linewidth', 5);\n\n\n", a, c, b, d); 00373 } 00374 } 00375 else 00376 if (dbg) 00377 { 00378 if (actually_crossing) 00379 printf("\n intersection_x too small: %e vs. x0=%e, inters_x-x0=%e", 00380 intersection_x, x0, intersection_x-x0); 00381 printf("\n\n\nhold on; plot(%f, %f, 'bd', 'markersize', 10, 'linewidth', 3); hold off;\n", x0, y0); 00382 printf("line([%f; %f], [%f; %f], 'color', 'b', 'linewidth', 5);\n\n\n", a, c, b, d); 00383 } 00384 if (dbg) printf("\n"); 00385 } 00386 00387 } 00388 if (dbg) printf(" Crossings in total: %d, i.e., %s\n", crossings, ((crossings & 1)==1) ? "INSIDE" : "OUTSIDE"); 00389 00390 return ((crossings & 1)==1); 00391 } 00392 00393 00394 00395 00396 00397 00398 //============================================================================================================== 00399 // 00400 // 090115: See also comment for 'point_inside_contour' above. 00401 // 090117: This must be (re)checked before being used... 00402 // 00403 //============================================================================================================== 00404 00405 bool point_on_contour_corner(const double x0, const double y0, 00406 const double * const vertices, const vector<int> &contour) 00407 { 00408 const double eps=1e-13; // 090115: Used for zero-tests for distances in the parameter domain. 00409 // Hmm... these should *really*, *really* be taken from some global variable 00410 // or something 00411 for (int ilim=contour.size(), i=0; i<ilim; i++) 00412 { 00413 const double &corner_u=vertices[contour[i]], &corner_v=vertices[contour[i]+1]; 00414 const double dist_squared = (corner_u-x0)*(corner_u-x0) + (corner_v-y0)*(corner_v-y0); 00415 00416 if (dist_squared<eps*eps) 00417 return true; 00418 } 00419 return false; 00420 } 00421 00422 00423 00424 00425 00426 00427 //============================================================================================================== 00428 // 00429 // 090115: See also comment for 'point_inside_contour' above. 00430 // For a first test, just doing brute force search. Note that this will negate the fast algorithm used 00431 // for the inside-tests. The current function should do something similar, if this turns out to be a 00432 // good idea. 00433 // 090129: This can probably be done in point_inside_contour at very little extra cost... 00434 // 100213: Seems to be a problem with this routine, for some special cases... investigating... 00435 // 00436 //============================================================================================================== 00437 00438 inline bool point_on_contour(const double x0, const double y0, 00439 const double * const vertices, const vector<int> &contour 00440 #ifdef DBG 00441 , const bool dbg=false 00442 #endif 00443 ) 00444 { 00445 #ifndef DBG 00446 const bool dbg=false; // This way, all tests on 'dbg' is removed compile-time when DBG is not defined. 00447 #endif 00448 00449 //const double eps=1e-12; // 090115: Used for zero-tests for degeneracy and parallellity tests. 00450 // 090203: Here it is also used for testing if a point is on a curve segment. 00451 //const double eps=1e-8; // 090203: Need 1e-8 for proper meshing of 'bin_p1_3.g2'. 00452 //const double eps=1e-6; // 100213: See comments below. 00453 const double eps=1e-5; // 100214: Reverting to 1e-8, think the need to have 1e-6 is really another problem 00454 // 100218: Keeping 1e-5 after discussion with Vibeke. 00455 00456 const double tau=1e-12; // 090115: Used for zero-tests for distances in the parameter domain. 00457 // 100214: Changed from 1e-14 to 1e-12, not needed right now, but 1e-14 seems small. 00458 00459 const int n=contour.size(); 00460 if (dbg) printf(" point_on_contour for (%f, %f), %d segments.\n", x0, y0, n); 00461 for (int i=0; i<n; i++) 00462 { 00463 const int j=(i+1)%n; 00464 const double &x2=vertices[contour[i]], &y2=vertices[contour[i]+1]; 00465 const double &x3=vertices[contour[j]], &y3=vertices[contour[j]+1]; 00466 const double contour_segment_length_squared = (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2); 00467 00468 if (contour_segment_length_squared>tau*tau) 00469 { 00470 // We project (x0, y0) onto the segment, and measure the distance. 00471 // The projection will be (x2, y2) + s*[(x3, y3)-(x2, y2)]. 00472 const double s = ( (x0-x2)*(x3-x2) + (y0-y2)*(y3-y2) ) / contour_segment_length_squared; 00473 00474 // 100213: Hmm... Checking for s in (-tau, 1+tau) is no good, if 1+tau==1. Better to use (-eps, 00475 // 1+eps) when eps is larger than tau, but this does not solve the problem! Which is that 00476 // 's' is not a measure of distance. Should rather use 00477 // s*sqrt(contour_segment_length_squared)! 00478 // 00479 // 100214: But, then, remember that the interval should be [0, 00480 // sqrt(contour_segment_length_squared)]. 00481 // 00482 // 100213: Hmm again. Seems that 1e-8 (and 1e-7!) for some reason is too restrictive. (How can this 00483 // be?!) Changing to 1e-6. 00484 // 00485 // 100214: Suspect that the reason for 1e-8 being to large is that something else is really wrong 00486 // with the contour of the test case. (Folding back onto itself or something strange?!) 00487 00488 const double s2 = s * sqrt(contour_segment_length_squared); 00489 00490 if ( (dbg) && (0) ) 00491 { 00492 printf(" %3d: (x2, y2) = (%f, %f), (x3, y3) = (%f, %f)\n", i, x2, y2, x3, y3); 00493 printf("line([%f; %f], [%f; %f], 'color', 'y', 'linewidth', 5);\n", x2, x3, y2, y3); 00494 printf(" s=%f, s2=%f\n", s, s2); 00495 } 00496 00497 // if ((s>-tau) && (s<1.0+tau)) 00498 // if ((s>=-eps) && (s<=1.0+eps)) 00499 // if ((s2>=-eps) && (s2<=1.0+eps)) 00500 if ( (s2>=-eps) && (s2<=sqrt(contour_segment_length_squared)+eps) ) 00501 { 00502 const double x = x2 + s*(x3-x2), y = y2 + s*(y3-y2); 00503 const double dist_squared = (x-x0)*(x-x0) + (y-y0)*(y-y0); 00504 if (dbg) printf("\t\t\t\t\t\t\t\t\t\tdist to segment = %e\n", sqrt(dist_squared)); 00505 00506 // 100214: That the tolerance required in the next test is as large as 1e-4 really shows that 00507 // the trimming curve is crappy, so by using 1-e4 we are really just hiding a symptom... 00508 00509 if (dist_squared<eps*eps) 00510 //if (dist_squared<1e-6*1e-6) 00511 //if (dist_squared<1e-4*1e-4) 00512 { 00513 if (dbg) printf("TRUE\n"); 00514 return true; 00515 } 00516 } 00517 } 00518 else 00519 { 00520 // 00521 // 090117: The contour segment is degenerate, and will be treated as a point. 00522 // 00523 // 100214: We should not get here, duplicate points should have been removed, but it's no harm in 00524 // keeping the branch I guess... 00525 // 00526 const double mx = 0.5*(x2+x3), my = 0.5*(y2+y3); 00527 const double dist_squared = (x0-mx)*(x0-mx) + (y0-my)*(y0-my); 00528 if (dist_squared<eps*eps) 00529 { 00530 if (dbg) printf("TRUE\n"); 00531 return true; 00532 } 00533 } 00534 } 00535 00536 if (dbg) printf("FALSE\n"); 00537 return false; 00538 } 00539 00540 00541 00542 00543 00544 00545 //============================================================================================================== 00546 // 00547 // 100219: See comments with corresponding dates in 'segment_contour_intersection_for_s2m()'. 00548 // 00549 //============================================================================================================== 00550 00551 bool point_on_segment(const double x0, const double y0, 00552 const double x2, const double y2, 00553 const double x3, const double y3, 00554 const double zero_eps, 00555 const double snap_eps 00556 #ifdef DBG 00557 , const bool dbg = false 00558 #endif 00559 ) 00560 { 00561 #ifndef DBG 00562 const bool dbg = false; // This way, all tests on 'dbg' is removed compile-time when DBG is not defined. 00563 #endif 00564 00565 const double t = ( (x0-x2)*(x3-x2) + (y0-y2)*(y3-y2) ) / ( (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2) ); 00566 const double x = x2 + t*(x3-x2) - x0; 00567 const double y = y2 + t*(y3-y2) - y0; 00568 const double dist = sqrt( x*x+y*y ); 00569 if ( (dist<snap_eps) && (t>=-zero_eps) && (t-zero_eps<=1.0) ) 00570 { 00571 if (dbg) 00572 { 00573 printf(" END TESTS 2: t=%f, dist=%f\n", t, dist); 00574 // printf("\nline([%f; %f], [%f; %f], 'color', 'k', 'linewidth', 6);\n", x0, x1, y0, y1); 00575 printf("line([%f; %f], [%f; %f], 'color', 'y', 'linewidth', 6);\n", x2, x3, y2, y3); 00576 } 00577 return true; 00578 } 00579 00580 return false; 00581 } 00582 00583 00584 00585 00586 00587 00588 //============================================================================================================== 00589 // 00590 // 081208: Comment added today. 00591 // 00592 // The contour consists of n segments, so we are finding an intersection between two straight line 00593 // pieces. The corners of the contour polygon are stored in 'vertices'. 00594 // The 'x0', 'y0', 'x1' and 'y1' seem to be points in the parameter domain of the surface/curve. 00595 // 00596 // 081209: There seems to be a problem with segments of length 0 in the contour... why haven't this happened 00597 // before? Or is it just a symptom of something else being the matter? 00598 // 00599 // 090115: The return value is true if an intersection between the segment and contour is found. 00600 // 00601 // 090117: Note that the contour may (of course) intersect the segment more than once. Just one intersection 00602 // is returned, and it is pretty indetermined which one it will be. 00603 // 00604 // 090117: If we modify the behaviour of "insideness-tests" in spline2mesh, this must be mirrored 00605 // here. Otherwise, 'split_quad' will not find the intersections it expect, and fail! 00606 // Or - we could do it in the caller, e.g., in 'split_quad'. 00607 // 00608 // 090203: There is a fundamental problem with the whole approach. In some cases, a quad (or it might as well 00609 // be any polygon, e.g., a triangle) has an edge that will intersect the trimming curve ('contour') in 00610 // more than one place. For example, a contour which follows an iso-parameter line, say a border of 00611 // the surface, and then turns sharpely back, making out a wedge. Somewhere near the tip of this edge, 00612 // this situation will always occur, independent of tesselation level etc. If the segment of a quad 00613 // (triangle) normal to the border is determined to intersect the contour in the end, i.e., at the 00614 // border of the surface, the wedge tip will be cut away wrongly... 00615 // How can this be solved?? 00616 // Not sure if this is a universal fix, but for now, we will favor inner intersections before 00617 // "end-of-segment" intersections... 00618 // This also means that we will not be able to return as soon as we find an intersection, slowing 00619 // things down somewhat. Considerably? 00620 // A compromise: When first intersection found is of type "end", continue until a "better" is found. 00621 // 00622 // 100218: This situation has again surfaced. Why does it not work, if it was fixed before? 00623 // Making the FIX090203 permanent. 00624 // Increasing tolerances. 00625 // 00626 // 100219: Hmm... I do not entirely like the way things have turned here. By favoring inner intersections, 00627 // we are in a way handling the case where there are many intersections on an edge, and choosing one 00628 // of them, as opposed to finding *all* and also treating *all*... 00629 // 00630 // 100219: There was a "bug" here: Originally the routine was not supposed to classify any intersection on 00631 // parallel segments as an intersection. However, this is needed when called from 'split_triangle' 00632 // (as opposed to from a ray-contour problem of a point-inside-contour context,) where we want 00633 // intersections to be found when the ends of the triangle edge lie *on* (or, actually, 1e-5 or so 00634 // from) such a parallel (to the edge) contour segment! Fixing this now, and adding a flag, so that 00635 // it will not be the default behaviour. ('snap_ends', see also comments in the code below.) 00636 // 00637 // 100225: Should 'snap_ends' be turned to default=true, and then removed altogether when proving 00638 // itself to be a good idea? 00639 // 00640 // 100222: The question is, then, who (and when) should turn this on? Following the conservative approach, 00641 // which is to make it default to false, and turn it (minimally) on in order to solve the current 00642 // (trim3.g2) case... 00643 // 00644 // 100222: One thing, which is not entirely satisfactory with the whole implementation: 00645 // 00646 // The "insidedness/onness" for points vs. contours used in routines for determination of which, 00647 // and how, triangles are to be clipped, is not the same, nor equivalent, to those used during the 00648 // actual clipping. The reason for this is, e.g., that these two operations due to their 00649 // implementation/organization require slightly different definitions of insidedness. This could 00650 // (if true, then it should) maybe be fixed. 00651 // 00652 //============================================================================================================== 00653 00654 bool segment_contour_intersection_for_s2m(const double x0, const double y0, 00655 const double x1, const double y1, 00656 const double * const vertices, 00657 const vector<int> &contour, 00658 double &x, double &y, // Resulting intersection 00659 double &s, // param for the segment containing the inters, [0, 1]. 00660 const bool snap_ends /* = false */ 00661 #ifdef DBG 00662 , const bool dbg /* = false */ 00663 #endif 00664 ) 00665 { 00666 #ifndef DBG 00667 const bool dbg=false; // This way, all tests on 'dbg' is removed compile-time when DBG is not defined. 00668 #endif 00669 00670 if (dbg) puts("\n\n # Entering segment_contour_intersection_for_s2m ########################################\n"); 00671 00672 // 00673 // Now we search for the (actually, just "a") piece of the contour which intersects the segment. Maybe it is 00674 // faster to just do this without the previous test, if we can use the 'sorted_segments' information. 00675 // 00676 00677 int i; 00678 //const double tau=1e-14; // 090115: Used for zero-tests for degeneracy and parallellity tests. 00679 const double tau=1e-12; // 100218: 00680 00681 const double parallellity_eps = 1e-10; // 100219 00682 00683 const double eps=1e-13; // 090115: Used for zero-tests for distances in the parameter domain. 00684 //const double eps2=1e-8; // 090203: This is used to determine whether or not an intersection is at one of the ends 00685 const double eps2=1e-5; // 100218: 00686 // of an edge. (In which case we continue to look for an "interior" intersection, 00687 // which will be favoured. Is this always sensible? Or is it possible to make a 00688 // case in which the opposite behaviour should be preferred? 00689 // 090204: We need this as large as 1e-8 for bin_p1_3.g2 to be meshed properly. 00690 bool intersection_found_at_the_end_of_the_segment = false; 00691 double x_saved=1e99, y_saved=1e99, s_saved=1e99; // Initial values to help spot unforeseen problems... 00692 const int n=contour.size(); 00693 for (i=0; i<n; i++) 00694 { 00695 const int ii=i, jj=(ii+1)%n; 00696 const double &x2=vertices[contour[ii]], &y2=vertices[contour[ii]+1]; 00697 const double &x3=vertices[contour[jj]], &y3=vertices[contour[jj]+1]; 00698 double t; 00699 00700 #if 0 00701 if (dbg) 00702 { 00703 const double d02 = sqrt( (x2-x0)*(x2-x0) + (y2-y0)*(y2-y0) ); 00704 const double d03 = sqrt( (x3-x0)*(x3-x0) + (y3-y0)*(y3-y0) ); 00705 const double d12 = sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) ); 00706 const double d13 = sqrt( (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1) ); 00707 00708 const double eps3 = 1e-3; 00709 00710 if ( (d02<eps3) || (d03<eps3) || (d12<eps3) || (d13<eps3) || (i==111169) ) 00711 { 00712 printf(" i=%4d: END TESTS: d02=%f, d03=%f, d12=%f, d13=%f\n", i, d02, d03, d12, d13); 00713 printf("\nline([%f; %f], [%f; %f], 'color', 'k', 'linewidth', 6);\n", x0, x1, y0, y1); 00714 printf("line([%f; %f], [%f; %f], 'color', 'y', 'linewidth', 6);\n", x2, x3, y2, y3); 00715 } 00716 00717 // testing if ends of first segment is in the other. 00718 00719 const double t0 = ( (x0-x2)*(x3-x2) + (y0-y2)*(y3-y2) ) / ( (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2) ); 00720 double p = x2+t0*(x3-x2) - x0; 00721 double q = y2+t0*(y3-y2) - y0; 00722 const double dist0 = sqrt(p*p+q*q); 00723 //if ( (dist0<1e-5) && (t0>=-eps) && (t0<=1.0+eps) ) 00724 if (i==69) 00725 { 00726 printf(" i=%4d: END TESTS 2: t0=%f, dist=%f\n", i, t0, dist0); 00727 printf("\nline([%f; %f], [%f; %f], 'color', 'k', 'linewidth', 6);\n", x0, x1, y0, y1); 00728 printf("line([%f; %f], [%f; %f], 'color', 'y', 'linewidth', 6);\n", x2, x3, y2, y3); 00729 } 00730 00731 p = x2+t0*(x3-x2) - x1; 00732 q = y2+t0*(y3-y2) - y1; 00733 const double dist1 = sqrt(p*p+q*q); 00734 if (i==69) 00735 //if ( (dist1<1e-5) && (t0>=-eps) && (t0<=1.0+eps) ) 00736 { 00737 printf(" i=%4d: END TESTS 3: t0=%f, dist=%f\n", i, t0, dist0); 00738 printf("\nline([%f; %f], [%f; %f], 'color', 'k', 'linewidth', 6);\n", x0, x1, y0, y1); 00739 printf("line([%f; %f], [%f; %f], 'color', 'y', 'linewidth', 6);\n", x2, x3, y2, y3); 00740 } 00741 00742 } 00743 #endif 00744 00745 // 081209: Adding this test... 00746 const double contour_segment_length_squared = (x2-x3)*(x2-x3)+(y2-y3)*(y2-y3); 00747 // if (dbg2) printf("i=%3d cont segm len=%f ", i, sqrt(contour_segment_length_squared)); 00748 00749 if (contour_segment_length_squared>tau*tau) 00750 { 00751 00752 if (snap_ends) 00753 { 00754 // 00755 // 100219: We check if the ends of the segment (the one given as input to the function) are 00756 // close to the contour piece in question. If 'snap_ends' is true, we could have that 00757 // situation even if the pieces are (close to) parallel. We do this testing now, 00758 // before the parallellity-test, since that test is too strict to determine if this 00759 // 'snap_ends'-test should be performed. (It is too strict because it is meant to trap 00760 // division with zero before they happen, mainly.) 00761 // 00762 // Note also that if one of the ends is on the contour, this is an "end intersection", 00763 // and this 'intersection_found_at_the_end_of_the_segment' is set and the search 00764 // continued, in case there is also some inner intersection. 00765 // 00766 00767 if (point_on_segment(x0, y0, x2, y2, x3, y3, eps, eps2)) // , prøve med true her og nedenfor, se om det alltid blir bra da.)) 00768 { 00769 intersection_found_at_the_end_of_the_segment = true; 00770 x_saved = x0; 00771 y_saved = y0; 00772 s_saved = 0.0; 00773 } 00774 else 00775 if (point_on_segment(x1, y1, x2, y2, x3, y3, eps, eps2 DBG_FLAG)) 00776 { 00777 intersection_found_at_the_end_of_the_segment = true; 00778 x_saved = x1; 00779 y_saved = y1; 00780 s_saved = 1.0; 00781 } 00782 } 00783 00784 // 090115: The next one should be zero for contour segment and segment parallel. If so, no intersection... 00785 const double parallellity = (x3-x2)*(y1-y0)-(y3-y2)*(x1-x0); 00786 if ( fabs(parallellity) > parallellity_eps ) // i.e., not parallel... 00787 { 00788 // 00789 // 100219: Let the (triangle) segment has ends A=(x0, y0) and B=(x1, y1), and the contour 00790 // piece has ends C=(x2, y2) and D=(x3, y3). With intersection S = A+s(B-A) = 00791 // C+t(D-C), we get the following: 00792 // 00793 t = ((x0-x2)*(y1-y0)-(y0-y2)*(x1-x0)) / parallellity; 00794 if ((t>=-eps) && (t<=1.0+eps)) 00795 { 00796 if (fabs(x1-x0)>fabs(y1-y0)) 00797 s = ((x2-x0)+t*(x3-x2))/(x1-x0); 00798 else 00799 s = ((y2-y0)+t*(y3-y2))/(y1-y0); 00800 // if (dbg2) printf("s=%f ", s); 00801 if ((s>=-eps) && (s<=1.0+eps)) 00802 { 00803 if (dbg) 00804 { 00805 const double ss = ((x0-x2)*(y3-y2)-(y0-y2)*(x3-x2)) / parallellity; 00806 printf(" i=%4d: s=%f (%f), t=%f, parallellity=%f\n", i, s, ss, t, parallellity); 00807 printf("\nline([%f; %f], [%f; %f], 'color', 'k', 'linewidth', 6);\n", x0, x1, y0, y1); 00808 printf("line([%f; %f], [%f; %f], 'color', 'y', 'linewidth', 6);\n", x2, x3, y2, y3); 00809 const double px = x0 + s*(x1-x0), py = y0 + s*(y1-y0); 00810 const double qx = x2 + t*(x3-x2), qy = y2 + t*(y3-y2); 00811 printf("hold on; plot(%f, %f, 'gs', 'markersize', 15, 'linewidth', 2);\n", px, py); 00812 printf("plot(%f, %f, 'ks', 'markersize', 15, 'linewidth', 2); hold off;\n\n", qx, qy); 00813 } 00814 s = std::max(std::min(s, 1.0), 0.0); 00815 x = x2+t*(x3-x2); 00816 y = y2+t*(y3-y2); 00817 //printf("***** s=%f t=%f x=%f y=%f\n", s, t, x, y); 00818 //printf("returning true from %s:%d...\n", __FILE__, __LINE__); 00819 if ((fabs(s)<eps2) || (fabs(s-1.0)<eps2)) 00820 { 00821 intersection_found_at_the_end_of_the_segment=true; 00822 if (dbg) printf("Intersection is at and end.\n"); 00823 x_saved=x; 00824 y_saved=y; 00825 s_saved=s; 00826 } 00827 else 00828 { 00829 if (dbg) puts(" # Returning true from segment_contour_intersection_for_s2m # (a) #####\n"); 00830 return true; 00831 } 00832 } 00833 else 00834 if (dbg) printf(" i=%4d: NO INTERS, s OUT OF RANGE s=%f, t=%f, p=%g\n", i, s, t, parallellity); 00835 } 00836 else 00837 if (dbg) printf(" i=%4d: no inters, t out of range s=%f, t=%f, p=%g\n", i, s, t, parallellity); 00838 } 00839 else 00840 if (dbg) 00841 { 00842 printf(" i=%4d: PARALLEL parallellity=%e\n", i, parallellity); 00843 // The egde of the triangle ("segment") and the piece of the contour: 00844 printf("\nline([%f; %f], [%f; %f], 'color', 'k', 'linewidth', 6);\n", x0, x1, y0, y1); 00845 printf("line([%f; %f], [%f; %f], 'color', 'y', 'linewidth', 6);\n\n", x2, x3, y2, y3); 00846 } 00847 } 00848 else 00849 { 00850 if (dbg) printf(" Hmm... degen contour segment... should that still be possible?!\n"); 00851 // 090117: The contour segment is degenerate, and will be treated as a point. 00852 const double mx = 0.5*(x2+x3), my = 0.5*(y2+y3); 00853 00854 // The projection will be (x0, y0) + s*[(x1, y1)-(x0, y0)]. 00855 s = ( (mx-x0)*(x1-x0)+(my-y0)*(y1-y0) ) / ( (x1-x0)*(x1-x0)+(y1-y0)*(y1-y0) ); 00856 // if (dbg2) printf("s=%f ", s); 00857 00858 if ((s>-eps) && (s<1.0+eps)) 00859 { 00860 x = x0 + s*(x1-x0); 00861 y = y0 + s*(y1-y0); 00862 00863 const double dist_squared = (x-mx)*(x-mx) + (y-my)*(y-my); 00864 if (dist_squared<eps*eps) 00865 { 00866 intersection_found_at_the_end_of_the_segment=true; 00867 x_saved=x; 00868 y_saved=y; 00869 s_saved=s; 00870 } 00871 } 00872 } 00873 // if (dbg2) printf("\n"); 00874 } 00875 00876 if (intersection_found_at_the_end_of_the_segment) 00877 { 00878 x=x_saved; 00879 y=y_saved; 00880 s=s_saved; 00881 if (dbg) puts(" # Returning true from segment_contour_intersection_for_s2m # (b) #######################\n"); 00882 return true; 00883 } 00884 00885 // 090205: This message is not longer as relevant (or correct!) Now (also before?) it is perfectly 00886 // acceptable that this returns 'false' for the "second intersection, first try". The problem is 00887 // that it should not happen in the "second try" done by 'split_triangle'... But the proper place to 00888 // catch that problem is most likely in 'split_triangle' itself, so I'll disable this old message... 00889 // MESSAGE("\nCouldn't find proper intersection between segment\n" 00890 // "and contour, while the segment's ends\n" 00891 // "are determined to be on opposite sides of the contour!!!\n" 00892 // "If you see this, the results are very likely going to look like crap, unfortunately.\n" 00893 // "It is currently unknown how this could happen.\n"); 00894 // printf("Segment length in the parameter domain: %g\n\n\n", sqrt((x0-x1)*(x0-x1) + (y0-y1)*(y0-y1))); 00895 00896 if (dbg) puts(" # Returning false from segment_contour_intersection_for_s2m # (c) ##########################\n"); 00897 00898 return false; 00899 } 00900 00901 00902 00903 00904 00905 00906 //============================================================================================================== 00907 // 00908 //============================================================================================================== 00909 00910 vector< short_list_short_list > sort_2dpoly_segments(const double * const vertices, 00911 const vector<int> &contour, 00912 const bool transposed /* =false */) 00913 { 00914 vector< short_list_short_list > sorted_segments; 00915 const int n=contour.size(); 00916 00917 const int debug=0; 00918 00919 // 020615: Introducing the 'transposed' for making a version in which x and y swap roles. 00920 00921 // 00922 // We cannot choose, it must contain indices into the 'contour' vector, if we were to copy that vector's 00923 // indices directly into 'vertices', we would lose information about one of the ends of all segments, as 00924 // soon as we started reordering the segments. Hence, initializing with 'contour[i]' was wrong. It must be 00925 // 'i'. 00926 sorted_segments.resize(n); 00927 for (int i=0; i<n; i++) 00928 sorted_segments[i].first=i; 00929 00930 // Now we sort this "outer" list wrt. largest x-value of the segments. 00931 { 00932 less_largest_x compare_functor(contour, vertices, 0.0, transposed); 00933 std::sort(sorted_segments.begin(), sorted_segments.end(), compare_functor); 00934 } 00935 00936 if ((debug) && (!transposed)) 00937 { 00938 //short_list_short_list *first_cut; 00939 vector<short_list_short_list>::iterator first_cut; 00940 00941 printf("\n"); 00942 for (first_cut=sorted_segments.begin(); first_cut<sorted_segments.end(); first_cut++) 00943 { 00944 const int a=contour[first_cut->first], b=contour[(first_cut->first+1)%n]; 00945 printf("%3d: %5d-%5d: (%10.3f, %10.3f, %10.3f)-" 00946 "(%10.3f, %10.3f, %10.3f) maxx: %10.3f\n", 00947 int(first_cut-sorted_segments.begin()), 00948 a, b, 00949 vertices[a+0], vertices[a+1], vertices[a+2], 00950 vertices[b+0], vertices[b+1], vertices[b+2], 00951 std::max(vertices[a+0], vertices[b+0])); 00952 } 00953 } 00954 00955 // Next, we fill inn all the lists of each pair of this "outer" list, i.e., the "middle" lists. 00956 if (transposed) 00957 { 00958 THROW("Huh?! This should not be in use, I think... J.O. 090129"); 00959 exit(0); 00960 } 00961 less_smallest_y lsy_comp(contour, vertices, 0.0); // , transposed); 00962 less_largest_y lly_comp(contour, vertices, 0.0); // , transposed); 00963 00964 for (int i=0; i<n; i++) 00965 { 00966 sorted_segments[i].second=new vector<short_list>; 00967 sorted_segments[i].second->resize(n-i); 00968 for (int j=0; j<n-i; j++) 00969 (*sorted_segments[i].second)[j].first=sorted_segments[i+j].first; 00970 00971 // Now we sort this "middle" list wrt. largest y-value of the segments. 00972 std::sort(sorted_segments[i].second->begin(), sorted_segments[i].second->end(), lly_comp); 00973 00974 if ((debug) && (!transposed)) 00975 { 00976 vector<short_list>::iterator second_cut; 00977 00978 printf("\n"); 00979 for (second_cut=sorted_segments[i].second->begin(); second_cut<sorted_segments[i].second->end(); second_cut++) 00980 { 00981 const int a=contour[second_cut->first], b=contour[(second_cut->first+1)%n]; 00982 printf(" %3d: %5d-%5d: (%10.3f, %10.3f, %10.3f)-" 00983 "(%10.3f, %10.3f, %10.3f) maxy: %10.3f\n", 00984 int(second_cut-sorted_segments[i].second->begin()), 00985 a, b, 00986 vertices[a+0], vertices[a+1], vertices[a+2], 00987 vertices[b+0], vertices[b+1], vertices[b+2], 00988 std::max(vertices[a+1], vertices[b+1])); 00989 } 00990 } 00991 00992 // Next, we fill inn all the lists of each pair of this "middle" list, i.e., the "inner" lists. 00993 for (int j=0; j<n-i; j++) 00994 { 00995 (*sorted_segments[i].second)[j].second=new vector<short>; 00996 (*sorted_segments[i].second)[j].second->resize(n-i-j); 00997 for (int k=0; k<n-i-j; k++) 00998 (*(*sorted_segments[i].second)[j].second)[k] = (*sorted_segments[i].second)[j+k].first; 00999 01000 // Now we sort this "inner" list wrt. smallest y-value of the segments. Default 'compare' should be 01001 // 'less'. 01002 // 01003 // 020614: Don't know why such a default comparison functor is mentioned. How could it ever work 01004 // with these compounded lists with the "strange" comparisons that are to be done?! 01005 // 020615: Probably I thought that as the inner list's entries consist of just an integer, there is 01006 // no need for the complex compare_functor, but then I forgot that it's not the integers (or 01007 // shorts) that are to be used as key in the sorting, but that which they point to! 01008 std::sort((*sorted_segments[i].second)[j].second->begin(), 01009 (*sorted_segments[i].second)[j].second->end(), 01010 lsy_comp); 01011 01012 if ((debug) && (!transposed)) 01013 { 01014 vector<short>::iterator third_cut; 01015 01016 printf("\n"); 01017 for (third_cut=(*sorted_segments[i].second)[j].second->begin(); 01018 third_cut<(*sorted_segments[i].second)[j].second->end(); 01019 third_cut++) 01020 { 01021 const int a=contour[*third_cut], b=contour[(*third_cut+1)%n]; 01022 printf(" %3d: %5d-%5d: (%10.3f, %10.3f, %10.3f)-" 01023 "(%10.3f, %10.3f, %10.3f) miny: %10.3f\n", 01024 int(third_cut- 01025 (*sorted_segments[i].second)[j].second->begin()), 01026 a, b, 01027 vertices[a+0], vertices[a+1], vertices[a+2], 01028 vertices[b+0], vertices[b+1], vertices[b+2], 01029 std::min(vertices[a+1], vertices[b+1])); 01030 } 01031 } 01032 01033 } 01034 } 01035 01036 // if (transposed) 01037 // { 01038 // printf("\n"); 01039 // print_sorted_segments2(sorted_segments, contour, vertices); 01040 // } 01041 01042 return sorted_segments; 01043 } 01044 01045 01046 01047 01048 01049 01050 //============================================================================================================== 01051 // 01052 // Function for testing of insidedness of point based on old piecewise linear contours. 01053 // 01054 // Note: The z-coo. is ignored by that routine. The curve is specified in the parameter domain. 01055 // 01056 // 081202: So what is the difference between 'trim_curve_p' and 'contour'?! 01057 // Ok, 'contour' is a list of indices into 'trim_curve_p', it seems. 01058 // (Suspecting that for usage here, 'contour' is simply {0, ..., n} for the proper 'n'...) 01059 // (Actually {0, 3, ..., 3*n}.) 01060 // 01061 // 100211: Removing traces of "sorted_segments". Also, making the 'contour' a reference. Why was it not?!?! 01062 // 01063 //============================================================================================================== 01064 01065 int is_inside(const vector< Vector3D > &trim_curve_p, 01066 const vector<int> &contour, 01067 const double u, const double v 01068 #ifdef DBG 01069 , const bool dbg /* =false */ 01070 #endif 01071 ) 01072 { 01073 #ifndef DBG 01074 // const bool dbg=false; // This way, all tests on 'dbg' is removed compile-time when DBG is not defined. 01075 #endif 01076 01077 01078 // test 100212 01079 01080 if (contour.size()==0) 01081 return 0; 01082 01083 #if 0 // 021012 01084 01085 const double * const vertices = &trim_curve_p[0][0]; 01086 01087 01088 vector<int> cont; 01089 cont.push_back(contour[0]); 01090 for (int i=0; i<int(contour.size()); i++) 01091 { 01092 const int j = (i+1)%int(contour.size()); 01093 const double &a=vertices[contour[i]], &b=vertices[contour[i]+1]; 01094 const double &c=vertices[contour[j]], &d=vertices[contour[j]+1]; 01095 const double len_squared = (a-c)*(a-c) + (b-d)*(b-d); 01096 const double eps = 1e-8; 01097 if ( len_squared < eps*eps ) 01098 { 01099 if (dbg) 01100 printf("is_inside preprocessing: HUH?! Segment with length %e, duplicate corner! REMOVING IT\n", 01101 sqrt(len_squared)); 01102 } 01103 else 01104 cont.push_back(contour[i]); 01105 } 01106 01107 01108 for (int s=2; s<=10; s++) 01109 { 01110 vector<int> cont2; 01111 cont2.push_back(cont[0]); 01112 for (int i=0; i<int(cont.size()); i++) 01113 { 01114 const int j = (i+1)%int(cont.size()); 01115 const double &a=vertices[cont[i]], &b=vertices[cont[i]+1]; 01116 const double &c=vertices[cont[j]], &d=vertices[cont[j]+1]; 01117 const double len_squared = (a-c)*(a-c) + (b-d)*(b-d); 01118 const double eps = 1e-8; 01119 if ( len_squared < eps*eps ) 01120 { 01121 if (dbg) 01122 printf("is_inside prep STAGE %d : HUH?! Segment %d with length %e, duplicate corner! REMOVING IT\n", 01123 s, i, sqrt(len_squared)); 01124 } 01125 else 01126 cont2.push_back(cont[i]); 01127 } 01128 cont=cont2; 01129 } 01130 01131 return point_inside_contour(u, v, &trim_curve_p[0][0], cont, dbg); 01132 #endif 01133 01134 01135 01136 01137 // 100212: The const_cast not needed any more. 01138 //return point_inside_contour(u, v, &trim_curve_p[0][0], const_cast<vector<int> &>(contour), dbg); 01139 return point_inside_contour(u, v, &trim_curve_p[0][0], contour DBG_FLAG); 01140 } 01141 01142 01143 01144 01145 01146 01147 //============================================================================================================== 01148 // 01149 // 090115: 01150 // 090117: This should be (re)checked before being used... Not currently in use. 01151 // 01152 //============================================================================================================== 01153 01154 bool is_on_corner(const vector< Vector3D > &trim_curve_p, const vector<int> &contour, const double u, const double v) 01155 { 01156 return point_on_contour_corner(u, v, &trim_curve_p[0][0], const_cast<vector<int> &>(contour)); 01157 } 01158 01159 01160 01161 01162 01163 01164 //============================================================================================================== 01165 // 01166 // 090117: 01167 // 01168 //============================================================================================================== 01169 01170 int is_on_contour(const vector< Vector3D > &trim_curve_p, const vector<int> &contour, const double u, const double v 01171 #ifdef DBG 01172 , const bool dbg /* = false */ 01173 #endif 01174 ) 01175 { 01176 return point_on_contour(u, v, &trim_curve_p[0][0], const_cast<vector<int> &>(contour) DBG_FLAG); 01177 } 01178 01179 01180 01181 01182 01183 01184 //============================================================================================================== 01185 // 01186 // 090204: 01187 // 100218: Eliminating unneeded test. Increasing the tolerance. 01188 // 01189 //============================================================================================================== 01190 01191 bool degenerate_triangle(const Vector2D &c1, const Vector2D &c2, const Vector2D &c3 01192 #ifdef DBG 01193 , const bool dbg /* = false */ 01194 #endif 01195 ) 01196 { 01197 #ifndef DBG 01198 const bool dbg = false; 01199 #endif 01200 01201 // const double eps=1e-14; 01202 const double eps=1e-9; 01203 const double eps_squared=eps*eps; 01204 01205 Vector2D e=c2-c1, f=c3-c1; 01206 const double e_len_squared=e*e, f_len_squared=f*f, g_len_squared=(c2-c3)*(c2-c3); 01207 01208 if (dbg) printf(" degen_tri: edge lengths: %f %f %f\n", sqrt(e*e), sqrt(f*f), sqrt(g_len_squared)); 01209 01210 if ( (e_len_squared<eps_squared) || (f_len_squared<eps_squared) || (g_len_squared<eps_squared) ) 01211 return true; 01212 01213 // Ok, even if all lengths are non-zero, it may still be degenerate. Then the area will be zero. 01214 // 01215 // 100219: Hmm... come to think of it, this test should be the only one needed... But on the other hand, 01216 // if we use the division below, we must still test on the size of the denominator, so we might as 01217 // well keep it as it is. (Or maybe exchange the whole thing for an area-test?!) 01218 01219 const double ef=(e*f)/sqrt(e_len_squared*f_len_squared); 01220 01221 if (dbg) printf(" degen_tri: ef=%e, fabs(ef)-1.0=%e\n", ef, fabs(ef)-1.0); 01222 01223 if (fabs(fabs(ef)-1.0) < eps) 01224 return true; 01225 01226 return false; 01227 } 01228 01229 01230 01231 01232 01233 01234 } // namespace Go
Generated on Tue Sep 21 15:44:24 2010 for GoTools Core by  doxygen 1.6.3