I have an x_monotone_curve_2 circular cuve in CGAL and I would like to find interpolated points along said curve. A MWE demonstrating my current approach is below, but it contains a lot of inexact math and angle calculations. It feels as though there's probably a more elegant way to do this with CGAL.
What's a good way of finding interpolated points along an x_monotone_curve_2 circular curve?
MWE
// Compile with: clang++ -DBOOST_ALL_NO_LIB -DCGAL_USE_GMPXX=1 -O2 -g -DNDEBUG -Wall -Wextra -pedantic -march=native -frounding-math main.cpp -lgmpxx -lmpfr -lgmp
#include <CGAL/Arr_circle_segment_traits_2.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Gps_circle_segment_traits_2.h>
using K = CGAL::Exact_predicates_exact_constructions_kernel;
using APoint_2 = CGAL::Arr_circle_segment_traits_2<K>::Point_2;
using Point_2 = K::Point_2;
using Circle_2 = K::Circle_2;
using Vector_2 = K::Vector_2;
using Traits_2 = CGAL::Gps_circle_segment_traits_2<K>;
using X_monotone_curve_2 = Traits_2::X_monotone_curve_2;
// Seems to be the only way to convert points between trait and non-trait systems.
Point_2 pc(const APoint_2 &pt){
return Point_2(CGAL::to_double(pt.x()), CGAL::to_double(pt.y()));
}
class CircularArcInterpolator {
public:
CircularArcInterpolator(const X_monotone_curve_2 &curve) : curve(curve) {
if(!curve.is_circular()){
throw std::runtime_error("CircularArcInterpolator requires a circular curve!");
}
}
// t is in the range [0,1] where 0 is the source of the curve and 1 is the
// target
Point_2 operator()(const double t) const {
std::cout<<"\n\n";
// Supporting circle
const auto& circle = curve.supporting_circle();
// Starting point of the circular arc
const auto start = pc(curve.source());
// Finishing point (counter-clockwise) of the circular arc
const auto end = pc(curve.target());
// If source==target this is a complete circle
// Current strategy: find the angle between start/end point and x-axis
// Get difference between angles and walk along this
const Vector_2 x_axis(1,0);
const auto start_angle = CGAL::angle(start - circle.center(), x_axis);
const auto end_angle = CGAL::angle(end - circle.center(), x_axis);
double ang_diff = end_angle - start_angle;
// Special case if the start and end point describe a full circle
if(ang_diff == 0){
ang_diff = 2*M_PI;
}
std::cout<<"ang_diff = "<<ang_diff<<std::endl;
const auto cx = circle.center().x();
const auto cy = circle.center().y();
const auto radius = std::sqrt(CGAL::to_double(circle.squared_radius()));
std::cout<<"cx = "<<cx<<std::endl;
std::cout<<"cy = "<<cy<<std::endl;
std::cout<<"radius = "<<radius<<std::endl;
std::cout<<"t = "<<t<<std::endl;
std::cout<<"inc angle = "<<(start_angle + ang_diff*t)<<std::endl;
return {
cx + radius * std::cos(start_angle + ang_diff*t),
cy + radius * std::sin(start_angle + ang_diff*t)
};
}
private:
const X_monotone_curve_2 curve;
};
int main(){
Point_2 center(0,0);
Circle_2 supporting_circle(center, 1);
APoint_2 start_end(1, 0);
X_monotone_curve_2 curve(supporting_circle, start_end, start_end, CGAL::COUNTERCLOCKWISE);
CircularArcInterpolator cai(curve);
const auto right = cai(0.00); // Should be Point_2(1, 0)
const auto up = cai(0.25); // Should be Point_2(0, 1)
const auto left = cai(0.50); // Should be Point_2(-1, 0)
const auto down = cai(0.75); // Should be Point_2(0, -1)
const auto end = cai(1.00); // Should be Point_2(1, 0)
std::cout<<right<<std::endl;
std::cout<<up <<std::endl;
std::cout<<left <<std::endl;
std::cout<<down <<std::endl;
std::cout<<end <<std::endl;
return 0;
}