1/*
2Open Asset Import Library (assimp)
3----------------------------------------------------------------------
4
5Copyright (c) 2006-2010, assimp team
6All rights reserved.
7
8Redistribution and use of this software in source and binary forms,
9with or without modification, are permitted provided that the
10following conditions are met:
11
12* Redistributions of source code must retain the above
13 copyright notice, this list of conditions and the
14 following disclaimer.
15
16* Redistributions in binary form must reproduce the above
17 copyright notice, this list of conditions and the
18 following disclaimer in the documentation and/or other
19 materials provided with the distribution.
20
21* Neither the name of the assimp team, nor the names of its
22 contributors may be used to endorse or promote products
23 derived from this software without specific prior
24 written permission of the assimp team.
25
26THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
27"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
29A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
30OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
31SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
32LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
33DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
34THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
35(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
36OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37
38----------------------------------------------------------------------
39*/
40
41/** @file IFCBoolean.cpp
42 * @brief Implements a subset of Ifc boolean operations
43 */
44
45
46#ifndef ASSIMP_BUILD_NO_IFC_IMPORTER
47#include "IFCUtil.h"
48#include "PolyTools.h"
49#include "ProcessHelper.h"
50#include <assimp/Defines.h>
51
52#include <iterator>
53#include <tuple>
54
55
56namespace Assimp {
57 namespace IFC {
58
59// ------------------------------------------------------------------------------------------------
60// Calculates intersection between line segment and plane. To catch corner cases, specify which side you prefer.
61// The function then generates a hit only if the end is beyond a certain margin in that direction, filtering out
62// "very close to plane" ghost hits as long as start and end stay directly on or within the given plane side.
63bool IntersectSegmentPlane(const IfcVector3& p,const IfcVector3& n, const IfcVector3& e0,
64 const IfcVector3& e1, bool assumeStartOnWhiteSide, IfcVector3& out)
65{
66 const IfcVector3 pdelta = e0 - p, seg = e1 - e0;
67 const IfcFloat dotOne = n*seg, dotTwo = -(n*pdelta);
68
69 // if segment ends on plane, do not report a hit. We stay on that side until a following segment starting at this
70 // point leaves the plane through the other side
71 if( std::abs(dotOne + dotTwo) < 1e-6 )
72 return false;
73
74 // if segment starts on the plane, report a hit only if the end lies on the *other* side
75 if( std::abs(dotTwo) < 1e-6 )
76 {
77 if( (assumeStartOnWhiteSide && dotOne + dotTwo < 1e-6) || (!assumeStartOnWhiteSide && dotOne + dotTwo > -1e-6) )
78 {
79 out = e0;
80 return true;
81 }
82 else
83 {
84 return false;
85 }
86 }
87
88 // ignore if segment is parallel to plane and far away from it on either side
89 // Warning: if there's a few thousand of such segments which slowly accumulate beyond the epsilon, no hit would be registered
90 if( std::abs(dotOne) < 1e-6 )
91 return false;
92
93 // t must be in [0..1] if the intersection point is within the given segment
94 const IfcFloat t = dotTwo / dotOne;
95 if( t > 1.0 || t < 0.0 )
96 return false;
97
98 out = e0 + t*seg;
99 return true;
100}
101
102// ------------------------------------------------------------------------------------------------
103void FilterPolygon(std::vector<IfcVector3>& resultpoly)
104{
105 if( resultpoly.size() < 3 )
106 {
107 resultpoly.clear();
108 return;
109 }
110
111 IfcVector3 vmin, vmax;
112 ArrayBounds(resultpoly.data(), static_cast<unsigned int>(resultpoly.size()), vmin, vmax);
113
114 // filter our IfcFloat points - those may happen if a point lies
115 // directly on the intersection line or directly on the clipping plane
116 const IfcFloat epsilon = (vmax - vmin).SquareLength() / 1e6f;
117 FuzzyVectorCompare fz(epsilon);
118 std::vector<IfcVector3>::iterator e = std::unique(resultpoly.begin(), resultpoly.end(), fz);
119
120 if( e != resultpoly.end() )
121 resultpoly.erase(e, resultpoly.end());
122
123 if( !resultpoly.empty() && fz(resultpoly.front(), resultpoly.back()) )
124 resultpoly.pop_back();
125}
126
127// ------------------------------------------------------------------------------------------------
128void WritePolygon(std::vector<IfcVector3>& resultpoly, TempMesh& result)
129{
130 FilterPolygon(resultpoly);
131
132 if( resultpoly.size() > 2 )
133 {
134 result.verts.insert(result.verts.end(), resultpoly.begin(), resultpoly.end());
135 result.vertcnt.push_back(static_cast<unsigned int>(resultpoly.size()));
136 }
137}
138
139
140// ------------------------------------------------------------------------------------------------
141void ProcessBooleanHalfSpaceDifference(const IfcHalfSpaceSolid* hs, TempMesh& result,
142 const TempMesh& first_operand,
143 ConversionData& /*conv*/)
144{
145 ai_assert(hs != NULL);
146
147 const IfcPlane* const plane = hs->BaseSurface->ToPtr<IfcPlane>();
148 if(!plane) {
149 IFCImporter::LogError("expected IfcPlane as base surface for the IfcHalfSpaceSolid");
150 return;
151 }
152
153 // extract plane base position vector and normal vector
154 IfcVector3 p,n(0.f,0.f,1.f);
155 if (plane->Position->Axis) {
156 ConvertDirection(n,plane->Position->Axis.Get());
157 }
158 ConvertCartesianPoint(p,plane->Position->Location);
159
160 if(!IsTrue(hs->AgreementFlag)) {
161 n *= -1.f;
162 }
163
164 // clip the current contents of `meshout` against the plane we obtained from the second operand
165 const std::vector<IfcVector3>& in = first_operand.verts;
166 std::vector<IfcVector3>& outvert = result.verts;
167
168 std::vector<unsigned int>::const_iterator begin = first_operand.vertcnt.begin(),
169 end = first_operand.vertcnt.end(), iit;
170
171 outvert.reserve(in.size());
172 result.vertcnt.reserve(first_operand.vertcnt.size());
173
174 unsigned int vidx = 0;
175 for(iit = begin; iit != end; vidx += *iit++) {
176
177 unsigned int newcount = 0;
178 bool isAtWhiteSide = (in[vidx] - p) * n > -1e-6;
179 for( unsigned int i = 0; i < *iit; ++i ) {
180 const IfcVector3& e0 = in[vidx + i], e1 = in[vidx + (i + 1) % *iit];
181
182 // does the next segment intersect the plane?
183 IfcVector3 isectpos;
184 if( IntersectSegmentPlane(p, n, e0, e1, isAtWhiteSide, isectpos) ) {
185 if( isAtWhiteSide ) {
186 // e0 is on the right side, so keep it
187 outvert.push_back(e0);
188 outvert.push_back(isectpos);
189 newcount += 2;
190 }
191 else {
192 // e0 is on the wrong side, so drop it and keep e1 instead
193 outvert.push_back(isectpos);
194 ++newcount;
195 }
196 isAtWhiteSide = !isAtWhiteSide;
197 }
198 else
199 {
200 if( isAtWhiteSide ) {
201 outvert.push_back(e0);
202 ++newcount;
203 }
204 }
205 }
206
207 if (!newcount) {
208 continue;
209 }
210
211 IfcVector3 vmin,vmax;
212 ArrayBounds(&*(outvert.end()-newcount),newcount,vmin,vmax);
213
214 // filter our IfcFloat points - those may happen if a point lies
215 // directly on the intersection line. However, due to IfcFloat
216 // precision a bitwise comparison is not feasible to detect
217 // this case.
218 const IfcFloat epsilon = (vmax-vmin).SquareLength() / 1e6f;
219 FuzzyVectorCompare fz(epsilon);
220
221 std::vector<IfcVector3>::iterator e = std::unique( outvert.end()-newcount, outvert.end(), fz );
222
223 if (e != outvert.end()) {
224 newcount -= static_cast<unsigned int>(std::distance(e,outvert.end()));
225 outvert.erase(e,outvert.end());
226 }
227 if (fz(*( outvert.end()-newcount),outvert.back())) {
228 outvert.pop_back();
229 --newcount;
230 }
231 if(newcount > 2) {
232 result.vertcnt.push_back(newcount);
233 }
234 else while(newcount-->0) {
235 result.verts.pop_back();
236 }
237
238 }
239 IFCImporter::LogDebug("generating CSG geometry by plane clipping (IfcBooleanClippingResult)");
240}
241
242// ------------------------------------------------------------------------------------------------
243// Check if e0-e1 intersects a sub-segment of the given boundary line.
244// note: this functions works on 3D vectors, but performs its intersection checks solely in xy.
245// New version takes the supposed inside/outside state as a parameter and treats corner cases as if
246// the line stays on that side. This should make corner cases more stable.
247// Two million assumptions! Boundary should have all z at 0.0, will be treated as closed, should not have
248// segments with length <1e-6, self-intersecting might break the corner case handling... just don't go there, ok?
249bool IntersectsBoundaryProfile(const IfcVector3& e0, const IfcVector3& e1, const std::vector<IfcVector3>& boundary,
250 const bool isStartAssumedInside, std::vector<std::pair<size_t, IfcVector3> >& intersect_results,
251 const bool halfOpen = false)
252{
253 ai_assert(intersect_results.empty());
254
255 // determine winding order - necessary to detect segments going "inwards" or "outwards" from a point directly on the border
256 // positive sum of angles means clockwise order when looking down the -Z axis
257 IfcFloat windingOrder = 0.0;
258 for( size_t i = 0, bcount = boundary.size(); i < bcount; ++i ) {
259 IfcVector3 b01 = boundary[(i + 1) % bcount] - boundary[i];
260 IfcVector3 b12 = boundary[(i + 2) % bcount] - boundary[(i + 1) % bcount];
261 IfcVector3 b1_side = IfcVector3(b01.y, -b01.x, 0.0); // rotated 90 clockwise in Z plane
262 // Warning: rough estimate only. A concave poly with lots of small segments each featuring a small counter rotation
263 // could fool the accumulation. Correct implementation would be sum( acos( b01 * b2) * sign( b12 * b1_side))
264 windingOrder += (b1_side.x*b12.x + b1_side.y*b12.y);
265 }
266 windingOrder = windingOrder > 0.0 ? 1.0 : -1.0;
267
268 const IfcVector3 e = e1 - e0;
269
270 for( size_t i = 0, bcount = boundary.size(); i < bcount; ++i ) {
271 // boundary segment i: b0-b1
272 const IfcVector3& b0 = boundary[i];
273 const IfcVector3& b1 = boundary[(i + 1) % bcount];
274 IfcVector3 b = b1 - b0;
275
276 // segment-segment intersection
277 // solve b0 + b*s = e0 + e*t for (s,t)
278 const IfcFloat det = (-b.x * e.y + e.x * b.y);
279 if( std::abs(det) < 1e-6 ) {
280 // no solutions (parallel lines)
281 continue;
282 }
283 IfcFloat b_sqlen_inv = 1.0 / b.SquareLength();
284
285 const IfcFloat x = b0.x - e0.x;
286 const IfcFloat y = b0.y - e0.y;
287 const IfcFloat s = (x*e.y - e.x*y) / det; // scale along boundary edge
288 const IfcFloat t = (x*b.y - b.x*y) / det; // scale along given segment
289 const IfcVector3 p = e0 + e*t;
290#ifdef ASSIMP_BUILD_DEBUG
291 const IfcVector3 check = b0 + b*s - p;
292 ai_assert((IfcVector2(check.x, check.y)).SquareLength() < 1e-5);
293#endif
294
295 // also calculate the distance of e0 and e1 to the segment. We need to detect the "starts directly on segment"
296 // and "ends directly at segment" cases
297 bool startsAtSegment, endsAtSegment;
298 {
299 // calculate closest point to each end on the segment, clamp that point to the segment's length, then check
300 // distance to that point. This approach is like testing if e0 is inside a capped cylinder.
301 IfcFloat et0 = (b.x*(e0.x - b0.x) + b.y*(e0.y - b0.y)) * b_sqlen_inv;
302 IfcVector3 closestPosToE0OnBoundary = b0 + std::max(IfcFloat(0.0), std::min(IfcFloat(1.0), et0)) * b;
303 startsAtSegment = (closestPosToE0OnBoundary - IfcVector3(e0.x, e0.y, 0.0)).SquareLength() < 1e-12;
304 IfcFloat et1 = (b.x*(e1.x - b0.x) + b.y*(e1.y - b0.y)) * b_sqlen_inv;
305 IfcVector3 closestPosToE1OnBoundary = b0 + std::max(IfcFloat(0.0), std::min(IfcFloat(1.0), et1)) * b;
306 endsAtSegment = (closestPosToE1OnBoundary - IfcVector3(e1.x, e1.y, 0.0)).SquareLength() < 1e-12;
307 }
308
309 // Line segment ends at boundary -> ignore any hit, it will be handled by possibly following segments
310 if( endsAtSegment && !halfOpen )
311 continue;
312
313 // Line segment starts at boundary -> generate a hit only if following that line would change the INSIDE/OUTSIDE
314 // state. This should catch the case where a connected set of segments has a point directly on the boundary,
315 // one segment not hitting it because it ends there and the next segment not hitting it because it starts there
316 // Should NOT generate a hit if the segment only touches the boundary but turns around and stays inside.
317 if( startsAtSegment )
318 {
319 IfcVector3 inside_dir = IfcVector3(b.y, -b.x, 0.0) * windingOrder;
320 bool isGoingInside = (inside_dir * e) > 0.0;
321 if( isGoingInside == isStartAssumedInside )
322 continue;
323
324 // only insert the point into the list if it is sufficiently far away from the previous intersection point.
325 // This way, we avoid duplicate detection if the intersection is directly on the vertex between two segments.
326 if( !intersect_results.empty() && intersect_results.back().first == i - 1 )
327 {
328 const IfcVector3 diff = intersect_results.back().second - e0;
329 if( IfcVector2(diff.x, diff.y).SquareLength() < 1e-10 )
330 continue;
331 }
332 intersect_results.push_back(std::make_pair(i, e0));
333 continue;
334 }
335
336 // for a valid intersection, s and t should be in range [0,1]. Including a bit of epsilon on s, potential double
337 // hits on two consecutive boundary segments are filtered
338 if( s >= -1e-6 * b_sqlen_inv && s <= 1.0 + 1e-6*b_sqlen_inv && t >= 0.0 && (t <= 1.0 || halfOpen) )
339 {
340 // only insert the point into the list if it is sufficiently far away from the previous intersection point.
341 // This way, we avoid duplicate detection if the intersection is directly on the vertex between two segments.
342 if( !intersect_results.empty() && intersect_results.back().first == i - 1 )
343 {
344 const IfcVector3 diff = intersect_results.back().second - p;
345 if( IfcVector2(diff.x, diff.y).SquareLength() < 1e-10 )
346 continue;
347 }
348 intersect_results.push_back(std::make_pair(i, p));
349 }
350 }
351
352 return !intersect_results.empty();
353}
354
355
356// ------------------------------------------------------------------------------------------------
357// note: this functions works on 3D vectors, but performs its intersection checks solely in xy.
358bool PointInPoly(const IfcVector3& p, const std::vector<IfcVector3>& boundary)
359{
360 // even-odd algorithm: take a random vector that extends from p to infinite
361 // and counts how many times it intersects edges of the boundary.
362 // because checking for segment intersections is prone to numeric inaccuracies
363 // or double detections (i.e. when hitting multiple adjacent segments at their
364 // shared vertices) we do it thrice with different rays and vote on it.
365
366 // the even-odd algorithm doesn't work for points which lie directly on
367 // the border of the polygon. If any of our attempts produces this result,
368 // we return false immediately.
369
370 std::vector<std::pair<size_t, IfcVector3> > intersected_boundary;
371 size_t votes = 0;
372
373 IntersectsBoundaryProfile(p, p + IfcVector3(1.0, 0, 0), boundary, true, intersected_boundary, true);
374 votes += intersected_boundary.size() % 2;
375
376 intersected_boundary.clear();
377 IntersectsBoundaryProfile(p, p + IfcVector3(0, 1.0, 0), boundary, true, intersected_boundary, true);
378 votes += intersected_boundary.size() % 2;
379
380 intersected_boundary.clear();
381 IntersectsBoundaryProfile(p, p + IfcVector3(0.6, -0.6, 0.0), boundary, true, intersected_boundary, true);
382 votes += intersected_boundary.size() % 2;
383
384 return votes > 1;
385}
386
387
388// ------------------------------------------------------------------------------------------------
389void ProcessPolygonalBoundedBooleanHalfSpaceDifference(const IfcPolygonalBoundedHalfSpace* hs, TempMesh& result,
390 const TempMesh& first_operand,
391 ConversionData& conv)
392{
393 ai_assert(hs != NULL);
394
395 const IfcPlane* const plane = hs->BaseSurface->ToPtr<IfcPlane>();
396 if(!plane) {
397 IFCImporter::LogError("expected IfcPlane as base surface for the IfcHalfSpaceSolid");
398 return;
399 }
400
401 // extract plane base position vector and normal vector
402 IfcVector3 p,n(0.f,0.f,1.f);
403 if (plane->Position->Axis) {
404 ConvertDirection(n,plane->Position->Axis.Get());
405 }
406 ConvertCartesianPoint(p,plane->Position->Location);
407
408 if(!IsTrue(hs->AgreementFlag)) {
409 n *= -1.f;
410 }
411
412 n.Normalize();
413
414 // obtain the polygonal bounding volume
415 std::shared_ptr<TempMesh> profile = std::shared_ptr<TempMesh>(new TempMesh());
416 if(!ProcessCurve(hs->PolygonalBoundary, *profile.get(), conv)) {
417 IFCImporter::LogError("expected valid polyline for boundary of boolean halfspace");
418 return;
419 }
420
421 // determine winding order by calculating the normal.
422 IfcVector3 profileNormal = TempMesh::ComputePolygonNormal(profile->verts.data(), profile->verts.size());
423
424 IfcMatrix4 proj_inv;
425 ConvertAxisPlacement(proj_inv,hs->Position);
426
427 // and map everything into a plane coordinate space so all intersection
428 // tests can be done in 2D space.
429 IfcMatrix4 proj = proj_inv;
430 proj.Inverse();
431
432 // clip the current contents of `meshout` against the plane we obtained from the second operand
433 const std::vector<IfcVector3>& in = first_operand.verts;
434 std::vector<IfcVector3>& outvert = result.verts;
435 std::vector<unsigned int>& outvertcnt = result.vertcnt;
436
437 outvert.reserve(in.size());
438 outvertcnt.reserve(first_operand.vertcnt.size());
439
440 unsigned int vidx = 0;
441 std::vector<unsigned int>::const_iterator begin = first_operand.vertcnt.begin();
442 std::vector<unsigned int>::const_iterator end = first_operand.vertcnt.end();
443 std::vector<unsigned int>::const_iterator iit;
444 for( iit = begin; iit != end; vidx += *iit++ )
445 {
446 // Our new approach: we cut the poly along the plane, then we intersect the part on the black side of the plane
447 // against the bounding polygon. All the white parts, and the black part outside the boundary polygon, are kept.
448 std::vector<IfcVector3> whiteside, blackside;
449
450 {
451 const IfcVector3* srcVertices = &in[vidx];
452 const size_t srcVtxCount = *iit;
453 if( srcVtxCount == 0 )
454 continue;
455
456 IfcVector3 polyNormal = TempMesh::ComputePolygonNormal(srcVertices, srcVtxCount, true);
457
458 // if the poly is parallel to the plane, put it completely on the black or white side
459 if( std::abs(polyNormal * n) > 0.9999 )
460 {
461 bool isOnWhiteSide = (srcVertices[0] - p) * n > -1e-6;
462 std::vector<IfcVector3>& targetSide = isOnWhiteSide ? whiteside : blackside;
463 targetSide.insert(targetSide.end(), srcVertices, srcVertices + srcVtxCount);
464 }
465 else
466 {
467 // otherwise start building one polygon for each side. Whenever the current line segment intersects the plane
468 // we put a point there as an end of the current segment. Then we switch to the other side, put a point there, too,
469 // as a beginning of the current segment, and simply continue accumulating vertices.
470 bool isCurrentlyOnWhiteSide = ((srcVertices[0]) - p) * n > -1e-6;
471 for( size_t a = 0; a < srcVtxCount; ++a )
472 {
473 IfcVector3 e0 = srcVertices[a];
474 IfcVector3 e1 = srcVertices[(a + 1) % srcVtxCount];
475 IfcVector3 ei;
476
477 // put starting point to the current mesh
478 std::vector<IfcVector3>& trgt = isCurrentlyOnWhiteSide ? whiteside : blackside;
479 trgt.push_back(srcVertices[a]);
480
481 // if there's an intersection, put an end vertex there, switch to the other side's mesh,
482 // and add a starting vertex there, too
483 bool isPlaneHit = IntersectSegmentPlane(p, n, e0, e1, isCurrentlyOnWhiteSide, ei);
484 if( isPlaneHit )
485 {
486 if( trgt.empty() || (trgt.back() - ei).SquareLength() > 1e-12 )
487 trgt.push_back(ei);
488 isCurrentlyOnWhiteSide = !isCurrentlyOnWhiteSide;
489 std::vector<IfcVector3>& newtrgt = isCurrentlyOnWhiteSide ? whiteside : blackside;
490 newtrgt.push_back(ei);
491 }
492 }
493 }
494 }
495
496 // the part on the white side can be written into the target mesh right away
497 WritePolygon(whiteside, result);
498
499 // The black part is the piece we need to get rid of, but only the part of it within the boundary polygon.
500 // So we now need to construct all the polygons that result from BlackSidePoly minus BoundaryPoly.
501 FilterPolygon(blackside);
502
503 // Complicated, II. We run along the polygon. a) When we're inside the boundary, we run on until we hit an
504 // intersection, which means we're leaving it. We then start a new out poly there. b) When we're outside the
505 // boundary, we start collecting vertices until we hit an intersection, then we run along the boundary until we hit
506 // an intersection, then we switch back to the poly and run on on this one again, and so on until we got a closed
507 // loop. Then we continue with the path we left to catch potential additional polys on the other side of the
508 // boundary as described in a)
509 if( !blackside.empty() )
510 {
511 // poly edge index, intersection point, edge index in boundary poly
512 std::vector<std::tuple<size_t, IfcVector3, size_t> > intersections;
513 bool startedInside = PointInPoly(proj * blackside.front(), profile->verts);
514 bool isCurrentlyInside = startedInside;
515
516 std::vector<std::pair<size_t, IfcVector3> > intersected_boundary;
517
518 for( size_t a = 0; a < blackside.size(); ++a )
519 {
520 const IfcVector3 e0 = proj * blackside[a];
521 const IfcVector3 e1 = proj * blackside[(a + 1) % blackside.size()];
522
523 intersected_boundary.clear();
524 IntersectsBoundaryProfile(e0, e1, profile->verts, isCurrentlyInside, intersected_boundary);
525 // sort the hits by distance from e0 to get the correct in/out/in sequence. Manually :-( I miss you, C++11.
526 if( intersected_boundary.size() > 1 )
527 {
528 bool keepSorting = true;
529 while( keepSorting )
530 {
531 keepSorting = false;
532 for( size_t b = 0; b < intersected_boundary.size() - 1; ++b )
533 {
534 if( (intersected_boundary[b + 1].second - e0).SquareLength() < (intersected_boundary[b].second - e0).SquareLength() )
535 {
536 keepSorting = true;
537 std::swap(intersected_boundary[b + 1], intersected_boundary[b]);
538 }
539 }
540 }
541 }
542 // now add them to the list of intersections
543 for( size_t b = 0; b < intersected_boundary.size(); ++b )
544 intersections.push_back(std::make_tuple(a, proj_inv * intersected_boundary[b].second, intersected_boundary[b].first));
545
546 // and calculate our new inside/outside state
547 if( intersected_boundary.size() & 1 )
548 isCurrentlyInside = !isCurrentlyInside;
549 }
550
551 // we got a list of in-out-combinations of intersections. That should be an even number of intersections, or
552 // we're fucked.
553 if( (intersections.size() & 1) != 0 )
554 {
555 IFCImporter::LogWarn("Odd number of intersections, can't work with that. Omitting half space boundary check.");
556 continue;
557 }
558
559 if( intersections.size() > 1 )
560 {
561 // If we started outside, the first intersection is a out->in intersection. Cycle them so that it
562 // starts with an intersection leaving the boundary
563 if( !startedInside )
564 for( size_t b = 0; b < intersections.size() - 1; ++b )
565 std::swap(intersections[b], intersections[(b + intersections.size() - 1) % intersections.size()]);
566
567 // Filter pairs of out->in->out that lie too close to each other.
568 for( size_t a = 0; intersections.size() > 0 && a < intersections.size() - 1; /**/ )
569 {
570 if( (std::get<1>(intersections[a]) - std::get<1>(intersections[(a + 1) % intersections.size()])).SquareLength() < 1e-10 )
571 intersections.erase(intersections.begin() + a, intersections.begin() + a + 2);
572 else
573 a++;
574 }
575 if( intersections.size() > 1 && (std::get<1>(intersections.back()) - std::get<1>(intersections.front())).SquareLength() < 1e-10 )
576 {
577 intersections.pop_back(); intersections.erase(intersections.begin());
578 }
579 }
580
581
582 // no intersections at all: either completely inside the boundary, so everything gets discarded, or completely outside.
583 // in the latter case we're implementional lost. I'm simply going to ignore this, so a large poly will not get any
584 // holes if the boundary is smaller and does not touch it anywhere.
585 if( intersections.empty() )
586 {
587 // starting point was outside -> everything is outside the boundary -> nothing is clipped -> add black side
588 // to result mesh unchanged
589 if( !startedInside )
590 {
591 outvertcnt.push_back(static_cast<unsigned int>(blackside.size()));
592 outvert.insert(outvert.end(), blackside.begin(), blackside.end());
593 continue;
594 }
595 else
596 {
597 // starting point was inside the boundary -> everything is inside the boundary -> nothing is spared from the
598 // clipping -> nothing left to add to the result mesh
599 continue;
600 }
601 }
602
603 // determine the direction in which we're marching along the boundary polygon. If the src poly is faced upwards
604 // and the boundary is also winded this way, we need to march *backwards* on the boundary.
605 const IfcVector3 polyNormal = IfcMatrix3(proj) * TempMesh::ComputePolygonNormal(blackside.data(), blackside.size());
606 bool marchBackwardsOnBoundary = (profileNormal * polyNormal) >= 0.0;
607
608 // Build closed loops from these intersections. Starting from an intersection leaving the boundary we
609 // walk along the polygon to the next intersection (which should be an IS entering the boundary poly).
610 // From there we walk along the boundary until we hit another intersection leaving the boundary,
611 // walk along the poly to the next IS and so on until we're back at the starting point.
612 // We remove every intersection we "used up", so any remaining intersection is the start of a new loop.
613 while( !intersections.empty() )
614 {
615 std::vector<IfcVector3> resultpoly;
616 size_t currentIntersecIdx = 0;
617
618 while( true )
619 {
620 ai_assert(intersections.size() > currentIntersecIdx + 1);
621 std::tuple<size_t, IfcVector3, size_t> currintsec = intersections[currentIntersecIdx + 0];
622 std::tuple<size_t, IfcVector3, size_t> nextintsec = intersections[currentIntersecIdx + 1];
623 intersections.erase(intersections.begin() + currentIntersecIdx, intersections.begin() + currentIntersecIdx + 2);
624
625 // we start with an in->out intersection
626 resultpoly.push_back(std::get<1>(currintsec));
627 // climb along the polygon to the next intersection, which should be an out->in
628 size_t numPolyPoints = (std::get<0>(currintsec) > std::get<0>(nextintsec) ? blackside.size() : 0)
629 + std::get<0>(nextintsec) - std::get<0>(currintsec);
630 for( size_t a = 1; a <= numPolyPoints; ++a )
631 resultpoly.push_back(blackside[(std::get<0>(currintsec) + a) % blackside.size()]);
632 // put the out->in intersection
633 resultpoly.push_back(std::get<1>(nextintsec));
634
635 // generate segments along the boundary polygon that lie in the poly's plane until we hit another intersection
636 IfcVector3 startingPoint = proj * std::get<1>(nextintsec);
637 size_t currentBoundaryEdgeIdx = (std::get<2>(nextintsec) + (marchBackwardsOnBoundary ? 1 : 0)) % profile->verts.size();
638 size_t nextIntsecIdx = SIZE_MAX;
639 while( nextIntsecIdx == SIZE_MAX )
640 {
641 IfcFloat t = 1e10;
642
643 size_t nextBoundaryEdgeIdx = marchBackwardsOnBoundary ? (currentBoundaryEdgeIdx + profile->verts.size() - 1) : currentBoundaryEdgeIdx + 1;
644 nextBoundaryEdgeIdx %= profile->verts.size();
645 // vertices of the current boundary segments
646 IfcVector3 currBoundaryPoint = profile->verts[currentBoundaryEdgeIdx];
647 IfcVector3 nextBoundaryPoint = profile->verts[nextBoundaryEdgeIdx];
648 // project the two onto the polygon
649 if( std::abs(polyNormal.z) > 1e-5 )
650 {
651 currBoundaryPoint.z = startingPoint.z + (currBoundaryPoint.x - startingPoint.x) * polyNormal.x/polyNormal.z + (currBoundaryPoint.y - startingPoint.y) * polyNormal.y/polyNormal.z;
652 nextBoundaryPoint.z = startingPoint.z + (nextBoundaryPoint.x - startingPoint.x) * polyNormal.x/polyNormal.z + (nextBoundaryPoint.y - startingPoint.y) * polyNormal.y/polyNormal.z;
653 }
654
655 // build a direction that goes along the boundary border but lies in the poly plane
656 IfcVector3 boundaryPlaneNormal = ((nextBoundaryPoint - currBoundaryPoint) ^ profileNormal).Normalize();
657 IfcVector3 dirAtPolyPlane = (boundaryPlaneNormal ^ polyNormal).Normalize() * (marchBackwardsOnBoundary ? -1.0 : 1.0);
658 // if we can project the direction to the plane, we can calculate a maximum marching distance along that dir
659 // until we finish that boundary segment and continue on the next
660 if( std::abs(polyNormal.z) > 1e-5 )
661 {
662 t = std::min(t, (nextBoundaryPoint - startingPoint).Length());
663 }
664
665 // check if the direction hits the loop start - if yes, we got a poly to output
666 IfcVector3 dirToThatPoint = proj * resultpoly.front() - startingPoint;
667 IfcFloat tpt = dirToThatPoint * dirAtPolyPlane;
668 if( tpt > -1e-6 && tpt <= t && (dirToThatPoint - tpt * dirAtPolyPlane).SquareLength() < 1e-10 )
669 {
670 nextIntsecIdx = intersections.size(); // dirty hack to end marching along the boundary and signal the end of the loop
671 t = tpt;
672 }
673
674 // also check if the direction hits any in->out intersections earlier. If we hit one, we can switch back
675 // to marching along the poly border from that intersection point
676 for( size_t a = 0; a < intersections.size(); a += 2 )
677 {
678 dirToThatPoint = proj * std::get<1>(intersections[a]) - startingPoint;
679 tpt = dirToThatPoint * dirAtPolyPlane;
680 if( tpt > -1e-6 && tpt <= t && (dirToThatPoint - tpt * dirAtPolyPlane).SquareLength() < 1e-10 )
681 {
682 nextIntsecIdx = a; // switch back to poly and march on from this in->out intersection
683 t = tpt;
684 }
685 }
686
687 // if we keep marching on the boundary, put the segment end point to the result poly and well... keep marching
688 if( nextIntsecIdx == SIZE_MAX )
689 {
690 resultpoly.push_back(proj_inv * nextBoundaryPoint);
691 currentBoundaryEdgeIdx = nextBoundaryEdgeIdx;
692 startingPoint = nextBoundaryPoint;
693 }
694
695 // quick endless loop check
696 if( resultpoly.size() > blackside.size() + profile->verts.size() )
697 {
698 IFCImporter::LogError("Encountered endless loop while clipping polygon against poly-bounded half space.");
699 break;
700 }
701 }
702
703 // we're back on the poly - if this is the intersection we started from, we got a closed loop.
704 if( nextIntsecIdx >= intersections.size() )
705 {
706 break;
707 }
708
709 // otherwise it's another intersection. Continue marching from there.
710 currentIntersecIdx = nextIntsecIdx;
711 }
712
713 WritePolygon(resultpoly, result);
714 }
715 }
716 }
717 IFCImporter::LogDebug("generating CSG geometry by plane clipping with polygonal bounding (IfcBooleanClippingResult)");
718}
719
720// ------------------------------------------------------------------------------------------------
721void ProcessBooleanExtrudedAreaSolidDifference(const IfcExtrudedAreaSolid* as, TempMesh& result,
722 const TempMesh& first_operand,
723 ConversionData& conv)
724{
725 ai_assert(as != NULL);
726
727 // This case is handled by reduction to an instance of the quadrify() algorithm.
728 // Obviously, this won't work for arbitrarily complex cases. In fact, the first
729 // operand should be near-planar. Luckily, this is usually the case in Ifc
730 // buildings.
731
732 std::shared_ptr<TempMesh> meshtmp = std::shared_ptr<TempMesh>(new TempMesh());
733 ProcessExtrudedAreaSolid(*as,*meshtmp,conv,false);
734
735 std::vector<TempOpening> openings(1, TempOpening(as,IfcVector3(0,0,0),meshtmp,std::shared_ptr<TempMesh>()));
736
737 result = first_operand;
738
739 TempMesh temp;
740
741 std::vector<IfcVector3>::const_iterator vit = first_operand.verts.begin();
742 for(unsigned int pcount : first_operand.vertcnt) {
743 temp.Clear();
744
745 temp.verts.insert(temp.verts.end(), vit, vit + pcount);
746 temp.vertcnt.push_back(pcount);
747
748 // The algorithms used to generate mesh geometry sometimes
749 // spit out lines or other degenerates which must be
750 // filtered to avoid running into assertions later on.
751
752 // ComputePolygonNormal returns the Newell normal, so the
753 // length of the normal is the area of the polygon.
754 const IfcVector3& normal = temp.ComputeLastPolygonNormal(false);
755 if (normal.SquareLength() < static_cast<IfcFloat>(1e-5)) {
756 IFCImporter::LogWarn("skipping degenerate polygon (ProcessBooleanExtrudedAreaSolidDifference)");
757 continue;
758 }
759
760 GenerateOpenings(openings, std::vector<IfcVector3>(1,IfcVector3(1,0,0)), temp, false, true);
761 result.Append(temp);
762
763 vit += pcount;
764 }
765
766 IFCImporter::LogDebug("generating CSG geometry by geometric difference to a solid (IfcExtrudedAreaSolid)");
767}
768
769// ------------------------------------------------------------------------------------------------
770void ProcessBoolean(const IfcBooleanResult& boolean, TempMesh& result, ConversionData& conv)
771{
772 // supported CSG operations:
773 // DIFFERENCE
774 if(const IfcBooleanResult* const clip = boolean.ToPtr<IfcBooleanResult>()) {
775 if(clip->Operator != "DIFFERENCE") {
776 IFCImporter::LogWarn("encountered unsupported boolean operator: " + (std::string)clip->Operator);
777 return;
778 }
779
780 // supported cases (1st operand):
781 // IfcBooleanResult -- call ProcessBoolean recursively
782 // IfcSweptAreaSolid -- obtain polygonal geometry first
783
784 // supported cases (2nd operand):
785 // IfcHalfSpaceSolid -- easy, clip against plane
786 // IfcExtrudedAreaSolid -- reduce to an instance of the quadrify() algorithm
787
788
789 const IfcHalfSpaceSolid* const hs = clip->SecondOperand->ResolveSelectPtr<IfcHalfSpaceSolid>(conv.db);
790 const IfcExtrudedAreaSolid* const as = clip->SecondOperand->ResolveSelectPtr<IfcExtrudedAreaSolid>(conv.db);
791 if(!hs && !as) {
792 IFCImporter::LogError("expected IfcHalfSpaceSolid or IfcExtrudedAreaSolid as second clipping operand");
793 return;
794 }
795
796 TempMesh first_operand;
797 if(const IfcBooleanResult* const op0 = clip->FirstOperand->ResolveSelectPtr<IfcBooleanResult>(conv.db)) {
798 ProcessBoolean(*op0,first_operand,conv);
799 }
800 else if (const IfcSweptAreaSolid* const swept = clip->FirstOperand->ResolveSelectPtr<IfcSweptAreaSolid>(conv.db)) {
801 ProcessSweptAreaSolid(*swept,first_operand,conv);
802 }
803 else {
804 IFCImporter::LogError("expected IfcSweptAreaSolid or IfcBooleanResult as first clipping operand");
805 return;
806 }
807
808 if(hs) {
809
810 const IfcPolygonalBoundedHalfSpace* const hs_bounded = clip->SecondOperand->ResolveSelectPtr<IfcPolygonalBoundedHalfSpace>(conv.db);
811 if (hs_bounded) {
812 ProcessPolygonalBoundedBooleanHalfSpaceDifference(hs_bounded, result, first_operand, conv);
813 }
814 else {
815 ProcessBooleanHalfSpaceDifference(hs, result, first_operand, conv);
816 }
817 }
818 else {
819 ProcessBooleanExtrudedAreaSolidDifference(as, result, first_operand, conv);
820 }
821 }
822 else {
823 IFCImporter::LogWarn("skipping unknown IfcBooleanResult entity, type is " + boolean.GetClassName());
824 }
825}
826
827} // ! IFC
828} // ! Assimp
829
830#endif
831
832