#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/geometries.hpp>
#include <boost/geometry/io/svg/svg_mapper.hpp>

#include <fstream>

namespace bg = boost::geometry;

const int point_count = 360; // points for a full circle

// Function to let buffer-distance depend on alpha, e.g.:
inline double corrected_distance(double distance, double alpha)
{
    return distance * 1.0 + 0.2 * sin(alpha * 6.0);
}

class buffer_point_strategy_sample
{
public :

    template
    <
        typename Point,
        typename OutputRange,
        typename DistanceStrategy
    >
    void apply(Point const& point,
            DistanceStrategy const& distance_strategy,
            OutputRange& output_range) const
    {
        double const distance = distance_strategy.apply(point, point,
                        bg::strategy::buffer::buffer_side_left);

        double const angle_increment = 2.0 * M_PI / double(point_count);

        double alpha = 0;
        for (std::size_t i = 0; i <= point_count; i++, alpha -= angle_increment)
        {
            double const cd = corrected_distance(distance, alpha);

            typename boost::range_value<OutputRange>::type output_point;
            bg::set<0>(output_point, bg::get<0>(point) + cd * cos(alpha));
            bg::set<1>(output_point, bg::get<1>(point) + cd * sin(alpha));
            output_range.push_back(output_point);
        }
    }
};

class buffer_join_strategy_sample
{
private :
    template
    <
        typename Point,
        typename DistanceType,
        typename RangeOut
    >
    inline void generate_points(Point const& vertex,
                Point const& perp1, Point const& perp2,
                DistanceType const& buffer_distance,
                RangeOut& range_out) const
    {
        double dx1 = bg::get<0>(perp1) - bg::get<0>(vertex);
        double dy1 = bg::get<1>(perp1) - bg::get<1>(vertex);
        double dx2 = bg::get<0>(perp2) - bg::get<0>(vertex);
        double dy2 = bg::get<1>(perp2) - bg::get<1>(vertex);

        // Assuming the corner is convex, angle2 < angle1
        double const angle1 = atan2(dy1, dx1);
        double angle2 = atan2(dy2, dx2);

        while (angle2 > angle1)
        {
            angle2 -= 2 * M_PI;
        }

        double const angle_increment = 2.0 * M_PI / double(point_count);
        double alpha = angle1 - angle_increment;

        for (int i = 0; alpha >= angle2 && i < point_count; i++, alpha -= angle_increment)
        {
            double cd = corrected_distance(buffer_distance, alpha);

            Point p;
            bg::set<0>(p, bg::get<0>(vertex) + cd * cos(alpha));
            bg::set<1>(p, bg::get<1>(vertex) + cd * sin(alpha));
            range_out.push_back(p);
        }
    }

public :

    template <typename Point, typename DistanceType, typename RangeOut>
    inline bool apply(Point const& ip, Point const& vertex,
                Point const& perp1, Point const& perp2,
                DistanceType const& buffer_distance,
                RangeOut& range_out) const
    {
        generate_points(vertex, perp1, perp2, buffer_distance, range_out);
        return true;
    }

    template <typename NumericType>
    static inline NumericType max_distance(NumericType const& distance)
    {
        return distance;
    }

};

class buffer_side_sample
{
public :
    template
    <
        typename Point,
        typename OutputRange,
        typename DistanceStrategy
    >
    static inline void apply(
                Point const& input_p1, Point const& input_p2,
                bg::strategy::buffer::buffer_side_selector side,
                DistanceStrategy const& distance,
                OutputRange& output_range)
    {
        // Generate a block along (left or right of) the segment

        double const dx = bg::get<0>(input_p2) - bg::get<0>(input_p1);
        double const dy = bg::get<1>(input_p2) - bg::get<1>(input_p1);

        // For normalization [0,1] (=dot product d.d, sqrt)
        double const length = bg::math::sqrt(dx * dx + dy * dy);

        if (bg::math::equals(length, 0))
        {
            return;
        }

        // Generate the normalized perpendicular p, to the left (ccw)
        double const px = -dy / length;
        double const py = dx / length;

        // Both vectors perpendicular to input p1 and input p2 have same angle
        double const alpha = atan2(py, px);

        double const d = distance.apply(input_p1, input_p2, side);

        double const cd = corrected_distance(d, alpha);

        output_range.resize(2);

        bg::set<0>(output_range.front(), bg::get<0>(input_p1) + px * cd);
        bg::set<1>(output_range.front(), bg::get<1>(input_p1) + py * cd);
        bg::set<0>(output_range.back(), bg::get<0>(input_p2) + px * cd);
        bg::set<1>(output_range.back(), bg::get<1>(input_p2) + py * cd);
    }
};


template <typename Geometry1, typename Geometry2>
void create_svg(std::string const& filename, Geometry1 const& original, Geometry2 const& buffer, std::string const& color)
{
    typedef typename bg::point_type<Geometry1>::type point_type;
    std::ofstream svg(filename.c_str());

    bg::svg_mapper<point_type> mapper(svg, 800, 800);
    mapper.add(buffer);

    mapper.map(original, "fill-opacity:0.3;fill:rgb(255,0,0);stroke:rgb(0,0,0);stroke-width:1");

    std::string style = "fill-opacity:0.3;fill:";
    style += color;
    style += ";stroke:rgb(0,0,0);stroke-width:1";
    mapper.map(buffer, style);
}

int main()
{
    typedef bg::model::d2::point_xy<double> point;
    typedef bg::model::polygon<point> polygon;

    // Predefined strategies
    bg::strategy::buffer::distance_symmetric<double> distance_strategy(1.3);
    bg::strategy::buffer::end_flat end_strategy; // not effectively used

    // Own strategies
    buffer_join_strategy_sample join_strategy;
    buffer_point_strategy_sample point_strategy;
    buffer_side_sample side_strategy;

    // Declare output
    bg::model::multi_polygon<polygon> result;

    // Declare/fill a multi-point
    bg::model::multi_point<point> mp;
    bg::read_wkt("MULTIPOINT((5 5),(7 5),(9 5),(7 10))", mp);

    // Create the buffer of a multi-point
    bg::buffer(mp, result,
                distance_strategy, side_strategy,
                join_strategy, end_strategy, point_strategy);

    create_svg("/tmp/buffer.svg", mp, result, "rgb(51,51,153)");

    // Polygon
    bg::model::multi_polygon<polygon> mpol;
    bg::read_wkt("MULTIPOLYGON(((0 10,1 15,6 11,0 10)),((2 0,3 8,4 0,2 0)))", mpol);
    result.clear();
    bg::buffer(mpol, result,
                distance_strategy, side_strategy,
                join_strategy, end_strategy, point_strategy);

    create_svg("/tmp/bufferpol.svg", mpol, result, "rgb(51,51,153)");

    return 0;
}