Вы просто хотите пересечь оригинал region_obj
с нижней половиной его envelope_box
.
Вместо этого, находя только верхнюю линию пересечения и принудительно удаляяуказывает на случайно внутри полигона, вы создавали недопустимый полигон:
Live On Coliru
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/linestring.hpp>
#include <boost/geometry/geometries/multi_polygon.hpp>
#include <boost/geometry/algorithms/envelope.hpp>
#include <boost/geometry/algorithms/intersection.hpp>
#include <boost/geometry/io/io.hpp>
#include <iostream>
namespace bg = boost::geometry;
using bg_point_type = bg::model::d2::point_xy<double>;
using bg_polygon_type = bg::model::ring<bg_point_type>;
template <typename G>
void diags(std::string name, G& geom) {
std::cout << name << ": " << bg::wkt(geom) << "\n";
std::string reason;
if (!bg::is_valid(geom, reason)) {
std::cout << name << " invalid: " << reason << "\n";
bg::correct(geom);
std::cout << name << " corrected: " << bg::wkt(geom) << "\n";
if (!bg::is_valid(geom, reason)) {
std::cout << name << " invalid: " << reason << "\n";
}
}
}
int main() {
double const ratio = 0.5;
//Find the envelope of the region of the object
bg_polygon_type obj_region;
bg::read_wkt("POLYGON((0.1 0.1,0.1 0.5,0.5 0.5,0.5 0.1,0.1 0.1))", obj_region);
diags("Input", obj_region);
bg::model::box<bg_point_type> envelope_box;
bg::envelope(obj_region, envelope_box);
auto const top_y = envelope_box.min_corner().y();
auto const bottom_y = envelope_box.max_corner().y();
auto const left = envelope_box.min_corner().x();
auto const right = envelope_box.max_corner().x();
//with min_corner and max_corner, we can know the top of the ring
auto const line_y = bottom_y - (bottom_y - top_y) * ratio;
//find the intersection of the line and the polygon
bg::model::linestring<bg_point_type> const line{{left, line_y}, {right, line_y}};
bg::model::multi_point<bg_point_type> intersect;
bg::intersection(obj_region, line, intersect);
std::cout << bg::wkt(intersect) << "\n";
//remove those points higher than top
bg_polygon_type result = obj_region;
auto it = std::remove_if(std::begin(result ), std::end(result ), [=](auto const &pt) { return pt.y() < line_y; });
result.erase(it, std::end(result ));
//insert the intersect points into the results
std::copy(std::begin(intersect), std::end(intersect), std::back_inserter(result));
diags("Result", result);
}
Печать
Input: POLYGON((0.1 0.1,0.1 0.5,0.5 0.5,0.5 0.1,0.1 0.1))
MULTIPOINT((0.1 0.3),(0.5 0.3))
Result: POLYGON((0.1 0.5,0.5 0.5,0.1 0.3,0.5 0.3))
Result invalid: Geometry is defined as closed but is open
Result corrected: POLYGON((0.1 0.5,0.5 0.5,0.1 0.3,0.5 0.3,0.1 0.5))
Result invalid: Geometry has invalid self-intersections. A self-intersection point was found at (0.3, 0.4); method: i; operations: u/i; segment IDs {source, multi, ring, segment}: {0, -1, -1, 1}/{0, -1, -1, 3}
Исправление
Просто спросив, что вы имеете в виду:
bg_polygon_type obj_region;
bg::read_wkt("POLYGON((0.1 0.1,0.1 0.5,0.5 0.5,0.5 0.1,0.1 0.1))", obj_region);
diags("Input", obj_region); // check validity/attempt correction
bg_polygon_type out;
bg::intersection(crop_box(obj_region, 0.5), obj_region, out);
std::cout << "Output: " << bg::wkt(out) << "\n";
Конечно, определите crop_box
аналогично тому, как вы это делали (только, получите конверт обрезки, а не многоточечный):
template <typename G>
bg_box_type crop_box(G const& geom, double ratio) {
bg_box_type env;
bg::envelope(geom, env);
auto miny = env.min_corner().y();
auto height = env.max_corner().y() - miny;
env.max_corner().set<1>(miny + height*ratio);
return env;
}
Демо
Добавление визуализации тоже:
Для визуализации мы масштабируем весь ввод в 100x
Live On Coliru
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/multi_polygon.hpp>
#include <boost/geometry/geometries/linestring.hpp>
#include <boost/geometry/geometries/multi_polygon.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/arithmetic/arithmetic.hpp>
#include <boost/geometry/algorithms/for_each.hpp>
#include <boost/geometry/algorithms/envelope.hpp>
#include <boost/geometry/algorithms/intersection.hpp>
#include <boost/geometry/io/io.hpp>
#include <iostream>
#include <fstream>
namespace bg = boost::geometry;
using bg_point_type = bg::model::d2::point_xy<double>;
using bg_ring_type = bg::model::ring<bg_point_type>;
using bg_polygon_type = bg::model::polygon<bg_point_type>;
using bg_multipolygon_type = bg::model::multi_polygon<bg_polygon_type>;
using bg_box_type = bg::model::box<bg_point_type>;
template <typename G>
bg_box_type crop_box(G const& geom, double ratio) {
bg_box_type env;
bg::envelope(geom, env);
auto const maxy = env.max_corner().y();
auto const height = maxy - env.min_corner().y() ;
env.min_corner().set<1>(maxy - height*ratio);
return env;
}
template <typename G>
void diags(std::string name, G& geom) {
std::cout << name << ": " << bg::wkt(geom) << "\n";
std::string reason;
if (!bg::is_valid(geom, reason)) {
std::cout << name << ": " << reason << "\n";
bg::correct(geom);
std::cout << bg::wkt(geom) << "\n";
if (!bg::is_valid(geom, reason)) {
std::cout << name << " corrected: " << reason << "\n";
}
}
}
int main(int argc, char** argv) {
bool const visualize = argc>1 && argv[1]==std::string("-v");
//Find the envelope of the region of the object
bg_polygon_type obj_region;
bg::read_wkt("POLYGON((0.1 0.1,0.1 0.5,0.5 0.5,0.5 0.1,0.1 0.1))", obj_region);
if (visualize) {
bg::for_each_point(obj_region, [](bg_point_type& p) { bg::multiply_value(p, 100.0); });
}
diags("Input", obj_region);
bg_multipolygon_type out;
bg::intersection(crop_box(obj_region, 0.5), obj_region, out);
diags("Output", out);
if (visualize) {
std::ofstream svg("svg.svg");
boost::geometry::svg_mapper<bg_point_type> mapper(svg, 600, 600);
mapper.add(obj_region);
mapper.add(out);
mapper.map(obj_region, "fill-opacity:0.5;fill:rgb(204,153,0);stroke:rgb(204,153,0);stroke-width:2");
mapper.map(out, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
}
}
Отпечатки
Input: POLYGON((0.1 0.1,0.1 0.5,0.5 0.5,0.5 0.1,0.1 0.1))
Output: MULTIPOLYGON(((0.1 0.5,0.5 0.5,0.5 0.3,0.1 0.3,0.1 0.5)))
Или когда Включение визуализации :
Input: POLYGON((10 10,10 50,50 50,50 10,10 10))
Output: MULTIPOLYGON(((10 50,50 50,50 30,10 30,10 50)))
И следующий вывод SVG: