How to increase speed/performance of boolean operations ?

[AlbertDeTerre]

EXPLAINATION

Hello,

I'm currently developping an algorithm that transforms a mesh into a primitive volume.

To do so, I start by dividing a bounding box around the mesh into 8 equal parts.
For each part, if the volume of the mesh contained within it is:

• Below a certain minimum threshold --> it is rejected.
• Over a specified maximum threshold --> we keep it.
• Between min and max threshold --> repeat the division process on each child, continuing until all remaining parts are retrieved

Finally, we assemble these parts to reconstruct the primitive shape, only composed of hexahedrons.

MY QUESTION

As this needs to do a lot of boolean operations, such as intersections to check volume and then compare to thresholds, I want to optimize the speed of these boolean operations. Is there a way to do that ?

OTHER INFORMATIONS

• Using exact_predicates_exact_constructions_kernel
• Using Surface_mesh
• Using CGAL 5.6.1

MY CODE

main function:

#include "iostream"
#include "vector"
#include "chrono"

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Lazy_exact_nt.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>

typedef CGAL::Exact_predicates_exact_constructions_kernel K;
typedef K::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Surface_mesh;

using namespace std;
namespace PMP = CGAL::Polygon_mesh_processing;

void calculate_min_max_coord(Surface_mesh, Point&, Point&);
Surface_mesh makeHexahedron(Point, Point);
vector<Surface_mesh> divide_equally(Surface_mesh);
vector<Surface_mesh> process_reshape(Surface_mesh, Surface_mesh &, double, double, int, int);

int main(int argc, char *argv[])
{
    // Declaring base variables
    Surface_mesh mesh_to_reshape, mesh_to_reshape_bbox, mergedMeshes;
    Point mesh_max, mesh_min;

    // Read file into mesh_to_reshape
    CGAL::IO::read_STL("./PipedPBase.stl", mesh_to_reshape);

    // Calcuate min and max coord to recreate the hexahedron bounding box of the mesh_to_reshape
    calculate_min_max_coord(mesh_to_reshape, mesh_min, mesh_max);

    // Create the bounding box
    mesh_to_reshape_bbox = makeHexahedron(mesh_min, mesh_max);

    // Calll the recursive function that will
    vector<Surface_mesh> list = process_reshape(mesh_to_reshape_bbox, mesh_to_reshape, 0.2, 0.8, 5, 0);

    // Build the primitive volume using all cubes
    for (auto mesh : list){
        PMP::corefine_and_compute_union(mergedMeshes, mesh, mergedMeshes);
    }

    // Write the result into a stl file
    CGAL::IO::write_STL("./output/merged_meshes.stl", mergedMeshes);

    return 0;
}

process_reshape function (Recursive function to keep or reject child):

vector<Surface_mesh> process_reshape(Surface_mesh mesh_to_divide, Surface_mesh original_mesh, double min_threshold, double max_threshold, int depth, int current_depth = 0)
{

    // Check if depth has been reached --> return the current mesh if it has
    if (current_depth == depth)
    {
        return vector<Surface_mesh>{mesh_to_divide};
    }

    // Compute the intersection_mesh and calculate the volume ratio
    Surface_mesh intersection_mesh;

    // Check if the meshes intersect
    bool has_intersection = CGAL::Polygon_mesh_processing::corefine_and_compute_intersection(original_mesh, mesh_to_divide, intersection_mesh, CGAL::parameters::throw_on_self_intersection(true));

    // If it has intersection -> compute the ratio bewteen volume of intersection mesh and the volume of the current cube
    auto volume_ratio = has_intersection ? (CGAL::to_double(CGAL::Polygon_mesh_processing::volume(intersection_mesh)) / CGAL::to_double(CGAL::Polygon_mesh_processing::volume(mesh_to_divide))) : 0.0;

    // If the volume ratio is between min and max threshold
    if (volume_ratio > min_threshold)
    {

        // If it's greater than max threshold --> return the current mesh
        if (volume_ratio > max_threshold)
        {
            return vector<Surface_mesh>{mesh_to_divide};
        }

        // If not --> divide the mesh into 8 equal children
        vector<Surface_mesh> created_children = divide_equally(mesh_to_divide);

        vector<Surface_mesh> to_return;

        // Iterate through children to call the recursion for each child
        for (Surface_mesh child : created_children)
        {
            vector<Surface_mesh> returned = process_reshape(child, intersection_mesh, min_threshold, max_threshold, depth, current_depth + 1);
            to_return.insert(to_return.end(), returned.begin(), returned.end());
        }

        return to_return;
    }
    else
    {
        // If the volume ratio is less than the min threshold, it means that we don't want to keep the current mesh
        return vector<Surface_mesh>{};
    }
}

divide_equally function (divide an hexahedron into 8 children):

vector<Surface_mesh> divide_equally(Surface_mesh mesh_to_divide)
{
    vector<Surface_mesh> result;

    // Calcuate min and max coordinates
    Point min_coord, max_coord;
    calculate_min_max_coord(mesh_to_divide, min_coord, max_coord);

    // Calculate length of each axis
    auto len_x = abs(max_coord.x() - min_coord.x()) / 2.0;
    auto len_y = abs(max_coord.y() - min_coord.y()) / 2.0;
    auto len_z = abs(max_coord.z() - min_coord.z()) / 2.0;

    auto x = min_coord.x(), y = min_coord.y(), z = min_coord.z();

    // Iterate 8 times to create 8 equal hexahedron
    for (int i = 0; i < 8; i++)
    {
        // Change x coordinate
        if ((i + 1) % 2 == 0)
            x = x == min_coord.x() ? min_coord.x() + len_x : min_coord.x();

        // Change y coordinate
        if (i != 0 && i % 2 == 0)
            y = y == min_coord.y() ? min_coord.y() + len_y : min_coord.y();

        // Change z coordinate
        if (i != 0 && i % 4 == 0 && i != 0)
            z = min_coord.z() + len_z;

        result.push_back(makeHexahedron(Point(x, y, z), Point(x + len_x, y + len_y, z + len_z)));
    }

    return result;
}

makeHexahedron function (create hexahedron using min and max coordinates):

vector<Surface_mesh> divide_equally(Surface_mesh mesh_to_divide)
{
    vector<Surface_mesh> result;

    // Calcuate min and max coordinates
    Point min_coord, max_coord;
    calculate_min_max_coord(mesh_to_divide, min_coord, max_coord);

    // Calculate length of each axis
    auto len_x = abs(max_coord.x() - min_coord.x()) / 2.0;
    auto len_y = abs(max_coord.y() - min_coord.y()) / 2.0;
    auto len_z = abs(max_coord.z() - min_coord.z()) / 2.0;

    auto x = min_coord.x(), y = min_coord.y(), z = min_coord.z();

    // Iterate 8 times to create 8 equal hexahedron
    for (int i = 0; i < 8; i++)
    {
        // Change x coordinate
        if ((i + 1) % 2 == 0)
            x = x == min_coord.x() ? min_coord.x() + len_x : min_coord.x();

        // Change y coordinate
        if (i != 0 && i % 2 == 0)
            y = y == min_coord.y() ? min_coord.y() + len_y : min_coord.y();

        // Change z coordinate
        if (i != 0 && i % 4 == 0 && i != 0)
            z = min_coord.z() + len_z;

        result.push_back(makeHexahedron(Point(x, y, z), Point(x + len_x, y + len_y, z + len_z)));
    }

    return result;
}

calculate_min_max_coord function (Compute min and max coordinates of a mesh):

void calculate_min_max_coord(Surface_mesh mesh, Point & min_coord, Point & max_coord)
{
    // Initialise max and min points
    Point pt_max = Point(-9999999999.0, -9999999999.0, -9999999999.0);
    Point pt_min = Point(9999999999.0, 9999999999.0, 9999999999.0);

    // Retrieve max and min from the mesh
    for (auto v : mesh.vertices())
    {
        Point p = mesh.point(v);
        pt_min = Point(CGAL::to_double(p.x() < pt_min.x() ? p.x() : pt_min.x()), CGAL::to_double(p.y() < pt_min.y() ? p.y() : pt_min.y()), CGAL::to_double(p.z() < pt_min.z() ? p.z() : pt_min.z()));
        pt_max = Point(CGAL::to_double(p.x() > pt_max.x() ? p.x() : pt_max.x()), CGAL::to_double(p.y() > pt_max.y() ? p.y() : pt_max.y()), CGAL::to_double(p.z() > pt_max.z() ? p.z() : pt_max.z()));
    }

    // Affect min and max to min_coord and max_coord
    min_coord = pt_min;
    max_coord = pt_max;
}

[sloriot]

It seems that one step could be replaced by what is done here, with an updated version of this class (currently in master branch, to be released in CGAL 6.0). You'll then be restricted to a smaller set of faces which can already improve the speed of the Boolean operation. I don't have to much time now but I'll have a closer look at your code later.

One thing that might be worth investigating is to use Monte Carlo to estimate the volume. Basically you can generate a sampling inside the volume. Then you can estimate the volume by counting the number of points inside. To do so, sample the cube using this class and using Side_of_triangle_mesh to select point inside the mesh. The contribution of each point to the volume being the volume of the cube/nb of sample points. That actually could be a nice example to be added to the orthtree package.

[AlbertDeTerre]

Thanks for your reply !

I found your first solution very interesting, but as it is not in a stable/official release of CGAL I don't know if it is really good for me to use it...

I'll have a look to Monte Carlo for the volume estimation but I don't think that it will be a game changer solution...

[sloriot]

I implemented the Monte Carlo to see how well it performs and it's pretty fast.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>

