Merge branch 'feature/improve-point-sampler-tests' into 'master'

Feature/improve point sampler tests

See merge request tools/frackit!199

frackit/geometryutilities/getboundingbox.hh

@@ -24,6 +24,7 @@
 #ifndef FRACKIT_GEOMETRY_GET_BOUNDING_BOX_HH
 #define FRACKIT_GEOMETRY_GET_BOUNDING_BOX_HH
 
+#include <frackit/geometry/hollowcylinder.hh>
 #include <frackit/occ/breputilities.hh>
 
 namespace Frackit {
@@ -41,6 +42,14 @@ template<class Geo>
 auto getBoundingBox(const Geo& geo)
 { return OCCUtilities::getBoundingBox( OCCUtilities::getShape(geo) ); }
 
+/*!
+ * \ingroup GeometryUtilities
+ * \brief Overload for hollow cylinders.
+ */
+template<class ctype>
+auto getBoundingBox(const HollowCylinder<ctype>& hollowCyl)
+{ return getBoundingBox( hollowCyl.fullCylinder() ); }
+
 } // end namespace Frackit
 
 #endif // FRACKIT_GEOMETRY_GET_BOUNDING_BOX_HH

frackit/occ/breputilities.hh

@@ -27,6 +27,7 @@
 #ifndef FRACKIT_BREP_UTILITIES_HH
 #define FRACKIT_BREP_UTILITIES_HH
 
+#include <cmath>
 #include <vector>
 #include <algorithm>
 #include <stdexcept>
@@ -39,6 +40,9 @@
 #include <gp_Pnt.hxx>
 #include <gp_Dir.hxx>
 #include <gp_Elips.hxx>
+#include <gp_Trsf.hxx>
+#include <gp_Ax1.hxx>
+#include <gp_Ax2.hxx>
 
 // shape classes from TopoDS package
 #include <TopoDS.hxx>
@@ -77,6 +81,8 @@
 // BRep primitives and operations
 #include <BRepPrimAPI_MakeCylinder.hxx>
 #include <BRepPrimAPI_MakeSphere.hxx>
+#include <BRepPrimAPI_MakeRevol.hxx>
+#include <BRepPrimAPI_MakePrism.hxx>
 #include <BRepPrimAPI_MakeBox.hxx>
 
 // internal geometry classes
@@ -88,9 +94,11 @@
 #include <frackit/geometry/quadrilateral.hh>
 #include <frackit/geometry/polygon.hh>
 #include <frackit/geometry/cylinder.hh>
+#include <frackit/geometry/hollowcylinder.hh>
 #include <frackit/geometry/cylindersurface.hh>
 #include <frackit/geometry/box.hh>
 #include <frackit/geometry/boundingbox.hh>
+#include <frackit/geometry/boundingbox.hh>
 #include <frackit/geometry/sphere.hh>
 
 #include <frackit/geometryutilities/name.hh>
