summaryrefslogtreecommitdiff
path: root/src/mesh/assimp-master/code/AssetLib/IFC/IFCBoolean.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/mesh/assimp-master/code/AssetLib/IFC/IFCBoolean.cpp')
-rw-r--r--src/mesh/assimp-master/code/AssetLib/IFC/IFCBoolean.cpp765
1 files changed, 765 insertions, 0 deletions
diff --git a/src/mesh/assimp-master/code/AssetLib/IFC/IFCBoolean.cpp b/src/mesh/assimp-master/code/AssetLib/IFC/IFCBoolean.cpp
new file mode 100644
index 0000000..36912a7
--- /dev/null
+++ b/src/mesh/assimp-master/code/AssetLib/IFC/IFCBoolean.cpp
@@ -0,0 +1,765 @@
+/*
+Open Asset Import Library (assimp)
+----------------------------------------------------------------------
+
+Copyright (c) 2006-2022, assimp team
+All rights reserved.
+
+Redistribution and use of this software in source and binary forms,
+with or without modification, are permitted provided that the
+following conditions are met:
+
+* Redistributions of source code must retain the above
+ copyright notice, this list of conditions and the
+ following disclaimer.
+
+* Redistributions in binary form must reproduce the above
+ copyright notice, this list of conditions and the
+ following disclaimer in the documentation and/or other
+ materials provided with the distribution.
+
+* Neither the name of the assimp team, nor the names of its
+ contributors may be used to endorse or promote products
+ derived from this software without specific prior
+ written permission of the assimp team.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+----------------------------------------------------------------------
+*/
+
+/** @file IFCBoolean.cpp
+ * @brief Implements a subset of Ifc boolean operations
+ */
+
+#ifndef ASSIMP_BUILD_NO_IFC_IMPORTER
+
+#include "AssetLib/IFC/IFCUtil.h"
+#include "Common/PolyTools.h"
+#include "PostProcessing/ProcessHelper.h"
+
+
+#include <iterator>
+#include <tuple>
+
+namespace Assimp {
+namespace IFC {
+
+// ------------------------------------------------------------------------------------------------
+// Calculates intersection between line segment and plane. To catch corner cases, specify which side you prefer.
+// The function then generates a hit only if the end is beyond a certain margin in that direction, filtering out
+// "very close to plane" ghost hits as long as start and end stay directly on or within the given plane side.
+bool IntersectSegmentPlane(const IfcVector3 &p, const IfcVector3 &n, const IfcVector3 &e0,
+ const IfcVector3 &e1, bool assumeStartOnWhiteSide, IfcVector3 &out) {
+ const IfcVector3 pdelta = e0 - p, seg = e1 - e0;
+ const IfcFloat dotOne = n * seg, dotTwo = -(n * pdelta);
+
+ // if segment ends on plane, do not report a hit. We stay on that side until a following segment starting at this
+ // point leaves the plane through the other side
+ if (std::abs(dotOne + dotTwo) < ai_epsilon)
+ return false;
+
+ // if segment starts on the plane, report a hit only if the end lies on the *other* side
+ if (std::abs(dotTwo) < ai_epsilon) {
+ if ((assumeStartOnWhiteSide && dotOne + dotTwo < ai_epsilon) || (!assumeStartOnWhiteSide && dotOne + dotTwo > -ai_epsilon)) {
+ out = e0;
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ // ignore if segment is parallel to plane and far away from it on either side
+ // Warning: if there's a few thousand of such segments which slowly accumulate beyond the epsilon, no hit would be registered
+ if (std::abs(dotOne) < ai_epsilon)
+ return false;
+
+ // t must be in [0..1] if the intersection point is within the given segment
+ const IfcFloat t = dotTwo / dotOne;
+ if (t > 1.0 || t < 0.0)
+ return false;
+
+ out = e0 + t * seg;
+ return true;
+}
+
+// ------------------------------------------------------------------------------------------------
+void FilterPolygon(std::vector<IfcVector3> &resultpoly) {
+ if (resultpoly.size() < 3) {
+ resultpoly.clear();
+ return;
+ }
+
+ IfcVector3 vmin, vmax;
+ ArrayBounds(resultpoly.data(), static_cast<unsigned int>(resultpoly.size()), vmin, vmax);
+
+ // filter our IfcFloat points - those may happen if a point lies
+ // directly on the intersection line or directly on the clipping plane
+ const IfcFloat epsilon = (vmax - vmin).SquareLength() / 1e6f;
+ FuzzyVectorCompare fz(epsilon);
+ std::vector<IfcVector3>::iterator e = std::unique(resultpoly.begin(), resultpoly.end(), fz);
+
+ if (e != resultpoly.end())
+ resultpoly.erase(e, resultpoly.end());
+
+ if (!resultpoly.empty() && fz(resultpoly.front(), resultpoly.back()))
+ resultpoly.pop_back();
+}
+
+// ------------------------------------------------------------------------------------------------
+void WritePolygon(std::vector<IfcVector3> &resultpoly, TempMesh &result) {
+ FilterPolygon(resultpoly);
+
+ if (resultpoly.size() > 2) {
+ result.mVerts.insert(result.mVerts.end(), resultpoly.begin(), resultpoly.end());
+ result.mVertcnt.push_back(static_cast<unsigned int>(resultpoly.size()));
+ }
+}
+
+// ------------------------------------------------------------------------------------------------
+void ProcessBooleanHalfSpaceDifference(const Schema_2x3::IfcHalfSpaceSolid *hs, TempMesh &result,
+ const TempMesh &first_operand,
+ ConversionData & /*conv*/) {
+ ai_assert(hs != nullptr);
+
+ const Schema_2x3::IfcPlane *const plane = hs->BaseSurface->ToPtr<Schema_2x3::IfcPlane>();
+ if (!plane) {
+ IFCImporter::LogError("expected IfcPlane as base surface for the IfcHalfSpaceSolid");
+ return;
+ }
+
+ // extract plane base position vector and normal vector
+ IfcVector3 p, n(0.f, 0.f, 1.f);
+ if (plane->Position->Axis) {
+ ConvertDirection(n, plane->Position->Axis.Get());
+ }
+ ConvertCartesianPoint(p, plane->Position->Location);
+
+ if (!IsTrue(hs->AgreementFlag)) {
+ n *= -1.f;
+ }
+
+ // clip the current contents of `meshout` against the plane we obtained from the second operand
+ const std::vector<IfcVector3> &in = first_operand.mVerts;
+ std::vector<IfcVector3> &outvert = result.mVerts;
+
+ std::vector<unsigned int>::const_iterator begin = first_operand.mVertcnt.begin(),
+ end = first_operand.mVertcnt.end(), iit;
+
+ outvert.reserve(in.size());
+ result.mVertcnt.reserve(first_operand.mVertcnt.size());
+
+ unsigned int vidx = 0;
+ for (iit = begin; iit != end; vidx += *iit++) {
+
+ unsigned int newcount = 0;
+ bool isAtWhiteSide = (in[vidx] - p) * n > -ai_epsilon;
+ for (unsigned int i = 0; i < *iit; ++i) {
+ const IfcVector3 &e0 = in[vidx + i], e1 = in[vidx + (i + 1) % *iit];
+
+ // does the next segment intersect the plane?
+ IfcVector3 isectpos;
+ if (IntersectSegmentPlane(p, n, e0, e1, isAtWhiteSide, isectpos)) {
+ if (isAtWhiteSide) {
+ // e0 is on the right side, so keep it
+ outvert.push_back(e0);
+ outvert.push_back(isectpos);
+ newcount += 2;
+ } else {
+ // e0 is on the wrong side, so drop it and keep e1 instead
+ outvert.push_back(isectpos);
+ ++newcount;
+ }
+ isAtWhiteSide = !isAtWhiteSide;
+ } else {
+ if (isAtWhiteSide) {
+ outvert.push_back(e0);
+ ++newcount;
+ }
+ }
+ }
+
+ if (!newcount) {
+ continue;
+ }
+
+ IfcVector3 vmin, vmax;
+ ArrayBounds(&*(outvert.end() - newcount), newcount, vmin, vmax);
+
+ // filter our IfcFloat points - those may happen if a point lies
+ // directly on the intersection line. However, due to IfcFloat
+ // precision a bitwise comparison is not feasible to detect
+ // this case.
+ const IfcFloat epsilon = (vmax - vmin).SquareLength() / 1e6f;
+ FuzzyVectorCompare fz(epsilon);
+
+ std::vector<IfcVector3>::iterator e = std::unique(outvert.end() - newcount, outvert.end(), fz);
+
+ if (e != outvert.end()) {
+ newcount -= static_cast<unsigned int>(std::distance(e, outvert.end()));
+ outvert.erase(e, outvert.end());
+ }
+ if (fz(*(outvert.end() - newcount), outvert.back())) {
+ outvert.pop_back();
+ --newcount;
+ }
+ if (newcount > 2) {
+ result.mVertcnt.push_back(newcount);
+ } else
+ while (newcount-- > 0) {
+ result.mVerts.pop_back();
+ }
+ }
+ IFCImporter::LogVerboseDebug("generating CSG geometry by plane clipping (IfcBooleanClippingResult)");
+}
+
+// ------------------------------------------------------------------------------------------------
+// Check if e0-e1 intersects a sub-segment of the given boundary line.
+// note: this functions works on 3D vectors, but performs its intersection checks solely in xy.
+// New version takes the supposed inside/outside state as a parameter and treats corner cases as if
+// the line stays on that side. This should make corner cases more stable.
+// Two million assumptions! Boundary should have all z at 0.0, will be treated as closed, should not have
+// segments with length <1e-6, self-intersecting might break the corner case handling... just don't go there, ok?
+bool IntersectsBoundaryProfile(const IfcVector3 &e0, const IfcVector3 &e1, const std::vector<IfcVector3> &boundary,
+ const bool isStartAssumedInside, std::vector<std::pair<size_t, IfcVector3>> &intersect_results,
+ const bool halfOpen = false) {
+ ai_assert(intersect_results.empty());
+
+ // determine winding order - necessary to detect segments going "inwards" or "outwards" from a point directly on the border
+ // positive sum of angles means clockwise order when looking down the -Z axis
+ IfcFloat windingOrder = 0.0;
+ for (size_t i = 0, bcount = boundary.size(); i < bcount; ++i) {
+ IfcVector3 b01 = boundary[(i + 1) % bcount] - boundary[i];
+ IfcVector3 b12 = boundary[(i + 2) % bcount] - boundary[(i + 1) % bcount];
+ IfcVector3 b1_side = IfcVector3(b01.y, -b01.x, 0.0); // rotated 90° clockwise in Z plane
+ // Warning: rough estimate only. A concave poly with lots of small segments each featuring a small counter rotation
+ // could fool the accumulation. Correct implementation would be sum( acos( b01 * b2) * sign( b12 * b1_side))
+ windingOrder += (b1_side.x * b12.x + b1_side.y * b12.y);
+ }
+ windingOrder = windingOrder > 0.0 ? 1.0 : -1.0;
+
+ const IfcVector3 e = e1 - e0;
+
+ for (size_t i = 0, bcount = boundary.size(); i < bcount; ++i) {
+ // boundary segment i: b0-b1
+ const IfcVector3 &b0 = boundary[i];
+ const IfcVector3 &b1 = boundary[(i + 1) % bcount];
+ IfcVector3 b = b1 - b0;
+
+ // segment-segment intersection
+ // solve b0 + b*s = e0 + e*t for (s,t)
+ const IfcFloat det = (-b.x * e.y + e.x * b.y);
+ if (std::abs(det) < ai_epsilon) {
+ // no solutions (parallel lines)
+ continue;
+ }
+ IfcFloat b_sqlen_inv = 1.0 / b.SquareLength();
+
+ const IfcFloat x = b0.x - e0.x;
+ const IfcFloat y = b0.y - e0.y;
+ const IfcFloat s = (x * e.y - e.x * y) / det; // scale along boundary edge
+ const IfcFloat t = (x * b.y - b.x * y) / det; // scale along given segment
+ const IfcVector3 p = e0 + e * t;
+#ifdef ASSIMP_BUILD_DEBUG
+ const IfcVector3 check = b0 + b * s - p;
+ ai_assert((IfcVector2(check.x, check.y)).SquareLength() < 1e-5);
+#endif
+
+ // also calculate the distance of e0 and e1 to the segment. We need to detect the "starts directly on segment"
+ // and "ends directly at segment" cases
+ bool startsAtSegment, endsAtSegment;
+ {
+ // calculate closest point to each end on the segment, clamp that point to the segment's length, then check
+ // distance to that point. This approach is like testing if e0 is inside a capped cylinder.
+ IfcFloat et0 = (b.x * (e0.x - b0.x) + b.y * (e0.y - b0.y)) * b_sqlen_inv;
+ IfcVector3 closestPosToE0OnBoundary = b0 + std::max(IfcFloat(0.0), std::min(IfcFloat(1.0), et0)) * b;
+ startsAtSegment = (closestPosToE0OnBoundary - IfcVector3(e0.x, e0.y, 0.0)).SquareLength() < 1e-12;
+ IfcFloat et1 = (b.x * (e1.x - b0.x) + b.y * (e1.y - b0.y)) * b_sqlen_inv;
+ IfcVector3 closestPosToE1OnBoundary = b0 + std::max(IfcFloat(0.0), std::min(IfcFloat(1.0), et1)) * b;
+ endsAtSegment = (closestPosToE1OnBoundary - IfcVector3(e1.x, e1.y, 0.0)).SquareLength() < 1e-12;
+ }
+
+ // Line segment ends at boundary -> ignore any hit, it will be handled by possibly following segments
+ if (endsAtSegment && !halfOpen)
+ continue;
+
+ // Line segment starts at boundary -> generate a hit only if following that line would change the INSIDE/OUTSIDE
+ // state. This should catch the case where a connected set of segments has a point directly on the boundary,
+ // one segment not hitting it because it ends there and the next segment not hitting it because it starts there
+ // Should NOT generate a hit if the segment only touches the boundary but turns around and stays inside.
+ if (startsAtSegment) {
+ IfcVector3 inside_dir = IfcVector3(b.y, -b.x, 0.0) * windingOrder;
+ bool isGoingInside = (inside_dir * e) > 0.0;
+ if (isGoingInside == isStartAssumedInside)
+ continue;
+
+ // only insert the point into the list if it is sufficiently far away from the previous intersection point.
+ // This way, we avoid duplicate detection if the intersection is directly on the vertex between two segments.
+ if (!intersect_results.empty() && intersect_results.back().first == i - 1) {
+ const IfcVector3 diff = intersect_results.back().second - e0;
+ if (IfcVector2(diff.x, diff.y).SquareLength() < 1e-10)
+ continue;
+ }
+ intersect_results.push_back(std::make_pair(i, e0));
+ continue;
+ }
+
+ // for a valid intersection, s and t should be in range [0,1]. Including a bit of epsilon on s, potential double
+ // hits on two consecutive boundary segments are filtered
+ if (s >= -ai_epsilon * b_sqlen_inv && s <= 1.0 + ai_epsilon * b_sqlen_inv && t >= 0.0 && (t <= 1.0 || halfOpen)) {
+ // only insert the point into the list if it is sufficiently far away from the previous intersection point.
+ // This way, we avoid duplicate detection if the intersection is directly on the vertex between two segments.
+ if (!intersect_results.empty() && intersect_results.back().first == i - 1) {
+ const IfcVector3 diff = intersect_results.back().second - p;
+ if (IfcVector2(diff.x, diff.y).SquareLength() < 1e-10)
+ continue;
+ }
+ intersect_results.push_back(std::make_pair(i, p));
+ }
+ }
+
+ return !intersect_results.empty();
+}
+
+// ------------------------------------------------------------------------------------------------
+// note: this functions works on 3D vectors, but performs its intersection checks solely in xy.
+bool PointInPoly(const IfcVector3 &p, const std::vector<IfcVector3> &boundary) {
+ // even-odd algorithm: take a random vector that extends from p to infinite
+ // and counts how many times it intersects edges of the boundary.
+ // because checking for segment intersections is prone to numeric inaccuracies
+ // or double detections (i.e. when hitting multiple adjacent segments at their
+ // shared vertices) we do it thrice with different rays and vote on it.
+
+ // the even-odd algorithm doesn't work for points which lie directly on
+ // the border of the polygon. If any of our attempts produces this result,
+ // we return false immediately.
+
+ std::vector<std::pair<size_t, IfcVector3>> intersected_boundary;
+ size_t votes = 0;
+
+ IntersectsBoundaryProfile(p, p + IfcVector3(1.0, 0, 0), boundary, true, intersected_boundary, true);
+ votes += intersected_boundary.size() % 2;
+
+ intersected_boundary.clear();
+ IntersectsBoundaryProfile(p, p + IfcVector3(0, 1.0, 0), boundary, true, intersected_boundary, true);
+ votes += intersected_boundary.size() % 2;
+
+ intersected_boundary.clear();
+ IntersectsBoundaryProfile(p, p + IfcVector3(0.6, -0.6, 0.0), boundary, true, intersected_boundary, true);
+ votes += intersected_boundary.size() % 2;
+
+ return votes > 1;
+}
+
+// ------------------------------------------------------------------------------------------------
+void ProcessPolygonalBoundedBooleanHalfSpaceDifference(const Schema_2x3::IfcPolygonalBoundedHalfSpace *hs, TempMesh &result,
+ const TempMesh &first_operand,
+ ConversionData &conv) {
+ ai_assert(hs != nullptr);
+
+ const Schema_2x3::IfcPlane *const plane = hs->BaseSurface->ToPtr<Schema_2x3::IfcPlane>();
+ if (!plane) {
+ IFCImporter::LogError("expected IfcPlane as base surface for the IfcHalfSpaceSolid");
+ return;
+ }
+
+ // extract plane base position vector and normal vector
+ IfcVector3 p, n(0.f, 0.f, 1.f);
+ if (plane->Position->Axis) {
+ ConvertDirection(n, plane->Position->Axis.Get());
+ }
+ ConvertCartesianPoint(p, plane->Position->Location);
+
+ if (!IsTrue(hs->AgreementFlag)) {
+ n *= -1.f;
+ }
+
+ n.Normalize();
+
+ // obtain the polygonal bounding volume
+ std::shared_ptr<TempMesh> profile = std::shared_ptr<TempMesh>(new TempMesh());
+ if (!ProcessCurve(hs->PolygonalBoundary, *profile.get(), conv)) {
+ IFCImporter::LogError("expected valid polyline for boundary of boolean halfspace");
+ return;
+ }
+
+ // determine winding order by calculating the normal.
+ IfcVector3 profileNormal = TempMesh::ComputePolygonNormal(profile->mVerts.data(), profile->mVerts.size());
+
+ IfcMatrix4 proj_inv;
+ ConvertAxisPlacement(proj_inv, hs->Position);
+
+ // and map everything into a plane coordinate space so all intersection
+ // tests can be done in 2D space.
+ IfcMatrix4 proj = proj_inv;
+ proj.Inverse();
+
+ // clip the current contents of `meshout` against the plane we obtained from the second operand
+ const std::vector<IfcVector3> &in = first_operand.mVerts;
+ std::vector<IfcVector3> &outvert = result.mVerts;
+ std::vector<unsigned int> &outvertcnt = result.mVertcnt;
+
+ outvert.reserve(in.size());
+ outvertcnt.reserve(first_operand.mVertcnt.size());
+
+ unsigned int vidx = 0;
+ std::vector<unsigned int>::const_iterator begin = first_operand.mVertcnt.begin();
+ std::vector<unsigned int>::const_iterator end = first_operand.mVertcnt.end();
+ std::vector<unsigned int>::const_iterator iit;
+ for (iit = begin; iit != end; vidx += *iit++) {
+ // Our new approach: we cut the poly along the plane, then we intersect the part on the black side of the plane
+ // against the bounding polygon. All the white parts, and the black part outside the boundary polygon, are kept.
+ std::vector<IfcVector3> whiteside, blackside;
+
+ {
+ const IfcVector3 *srcVertices = &in[vidx];
+ const size_t srcVtxCount = *iit;
+ if (srcVtxCount == 0)
+ continue;
+
+ IfcVector3 polyNormal = TempMesh::ComputePolygonNormal(srcVertices, srcVtxCount, true);
+
+ // if the poly is parallel to the plane, put it completely on the black or white side
+ if (std::abs(polyNormal * n) > 0.9999) {
+ bool isOnWhiteSide = (srcVertices[0] - p) * n > -ai_epsilon;
+ std::vector<IfcVector3> &targetSide = isOnWhiteSide ? whiteside : blackside;
+ targetSide.insert(targetSide.end(), srcVertices, srcVertices + srcVtxCount);
+ } else {
+ // otherwise start building one polygon for each side. Whenever the current line segment intersects the plane
+ // we put a point there as an end of the current segment. Then we switch to the other side, put a point there, too,
+ // as a beginning of the current segment, and simply continue accumulating vertices.
+ bool isCurrentlyOnWhiteSide = ((srcVertices[0]) - p) * n > -ai_epsilon;
+ for (size_t a = 0; a < srcVtxCount; ++a) {
+ IfcVector3 e0 = srcVertices[a];
+ IfcVector3 e1 = srcVertices[(a + 1) % srcVtxCount];
+ IfcVector3 ei;
+
+ // put starting point to the current mesh
+ std::vector<IfcVector3> &trgt = isCurrentlyOnWhiteSide ? whiteside : blackside;
+ trgt.push_back(srcVertices[a]);
+
+ // if there's an intersection, put an end vertex there, switch to the other side's mesh,
+ // and add a starting vertex there, too
+ bool isPlaneHit = IntersectSegmentPlane(p, n, e0, e1, isCurrentlyOnWhiteSide, ei);
+ if (isPlaneHit) {
+ if (trgt.empty() || (trgt.back() - ei).SquareLength() > 1e-12)
+ trgt.push_back(ei);
+ isCurrentlyOnWhiteSide = !isCurrentlyOnWhiteSide;
+ std::vector<IfcVector3> &newtrgt = isCurrentlyOnWhiteSide ? whiteside : blackside;
+ newtrgt.push_back(ei);
+ }
+ }
+ }
+ }
+
+ // the part on the white side can be written into the target mesh right away
+ WritePolygon(whiteside, result);
+
+ // The black part is the piece we need to get rid of, but only the part of it within the boundary polygon.
+ // So we now need to construct all the polygons that result from BlackSidePoly minus BoundaryPoly.
+ FilterPolygon(blackside);
+
+ // Complicated, II. We run along the polygon. a) When we're inside the boundary, we run on until we hit an
+ // intersection, which means we're leaving it. We then start a new out poly there. b) When we're outside the
+ // boundary, we start collecting vertices until we hit an intersection, then we run along the boundary until we hit
+ // an intersection, then we switch back to the poly and run on on this one again, and so on until we got a closed
+ // loop. Then we continue with the path we left to catch potential additional polys on the other side of the
+ // boundary as described in a)
+ if (!blackside.empty()) {
+ // poly edge index, intersection point, edge index in boundary poly
+ std::vector<std::tuple<size_t, IfcVector3, size_t>> intersections;
+ bool startedInside = PointInPoly(proj * blackside.front(), profile->mVerts);
+ bool isCurrentlyInside = startedInside;
+
+ std::vector<std::pair<size_t, IfcVector3>> intersected_boundary;
+
+ for (size_t a = 0; a < blackside.size(); ++a) {
+ const IfcVector3 e0 = proj * blackside[a];
+ const IfcVector3 e1 = proj * blackside[(a + 1) % blackside.size()];
+
+ intersected_boundary.clear();
+ IntersectsBoundaryProfile(e0, e1, profile->mVerts, isCurrentlyInside, intersected_boundary);
+ // sort the hits by distance from e0 to get the correct in/out/in sequence. Manually :-( I miss you, C++11.
+ if (intersected_boundary.size() > 1) {
+ bool keepSorting = true;
+ while (keepSorting) {
+ keepSorting = false;
+ for (size_t b = 0; b < intersected_boundary.size() - 1; ++b) {
+ if ((intersected_boundary[b + 1].second - e0).SquareLength() < (intersected_boundary[b].second - e0).SquareLength()) {
+ keepSorting = true;
+ std::swap(intersected_boundary[b + 1], intersected_boundary[b]);
+ }
+ }
+ }
+ }
+ // now add them to the list of intersections
+ for (size_t b = 0; b < intersected_boundary.size(); ++b)
+ intersections.push_back(std::make_tuple(a, proj_inv * intersected_boundary[b].second, intersected_boundary[b].first));
+
+ // and calculate our new inside/outside state
+ if (intersected_boundary.size() & 1)
+ isCurrentlyInside = !isCurrentlyInside;
+ }
+
+ // we got a list of in-out-combinations of intersections. That should be an even number of intersections, or
+ // we are facing a non-recoverable error.
+ if ((intersections.size() & 1) != 0) {
+ IFCImporter::LogWarn("Odd number of intersections, can't work with that. Omitting half space boundary check.");
+ continue;
+ }
+
+ if (intersections.size() > 1) {
+ // If we started outside, the first intersection is a out->in intersection. Cycle them so that it
+ // starts with an intersection leaving the boundary
+ if (!startedInside)
+ for (size_t b = 0; b < intersections.size() - 1; ++b)
+ std::swap(intersections[b], intersections[(b + intersections.size() - 1) % intersections.size()]);
+
+ // Filter pairs of out->in->out that lie too close to each other.
+ for (size_t a = 0; intersections.size() > 0 && a < intersections.size() - 1; /**/) {
+ if ((std::get<1>(intersections[a]) - std::get<1>(intersections[(a + 1) % intersections.size()])).SquareLength() < 1e-10)
+ intersections.erase(intersections.begin() + a, intersections.begin() + a + 2);
+ else
+ a++;
+ }
+ if (intersections.size() > 1 && (std::get<1>(intersections.back()) - std::get<1>(intersections.front())).SquareLength() < 1e-10) {
+ intersections.pop_back();
+ intersections.erase(intersections.begin());
+ }
+ }
+
+ // no intersections at all: either completely inside the boundary, so everything gets discarded, or completely outside.
+ // in the latter case we're implementional lost. I'm simply going to ignore this, so a large poly will not get any
+ // holes if the boundary is smaller and does not touch it anywhere.
+ if (intersections.empty()) {
+ // starting point was outside -> everything is outside the boundary -> nothing is clipped -> add black side
+ // to result mesh unchanged
+ if (!startedInside) {
+ outvertcnt.push_back(static_cast<unsigned int>(blackside.size()));
+ outvert.insert(outvert.end(), blackside.begin(), blackside.end());
+ continue;
+ } else {
+ // starting point was inside the boundary -> everything is inside the boundary -> nothing is spared from the
+ // clipping -> nothing left to add to the result mesh
+ continue;
+ }
+ }
+
+ // determine the direction in which we're marching along the boundary polygon. If the src poly is faced upwards
+ // and the boundary is also winded this way, we need to march *backwards* on the boundary.
+ const IfcVector3 polyNormal = IfcMatrix3(proj) * TempMesh::ComputePolygonNormal(blackside.data(), blackside.size());
+ bool marchBackwardsOnBoundary = (profileNormal * polyNormal) >= 0.0;
+
+ // Build closed loops from these intersections. Starting from an intersection leaving the boundary we
+ // walk along the polygon to the next intersection (which should be an IS entering the boundary poly).
+ // From there we walk along the boundary until we hit another intersection leaving the boundary,
+ // walk along the poly to the next IS and so on until we're back at the starting point.
+ // We remove every intersection we "used up", so any remaining intersection is the start of a new loop.
+ while (!intersections.empty()) {
+ std::vector<IfcVector3> resultpoly;
+ size_t currentIntersecIdx = 0;
+
+ while (true) {
+ ai_assert(intersections.size() > currentIntersecIdx + 1);
+ std::tuple<size_t, IfcVector3, size_t> currintsec = intersections[currentIntersecIdx + 0];
+ std::tuple<size_t, IfcVector3, size_t> nextintsec = intersections[currentIntersecIdx + 1];
+ intersections.erase(intersections.begin() + currentIntersecIdx, intersections.begin() + currentIntersecIdx + 2);
+
+ // we start with an in->out intersection
+ resultpoly.push_back(std::get<1>(currintsec));
+ // climb along the polygon to the next intersection, which should be an out->in
+ size_t numPolyPoints = (std::get<0>(currintsec) > std::get<0>(nextintsec) ? blackside.size() : 0) + std::get<0>(nextintsec) - std::get<0>(currintsec);
+ for (size_t a = 1; a <= numPolyPoints; ++a)
+ resultpoly.push_back(blackside[(std::get<0>(currintsec) + a) % blackside.size()]);
+ // put the out->in intersection
+ resultpoly.push_back(std::get<1>(nextintsec));
+
+ // generate segments along the boundary polygon that lie in the poly's plane until we hit another intersection
+ IfcVector3 startingPoint = proj * std::get<1>(nextintsec);
+ size_t currentBoundaryEdgeIdx = (std::get<2>(nextintsec) + (marchBackwardsOnBoundary ? 1 : 0)) % profile->mVerts.size();
+ size_t nextIntsecIdx = SIZE_MAX;
+ while (nextIntsecIdx == SIZE_MAX) {
+ IfcFloat t = 1e10;
+
+ size_t nextBoundaryEdgeIdx = marchBackwardsOnBoundary ? (currentBoundaryEdgeIdx + profile->mVerts.size() - 1) : currentBoundaryEdgeIdx + 1;
+ nextBoundaryEdgeIdx %= profile->mVerts.size();
+ // vertices of the current boundary segments
+ IfcVector3 currBoundaryPoint = profile->mVerts[currentBoundaryEdgeIdx];
+ IfcVector3 nextBoundaryPoint = profile->mVerts[nextBoundaryEdgeIdx];
+ // project the two onto the polygon
+ if (std::abs(polyNormal.z) > 1e-5) {
+ currBoundaryPoint.z = startingPoint.z + (currBoundaryPoint.x - startingPoint.x) * polyNormal.x / polyNormal.z + (currBoundaryPoint.y - startingPoint.y) * polyNormal.y / polyNormal.z;
+ nextBoundaryPoint.z = startingPoint.z + (nextBoundaryPoint.x - startingPoint.x) * polyNormal.x / polyNormal.z + (nextBoundaryPoint.y - startingPoint.y) * polyNormal.y / polyNormal.z;
+ }
+
+ // build a direction that goes along the boundary border but lies in the poly plane
+ IfcVector3 boundaryPlaneNormal = ((nextBoundaryPoint - currBoundaryPoint) ^ profileNormal).Normalize();
+ IfcVector3 dirAtPolyPlane = (boundaryPlaneNormal ^ polyNormal).Normalize() * (marchBackwardsOnBoundary ? -1.0 : 1.0);
+ // if we can project the direction to the plane, we can calculate a maximum marching distance along that dir
+ // until we finish that boundary segment and continue on the next
+ if (std::abs(polyNormal.z) > 1e-5) {
+ t = std::min(t, (nextBoundaryPoint - startingPoint).Length());
+ }
+
+ // check if the direction hits the loop start - if yes, we got a poly to output
+ IfcVector3 dirToThatPoint = proj * resultpoly.front() - startingPoint;
+ IfcFloat tpt = dirToThatPoint * dirAtPolyPlane;
+ if (tpt > -1e-6 && tpt <= t && (dirToThatPoint - tpt * dirAtPolyPlane).SquareLength() < 1e-10) {
+ nextIntsecIdx = intersections.size(); // dirty hack to end marching along the boundary and signal the end of the loop
+ t = tpt;
+ }
+
+ // also check if the direction hits any in->out intersections earlier. If we hit one, we can switch back
+ // to marching along the poly border from that intersection point
+ for (size_t a = 0; a < intersections.size(); a += 2) {
+ dirToThatPoint = proj * std::get<1>(intersections[a]) - startingPoint;
+ tpt = dirToThatPoint * dirAtPolyPlane;
+ if (tpt > -1e-6 && tpt <= t && (dirToThatPoint - tpt * dirAtPolyPlane).SquareLength() < 1e-10) {
+ nextIntsecIdx = a; // switch back to poly and march on from this in->out intersection
+ t = tpt;
+ }
+ }
+
+ // if we keep marching on the boundary, put the segment end point to the result poly and well... keep marching
+ if (nextIntsecIdx == SIZE_MAX) {
+ resultpoly.push_back(proj_inv * nextBoundaryPoint);
+ currentBoundaryEdgeIdx = nextBoundaryEdgeIdx;
+ startingPoint = nextBoundaryPoint;
+ }
+
+ // quick endless loop check
+ if (resultpoly.size() > blackside.size() + profile->mVerts.size()) {
+ IFCImporter::LogError("Encountered endless loop while clipping polygon against poly-bounded half space.");
+ break;
+ }
+ }
+
+ // we're back on the poly - if this is the intersection we started from, we got a closed loop.
+ if (nextIntsecIdx >= intersections.size()) {
+ break;
+ }
+
+ // otherwise it's another intersection. Continue marching from there.
+ currentIntersecIdx = nextIntsecIdx;
+ }
+
+ WritePolygon(resultpoly, result);
+ }
+ }
+ }
+ IFCImporter::LogVerboseDebug("generating CSG geometry by plane clipping with polygonal bounding (IfcBooleanClippingResult)");
+}
+
+// ------------------------------------------------------------------------------------------------
+void ProcessBooleanExtrudedAreaSolidDifference(const Schema_2x3::IfcExtrudedAreaSolid *as, TempMesh &result,
+ const TempMesh &first_operand,
+ ConversionData &conv) {
+ ai_assert(as != nullptr);
+
+ // This case is handled by reduction to an instance of the quadrify() algorithm.
+ // Obviously, this won't work for arbitrarily complex cases. In fact, the first
+ // operand should be near-planar. Luckily, this is usually the case in Ifc
+ // buildings.
+
+ std::shared_ptr<TempMesh> meshtmp = std::shared_ptr<TempMesh>(new TempMesh());
+ ProcessExtrudedAreaSolid(*as, *meshtmp, conv, false);
+
+ std::vector<TempOpening> openings(1, TempOpening(as, IfcVector3(0, 0, 0), meshtmp, std::shared_ptr<TempMesh>()));
+
+ result = first_operand;
+
+ TempMesh temp;
+
+ std::vector<IfcVector3>::const_iterator vit = first_operand.mVerts.begin();
+ for (unsigned int pcount : first_operand.mVertcnt) {
+ temp.Clear();
+
+ temp.mVerts.insert(temp.mVerts.end(), vit, vit + pcount);
+ temp.mVertcnt.push_back(pcount);
+
+ // The algorithms used to generate mesh geometry sometimes
+ // spit out lines or other degenerates which must be
+ // filtered to avoid running into assertions later on.
+
+ // ComputePolygonNormal returns the Newell normal, so the
+ // length of the normal is the area of the polygon.
+ const IfcVector3 &normal = temp.ComputeLastPolygonNormal(false);
+ if (normal.SquareLength() < static_cast<IfcFloat>(1e-5)) {
+ IFCImporter::LogWarn("skipping degenerate polygon (ProcessBooleanExtrudedAreaSolidDifference)");
+ continue;
+ }
+
+ GenerateOpenings(openings, temp, false, true);
+ result.Append(temp);
+
+ vit += pcount;
+ }
+
+ IFCImporter::LogVerboseDebug("generating CSG geometry by geometric difference to a solid (IfcExtrudedAreaSolid)");
+}
+
+// ------------------------------------------------------------------------------------------------
+void ProcessBoolean(const Schema_2x3::IfcBooleanResult &boolean, TempMesh &result, ConversionData &conv) {
+ // supported CSG operations:
+ // DIFFERENCE
+ if (const Schema_2x3::IfcBooleanResult *const clip = boolean.ToPtr<Schema_2x3::IfcBooleanResult>()) {
+ if (clip->Operator != "DIFFERENCE") {
+ IFCImporter::LogWarn("encountered unsupported boolean operator: ", (std::string)clip->Operator);
+ return;
+ }
+
+ // supported cases (1st operand):
+ // IfcBooleanResult -- call ProcessBoolean recursively
+ // IfcSweptAreaSolid -- obtain polygonal geometry first
+
+ // supported cases (2nd operand):
+ // IfcHalfSpaceSolid -- easy, clip against plane
+ // IfcExtrudedAreaSolid -- reduce to an instance of the quadrify() algorithm
+
+ const Schema_2x3::IfcHalfSpaceSolid *const hs = clip->SecondOperand->ResolveSelectPtr<Schema_2x3::IfcHalfSpaceSolid>(conv.db);
+ const Schema_2x3::IfcExtrudedAreaSolid *const as = clip->SecondOperand->ResolveSelectPtr<Schema_2x3::IfcExtrudedAreaSolid>(conv.db);
+ if (!hs && !as) {
+ IFCImporter::LogError("expected IfcHalfSpaceSolid or IfcExtrudedAreaSolid as second clipping operand");
+ return;
+ }
+
+ TempMesh first_operand;
+ if (const Schema_2x3::IfcBooleanResult *const op0 = clip->FirstOperand->ResolveSelectPtr<Schema_2x3::IfcBooleanResult>(conv.db)) {
+ ProcessBoolean(*op0, first_operand, conv);
+ } else if (const Schema_2x3::IfcSweptAreaSolid *const swept = clip->FirstOperand->ResolveSelectPtr<Schema_2x3::IfcSweptAreaSolid>(conv.db)) {
+ ProcessSweptAreaSolid(*swept, first_operand, conv);
+ } else {
+ IFCImporter::LogError("expected IfcSweptAreaSolid or IfcBooleanResult as first clipping operand");
+ return;
+ }
+
+ if (hs) {
+
+ const Schema_2x3::IfcPolygonalBoundedHalfSpace *const hs_bounded = clip->SecondOperand->ResolveSelectPtr<Schema_2x3::IfcPolygonalBoundedHalfSpace>(conv.db);
+ if (hs_bounded) {
+ ProcessPolygonalBoundedBooleanHalfSpaceDifference(hs_bounded, result, first_operand, conv);
+ } else {
+ ProcessBooleanHalfSpaceDifference(hs, result, first_operand, conv);
+ }
+ } else {
+ ProcessBooleanExtrudedAreaSolidDifference(as, result, first_operand, conv);
+ }
+ } else {
+ IFCImporter::LogWarn("skipping unknown IfcBooleanResult entity, type is ", boolean.GetClassName());
+ }
+}
+
+} // namespace IFC
+} // namespace Assimp
+
+#endif