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 <assimp/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. |

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 | // ------------------------------------------------------------------------------------------------ |

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(), 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 | // ------------------------------------------------------------------------------------------------ |

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(static_cast<unsigned int>(resultpoly.size())); |

136 | } |

137 | } |

138 | |

139 | |

140 | // ------------------------------------------------------------------------------------------------ |

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? |

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 | |

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. |

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 | return votes > 1; |

385 | } |

386 | |

387 | |

388 | // ------------------------------------------------------------------------------------------------ |

389 | void 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 | // ------------------------------------------------------------------------------------------------ |

721 | void 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 | // ------------------------------------------------------------------------------------------------ |

770 | void 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 |