src/INTERP_KERNEL/TransformedTriangleIntersect.cxx

   1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "TransformedTriangle.hxx"
  21 #include <iostream>
  22 #include <fstream>
  23 #include <cassert>
  24 #include <cmath>
  25 #include "VectorUtils.hxx"
  26
  27 namespace INTERP_KERNEL
  28 {
  29
  30   // ----------------------------------------------------------------------------------
  31   //  Correspondance tables describing all the variations of formulas.
  32   // ----------------------------------------------------------------------------------
  33
  34   /// \brief Correspondance between facets and double products.
  35   ///
  36   /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
  37   const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_SEG_FACET_INTERSECTION[12] =
  38     {
  39       C_XH, C_XY, C_ZX, // OYZ
  40       C_YH, C_YZ, C_XY, // OZX
  41       C_ZH, C_ZX, C_YZ, // OXY
  42       C_XH, C_YH, C_ZH  // XYZ
  43     };
  44
  45   /// \brief Signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION.
  46   ///
  47   /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
  48   const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] =
  49     {
  50       1.0, 1.0, -1.0,
  51       1.0, 1.0, -1.0,
  52       1.0, 1.0, -1.0,
  53       1.0, 1.0,  1.0
  54     };
  55
  56   /// \brief Coordinates of corners of tetrahedron.
  57   ///
  58   /// Use 3*Corner + coordinate as index
  59   const double TransformedTriangle::COORDS_TET_CORNER[12] =
  60     {
  61       0.0, 0.0, 0.0,
  62       1.0, 0.0, 0.0,
  63       0.0, 1.0, 0.0,
  64       0.0, 0.0, 1.0
  65     };
  66
  67   /// \brief Indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
  68   /// for the calculation of the coordinates (x,y,z) of the intersection points
  69   /// for Segment-Facet and Segment-Edge intersections.
  70   ///
  71   /// Use 3*facet + coordinate as index. -1 indicates that the coordinate is 0.
  72   const int TransformedTriangle::DP_INDEX[12] =
  73     {
  74       // x, y, z
  75       -1, 1, 2,  // OYZ
  76       5, -1, 4,  // OZX
  77       7, 8, -1,  // OXY
  78       9, 10, 11  // XYZ
  79     };
  80
  81   /// \brief Correspondance edge - corners.
  82   ///
  83   /// Gives the two corners associated with each edge
  84   /// Use 2*edge + {0, 1} as index
  85   const TransformedTriangle::TetraCorner TransformedTriangle::CORNERS_FOR_EDGE[12] =
  86     {
  87       O, X, // OX
  88       O, Y, // OY
  89       O, Z, // OZ
  90       X, Y, // XY
  91       Y, Z, // YZ
  92       Z, X  // ZX
  93     };
  94
  95   /// \brief Correspondance edge - facets.
  96   ///
  97   /// Gives the two facets shared by and edge. Use 2*facet + {0, 1} as index
  98   const TransformedTriangle::TetraFacet TransformedTriangle::FACET_FOR_EDGE[12] =
  99     {
 100       OXY, OZX, // OX
 101       OXY, OYZ, // OY
 102       OZX, OYZ, // OZ
 103       OXY, XYZ, // XY
 104       OYZ, XYZ, // YZ
 105       OZX, XYZ  // ZX
 106     };
 107
 108   /// \brief Correspondance corners - edges.
 109   ///
 110   /// Gives edges meeting at a given corner. Use 3*corner + {0,1,2} as index
 111   const TransformedTriangle::TetraEdge TransformedTriangle::EDGES_FOR_CORNER[12] =
 112     {
 113       OX, OY, OZ, // O
 114       OX, XY, ZX, // X
 115       OY, XY, YZ, // Y
 116       OZ, ZX, YZ  // Z
 117     };
 118
 119   /// \brief Double products to use in halfstrip intersection tests.
 120   ///
 121   /// Use 4*(offset_edge) + {0,1,2,3} as index. offset_edge = edge - 3  (so that XY -> 0, YZ -> 1, ZX -> 2)
 122   /// Entries with offset 0 and 1 are for the first condition (positive product)
 123   /// and those with offset 2 and 3 are for the second condition (negative product).
 124   const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_HALFSTRIP_INTERSECTION[12] =
 125     {
 126       C_10, C_01, C_ZH, C_10, // XY
 127       C_01, C_XY, C_XH, C_01, // YZ
 128       C_XY, C_10, C_YH, C_XY  // ZX
 129     };
 130
 131   /// \brief Double products to use in segment-ray test.
 132   ///
 133   /// Use 7*corner_offset + {0,1,2,3,4,5,6} as index. corner_offset = corner - 1 (so that X -> 0, Y-> 1, Z->2)
 134   /// Entries with offset 0 are for first condition (zero double product) and the rest are for condition 3 (in the same
 135   /// order as in the article)
 136   const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] =
 137     {
 138       C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
 139       C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y
 140       C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01  // Z
 141     };
 142
 143   /**
 144    * Calculates the point of intersection between the given edge of the tetrahedron and the
 145    * triangle PQR. (Grandy, eq [22])
 146    *
 147    * @pre   testSurfaceEdgeIntersection(edge) returns true
 148    * @param edge   edge of tetrahedron
 149    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 150    */
 151   void TransformedTriangle::calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const
 152   {
 153     assert(edge < H01);
 154
 155     // barycentric interpolation between points A and B
 156     // : (x,y,z)* = (1-alpha)*A + alpha*B where
 157     // alpha = t_A / (t_A - t_B)
 158
 159     const TetraCorner corners[2] =
 160       {
 161         CORNERS_FOR_EDGE[2*edge],
 162         CORNERS_FOR_EDGE[2*edge + 1]
 163       };
 164
 165     // calculate alpha
 166     const double tA = calcStableT(corners[0]);
 167     const double tB = calcStableT(corners[1]);
 168     const double alpha = tA / (tA - tB);
 169
 170     // calculate point
 171     LOG(4, "corner A = " << corners[0] << " corner B = " << corners[1] );
 172     LOG(4, "tA = " << tA << " tB = " << tB << " alpha= " << alpha );
 173     for(int i = 0; i < 3; ++i)
 174       {
 175
 176         pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] +
 177           alpha * COORDS_TET_CORNER[3*corners[1] + i];
 178 #if 0
 179         pt[i] = (1 - alpha) * getCoordinateForTetCorner<corners[0], i>() +
 180           alpha * getCoordinateForTetCorner<corners[0], i>();
 181 #endif
 182         LOG(6, pt[i] );
 183         assert(pt[i] >= 0.0);
 184         assert(pt[i] <= 1.0);
 185       }
 186   }
 187
 188   /**
 189    * Calculates the point of intersection between the given segment of the triangle
 190    * and the given facet of the tetrahedron. (Grandy, eq. [23])
 191    *
 192    * @pre   testSurfaceEdgeIntersection(seg, facet) returns true
 193    *
 194    * @param seg    segment of the triangle
 195    * @param facet  facet of the tetrahedron
 196    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 197    */
 198   void TransformedTriangle::calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const
 199   {
 200     // calculate s
 201     double s = 0.0;
 202     for(int i = 0; i < 3; ++i)
 203       {
 204         const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[3*facet + i];
 205         const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + i];
 206         s -= sign * calcStableC(seg, dp);
 207       }
 208
 209     assert(s != 0.0);
 210
 211     // calculate coordinates of intersection point
 212     for(int i = 0 ; i < 3; ++i)
 213       {
 214         const int dpIdx = DP_INDEX[3*facet + i];
 215
 216         if(dpIdx < 0)
 217           {
 218             pt[i] = 0.0;
 219           }
 220         else
 221           {
 222             const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
 223             const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
 224             pt[i] = -( sign * calcStableC(seg, dp) ) / s;
 225
 226             LOG(4, "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i]  );
 227             LOG(4, "c(" << seg << ", " << dp << ") = " <<  sign * calcStableC(seg, dp) );
 228             assert(pt[i] >= 0.0);
 229             assert(pt[i] <= 1.0);
 230           }
 231       }
 232
 233   }
 234
 235   /**
 236    * Tests if the given segment of the triangle intersects the given edge of the tetrahedron (Grandy, eq. [20]
 237    * If the OPTIMIZE is defined, it does not do the test the double product that should be zero.
 238    * @param seg    segment of the triangle
 239    * @param edge   edge of tetrahedron
 240    * @return      true if the segment intersects the edge
 241    */
 242   bool TransformedTriangle::testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const
 243   {
 244       {
 245         // check condition that the double products for one of the two
 246         // facets adjacent to the edge has a positive product
 247         bool isFacetCondVerified = false;
 248         TetraFacet facet[2];
 249         for(int i = 0 ; i < 2 ; ++i)
 250           {
 251             facet[i] = FACET_FOR_EDGE[2*edge + i];
 252
 253             // find the two c-values -> the two for the other edges of the facet
 254             int idx1 = 0 ;
 255             int idx2 = 1;
 256             DoubleProduct dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1];
 257             DoubleProduct dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
 258
 259             if(dp1 == DoubleProduct( edge ))
 260               {
 261                 idx1 = 2;
 262                 dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1];
 263               }
 264             else if(dp2 == DoubleProduct( edge ))
 265               {
 266                 idx2 = 2;
 267                 dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
 268               }
 269
 270             const double c1 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1]*calcStableC(seg, dp1);
 271             const double c2 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2]*calcStableC(seg, dp2);
 272
 273             //isFacetCondVerified = isFacetCondVerified || c1*c2 > 0.0;
 274             if(c1*c2 > 0.0)
 275               {
 276                 isFacetCondVerified = true;
 277               }
 278           }
 279
 280         if(!isFacetCondVerified)
 281           {
 282             return false;
 283           }
 284         else
 285           {
 286             return testSegmentIntersectsFacet(seg, facet[0]) || testSegmentIntersectsFacet(seg, facet[1]);
 287           }
 288       }
 289   }
 290
 291   /**
 292    * Calculates the point of intersection between the given segment of the triangle
 293    * and the given edge of the tetrahedron. (Grandy, eq. [25])
 294    *
 295    * @pre   testSegmentEdgeIntersection(seg, edge) returns true
 296    *
 297    * @param seg    segment of the triangle
 298    * @param edge   edge of the tetrahedron
 299    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 300    */
 301   void TransformedTriangle::calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const
 302   {
 303     assert(edge < H01);
 304
 305     // get the two facets associated with the edge
 306     static const TetraFacet FACETS_FOR_EDGE[12] =
 307       {
 308         OXY, OZX, // OX
 309         OXY, OYZ, // OY
 310         OZX, OYZ, // OZ
 311         OXY, XYZ, // XY
 312         OYZ, XYZ, // YZ
 313         OZX, XYZ  // ZX
 314       };
 315
 316     const TetraFacet facets[2] =
 317       {
 318         FACETS_FOR_EDGE[2*edge],
 319         FACETS_FOR_EDGE[2*edge + 1]
 320       };
 321
 322     // calculate s for the two edges
 323     double s[2];
 324     for(int i = 0; i < 2; ++i)
 325       {
 326         s[i] = 0.0;
 327         for(int j = 0; j < 3; ++j)
 328           {
 329             const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[3*facets[i] + j];
 330             const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[3*facets[i] + j];
 331             s[i] += sign * calcStableC(seg, dp);
 332           }
 333       }
 334
 335     // calculate denominator
 336     const double denominator = s[0]*s[0] + s[1]*s[1];
 337
 338     // calculate intersection point
 339     for(int i = 0; i < 3; ++i)
 340       {
 341         // calculate double product values for the two faces
 342         double c[2];
 343         for(int j = 0 ; j < 2; ++j)
 344           {
 345             const int dpIdx = DP_INDEX[3*facets[j] + i];
 346             const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
 347             const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
 348             c[j] = dpIdx < 0.0 ? 0.0 : sign * calcStableC(seg, dp);
 349           }
 350
 351         // pt[i] = (c1*s1 + c2*s2) / (s1^2 + s2^2)
 352
 353         pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator;
 354
 355         // strange bug with -O2 enabled : assertion fails when we don't have the following
 356         // trace - line
 357         //std::cout << "pt[i] = " << pt[i] << std::endl;
 358         //assert(pt[i] >= 0.0); // check we are in tetraeder
 359         //assert(pt[i] <= 1.0);
 360
 361       }
 362   }
 363
 364
 365   /**
 366    * Tests if the given segment of the triangle intersects the given corner of the tetrahedron.
 367    * (Grandy, eq. [21]). If OPTIMIZE is defined, the double products that should be zero are not verified.
 368    *
 369    * @param seg    segment of the triangle
 370    * @param corner corner of the tetrahedron
 371    * @return      true if the segment intersects the corner
 372    */
 373   bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
 374   {
 375
 376
 377     // facets meeting at a given corner
 378     static const TetraFacet FACETS_FOR_CORNER[12] =
 379       {
 380         OXY, OYZ, OZX, // O
 381         OZX, OXY, XYZ, // X
 382         OYZ, XYZ, OXY, // Y
 383         OZX, XYZ, OYZ  // Z
 384       };
 385
 386     // check segment intersect a facet
 387     for(int i = 0 ; i < 3 ; ++i)
 388       {
 389         const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i];
 390         if(testSegmentIntersectsFacet(seg, facet))
 391           {
 392             return true;
 393           }
 394       }
 395
 396     return false;
 397   }
 398
 399   /**
 400    * Tests if the given segment of the triangle intersects the half-strip above the
 401    * given edge of the h = 0 plane. (Grandy, eq. [30])
 402    *
 403    * @param seg    segment of the triangle
 404    * @param edge   edge of the h = 0 plane of the tetrahedron (XY, YZ, ZX)
 405    * @return      true if the upwards ray from the corner intersects the triangle.
 406    */
 407   bool TransformedTriangle::testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edge)
 408   {
 409     // get right index here to avoid "filling out" array
 410     const int edgeIndex = static_cast<int>(edge) - 3;
 411
 412     // double products used in test
 413     // products 1 and 2 for each edge -> first condition in Grandy [30]
 414     // products 3 and 4 for each edge -> third condition
 415     // NB : some uncertainty whether these last are correct
 416     // DP_FOR_HALFSTRIP_INTERSECTION
 417
 418     // facets to use in second condition (S_m)
 419     static const TetraFacet FACET_FOR_HALFSTRIP_INTERSECTION[3] =
 420       {
 421         NO_TET_FACET, // XY -> special case : test with plane H = 0
 422         OYZ, // YZ
 423         OZX  // ZX
 424       };
 425
 426     const double cVals[4] =
 427       {
 428         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]),
 429         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]),
 430         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 2]),
 431         calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 3])
 432       };
 433
 434     const TetraFacet facet =  FACET_FOR_HALFSTRIP_INTERSECTION[edgeIndex];
 435
 436
 437     // special case : facet H = 0
 438     const bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
 439     LOG(4, "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) );
 440     LOG(4, "c2 = " << cVals[2] << ", c3 = " << cVals[3] );
 441
 442     return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0);
 443   }
 444
 445   /**
 446    * Calculates the point of intersection between the given segment of the triangle
 447    * and the halfstrip above the given edge of the tetrahedron. (Grandy, eq. [31])
 448    *
 449    * @pre   testSegmentHalfstripIntersection(seg, edge) returns true
 450    *
 451    * @param seg    segment of the triangle
 452    * @param edge   edge of the tetrahedron defining the halfstrip
 453    * @param pt     array of three doubles in which to store the coordinates of the intersection point
 454    */
 455   void TransformedTriangle::calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const
 456   {
 457     assert(edge > OZ);
 458     assert(edge < H01);
 459
 460     // get right index here to avoid "filling out" array
 461     const int edgeIndex = static_cast<int>(edge) - 3;
 462     assert(edgeIndex >= 0);
 463     assert(edgeIndex < 3);
 464
 465     // Barycentric interpolation on the edge
 466     // for edge AB : (x,y,z)* = (1-alpha) * A + alpha * B
 467     // where alpha = cB / (cB - cA)
 468
 469     const double cA = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]);
 470     const double cB = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]);
 471     assert(cA != cB);
 472
 473     const double alpha = cA / (cA - cB);
 474
 475     for(int i = 0; i < 3; ++i)
 476       {
 477         const TetraCorner corners[2] =
 478           {
 479             CORNERS_FOR_EDGE[2*edge],
 480             CORNERS_FOR_EDGE[2*edge + 1]
 481           };
 482
 483         const double cornerCoords[2] =
 484           {
 485             COORDS_TET_CORNER[3*corners[0] + i],
 486             COORDS_TET_CORNER[3*corners[1] + i]
 487           };
 488
 489         pt[i] = (1 - alpha) * cornerCoords[0] + alpha * cornerCoords[1];
 490         LOG(6, pt[i] );
 491         assert(pt[i] >= 0.0);
 492         assert(pt[i] <= 1.0);
 493       }
 494     assert(epsilonEqualRelative(pt[0] + pt[1] + pt[2], 1.0));
 495   }
 496
 497   /**
 498    * Tests if the given segment of triangle PQR intersects the ray pointing
 499    * in the upwards z - direction from the given corner of the tetrahedron. (Grandy eq. [29])
 500    * If OPTIMIZE is defined, the double product that should be zero is not verified.
 501    *
 502    * @param seg    segment of the triangle PQR
 503    * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z)
 504    * @return      true if the upwards ray from the corner intersects the segment.
 505    */
 506   bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const
 507   {
 508     assert(corner == X || corner == Y || corner == Z);
 509     LOG(4, "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner );
 510
 511     // readjust index since O is not used
 512     const int cornerIdx = static_cast<int>(corner) - 1;
 513
 514     // facets to use
 515     //? not sure this is correct
 516     static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] =
 517       {
 518         OZX, // X
 519         OYZ, // Y
 520         OZX, // Z
 521       };
 522
 523     // cond 2
 524     const bool cond21 = testSegmentIntersectsFacet(seg, FIRST_FACET_SEGMENT_RAY_INTERSECTION[cornerIdx]);
 525     const bool cond22  = (corner == Z) ? testSegmentIntersectsFacet(seg, OYZ) : testSegmentIntersectsHPlane(seg);
 526
 527     if(!(cond21 || cond22))
 528       {
 529         LOG(4, "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 );
 530         return false;
 531       }
 532
 533     // cond 3
 534     const double cVals[6] =
 535       {
 536         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 1]),
 537         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 2]),
 538         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 3]),
 539         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 4]),
 540         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 5]),
 541         calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 6]),
 542       };
 543
 544     // cond. 3
 545     if(( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) >= 0.0)
 546       {
 547         LOG(4, "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5]  );
 548       }
 549     return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0;
 550
 551   }
 552
 553   // /////////////////////////////////////////////////////////////////////////////////
 554   //  Utility methods used in intersection tests                       ///////////////
 555   // /////////////////////////////////////////////////////////////////////////////////
 556   /**
 557    * Tests if the triangle PQR surrounds the axis on which the
 558    * given edge of the tetrahedron lies.
 559    *
 560    * @param edge   edge of tetrahedron
 561    * @return      true if PQR surrounds edge, false if not (see Grandy, eq. [53])
 562    */
 563   bool TransformedTriangle::testTriangleSurroundsEdge(const TetraEdge edge) const
 564   {
 565     // NB DoubleProduct enum corresponds to TetraEdge enum according to Grandy, table III
 566     // so we can use the edge directly
 567
 568     const double cPQ = calcStableC(PQ, DoubleProduct(edge));
 569     const double cQR = calcStableC(QR, DoubleProduct(edge));
 570     const double cRP = calcStableC(RP, DoubleProduct(edge));
 571
 572     LOG(5, "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" );
 573
 574     // if two or more c-values are zero we disallow x-edge intersection
 575     // Grandy, p.446
 576     const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
 577
 578     if(numZeros >= 2 )
 579       {
 580         LOG(5, "TriangleSurroundsEdge test fails due to too many 0 dp" );
 581       }
 582
 583     return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2;
 584   }
 585
 586 } // NAMESPACE