4

I am currently working on a project which deals with geometric problems. Since this project will be used commercially I cannot use libraries like CGAL.

I am currently using boost::geometry with inexact types but I encountered numeric issues. I tried to simply use an exact point type from boost::multiprecision but it doesn't compile when I call boost::geometry functions.

I found this page which shows how to use a numeric_adaptor to use boost::geometry with exact number types. However, it seems outdated and I wasn't able to make it work.

Can boost::geometry be used with exact number types ? How ?

#include <vector>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp> 
#include <boost/geometry/geometries/segment.hpp>
#include <boost/geometry/algorithms/intersection.hpp>

#include <boost/multiprecision/gmp.hpp>

namespace bg = boost::geometry;
namespace bm = boost::multiprecision;

typedef bg::model::d2::point_xy<bm::mpq_rational> point;
typedef boost::geometry::model::segment<point> segment;

int main(void)
{
    point a(0,0);
    point b(1,0);
    point c(1,1);
    point d(0,1);

    segment s1(a,c);
    segment s2(b,d);

    std::vector<point> ip;
    bg::intersection(s1, s2, ip); // Doesn't compile

    return 0;
}

clang++3.4.2 output:

In file included from boost_geom_intersect.cpp:3:
In file included from /usr/include/boost/geometry.hpp:17:
In file included from /usr/include/boost/geometry/geometry.hpp:36:
In file included from /usr/include/boost/geometry/core/radian_access.hpp:21:
In file included from /usr/include/boost/numeric/conversion/cast.hpp:33:
In file included from /usr/include/boost/numeric/conversion/converter.hpp:14:
/usr/include/boost/numeric/conversion/converter_policies.hpp:187:69: error: cannot convert 'const
      boost::multiprecision::detail::expression<boost::multiprecision::detail::divides,
      boost::multiprecision::detail::expression<boost::multiprecision::detail::multiply_immediates,
      boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      void, void>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, void, void>' to 'result_type' (aka 'double') without a conversion
      operator
  static result_type low_level_convert ( argument_type s ) { return static_cast<result_type>(s) ; }
                                                                    ^~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/boost/numeric/conversion/detail/converter.hpp:524:32: note: in instantiation of member function
      'boost::numeric::raw_converter<boost::numeric::conversion_traits<double, boost::multiprecision::detail::expression<boost::multiprecision::detail::divides,
      boost::multiprecision::detail::expression<boost::multiprecision::detail::multiply_immediates,
      boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      void, void>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, void, void> > >::low_level_convert' requested here
      return RawConverterBase::low_level_convert(s);
                               ^
/usr/include/boost/numeric/conversion/cast.hpp:53:27: note: in instantiation of member function
      'boost::numeric::convdetail::non_rounding_converter<boost::numeric::conversion_traits<double,
      boost::multiprecision::detail::expression<boost::multiprecision::detail::divides,
      boost::multiprecision::detail::expression<boost::multiprecision::detail::multiply_immediates,
      boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      void, void>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, void, void> >,
      boost::numeric::convdetail::dummy_range_checker<boost::numeric::conversion_traits<double,
      boost::multiprecision::detail::expression<boost::multiprecision::detail::divides,
      boost::multiprecision::detail::expression<boost::multiprecision::detail::multiply_immediates,
      boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      void, void>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, void, void> > >,
      boost::numeric::raw_converter<boost::numeric::conversion_traits<double, boost::multiprecision::detail::expression<boost::multiprecision::detail::divides,
      boost::multiprecision::detail::expression<boost::multiprecision::detail::multiply_immediates,
      boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      void, void>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, void, void> > > >::convert' requested here
        return converter::convert(arg);
                          ^
/usr/include/boost/geometry/policies/robustness/segment_ratio.hpp:146:22: note: in instantiation of function template specialization 'boost::numeric_cast<double,
      boost::multiprecision::detail::expression<boost::multiprecision::detail::divides,
      boost::multiprecision::detail::expression<boost::multiprecision::detail::multiply_immediates,
      boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      void, void>, boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, void, void> >' requested here
            : boost::numeric_cast<double>
                     ^
/usr/include/boost/geometry/policies/robustness/segment_ratio.hpp:129:9: note: in instantiation of member function
      'boost::geometry::segment_ratio<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1> >::initialize' requested here
        initialize();
        ^
/usr/include/boost/geometry/strategies/cartesian/cart_intersect.hpp:207:33: note: in instantiation of member function
      'boost::geometry::segment_ratio<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1> >::assign' requested here
                sinfo.robust_ra.assign(robust_da, robust_da0);
                                ^
