Skip to content

Commit

Permalink
Optimize rebars (sweeps)
Browse files Browse the repository at this point in the history
  • Loading branch information
QuimMoya committed Aug 28, 2023
1 parent 1f1eace commit 5e5443a
Show file tree
Hide file tree
Showing 3 changed files with 277 additions and 19 deletions.
14 changes: 11 additions & 3 deletions src/wasm/geometry/IfcGeometryProcessor.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -666,6 +666,9 @@ namespace webifc::geometry
}
case schema::IFCSWEPTDISKSOLID:
{
_expressIDToGeometry[1] = predefinedCylinder;
_expressIDToGeometry[2] = predefinedCube;

// TODO: prevent self intersections in Sweep function still not working properly

bool closed = false;
Expand Down Expand Up @@ -703,7 +706,7 @@ namespace webifc::geometry
IfcProfile profile;
profile.curve = GetCircleCurve(radius, _circleSegments);

IfcGeometry geom = Sweep(_geometryLoader.GetLinearScalingFactor(), closed, profile, directrix);
IfcGeometry geom = SweepCircular(_geometryLoader.GetLinearScalingFactor(), mesh, _optimize_profiles, closed, profile, radius, directrix);

_expressIDToGeometry[expressID] = geom;
mesh.expressID = expressID;
Expand Down Expand Up @@ -760,6 +763,8 @@ namespace webifc::geometry
}
case schema::IFCEXTRUDEDAREASOLID:
{
_expressIDToGeometry[1] = predefinedCylinder;
_expressIDToGeometry[2] = predefinedCube;
_loader.MoveToArgumentOffset(expressID, 0);
uint32_t profileID = _loader.GetRefArgument();
uint32_t placementID = _loader.GetOptionalRefArgument();
Expand All @@ -776,6 +781,11 @@ namespace webifc::geometry
_loader.MoveToArgumentOffset(profileID, 2);
uint32_t profilePlacementID = _loader.GetRefArgument();
double radius = _loader.GetDoubleArgument();

// std::cout << radius << std::endl;
// std::cout << depth << std::endl;
// std::cout << profilePlacementID << std::endl;

// double thickness = _loader.GetDoubleArgument(); // Read this property only in hollow profiles

if (placementID)
Expand Down Expand Up @@ -819,7 +829,6 @@ namespace webifc::geometry
extrusionScale *= profileTransform;
mesh.transformation *= extrusionScale;

_expressIDToGeometry[1] = predefinedCylinder;
mesh.expressID = 1;
mesh.hasGeometry = true;
return mesh;
Expand Down Expand Up @@ -877,7 +886,6 @@ namespace webifc::geometry
extrusionScale *= profileTransform;
mesh.transformation *= extrusionScale;

_expressIDToGeometry[2] = predefinedCube;
mesh.expressID = 2;
mesh.hasGeometry = true;
return mesh;
Expand Down
259 changes: 255 additions & 4 deletions src/wasm/geometry/operations/geometryutils.h
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,7 @@ namespace webifc::geometry {

// compute curve for each part of the directrix
std::vector<IfcCurve> curves;
std::vector<glm::dmat4> transforms;

for (size_t i = 0; i < dpts.size(); i++)
{
Expand Down Expand Up @@ -205,7 +206,7 @@ namespace webifc::geometry {
pt = -pt2D.x * right - pt2D.y * left + planeOrigin;
}
glm::dvec3 proj = projectOntoPlane(planeOrigin, planeNormal, pt, directrixSegmentNormal);

segmentForCurve.Add(proj);
}
}
Expand Down Expand Up @@ -245,8 +246,6 @@ namespace webifc::geometry {

//Only segments smaller than 10 cm will be represented, those that are bigger will be standardized

// if(di < 0.1 / scaling)
// {
const auto &c1 = curves[i - 1].points;
const auto &c2 = curves[i].points;

Expand All @@ -262,7 +261,259 @@ namespace webifc::geometry {
geom.AddFace(tl, br, bl);
geom.AddFace(tl, tr, br);
}
// }

}

// DumpSVGCurve(directrix.points, glm::dvec3(), "directrix.html");
// DumpIfcGeometry(geom, "sweep.obj");

return geom;
}

