通过增强几何从多边形(环)裁剪多边形(环)的一部分
Posted
技术标签:
【中文标题】通过增强几何从多边形(环)裁剪多边形(环)的一部分【英文标题】:Crop part of a polygon(ring) from a polygon(ring) by boost geometry 【发布时间】:2018-05-25 16:28:22 【问题描述】:尝试从一个环从下到上裁剪一个环(50%),但结果不像我预期的那样工作。
我的解决办法是
using bg_point_type = boost::geometry::model::d2::point_xy<double>;
using bg_polygon_type = boost::geometry::model::ring<bg_point_type>;
//Find the envelope of the region of the object
bg_polygon_type obj_region = //... it is a rectangle
boost::geometry::model::box<bg_point_type> envelope_box;
boost::geometry::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
boost::geometry::model::linestring<bg_point_type> const lineleft, line_y, right, line_y;
std::vector<bg_point_type> intersect;
boost::geometry::intersection(obj_region, line, intersect);
//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));
boost::geometry::correct(result);
这个解决方案给了我需要的点,但它没有给我想要的形状。
点为“0.1,0.1,0.1,0.5,0.5,0.5,0.5,0.1,0.1,0.1”,我有一个这样的环(为简单起见,我选择矩形为例)
在我裁剪戒指之后,我希望它可以(有点,0.1,0.5,0.1,0.3,0.5,0.3,0.5,0.5,0.1,0.5)
但它给了我 0.1,0.5,0.5,0.5,0.1,0.3,0.5,0.3,0.1,0.5
问题是,“正确”的多边形不是唯一的,我怎么能得到我想要的结果,如果在这里说太长/太复杂,我应该搜索什么样的关键字?我试过“裁剪多边形”、“裁剪多边形的一部分”、“从多边形裁剪多边形”等。谢谢
编辑 1:obj_region 是简单的多边形
编辑2:sehe的解决问题
我将crop_box 函数更改为从下到上裁剪并期望 交叉点 api 会给我相同的结果,但它会返回一个奇怪的结果。
template <typename G>
bg_box_type crop_box(G const& geom, double ratio)
bg_box_type env;
bg::envelope(geom, env);
auto const miny = env.min_corner().y();
auto const height = env.max_corner().y() - miny;
env.min_corner().set<1>(env.max_corner().y() - height*ratio);
return env;
输入:多边形((100.0 100.0,100.0 200.0,200.0 200.0,200.0 100.0,100.0 100.0)) 输出:POLYGON((100 200,100 150,200 200,200 150,100 200))--无效 正确后的输出:POLYGON((100 200,100 150,200 200,200 150,100 200))--无效
完整的源代码
#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_polygon_type = bg::model::ring<bg_point_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 miny = env.min_corner().y();
auto const height = env.max_corner().y() - miny;
env.min_corner().set<1>(env.max_corner().y() - 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()
bool const visualize = true;
//Find the envelope of the region of the object
bg_polygon_type obj_region;
bg::read_wkt("POLYGON((100.0 100.0,100.0 200.0,200.0 200.0,200.0 100.0,100.0 100.0))", obj_region);
diags("Input", obj_region);
bg_polygon_type out;
bg::intersection(crop_box(obj_region, 0.5), obj_region, out);
diags("output", out);
std::cout << "Output: " << bg::wkt(out) << "\n";
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");
【问题讨论】:
【参考方案1】:您只想将原始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 lineleft, 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;
演示
也添加可视化:
为了可视化,我们将所有输入缩放 100 倍
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)))
或者当enabling visualization:
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 输出:
【讨论】:
我在使用您的解决方案时发现了一个问题,如果我更改了 min_coordinate 它将不起作用,如果您有时间,请查看我的帖子,我会编辑并添加新信息,谢谢。 是的,哎呀,我的错。在这种情况下,GeometryOut
的 docs state 应该是多面体。修正了我的答案,现在输出i.imgur.com/KonUEkj.png
谢谢,它有效。你能指出我的文档的链接吗?我想知道为什么它应该是 multi_polygon。
谢谢,没注意到。我以前读过它,但没有完全理解文档,根据你的代码和文档,我猜它是在说“如果你想找到交叉多边形,输出类型应该是 multi_polygon”。以上是关于通过增强几何从多边形(环)裁剪多边形(环)的一部分的主要内容,如果未能解决你的问题,请参考以下文章
有一个整数矩阵 MxN 如何将它们分组为具有增强几何形状的多边形?