1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 #include "TransformedTriangle.hxx"
24 #include "VectorUtils.hxx"
26 namespace INTERP_KERNEL
29 // ----------------------------------------------------------------------------------
30 // Correspondance tables describing all the variations of formulas.
31 // ----------------------------------------------------------------------------------
33 /// \brief Correspondance between facets and double products.
35 /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
36 const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_SEG_FACET_INTERSECTION[12] =
38 C_XH, C_XY, C_ZX, // OYZ
39 C_YH, C_YZ, C_XY, // OZX
40 C_ZH, C_ZX, C_YZ, // OXY
41 C_XH, C_YH, C_ZH // XYZ
44 /// \brief Signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION.
46 /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
47 const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] =
55 /// \brief Coordinates of corners of tetrahedron.
57 /// Use 3*Corner + coordinate as index
58 const double TransformedTriangle::COORDS_TET_CORNER[12] =
66 /// \brief Indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
67 /// for the calculation of the coordinates (x,y,z) of the intersection points
68 /// for Segment-Facet and Segment-Edge intersections.
70 /// Use 3*facet + coordinate as index. -1 indicates that the coordinate is 0.
71 const int TransformedTriangle::DP_INDEX[12] =
80 /// \brief Correspondance edge - corners.
82 /// Gives the two corners associated with each edge
83 /// Use 2*edge + {0, 1} as index
84 const TransformedTriangle::TetraCorner TransformedTriangle::CORNERS_FOR_EDGE[12] =
94 /// \brief Correspondance edge - facets.
96 /// Gives the two facets shared by and edge. Use 2*facet + {0, 1} as index
97 const TransformedTriangle::TetraFacet TransformedTriangle::FACET_FOR_EDGE[12] =
107 /// \brief Correspondance corners - edges.
109 /// Gives edges meeting at a given corner. Use 3*corner + {0,1,2} as index
110 const TransformedTriangle::TetraEdge TransformedTriangle::EDGES_FOR_CORNER[12] =
118 /// \brief Double products to use in halfstrip intersection tests.
120 /// Use 4*(offset_edge) + {0,1,2,3} as index. offset_edge = edge - 3 (so that XY -> 0, YZ -> 1, ZX -> 2)
121 /// Entries with offset 0 and 1 are for the first condition (positive product)
122 /// and those with offset 2 and 3 are for the second condition (negative product).
123 const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_HALFSTRIP_INTERSECTION[12] =
125 C_10, C_01, C_ZH, C_10, // XY
126 C_01, C_XY, C_XH, C_01, // YZ
127 C_XY, C_10, C_YH, C_XY // ZX
130 /// \brief Double products to use in segment-ray test.
132 /// 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)
133 /// Entries with offset 0 are for first condition (zero double product) and the rest are for condition 3 (in the same
134 /// order as in the article)
135 const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] =
137 C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
138 C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y
139 C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01 // Z
143 * Calculates the point of intersection between the given edge of the tetrahedron and the
144 * triangle PQR. (Grandy, eq [22])
146 * @pre testSurfaceEdgeIntersection(edge) returns true
147 * @param edge edge of tetrahedron
148 * @param pt array of three doubles in which to store the coordinates of the intersection point
150 void TransformedTriangle::calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const
154 // barycentric interpolation between points A and B
155 // : (x,y,z)* = (1-alpha)*A + alpha*B where
156 // alpha = t_A / (t_A - t_B)
158 const TetraCorner corners[2] =
160 CORNERS_FOR_EDGE[2*edge],
161 CORNERS_FOR_EDGE[2*edge + 1]
165 const double tA = calcStableT(corners[0]);
166 const double tB = calcStableT(corners[1]);
167 const double alpha = tA / (tA - tB);
170 LOG(4, "corner A = " << corners[0] << " corner B = " << corners[1] );
171 LOG(4, "tA = " << tA << " tB = " << tB << " alpha= " << alpha );
172 for(int i = 0; i < 3; ++i)
175 pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] +
176 alpha * COORDS_TET_CORNER[3*corners[1] + i];
178 pt[i] = (1 - alpha) * getCoordinateForTetCorner<corners[0], i>() +
179 alpha * getCoordinateForTetCorner<corners[0], i>();
182 assert(pt[i] >= 0.0);
183 assert(pt[i] <= 1.0);
188 * Calculates the point of intersection between the given segment of the triangle
189 * and the given facet of the tetrahedron. (Grandy, eq. [23])
191 * @pre testSurfaceEdgeIntersection(seg, facet) returns true
193 * @param seg segment of the triangle
194 * @param facet facet of the tetrahedron
195 * @param pt array of three doubles in which to store the coordinates of the intersection point
197 void TransformedTriangle::calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const
201 for(int i = 0; i < 3; ++i)
203 const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[3*facet + i];
204 const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + i];
205 s -= sign * calcStableC(seg, dp);
210 // calculate coordinates of intersection point
211 for(int i = 0 ; i < 3; ++i)
213 const int dpIdx = DP_INDEX[3*facet + i];
221 const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
222 const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
223 pt[i] = -( sign * calcStableC(seg, dp) ) / s;
225 LOG(4, "SegmentFacetIntPtCalc : pt[" << i << "] = " << pt[i] );
226 LOG(4, "c(" << seg << ", " << dp << ") = " << sign * calcStableC(seg, dp) );
227 assert(pt[i] >= 0.0);
228 assert(pt[i] <= 1.0);
235 * Tests if the given segment of the triangle intersects the given edge of the tetrahedron (Grandy, eq. [20]
236 * If the OPTIMIZE is defined, it does not do the test the double product that should be zero.
237 * @param seg segment of the triangle
238 * @param edge edge of tetrahedron
239 * @return true if the segment intersects the edge
241 bool TransformedTriangle::testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const
244 // check condition that the double products for one of the two
245 // facets adjacent to the edge has a positive product
246 bool isFacetCondVerified = false;
248 for(int i = 0 ; i < 2 ; ++i)
250 facet[i] = FACET_FOR_EDGE[2*edge + i];
252 // find the two c-values -> the two for the other edges of the facet
255 DoubleProduct dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1];
256 DoubleProduct dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
258 if(dp1 == DoubleProduct( edge ))
261 dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1];
263 else if(dp2 == DoubleProduct( edge ))
266 dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
269 const double c1 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1]*calcStableC(seg, dp1);
270 const double c2 = SIGN_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2]*calcStableC(seg, dp2);
272 //isFacetCondVerified = isFacetCondVerified || c1*c2 > 0.0;
275 isFacetCondVerified = true;
279 if(!isFacetCondVerified)
285 return testSegmentIntersectsFacet(seg, facet[0]) || testSegmentIntersectsFacet(seg, facet[1]);
291 * Calculates the point of intersection between the given segment of the triangle
292 * and the given edge of the tetrahedron. (Grandy, eq. [25])
294 * @pre testSegmentEdgeIntersection(seg, edge) returns true
296 * @param seg segment of the triangle
297 * @param edge edge of the tetrahedron
298 * @param pt array of three doubles in which to store the coordinates of the intersection point
300 void TransformedTriangle::calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const
304 // get the two facets associated with the edge
305 static const TetraFacet FACETS_FOR_EDGE[12] =
315 const TetraFacet facets[2] =
317 FACETS_FOR_EDGE[2*edge],
318 FACETS_FOR_EDGE[2*edge + 1]
321 // calculate s for the two edges
323 for(int i = 0; i < 2; ++i)
326 for(int j = 0; j < 3; ++j)
328 const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[3*facets[i] + j];
329 const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[3*facets[i] + j];
330 s[i] += sign * calcStableC(seg, dp);
334 // calculate denominator
335 const double denominator = s[0]*s[0] + s[1]*s[1];
337 // calculate intersection point
338 for(int i = 0; i < 3; ++i)
340 // calculate double product values for the two faces
342 for(int j = 0 ; j < 2; ++j)
344 const int dpIdx = DP_INDEX[3*facets[j] + i];
345 const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
346 const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
347 c[j] = dpIdx < 0.0 ? 0.0 : sign * calcStableC(seg, dp);
350 // pt[i] = (c1*s1 + c2*s2) / (s1^2 + s2^2)
352 pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator;
354 // strange bug with -O2 enabled : assertion fails when we don't have the following
356 //std::cout << "pt[i] = " << pt[i] << std::endl;
357 //assert(pt[i] >= 0.0); // check we are in tetraeder
358 //assert(pt[i] <= 1.0);
365 * Tests if the given segment of the triangle intersects the given corner of the tetrahedron.
366 * (Grandy, eq. [21]). If OPTIMIZE is defined, the double products that should be zero are not verified.
368 * @param seg segment of the triangle
369 * @param corner corner of the tetrahedron
370 * @return true if the segment intersects the corner
372 bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
376 // facets meeting at a given corner
377 static const TetraFacet FACETS_FOR_CORNER[12] =
385 // check segment intersect a facet
386 for(int i = 0 ; i < 3 ; ++i)
388 const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i];
389 if(testSegmentIntersectsFacet(seg, facet))
399 * Tests if the given segment of the triangle intersects the half-strip above the
400 * given edge of the h = 0 plane. (Grandy, eq. [30])
402 * @param seg segment of the triangle
403 * @param edge edge of the h = 0 plane of the tetrahedron (XY, YZ, ZX)
404 * @return true if the upwards ray from the corner intersects the triangle.
406 bool TransformedTriangle::testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edge)
408 // get right index here to avoid "filling out" array
409 const int edgeIndex = static_cast<int>(edge) - 3;
411 // double products used in test
412 // products 1 and 2 for each edge -> first condition in Grandy [30]
413 // products 3 and 4 for each edge -> third condition
414 // NB : some uncertainty whether these last are correct
415 static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12] =
417 C_10, C_01, C_ZH, C_10, // XY
418 C_01, C_XY, C_XH, C_01, // YZ
419 C_XY, C_10, C_YH, C_XY // ZX
422 // facets to use in second condition (S_m)
423 static const TetraFacet FACET_FOR_HALFSTRIP_INTERSECTION[3] =
425 NO_TET_FACET, // XY -> special case : test with plane H = 0
430 const double cVals[4] =
432 calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]),
433 calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]),
434 calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 2]),
435 calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 3])
438 const TetraFacet facet = FACET_FOR_HALFSTRIP_INTERSECTION[edgeIndex];
441 // special case : facet H = 0
442 const bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
443 LOG(4, "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) );
444 LOG(4, "c2 = " << cVals[2] << ", c3 = " << cVals[3] );
446 return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0);
450 * Calculates the point of intersection between the given segment of the triangle
451 * and the halfstrip above the given edge of the tetrahedron. (Grandy, eq. [31])
453 * @pre testSegmentHalfstripIntersection(seg, edge) returns true
455 * @param seg segment of the triangle
456 * @param edge edge of the tetrahedron defining the halfstrip
457 * @param pt array of three doubles in which to store the coordinates of the intersection point
459 void TransformedTriangle::calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const
464 // get right index here to avoid "filling out" array
465 const int edgeIndex = static_cast<int>(edge) - 3;
466 assert(edgeIndex >= 0);
467 assert(edgeIndex < 3);
469 // Barycentric interpolation on the edge
470 // for edge AB : (x,y,z)* = (1-alpha) * A + alpha * B
471 // where alpha = cB / (cB - cA)
473 const double cA = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]);
474 const double cB = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]);
477 const double alpha = cA / (cA - cB);
479 for(int i = 0; i < 3; ++i)
481 const TetraCorner corners[2] =
483 CORNERS_FOR_EDGE[2*edge],
484 CORNERS_FOR_EDGE[2*edge + 1]
487 const double cornerCoords[2] =
489 COORDS_TET_CORNER[3*corners[0] + i],
490 COORDS_TET_CORNER[3*corners[1] + i]
493 pt[i] = (1 - alpha) * cornerCoords[0] + alpha * cornerCoords[1];
495 assert(pt[i] >= 0.0);
496 assert(pt[i] <= 1.0);
498 assert(epsilonEqualRelative(pt[0] + pt[1] + pt[2], 1.0));
502 * Tests if the given segment of triangle PQR intersects the ray pointing
503 * in the upwards z - direction from the given corner of the tetrahedron. (Grandy eq. [29])
504 * If OPTIMIZE is defined, the double product that should be zero is not verified.
506 * @param seg segment of the triangle PQR
507 * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z)
508 * @return true if the upwards ray from the corner intersects the segment.
510 bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const
512 assert(corner == X || corner == Y || corner == Z);
513 LOG(4, "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner );
515 // readjust index since O is not used
516 const int cornerIdx = static_cast<int>(corner) - 1;
519 //? not sure this is correct
520 static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] =
528 const bool cond21 = testSegmentIntersectsFacet(seg, FIRST_FACET_SEGMENT_RAY_INTERSECTION[cornerIdx]);
529 const bool cond22 = (corner == Z) ? testSegmentIntersectsFacet(seg, OYZ) : testSegmentIntersectsHPlane(seg);
531 if(!(cond21 || cond22))
533 LOG(4, "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 );
538 const double cVals[6] =
540 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 1]),
541 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 2]),
542 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 3]),
543 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 4]),
544 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 5]),
545 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 6]),
549 if(( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) >= 0.0)
551 LOG(4, "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] );
553 return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0;
557 // /////////////////////////////////////////////////////////////////////////////////
558 // Utility methods used in intersection tests ///////////////
559 // /////////////////////////////////////////////////////////////////////////////////
561 * Tests if the triangle PQR surrounds the axis on which the
562 * given edge of the tetrahedron lies.
564 * @param edge edge of tetrahedron
565 * @return true if PQR surrounds edge, false if not (see Grandy, eq. [53])
567 bool TransformedTriangle::testTriangleSurroundsEdge(const TetraEdge edge) const
569 // NB DoubleProduct enum corresponds to TetraEdge enum according to Grandy, table III
570 // so we can use the edge directly
572 const double cPQ = calcStableC(PQ, DoubleProduct(edge));
573 const double cQR = calcStableC(QR, DoubleProduct(edge));
574 const double cRP = calcStableC(RP, DoubleProduct(edge));
576 LOG(5, "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" );
578 // if two or more c-values are zero we disallow x-edge intersection
580 const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
584 LOG(5, "TriangleSurroundsEdge test fails due to too many 0 dp" );
587 return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2;