#include "VectorUtils.hxx"
-
-/**
- * Class representing a circular order of a set of points around their barycenter.
- * It is used with the STL sort() algorithm to sort the point of the two polygons
- *
- */
-class ProjectedCentralCircularSortOrder
+namespace INTERP_UTILS
{
-public:
-
- /// Enumeration of different planes to project on when calculating order
- enum CoordType { XY, XZ, YZ };
-
- /**
- * Constructor
- *
- * @param barycenter double[3] containing the barycenter of the points to be compared
- * @param type plane to project on when comparing. The comparison will not work if all the points are in a plane perpendicular
- * to the plane being projected on
- */
- ProjectedCentralCircularSortOrder(const double* barycenter, const CoordType type)
- : _aIdx((type == YZ) ? 2 : 0),
- _bIdx((type == XY) ? 1 : 2),
- _a(barycenter[_aIdx]),
- _b(barycenter[_bIdx])
- {
- }
/**
- * Comparison operator.
- * Compares the relative position between two points in their ordering around the barycenter.
+ * Class representing a circular order of a set of points around their barycenter.
+ * It is used with the STL sort() algorithm to sort the point of the two polygons
*
- * @param pt1 a double[3] representing a point
- * @param pt2 a double[3] representing a point
- * @return true if the angle of the difference vector between pt1 and the barycenter is greater than that
- * of the difference vector between pt2 and the barycenter.
*/
- bool operator()(const double* pt1, const double* pt2)
+ class ProjectedCentralCircularSortOrder
{
- // calculate angles with the axis
- const double ang1 = atan2(pt1[_aIdx] - _a, pt1[_bIdx] - _b);
- const double ang2 = atan2(pt2[_aIdx] - _a, pt2[_bIdx] - _b);
+ public:
- return ang1 > ang2;
- }
-
-private:
- /// index corresponding to first coordinate of plane on which points are projected
- const int _aIdx;
+ /// Enumeration of different planes to project on when calculating order
+ enum CoordType { XY, XZ, YZ };
- /// index corresponding to second coordinate of plane on which points are projected
- const int _bIdx;
+ /**
+ * Constructor
+ *
+ * @param barycenter double[3] containing the barycenter of the points to be compared
+ * @param type plane to project on when comparing. The comparison will not work if all the points are in a plane perpendicular
+ * to the plane being projected on
+ */
+ ProjectedCentralCircularSortOrder(const double* barycenter, const CoordType type)
+ : _aIdx((type == YZ) ? 2 : 0),
+ _bIdx((type == XY) ? 1 : 2),
+ _a(barycenter[_aIdx]),
+ _b(barycenter[_bIdx])
+ {
+ }
- /// value of first projected coordinate of the barycenter
- const double _a;
-
- /// value of second projected coordinate of the barycenter
- const double _b;
-};
+ /**
+ * Comparison operator.
+ * Compares the relative position between two points in their ordering around the barycenter.
+ *
+ * @param pt1 a double[3] representing a point
+ * @param pt2 a double[3] representing a point
+ * @return true if the angle of the difference vector between pt1 and the barycenter is greater than that
+ * of the difference vector between pt2 and the barycenter.
+ */
+ bool operator()(const double* pt1, const double* pt2)
+ {
+ // calculate angles with the axis
+ const double ang1 = atan2(pt1[_aIdx] - _a, pt1[_bIdx] - _b);
+ const double ang2 = atan2(pt2[_aIdx] - _a, pt2[_bIdx] - _b);
+ return ang1 > ang2;
+ }
-//class Vector3Cmp
-//{
-// public:
-// bool operator()(double* const& pt1, double* const& pt2)
-// {
-// LOG(6, "points are equal ? : " << int((pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2])));
-// return (pt1[0] == pt2[0]) && (pt1[1] == pt2[1]) && (pt1[2] == pt2[2]);
-// }
-// };
+ private:
+ /// index corresponding to first coordinate of plane on which points are projected
+ const int _aIdx;
+
+ /// index corresponding to second coordinate of plane on which points are projected
+ const int _bIdx;
-namespace INTERP_UTILS
-{
- ////////////////////////////////////////////////////////////////////////////
- /// PUBLIC ////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////
+ /// value of first projected coordinate of the barycenter
+ const double _a;
+
+ /// value of second projected coordinate of the barycenter
+ const double _b;
+ };
+ // ----------------------------------------------------------------------------------
+ // TransformedTriangle PUBLIC
+ // ----------------------------------------------------------------------------------
+
/**
* Constructor
*
double volB = 0.0;
// if triangle is not in h = 0 plane, calculate volume under B
- if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ))
- {
- LOG(2, "---- Treating polygon B ... ");
+ if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ))
+ {
+ LOG(2, "---- Treating polygon B ... ");
- calculatePolygonBarycenter(B, barycenter);
- sortIntersectionPolygon(B, barycenter);
- volB = calculateVolumeUnderPolygon(B, barycenter);
- LOG(2, "Volume is " << sign * volB);
- }
+ calculatePolygonBarycenter(B, barycenter);
+ sortIntersectionPolygon(B, barycenter);
+ volB = calculateVolumeUnderPolygon(B, barycenter);
+ LOG(2, "Volume is " << sign * volB);
+ }
LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********");
return sign * (volA + volB);
}
-
- // //////////////////////////////////////////////////////////////////////////
- // PRIVATE /////////////////////////////////////////////////////////////////
- // //////////////////////////////////////////////////////////////////////////
-
- // //////////////////////////////////////////////////////////////////////////////////
- // High-level methods called directly by calculateIntersectionVolume()
- // //////////////////////////////////////////////////////////////////////////////////
+
+ // ----------------------------------------------------------------------------------
+ // TransformedTriangle PRIVATE
+ // ----------------------------------------------------------------------------------
/**
* Calculates the intersection polygons A and B, performing the intersection tests
#else
- // -- segment intersections
- for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
- {
-
- // segment - facet
- for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
+ // -- segment intersections
+ for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
{
- if(testSegmentFacetIntersection(seg, facet))
+
+ // segment - facet
+ for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
{
- double* ptA = new double[3];
- calcIntersectionPtSegmentFacet(seg, facet, ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
- if(facet == XYZ)
+ if(testSegmentFacetIntersection(seg, facet))
{
- double* ptB = new double[3];
- copyVector3(ptA, ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
- }
+ double* ptA = new double[3];
+ calcIntersectionPtSegmentFacet(seg, facet, ptA);
+ _polygonA.push_back(ptA);
+ LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
+ if(facet == XYZ)
+ {
+ double* ptB = new double[3];
+ copyVector3(ptA, ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
+ }
+ }
}
- }
- // segment - edge
- for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
- {
- if(testSegmentEdgeIntersection(seg, edge))
+ // segment - edge
+ for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
{
- double* ptA = new double[3];
- calcIntersectionPtSegmentEdge(seg, edge, ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
- if(edge >= XY)
+ if(testSegmentEdgeIntersection(seg, edge))
{
- double* ptB = new double[3];
- copyVector3(ptA, ptB);
- _polygonB.push_back(ptB);
+ double* ptA = new double[3];
+ calcIntersectionPtSegmentEdge(seg, edge, ptA);
+ _polygonA.push_back(ptA);
+ LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
+ if(edge >= XY)
+ {
+ double* ptB = new double[3];
+ copyVector3(ptA, ptB);
+ _polygonB.push_back(ptB);
+ }
}
}
- }
- // segment - corner
- for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
- {
- if(testSegmentCornerIntersection(seg, corner))
+ // segment - corner
+ for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
{
- double* ptA = new double[3];
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
- if(corner != O)
+ if(testSegmentCornerIntersection(seg, corner))
{
- double* ptB = new double[3];
- _polygonB.push_back(ptB);
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
- LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
+ double* ptA = new double[3];
+ copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
+ _polygonA.push_back(ptA);
+ LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
+ if(corner != O)
+ {
+ double* ptB = new double[3];
+ _polygonB.push_back(ptB);
+ copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
+ LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
+ }
}
}
- }
#endif
#ifdef OPTIMIZE
- // segment - ray
- for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
- {
- if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
+ // segment - ray
+ for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
{
- double* ptB = new double[3];
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
+ if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
+ {
+ double* ptB = new double[3];
+ copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
+ }
}
- }
- // segment - halfstrip
- for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
- {
+ // segment - halfstrip
+ for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
+ {
#if 0
- const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
- const bool doTest =
- !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
- !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
+ const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
+ const bool doTest =
+ !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
+ !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
- if(doTest && testSegmentHalfstripIntersection(seg, edge))
+ if(doTest && testSegmentHalfstripIntersection(seg, edge))
#endif
- if(testSegmentHalfstripIntersection(seg, edge))
- {
- double* ptB = new double[3];
- calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
+ if(testSegmentHalfstripIntersection(seg, edge))
+ {
+ double* ptB = new double[3];
+ calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
+ }
}
- }
#else
- // segment - ray
- for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
- {
- if(testSegmentRayIntersection(seg, corner))
+ // segment - ray
+ for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
{
- double* ptB = new double[3];
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
+ if(testSegmentRayIntersection(seg, corner))
+ {
+ double* ptB = new double[3];
+ copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
+ }
}
- }
- // segment - halfstrip
- for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
- {
- if(testSegmentHalfstripIntersection(seg, edge))
+ // segment - halfstrip
+ for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
{
- double* ptB = new double[3];
- calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
+ if(testSegmentHalfstripIntersection(seg, edge))
+ {
+ double* ptB = new double[3];
+ calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
+ }
}
- }
#endif
- }
+ }
- // inclusion tests
- for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
- {
- // { XYZ - inclusion only possible if in Tetrahedron?
- // tetrahedron
- if(testCornerInTetrahedron(corner))
+ // inclusion tests
+ for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
{
- double* ptA = new double[3];
- copyVector3(&_coords[5*corner], ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
- }
+ // { XYZ - inclusion only possible if in Tetrahedron?
+ // tetrahedron
+ if(testCornerInTetrahedron(corner))
+ {
+ double* ptA = new double[3];
+ copyVector3(&_coords[5*corner], ptA);
+ _polygonA.push_back(ptA);
+ LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
+ }
- // XYZ - plane
- if(testCornerOnXYZFacet(corner))
- {
- double* ptB = new double[3];
- copyVector3(&_coords[5*corner], ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
- }
+ // XYZ - plane
+ if(testCornerOnXYZFacet(corner))
+ {
+ double* ptB = new double[3];
+ copyVector3(&_coords[5*corner], ptB);
+ _polygonB.push_back(ptB);
+ LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
+ }
+
+ // projection on XYZ - facet
+ if(testCornerAboveXYZFacet(corner))
+ {
+ double* ptB = new double[3];
+ copyVector3(&_coords[5*corner], ptB);
+ ptB[2] = 1 - ptB[0] - ptB[1];
+ assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
+ _polygonB.push_back(ptB);
+ LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
+ }
- // projection on XYZ - facet
- if(testCornerAboveXYZFacet(corner))
- {
- double* ptB = new double[3];
- copyVector3(&_coords[5*corner], ptB);
- ptB[2] = 1 - ptB[0] - ptB[1];
- assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
- _polygonB.push_back(ptB);
- LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
}
}
- }
-
- /**
- * Calculates the barycenters of the given intersection polygon.
- *
- * @pre the intersection polygons have been calculated with calculateIntersectionPolygons()
- *
- * @param poly one of the two intersection polygons
- * @param barycenter array of three doubles where barycenter is stored
- *
- */
- void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
- {
- LOG(3,"--- Calculating polygon barycenter");
+ /**
+ * Calculates the barycenters of the given intersection polygon.
+ *
+ * @pre the intersection polygons have been calculated with calculateIntersectionPolygons()
+ *
+ * @param poly one of the two intersection polygons
+ * @param barycenter array of three doubles where barycenter is stored
+ *
+ */
+ void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
+ {
+ LOG(3,"--- Calculating polygon barycenter");
- // get the polygon points
- std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
+ // get the polygon points
+ std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
- // calculate barycenter
- const int m = polygon.size();
+ // calculate barycenter
+ const int m = polygon.size();
- for(int j = 0 ; j < 3 ; ++j)
- {
- barycenter[j] = 0.0;
- }
+ for(int j = 0 ; j < 3 ; ++j)
+ {
+ barycenter[j] = 0.0;
+ }
- if(m != 0)
- {
- for(int i = 0 ; i < m ; ++i)
+ if(m != 0)
{
- const double* pt = polygon[i];
- for(int j = 0 ; j < 3 ; ++j)
+ for(int i = 0 ; i < m ; ++i)
{
- barycenter[j] += pt[j] / double(m);
+ const double* pt = polygon[i];
+ for(int j = 0 ; j < 3 ; ++j)
+ {
+ barycenter[j] += pt[j] / double(m);
+ }
}
}
+ LOG(3,"Barycenter is " << vToStr(barycenter));
}
- LOG(3,"Barycenter is " << vToStr(barycenter));
- }
- /**
- * Sorts the given intersection polygon in circular order around its barycenter.
- * @pre the intersection polygons have been calculated with calculateIntersectionPolygons()
- * @post the vertices in _polygonA and _polygonB are sorted in circular order around their
- * respective barycenters
- *
- * @param poly one of the two intersection polygons
- * @param barycenter array of three doubles with the coordinates of the barycenter
- *
- */
- void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter)
- {
- LOG(3,"--- Sorting polygon ...");
-
- using ::ProjectedCentralCircularSortOrder;
- typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor
- typedef SortOrder::CoordType CoordType;
+ /**
+ * Sorts the given intersection polygon in circular order around its barycenter.
+ * @pre the intersection polygons have been calculated with calculateIntersectionPolygons()
+ * @post the vertices in _polygonA and _polygonB are sorted in circular order around their
+ * respective barycenters
+ *
+ * @param poly one of the two intersection polygons
+ * @param barycenter array of three doubles with the coordinates of the barycenter
+ *
+ */
+ void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter)
+ {
+ LOG(3,"--- Sorting polygon ...");
- // get the polygon points
- std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
+ using ::ProjectedCentralCircularSortOrder;
+ typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor
+ typedef SortOrder::CoordType CoordType;
- if(polygon.size() == 0)
- return;
+ // get the polygon points
+ std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
- // determine type of sorting
- CoordType type = SortOrder::XY;
- if(poly == A) // B is on h = 0 plane -> ok
- {
- // NB : the following test is never true if we have eliminated the
- // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume().
- // We keep the test here anyway, to avoid interdependency.
+ if(polygon.size() == 0)
+ return;
- // is triangle parallel to x == 0 ?
- if(isTriangleParallelToFacet(OZX))
+ // determine type of sorting
+ CoordType type = SortOrder::XY;
+ if(poly == A) // B is on h = 0 plane -> ok
{
- type = SortOrder::YZ;
- }
- else if(isTriangleParallelToFacet(OYZ))
- {
- type = SortOrder::XZ;
+ // NB : the following test is never true if we have eliminated the
+ // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume().
+ // We keep the test here anyway, to avoid interdependency.
+
+ // is triangle parallel to x == 0 ?
+ if(isTriangleParallelToFacet(OZX))
+ {
+ type = SortOrder::YZ;
+ }
+ else if(isTriangleParallelToFacet(OYZ))
+ {
+ type = SortOrder::XZ;
+ }
}
- }
- // create order object
- SortOrder order(barycenter, type);
+ // create order object
+ SortOrder order(barycenter, type);
- // sort vector with this object
- // NB : do not change place of first object, with respect to which the order
- // is defined
- sort((polygon.begin()), polygon.end(), order);
+ // sort vector with this object
+ // NB : do not change place of first object, with respect to which the order
+ // is defined
+ sort((polygon.begin()), polygon.end(), order);
- LOG(3,"Sorted polygon is ");
- for(int i = 0 ; i < polygon.size() ; ++i)
- {
- LOG(3,vToStr(polygon[i]));
- }
+ LOG(3,"Sorted polygon is ");
+ for(int i = 0 ; i < polygon.size() ; ++i)
+ {
+ LOG(3,vToStr(polygon[i]));
+ }
- }
+ }
- /**
- * Calculates the volume between the given polygon and the z = 0 plane.
- *
- * @pre the intersection polygones have been calculated with calculateIntersectionPolygons(),
- * and they have been sorted in circular order with sortIntersectionPolygons(void)
- *
- * @param poly one of the two intersection polygons
- * @param barycenter array of three doubles with the coordinates of the barycenter
- * @return the volume between the polygon and the z = 0 plane
- *
- */
- double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter)
- {
- LOG(2,"--- Calculating volume under polygon");
+ /**
+ * Calculates the volume between the given polygon and the z = 0 plane.
+ *
+ * @pre the intersection polygones have been calculated with calculateIntersectionPolygons(),
+ * and they have been sorted in circular order with sortIntersectionPolygons(void)
+ *
+ * @param poly one of the two intersection polygons
+ * @param barycenter array of three doubles with the coordinates of the barycenter
+ * @return the volume between the polygon and the z = 0 plane
+ *
+ */
+ double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter)
+ {
+ LOG(2,"--- Calculating volume under polygon");
- // get the polygon points
- std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
+ // get the polygon points
+ std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
- double vol = 0.0;
- const int m = polygon.size();
+ double vol = 0.0;
+ const int m = polygon.size();
- for(int i = 0 ; i < m ; ++i)
- {
- const double* ptCurr = polygon[i]; // pt "i"
- const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0)
+ for(int i = 0 ; i < m ; ++i)
+ {
+ const double* ptCurr = polygon[i]; // pt "i"
+ const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0)
- const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2];
- const double factor2 =
- ptCurr[0]*(ptNext[1] - barycenter[1])
- + ptNext[0]*(barycenter[1] - ptCurr[1])
- + barycenter[0]*(ptCurr[1] - ptNext[1]);
- vol += (factor1 * factor2) / 6.0;
- }
+ const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2];
+ const double factor2 =
+ ptCurr[0]*(ptNext[1] - barycenter[1])
+ + ptNext[0]*(barycenter[1] - ptCurr[1])
+ + barycenter[0]*(ptCurr[1] - ptNext[1]);
+ vol += (factor1 * factor2) / 6.0;
+ }
- LOG(2,"Abs. Volume is " << vol);
- return vol;
- }
+ LOG(2,"Abs. Volume is " << vol);
+ return vol;
+ }
- ////////////////////////////////////////////////////////////////////////////////////
- // Detection of (very) degenerate cases /////////////
- ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Detection of (very) degenerate cases /////////////
+ ////////////////////////////////////////////////////////////////////////////////////
- /**
- * Checks if the triangle lies in the plane of a given facet
- *
- * @param facet one of the facets of the tetrahedron
- * @return true if PQR lies in the plane of the facet, false if not
- */
- bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const
- {
+ /**
+ * Checks if the triangle lies in the plane of a given facet
+ *
+ * @param facet one of the facets of the tetrahedron
+ * @return true if PQR lies in the plane of the facet, false if not
+ */
+ bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const
+ {
- // coordinate to check
- const int coord = static_cast<int>(facet);
+ // coordinate to check
+ const int coord = static_cast<int>(facet);
- for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
- {
- if(_coords[5*c + coord] != 0.0)
+ for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
{
- return false;
+ if(_coords[5*c + coord] != 0.0)
+ {
+ return false;
+ }
}
- }
- return true;
- }
+ return true;
+ }
- /**
- * Checks if the triangle is parallel to the given facet
- *
- * @param facet one of the facets of the unit tetrahedron
- * @return true if triangle is parallel to facet, false if not
- */
- bool TransformedTriangle::isTriangleParallelToFacet(const TetraFacet facet) const
- {
- // coordinate to check
- const int coord = static_cast<int>(facet);
- return (_coords[5*P + coord] == _coords[5*Q + coord]) && (_coords[5*P + coord] == _coords[5*R + coord]);
- }
+ /**
+ * Checks if the triangle is parallel to the given facet
+ *
+ * @param facet one of the facets of the unit tetrahedron
+ * @return true if triangle is parallel to facet, false if not
+ */
+ bool TransformedTriangle::isTriangleParallelToFacet(const TetraFacet facet) const
+ {
+ // coordinate to check
+ const int coord = static_cast<int>(facet);
+ return (_coords[5*P + coord] == _coords[5*Q + coord]) && (_coords[5*P + coord] == _coords[5*R + coord]);
+ }
- /**
- * Determines whether the triangle is below the z-plane.
- *
- * @return true if the z-coordinate of the three corners of the triangle are all less than 0, false otherwise.
- */
- bool TransformedTriangle::isTriangleBelowTetraeder() const
- {
- for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
+ /**
+ * Determines whether the triangle is below the z-plane.
+ *
+ * @return true if the z-coordinate of the three corners of the triangle are all less than 0, false otherwise.
+ */
+ bool TransformedTriangle::isTriangleBelowTetraeder() const
{
- // check z-coords for all points
- if(_coords[5*c + 2] >= 0.0)
+ for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
{
- return false;
+ // check z-coords for all points
+ if(_coords[5*c + 2] >= 0.0)
+ {
+ return false;
+ }
}
+ return true;
}
- return true;
- }
- /**
- * Prints the coordinates of the triangle to std::cout
- *
- */
- void TransformedTriangle::dumpCoords() const
- {
- std::cout << "Coords : ";
- for(int i = 0 ; i < 3; ++i)
+ /**
+ * Prints the coordinates of the triangle to std::cout
+ *
+ */
+ void TransformedTriangle::dumpCoords() const
{
- std::cout << vToStr(&_coords[5*i]) << ",";
+ std::cout << "Coords : ";
+ for(int i = 0 ; i < 3; ++i)
+ {
+ std::cout << vToStr(&_coords[5*i]) << ",";
+ }
+ std::cout << std::endl;
}
- std::cout << std::endl;
- }
-}; // NAMESPACE
+ }; // NAMESPACE
#include <vector>
+
+
#ifdef OPTIMIZE
+/// OPT_INLINE will be replaced by "inline" if OPTIMIZE is defined and else nothing
#define OPT_INLINE inline
#else
+/// OPT_INLINE will be replaced by "inline" if OPTIMIZE is defined and else nothing
#define OPT_INLINE
#endif
* When an intersection point has been detected it is calculated with a corresponding calc* - method in the cases where it
* is not known directly. It is then added to the polygon A and/or B as necessary.
*
- *
+ * OPTIMIZE :
+ * If OPTIMIZE is defined, a large number of methods will be prefixed with inline and some optimizations concerning the tests
+ * with zero double products will be used.
*/
class TransformedTriangle
{
private:
+
+ // ----------------------------------------------------------------------------------
+ // High-level methods called directly by calculateIntersectionVolume()
+ // ----------------------------------------------------------------------------------
void calculateIntersectionPolygons();
void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter);
double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter);
- ////////////////////////////////////////////////////////////////////////////////////
- /// Detection of degenerate triangles ////////////
- ////////////////////////////////////////////////////////////////////////////////////
+ // ----------------------------------------------------------------------------------
+ // Detection of degenerate triangles
+ // ----------------------------------------------------------------------------------
bool isTriangleInPlaneOfFacet(const TetraFacet facet) const;
bool isTriangleBelowTetraeder() const;
- ////////////////////////////////////////////////////////////////////////////////////
- /// Intersection test methods and intersection point calculations ////////
- ////////////////////////////////////////////////////////////////////////////////////
+ // ----------------------------------------------------------------------------------
+ // Intersection test methods and intersection point calculations
+ // ----------------------------------------------------------------------------------
OPT_INLINE bool testSurfaceEdgeIntersection(const TetraEdge edge) const;
OPT_INLINE bool testCornerAboveXYZFacet(const TriCorner corner) const;
- ////////////////////////////////////////////////////////////////////////////////////
- /// Utility methods used in intersection tests ///////////////
- ////////////////////////////////////////////////////////////////////////////////////
+ // ----------------------------------------------------------------------------------
+ // Utility methods used in intersection tests
+ // ----------------------------------------------------------------------------------
bool testTriangleSurroundsEdge(const TetraEdge edge) const;
bool testTriangleSurroundsRay(const TetraCorner corner) const;
- ////////////////////////////////////////////////////////////////////////////////////
- /// Double and triple product calculations ///////////////
- ////////////////////////////////////////////////////////////////////////////////////
+ // ----------------------------------------------------------------------------------
+ // Double and triple product calculations
+ // ----------------------------------------------------------------------------------
- void preCalculateDoubleProducts(void);
+
bool areDoubleProductsConsistent(const TriSegment seg) const;
+ void preCalculateDoubleProducts(void);
+
OPT_INLINE void resetDoubleProducts(const TriSegment seg, const TetraCorner corner);
double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const;
double calcTByDevelopingRow(const TetraCorner corner, const int row = 1, const bool project = false) const;
- ////////////////////////////////////////////////////////////////////////////////////
- /// Member variables ///////////////
- ////////////////////////////////////////////////////////////////////////////////////
+ // ----------------------------------------------------------------------------------
+ // Member variables
+ // ----------------------------------------------------------------------------------
private:
- // order :
- // [ p_x, p_y, p_z, p_h, p_H, q_x, q_y, q_z, q_h, q_H, r_x, r_y, r_z, r_h, r_H ]
+ /// Array holding the coordinates of the triangle's three corners
+ /// order :
+ /// [ p_x, p_y, p_z, p_h, p_H, q_x, q_y, q_z, q_h, q_H, r_x, r_y, r_z, r_h, r_H ]
double _coords[15];
- /// flags showing whether the double and triple products have been precalculated for this class
- bool _isDoubleProductsCalculated, _isTripleProductsCalculated;
+ /// Flag showing whether the double products have been calculated yet
+ bool _isDoubleProductsCalculated;
+
+ /// Flag showing whether the triple products have been calculated yet
+ bool _isTripleProductsCalculated;
- /// array containing the 24 double products
+ /// Array containing the 24 double products.
/// order : c^PQ_YZ, ... ,cPQ_10, ... c^QR_YZ, ... c^RP_YZ
/// following order in enumeration DoubleProduct
double _doubleProducts[24];
- /// array containing the 4 triple products
+ /// Array containing the 4 triple products.
/// order : t_O, t_X, t_Y, t_Z
double _tripleProducts[4];
- /// arrays holding the points in the two intersection polygons A and B
+ /// Vector holding the points of the intersection polygon A.
/// these points are allocated in calculateIntersectionPolygons() and liberated in the destructor
- std::vector<double*> _polygonA, _polygonB;
+ std::vector<double*> _polygonA;
- /// vectors holding the coordinates of the barycenters of the polygons A and B
- /// these points are calculated in calculatePolygonBarycenter
- double _barycenterA[3], _barycenterB[3];
+ /// Vector holding the points of the intersection polygon B.
+ /// These points are allocated in calculateIntersectionPolygons() and liberated in the destructor
+ std::vector<double*> _polygonB;
+
+ /// Array holding the coordinates of the barycenter of the polygon A
+ /// This point is calculated in calculatePolygonBarycenter
+ double _barycenterA[3];
+
+ /// Array holding the coordinates of the barycenter of the polygon B
+ /// This point is calculated in calculatePolygonBarycenter
+ double _barycenterB[3];
- // used for debugging
+ /// Array of flags indicating which of the four triple products have been correctly calculated.
+ /// Used for asserts in debug mode
bool _validTP[4];
#ifdef OPTIMIZE
- void preCalculateTriangleSurroundsEdge();
+ void preCalculateTriangleSurroundsEdge();
+
+ /// Array holding results of the test testTriangleSurroundsEdge() for all the edges.
+ /// These are calculated in preCalculateTriangleSurroundsEdge().
bool _triangleSurroundsEdgeCache[NO_TET_EDGE];
- bool _isOutsideTetra;
#endif
- ////////////////////////////////////////////////////////////////////////////////////
- /// Constants /////////////////
- ////////////////////////////////////////////////////////////////////////////////////
+ // ----------------------------------------------------------------------------------
+ // Constants
+ // ----------------------------------------------------------------------------------
// offsets : 0 -> x, 1 -> y, 2 -> z, 3 -> h, 4 -> H
// corresponds to order of double products in DoubleProduct
// double products used in segment - ray test
static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21];
-
- inline bool isTriangleOutsideTetra(void) const;
-
};
#include <cmath>
#include "VectorUtils.hxx"
-#define SEG_RAY_TABLE 1 // seems correct
-
namespace INTERP_UTILS
{
- ////////////////////////////////////////////////////////////////////////////////////
- /// Correspondance tables describing all the variations of formulas. //////////////
- ////////////////////////////////////////////////////////////////////////////////////
+ // ----------------------------------------------------------------------------------
+ // Correspondance tables describing all the variations of formulas.
+ // ----------------------------------------------------------------------------------
- // correspondance facet - double product
- // Grandy, table IV
+ /// Correspondance between facets and double products.
+ /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_SEG_FACET_INTERSECTION[12] =
{
C_XH, C_XY, C_ZX, // OYZ
C_XH, C_YH, C_ZH // XYZ
};
- // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION
+ /// Signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION
+ /// This table encodes Grandy, table IV. Use 3*facet + {0,1,2} as index
const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] =
{
1.0, 1.0, -1.0,
1.0, 1.0, 1.0
};
- // coordinates of corners of tetrahedron
+ /// Coordinates of corners of tetrahedron.
+ /// Use 3*Corner + coordinate as index
const double TransformedTriangle::COORDS_TET_CORNER[12] =
{
0.0, 0.0, 0.0,
0.0, 0.0, 1.0
};
- // indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
- // for the calculation of the coordinates (x,y,z) of the intersection points
- // for Segment-Facet and Segment-Edge intersections
- // -1 indicates that the coordinate is 0
+ /// Indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
+ /// for the calculation of the coordinates (x,y,z) of the intersection points
+ /// for Segment-Facet and Segment-Edge intersections.
+ /// Use 3*facet + coordinate as index. -1 indicates that the coordinate is 0.
const int TransformedTriangle::DP_INDEX[12] =
{
// x, y, z
9, 10, 11 // XYZ
};
- // correspondance edge - corners
+ /// Correspondance edge - corners
+ /// Gives the two corners associated with each edge
+ /// Use 2*edge + {0, 1} as index
const TransformedTriangle::TetraCorner TransformedTriangle::CORNERS_FOR_EDGE[12] =
{
O, X, // OX
Z, X // ZX
};
- // correspondance edge - facets
- // facets shared by each edge
+ /// Correspondance edge - facets.
+ /// Gives the two facets shared by and edge. Use 2*facet + {0, 1} as index
const TransformedTriangle::TetraFacet TransformedTriangle::FACET_FOR_EDGE[12] =
{
OXY, OZX, // OX
OZX, XYZ // ZX
};
- // edges meeting at a given corner
+ /// Correspondance corners - edges
+ /// Gives edges meeting at a given corner. Use 3*corner + {0,1,2} as index
const TransformedTriangle::TetraEdge TransformedTriangle::EDGES_FOR_CORNER[12] =
{
OX, OY, OZ, // O
OZ, ZX, YZ // Z
};
- // NB : some uncertainty whether these last are correct
+ /// Double products to use in halfstrip intersection tests
+ /// Use 4*(offset_edge) + {0,1,2,3} as index. offset_edge = edge - 3 (so that XY -> 0, YZ -> 1, ZX -> 2)
+ /// Entries with offset 0 and 1 are for the first condition (positive product)
+ /// and those with offset 2 and 3 are for the second condition (negative product).
const TransformedTriangle::DoubleProduct TransformedTriangle::DP_FOR_HALFSTRIP_INTERSECTION[12] =
{
C_10, C_01, C_ZH, C_10, // XY
C_XY, C_10, C_YH, C_XY // ZX
};
- // double products to use in segment-ray test
- // dp 1 -> cond 1
- // dp 2-7 -> cond 3
-#if SEG_RAY_TABLE==1
+ /// Double products to use in segment-ray test
+ /// 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)
+ /// Entries with offset 0 are for first condition (zero double product) and the rest are for condition 3 (in the same
+ /// order as in the article)
const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] =
{
C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_10, // Y
C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_01 // Z
};
-#else
-
- const TransformedTriangle::DoubleProduct TransformedTriangle::DP_SEGMENT_RAY_INTERSECTION[21] =
- {
- C_10, C_YH, C_ZH, C_01, C_XY, C_YH, C_XY, // X
- C_01, C_XH, C_ZH, C_XY, C_10, C_ZH, C_YZ, // Y
- C_XY, C_YH, C_XH, C_10, C_01, C_XH, C_ZX // Z
- };
-#endif
-
- ////////////////////////////////////////////////////////////////////////////////////
- /// Intersection test methods and intersection point calculations ////////
- ////////////////////////////////////////////////////////////////////////////////////
+ // ----------------------------------------------------------------------------------
+ // Intersection test methods and intersection point calculations
+ // ----------------------------------------------------------------------------------
#ifndef OPTIMIZE // inlined otherwise -> see TransformedTriangle_inline.hxx
/**
* Tests if the given edge of the tetrahedron intersects the triangle PQR. (Grandy, eq [17])
#endif
- ////////////////////////////////////////////////////////////////////////////////////
- /// Utility methods used in intersection tests ///////////////
- ////////////////////////////////////////////////////////////////////////////////////
+ // /////////////////////////////////////////////////////////////////////////////////
+ // Utility methods used in intersection tests ///////////////
+ // /////////////////////////////////////////////////////////////////////////////////
/**
* Tests if the triangle PQR surrounds the axis on which the
* given edge of the tetrahedron lies.
