#include <iostream>
#include <vector>

#define BOOST_NO_CXX11_NOEXCEPT

#include <boost/multiprecision/gmp.hpp>
//#include <boost/multiprecision/cpp_int.hpp>
//#include <boost/multiprecision/number.hpp>
//#include <boost/numeric/conversion/converter.hpp>

//typedef boost::multiprecision::mpq_rational GeometryCoord;
//typedef boost::multiprecision::cpp_rational GeometryCoord;
//typedef boost::numeric::conversion_traits<double , GeometryCoord> MyFloat_to_double_Traits;

/*namespace boost {
  namespace numeric {
template<>
struct raw_converter<MyFloat_to_double_Traits>
{
  static double low_level_convert ( GeometryCoord const&  s ) { return( 0.0 ); } // Must be real conversion
};
}
}*/

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/io/wkt/wkt.hpp>

#include <boost/foreach.hpp>

typedef boost::multiprecision::number<boost::multiprecision::gmp_rational, boost::multiprecision::et_off> GeometryCoord;
//typedef boost::multiprecision::number<boost::multiprecision::gmp_rational> GeometryCoord;
//typedef boost::multiprecision::mpq_rational GeometryCoord;
//typedef boost::multiprecision::number<boost::multiprecision::cpp_rational, boost::multiprecision::et_off> GeometryCoord;

typedef boost::geometry::model::d2::point_xy<GeometryCoord> GeometryPoint;
typedef boost::geometry::model::multi_point<GeometryPoint> GeometryMultiPoint;
typedef boost::geometry::model::segment<GeometryPoint> GeometrySegment;
typedef boost::geometry::model::box<GeometryPoint> GeometryBox;
typedef boost::geometry::model::polygon<GeometryPoint> GeometryPolygon;

