1 // Copyright (C) 2007-2024 CEA, EDF
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.
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
20 #include "TransformedTriangle.hxx"
25 #include "VectorUtils.hxx"
27 namespace INTERP_KERNEL
30 // ----------------------------------------------------------------------------------
31 // Correspondence tables describing all the variations of formulas.
32 // ----------------------------------------------------------------------------------
34 /// \brief Correspondence between facets and double products.
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] =
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
45 /// \brief Signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION.
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] =
56 /// \brief Coordinates of corners of tetrahedron.
58 /// Use 3*Corner + coordinate as index
59 const double TransformedTriangle::COORDS_TET_CORNER[12] =
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.
71 /// Use 3*facet + coordinate as index. -1 indicates that the coordinate is 0.
72 const int TransformedTriangle::DP_INDEX[12] =
81 /// \brief Correspondence edge - corners.
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] =
95 /// \brief Correspondence edge - facets.
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] =
108 /// \brief Correspondence corners - edges.
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] =
119 /// \brief Double products to use in halfstrip intersection tests.
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] =
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
131 /// \brief Double products to use in segment-ray test.
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] =
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
144 * Calculates the point of intersection between the given edge of the tetrahedron and the
145 * triangle PQR. (Grandy, eq [22])
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
151 void TransformedTriangle::calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const
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)
159 const TetraCorner corners[2] =
161 CORNERS_FOR_EDGE[2*edge],
162 CORNERS_FOR_EDGE[2*edge + 1]
166 const double tA = calcStableT(corners[0]);
167 const double tB = calcStableT(corners[1]);
168 const double alpha = tA / (tA - tB);
171 LOG(4, "corner A = " << strTC(corners[0]) << " corner B = " << strTC(corners[1]) );
172 LOG(4, "tA = " << tA << " tB = " << tB << " alpha= " << alpha );
173 for(int i = 0; i < 3; ++i)
176 pt[i] = (1 - alpha)*COORDS_TET_CORNER[3*corners[0] + i] + 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]
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];
347 const DoubleProduct dp = DP_FOR_SEG_FACET_INTERSECTION[dpIdx];
348 const double sign = SIGN_FOR_SEG_FACET_INTERSECTION[dpIdx];
349 c[j] = sign * calcStableC(seg, dp);
355 pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator;
361 * Tests if the given segment of the triangle intersects the given corner of the tetrahedron.
362 * (Grandy, eq. [21]).
364 * @param seg segment of the triangle
365 * @param corner corner of the tetrahedron
366 * @return true if the segment intersects the corner
368 bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
370 // facets meeting at a given corner
371 static const TetraFacet FACETS_FOR_CORNER[12] =
379 // check segment intersect a facet
380 for(int i = 0 ; i < 3 ; ++i)
382 const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i];
383 if(testSegmentIntersectsFacet(seg, facet))
391 * Tests if the given segment of the triangle intersects the half-strip above the
392 * given edge of the h = 0 plane. (Grandy, eq. [30])
394 * @param seg segment of the triangle
395 * @param edge edge of the h = 0 plane of the tetrahedron (XY, YZ, ZX)
396 * @return true if the upwards ray from the corner intersects the triangle.
398 bool TransformedTriangle::testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edge)
400 // get right index here to avoid "filling out" array
401 const int edgeIndex = static_cast<int>(edge) - 3;
403 // double products used in test
404 // products 1 and 2 for each edge -> first condition in Grandy [30]
405 // products 3 and 4 for each edge -> third condition
406 // NB : some uncertainty whether these last are correct
407 // DP_FOR_HALFSTRIP_INTERSECTION
409 // facets to use in second condition (S_m)
410 static const TetraFacet FACET_FOR_HALFSTRIP_INTERSECTION[3] =
412 NO_TET_FACET, // XY -> special case : test with plane H = 0
417 const double cVals[4] =
419 calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]),
420 calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]),
421 calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 2]),
422 calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 3])
425 const TetraFacet facet = FACET_FOR_HALFSTRIP_INTERSECTION[edgeIndex];
427 // special case : facet H = 0
428 const bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
429 LOG(4, "Halfstrip tests (" << strTriS(seg) << ", " << strTE(edge) << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) );
430 LOG(4, "c2 = " << cVals[2] << ", c3 = " << cVals[3] );
432 return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0);
436 * Calculates the point of intersection between the given segment of the triangle
437 * and the halfstrip above the given edge of the tetrahedron. (Grandy, eq. [31])
439 * @pre testSegmentHalfstripIntersection(seg, edge) returns true
441 * @param seg segment of the triangle
442 * @param edge edge of the tetrahedron defining the halfstrip
443 * @param pt array of three doubles in which to store the coordinates of the intersection point
445 void TransformedTriangle::calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const
450 // get right index here to avoid "filling out" array
451 const int edgeIndex = static_cast<int>(edge) - 3;
452 assert(edgeIndex >= 0);
453 assert(edgeIndex < 3);
455 // Barycentric interpolation on the edge
456 // for edge AB : (x,y,z)* = (1-alpha) * A + alpha * B
457 // where alpha = cB / (cB - cA)
459 const double cA = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex]);
460 const double cB = calcStableC(seg, DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIndex + 1]);
463 const double alpha = cA / (cA - cB);
465 for(int i = 0; i < 3; ++i)
467 const TetraCorner corners[2] =
469 CORNERS_FOR_EDGE[2*edge],
470 CORNERS_FOR_EDGE[2*edge + 1]
473 const double cornerCoords[2] =
475 COORDS_TET_CORNER[3*corners[0] + i],
476 COORDS_TET_CORNER[3*corners[1] + i]
479 pt[i] = (1 - alpha) * cornerCoords[0] + alpha * cornerCoords[1];
481 assert(pt[i] >= 0.0);
482 assert(pt[i] <= 1.0);
484 assert(epsilonEqualRelative(pt[0] + pt[1] + pt[2], 1.0));
488 * Tests if the given segment of triangle PQR intersects the ray pointing
489 * in the upwards z - direction from the given corner of the tetrahedron. (Grandy eq. [29])
491 * @param seg segment of the triangle PQR
492 * @param corner corner of the tetrahedron on the h = 0 facet (X, Y, or Z)
493 * @return true if the upwards ray from the corner intersects the segment.
495 bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const
497 assert(corner == X || corner == Y || corner == Z);
498 LOG(4, "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner );
500 // readjust index since O is not used
501 const int cornerIdx = static_cast<int>(corner) - 1;
504 //? not sure this is correct
505 static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] =
513 const bool cond21 = testSegmentIntersectsFacet(seg, FIRST_FACET_SEGMENT_RAY_INTERSECTION[cornerIdx]);
514 const bool cond22 = (corner == Z) ? testSegmentIntersectsFacet(seg, OYZ) : testSegmentIntersectsHPlane(seg);
516 if(!(cond21 || cond22))
518 LOG(4, "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 );
523 const double cVals[6] =
525 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 1]),
526 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 2]),
527 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 3]),
528 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 4]),
529 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 5]),
530 calcStableC(seg, DP_SEGMENT_RAY_INTERSECTION[7*cornerIdx + 6]),
534 if(( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) >= 0.0)
536 LOG(4, "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] );
538 return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0;
542 // /////////////////////////////////////////////////////////////////////////////////
543 // Utility methods used in intersection tests ///////////////
544 // /////////////////////////////////////////////////////////////////////////////////
546 * Tests if the triangle PQR surrounds the axis on which the
547 * given edge of the tetrahedron lies.
549 * @param edge edge of tetrahedron
550 * @return true if PQR surrounds edge, false if not (see Grandy, eq. [53])
552 bool TransformedTriangle::testTriangleSurroundsEdge(const TetraEdge edge) const
554 // NB DoubleProduct enum corresponds to TetraEdge enum according to Grandy, table III
555 // so we can use the edge directly
557 const double cPQ = calcStableC(PQ, DoubleProduct(edge));
558 const double cQR = calcStableC(QR, DoubleProduct(edge));
559 const double cRP = calcStableC(RP, DoubleProduct(edge));
561 LOG(5, "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" );
563 // if two or more c-values are zero we disallow x-edge intersection
565 // test for == 0.0 is OK here, if you look at the way double products were pre-comuted:
566 const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
570 LOG(5, "TriangleSurroundsEdge test fails due to too many 0 dp" );
573 return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2;