c++cgal
How does one obtain the area of a general polygon set in CGAL?
In the following CGAL example, drawn from here, I have a union of circles and rectangles. It looks like so:
How do I get the area of the resulting union using CGAL?
// Compile with: clang++ -DBOOST_ALL_NO_LIB -DCGAL_USE_GMPXX=1 -O2 -g -DNDEBUG -Wall -Wextra -pedantic -march=native -frounding-math bob.cpp -lgmpxx -lmpfr -lgmp
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Gps_circle_segment_traits_2.h>
#include <CGAL/General_polygon_set_2.h>
#include <CGAL/Lazy_exact_nt.h>
#include <list>
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_2 Point_2;
typedef Kernel::Circle_2 Circle_2;
typedef CGAL::Gps_circle_segment_traits_2<Kernel> Traits_2;
typedef CGAL::General_polygon_set_2<Traits_2> Polygon_set_2;
typedef Traits_2::General_polygon_2 Polygon_2;
typedef Traits_2::General_polygon_with_holes_2 Polygon_with_holes_2;
typedef Traits_2::Curve_2 Curve_2;
typedef Traits_2::X_monotone_curve_2 X_monotone_curve_2;
// Construct a polygon from a circle.
Polygon_2 construct_polygon (const Circle_2& circle)
{
// Subdivide the circle into two x-monotone arcs.
Traits_2 traits;
Curve_2 curve (circle);
std::list<CGAL::Object> objects;
traits.make_x_monotone_2_object() (curve, std::back_inserter(objects));
CGAL_assertion (objects.size() == 2);
// Construct the polygon.
Polygon_2 pgn;
X_monotone_curve_2 arc;
std::list<CGAL::Object>::iterator iter;
for (iter = objects.begin(); iter != objects.end(); ++iter) {
CGAL::assign (arc, *iter);
pgn.push_back (arc);
}
return pgn;
}
// Construct a polygon from a rectangle.
Polygon_2 construct_polygon (const Point_2& p1, const Point_2& p2,
const Point_2& p3, const Point_2& p4)
{
Polygon_2 pgn;
X_monotone_curve_2 s1(p1, p2); pgn.push_back(s1);
X_monotone_curve_2 s2(p2, p3); pgn.push_back(s2);
X_monotone_curve_2 s3(p3, p4); pgn.push_back(s3);
X_monotone_curve_2 s4(p4, p1); pgn.push_back(s4);
return pgn;
}
// The main program:
int main ()
{
// Insert four non-intersecting circles.
Polygon_set_2 S;
Polygon_2 circ1, circ2, circ3, circ4;
circ1 = construct_polygon(Circle_2(Point_2(1, 1), 1)); S.insert(circ1);
circ2 = construct_polygon(Circle_2(Point_2(5, 1), 1)); S.insert(circ2);
circ3 = construct_polygon(Circle_2(Point_2(5, 5), 1)); S.insert(circ3);
circ4 = construct_polygon(Circle_2(Point_2(1, 5), 1)); S.insert(circ4);
// Compute the union with four rectangles incrementally.
Polygon_2 rect1, rect2, rect3, rect4;
rect1 = construct_polygon(Point_2(1, 0), Point_2(5, 0),
Point_2(5, 2), Point_2(1, 2));
S.join (rect1);
rect2 = construct_polygon(Point_2(1, 4), Point_2(5, 4),
Point_2(5, 6), Point_2(1, 6));
S.join (rect2);
rect3 = construct_polygon(Point_2(0, 1), Point_2(2, 1),
Point_2(2, 5), Point_2(0, 5));
S.join (rect3);
rect4 = construct_polygon(Point_2(4, 1), Point_2(6, 1),
Point_2(6, 5), Point_2(4, 5));
S.join (rect4);
// Print the output.
std::list<Polygon_with_holes_2> res;
S.polygons_with_holes (std::back_inserter (res));
std::copy (res.begin(), res.end(),
std::ostream_iterator<Polygon_with_holes_2>(std::cout, "\n"));
std::cout << std::endl;
return 0;
}
Solution
General polygons with possibly circular edges don't have the area function - so we need to get into details of their representation and write this kind of functions ourselves (sigh!!). An example of such an area function is below:
// ------ return signed area under the linear segment (P1, P2)
auto area(Traits_2::Point_2 const& P1, Traits_2::Point_2 const& P2)
{
auto const dx = CGAL::to_double(P1.x()) - CGAL::to_double(P2.x());
auto const sy = CGAL::to_double(P1.y()) + CGAL::to_double(P2.y());
return dx * sy / 2;
}
// ------ return signed area under the circular segment (P1, P2, C)
auto area(Traits_2::Point_2 const& P1, Traits_2::Point_2 const& P2, Circle const& C)
{
auto const dx = CGAL::to_double(P1.x()) - CGAL::to_double(P2.x());
auto const dy = CGAL::to_double(P1.y()) - CGAL::to_double(P2.y());
auto const squaredChord = dx * dx + dy * dy;
auto const chord = std::sqrt(squaredChord);
auto const squaredRadius = CGAL::to_double(C.squared_radius());
auto const areaSector = squaredRadius * std::asin(std::min(1.0, chord / (std::sqrt(squaredRadius) * 2)));
auto const areaTriangle = chord * std::sqrt(std::max(0.0, squaredRadius * 4 - squaredChord)) / 4;
auto const areaCircularSegment = areaSector - areaTriangle;
return area(P1, P2) + C.orientation() * areaCircularSegment;
}
// ------ return signed area under the X-monotone curve
auto area(X_monotone_curve_2 const& XCV)
{
if (XCV.is_linear())
{
return area(XCV.source(), XCV.target());
}
else if (XCV.is_circular())
{
return area(XCV.source(), XCV.target(), XCV.supporting_circle());
}
else
{
return 0.0;
}
}
// ------ return area of the simple polygon
auto area(Polygon_2 const& P)
{
auto res = 0.0;
for (auto it = P.curves_begin(); it != P.curves_end(); ++it) res += area(*it);
return res;
}
// ------ return area of the polygon with (optional) holes
auto area(Polygon_with_holes_2 const& P)
{
auto res = area(P.outer_boundary());
for (auto it = P.holes_begin(); it != P.holes_end(); ++it) res += area(*it);
return res;
}
As you can see I used here Polygon_with_holes_2 access functions outer_boundary, holes_begin and holes_end. For Polygon_2 I used only curves_begin and curves_end access functions. Much more efforts are required to process edges of this polygon, which are of the X_monotone_curve_2 type. All these access functions will be necessary for any other functions you'll need to develop for general polygons.
A full MWE incorporating your example is below:
// Compile with: clang++ -DBOOST_ALL_NO_LIB -DCGAL_USE_GMPXX=1 -O2 -g -DNDEBUG -Wall -Wextra -pedantic -march=native -frounding-math bob.cpp -lgmpxx -lmpfr -lgmp
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Gps_circle_segment_traits_2.h>
#include <CGAL/General_polygon_set_2.h>
#include <CGAL/Lazy_exact_nt.h>
#include <list>
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_2 Point_2;
typedef Kernel::Circle_2 Circle_2;
typedef CGAL::Gps_circle_segment_traits_2<Kernel> Traits_2;
typedef CGAL::General_polygon_set_2<Traits_2> Polygon_set_2;
typedef Traits_2::General_polygon_2 Polygon_2;
typedef Traits_2::General_polygon_with_holes_2 Polygon_with_holes_2;
typedef Traits_2::Curve_2 Curve_2;
typedef Traits_2::X_monotone_curve_2 X_monotone_curve_2;
//For two circles of radii R and r and centered at (0,0) and (d,0) intersecting
//in a region shaped like an asymmetric lens.
constexpr double lens_area(const double r, const double R, const double d){
return r*r*std::acos((d*d+r*r-R*R)/2/d/r) + R*R*std::acos((d*d+R*R-r*r)/2/d/R) - 0.5*std::sqrt((-d+r+R)*(d+r-R)*(d-r+R)*(d+r+R));
}
// ------ return signed area under the linear segment (P1, P2)
auto area(const Traits_2::Point_2& P1, const Traits_2::Point_2& P2){
auto const dx = CGAL::to_double(P1.x()) - CGAL::to_double(P2.x());
auto const sy = CGAL::to_double(P1.y()) + CGAL::to_double(P2.y());
return dx * sy / 2;
}
// ------ return signed area under the circular segment (P1, P2, C)
auto area(const Traits_2::Point_2& P1, const Traits_2::Point_2& P2, const Circle_2& C){
auto const dx = CGAL::to_double(P1.x()) - CGAL::to_double(P2.x());
auto const dy = CGAL::to_double(P1.y()) - CGAL::to_double(P2.y());
auto const squaredChord = dx * dx + dy * dy;
auto const chord = std::sqrt(squaredChord);
auto const squaredRadius = CGAL::to_double(C.squared_radius());
auto const areaSector = squaredRadius * std::asin(std::min(1.0, chord / (std::sqrt(squaredRadius) * 2)));
auto const areaTriangle = chord * std::sqrt(std::max(0.0, squaredRadius * 4 - squaredChord)) / 4;
auto const areaCircularSegment = areaSector - areaTriangle;
return area(P1, P2) + C.orientation() * areaCircularSegment;
}
// ------ return signed area under the X-monotone curve
auto area(const X_monotone_curve_2& XCV){
if (XCV.is_linear())
{
return area(XCV.source(), XCV.target());
}
else if (XCV.is_circular())
{
return area(XCV.source(), XCV.target(), XCV.supporting_circle());
}
else
{
return 0.0;
}
}
// ------ return area of the simple polygon
auto area(const Polygon_2 P){
auto res = 0.0;
for (auto it = P.curves_begin(); it != P.curves_end(); ++it) {
res += area(*it);
}
return res;
}
// ------ return area of the polygon with (optional) holes
auto area(const Polygon_with_holes_2& P){
auto res = area(P.outer_boundary());
for (auto it = P.holes_begin(); it != P.holes_end(); ++it) {
res += area(*it);
}
return res;
}
// Construct a polygon from a circle.
Polygon_2 construct_polygon (const Circle_2& circle){
// Subdivide the circle into two x-monotone arcs.
Traits_2 traits;
Curve_2 curve (circle);
std::list<CGAL::Object> objects;
traits.make_x_monotone_2_object() (curve, std::back_inserter(objects));
CGAL_assertion (objects.size() == 2);
// Construct the polygon.
Polygon_2 pgn;
X_monotone_curve_2 arc;
std::list<CGAL::Object>::iterator iter;
for (iter = objects.begin(); iter != objects.end(); ++iter) {
CGAL::assign (arc, *iter);
pgn.push_back (arc);
}
return pgn;
}
// Construct a polygon from a rectangle.
Polygon_2 construct_polygon (
const Point_2& p1,
const Point_2& p2,
const Point_2& p3,
const Point_2& p4
){
Polygon_2 pgn;
X_monotone_curve_2 s1(p1, p2); pgn.push_back(s1);
X_monotone_curve_2 s2(p2, p3); pgn.push_back(s2);
X_monotone_curve_2 s3(p3, p4); pgn.push_back(s3);
X_monotone_curve_2 s4(p4, p1); pgn.push_back(s4);
return pgn;
}
double area(const Polygon_set_2& Pset){
std::list<Polygon_with_holes_2> res;
Pset.polygons_with_holes(std::back_inserter(res));
double ret = 0;
for(const auto &x : res){
const auto temp = area(x);
if(!std::isnan(temp)){ //TODO: This seems to be necessary, but I'm not sure why
ret += area(x);
}
}
return ret;
}
void test1(){
// Insert four non-intersecting circles.
Polygon_set_2 S;
const auto circ1 = construct_polygon(Circle_2(Point_2(1, 1), 1));
const auto circ2 = construct_polygon(Circle_2(Point_2(5, 1), 1));
const auto circ3 = construct_polygon(Circle_2(Point_2(5, 5), 1));
const auto circ4 = construct_polygon(Circle_2(Point_2(1, 5), 1));
S.insert(circ1);
S.insert(circ2);
S.insert(circ3);
S.insert(circ4);
// Compute the union with four rectangles incrementally.
const auto rect1 = construct_polygon(Point_2(1, 0), Point_2(5, 0),
Point_2(5, 2), Point_2(1, 2));
const auto rect2 = construct_polygon(Point_2(1, 4), Point_2(5, 4),
Point_2(5, 6), Point_2(1, 6));
const auto rect3 = construct_polygon(Point_2(0, 1), Point_2(2, 1),
Point_2(2, 5), Point_2(0, 5));
const auto rect4 = construct_polygon(Point_2(4, 1), Point_2(6, 1),
Point_2(6, 5), Point_2(4, 5));
S.join(rect1);
S.join(rect2);
S.join(rect3);
S.join(rect4);
// Print the output.
std::list<Polygon_with_holes_2> res;
S.polygons_with_holes (std::back_inserter (res));
std::copy (res.begin(), res.end(),
std::ostream_iterator<Polygon_with_holes_2>(std::cout, "\n"));
std::cout << std::endl;
for(const auto &x: res){
std::cout << "Area = "<< area(x)<< " and should be "<<(4*8-4+M_PI)<<std::endl;
}
}
void test2(){
Polygon_set_2 S;
const auto circ1 = construct_polygon(Circle_2(Point_2(0, 0), 1));
const auto circ2 = construct_polygon(Circle_2(Point_2(0.5, 0), 1));
S.join(circ1);
S.join(circ2);
std::cout<<"Number of holes was "<<(S.number_of_polygons_with_holes()==1)<<" should be 1"<<std::endl;
std::cout<<"Area was "<<CGAL::to_double(area(S))<<" should be "<<(2*M_PI - lens_area(1, 1, 0.5))<<std::endl;
}
void test3(){
Polygon_set_2 S;
const auto circ1 = construct_polygon(Circle_2(Point_2(0, 0), 1));
const auto circ2 = construct_polygon(Circle_2(Point_2(0.5, 0), 1));
const auto circ3 = construct_polygon(Circle_2(Point_2(-1.8, 0), 1));
S.join(circ1);
S.join(circ2);
S.join(circ3);
std::cout<<"Number of holes was "<<(S.number_of_polygons_with_holes()==1)<<" should be 1"<<std::endl;
std::cout<<"Area was "<<CGAL::to_double(area(S))<<" should be "<<(3*M_PI - lens_area(1, 1, 0.5) - lens_area(1, 1, 1.8))<<std::endl;
}
// The main program:
int main (){
std::cout<<"test1"<<std::endl;
test1();
std::cout<<"test2"<<std::endl;
test2();
std::cout<<"test3"<<std::endl;
test3();
return 0;
}
Running the test1 gives an output of 31.1416. You have four rectangles of area 8 each with four overlapping regions of area 1 and 4 quarter circles giving: 4*8-4+π=31.141592653589793, so the results match the theory. Other two tests also work OK.
- How memory address for pointer to arrays is same as an element in 2D array?
- How to use ellipsis in c's case statement?
- Fast & accurate atan/arctan approximation algorithm
- How can I exclude non-numeric keys? CS50 Caesar Pset2
- Fast ceiling of an integer division in C / C++
- Is there an invalid pthread_t id?
- How does SIMD (avx) processing work? for example, if I want 10 32 bit floats how do i fit in a 256 bit avx vector?
- FDCAN problems on STM32G4
- How does the call macro enable mutual recursion between functions f and g in this Hanoi Tower implementation?
- Running test on Rocket core CPU - global variable initialized to 0 is unsuccessful, output wrong value instead
- Interacting with C arrays without knowing the size
- Combination of two strings
- Avoiding strcpy overflow destination warning
- carriage return by fgets
- How to use special characters in C?
- Why does 1.0/100.0 == 0.1/10.0 give True?
- Is it correct to compare pointers in C?
- Force free() to return malloc memory back to OS
- How can I print to standard error in C with 'printf'?
- What is the standard behavior of fread in C on Windows?
- How is strtok removing lines it shouldn't have access to?
- Using array as smart point in C
- Assigning string to malloced 2d char array not working as intended
- How to refactor repetition inside a Makefile?
- Why does an empty preprocessor command still evaluate to something?
- How to implement variable sized array within C struct
- Character array typecasting to integer
- Handling HTTP Headers in a Minimal C HTTP Server
- How to get the sign, mantissa and exponent of a floating point number
- Why do MCU libraries use logic operations instead of bitfield structs?