assimp/code/IFCBoolean.cpp

/*
Open Asset Import Library (assimp)
----------------------------------------------------------------------

Copyright (c) 2006-2010, 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 "IFCUtil.h"
#include "PolyTools.h"
#include "ProcessHelper.h"
#include <assimp/Defines.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) < 1e-6 )
        return false;

    // if segment starts on the plane, report a hit only if the end lies on the *other* side
    if( std::abs(dotTwo) < 1e-6 )
    {
        if( (assumeStartOnWhiteSide && dotOne + dotTwo < 1e-6) || (!assumeStartOnWhiteSide && dotOne + dotTwo > -1e-6) )
        {
            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) < 1e-6 )
        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.verts.insert(result.verts.end(), resultpoly.begin(), resultpoly.end());
        result.vertcnt.push_back(static_cast<unsigned int>(resultpoly.size()));
    }
}


// ------------------------------------------------------------------------------------------------
void ProcessBooleanHalfSpaceDifference(const IfcHalfSpaceSolid* hs, TempMesh& result,
    const TempMesh& first_operand,
    ConversionData& /*conv*/)
{
    ai_assert(hs != NULL);

    const IfcPlane* const plane = hs->BaseSurface->ToPtr<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.verts;
    std::vector<IfcVector3>& outvert = result.verts;

    std::vector<unsigned int>::const_iterator begin = first_operand.vertcnt.begin(),
        end = first_operand.vertcnt.end(), iit;

    outvert.reserve(in.size());
    result.vertcnt.reserve(first_operand.vertcnt.size());

    unsigned int vidx = 0;
    for(iit = begin; iit != end; vidx += *iit++) {

        unsigned int newcount = 0;
        bool isAtWhiteSide = (in[vidx] - p) * n > -1e-6;
        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.vertcnt.push_back(newcount);
        }
        else while(newcount-->0) {
            result.verts.pop_back();
        }

    }
    IFCImporter::LogDebug("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<39> 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;
        IfcFloat b_sqlen_inv = 1.0 / b.SquareLength();

        // 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) < 1e-6 ) {
            // no solutions (parallel lines)
            continue;
        }

        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 >= -1e-6 * b_sqlen_inv && s <= 1.0 + 1e-6*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;

//  ai_assert(votes == 3 || votes == 0);
    return votes > 1;
}


// ------------------------------------------------------------------------------------------------
void ProcessPolygonalBoundedBooleanHalfSpaceDifference(const IfcPolygonalBoundedHalfSpace* hs, TempMesh& result,
                                                       const TempMesh& first_operand,
                                                       ConversionData& conv)
{
    ai_assert(hs != NULL);

    const IfcPlane* const plane = hs->BaseSurface->ToPtr<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->verts.data(), profile->verts.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.verts;
    std::vector<IfcVector3>& outvert = result.verts;
    std::vector<unsigned int>& outvertcnt = result.vertcnt;

    outvert.reserve(in.size());
    outvertcnt.reserve(first_operand.vertcnt.size());

    unsigned int vidx = 0;
    std::vector<unsigned int>::const_iterator begin = first_operand.vertcnt.begin();
    std::vector<unsigned int>::const_iterator end = first_operand.vertcnt.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 > -1e-6;
                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 > -1e-6;
                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->verts);
            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->verts, 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're fucked.
            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->verts.size();
                    size_t nextIntsecIdx = SIZE_MAX;
                    while( nextIntsecIdx == SIZE_MAX )
                    {
                        IfcFloat t = 1e10;

                        size_t nextBoundaryEdgeIdx = marchBackwardsOnBoundary ? (currentBoundaryEdgeIdx + profile->verts.size() - 1) : currentBoundaryEdgeIdx + 1;
                        nextBoundaryEdgeIdx %= profile->verts.size();
                        // vertices of the current boundary segments
                        IfcVector3 currBoundaryPoint = profile->verts[currentBoundaryEdgeIdx];
                        IfcVector3 nextBoundaryPoint = profile->verts[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->verts.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::LogDebug("generating CSG geometry by plane clipping with polygonal bounding (IfcBooleanClippingResult)");
}

// ------------------------------------------------------------------------------------------------
void ProcessBooleanExtrudedAreaSolidDifference(const IfcExtrudedAreaSolid* as, TempMesh& result,
                                               const TempMesh& first_operand,
                                               ConversionData& conv)
{
    ai_assert(as != NULL);

    // 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.verts.begin();
    for(unsigned int pcount : first_operand.vertcnt) {
        temp.Clear();

        temp.verts.insert(temp.verts.end(), vit, vit + pcount);
        temp.vertcnt.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, std::vector<IfcVector3>(1,IfcVector3(1,0,0)), temp, false, true);
        result.Append(temp);

        vit += pcount;
    }

    IFCImporter::LogDebug("generating CSG geometry by geometric difference to a solid (IfcExtrudedAreaSolid)");
}

// ------------------------------------------------------------------------------------------------
void ProcessBoolean(const IfcBooleanResult& boolean, TempMesh& result, ConversionData& conv)
{
    // supported CSG operations:
    //   DIFFERENCE
    if(const IfcBooleanResult* const clip = boolean.ToPtr<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 IfcHalfSpaceSolid* const hs = clip->SecondOperand->ResolveSelectPtr<IfcHalfSpaceSolid>(conv.db);
        const IfcExtrudedAreaSolid* const as = clip->SecondOperand->ResolveSelectPtr<IfcExtrudedAreaSolid>(conv.db);
        if(!hs && !as) {
            IFCImporter::LogError("expected IfcHalfSpaceSolid or IfcExtrudedAreaSolid as second clipping operand");
            return;
        }

        TempMesh first_operand;
        if(const IfcBooleanResult* const op0 = clip->FirstOperand->ResolveSelectPtr<IfcBooleanResult>(conv.db)) {
            ProcessBoolean(*op0,first_operand,conv);
        }
        else if (const IfcSweptAreaSolid* const swept = clip->FirstOperand->ResolveSelectPtr<IfcSweptAreaSolid>(conv.db)) {
            ProcessSweptAreaSolid(*swept,first_operand,conv);
        }
        else {
            IFCImporter::LogError("expected IfcSweptAreaSolid or IfcBooleanResult as first clipping operand");
            return;
        }

        if(hs) {

            const IfcPolygonalBoundedHalfSpace* const hs_bounded = clip->SecondOperand->ResolveSelectPtr<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());
    }
}

} // ! IFC
} // ! Assimp

#endif