Буст, геометрия

Вы должны использовать интегральные координаты.

Из документации: http://www.boost.org/doc/libs/1_58_0/libs/polygon/doc/index.htm

Тип данных координат является параметром шаблона всех типов данных и алгоритмов, предоставляемых библиотекой, и ожидается, что он будет интегральным. Типы данных координат с плавающей запятой не поддерживаются алгоритмами, реализованными в библиотеке, из-за того, что достижение надежности с плавающей запятой подразумевает другой набор алгоритмов и, как правило, специфичные для платформы предположения о представлениях с плавающей запятой.

То же самое относится и к более ранним версиям.

В вашем случае вывод Boost(1_55 или 1_52) правильный (случайно).

Вывод может показаться правильным на первый взгляд, но на самом деле это не так, если присмотреться:

Слегка измененный пример: Live On Coliru

#include <boost/geometry.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/array.hpp>
#include <vector>
#include <iomanip>

namespace bg = boost::geometry;

template <typename C, typename T = typename C::value_type>
    void append(C& container, std::initializer_list<T> init) {
        container.insert(container.end(), init);
    }

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

    std::vector<TBoostPoly> squares;

    using Eps = boost::array<double, 4>;

    for (auto const& eps : { 
                Eps {{     0,     0,     0,     0 }},
                Eps {{ 1e-15, 1e-15, 2e-15, 2e-15 }},
            })
    {
        TBoostPoly square;
        append(square.outer(), {
             { 0.5,  4.25 + eps[0] },
             { 0.5,   4.5 + eps[1] },
             { 1.0,   4.5 + eps[2] },
             { 1.0,  4.25 + eps[3] },
             { 0.5,  4.25 + eps[0] }
            });

        squares.push_back(std::move(square));
    }

    for (auto& p : squares)
        bg::correct(p);

    std::vector<TBoostPoly> output;
    bg::intersection(squares[0], squares[1], output);

    for (auto& p : output) std::cout << "Output: " << bg::wkt(p) << "\n";
    std::cout << std::fixed << std::setprecision(std::numeric_limits<bg::coordinate_type<TBoostPoint>::type >::max_digits10);
    for (auto& p : output) std::cout << "Output: " << bg::wkt(p) << "\n";
}

Какие печатает

Output: POLYGON((0.5 4.5,1 4.5,1 4.25,0.5 4.25,0.5 4.5))
Output: POLYGON((0.50000000000000000 4.50000000000000000,1.00000000000000000 4.50000000000000000,1.00000000000000000 4.25000000000000178,0.50000000000000000 4.25000004999999970,0.50000000000000000 4.50000000000000000))

Как видите, наивный, естественный вывод может казаться равным 4.25 в какой-то момент, но фактическое сохраненное значение равно 4.25000000000000178 именно в этот момент.

В зависимости от того, какие у вас есть требования, вы можете быть счастливее с некоторыми типами десятичного представления произвольной точности. В качестве доказательства концепции, вот та же программа, параметризованная для использования 50-значных десятичных чисел с плавающей запятой:

Жить на Coliru

#include <boost/geometry.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/array.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <vector>
#include <iomanip>

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

template <typename C, typename T = typename C::value_type>
    void append(C& container, std::initializer_list<T> init) {
        container.insert(container.end(), init);
    }

int main() {
    typedef bmp::number<bmp::cpp_dec_float<50>, bmp::et_off> Decimal;
    typedef bg::model::d2::point_xy<Decimal> TBoostPoint;
    typedef bg::model::polygon<TBoostPoint> TBoostPoly;

    std::vector<TBoostPoly> squares;

    using Eps = boost::array<Decimal, 4>;

    for (auto const& eps : { 
                Eps {{     0,     0,     0,     0 }},
                Eps {{ 1e-15, 1e-15, 2e-15, 2e-15 }},
            })
    {
        TBoostPoly square;
        append(square.outer(), {
             { 0.5,  4.25 + eps[0] },
             { 0.5,   4.5 + eps[1] },
             { 1.0,   4.5 + eps[2] },
             { 1.0,  4.25 + eps[3] },
             { 0.5,  4.25 + eps[0] }
            });

        squares.push_back(std::move(square));
    }

    for (auto& p : squares)
        bg::correct(p);

    std::vector<TBoostPoly> output;
    bg::intersection(squares[0], squares[1], output);

    for (auto& p : output) std::cout << "Output: " << bg::wkt(p) << "\n";
    std::cout << std::fixed << std::setprecision(std::numeric_limits<bg::coordinate_type<TBoostPoint>::type >::max_digits10);
    for (auto& p : output) std::cout << "Output: " << bg::wkt(p) << "\n";
}

Что печатает:

Output: POLYGON((0.5 4.5,1 4.5,1 4.25,0.5 4.25,0.5 4.5))
Output: POLYGON((0.50000000000000000000000000000000000000000000000000000000000000000000000000000000 4.50000000000000000000000000000000000000000000000000000000000000000000000000000000,1.00000000000000000000000000000000000000000000000000000000000000000000000000000000 4.50000000000000000000000000000000000000000000000000000000000000000000000000000000,1.00000000000000000000000000000000000000000000000000000000000000000000000000000000 4.25000000000000200000000000000015541079975332215847661437120239003029098500000000,0.50000000000000000000000000000000000000000000000000000000000000000000000000000000 4.25000000000000100000000000000007770539987666107923830718560119501514549200000000,0.50000000000000000000000000000000000000000000000000000000000000000000000000000000 4.50000000000000000000000000000000000000000000000000000000000000000000000000000000))