int main()
{
    /*std::vector<GeometryPoint> points;

    GeometrySegment Segment1(GeometryPoint(0, 3), GeometryPoint(3, 2));
    GeometrySegment Segment2(GeometryPoint(2, 0), GeometryPoint(0, 5));

    boost::geometry::intersection(Segment1, Segment2, points);*/
/*
    for (std::vector<point>::const_iterator i = points.begin(); i != points.end(); ++i) {
        point p = *i;

	std::cout << "x = " << p.x() << " y = " << p.y() << std::endl;
    }

    std::cout << std::endl;*/

    GeometryCoord c(3, 4);

    std::cout << " c = " << c << std::endl;

    GeometryPolygon green, blue;

    /*boost::geometry::read_wkt(
        "POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))", green);   

    boost::geometry::read_wkt(
    "POLYGON((10 10, 10 20, 20 20, 20 10, 10 10))", blue);*/




    /*boost::geometry::read_wkt(
        "POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))", green);   

    boost::geometry::read_wkt(
    "POLYGON((0 0, 0 20, 20 20, 20 0, 0 0))", blue);*/

    /*boost::geometry::read_wkt(
        "POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", green);   
	
    boost::geometry::read_wkt(
    "POLYGON((1 1, 1 2, 2 2, 2 1, 1 1))", blue);*/

    green.outer().push_back(GeometryPoint(GeometryCoord(0), GeometryCoord(0)));
    green.outer().push_back(GeometryPoint(GeometryCoord(0), GeometryCoord(101,10)));
    green.outer().push_back(GeometryPoint(GeometryCoord(101, 10), GeometryCoord(101, 10)));
    green.outer().push_back(GeometryPoint(GeometryCoord(101,10), GeometryCoord(0)));
    green.outer().push_back(GeometryPoint(GeometryCoord(0), GeometryCoord(0)));

    {
      const boost::geometry::model::ring<GeometryPoint> &ring = green.outer();

      for(boost::geometry::model::ring<GeometryPoint>::const_iterator i = ring.begin(); i != ring.end(); ++i ) {
        GeometryPoint Point = *i;

	std::cout << "x = " << Point.x() << " y = " << Point.y() << std::endl;
      }
    }

    std::cout << std::endl;

    blue.outer().push_back(GeometryPoint(GeometryCoord(10), GeometryCoord(10)));
    blue.outer().push_back(GeometryPoint(GeometryCoord(10), GeometryCoord(20)));
    blue.outer().push_back(GeometryPoint(GeometryCoord(20), GeometryCoord(20)));
    blue.outer().push_back(GeometryPoint(GeometryCoord(20), GeometryCoord(10)));
    blue.outer().push_back(GeometryPoint(GeometryCoord(10), GeometryCoord(10)));

    {
      const boost::geometry::model::ring<GeometryPoint> &ring = blue.outer();

      for(boost::geometry::model::ring<GeometryPoint>::const_iterator i = ring.begin(); i != ring.end(); ++i ) {
        GeometryPoint Point = *i;

	std::cout << "x = " << Point.x() << " y = " << Point.y() << std::endl;
      }
    }

    std::cout << std::endl;

    std::vector<GeometryPolygon> output;
    //boost::geometry::union_(green, blue, output);

    boost::geometry::intersection(green, blue, output);

    //boost::geometry::correct(green);

    std::cout << "size = " << output.size() << std::endl;

    //std::cout << "Intersects = " << boost::geometry::intersects(green, blue) << std::endl;

    BOOST_FOREACH(GeometryPolygon const& p, output) {
      //std::cout << boost::geometry::wkt<polygon>(p) << std::endl;

      const boost::geometry::model::ring<GeometryPoint> &ring = p.outer();

      for(boost::geometry::model::ring<GeometryPoint>::const_iterator i = ring.begin(); i != ring.end(); ++i ) {
        GeometryPoint Point = *i;

	std::cout << "x = " << Point.x() << " y = " << Point.y() << std::endl;
      }

      //std::cout << "Point inside: " << boost::geometry::within( point(0, 0), p ) << std::endl;
    }

    /*GeometryPoint p1(GeometryCoord(-1), GeometryCoord(-1));
    GeometryPoint p2(GeometryCoord(1), GeometryCoord(1));

    boost::geometry::distance(p1, p2);*/

    /*boost::geometry::strategy::buffer::distance_symmetric<GeometryCoord> distance_strategy(10);
    boost::geometry::strategy::buffer::join_round join_strategy;
    boost::geometry::strategy::buffer::end_round end_strategy;
    boost::geometry::strategy::buffer::point_circle point_strategy;
    boost::geometry::strategy::buffer::side_straight side_strategy;

    boost::geometry::model::multi_polygon<GeometryPolygon> result;

    boost::geometry::buffer(green, result, distance_strategy, side_strategy, join_strategy, end_strategy, point_strategy);

    if (result.size() != 1) {
      std::cout << "result.size() != 1" << std::endl;
      }*/

    GeometryMultiPoint MultiPoint;

    const boost::geometry::model::ring<GeometryPoint> &ring = blue.outer();

    for(boost::geometry::model::ring<GeometryPoint>::const_iterator i = ring.begin(); i != ring.end(); ++i ) {
      GeometryPoint Point = *i;

      boost::geometry::append(MultiPoint, Point);
    }

    GeometryPolygon Hull;

    boost::geometry::convex_hull(MultiPoint, Hull);

    boost::geometry::index::rtree< std::pair< GeometryBox, const GeometryPolygon * >, boost::geometry::index::quadratic<16, 4> > rtree;

    GeometryBox Box( GeometryPoint( 10, 10 ), GeometryPoint( 11, 11 ) );
    boost::geometry::envelope( green, Box );
    rtree.insert(std::make_pair(Box, &green));

    GeometryPoint center;

    boost::geometry::centroid( green, center );

    std::vector< std::pair< GeometryBox, const GeometryPolygon * > > result;

    GeometryBox searchBox( GeometryPoint( GeometryCoord( center.x() - 1 ), GeometryCoord( center.y() - 1 ) ),
                           GeometryPoint( GeometryCoord( center.x() + 1 ), GeometryCoord( center.y() + 1 ) ) );

    rtree.query( boost::geometry::index::intersects( searchBox ), std::back_inserter( result ) );

    std::cout << "Found " << result.size() << " polygone(s)" << std::endl;

    return 0;
}