/usr/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp:114:47: note: in instantiation of function template specialization
      'boost::geometry::strategy::intersection::relate_cartesian_segments<boost::geometry::policies::relate::segments_intersection_points<boost::geometry::segment_intersection_points<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational,
      1>, boost::geometry::cs::cartesian>, boost::geometry::segment_ratio<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1> > > >,
      void>::apply<boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >,
      boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >, boost::geometry::detail::no_rescale_policy,
      boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, boost::geometry::cs::cartesian> >' requested
      here
        intersection_return_type is = policy::apply(segment1, segment2,
                                              ^
/usr/include/boost/geometry/algorithms/intersection.hpp:51:12: note: in instantiation of function template specialization
      'boost::geometry::detail::intersection::intersection_segment_segment_point<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational,
      1>, boost::geometry::cs::cartesian>
      >::apply<boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >,
      boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >, boost::geometry::detail::no_rescale_policy,
      std::back_insert_iterator<std::vector<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian>, std::allocator<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational,
      1>, boost::geometry::cs::cartesian> > > >, boost::geometry::strategy_intersection<boost::geometry::cartesian_tag,
      boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >,
      boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >, boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian>, boost::geometry::detail::no_rescale_policy, void> >' requested here
        >::apply(geometry1, geometry2, robust_policy, std::back_inserter(geometry_out), strategy);
           ^
/usr/include/boost/geometry/algorithms/intersection.hpp:148:12: note: in instantiation of function template specialization
      'boost::geometry::dispatch::intersection<boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational,
      1>, boost::geometry::cs::cartesian> >,
      boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >, boost::geometry::segment_tag, boost::geometry::segment_tag, false>::apply<boost::geometry::detail::no_rescale_policy,
      std::vector<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>, boost::geometry::cs::cartesian>,
      std::allocator<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> > >, boost::geometry::strategy_intersection<boost::geometry::cartesian_tag,
      boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >,
      boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >, boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian>, boost::geometry::detail::no_rescale_policy, void> >' requested here
        >::apply(geometry1, geometry2, robust_policy, geometry_out, strategy());
           ^
/usr/include/boost/geometry/algorithms/intersection.hpp:308:21: note: in instantiation of function template specialization
      'boost::geometry::resolve_variant::intersection<boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational,
      1>, boost::geometry::cs::cartesian> >,
      boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >
      >::apply<std::vector<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian>, std::allocator<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational,
      1>, boost::geometry::cs::cartesian> > > >' requested here
        >::template apply
                    ^
boost_geom_intersect.cpp:27:9: note: in instantiation of function template specialization
      'boost::geometry::intersection<boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational,
      1>, boost::geometry::cs::cartesian> >,
      boost::geometry::model::segment<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> >, std::vector<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational,
      1>, boost::geometry::cs::cartesian>,
      std::allocator<boost::geometry::model::d2::point_xy<boost::multiprecision::number<boost::multiprecision::backends::gmp_rational, 1>,
      boost::geometry::cs::cartesian> > > >' requested here
    bg::intersection(s1, s2, ip); // Doesn't compile
        ^
1 error generated.
deadlock
  • 93
  • 1
  • 5
  • Side note: you can totally use CGAL in a commercial project, you should discuss with GeometryFactory, the commercial licenses are affordable. Though if Boost.Geometry works for you, that's even cheaper indeed. – Marc Glisse Sep 09 '14 at 11:43

1 Answers1

3

Boost Geometry is getting confused by the proxy types returned from expression templates, where it is expecting the concrete numeric results: documentation

The Multiprecision library comes in two distinct parts:

  • An expression-template-enabled front-end number that handles all the operator overloading, expression evaluation optimization, and code reduction.
  • A selection of back-ends that implement the actual arithmetic operations, and need conform only to the reduced interface requirements of the front-end.

The meta-programming grinds to a halt there.

Fortunately, you can simply use the mpq_rational modified to disable expressions templates:

typedef bm::number<bm::gmp_rational, bm::et_off> my_rational;

This will compile without problems.

Coliru chokes on it, but here it is: http://coliru.stacked-crooked.com/a/232d98bfbb430468

#include <vector>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp> 
#include <boost/geometry/geometries/segment.hpp>
#include <boost/geometry/algorithms/intersection.hpp>

#include <boost/multiprecision/gmp.hpp>
#include <boost/multiprecision/number.hpp>

namespace bg = boost::geometry;
namespace bm = boost::multiprecision;

typedef bm::number<bm::gmp_rational, bm::et_off> my_rational;
typedef bg::model::d2::point_xy<my_rational > point;
typedef boost::geometry::model::segment<point> segment;

int main(void)
{
    point a(0,0);
    point b(1,0);
    point c(1,1);
    point d(0,1);

    segment s1(a,c);
    segment s2(b,d);

    std::vector<point> ip;
    bg::intersection(s1, s2, ip); // Doesn't compile

    return 0;
}
sehe
  • 374,641
  • 47
  • 450
  • 633
  • It works, thanks =). My running time is now 73 minutes (inexact: 10 seconds). To speed this up, I would try to use a lazy exact type (which avoids unneeded exact computations) but I didn't find any in boost. CGAL handles these types. However it looks like CGAL rational types use expression templates that cannot be deactivated (I get the same error as above when using CGAL exact number types). On the other hand, using CGAL::Lazy_exact_nt of boost::multiprecision types doesn't seem to be simple (I didn't find examples of custom number type in CGAL). Do you have any idea on how to do this? – deadlock Aug 28 '14 at 14:40
  • The original `mpq_rational` is exactly that: a lazy exact type. The problem is, Geometry does not support unevaluated expressions, so that's the end of the road there. The problem not being that there is no lazily-evaluated type, but such a type is not being supported in Boost Geo. You could get away with wrapping your own type that hides it, but this would require a fair bit of type-erasure-fu – sehe Aug 28 '14 at 14:43
  • @deadlock CGAL number types do not use expression templates (mpq_class from GMP does, but it is not used by default in CGAL), so the errors you get are probably different. – Marc Glisse Sep 09 '14 at 11:48