I have a base multipolygon B and a number of other multipolygons p_0, p_1, p_2, ... each of which has a weight w_0, w_1, w_2, ....
I would like to:
B and p_0 and multiply this area by w_0.B and p_0 to form a new multipolygon B_1 which is B with p_0 cut out of it.B_1 and p_1 and multiply this area by w_1.B_1 and p_1 to form a new multipolygon B_2 which is B_1 with p_1 cut out of it.I am trying to do this using the boost::polygon library. They have an example of how to do intersections here and differences here.
The intersection function is defined as follows:
bool intersection(Geometry1 const & geometry1, Geometry2 const & geometry2, GeometryOut & geometry_out)
In order to use intersection in a loop, and measure areas, I think I need to convert the GeometryOut type to the Geometry1 type. It is not clear how to do this.
A reduced version of my code thus far is:
#include <boost/polygon/polygon.hpp>
#include <boost/geometry/geometry.hpp>
#include <boost/geometry/multi/geometries/multi_polygon.hpp>
#include <iostream>
#include <vector>
#include <cstdlib>
namespace gtl = boost::polygon;
typedef gtl::polygon_data<float> Polygon;
Polygon MakeBox(float xmin, float ymin, float xmax, float ymax){
typedef gtl::polygon_traits<Polygon>::point_type Point;
Point pts[] = {
gtl::construct<Point>(xmin, ymin),
gtl::construct<Point>(xmin, ymax),
gtl::construct<Point>(xmax, ymax),
gtl::construct<Point>(xmax, ymin)
};
Polygon poly;
gtl::set_points(poly, pts, pts+4);
return poly;
}
double GetValue(const Polygon &base, const std::vector<Polygon> &polys, const std::vector<double> &vals){
double value = 0;
double base_area = gtl::area(base);
boost::geometry::model::multi_polygon<Polygon> multipoly, isect, multipolynew;
//I've also tried the following line in place of the line above: doesn't work.
//std::deque<Polygon> multipoly, isect, multipolynew;
multipoly.push_back(base);
for(int i=0;i<polys.size();i++){
boost::geometry::intersection(multipoly, polys[i], isect);
value += gtl::area(isect)/base_area*vals[i];
boost::geometry::difference(multipoly,polys[i],multipolynew);
multipoly = multipolynew;
}
return value;
}
int main() {
Polygon base = MakeBox(0,0,10,10); //Base polygon
std::vector<Polygon> polys;
std::vector<double> vals;
//Set up some random input
for(int i=0;i<3;i++){
int xmin=rand()%10;
int xmax=xmin+rand()%10;
int ymin=rand()%10;
int ymax=ymin+rand()%10;
polys.push_back(MakeBox(xmin,ymin,xmax,ymax));
vals.push_back(rand()%100);
}
std::cout<<GetValue(base,polys,vals)<<std::endl;
}
Looking at this from all angles I would conclude that this is a limitation of Boost Geometry adaptation of Boost Polygon types.
Even when replacing the "vanilla" ring type (boost::polygon_data<>) by the more full featured polygon concept implementation from Boost Polygon (boost::polygon_with_holes_data<>) doesn't allow the adapted Polygon to be used in a boost::geometry::multi_poygon<T> instantiation.
In case you are interested I have some more information uncovered during my searches:
You can do what you wish when using Boost Geometry models instead of Boost Polygon types. There is no problem using the multi_polygon as a Geometry directly there:
#include <fstream>
#include <boost/geometry/geometry.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/geometries/geometries.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/foreach.hpp>
typedef boost::geometry::model::d2::point_xy<double> PointType;
using Polygon = boost::geometry::model::polygon<PointType>;
Polygon MakeBox(float xmin, float ymin, float xmax, float ymax){
Polygon poly;
poly.outer().assign({
{xmin, ymin},
{xmin, ymax},
{xmax, ymax},
{xmax, ymin}
});
return poly;
}
double GetValue(const Polygon &base, const std::vector<Polygon> &polys, const std::vector<double> &vals){
double value = 0;
double base_area = boost::geometry::area(base);
boost::geometry::model::multi_polygon<Polygon> multipoly, isect, multipolynew;
multipoly.push_back(base);
for(size_t i=0;i<polys.size();i++){
boost::geometry::intersection(multipoly, polys[i], isect);
value += boost::geometry::area(isect)/base_area*vals[i];
boost::geometry::difference(multipoly,polys[i],multipolynew);
multipoly = multipolynew;
}
return value;
}
int main()
{
Polygon base = MakeBox(0,0,10,10); //Base polygon
std::vector<Polygon> polys;
std::vector<double> vals;
//Set up some random input
for(int i=0;i<3;i++){
int xmin=rand()%10;
int xmax=xmin+rand()%10;
int ymin=rand()%10;
int ymax=ymin+rand()%10;
polys.push_back(MakeBox(xmin,ymin,xmax,ymax));
vals.push_back(rand()%100);
}
std::cout<<GetValue(base,polys,vals)<<std::endl;
}
Alternatively you could implement a O(n log n) visitation version of intersect that intersects two collections (e.g. vectors) of Polygons even though they're not a single Geometry.
This mailing list exchange details precisely this http://boost-geometry.203548.n3.nabble.com/intersection-of-two-vectors-of-polygons-tp2875513p2880997.html (even though the use case here is so that the individual polygon can be associated with unique bits of information (ids) which is hard to do with multi_polygons).
Here's a simplified adaptation of that code that does work for intersection. I suppose you'd need to do the similar approach for the difference call (?). I'll leave that as an exercise for the reader.
Another relevant thread might be http://boost-geometry.203548.n3.nabble.com/Unioning-many-polygons-td3405482.html