#include <CGAL/Octree.h>
#include <CGAL/Orthtree_traits_face_graph.h>
#include <CGAL/Side_of_triangle_mesh.h>
#include <CGAL/Random.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Mesh = CGAL::Surface_mesh<K::Point_3>;

using OTraits = CGAL::Orthtree_traits_face_graph<Mesh, Mesh::Property_map<Mesh::Vertex_index, K::Point_3>>;

namespace PMP = CGAL::Polygon_mesh_processing;

class Split_predicate_min_volume {

  double m_min_vol;
  std::shared_ptr<CGAL::Side_of_triangle_mesh<Mesh, K>> m_side_of_ptr;

  double estimate_volume(const std::vector<K::Point_3>& samples, double cube_vol) const
  {
    std::size_t inside=0;
    for (const K::Point_3& p : samples)
      if ((*m_side_of_ptr)(p)==CGAL::ON_BOUNDED_SIDE)
        ++inside;
    return cube_vol/samples.size() * inside;
  }

  void fill_bbox(const OTraits::Bbox_d& bb, std::vector<K::Point_3>& samples, std::size_t nb_samples) const
  {
    samples.reserve(nb_samples);

    CGAL::Random& rg = CGAL::get_default_random();
    double grid_dx = bb.xmax() - bb.xmin();
    double grid_dy = bb.ymax() - bb.ymin();
    double grid_dz = bb.zmax() - bb.zmin();

    for(std::size_t i=0; i < nb_samples; ++i){
      samples.push_back(
        K::Point_3(bb.xmin ()+ rg.get_double()* grid_dx,
                   bb.ymin() + rg.get_double()* grid_dy,
                   bb.zmin() + rg.get_double()* grid_dz)
      );
    }
  }

public:
  Split_predicate_min_volume(){}

  Split_predicate_min_volume(const double& mv, const Mesh& tmesh)
    : m_min_vol(mv)
    , m_side_of_ptr(new CGAL::Side_of_triangle_mesh<Mesh, K>(tmesh))
  {}

  /*!
    \brief returns `true` if `ni` should be split, `false` otherwise.
   */
  template <typename NodeIndex, typename Tree>
  bool operator()(NodeIndex ni, const Tree& tree) const {
    if (tree.data(ni).empty()) return false;

    OTraits::Bbox_d bb = tree.bbox(ni);
    std::vector<K::Point_3> samples;
    fill_bbox(bb, samples, 1000);

    return (estimate_volume(samples, bb.volume()) > m_min_vol);
  }
};


using Octree = CGAL::Orthtree<OTraits>;

void dump_as_polylines(const Octree& ot)
{
    std::ofstream out("octree.polylines.txt");
    for (Octree::Node_index node : ot.traverse(CGAL::Orthtrees::Leaves_traversal<Octree>(ot)))
    {
      if (!ot.is_leaf(node) || ot.data(node).empty())
        continue;
      auto bb = ot.bbox(node);
      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmin()
          << "  " << bb.xmax() << " " << bb.ymin() << " " << bb.zmin() << "\n";
      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmin()
          << "  " << bb.xmin() << " " << bb.ymax() << " " << bb.zmin() << "\n";
      out << "2 " << bb.xmax() << " " << bb.ymin() << " " << bb.zmin()
          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmin() << "\n";
      out << "2 " << bb.xmin() << " " << bb.ymax() << " " << bb.zmin()
          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmin() << "\n";
//
      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmin()
          << "  " << bb.xmin() << " " << bb.ymin() << " " << bb.zmax() << "\n";
      out << "2 " << bb.xmax() << " " << bb.ymin() << " " << bb.zmin()
          << "  " << bb.xmax() << " " << bb.ymin() << " " << bb.zmax() << "\n";
      out << "2 " << bb.xmin() << " " << bb.ymax() << " " << bb.zmin()
          << "  " << bb.xmin() << " " << bb.ymax() << " " << bb.zmax() << "\n";
      out << "2 " << bb.xmax() << " " << bb.ymax() << " " << bb.zmin()
          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmax() << "\n";
//
      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmax()
          << "  " << bb.xmax() << " " << bb.ymin() << " " << bb.zmax() << "\n";
      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmax()
          << "  " << bb.xmin() << " " << bb.ymax() << " " << bb.zmax() << "\n";
      out << "2 " << bb.xmax() << " " << bb.ymin() << " " << bb.zmax()
          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmax() << "\n";
      out << "2 " << bb.xmin() << " " << bb.ymax() << " " << bb.zmax()
          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmax() << "\n";
    }
}

int main(int argc, char** argv)
{
  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off");
  const double mv = argc>2? atof(argv[2]) : 0.00004;
  Mesh mesh;
  if(!CGAL::IO::read_polygon_mesh(filename, mesh))
  {
    std::cerr << "Error: cannot read file" << std::endl;
    return EXIT_FAILURE;
  }

  Octree tree(mesh, mesh.points());
  Split_predicate_min_volume sp(mv, mesh);
  tree.refine(sp);

  dump_as_polylines(tree);
}