inline IfcGeometry SweepCircular(const double scaling, IfcComposedMesh &mesh, const bool optimizeProfiles, const bool closed, const IfcProfile &profile, const double radius, const IfcCurve &directrix, const glm::dvec3 &initialDirectrixNormal = glm::dvec3(0), const bool rotate90 = false)
{
IfcGeometry geom;

std::vector<glm::vec<3, glm::f64>> dpts;

// Remove repeated points
for (size_t i = 0; i < directrix.points.size(); i++)
{
if (i < directrix.points.size() - 1)
{
if (glm::distance(directrix.points[i], directrix.points[i + 1]) > EPS_BIG2 / scaling)
{
dpts.push_back(directrix.points[i]);
}
}
else
{
dpts.push_back(directrix.points[i]);
}
}

if (closed)
{
glm::vec<3, glm::f64> dirStart = dpts[dpts.size() - 2] - dpts[dpts.size() - 1];
glm::vec<3, glm::f64> dirEnd = dpts[1] - dpts[0];
std::vector<glm::vec<3, glm::f64>> newDpts;
newDpts.push_back(dpts[0] + dirStart);
for (size_t i = 0; i < dpts.size(); i++)
{
newDpts.push_back(dpts[i]);
}
newDpts.push_back(dpts[dpts.size() - 1] + dirEnd);
dpts = newDpts;
}

if (dpts.size() <= 1)
{
// nothing to sweep
return geom;
}

// compute curve for each part of the directrix
std::vector<IfcCurve> curves;
std::vector<glm::dmat4> transforms;

for (size_t i = 0; i < dpts.size(); i++)
{
IfcCurve segmentForCurve;

glm::dvec3 directrix2;
glm::dvec3 planeNormal;
glm::dvec3 directrixSegmentNormal;
glm::dvec3 planeOrigin;

if (i == 0) // start
{
planeNormal = glm::normalize(dpts[1] - dpts[0]);
directrixSegmentNormal = planeNormal;
planeOrigin = dpts[0];
directrix2 = planeNormal;
}
else if (i == dpts.size() - 1) // end
{
planeNormal = glm::normalize(dpts[i] - dpts[i - 1]);
directrixSegmentNormal = planeNormal;
planeOrigin = dpts[i];
directrix2 = planeNormal;
}
else // middle
{
// possibly the directrix is bad
glm::dvec3 n1 = glm::normalize(dpts[i] - dpts[i - 1]);
glm::dvec3 n2 = glm::normalize(dpts[i + 1] - dpts[i]);
glm::dvec3 p = glm::normalize(glm::cross(n1, n2));
directrix2 = -n1;

// double prod = glm::dot(n1, n2);

if (std::isnan(p.x))
{
// TODO: sometimes outliers cause the perp to become NaN!
// this is bad news, as it nans the points added to the final mesh
// also, it's hard to bail out now :/
// see curve.add() for more info on how this is currently "solved"
printf("NaN perp!\n");
}

glm::dvec3 u1 = glm::normalize(glm::cross(n1, p));
glm::dvec3 u2 = glm::normalize(glm::cross(n2, p));

// TODO: When n1 and n2 have similar direction but opposite side...
// ... projection tend to infinity. -> glm::dot(n1, n2)
// I implemented a bad solution to prevent projection to infinity
if (glm::dot(n1, n2) < -0.9)
{
n2 = -n2;
u2 = -u2;
}

glm::dvec3 au = glm::normalize(u1 + u2);
planeNormal = glm::normalize(glm::cross(au, p));
directrixSegmentNormal = n1; // n1 or n2 doesn't matter

planeOrigin = dpts[i];
}

glm::dvec3 dz = glm::normalize(directrix2);
glm::dvec3 dx = glm::dvec3(1, 0, 0);
glm::dvec3 dy = glm::dvec3(0, 1, 0);

double parallelZ = glm::abs(glm::dot(dz, glm::dvec3(0, 0, 1)));

if(parallelZ > 1 - EPS_BIG2)
{
dx = glm::normalize(glm::cross(dz, glm::dvec3(0, 1, 0)));
} else {
dx = glm::normalize(glm::cross(dz, glm::dvec3(0, 0, 1)));
}

dy = glm::normalize(glm::cross(dz, dx));

glm::dmat4 profileScale = glm::dmat4(
glm::dvec4(dx * radius, 0),
glm::dvec4(dy * radius, 0),
glm::dvec4(dz, 0),
glm::dvec4(planeOrigin, 1));

transforms.push_back(profileScale);

if (curves.empty())
{
// construct initial curve
glm::dvec3 left;
glm::dvec3 right;
if (initialDirectrixNormal == glm::dvec3(0))
{
left = glm::cross(directrixSegmentNormal, glm::dvec3(directrixSegmentNormal.y, directrixSegmentNormal.x, directrixSegmentNormal.z));
if (left == glm::dvec3(0, 0, 0))
{
left = glm::cross(directrixSegmentNormal, glm::dvec3(directrixSegmentNormal.x, directrixSegmentNormal.z, directrixSegmentNormal.y));
}
if (left == glm::dvec3(0, 0, 0))
{
left = glm::cross(directrixSegmentNormal, glm::dvec3(directrixSegmentNormal.z, directrixSegmentNormal.y, directrixSegmentNormal.x));
}
right = glm::normalize(glm::cross(directrixSegmentNormal, left));
left = glm::normalize(glm::cross(directrixSegmentNormal, right));
}
else
{
left = glm::cross(directrixSegmentNormal, initialDirectrixNormal);
glm::dvec3 side = glm::normalize(initialDirectrixNormal);
right = glm::normalize(glm::cross(directrixSegmentNormal, left));
left = glm::normalize(glm::cross(directrixSegmentNormal, right));
right *= side;
}

if (left == glm::dvec3(0, 0, 0))
{
printf("0 left vec in sweep!\n");
}

// project profile onto planeNormal, place on planeOrigin
// TODO: look at holes
auto &ppts = profile.curve.points;
for (auto &pt2D : ppts)
{
glm::dvec3 pt = -pt2D.x * left + -pt2D.y * right + planeOrigin;
if(rotate90)
{
pt = -pt2D.x * right - pt2D.y * left + planeOrigin;
}
glm::dvec3 proj = projectOntoPlane(planeOrigin, planeNormal, pt, directrixSegmentNormal);

segmentForCurve.Add(proj);
}
}
else
{
// project previous curve onto the normal
const IfcCurve &prevCurve = curves.back();

auto &ppts = prevCurve.points;
for (auto &pt : ppts)
{
glm::dvec3 proj = projectOntoPlane(planeOrigin, planeNormal, pt, directrixSegmentNormal);

segmentForCurve.Add(proj);
}
}

if (!closed || (i != 0 && i != dpts.size() - 1))
{
curves.push_back(segmentForCurve);
}
}

if (closed)
{
dpts.pop_back();
dpts.erase(dpts.begin());
}

// connect the curves
for (size_t i = 1; i < dpts.size(); i++)
{
glm::dvec3 p1 = dpts[i - 1];
glm::dvec3 p2 = dpts[i];
glm::dvec3 dir = p1 - p2;
glm::dvec4 ddir = glm::dvec4(dir, 0);
const double di = glm::distance(p1, p2);

//Only segments smaller than 10 cm will be represented, those that are bigger will be standardized

if(!optimizeProfiles || di < 0.5 / scaling)
{
const auto &c1 = curves[i - 1].points;
const auto &c2 = curves[i].points;

uint32_t capSize = c1.size();
for (size_t j = 1; j < capSize; j++)
{
glm::dvec3 bl = c1[j - 1];
glm::dvec3 br = c1[j - 0];

glm::dvec3 tl = c2[j - 1];
glm::dvec3 tr = c2[j - 0];

geom.AddFace(tl, br, bl);
geom.AddFace(tl, tr, br);
}
}
else
{
transforms[i] = glm::dmat4(transforms[i][0], transforms[i][1], ddir, transforms[i][3]);
IfcComposedMesh newMesh;
newMesh.expressID = 1;
newMesh.hasColor = mesh.hasColor;
newMesh.color = mesh.color;
newMesh.hasGeometry = true;
newMesh.transformation = transforms[i];
mesh.children.push_back(newMesh);
}
}

// DumpSVGCurve(directrix.points, glm::dvec3(), "directrix.html");
Expand Down
Loading

0 comments on commit 5e5443a

Please sign in to comment.