1 /*
2 Open Asset Import Library (assimp)
3 ----------------------------------------------------------------------
4
5 Copyright (c) 2006-2010, assimp team
6 All rights reserved.
7
8 Redistribution and use of this software in source and binary forms,
9 with or without modification, are permitted provided that the
10 following 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
26 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
27 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
29 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
30 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
31 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
32 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
33 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
34 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
35 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
36 OF 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 "Defines.h"
51
52 #include <iterator>
53 #include <tuple>
54
55
56 namespace 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.
IntersectSegmentPlane(const IfcVector3 & p,const IfcVector3 & n,const IfcVector3 & e0,const IfcVector3 & e1,bool assumeStartOnWhiteSide,IfcVector3 & out)63 bool 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 // ------------------------------------------------------------------------------------------------
FilterPolygon(std::vector<IfcVector3> & resultpoly)103 void 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(), 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 // ------------------------------------------------------------------------------------------------
WritePolygon(std::vector<IfcVector3> & resultpoly,TempMesh & result)128 void 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(resultpoly.size());
136 }
137 }
138
139
140 // ------------------------------------------------------------------------------------------------
ProcessBooleanHalfSpaceDifference(const IfcHalfSpaceSolid * hs,TempMesh & result,const TempMesh & first_operand,ConversionData &)141 void 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?
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)249 bool 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 IfcFloat b_sqlen_inv = 1.0 / b.SquareLength();
276
277 // segment-segment intersection
278 // solve b0 + b*s = e0 + e*t for (s,t)
279 const IfcFloat det = (-b.x * e.y + e.x * b.y);
280 if( std::abs(det) < 1e-6 ) {
281 // no solutions (parallel lines)
282 continue;
283 }
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.
PointInPoly(const IfcVector3 & p,const std::vector<IfcVector3> & boundary)358 bool 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 // ai_assert(votes == 3 || votes == 0);
385 return votes > 1;
386 }
387
388
389 // ------------------------------------------------------------------------------------------------
ProcessPolygonalBoundedBooleanHalfSpaceDifference(const IfcPolygonalBoundedHalfSpace * hs,TempMesh & result,const TempMesh & first_operand,ConversionData & conv)390 void ProcessPolygonalBoundedBooleanHalfSpaceDifference(const IfcPolygonalBoundedHalfSpace* hs, TempMesh& result,
391 const TempMesh& first_operand,
392 ConversionData& conv)
393 {
394 ai_assert(hs != NULL);
395
396 const IfcPlane* const plane = hs->BaseSurface->ToPtr<IfcPlane>();
397 if(!plane) {
398 IFCImporter::LogError("expected IfcPlane as base surface for the IfcHalfSpaceSolid");
399 return;
400 }
401
402 // extract plane base position vector and normal vector
403 IfcVector3 p,n(0.f,0.f,1.f);
404 if (plane->Position->Axis) {
405 ConvertDirection(n,plane->Position->Axis.Get());
406 }
407 ConvertCartesianPoint(p,plane->Position->Location);
408
409 if(!IsTrue(hs->AgreementFlag)) {
410 n *= -1.f;
411 }
412
413 n.Normalize();
414
415 // obtain the polygonal bounding volume
416 std::shared_ptr<TempMesh> profile = std::shared_ptr<TempMesh>(new TempMesh());
417 if(!ProcessCurve(hs->PolygonalBoundary, *profile.get(), conv)) {
418 IFCImporter::LogError("expected valid polyline for boundary of boolean halfspace");
419 return;
420 }
421
422 // determine winding order by calculating the normal.
423 IfcVector3 profileNormal = TempMesh::ComputePolygonNormal(profile->verts.data(), profile->verts.size());
424
425 IfcMatrix4 proj_inv;
426 ConvertAxisPlacement(proj_inv,hs->Position);
427
428 // and map everything into a plane coordinate space so all intersection
429 // tests can be done in 2D space.
430 IfcMatrix4 proj = proj_inv;
431 proj.Inverse();
432
433 // clip the current contents of `meshout` against the plane we obtained from the second operand
434 const std::vector<IfcVector3>& in = first_operand.verts;
435 std::vector<IfcVector3>& outvert = result.verts;
436 std::vector<unsigned int>& outvertcnt = result.vertcnt;
437
438 outvert.reserve(in.size());
439 outvertcnt.reserve(first_operand.vertcnt.size());
440
441 unsigned int vidx = 0;
442 std::vector<unsigned int>::const_iterator begin = first_operand.vertcnt.begin();
443 std::vector<unsigned int>::const_iterator end = first_operand.vertcnt.end();
444 std::vector<unsigned int>::const_iterator iit;
445 for( iit = begin; iit != end; vidx += *iit++ )
446 {
447 // Our new approach: we cut the poly along the plane, then we intersect the part on the black side of the plane
448 // against the bounding polygon. All the white parts, and the black part outside the boundary polygon, are kept.
449 std::vector<IfcVector3> whiteside, blackside;
450
451 {
452 const IfcVector3* srcVertices = &in[vidx];
453 const size_t srcVtxCount = *iit;
454 if( srcVtxCount == 0 )
455 continue;
456
457 IfcVector3 polyNormal = TempMesh::ComputePolygonNormal(srcVertices, srcVtxCount, true);
458
459 // if the poly is parallel to the plane, put it completely on the black or white side
460 if( std::abs(polyNormal * n) > 0.9999 )
461 {
462 bool isOnWhiteSide = (srcVertices[0] - p) * n > -1e-6;
463 std::vector<IfcVector3>& targetSide = isOnWhiteSide ? whiteside : blackside;
464 targetSide.insert(targetSide.end(), srcVertices, srcVertices + srcVtxCount);
465 }
466 else
467 {
468 // otherwise start building one polygon for each side. Whenever the current line segment intersects the plane
469 // we put a point there as an end of the current segment. Then we switch to the other side, put a point there, too,
470 // as a beginning of the current segment, and simply continue accumulating vertices.
471 bool isCurrentlyOnWhiteSide = ((srcVertices[0]) - p) * n > -1e-6;
472 for( size_t a = 0; a < srcVtxCount; ++a )
473 {
474 IfcVector3 e0 = srcVertices[a];
475 IfcVector3 e1 = srcVertices[(a + 1) % srcVtxCount];
476 IfcVector3 ei;
477
478 // put starting point to the current mesh
479 std::vector<IfcVector3>& trgt = isCurrentlyOnWhiteSide ? whiteside : blackside;
480 trgt.push_back(srcVertices[a]);
481
482 // if there's an intersection, put an end vertex there, switch to the other side's mesh,
483 // and add a starting vertex there, too
484 bool isPlaneHit = IntersectSegmentPlane(p, n, e0, e1, isCurrentlyOnWhiteSide, ei);
485 if( isPlaneHit )
486 {
487 if( trgt.empty() || (trgt.back() - ei).SquareLength() > 1e-12 )
488 trgt.push_back(ei);
489 isCurrentlyOnWhiteSide = !isCurrentlyOnWhiteSide;
490 std::vector<IfcVector3>& newtrgt = isCurrentlyOnWhiteSide ? whiteside : blackside;
491 newtrgt.push_back(ei);
492 }
493 }
494 }
495 }
496
497 // the part on the white side can be written into the target mesh right away
498 WritePolygon(whiteside, result);
499
500 // The black part is the piece we need to get rid of, but only the part of it within the boundary polygon.
501 // So we now need to construct all the polygons that result from BlackSidePoly minus BoundaryPoly.
502 FilterPolygon(blackside);
503
504 // Complicated, II. We run along the polygon. a) When we're inside the boundary, we run on until we hit an
505 // intersection, which means we're leaving it. We then start a new out poly there. b) When we're outside the
506 // boundary, we start collecting vertices until we hit an intersection, then we run along the boundary until we hit
507 // an intersection, then we switch back to the poly and run on on this one again, and so on until we got a closed
508 // loop. Then we continue with the path we left to catch potential additional polys on the other side of the
509 // boundary as described in a)
510 if( !blackside.empty() )
511 {
512 // poly edge index, intersection point, edge index in boundary poly
513 std::vector<std::tuple<size_t, IfcVector3, size_t> > intersections;
514 bool startedInside = PointInPoly(proj * blackside.front(), profile->verts);
515 bool isCurrentlyInside = startedInside;
516
517 std::vector<std::pair<size_t, IfcVector3> > intersected_boundary;
518
519 for( size_t a = 0; a < blackside.size(); ++a )
520 {
521 const IfcVector3 e0 = proj * blackside[a];
522 const IfcVector3 e1 = proj * blackside[(a + 1) % blackside.size()];
523
524 intersected_boundary.clear();
525 IntersectsBoundaryProfile(e0, e1, profile->verts, isCurrentlyInside, intersected_boundary);
526 // sort the hits by distance from e0 to get the correct in/out/in sequence. Manually :-( I miss you, C++11.
527 if( intersected_boundary.size() > 1 )
528 {
529 bool keepSorting = true;
530 while( keepSorting )
531 {
532 keepSorting = false;
533 for( size_t b = 0; b < intersected_boundary.size() - 1; ++b )
534 {
535 if( (intersected_boundary[b + 1].second - e0).SquareLength() < (intersected_boundary[b].second - e0).SquareLength() )
536 {
537 keepSorting = true;
538 std::swap(intersected_boundary[b + 1], intersected_boundary[b]);
539 }
540 }
541 }
542 }
543 // now add them to the list of intersections
544 for( size_t b = 0; b < intersected_boundary.size(); ++b )
545 intersections.push_back(std::make_tuple(a, proj_inv * intersected_boundary[b].second, intersected_boundary[b].first));
546
547 // and calculate our new inside/outside state
548 if( intersected_boundary.size() & 1 )
549 isCurrentlyInside = !isCurrentlyInside;
550 }
551
552 // we got a list of in-out-combinations of intersections. That should be an even number of intersections, or
553 // we're fucked.
554 if( (intersections.size() & 1) != 0 )
555 {
556 IFCImporter::LogWarn("Odd number of intersections, can't work with that. Omitting half space boundary check.");
557 continue;
558 }
559
560 if( intersections.size() > 1 )
561 {
562 // If we started outside, the first intersection is a out->in intersection. Cycle them so that it
563 // starts with an intersection leaving the boundary
564 if( !startedInside )
565 for( size_t b = 0; b < intersections.size() - 1; ++b )
566 std::swap(intersections[b], intersections[(b + intersections.size() - 1) % intersections.size()]);
567
568 // Filter pairs of out->in->out that lie too close to each other.
569 for( size_t a = 0; intersections.size() > 0 && a < intersections.size() - 1; /**/ )
570 {
571 if( (std::get<1>(intersections[a]) - std::get<1>(intersections[(a + 1) % intersections.size()])).SquareLength() < 1e-10 )
572 intersections.erase(intersections.begin() + a, intersections.begin() + a + 2);
573 else
574 a++;
575 }
576 if( intersections.size() > 1 && (std::get<1>(intersections.back()) - std::get<1>(intersections.front())).SquareLength() < 1e-10 )
577 {
578 intersections.pop_back(); intersections.erase(intersections.begin());
579 }
580 }
581
582
583 // no intersections at all: either completely inside the boundary, so everything gets discarded, or completely outside.
584 // in the latter case we're implementional lost. I'm simply going to ignore this, so a large poly will not get any
585 // holes if the boundary is smaller and does not touch it anywhere.
586 if( intersections.empty() )
587 {
588 // starting point was outside -> everything is outside the boundary -> nothing is clipped -> add black side
589 // to result mesh unchanged
590 if( !startedInside )
591 {
592 outvertcnt.push_back(blackside.size());
593 outvert.insert(outvert.end(), blackside.begin(), blackside.end());
594 continue;
595 }
596 else
597 {
598 // starting point was inside the boundary -> everything is inside the boundary -> nothing is spared from the
599 // clipping -> nothing left to add to the result mesh
600 continue;
601 }
602 }
603
604 // determine the direction in which we're marching along the boundary polygon. If the src poly is faced upwards
605 // and the boundary is also winded this way, we need to march *backwards* on the boundary.
606 const IfcVector3 polyNormal = IfcMatrix3(proj) * TempMesh::ComputePolygonNormal(blackside.data(), blackside.size());
607 bool marchBackwardsOnBoundary = (profileNormal * polyNormal) >= 0.0;
608
609 // Build closed loops from these intersections. Starting from an intersection leaving the boundary we
610 // walk along the polygon to the next intersection (which should be an IS entering the boundary poly).
611 // From there we walk along the boundary until we hit another intersection leaving the boundary,
612 // walk along the poly to the next IS and so on until we're back at the starting point.
613 // We remove every intersection we "used up", so any remaining intersection is the start of a new loop.
614 while( !intersections.empty() )
615 {
616 std::vector<IfcVector3> resultpoly;
617 size_t currentIntersecIdx = 0;
618
619 while( true )
620 {
621 ai_assert(intersections.size() > currentIntersecIdx + 1);
622 std::tuple<size_t, IfcVector3, size_t> currintsec = intersections[currentIntersecIdx + 0];
623 std::tuple<size_t, IfcVector3, size_t> nextintsec = intersections[currentIntersecIdx + 1];
624 intersections.erase(intersections.begin() + currentIntersecIdx, intersections.begin() + currentIntersecIdx + 2);
625
626 // we start with an in->out intersection
627 resultpoly.push_back(std::get<1>(currintsec));
628 // climb along the polygon to the next intersection, which should be an out->in
629 size_t numPolyPoints = (std::get<0>(currintsec) > std::get<0>(nextintsec) ? blackside.size() : 0)
630 + std::get<0>(nextintsec) - std::get<0>(currintsec);
631 for( size_t a = 1; a <= numPolyPoints; ++a )
632 resultpoly.push_back(blackside[(std::get<0>(currintsec) + a) % blackside.size()]);
633 // put the out->in intersection
634 resultpoly.push_back(std::get<1>(nextintsec));
635
636 // generate segments along the boundary polygon that lie in the poly's plane until we hit another intersection
637 IfcVector3 startingPoint = proj * std::get<1>(nextintsec);
638 size_t currentBoundaryEdgeIdx = (std::get<2>(nextintsec) + (marchBackwardsOnBoundary ? 1 : 0)) % profile->verts.size();
639 size_t nextIntsecIdx = SIZE_MAX;
640 while( nextIntsecIdx == SIZE_MAX )
641 {
642 IfcFloat t = 1e10;
643
644 size_t nextBoundaryEdgeIdx = marchBackwardsOnBoundary ? (currentBoundaryEdgeIdx + profile->verts.size() - 1) : currentBoundaryEdgeIdx + 1;
645 nextBoundaryEdgeIdx %= profile->verts.size();
646 // vertices of the current boundary segments
647 IfcVector3 currBoundaryPoint = profile->verts[currentBoundaryEdgeIdx];
648 IfcVector3 nextBoundaryPoint = profile->verts[nextBoundaryEdgeIdx];
649 // project the two onto the polygon
650 if( std::abs(polyNormal.z) > 1e-5 )
651 {
652 currBoundaryPoint.z = startingPoint.z + (currBoundaryPoint.x - startingPoint.x) * polyNormal.x/polyNormal.z + (currBoundaryPoint.y - startingPoint.y) * polyNormal.y/polyNormal.z;
653 nextBoundaryPoint.z = startingPoint.z + (nextBoundaryPoint.x - startingPoint.x) * polyNormal.x/polyNormal.z + (nextBoundaryPoint.y - startingPoint.y) * polyNormal.y/polyNormal.z;
654 }
655
656 // build a direction that goes along the boundary border but lies in the poly plane
657 IfcVector3 boundaryPlaneNormal = ((nextBoundaryPoint - currBoundaryPoint) ^ profileNormal).Normalize();
658 IfcVector3 dirAtPolyPlane = (boundaryPlaneNormal ^ polyNormal).Normalize() * (marchBackwardsOnBoundary ? -1.0 : 1.0);
659 // if we can project the direction to the plane, we can calculate a maximum marching distance along that dir
660 // until we finish that boundary segment and continue on the next
661 if( std::abs(polyNormal.z) > 1e-5 )
662 {
663 t = std::min(t, (nextBoundaryPoint - startingPoint).Length());
664 }
665
666 // check if the direction hits the loop start - if yes, we got a poly to output
667 IfcVector3 dirToThatPoint = proj * resultpoly.front() - startingPoint;
668 IfcFloat tpt = dirToThatPoint * dirAtPolyPlane;
669 if( tpt > -1e-6 && tpt <= t && (dirToThatPoint - tpt * dirAtPolyPlane).SquareLength() < 1e-10 )
670 {
671 nextIntsecIdx = intersections.size(); // dirty hack to end marching along the boundary and signal the end of the loop
672 t = tpt;
673 }
674
675 // also check if the direction hits any in->out intersections earlier. If we hit one, we can switch back
676 // to marching along the poly border from that intersection point
677 for( size_t a = 0; a < intersections.size(); a += 2 )
678 {
679 dirToThatPoint = proj * std::get<1>(intersections[a]) - startingPoint;
680 tpt = dirToThatPoint * dirAtPolyPlane;
681 if( tpt > -1e-6 && tpt <= t && (dirToThatPoint - tpt * dirAtPolyPlane).SquareLength() < 1e-10 )
682 {
683 nextIntsecIdx = a; // switch back to poly and march on from this in->out intersection
684 t = tpt;
685 }
686 }
687
688 // if we keep marching on the boundary, put the segment end point to the result poly and well... keep marching
689 if( nextIntsecIdx == SIZE_MAX )
690 {
691 resultpoly.push_back(proj_inv * nextBoundaryPoint);
692 currentBoundaryEdgeIdx = nextBoundaryEdgeIdx;
693 startingPoint = nextBoundaryPoint;
694 }
695
696 // quick endless loop check
697 if( resultpoly.size() > blackside.size() + profile->verts.size() )
698 {
699 IFCImporter::LogError("Encountered endless loop while clipping polygon against poly-bounded half space.");
700 break;
701 }
702 }
703
704 // we're back on the poly - if this is the intersection we started from, we got a closed loop.
705 if( nextIntsecIdx >= intersections.size() )
706 {
707 break;
708 }
709
710 // otherwise it's another intersection. Continue marching from there.
711 currentIntersecIdx = nextIntsecIdx;
712 }
713
714 WritePolygon(resultpoly, result);
715 }
716 }
717 }
718 IFCImporter::LogDebug("generating CSG geometry by plane clipping with polygonal bounding (IfcBooleanClippingResult)");
719 }
720
721 // ------------------------------------------------------------------------------------------------
ProcessBooleanExtrudedAreaSolidDifference(const IfcExtrudedAreaSolid * as,TempMesh & result,const TempMesh & first_operand,ConversionData & conv)722 void ProcessBooleanExtrudedAreaSolidDifference(const IfcExtrudedAreaSolid* as, TempMesh& result,
723 const TempMesh& first_operand,
724 ConversionData& conv)
725 {
726 ai_assert(as != NULL);
727
728 // This case is handled by reduction to an instance of the quadrify() algorithm.
729 // Obviously, this won't work for arbitrarily complex cases. In fact, the first
730 // operand should be near-planar. Luckily, this is usually the case in Ifc
731 // buildings.
732
733 std::shared_ptr<TempMesh> meshtmp = std::shared_ptr<TempMesh>(new TempMesh());
734 ProcessExtrudedAreaSolid(*as,*meshtmp,conv,false);
735
736 std::vector<TempOpening> openings(1, TempOpening(as,IfcVector3(0,0,0),meshtmp,std::shared_ptr<TempMesh>()));
737
738 result = first_operand;
739
740 TempMesh temp;
741
742 std::vector<IfcVector3>::const_iterator vit = first_operand.verts.begin();
743 for(unsigned int pcount : first_operand.vertcnt) {
744 temp.Clear();
745
746 temp.verts.insert(temp.verts.end(), vit, vit + pcount);
747 temp.vertcnt.push_back(pcount);
748
749 // The algorithms used to generate mesh geometry sometimes
750 // spit out lines or other degenerates which must be
751 // filtered to avoid running into assertions later on.
752
753 // ComputePolygonNormal returns the Newell normal, so the
754 // length of the normal is the area of the polygon.
755 const IfcVector3& normal = temp.ComputeLastPolygonNormal(false);
756 if (normal.SquareLength() < static_cast<IfcFloat>(1e-5)) {
757 IFCImporter::LogWarn("skipping degenerate polygon (ProcessBooleanExtrudedAreaSolidDifference)");
758 continue;
759 }
760
761 GenerateOpenings(openings, std::vector<IfcVector3>(1,IfcVector3(1,0,0)), temp, false, true);
762 result.Append(temp);
763
764 vit += pcount;
765 }
766
767 IFCImporter::LogDebug("generating CSG geometry by geometric difference to a solid (IfcExtrudedAreaSolid)");
768 }
769
770 // ------------------------------------------------------------------------------------------------
ProcessBoolean(const IfcBooleanResult & boolean,TempMesh & result,ConversionData & conv)771 void ProcessBoolean(const IfcBooleanResult& boolean, TempMesh& result, ConversionData& conv)
772 {
773 // supported CSG operations:
774 // DIFFERENCE
775 if(const IfcBooleanResult* const clip = boolean.ToPtr<IfcBooleanResult>()) {
776 if(clip->Operator != "DIFFERENCE") {
777 IFCImporter::LogWarn("encountered unsupported boolean operator: " + (std::string)clip->Operator);
778 return;
779 }
780
781 // supported cases (1st operand):
782 // IfcBooleanResult -- call ProcessBoolean recursively
783 // IfcSweptAreaSolid -- obtain polygonal geometry first
784
785 // supported cases (2nd operand):
786 // IfcHalfSpaceSolid -- easy, clip against plane
787 // IfcExtrudedAreaSolid -- reduce to an instance of the quadrify() algorithm
788
789
790 const IfcHalfSpaceSolid* const hs = clip->SecondOperand->ResolveSelectPtr<IfcHalfSpaceSolid>(conv.db);
791 const IfcExtrudedAreaSolid* const as = clip->SecondOperand->ResolveSelectPtr<IfcExtrudedAreaSolid>(conv.db);
792 if(!hs && !as) {
793 IFCImporter::LogError("expected IfcHalfSpaceSolid or IfcExtrudedAreaSolid as second clipping operand");
794 return;
795 }
796
797 TempMesh first_operand;
798 if(const IfcBooleanResult* const op0 = clip->FirstOperand->ResolveSelectPtr<IfcBooleanResult>(conv.db)) {
799 ProcessBoolean(*op0,first_operand,conv);
800 }
801 else if (const IfcSweptAreaSolid* const swept = clip->FirstOperand->ResolveSelectPtr<IfcSweptAreaSolid>(conv.db)) {
802 ProcessSweptAreaSolid(*swept,first_operand,conv);
803 }
804 else {
805 IFCImporter::LogError("expected IfcSweptAreaSolid or IfcBooleanResult as first clipping operand");
806 return;
807 }
808
809 if(hs) {
810
811 const IfcPolygonalBoundedHalfSpace* const hs_bounded = clip->SecondOperand->ResolveSelectPtr<IfcPolygonalBoundedHalfSpace>(conv.db);
812 if (hs_bounded) {
813 ProcessPolygonalBoundedBooleanHalfSpaceDifference(hs_bounded, result, first_operand, conv);
814 }
815 else {
816 ProcessBooleanHalfSpaceDifference(hs, result, first_operand, conv);
817 }
818 }
819 else {
820 ProcessBooleanExtrudedAreaSolidDifference(as, result, first_operand, conv);
821 }
822 }
823 else {
824 IFCImporter::LogWarn("skipping unknown IfcBooleanResult entity, type is " + boolean.GetClassName());
825 }
826 }
827
828 } // ! IFC
829 } // ! Assimp
830
831 #endif
832
833