@@ -170,49 +178,6 @@ namespace OCCUtilities {
     TopoDS_Edge getShape(const Segment<ctype, worldDim>& segment)
     { return makeEdge(segment.source(), segment.target()); }
 
-    /*!
-     * \ingroup OpenCascade
-     * \brief Get the BRep of a cylinder.
-     */
-    template<class ctype>
-    TopoDS_Solid getShape(const Cylinder<ctype>& cylinder)
-    {
-        const auto& lateral = cylinder.lateralFace();
-        const auto& bottom = lateral.lowerBoundingCircle();
-        auto axis = direction(bottom.normal());
-        auto base1 = direction(bottom.base1());
-        auto center = point(bottom.center());
-        BRepPrimAPI_MakeCylinder makeCylinder(gp_Ax2(center, axis, base1),
-                                              bottom.radius(),
-                                              lateral.height(),
-                                              2.0*M_PI);
-        makeCylinder.Build();
-        if (!makeCylinder.IsDone())
-            throw std::runtime_error("Could not build cylinder");
-        return TopoDS::Solid(makeCylinder.Shape());
-    }
-
-    /*!
-     * \ingroup OpenCascade
-     * \brief Get the BRep of a cylinder surface.
-     */
-    template<class ctype>
-    TopoDS_Face getShape(const CylinderSurface<ctype>& cylSurface)
-    {
-        const auto& bottom = cylSurface.lowerBoundingCircle();
-        auto axis = direction(bottom.normal());
-        auto base1 = direction(bottom.base1());
-        auto center = point(bottom.center());
-        BRepPrimAPI_MakeCylinder makeCylinder(gp_Ax2(center, axis, base1),
-                                              bottom.radius(),
-                                              cylSurface.height(),
-                                              2.0*M_PI);
-        makeCylinder.Build();
-        if (!makeCylinder.IsDone())
-            throw std::runtime_error("Could not build cylinder");
-        return makeCylinder.Cylinder().LateralFace();
-    }
-
     /*!
      * \ingroup OpenCascade
      * \brief Get the BRep of an ellipse in 3d space.
@@ -320,9 +285,32 @@ namespace OCCUtilities {
         return BRepBuilderAPI_MakeFace(wire);
     }
 
+    /*!
+     * \ingroup OpenCascade
+     * \brief Get the BRep of a cylinder surface.
+     */
+    template<class ctype>
+    TopoDS_Face getShape(const CylinderSurface<ctype>& cylSurface)
+    {
+        const auto& bottom = cylSurface.lowerBoundingCircle();
+        auto axis = direction(bottom.normal());
+        auto base1 = direction(bottom.base1());
+        auto center = point(bottom.center());
+        BRepPrimAPI_MakeCylinder makeCylinder(gp_Ax2(center, axis, base1),
+                                              bottom.radius(),
+                                              cylSurface.height(),
+                                              2.0*M_PI);
+        makeCylinder.Build();
+        if (!makeCylinder.IsDone())
+            throw std::runtime_error("Could not build cylinder");
+        return makeCylinder.Cylinder().LateralFace();
+    }
+
     /*!
      * \ingroup OpenCascade
      * \brief Get the BRep of a box.
+     * \todo TODO: Once boxes are more generic (than axis-aligned
+     *             bounding boxes), this must be adapted
      */
     template<class ctype>
     TopoDS_Solid getShape(const Box<ctype>& box)
@@ -336,6 +324,75 @@ namespace OCCUtilities {
         return TopoDS::Solid(makeBox.Shape());
     }
 
+    /*!
+     * \ingroup OpenCascade
+     * \brief Get the BRep of a bounding box.
+     */
+    template<class ctype>
+    TopoDS_Solid getShape(const BoundingBox<ctype>& bbox)
+    {
+        const gp_Pnt pMin(bbox.xMin(), bbox.yMin(), bbox.zMin());
+        const gp_Pnt pMax(bbox.xMax(), bbox.yMax(), bbox.zMax());
+        BRepPrimAPI_MakeBox makeBox(pMin, pMax);
+        makeBox.Build();
+        if (!makeBox.IsDone())
+            throw std::runtime_error("Could not build box");
+        return TopoDS::Solid(makeBox.Shape());
+    }
+
+    /*!
+     * \ingroup OpenCascade
+     * \brief Get the BRep of a cylinder.
+     */
+    template<class ctype>
+    TopoDS_Solid getShape(const Cylinder<ctype>& cylinder)
+    {
+        const auto& lateral = cylinder.lateralFace();
+        const auto& bottom = lateral.lowerBoundingCircle();
+        auto axis = direction(bottom.normal());
+        auto base1 = direction(bottom.base1());
+        auto center = point(bottom.center());
+        BRepPrimAPI_MakeCylinder makeCylinder(gp_Ax2(center, axis, base1),
+                                              bottom.radius(),
+                                              lateral.height(),
+                                              2.0*M_PI);
+        makeCylinder.Build();
+        if (!makeCylinder.IsDone())
+            throw std::runtime_error("Could not build cylinder");
+        return TopoDS::Solid(makeCylinder.Shape());
+    }
+
+    /*!
+     * \ingroup OpenCascade
+     * \brief Get the BRep of a hollow cylinder.
+     */
+    template<class ctype>
+    TopoDS_Solid getShape(const HollowCylinder<ctype>& hollowCyl)
+    {
+        const auto bottomCenter = hollowCyl.bottomInnerCircle().center();
+        auto p1 = bottomCenter; p1 += hollowCyl.base1()*hollowCyl.innerRadius();
+        auto p2 = p1; p2 += hollowCyl.base1()*(hollowCyl.outerRadius() - hollowCyl.innerRadius());
+
+        // make side segment and rotate around axis
+        const auto& segShape = getShape(Segment<ctype, 3>{p1, p2});
+        const gp_Ax1 cylAxis(point(bottomCenter), direction(hollowCyl.axis()));
+        BRepPrimAPI_MakeRevol makeBottom(segShape, cylAxis, 2.0*M_PI);
+        makeBottom.Build();
+        if (!makeBottom.IsDone())
+            throw std::runtime_error("Could not build hollow cylinder bottom");
+        const auto& bottomShape = makeBottom.Shape();
+
+        // extrude bottom
+        const auto& cylAxisSeg = hollowCyl.centerSegment();
+        const auto& v = vector(Vector<ctype, 3>{cylAxisSeg.source(), cylAxisSeg.target()});
+        BRepPrimAPI_MakePrism makeExtrusion(bottomShape, v);
+        makeExtrusion.Build();
+        if (!makeExtrusion.IsDone())
+            throw std::runtime_error("Could not extrude hollow cylinder bottom");
+
+        return TopoDS::Solid(makeExtrusion.Shape());
+    }
+
     /*!
      * \ingroup OpenCascade
      * \brief Get the BRep of a sphere.

frackit/python/geometry/boundingbox.hh

@@ -44,6 +44,7 @@ void registerBoundingBox(py::module& module)
             "xmin"_a, "ymin"_a, "zmin"_a, "xmax"_a, "ymax"_a, "zmax"_a);
     cls.def(py::init<const Point&, const Point&>(), "firstCorner"_a, "secondCorner"_a);
     cls.def(py::init<const Box&>(), "box"_a);
+    cls.def(py::init<>());
 
     // dimensionality properties
     registerDimensionProperties(cls);

frackit/python/geometry/getboundingbox.hh

@@ -27,6 +27,7 @@
 #include <frackit/geometry/boundingbox.hh>
 #include <frackit/geometry/circle.hh>
 #include <frackit/geometry/cylinder.hh>
+#include <frackit/geometry/hollowcylinder.hh>
 #include <frackit/geometry/cylindersurface.hh>
 #include <frackit/geometry/disk.hh>
 #include <frackit/geometry/ellipse.hh>
@@ -35,6 +36,7 @@
 #include <frackit/geometry/quadrilateral.hh>
 #include <frackit/geometry/segment.hh>
 #include <frackit/geometry/sphere.hh>
+#include <frackit/geometry/point.hh>
 
 #include <frackit/python/geometry/ctype.hh>
 #include <frackit/python/geometry/brepwrapper.hh>
@@ -66,6 +68,7 @@ void registerGetBoundingBox(py::module& module)
 {
     Detail::registerGetBoundingBox<Box<ctype>>(module, "box.");
     Detail::registerGetBoundingBox<Cylinder<ctype>>(module, "cylinder.");
+    Detail::registerGetBoundingBox<HollowCylinder<ctype>>(module, "hollow cylinder.");
     Detail::registerGetBoundingBox<CylinderSurface<ctype>>(module, "cylinder surface.");
     Detail::registerGetBoundingBox<Disk<ctype>>(module, "disk.");
     Detail::registerGetBoundingBox<Sphere<ctype>>(module, "sphere.");
@@ -76,6 +79,7 @@ void registerGetBoundingBox(py::module& module)
     Detail::registerGetBoundingBox<Polygon<ctype, 3>>(module, "polygon in 3d space.");
     Detail::registerGetBoundingBox<Quadrilateral<ctype, 3>>(module, "quadrilateral in 3d space.");
     Detail::registerGetBoundingBox<Segment<ctype, 3>>(module, "segment in 3d space.");
+    Detail::registerGetBoundingBox<Point<ctype, 3>>(module, "point in 3d space.");
 
     Detail::registerGetBoundingBox<ShapeWrapper>(module, "(wrapped) OpenCascade shape.");
     Detail::registerGetBoundingBox<VertexWrapper>(module, "(wrapped) OpenCascade vertex.");

frackit/python/occutilities/breputilities.hh

@@ -39,6 +39,7 @@
 #include <frackit/geometry/disk.hh>
 #include <frackit/geometry/box.hh>
 #include <frackit/geometry/cylinder.hh>
+#include <frackit/geometry/hollowcylinder.hh>
 #include <frackit/geometry/cylindersurface.hh>
 #include <frackit/geometry/sphere.hh>
 #include <frackit/python/geometry/brepwrapper.hh>
@@ -333,6 +334,7 @@ void registerBRepUtilities(pybind11::module& module)
     module.def("getShape", &OCCUtilities::getShape<Disk<ctype>>, "Returns the OCC BRep of a disk");
     module.def("getShape", &OCCUtilities::getShape<Box<ctype>>, "Returns the OCC BRep of a box");
     module.def("getShape", &OCCUtilities::getShape<Cylinder<ctype>>, "Returns the OCC BRep of a cylinder");
+    module.def("getShape", &OCCUtilities::getShape<HollowCylinder<ctype>>, "Returns the OCC BRep of a hollow cylinder");
     module.def("getShape", &OCCUtilities::getShape<CylinderSurface<ctype>>, "Returns the OCC BRep of the lateral surface of a cylinder");
     module.def("getShape", &OCCUtilities::getShape<Sphere<ctype>>, "Returns the OCC BRep of a sphere");
 

test/sampling/test_pointsampling.cc

 #include <stdexcept>
 #include <string>
-#include <array>
+#include <vector>
+#include <cmath>
 
 #include <BRepTools.hxx>
 #include <BRep_Builder.hxx>
@@ -8,51 +9,89 @@
 
 #include <frackit/sampling/makeuniformpointsampler.hh>
 #include <frackit/geometry/point.hh>
+#include <frackit/geometry/box.hh>
+#include <frackit/geometry/boundingbox.hh>
 #include <frackit/geometry/cylinder.hh>
 #include <frackit/geometry/hollowcylinder.hh>
+#include <frackit/geometryutilities/getboundingbox.hh>
 #include <frackit/occ/breputilities.hh>
 
-//! test random sampling of points on geometries
-int main()
+//! sample points on a geometry
+//! if a brep file name is provided, the result is written
+//! into a brep file, such that the result can be checked visually
+template<class G>
+void testPointSampling(const G& geometry,
+                       std::size_t numSamples,
+                       const std::string& brepFileName = "")
 {
-    using ctype = double;
-
     using namespace Frackit;
-    using Cylinder = Cylinder<ctype>;
+    using ctype = typename G::ctype;
+    static constexpr int worldDim = G::worldDimension();
 
-    Cylinder cylinder(0.5, 1.0);
-    auto cylPointSampler = makeUniformPointSampler(cylinder);
+    // sample points
+    auto pointSampler = makeUniformPointSampler(geometry);
+    std::vector<Frackit::Point<ctype, worldDim>> points(numSamples);
+    for (unsigned int i = 0; i < numSamples; ++i) points[i] = pointSampler();
 
-    // sample 500 points
-    std::array<Frackit::Point<ctype, 3>, 500> points;
-    for (unsigned int i = 0; i < 500; ++i)
-        points[i] = cylPointSampler();
+    // check if all of the points are contained in the geometry
+    if ( std::any_of(points.begin(), points.end(),
+                     [&geometry] (const auto& p)
+                     { return !geometry.contains(p); }) )
+        throw std::runtime_error("Point not contained in geometry");
 
-    // create a single compound and write to .brep file
-    // build a single compound shape
-    BRep_Builder b;
-    TopoDS_Compound c;
-    b.MakeCompound(c);
+    // compute bounding box of the point cloud
+    BoundingBox<ctype> bbpc;
     for (const auto& p : points)
-        b.Add(c, OCCUtilities::getShape(p));
-    b.Add(c, OCCUtilities::getShape(cylinder));
+        bbpc += getBoundingBox(p);
 
-    BRepTools::Write(c, "cylinderpoints.brep");
+    // we expect less than 5% deviation to the geometry's bbox
+    const auto bbg = getBoundingBox(geometry);
+    const auto dx = bbg.xMax() - bbg.xMin();
+    const auto dy = bbg.yMax() - bbg.yMin();
+    const auto dz = bbg.zMax() - bbg.zMin();
 
-    // test sampling on a hollow cylinder
-    HollowCylinder<ctype> hollowCylinder(0.4, 0.5, 1.0);
-    cylPointSampler = makeUniformPointSampler(hollowCylinder);
-    for (unsigned int i = 0; i < 500; ++i)
-        points[i] = cylPointSampler();
+    using std::abs;
+    if ( abs(bbpc.xMin() - bbg.xMin()) > 0.05*dx ) throw std::runtime_error("Bounding box deviation > 5%");
+    if ( abs(bbpc.xMax() - bbg.xMax()) > 0.05*dx ) throw std::runtime_error("Bounding box deviation > 5%");
+    if ( abs(bbpc.yMin() - bbg.yMin()) > 0.05*dy ) throw std::runtime_error("Bounding box deviation > 5%");
+    if ( abs(bbpc.yMax() - bbg.yMax()) > 0.05*dy ) throw std::runtime_error("Bounding box deviation > 5%");
+    if ( abs(bbpc.zMin() - bbg.zMin()) > 0.05*dz ) throw std::runtime_error("Bounding box deviation > 5%");
+    if ( abs(bbpc.zMax() - bbg.zMax()) > 0.05*dz ) throw std::runtime_error("Bounding box deviation > 5%");
 
-    // create a single compound and write to .brep file
-    // build a single compound shape
-    TopoDS_Compound c2;
-    b.MakeCompound(c2);
-    for (const auto& p : points)
-        b.Add(c2, OCCUtilities::getShape(p));
-    b.Add(c2, OCCUtilities::getShape(cylinder));
-    BRepTools::Write(c2, "cylinderpoints_hollow.brep");
+    std::cout << "Successfully tested geometry: " << geometry.name() << std::endl;
+    std::cout << "\t- Deviation in x-direction: " << abs(dx - bbpc.xMax() + bbpc.xMin())/dx*100.0 << "%" << std::endl;
+    std::cout << "\t- Deviation in y-direction: " << abs(dy - bbpc.yMax() + bbpc.yMin())/dy*100.0 << "%" << std::endl;
+    std::cout << "\t- Deviation in z-direction: " << abs(dz - bbpc.zMax() + bbpc.zMin())/dz*100.0 << "%" << std::endl;
+
+    if (brepFileName != "")
+    {
+        // create a single compound and write to .brep file
+        BRep_Builder b;
+        TopoDS_Compound c;
+        b.MakeCompound(c);
+        for (const auto& p : points)
+            b.Add(c, OCCUtilities::getShape(p));
+        b.Add(c, OCCUtilities::getShape(geometry));
+
+        BRepTools::Write(c, brepFileName.c_str());
+    }
+}
+
+//! test random sampling of points on geometries
+int main()
+{
+    using ctype = double;
+
+    using namespace Frackit;
+    using Box = Box<ctype>;
+    using BBox = BoundingBox<ctype>;
+    using Cylinder = Cylinder<ctype>;
+    using HollowCyl = HollowCylinder<ctype>;
+
+    testPointSampling(Box(0.,0.,0.,1.,1.,1.), 1000, "pointsamples_box.brep");
+    testPointSampling(BBox(0.,0.,0.,1.,1.,1.), 1000, "pointsamples_bbox.brep");
+    testPointSampling(Cylinder(0.5, 1.0), 1000, "pointsamples_cyl.brep");
+    testPointSampling(HollowCyl(0.4, 0.5, 1.0), 1000, "pointsamples_hollowcyl.brep");
 
     std::cout << "All tests passed" << std::endl;
     return 0;