Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members  

guar_intersect.cpp

Go to the documentation of this file.
00001 /* LIBGUL - Geometry Utility Library
00002  * Copyright (C) 1998-1999 Norbert Irmer
00003  *
00004  * This library is free software; you can redistribute it and/or
00005  * modify it under the terms of the GNU Library General Public
00006  * License as published by the Free Software Foundation; either
00007  * version 2 of the License, or (at your option) any later version.
00008  *
00009  * This library is distributed in the hope that it will be useful,
00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * Library General Public License for more details.
00013  *
00014  * You should have received a copy of the GNU Library General Public
00015  * License along with this library; if not, write to the
00016  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
00017  * Boston, MA 02111-1307, USA.
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   calculate the intersection point of two lines (in 3D). It must exist,
00044   or the result will be wrong.
00045 ----------------------------------------------------------------------- */
00046 bool RegularIntersectLines(
00047   const point<rational>& A,     // point on line 1
00048   const point<rational>& B,     // direction vector of line 1
00049   const point<rational>& a,     // point on line 2
00050   const point<rational>& b,     // direction vector of line 2
00051   rational *lambda,  // parameter value of intersect. point (for line 1)
00052   rational *mu       // parameter value of intersect. point (for line 2)
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   // more convenient notation
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) ) // calculation without z (for exact arithmetic the
00068   {              // goodness of N1 doesn't matter)
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) )  // calculation without y (for exact arithmetic the
00077     {               // goodness of N2 doesn't matter)
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) ) // calculation without y (for exact arithmetic the 
00086       {              // goodness of N3 doesn't matter)
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;  // no intersection point or lines colinear
00093       }      
00094     }
00095   }  
00096   return true;
00097 }
00098 
00099 /*-----------------------------------------------------------------------
00100   Estimate the (u,v) parameter values of a point P for a NURBS surface F,
00101   assuming that A = F(0,0), B = F(1,0), C = F(0,1)
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   // calculate barycentric coordinates in xy-plane
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) ) // calculate barycentric coordinates in xz-plane
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 // calculate barycentric coordinates in yz-plane
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 // template instantiation
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   Estimate the (u,v) parameter values of a point P for a NURBS surface F,
00217   assuming that A = F(0,0), B = F(1,0), C = F(0,1)
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   // calculate barycentric coordinates in xy-plane
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) ) // calculate barycentric coordinates in xz-plane
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 // calculate barycentric coordinates in yz-plane
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   cout << "N = (" << N << ")\n";
00308   cout << "P12 = (" << P12 << "), ";
00309   cout << "WU = (" << WU << ")\n";
00310   cout << "P13 = (" << P13 << "), ";
00311   cout << "WV = (" << WV << ")\n";
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 // template instantiation
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   calculates the endpoints of the intersection line segment of two
00349   triangles (the cases with both triangles being coplanar is not 
00350   handled).
00351   it also checks if the two halves of the triangles left and right
00352   of the intersection line are lying above or below the other
00353   triangle
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   // calculate plane equation for 1.triangle
00369   // (Ax + By + Cz = R)
00370 
00371   v11 = tri0.P2 - tri0.P1;       // 1. direction vector
00372   v12 = tri0.P3 - tri0.P1;       // 2. direction vector
00373   v13 = tri0.P3 - tri0.P2;       // (needed later)
00374   n = cross_product( v11, v12 ); // normal vector of the plane
00375   A = n.x;  B = n.y;  C = n.z;
00376   R = n * tri0.P1;               // right side of plane equation
00377 
00378   // calculate plane equation for 2.triangle
00379   // (ax + by + cz = r)
00380   v21 = tri1.P2 - tri1.P1;       // 1. direction vector
00381   v22 = tri1.P3 - tri1.P1;       // 2. direction vector
00382   v23 = tri1.P3 - tri1.P2;       // (needed later)
00383   n = cross_product( v21, v22 ); // normal vector of the plane
00384   a = n.x;  b = n.y;  c = n.z;
00385   r = n * tri1.P1;               // right side of plane equation
00386 
00387   // calculate intersection line
00388   // (G : x = g + lambda * t)
00389 
00390   // if intersection line is not parallel to the X-axis,
00391   // set g.x = 0 and t.x = 1
00392   // ==> g.y, g.z and t.x, t.z
00393    
00394   N = B*c - b*C;    // determinant
00395   dim = 1;          // assume intersection has dimension 1
00396 
00397   if( test(N) )     // if not parallel to the X-axis 
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) )   // if not parallel to the Y-axis
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) )   // if not parallel to the Z-axis
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        // normal vectors of the two planes are linear dependent
00433       {
00434         // check if one point of the second triangle lies within the plane
00435         // of the first triangle
00436         if( test(A*tri1.P1.x + B*tri1.P1.y + C*tri1.P1.z - R) )
00437           return 0; // no intersection
00438         else
00439           dim = 2;
00440       }  
00441     }
00442   }
00443 
00444   if( dim == 1 )
00445   {
00446     // intersect intersection line of the planes with the border of the first
00447     // triangle. the result is an interval [l1,r1] in the parametric domain
00448     // of the intersection line.
00449 
00450     // convenient notation
00451     P1 = tri0.P1;
00452     P2 = tri0.P2;
00453     P3 = tri0.P3;
00454 
00455     // intersection line: X = g + mu * t
00456     // line through triangle edge P1P2:  X = P1 + lambda * v11
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     // line through triangle edge P1P3:  X = P1 + lambda * v12
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     // line through triangle edge P2P3:  X = P2 + lambda * v13
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     // intersect intersection line of the planes with the border of the second
00508     // triangle. the result is an interval [l2,r2] in the parametric domain
00509     // of the intersection line.
00510 
00511     // convenient notation
00512     P1 = tri1.P1;
00513     P2 = tri1.P2;
00514     P3 = tri1.P3;
00515 
00516     // intersection line: X = g + mu * t
00517     // line through triangle edge P1P2:  X = P1 + lambda * v21
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     // line through triangle edge P1P3:  X = P1 + lambda * v22
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     // line through triangle edge P2P3:  X = P2 + lambda * v23
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     // intersection of the intervals [l1,r1] and [l2,r2]  delivers the 
00569     // resulting segment of the intersection line 
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     // calculate the endpoints of the intersection segment
00577     retP[0] = g + l * t;
00578     retP[1] = g + r * t;
00579     return 2;
00580   }
00581 
00582   return 0;
00583 }
00584 
00585 /*----------------------------------------------------------------
00586   calculates the endpoints of the intersection line segment of two
00587   triangles (the cases with both triangles being coplanar or just
00588   touching in one point is not handled).
00589   it also checks if the two halves of the triangles left and right
00590   of the intersection line are lying above or below the other
00591   triangle
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   // calculate plane equation for 1.triangle
00615   // (Ax + By + Cz = R)
00616 
00617   v11 = tri1.P2 - tri1.P1;       // 1. direction vector
00618   v12 = tri1.P3 - tri1.P1;       // 2. direction vector
00619   v13 = tri1.P3 - tri1.P2;       // (needed later)
00620   N = cross_product( v11, v12 ); // normal vector of the plane
00621   A = N.x;  B = N.y;  C = N.z;
00622   R = N * tri1.P1;               // right side of plane equation
00623   // cout << "N = (" << N << ")\n";
00624 
00625   // calculate plane equation for 2.triangle
00626   // (ax + by + cz = r)
00627   v21 = tri2.P2 - tri2.P1;       // 1. direction vector
00628   v22 = tri2.P3 - tri2.P1;       // 2. direction vector
00629   v23 = tri2.P3 - tri2.P2;       // (needed later)
00630   n = cross_product( v21, v22 ); // normal vector of the plane
00631   a = n.x;  b = n.y;  c = n.z;
00632   r = n * tri2.P1;               // right side of plane equation
00633 
00634   // calculate intersection line
00635   // (G : x = g + lambda * t)
00636 
00637   // if intersection line is not parallel to the X-axis,
00638   // set g.x = 0 and t.x = 1
00639   // ==> g.y, g.z and t.x, t.z
00640    
00641   T = B*c - b*C;    // determinant
00642   dim = 1;          // assume intersection has dimension 1
00643 
00644   if( test(T) )     // if not parallel to the X-axis 
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) )   // if not parallel to the Y-axis
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) )   // if not parallel to the Z-axis
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        // normal vectors of the two planes are linear dependent
00680       {
00681         // check if one point of the second triangle lies within the plane
00682         // of the first triangle
00683         if( test(A*tri2.P1.x + B*tri2.P1.y + C*tri2.P1.z - R) )
00684           return false; // no intersection
00685         else
00686           dim = 2;
00687       }  
00688     }
00689   }
00690 
00691   if( dim == 1 )
00692   {
00693     // intersect intersection line of the planes with the border of the first
00694     // triangle. the result is an interval [l1,r1] in the parametric domain
00695     // of the intersection line.
00696 
00697     mid1 = false; 
00698     inter1_12 = false;
00699     inter1_13 = false;
00700     inter1_23 = false;
00701 
00702     // convenient notation
00703     P1 = tri1.P1;
00704     P2 = tri1.P2;
00705     P3 = tri1.P3;
00706 
00707     // intersection line: X = g + mu * t
00708     // line through triangle edge P1P2:  X = P1 + lambda * v11
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;  // will be used to classify the two halves of
00730                             // tri 0 left and right of the intersection line
00731                             // as lying above or below tri2
00732           }
00733         }
00734       }
00735     }
00736 
00737     // line through triangle edge P1P3:  X = P1 + lambda * v12
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;  // will be used to classify the two halves of
00766                             // tri 0 left and right of the intersection line
00767                             // as lying above or below tri2
00768           }
00769         }
00770       }
00771     }
00772 
00773     // line through triangle edge P2P3:  X = P2 + lambda * v13
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;  // will be used to classify the two halves of
00802                             // tri 0 left and right of the intersection line
00803                             // as lying above or below tri2
00804           }
00805         }
00806       }
00807     }
00808 
00809     if( !init )
00810       return 0;
00811     // cout << "l1 = " << l1.dump() << ", r1 = " << r1.dump() << "\n";
00812 
00813     // intersect intersection line of the planes with the border of the second
00814     // triangle. the result is an interval [l2,r2] in the parametric domain
00815     // of the intersection line.
00816 
00817     mid2 = false;
00818     inter2_12 = false;
00819     inter2_13 = false;
00820     inter2_23 = false;
00821 
00822     // convenient notation
00823     P1 = tri2.P1;
00824     P2 = tri2.P2;
00825     P3 = tri2.P3;
00826 
00827     // intersection line: X = g + mu * t
00828     // line through triangle edge P1P2:  X = P1 + lambda * v21
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;   // will be used to classify the two halves of
00850                             // tri 0 left and right of the intersection line
00851                             // as lying above or below tri2
00852           }
00853         }
00854       }
00855     }
00856   
00857     // line through triangle edge P1P3:  X = P1 + lambda * v22
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;  // will be used to classify the two halves of
00886                             // tri 0 left and right of the intersection line
00887                             // as lying above or below tri2
00888           }
00889         }
00890       }
00891     }
00892 
00893     // line through triangle edge P2P3:  X = P2 + lambda * v23
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;  // will be used to classify the two halves of
00922                             // tri 0 left and right of the intersection line
00923                             // as lying above or below tri2
00924           }
00925         }
00926       }
00927     }
00928 
00929     if( !init )
00930       return 0;
00931     // cout << "l2 = " << l2.dump() << ", r2 = " << r2.dump() << "\n";
00932 
00933     // intersection of the intervals [l1,r1] and [l2,r2]  delivers the 
00934     // resulting segment of the intersection line 
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 ) // single point intersection
00941     {
00942       // retP[0] = g + l * t;
00943       return false;
00944     }
00945 
00946     // calculate the endpoints of the intersection segment
00947     S[0] = g + l * t;
00948     S[1] = g + r * t;
00949 
00950     // in triangle coordinates of triangle 1
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     // calculate the edge/side pairs of the intersection segment in tri1
00960     if( mid1 )
00961     {
00962       sign = test(n * Pc1 - r);
00963 
00964       // this is a bit of a overkill since the coordinates of 'Pc1' in the 
00965       // parametric domain are known
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       // this is a bit of a overkill since the coordinates of 'Pc1' in the 
00993       // parametric domain are known
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;  // undetermined
01003     }
01004 
01005     // in triangle coordinates of triangle 2
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     // calculate the edge/side pairs of the intersection segment in tri2
01015     if( mid2 )
01016     {
01017       sign = test(N * Pc2 - R);
01018 
01019       // this is a bit of a overkill since the coordinates of 'Pc1' in the 
01020       // parametric domain are known
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       // this is a bit of a overkill since the coordinates of 'Pc1' in the 
01048       // parametric domain are known
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;  // undetermined
01058     }
01059 
01060     // two intersection points
01061     return true;
01062   }
01063 
01064   return false;
01065 }
01066 
01067 }

Generated on Mon Jan 21 04:17:18 2002 for GUL 0.6 - Geometry Utility Library by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001