diff --git a/code/IFCUtil.cpp b/code/IFCUtil.cpp index e7a9e1b9c..483410ff3 100644 --- a/code/IFCUtil.cpp +++ b/code/IFCUtil.cpp @@ -570,6 +570,7 @@ void ConvertTransformOperator(IfcMatrix4& out, const IfcCartesianTransformationO out = locm * out * s; } + } // ! IFC } // ! Assimp diff --git a/code/IFCUtil.h b/code/IFCUtil.h index 9c434a490..9e2e15d79 100644 --- a/code/IFCUtil.h +++ b/code/IFCUtil.h @@ -299,6 +299,7 @@ bool GenerateOpenings(std::vector& openings, const IfcVector3& wall_extrusion_axis = IfcVector3(0,1,0)); + // IFCCurve.cpp // ------------------------------------------------------------------------------------------------ @@ -403,7 +404,6 @@ public: // IfcProfile.cpp bool ProcessCurve(const IfcCurve& curve, TempMesh& meshout, ConversionData& conv); - } } diff --git a/code/IfcBoolean.cpp b/code/IfcBoolean.cpp index 1d3226700..ce272ca1b 100644 --- a/code/IfcBoolean.cpp +++ b/code/IfcBoolean.cpp @@ -182,6 +182,108 @@ void ProcessBooleanHalfSpaceDifference(const IfcHalfSpaceSolid* hs, TempMesh& re IFCImporter::LogDebug("generating CSG geometry by plane clipping (IfcBooleanClippingResult)"); } +// ------------------------------------------------------------------------------------------------ +// Check if e0-e1 intersects a sub-segment of the given boundary line. +// note: this method works on 3D vectors, but performs its intersection checks solely in xy. +bool IntersectsBoundaryProfile( const IfcVector3& e0, const IfcVector3& e1, const std::vector& boundary, + std::vector& intersected_boundary_segments, + std::vector& intersected_boundary_points, + bool half_open = false) +{ + ai_assert(intersected_boundary_segments.empty()); + ai_assert(intersected_boundary_points.empty()); + + const IfcVector3& e = e1 - e0; + + for (size_t i = 0, bcount = boundary.size(); i < bcount; ++i) { + // boundary segment i: b0-b1 + const IfcVector3& b0 = boundary[i]; + const IfcVector3& b1 = boundary[(i+1) % bcount]; + + const IfcVector3& b = b1 - b0; + + // segment-segment intersection + // solve b0 + b*s = e0 + e*s for (s,t) + const IfcFloat det = (-b.x * e.y + e.x * b.y); + if(fabs(det) < 1e-6) { + // no solutions (parallel lines) + continue; + } + + const IfcFloat x = b0.x - e0.x; + const IfcFloat y = b0.y - e0.y; + + const IfcFloat s = (x*e.y - e.x*y)/det; + const IfcFloat t = (x*b.y - b.x*y)/det; + +#ifdef _DEBUG + const IfcVector3 check = b0 + b*s - (e0 + e*t); + ai_assert((IfcVector2(check.x,check.y)).SquareLength() < 1e-5); +#endif + + // for a valid intersection, s-t should be in range [0,1] + if (s >= 0.0 && (s <= 1.0 || half_open) && t >= 0.0 && t <= 1.0) { + + const IfcVector3& p = b0 + b*s; + + // only insert the point into the list if it is sufficiently + // far away from the previous intersection point. This way, + // we avoid duplicate detection if the intersection is + // directly on the vertex between two segments. + if (!intersected_boundary_points.empty() && intersected_boundary_segments.back()==(i==0?bcount-1:i-1) ) { + if((intersected_boundary_points.back() - p).SquareLength() < 1e-5) { + continue; + } + } + intersected_boundary_segments.push_back(i); + intersected_boundary_points.push_back(p); + } + } + + return false; +} + + +// ------------------------------------------------------------------------------------------------ +bool PointInPoly(const IfcVector3& p, const std::vector& boundary) +{ + // even-odd algorithm: take a random vector that extends from p to infinite + // and counts how many times it intersects edges of the boundary. + // because checking for segment intersections is prone to numeric inaccuracies + // or double detections (i.e. when hitting multiple adjacent segments at their + // shared vertices) we do it thrice with different rays and vote on it. + + std::vector intersected_boundary_segments; + std::vector intersected_boundary_points; + size_t votes = 0; + + IntersectsBoundaryProfile(p, p + IfcVector3(1.0,0,0), boundary, + intersected_boundary_segments, + intersected_boundary_points, true); + + votes += intersected_boundary_segments.size() % 2; + + intersected_boundary_segments.clear(); + intersected_boundary_points.clear(); + + IntersectsBoundaryProfile(p, p + IfcVector3(0,1.0,0), boundary, + intersected_boundary_segments, + intersected_boundary_points, true); + + votes += intersected_boundary_segments.size() % 2; + + intersected_boundary_segments.clear(); + intersected_boundary_points.clear(); + + IntersectsBoundaryProfile(p, p + IfcVector3(0,0,1.0), boundary, + intersected_boundary_segments, + intersected_boundary_points, true); + + votes += intersected_boundary_segments.size() % 2; + return votes > 1; +} + + // ------------------------------------------------------------------------------------------------ void ProcessPolygonalBoundedBooleanHalfSpaceDifference(const IfcPolygonalBoundedHalfSpace* hs, TempMesh& result, const TempMesh& first_operand, @@ -189,9 +291,249 @@ void ProcessPolygonalBoundedBooleanHalfSpaceDifference(const IfcPolygonalBounded { ai_assert(hs != NULL); - return; // niy + const IfcPlane* const plane = hs->BaseSurface->ToPtr(); + if(!plane) { + IFCImporter::LogError("expected IfcPlane as base surface for the IfcHalfSpaceSolid"); + return; + } + // extract plane base position vector and normal vector + IfcVector3 p,n(0.f,0.f,1.f); + if (plane->Position->Axis) { + ConvertDirection(n,plane->Position->Axis.Get()); + } + ConvertCartesianPoint(p,plane->Position->Location); + + if(!IsTrue(hs->AgreementFlag)) { + n *= -1.f; + } + + n.Normalize(); + + // obtain the polygonal bounding volume + boost::shared_ptr profile = boost::shared_ptr(new TempMesh()); + if(!ProcessCurve(hs->PolygonalBoundary, *profile.get(), conv)) { + IFCImporter::LogError("expected valid polyline for boundary of boolean halfspace"); + return; + } + + IfcMatrix4 mat; + ConvertAxisPlacement(mat,hs->Position); + profile->Transform(mat); + + // project the profile onto the plane (orthogonally along the plane normal) + IfcVector3 r; + bool have_r = false; + BOOST_FOREACH(IfcVector3& vec, profile->verts) { + vec = vec + ((p - vec) * n) * n; + ai_assert(fabs((vec-p) * n) < 1e-6); + + if (!have_r && (vec-p).SquareLength() > 1e-8) { + r = vec-p; + have_r = true; + } + } + + if (!have_r) { + IFCImporter::LogError("polyline for boundary of boolean halfspace is degenerate"); + return; + } + + // and map everything into a plane coordinate space so all intersection + // tests can be done in 2D space. + IfcMatrix4 proj; + + r.Normalize(); + + IfcVector3 u = n ^ r; + u.Normalize(); + + proj.a1 = r.x; + proj.a2 = r.y; + proj.a3 = r.z; + proj.b1 = u.x; + proj.b2 = u.y; + proj.b3 = u.z; + proj.c1 = n.x; + proj.c2 = n.y; + proj.c3 = n.z; + BOOST_FOREACH(IfcVector3& vec, profile->verts) { + vec *= proj; + } + + const IfcMatrix4 proj_inv = IfcMatrix4(proj).Inverse(); + + // clip the current contents of `meshout` against the plane we obtained from the second operand + const std::vector& in = first_operand.verts; + std::vector& outvert = result.verts; + + std::vector::const_iterator begin = first_operand.vertcnt.begin(), + end = first_operand.vertcnt.end(), iit; + + outvert.reserve(in.size()); + result.vertcnt.reserve(first_operand.vertcnt.size()); + + std::vector intersected_boundary_segments; + std::vector intersected_boundary_points; + + unsigned int vidx = 0; + for(iit = begin; iit != end; vidx += *iit++) { + if (!*iit) { + continue; + } + + unsigned int newcount = 0; + bool was_outside_boundary = !PointInPoly(in[vidx], profile->verts); + + size_t last_intersected_boundary_segment; + IfcVector3 last_intersected_boundary_point; + + bool extra_point_flag = false; + IfcVector3 extra_point; + + for(unsigned int i = 0; i < *iit; ++i) { + // current segment: [i,i+1 mod size] or [*extra_point,i] if extra_point_flag is set + const IfcVector3& e0 = extra_point_flag ? extra_point : in[vidx+i]; + const IfcVector3& e1 = extra_point_flag ? in[vidx+i] : in[vidx+(i+1)%*iit]; + + // does the current segment intersect the polygonal boundary? + const IfcVector3& e0_plane = proj * e0; + const IfcVector3& e1_plane = proj * e1; + + intersected_boundary_segments.clear(); + intersected_boundary_points.clear(); + + const bool is_boundary_intersection = IntersectsBoundaryProfile(e0_plane, e1_plane, profile->verts, + intersected_boundary_segments, + intersected_boundary_points); + + const bool is_outside_boundary = is_boundary_intersection ? !was_outside_boundary : was_outside_boundary; + + // does the current segment intersect the plane? + // (no extra check if this is an extra point) + IfcVector3 isectpos; + const Intersect isect = extra_point_flag ? Intersect_No : IntersectSegmentPlane(p,n,e0,e1,isectpos); + + // is it on the side of the plane that we keep? + const bool is_white_side =(e0-p).Normalize()*n > 0; + + // e0 on good side of plane? (i.e. we should keep geometry on this side) + if (is_white_side) { + // but is there an intersection in e0-e1 and is e1 in the clipping + // boundary? In this case, generate a line that only goes to the + // intersection point. + if (isect == Intersect_Yes && PointInPoly(e1, profile->verts)) { + outvert.push_back(e0); + ++newcount; + + outvert.push_back(isectpos); + ++newcount; + + // this is, however, only a line that goes to the plane, but not + // necessarily to the point where the bounding volume on the + // black side of the plane is hit. So basically, we need another + // check for [isectpos-e1], which should give an intersection + // point and also set the last_intersected_boundary_*'s. + extra_point_flag = true; + extra_point = isectpos; + continue; + } + else { + outvert.push_back(e0); + ++newcount; + } + } + // e0 on bad side of plane (i.e. we should remove geometry on this side, + // but only if it is within the bounding volume). + else if (isect == Intersect_Yes) { + if (is_boundary_intersection) {} + // drop it and keep e1 instead + outvert.push_back(isectpos); + ++newcount; + } + else { + + // did we just pass the boundary line? + if (is_boundary_intersection) { + + // and are now outside the clipping boundary? + if (is_outside_boundary) { + // in this case, get the point where the clipping boundary + // was entered first. Then, get the point where the clipping + // boundary volume was left! These two points with the plane + // normal form another plane that intersects the clipping + // volume. There are two ways to get from the first to the + // second point along the intersection curve, try to pick the + // one that lies within the current polygon. + + // TODO this approach doesn't handle all cases + + // ... + + outvert.push_back(proj_inv * intersected_boundary_points.back()); + ++newcount; + + outvert.push_back(e1); + ++newcount; + } + else { + // we just entered the clipping boundary. Record the point + // and the segment where we entered and also generate this point. + last_intersected_boundary_segment = intersected_boundary_segments.front(); + last_intersected_boundary_point = intersected_boundary_points.front(); + + outvert.push_back(e0); + ++newcount; + + outvert.push_back(proj_inv * last_intersected_boundary_point); + ++newcount; + } + } + // if not, we just keep the vertex + else if (is_outside_boundary) { + outvert.push_back(e0); + ++newcount; + } + } + + was_outside_boundary = is_outside_boundary; + extra_point_flag = false; + } + + if (!newcount) { + continue; + } + + IfcVector3 vmin,vmax; + ArrayBounds(&*(outvert.end()-newcount),newcount,vmin,vmax); + + // filter our IfcFloat points - those may happen if a point lies + // directly on the intersection line. However, due to IfcFloat + // precision a bitwise comparison is not feasible to detect + // this case. + const IfcFloat epsilon = (vmax-vmin).SquareLength() / 1e6f; + FuzzyVectorCompare fz(epsilon); + + std::vector::iterator e = std::unique( outvert.end()-newcount, outvert.end(), fz ); + + if (e != outvert.end()) { + newcount -= static_cast(std::distance(e,outvert.end())); + outvert.erase(e,outvert.end()); + } + if (fz(*( outvert.end()-newcount),outvert.back())) { + outvert.pop_back(); + --newcount; + } + if(newcount > 2) { + result.vertcnt.push_back(newcount); + } + else while(newcount-->0) { + result.verts.pop_back(); + } + + } + IFCImporter::LogDebug("generating CSG geometry by plane clipping with polygonal bounding (IfcBooleanClippingResult)"); } // ------------------------------------------------------------------------------------------------ @@ -283,6 +625,7 @@ void ProcessBoolean(const IfcBooleanResult& boolean, TempMesh& result, Conversio } if(hs) { + const IfcPolygonalBoundedHalfSpace* const hs_bounded = clip->SecondOperand->ResolveSelectPtr(conv.db); if (hs_bounded) { ProcessPolygonalBoundedBooleanHalfSpaceDifference(hs_bounded, result, first_operand, conv);