00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "stdafx.h"
00021
00022 #include <iostream>
00023
00024 #include "gul_types.h"
00025 #include "gul_vector.h"
00026 #include "guar_exact.h"
00027 #include "guar_intersect.h"
00028 #include "gul_io.h"
00029
00030 namespace guar {
00031
00032 using gul::Ptr;
00033 using gul::point;
00034 using gul::point2;
00035 using gul::triangle;
00036 using gul::triangle2;
00037 using gul::cross_product;
00038 using gul::test;
00039 using gul::compare;
00040 using std::cout;
00041
00042
00043
00044
00045
00046 bool RegularIntersectLines(
00047 const point<rational>& A,
00048 const point<rational>& B,
00049 const point<rational>& a,
00050 const point<rational>& b,
00051 rational *lambda,
00052 rational *mu
00053 )
00054 {
00055 rational A0,A1,A2,B0,B1,B2;
00056 rational a0,a1,a2,b0,b1,b2;
00057 rational N1,N2,N3;
00058
00059
00060 A0 = A.x; A1 = A.y; A2 = A.z;
00061 B0 = B.x; B1 = B.y; B2 = B.z;
00062 a0 = a.x; a1 = a.y; a2 = a.z;
00063 b0 = b.x; b1 = b.y; b2 = b.z;
00064
00065 N1 = B1*b0 - b1*B0;
00066
00067 if( test(N1) )
00068 {
00069 *lambda = (b0*(a1-A1) - b1*(a0-A0)) / N1;
00070 *mu = (B0*(a1-A1) - B1*(a0-A0)) / N1;
00071 }
00072 else
00073 {
00074 N2 = B2*b0 - b2*B0;
00075
00076 if( test(N2) )
00077 {
00078 *lambda = (b0*(a2-A2) - b2*(a0-A0)) / N2;
00079 *mu = (B0*(a2-A2) - B2*(a0-A0)) / N2;
00080 }
00081 else
00082 {
00083 N3 = B2*b1 - b2*B1;
00084
00085 if( test(N3) )
00086 {
00087 *lambda = (b1*(a2-A2) - b2*(a1-A1)) / N3;
00088 *mu = (B1*(a2-A2) - B2*(a1-A1)) / N3;
00089 }
00090 else
00091 {
00092 return false;
00093 }
00094 }
00095 }
00096 return true;
00097 }
00098
00099
00100
00101
00102
00103 template< class T >
00104 GULAPI void BarycentricCoordinates(
00105 const point<T>& P1, const point<T>& P2,
00106 const point<T>& P3, const point<T>& Pin,
00107 T *w, T *u, T *v )
00108 {
00109 point2<T> W,U,V,P;
00110 point2<T> WU,WV,UV,WP,UP;
00111 point<T> UVW,PVW,PWU,PUV,P12,P13,N;
00112 T aUVW,aPVW,aPWU,aPUV;
00113 rational one_neg(ULong(1),-1);
00114
00115 P12 = P2 - P1;
00116 P13 = P3 - P1;
00117 N = cross_product( P12, P13 );
00118
00119
00120 if( test(N.z) )
00121 {
00122 W.x = P1.x;
00123 W.y = P1.y;
00124
00125 U.x = P2.x;
00126 U.y = P2.y;
00127
00128 V.x = P3.x;
00129 V.y = P3.y;
00130
00131 WU.x = P12.x;
00132 WU.y = P12.y;
00133
00134 WV.x = P13.x;
00135 WV.y = P13.y;
00136
00137 P.x = Pin.x;
00138 P.y = Pin.y;
00139 }
00140 else if( test(N.y) )
00141 {
00142 W.x = P1.x;
00143 W.y = P1.z;
00144
00145 U.x = P2.x;
00146 U.y = P2.z;
00147
00148 V.x = P3.x;
00149 V.y = P3.z;
00150
00151 WU.x = P12.x;
00152 WU.y = P12.z;
00153
00154 WV.x = P13.x;
00155 WV.y = P13.z;
00156
00157 P.x = Pin.x;
00158 P.y = Pin.z;
00159 }
00160 else
00161 {
00162 gul::Assert<InternalError>( ndebug || test(N.x) );
00163
00164 W.x = P1.y;
00165 W.y = P1.z;
00166
00167 U.x = P2.y;
00168 U.y = P2.z;
00169
00170 V.x = P3.y;
00171 V.y = P3.z;
00172
00173 WU.x = P12.y;
00174 WU.y = P12.z;
00175
00176 WV.x = P13.y;
00177 WV.y = P13.z;
00178
00179 P.x = Pin.y;
00180 P.y = Pin.z;
00181 }
00182
00183 UVW = cross_product(WU,WV);
00184 aUVW = UVW.z;
00185 if( test(aUVW) < 0 ) aUVW = aUVW * one_neg;
00186
00187 WP = P-W;
00188 PVW = cross_product(WP,WV);
00189 aPVW = PVW.z;
00190 if( test(aPVW) < 0 ) aPVW = aPVW * one_neg;
00191
00192 *u = aPVW/aUVW;
00193
00194 PWU = cross_product(WP,WU);
00195 aPWU = PWU.z;
00196 if( test(aPWU) < 0 ) aPWU = aPWU *one_neg;
00197
00198 *v = aPWU/aUVW;
00199
00200 UV = V-U;
00201 UP = P-U;
00202
00203 PUV = cross_product(UV,UP);
00204 aPUV = PUV.z;
00205 if( test(aPUV) < 0 ) aPUV = aPUV * one_neg;
00206
00207 *w = aPUV/aUVW;
00208 }
00209
00210 template GULAPI void BarycentricCoordinates(
00211 const point<rational>& P1, const point<rational>& P2,
00212 const point<rational>& P3, const point<rational>& Pin,
00213 rational *w, rational *u, rational *v );
00214
00215
00216
00217
00218
00219 template< class T >
00220 GULAPI void BarycentricCoordinates(
00221 const point<T>& P1, const point<T>& P2,
00222 const point<T>& P3, const point<T>& P12,
00223 const point<T>& P13, const point<T>& P23,
00224 const point<T>& N, const point<T>& Pin,
00225 T *w, T *u, T *v )
00226 {
00227 point2<T> W,U,V,P;
00228 point2<T> WU,WV,UV,WP,UP;
00229 point<T> UVW,PVW,PWU,PUV;
00230 T aUVW,aPVW,aPWU,aPUV;
00231 rational one_neg(ULong(1),-1);
00232
00233
00234 if( test(N.z) )
00235 {
00236 W.x = P1.x;
00237 W.y = P1.y;
00238
00239 U.x = P2.x;
00240 U.y = P2.y;
00241
00242 V.x = P3.x;
00243 V.y = P3.y;
00244
00245 WU.x = P12.x;
00246 WU.y = P12.y;
00247
00248 WV.x = P13.x;
00249 WV.y = P13.y;
00250
00251 UV.x = P23.x;
00252 UV.y = P23.y;
00253
00254 P.x = Pin.x;
00255 P.y = Pin.y;
00256 }
00257 else if( test(N.y) )
00258 {
00259 W.x = P1.x;
00260 W.y = P1.z;
00261
00262 U.x = P2.x;
00263 U.y = P2.z;
00264
00265 V.x = P3.x;
00266 V.y = P3.z;
00267
00268 WU.x = P12.x;
00269 WU.y = P12.z;
00270
00271 WV.x = P13.x;
00272 WV.y = P13.z;
00273
00274 UV.x = P23.x;
00275 UV.y = P23.z;
00276
00277 P.x = Pin.x;
00278 P.y = Pin.z;
00279 }
00280 else
00281 {
00282 gul::Assert<InternalError>( ndebug || test(N.x) );
00283
00284 W.x = P1.y;
00285 W.y = P1.z;
00286
00287 U.x = P2.y;
00288 U.y = P2.z;
00289
00290 V.x = P3.y;
00291 V.y = P3.z;
00292
00293 WU.x = P12.y;
00294 WU.y = P12.z;
00295
00296 WV.x = P13.y;
00297 WV.y = P13.z;
00298
00299 UV.x = P23.y;
00300 UV.y = P23.z;
00301
00302 P.x = Pin.y;
00303 P.y = Pin.z;
00304 }
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314 UVW = cross_product(WU,WV);
00315 aUVW = UVW.z;
00316 if( test(aUVW) < 0 ) aUVW = aUVW * one_neg;
00317
00318 WP = P-W;
00319 PVW = cross_product(WP,WV);
00320 aPVW = PVW.z;
00321 if( test(aPVW) < 0 ) aPVW = aPVW * one_neg;
00322
00323 *u = aPVW/aUVW;
00324
00325 PWU = cross_product(WP,WU);
00326 aPWU = PWU.z;
00327 if( test(aPWU) < 0 ) aPWU = aPWU *one_neg;
00328
00329 *v = aPWU/aUVW;
00330
00331 UP = P-U;
00332 PUV = cross_product(UV,UP);
00333 aPUV = PUV.z;
00334 if( test(aPUV) < 0 ) aPUV = aPUV * one_neg;
00335
00336 *w = aPUV/aUVW;
00337 }
00338
00339 template GULAPI void BarycentricCoordinates(
00340 const point<rational>& P1, const point<rational>& P2,
00341 const point<rational>& P3, const point<rational>& P12,
00342 const point<rational>& P13, const point<rational>& P23,
00343 const point<rational>& N, const point<rational>& Pin,
00344 rational *w, rational *u, rational *v );
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 GULAPI int IntersectTriangles(
00356 const triangle<rational>& tri0,
00357 const triangle<rational>& tri1,
00358 Ptr<point<rational> >& retP )
00359 {
00360 point<rational> v11,v12,v13,n,v21,v22,v23,g,t;
00361 point<rational> P1,P2,P3;
00362 rational A,B,C,R,a,b,c,r,N;
00363 rational one(ULong(1)),zero;
00364 rational lambda,mu,l1,r1,l2,r2,l;
00365 bool result,init;
00366 int dim;
00367
00368
00369
00370
00371 v11 = tri0.P2 - tri0.P1;
00372 v12 = tri0.P3 - tri0.P1;
00373 v13 = tri0.P3 - tri0.P2;
00374 n = cross_product( v11, v12 );
00375 A = n.x; B = n.y; C = n.z;
00376 R = n * tri0.P1;
00377
00378
00379
00380 v21 = tri1.P2 - tri1.P1;
00381 v22 = tri1.P3 - tri1.P1;
00382 v23 = tri1.P3 - tri1.P2;
00383 n = cross_product( v21, v22 );
00384 a = n.x; b = n.y; c = n.z;
00385 r = n * tri1.P1;
00386
00387
00388
00389
00390
00391
00392
00393
00394 N = B*c - b*C;
00395 dim = 1;
00396
00397 if( test(N) )
00398 {
00399 g.x = zero;
00400 g.y = (R*c - r*C)/N;
00401 g.z = (B*r - b*R)/N;
00402 t.x = one;
00403 t.y = (C*a - c*A)/N;
00404 t.z = (A*b - a*B)/N;
00405 }
00406 else
00407 {
00408 N = A*c - a*C;
00409
00410 if( test(N) )
00411 {
00412 g.x = (R*c - r*C) / N;
00413 g.y = zero;
00414 g.z = (A*r - a*R) / N;
00415 t.x = (C*b - c*B) / N;
00416 t.y = one;
00417 t.z = (B*a - b*A) / N;
00418 }
00419 else
00420 {
00421 N = A*b - a*B;
00422
00423 if( test(N) )
00424 {
00425 g.x = (R*b - r*B) / N;
00426 g.y = (A*r - a*R) / N;
00427 g.z = zero;
00428 t.x = (B*c - b*C) / N;
00429 t.y = (C*a - c*A) / N;
00430 t.z = one;
00431 }
00432 else
00433 {
00434
00435
00436 if( test(A*tri1.P1.x + B*tri1.P1.y + C*tri1.P1.z - R) )
00437 return 0;
00438 else
00439 dim = 2;
00440 }
00441 }
00442 }
00443
00444 if( dim == 1 )
00445 {
00446
00447
00448
00449
00450
00451 P1 = tri0.P1;
00452 P2 = tri0.P2;
00453 P3 = tri0.P3;
00454
00455
00456
00457
00458 init = false;
00459 result = RegularIntersectLines( P1, v11, g, t, &lambda, &mu );
00460 if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00461 {
00462 init = true;
00463 l1 = mu;
00464 r1 = mu;
00465 }
00466
00467
00468
00469 result = RegularIntersectLines( P1, v12, g, t, &lambda, &mu );
00470 if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00471 {
00472 if( init )
00473 {
00474 if( compare(mu,l1) < 0 ) l1 = mu;
00475 else if( compare(mu,r1) > 0 ) r1 = mu;
00476 }
00477 else
00478 {
00479 init = true;
00480 l1 = mu;
00481 r1 = mu;
00482 }
00483 }
00484
00485
00486
00487 result = RegularIntersectLines( P2, v13, g, t, &lambda, &mu );
00488 if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00489 {
00490 if( init )
00491 {
00492 if( compare(mu,l1) < 0 ) l1 = mu;
00493 else if( compare(mu,r1) > 0 ) r1 = mu;
00494 }
00495 else
00496 {
00497 init = true;
00498 l1 = mu;
00499 r1 = mu;
00500 }
00501 }
00502
00503 if( !init )
00504 return 0;
00505 cout << "l1 = " << l1.dump() << ", r1 = " << r1.dump() << "\n";
00506
00507
00508
00509
00510
00511
00512 P1 = tri1.P1;
00513 P2 = tri1.P2;
00514 P3 = tri1.P3;
00515
00516
00517
00518
00519 init = false;
00520 result = RegularIntersectLines( P1, v21, g, t, &lambda, &mu );
00521 if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00522 {
00523 init = true;
00524 l2 = mu;
00525 r2 = mu;
00526 }
00527
00528
00529
00530 result = RegularIntersectLines( P1, v22, g, t, &lambda, &mu );
00531 if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00532 {
00533 if( init )
00534 {
00535 if( compare(mu,l2) < 0 ) l2 = mu;
00536 else if( compare(mu,r2) > 0 ) r2 = mu;
00537 }
00538 else
00539 {
00540 init = true;
00541 l2 = mu;
00542 r2 = mu;
00543 }
00544 }
00545
00546
00547
00548 result = RegularIntersectLines( P2, v23, g, t, &lambda, &mu );
00549 if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00550 {
00551 if( init )
00552 {
00553 if( compare(mu,l2) < 0 ) l2 = mu;
00554 else if( compare(mu,r2) > 0 ) r2 = mu;
00555 }
00556 else
00557 {
00558 init = true;
00559 l2 = mu;
00560 r2 = mu;
00561 }
00562 }
00563
00564 if( !init )
00565 return 0;
00566 cout << "l2 = " << l2.dump() << ", r2 = " << r2.dump() << "\n";
00567
00568
00569
00570
00571 if(compare(l1,l2)>0) l = l1; else l = l2;
00572 if(compare(r1,r2)<0) r = r1; else r = r2;
00573 if( compare(r,l) < 0 )
00574 return 0;
00575
00576
00577 retP[0] = g + l * t;
00578 retP[1] = g + r * t;
00579 return 2;
00580 }
00581
00582 return 0;
00583 }
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593 GULAPI bool IntersectTrianglesUV(
00594 const triangle<rational>& tri1, const triangle2<rational>& tri1uv,
00595 const triangle<rational>& tri2, const triangle2<rational>& tri2uv,
00596 point2<rational> *S1uv, int *S1flag,
00597 point2<rational> *S2uv, int *S2flag,
00598 point<rational> *S )
00599 {
00600 point<rational> v11,v12,v13,n,N,v21,v22,v23,g,t;
00601 point<rational> P1,P2,P3;
00602 rational A,B,C,R,a,b,c,r,T;
00603 rational one(ULong(1)),zero;
00604 rational lambda,mu,l1,r1,l2,r2,l;
00605 bool result,init;
00606 int dim;
00607 bool inter1_12, inter1_13, inter1_23;
00608 bool inter2_12, inter2_13, inter2_23;
00609 int mid1,mid2,t0,t1,sign,side;
00610 point<rational> Pc1,Pc2;
00611 rational u,v,w;
00612 point2<rational> Pc1uv, Pc2uv;
00613
00614
00615
00616
00617 v11 = tri1.P2 - tri1.P1;
00618 v12 = tri1.P3 - tri1.P1;
00619 v13 = tri1.P3 - tri1.P2;
00620 N = cross_product( v11, v12 );
00621 A = N.x; B = N.y; C = N.z;
00622 R = N * tri1.P1;
00623
00624
00625
00626
00627 v21 = tri2.P2 - tri2.P1;
00628 v22 = tri2.P3 - tri2.P1;
00629 v23 = tri2.P3 - tri2.P2;
00630 n = cross_product( v21, v22 );
00631 a = n.x; b = n.y; c = n.z;
00632 r = n * tri2.P1;
00633
00634
00635
00636
00637
00638
00639
00640
00641 T = B*c - b*C;
00642 dim = 1;
00643
00644 if( test(T) )
00645 {
00646 g.x = zero;
00647 g.y = (R*c - r*C)/T;
00648 g.z = (B*r - b*R)/T;
00649 t.x = one;
00650 t.y = (C*a - c*A)/T;
00651 t.z = (A*b - a*B)/T;
00652 }
00653 else
00654 {
00655 T = A*c - a*C;
00656
00657 if( test(T) )
00658 {
00659 g.x = (R*c - r*C) / T;
00660 g.y = zero;
00661 g.z = (A*r - a*R) / T;
00662 t.x = (C*b - c*B) / T;
00663 t.y = one;
00664 t.z = (B*a - b*A) / T;
00665 }
00666 else
00667 {
00668 T = A*b - a*B;
00669
00670 if( test(T) )
00671 {
00672 g.x = (R*b - r*B) / T;
00673 g.y = (A*r - a*R) / T;
00674 g.z = zero;
00675 t.x = (B*c - b*C) / T;
00676 t.y = (C*a - c*A) / T;
00677 t.z = one;
00678 }
00679 else
00680 {
00681
00682
00683 if( test(A*tri2.P1.x + B*tri2.P1.y + C*tri2.P1.z - R) )
00684 return false;
00685 else
00686 dim = 2;
00687 }
00688 }
00689 }
00690
00691 if( dim == 1 )
00692 {
00693
00694
00695
00696
00697 mid1 = false;
00698 inter1_12 = false;
00699 inter1_13 = false;
00700 inter1_23 = false;
00701
00702
00703 P1 = tri1.P1;
00704 P2 = tri1.P2;
00705 P3 = tri1.P3;
00706
00707
00708
00709
00710 init = false;
00711 result = RegularIntersectLines( P1, v11, g, t, &lambda, &mu );
00712 if( result )
00713 {
00714 t0 = test(lambda);
00715 if( t0 >= 0 )
00716 {
00717 t1 = compare(lambda,one);
00718 if( t1 <= 0 )
00719 {
00720 inter1_12 = true;
00721
00722 init = true;
00723 l1 = mu;
00724 r1 = mu;
00725
00726 if( (t0 > 0) && (t1 < 0) )
00727 {
00728 mid1 = true;
00729 Pc1 = tri1.P1;
00730
00731
00732 }
00733 }
00734 }
00735 }
00736
00737
00738
00739 result = RegularIntersectLines( P1, v12, g, t, &lambda, &mu );
00740 if( result )
00741 {
00742 t0 = test(lambda);
00743 if( t0 >= 0 )
00744 {
00745 t1 = compare(lambda,one);
00746 if( t1 <= 0 )
00747 {
00748 inter1_13 = true;
00749
00750 if( init )
00751 {
00752 if( compare(mu,l1) < 0 ) l1 = mu;
00753 else if( compare(mu,r1) > 0 ) r1 = mu;
00754 }
00755 else
00756 {
00757 init = true;
00758 l1 = mu;
00759 r1 = mu;
00760 }
00761
00762 if( (t0 > 0) && (t1 < 0) )
00763 {
00764 mid1 = true;
00765 Pc1 = tri1.P1;
00766
00767
00768 }
00769 }
00770 }
00771 }
00772
00773
00774
00775 result = RegularIntersectLines( P2, v13, g, t, &lambda, &mu );
00776 if( result )
00777 {
00778 t0 = test(lambda);
00779 if( t0 >= 0 )
00780 {
00781 t1 = compare(lambda,one);
00782 if( t1 <= 0 )
00783 {
00784 inter1_23 = true;
00785
00786 if( init )
00787 {
00788 if( compare(mu,l1) < 0 ) l1 = mu;
00789 else if( compare(mu,r1) > 0 ) r1 = mu;
00790 }
00791 else
00792 {
00793 init = true;
00794 l1 = mu;
00795 r1 = mu;
00796 }
00797
00798 if( (t0 > 0) && (t1 < 0) )
00799 {
00800 mid1 = true;
00801 Pc1 = tri1.P2;
00802
00803
00804 }
00805 }
00806 }
00807 }
00808
00809 if( !init )
00810 return 0;
00811
00812
00813
00814
00815
00816
00817 mid2 = false;
00818 inter2_12 = false;
00819 inter2_13 = false;
00820 inter2_23 = false;
00821
00822
00823 P1 = tri2.P1;
00824 P2 = tri2.P2;
00825 P3 = tri2.P3;
00826
00827
00828
00829
00830 init = false;
00831 result = RegularIntersectLines( P1, v21, g, t, &lambda, &mu );
00832 if( result )
00833 {
00834 t0 = test(lambda);
00835 if( t0 >= 0 )
00836 {
00837 t1 = compare(lambda,one);
00838 if( t1 <= 0 )
00839 {
00840 inter2_12 = true;
00841
00842 init = true;
00843 l2 = mu;
00844 r2 = mu;
00845
00846 if( (t0 > 0) && (t1 < 0) )
00847 {
00848 mid2 = true;
00849 Pc2 = tri2.P1;
00850
00851
00852 }
00853 }
00854 }
00855 }
00856
00857
00858
00859 result = RegularIntersectLines( P1, v22, g, t, &lambda, &mu );
00860 if( result )
00861 {
00862 t0 = test(lambda);
00863 if( t0 >= 0 )
00864 {
00865 t1 = compare(lambda,one);
00866 if( t1 <= 0 )
00867 {
00868 inter2_13 = true;
00869
00870 if( init )
00871 {
00872 if( compare(mu,l2) < 0 ) l2 = mu;
00873 else if( compare(mu,r2) > 0 ) r2 = mu;
00874 }
00875 else
00876 {
00877 init = true;
00878 l2 = mu;
00879 r2 = mu;
00880 }
00881
00882 if( (t0 > 0) && (t1 < 0) )
00883 {
00884 mid2 = true;
00885 Pc2 = tri2.P1;
00886
00887
00888 }
00889 }
00890 }
00891 }
00892
00893
00894
00895 result = RegularIntersectLines( P2, v23, g, t, &lambda, &mu );
00896 if( result )
00897 {
00898 t0 = test(lambda);
00899 if( t0 >= 0 )
00900 {
00901 t1 = compare(lambda,one);
00902 if( t1 <= 0 )
00903 {
00904 inter2_23 = true;
00905
00906 if( init )
00907 {
00908 if( compare(mu,l2) < 0 ) l2 = mu;
00909 else if( compare(mu,r2) > 0 ) r2 = mu;
00910 }
00911 else
00912 {
00913 init = true;
00914 l2 = mu;
00915 r2 = mu;
00916 }
00917
00918 if( (t0 > 0) && (t1 < 0) )
00919 {
00920 mid2 = true;
00921 Pc2 = tri2.P2;
00922
00923
00924 }
00925 }
00926 }
00927 }
00928
00929 if( !init )
00930 return 0;
00931
00932
00933
00934
00935
00936 if(compare(l1,l2)>0) l = l1; else l = l2;
00937 if(compare(r1,r2)<0) r = r1; else r = r2;
00938 t1 = compare(r,l);
00939 if( t1 < 0 ) return false;
00940 else if( t1 == 0 )
00941 {
00942
00943 return false;
00944 }
00945
00946
00947 S[0] = g + l * t;
00948 S[1] = g + r * t;
00949
00950
00951 BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N,
00952 S[0], &u, &v, &w );
00953 S1uv[0] = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;
00954
00955 BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N,
00956 S[1], &u, &v, &w );
00957 S1uv[1] = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;
00958
00959
00960 if( mid1 )
00961 {
00962 sign = test(n * Pc1 - r);
00963
00964
00965
00966 BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N,
00967 Pc1, &u, &v, &w );
00968 Pc1uv = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;
00969
00970 side = test( cross_product(S1uv[1]-S1uv[0],Pc1uv-S1uv[0]).z );
00971
00972 gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
00973 if( side < 0 ) side = 0;
00974 else if( side > 0 ) side = 1;
00975 S1flag[side] = sign;
00976 S1flag[!side] = -sign;
00977 }
00978 else
00979 {
00980 if( inter1_12 && inter1_13 )
00981 Pc1 = tri1.P1;
00982 else if( inter1_12 && inter1_23 )
00983 Pc1 = tri1.P2;
00984 else
00985 {
00986 gul::Assert<InternalError>( ndebug || (inter1_13 && inter1_23));
00987 Pc1 = tri1.P3;
00988 }
00989
00990 sign = test(n * Pc1 - r);
00991
00992
00993
00994 BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N,
00995 Pc1, &u, &v, &w );
00996 Pc1uv = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;
00997
00998 side = test( cross_product(S1uv[1]-S1uv[0],Pc1uv-S1uv[0]).z );
00999
01000 gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
01001 S1flag[side] = sign;
01002 S1flag[!side] = 0;
01003 }
01004
01005
01006 BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n,
01007 S[0], &u, &v, &w );
01008 S2uv[0] = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;
01009
01010 BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n,
01011 S[1], &u, &v, &w );
01012 S2uv[1] = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;
01013
01014
01015 if( mid2 )
01016 {
01017 sign = test(N * Pc2 - R);
01018
01019
01020
01021 BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n,
01022 Pc2, &u, &v, &w );
01023 Pc2uv = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;
01024
01025 side = test( cross_product(S2uv[1]-S2uv[0],Pc2uv-S2uv[0]).z );
01026
01027 gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
01028 if( side < 0 ) side = 0;
01029 else if( side > 0 ) side = 1;
01030 S2flag[side] = sign;
01031 S2flag[!side] = -sign;
01032 }
01033 else
01034 {
01035 if( inter2_12 && inter2_13 )
01036 Pc2 = tri2.P1;
01037 else if( inter2_12 && inter2_23 )
01038 Pc2 = tri2.P2;
01039 else
01040 {
01041 gul::Assert<InternalError>( ndebug || (inter2_13 && inter2_23));
01042 Pc2 = tri1.P3;
01043 }
01044
01045 sign = test(N * Pc2 - R);
01046
01047
01048
01049 BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n,
01050 Pc2, &u, &v, &w );
01051 Pc2uv = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;
01052
01053 side = test( cross_product(S2uv[1]-S2uv[0],Pc2uv-S2uv[0]).z );
01054
01055 gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
01056 S2flag[side] = sign;
01057 S2flag[!side] = 0;
01058 }
01059
01060
01061 return true;
01062 }
01063
01064 return false;
01065 }
01066
01067 }