1 // Copyright (C) 2007-2019 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, 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 = " << corners[0] << " corner B = " << 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] +
177 alpha * COORDS_TET_CORNER[3*corners[1] + i];
179 pt[i] = (1 - alpha) * getCoordinateForTetCorner<corners[0], i>() +
180 alpha * getCoordinateForTetCorner<corners[0], i>();
183 assert(pt[i] >= 0.0);
184 assert(pt[i] <= 1.0);
189 * Calculates the point of intersection between the given segment of the triangle
190 * and the given facet of the tetrahedron. (Grandy, eq. [23])
192 * @pre testSurfaceEdgeIntersection(seg, facet) returns true
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
198 void TransformedTriangle::calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const
202 for(int i = 0; i < 3; ++i)
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);
211 // calculate coordinates of intersection point
212 for(int i = 0 ; i < 3; ++i)
214 const int dpIdx = DP_INDEX[3*facet + i];
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;
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);
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
242 bool TransformedTriangle::testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const
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;
249 for(int i = 0 ; i < 2 ; ++i)
251 facet[i] = FACET_FOR_EDGE[2*edge + i];
253 // find the two c-values -> the two for the other edges of the facet
256 DoubleProduct dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1];
257 DoubleProduct dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
259 if(dp1 == DoubleProduct( edge ))
262 dp1 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx1];
264 else if(dp2 == DoubleProduct( edge ))
267 dp2 = DP_FOR_SEG_FACET_INTERSECTION[3*facet[i] + idx2];
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);
273 //isFacetCondVerified = isFacetCondVerified || c1*c2 > 0.0;
276 isFacetCondVerified = true;
280 if(!isFacetCondVerified)
286 return testSegmentIntersectsFacet(seg, facet[0]) || testSegmentIntersectsFacet(seg, facet[1]);
292 * Calculates the point of intersection between the given segment of the triangle
293 * and the given edge of the tetrahedron. (Grandy, eq. [25])
295 * @pre testSegmentEdgeIntersection(seg, edge) returns true
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
301 void TransformedTriangle::calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const
305 // get the two facets associated with the edge
306 static const TetraFacet FACETS_FOR_EDGE[12] =
316 const TetraFacet facets[2] =
318 FACETS_FOR_EDGE[2*edge],
319 FACETS_FOR_EDGE[2*edge + 1]
322 // calculate s for the two edges
324 for(int i = 0; i < 2; ++i)
327 for(int j = 0; j < 3; ++j)
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);
335 // calculate denominator
336 const double denominator = s[0]*s[0] + s[1]*s[1];
338 // calculate intersection point
339 for(int i = 0; i < 3; ++i)
341 // calculate double product values for the two faces
343 for(int j = 0 ; j < 2; ++j)
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);
351 // pt[i] = (c1*s1 + c2*s2) / (s1^2 + s2^2)
353 pt[i] = (c[0] * s[0] + c[1] * s[1]) / denominator;
355 // strange bug with -O2 enabled : assertion fails when we don't have the following
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);
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.
369 * @param seg segment of the triangle
370 * @param corner corner of the tetrahedron
371 * @return true if the segment intersects the corner
373 bool TransformedTriangle::testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const
377 // facets meeting at a given corner
378 static const TetraFacet FACETS_FOR_CORNER[12] =
386 // check segment intersect a facet
387 for(int i = 0 ; i < 3 ; ++i)
389 const TetraFacet facet = FACETS_FOR_CORNER[3*corner + i];
390 if(testSegmentIntersectsFacet(seg, facet))
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])
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.
407 bool TransformedTriangle::testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edge)
409 // get right index here to avoid "filling out" array
410 const int edgeIndex = static_cast<int>(edge) - 3;
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
418 // facets to use in second condition (S_m)
419 static const TetraFacet FACET_FOR_HALFSTRIP_INTERSECTION[3] =
421 NO_TET_FACET, // XY -> special case : test with plane H = 0
426 const double cVals[4] =
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])
434 const TetraFacet facet = FACET_FOR_HALFSTRIP_INTERSECTION[edgeIndex];
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] );
442 return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0);
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])
449 * @pre testSegmentHalfstripIntersection(seg, edge) returns true
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
455 void TransformedTriangle::calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const
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);
465 // Barycentric interpolation on the edge
466 // for edge AB : (x,y,z)* = (1-alpha) * A + alpha * B
467 // where alpha = cB / (cB - cA)
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]);
473 const double alpha = cA / (cA - cB);
475 for(int i = 0; i < 3; ++i)
477 const TetraCorner corners[2] =
479 CORNERS_FOR_EDGE[2*edge],
480 CORNERS_FOR_EDGE[2*edge + 1]
483 const double cornerCoords[2] =
485 COORDS_TET_CORNER[3*corners[0] + i],
486 COORDS_TET_CORNER[3*corners[1] + i]
489 pt[i] = (1 - alpha) * cornerCoords[0] + alpha * cornerCoords[1];
491 assert(pt[i] >= 0.0);
492 assert(pt[i] <= 1.0);
494 assert(epsilonEqualRelative(pt[0] + pt[1] + pt[2], 1.0));
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.
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.
506 bool TransformedTriangle::testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const
508 assert(corner == X || corner == Y || corner == Z);
509 LOG(4, "Testing seg - ray intersection for seg = " << seg << ", corner = " << corner );
511 // readjust index since O is not used
512 const int cornerIdx = static_cast<int>(corner) - 1;
515 //? not sure this is correct
516 static const TetraFacet FIRST_FACET_SEGMENT_RAY_INTERSECTION[3] =
524 const bool cond21 = testSegmentIntersectsFacet(seg, FIRST_FACET_SEGMENT_RAY_INTERSECTION[cornerIdx]);
525 const bool cond22 = (corner == Z) ? testSegmentIntersectsFacet(seg, OYZ) : testSegmentIntersectsHPlane(seg);
527 if(!(cond21 || cond22))
529 LOG(4, "SR fails at cond 2 : cond21 = " << cond21 << ", cond22 = " << cond22 );
534 const double cVals[6] =
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]),
545 if(( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) >= 0.0)
547 LOG(4, "SR fails at cond 3 : " << (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] );
549 return ( (cVals[0] + cVals[1])*(cVals[2] - cVals[3]) - cVals[4]*cVals[5] ) < 0.0;
553 // /////////////////////////////////////////////////////////////////////////////////
554 // Utility methods used in intersection tests ///////////////
555 // /////////////////////////////////////////////////////////////////////////////////
557 * Tests if the triangle PQR surrounds the axis on which the
558 * given edge of the tetrahedron lies.
560 * @param edge edge of tetrahedron
561 * @return true if PQR surrounds edge, false if not (see Grandy, eq. [53])
563 bool TransformedTriangle::testTriangleSurroundsEdge(const TetraEdge edge) const
565 // NB DoubleProduct enum corresponds to TetraEdge enum according to Grandy, table III
566 // so we can use the edge directly
568 const double cPQ = calcStableC(PQ, DoubleProduct(edge));
569 const double cQR = calcStableC(QR, DoubleProduct(edge));
570 const double cRP = calcStableC(RP, DoubleProduct(edge));
572 LOG(5, "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" );
574 // if two or more c-values are zero we disallow x-edge intersection
576 const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
580 LOG(5, "TriangleSurroundsEdge test fails due to too many 0 dp" );
583 return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2;