CGAL / cgal

The public CGAL repository, see the README below
https://github.com/CGAL/cgal#readme
Other
4.86k stars 1.38k forks source link

Can Nef_polyhedron_3 (and minkowski_sum_3 ) be used with inexact kernels? #5490

Open ochafik opened 3 years ago

ochafik commented 3 years ago

Hi CGAL team!

I've been struggling to use Epick (or Cartesian<double>) with Nef_polyhedron_3: getting failures with even quite simple polyhedra, either during creation of the nef itself or during calls to minkowski_sum_3. The same models work great w/ Epeck and other exact kernels. Is that expected?

(also, not sure if relevant, but Polyhedron_3::is_valid(false, /*level*/ 1) says these models aren't valid but level 0 says they are?)

Context: I'm trying to find ways to speed up OpenSCAD's minkowski operations, which only use CGAL's minkowski with an exact kernel as a fallback. Their default algorithm is to decompose each operand in convex parts and union the hulls of the cross products of these respective parts. Since these hulls operate with Epick, I figured, why not be inexact all the way and give minkowski_sum_3 another try? But it doesn't seem to work out of the box :-(

I've pasted a self-contained reproduction case below, which can be compiled & run with:

# Add -DEXACT to use Epeck instead of Epick: it then all runs fine!
# Add -DMESH to use Surface_mesh instead of Polyhedron_3: same issue.
# Also tried the various CGAL_NO_ASSERTIONS and such: doesn't do any good.

g++ -I../CGAL-5.2/include -stdlib=libc++ -std=c++1y -lgmp invalid_poly.cc -o invalid_poly

./invalid_poly 1
./invalid_poly 2

Resulting output:

# Weird: polyhedron is 'invalid' as per Polyhedron_3::is_valid(level=1)!
# Creating nef...
CGAL error: assertion violation!
Expression : ss_circle.has_on(sp)
File       : ../CGAL-5.2/include/CGAL/Nef_3/polygon_mesh_to_nef_3.h
Line       : 255
Explanation: 
Refer to the bug-reporting instructions at https://www.cgal.org/bug_report.html
libc++abi.dylib: terminating with uncaught exception of type CGAL::Assertion_exception: CGAL ERROR: assertion violation!
Expr: ss_circle.has_on(sp)
File: ../CGAL-5.2/include/CGAL/Nef_3/polygon_mesh_to_nef_3.h
Line: 255
# Weird: polyhedron is 'invalid' as per Polyhedron_3::is_valid(level=1)!
# Running minkowski_sum_3...
CGAL error: assertion violation!
Expression : orientation(pred, p_sweep) == 0
File       : ../CGAL-5.2/include/CGAL/Nef_2/Segment_overlay_traits.h
Line       : 735
Explanation: 
Refer to the bug-reporting instructions at https://www.cgal.org/bug_report.html
libc++abi.dylib: terminating with uncaught exception of type CGAL::Assertion_exception: CGAL ERROR: assertion violation!
Expr: orientation(pred, p_sweep) == 0
File: ../CGAL-5.2/include/CGAL/Nef_2/Segment_overlay_traits.h
Line: 735

Source Code

invalid_poly.cc:

// Copyright 2021 Google LLC.
// SPDX-License-Identifier: Apache-2.0
/*
  g++ -I../CGAL-5.2/include -stdlib=libc++ -std=c++1y -lgmp invalid_poly.cc -o invalid_poly && \
      ./invalid_poly 1
*/

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Nef_polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/minkowski_sum_3.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/Surface_mesh.h>
#include <fstream>

extern std::string example1; // see bottom of file
extern std::string example2;

int main(int argc, char *argv[])
{
#ifdef EXACT
  typedef CGAL::Epeck Kernel;
#else
  typedef CGAL::Epick Kernel;
#endif
  typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
  typedef CGAL::Surface_mesh<CGAL::Point_3<Kernel>> Mesh;
  typedef CGAL::Nef_polyhedron_3<Kernel> Nef;

#ifdef MESH
  Mesh polyhedron;
#else
  Polyhedron polyhedron;
#endif

  std::stringstream ss;
  ss << ((argc > 1 && std::string(argv[1]) == "2") ? example2 : example1);
  ss >> polyhedron;

  assert(CGAL::is_closed(polyhedron)); // OK
  assert(CGAL::is_valid_halfedge_graph(polyhedron, /* verb */ false)); // OK

#ifndef MESH
  assert(polyhedron.is_valid(/* verbose */ false, 0)); // OK
  if (!polyhedron.is_valid(/* verbose */ false, 1)) { // NOT OK?
    std::cerr << "# Weird: polyhedron is 'invalid' as per Polyhedron_3::is_valid(level=1)!\n";
  }
#endif

  std::cout << "# Creating nef...\n";
  Nef lhs(polyhedron); // Chokes on example 1
  Nef rhs = lhs;
  std::cout << "# nef has " << lhs.number_of_facets() << " facets\n";

  std::cout << "# Running minkowski_sum_3...\n";
  auto result = CGAL::minkowski_sum_3(lhs, rhs); // Chokes on example 2

  // Polyhedron resultPoly;
  // result.convert_to_polyhedron(resultPoly);
  // std::ofstream("out.off") << resultPoly;
  // std::cout << "out.off: " << result.number_of_facets() << " facets\n";

  std::cout << "# All gooood!\n"; // Never getting here
  return EXIT_SUCCESS;
}

// `sphere(1, $fn=8);` in OpenSCAD
std::string example1(
"OFF\n\
32 60 0\n\
0.270598 0.270598 0.92388\n\
0 0.382683 0.92388\n\
-0.270598 0.270598 0.92388\n\
-0.382683 0 0.92388\n\
-0.270598 -0.270598 0.92388\n\
0 -0.382683 0.92388\n\
0.270598 -0.270598 0.92388\n\
0.382683 0 0.92388\n\
0.653281 0.653281 0.382683\n\
0 0.92388 0.382683\n\
-0.653281 0.653281 0.382683\n\
-0.92388 0 0.382683\n\
-0.653281 -0.653281 0.382683\n\
0 -0.92388 0.382683\n\
0.653281 -0.653281 0.382683\n\
0.92388 0 0.382683\n\
0.653281 0.653281 -0.382683\n\
0 0.92388 -0.382683\n\
-0.653281 0.653281 -0.382683\n\
-0.92388 0 -0.382683\n\
-0.653281 -0.653281 -0.382683\n\
0 -0.92388 -0.382683\n\
0.653281 -0.653281 -0.382683\n\
0.382683 0 -0.92388\n\
0.92388 0 -0.382683\n\
0.270598 0.270598 -0.92388\n\
0 0.382683 -0.92388\n\
-0.270598 0.270598 -0.92388\n\
-0.382683 0 -0.92388\n\
-0.270598 -0.270598 -0.92388\n\
0 -0.382683 -0.92388\n\
0.270598 -0.270598 -0.92388\n\
3  7 15 8\n\
3  7 8 0\n\
3  0 8 9\n\
3  0 9 1\n\
3  1 9 10\n\
3  1 10 2\n\
3  2 10 11\n\
3  2 11 3\n\
3  3 11 12\n\
3  3 12 4\n\
3  4 12 13\n\
3  4 13 5\n\
3  5 13 14\n\
3  5 14 6\n\
3  6 14 15\n\
3  6 15 7\n\
3  15 24 16\n\
3  15 16 8\n\
3  8 16 17\n\
3  8 17 9\n\
3  9 17 18\n\
3  9 18 10\n\
3  10 18 19\n\
3  10 19 11\n\
3  11 19 20\n\
3  11 20 12\n\
3  12 20 21\n\
3  12 21 13\n\
3  13 21 22\n\
3  13 22 14\n\
3  14 22 24\n\
3  14 24 15\n\
3  24 23 25\n\
3  24 25 16\n\
3  16 25 26\n\
3  16 26 17\n\
3  17 26 27\n\
3  17 27 18\n\
3  18 27 28\n\
3  18 28 19\n\
3  19 28 29\n\
3  19 29 20\n\
3  20 29 30\n\
3  20 30 21\n\
3  21 30 31\n\
3  21 31 22\n\
3  22 31 23\n\
3  22 23 24\n\
3  0 1 7\n\
3  5 1 3\n\
3  3 1 2\n\
3  5 3 4\n\
3  7 1 5\n\
3  7 5 6\n\
3  30 29 28\n\
3  26 30 28\n\
3  27 26 28\n\
3  25 23 26\n\
3  23 30 26\n\
3  31 30 23\n\
");

// Derived from `sphere(1, $fn=3);` in OpenSCAD (edited to have integer coordinates, but same issue as the original + same general shape)
std::string example2(
"OFF\n\
6 8 0\n\
-3.0 6.0 7.0\n\
-3.0 -6.0 7.0\n\
7.0 0 -7.0\n\
7.0 0 7.0\n\
-3.0 6.0 -7.0\n\
-3.0 -6.0 -7.0\n\
3  1 3 0\n\
3  3 2 4\n\
3  3 4 0\n\
3  0 4 5\n\
3  0 5 1\n\
3  1 5 2\n\
3  1 2 3\n\
3  2 5 4\n\
");

Environment

sloriot commented 3 years ago

Hi Olivier,

Nef (and Minkowski) requires exact constructions so it cannot work with floating point based number type.

ochafik commented 3 years ago

Hey Sebastien, thanks for confirming!

It's a bit far-fetched, but wondering how far I could go writing my own approximate kernel where a == b allows for some epsilon difference... do you know if anything of the sort has been attempted, or how doomed to fail it would be?

Also, exploring other leads, I was wondering if I could maybe use some approximate convex decomposition to limit the worst cases of many convex parts in OpenSCAD's custom minkowski algorithm. I came across https://github.com/CGAL/cgal/issues/3270, and I'm not sure whether the project was integrated into CGAL proper: the Triangulated Surface Mesh Segmentation looks close but not quite the same. Would you know how it compares / if there is hope to get the GSOC version for general use?

sloriot commented 3 years ago

Approximate won't work because internally it assumes that all the vertices of a face are coplanar which is no longer the case as soon as you are using a kernel with inexact constructions.

About the approximate convex decomposition, the PR is still opened https://github.com/CGAL/cgal/pull/3302 as I have to investigate why it is that slow (compared to the author's version).

ochafik commented 3 years ago

Approximate won't work because internally it assumes that all the vertices of a face are coplanar which is no longer the case as soon as you are using a kernel with inexact constructions.

I would be okay with lots of inaccuracies piling up and some failed assertions down the line forcing me to fall back to an exact kernel. Do you think the lack of assertions being thrown would be a good proxy to evaluate whether the nef / minkowski sum succeeded against all odds? A fast (best effort), approximate minkowski could potentially be great to have in OpenSCAD, if only for previews maybe.

the PR is still opened #3302 as I have to investigate why it is that slow (compared to the author's version).

Hah, cool, good luck! Also, even a suboptimal version could have great practical use, hope it lands!

ochafik commented 3 years ago

approximate kernel where a == b allows for some epsilon difference.

FWIW, I've taken a stab at making such a kernel, just to convince myself it was a bad idea. It's not pretty, it builds nefs alright on a handful of trivial models (Cartesian<FuzzyDouble>), even agrees to do nef csg ops, but... explodes in Convex_decomposition_3. Oh well, you warned me 🤷‍♂️.

// Copyright 2021 Google LLC.
// SPDX-License-Identifier: Apache-2.0
#pragma once

#include <cmath>
#include <type_traits>

// Field type for an Ipick: inexact predicates (epsilon), inexact construction (double) kernel
class FuzzyDouble {
public:
  static constexpr double epsilon = 0.00000001;
  double value;

  template <typename T, std::enable_if_t<std::is_arithmetic<T>::value, bool> = true>
  FuzzyDouble(T value) : value(static_cast<double>(value)) {}
  FuzzyDouble(const FuzzyDouble& other) { *this = other; }
  FuzzyDouble() : value(0.0) {}

  operator double() const { return value; }
  FuzzyDouble &operator=(const FuzzyDouble& other) { value = other.value; return *this; }
  FuzzyDouble operator-() const { return -value; }

#define OP_TEMPLATE(T) \
  template <typename T, std::enable_if_t<(std::is_arithmetic<T>::value || std::is_assignable<FuzzyDouble, T>::value), bool> = true>

  OP_TEMPLATE(T) bool operator==(const T& x) const { return std::abs(value - static_cast<double>(x)) <= FuzzyDouble::epsilon;}
  OP_TEMPLATE(T) bool operator!=(const T& x) const { return !(*this == x); }
  OP_TEMPLATE(T) bool operator<(const T& x) const { return value + FuzzyDouble::epsilon < static_cast<double>(x); }
  OP_TEMPLATE(T) bool operator<=(const T& x) const { return value <= static_cast<double>(x) + FuzzyDouble::epsilon; }
  OP_TEMPLATE(T) bool operator>(const T& x) const { return value > static_cast<double>(x) + FuzzyDouble::epsilon; }
  OP_TEMPLATE(T) bool operator>=(const T& x) const { return value + FuzzyDouble::epsilon >= static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble operator*(const T& x) const { return value * static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble operator/(const T& x) const { return value / static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble operator+(const T& x) const { return value + static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble operator-(const T& x) const { return value - static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble& operator*=(const T& x) { value *= static_cast<double>(x); return *this; }
  OP_TEMPLATE(T) FuzzyDouble& operator/=(const T& x) { value /= static_cast<double>(x); return *this; }
  OP_TEMPLATE(T) FuzzyDouble& operator+=(const T& x) { value += static_cast<double>(x); return *this; }
  OP_TEMPLATE(T) FuzzyDouble& operator-=(const T& x) { value -= static_cast<double>(x); return *this; }
};

std::ostream& operator<<(std::ostream& out, const FuzzyDouble& n) { out << n.value; return out; }
std::istream& operator>>(std::istream& in, FuzzyDouble& n) { in >> n.value; return in; }

namespace std { FuzzyDouble sqrt(const FuzzyDouble& x) { return std::sqrt(x.value); } }

namespace CGAL {
  inline FuzzyDouble min BOOST_PREVENT_MACRO_SUBSTITUTION(const FuzzyDouble& x,const FuzzyDouble& y){ return std::min(x.value, y.value); }
  inline FuzzyDouble max BOOST_PREVENT_MACRO_SUBSTITUTION(const FuzzyDouble& x,const FuzzyDouble& y){ return std::max(x.value, y.value); }

  CGAL_DEFINE_COERCION_TRAITS_FOR_SELF(FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(short    ,FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(int      ,FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(long     ,FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(float    ,FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(double   ,FuzzyDouble)

  template <>
  class Algebraic_structure_traits<FuzzyDouble> : public Algebraic_structure_traits_base<FuzzyDouble,Field_with_sqrt_tag> {
  public:
    typedef Tag_false       Is_exact;
    typedef Tag_true        Is_numerical_sensitive;
  };

  template <>
  class Real_embeddable_traits<FuzzyDouble>: public INTERN_RET::Real_embeddable_traits_base<FuzzyDouble,CGAL::Tag_true>
  {
    typedef Algebraic_structure_traits<Type>        AST;
  public:

    typedef Tag_true                                Is_real_embeddable;
    typedef Uncertain<bool>                         Boolean;
    typedef Uncertain<CGAL::Comparison_result>      Comparison_result;
    typedef Uncertain<CGAL::Sign>                   Sign;

    typedef AST::Is_zero    Is_zero;

    struct Is_finite: public CGAL::cpp98::unary_function<Type,Boolean>{
      inline Boolean operator()(const Type &x) const { return !std::isinf(x.value) && !std::isnan(x.value); };
    };

    struct Abs: public CGAL::cpp98::unary_function<Type,Type>{
      inline Type operator()(const Type &x) const { return std::abs(x.value); };
    };

    struct To_interval: public CGAL::cpp98::unary_function<Type,std::pair<double,double> >{
      inline std::pair<double,double> operator()(const Type &x)const{
        return std::make_pair(x.value - FuzzyDouble::epsilon, x.value + FuzzyDouble::epsilon);
      };
    };
  };

} //namespace CGAL