Skip to content

Commit

Permalink
Merge branch 'gsoc2023-aos_sphere_demo-denizdiktas' of github.com:CGA…
Browse files Browse the repository at this point in the history
…L/cgal-public-dev into gsoc2023-aos_sphere_demo-denizdiktas
  • Loading branch information
efifogel committed Jun 1, 2023
2 parents e55287a + 34b9a7d commit d565320
Showing 1 changed file with 82 additions and 66 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -2916,8 +2916,8 @@ class Arr_geodesic_arc_on_sphere_traits_2 : public Kernel_ {
typedef double Approximate_number_type;
typedef CGAL::Cartesian<Approximate_number_type> Approximate_kernel;
typedef Arr_extended_direction_3<Approximate_kernel> Approximate_point_2;
typedef CGAL::Vector_3<Approximate_kernel> Approximate_Vector_3;
typedef Approximate_kernel::Direction_3 Approximate_Direction_3;
typedef CGAL::Vector_3<Approximate_kernel> Approximate_vector_3;
typedef Approximate_kernel::Direction_3 Approximate_direction_3;

class Approximate_2 {
public:
Expand Down Expand Up @@ -2946,73 +2946,89 @@ class Arr_geodesic_arc_on_sphere_traits_2 : public Kernel_ {
/*! Obtain an approximation of an \f$x\f$-monotone curve.
*/
template <typename OutputIterator>
OutputIterator operator()(const X_monotone_curve_2& xcv, Approximate_number_type error,
OutputIterator operator()(const X_monotone_curve_2& xcv,
Approximate_number_type error,
OutputIterator oi, bool l2r = true) const {
auto s = xcv.source();
auto t = xcv.target();
auto n = xcv.normal();

// get the approximate points;
auto as = (*this)(s);
auto at = (*this)(t);

// convert the approximate points to vectors with approximate-kernel
auto normalize = [](auto& x) { x /= sqrt(x.squared_length()); };
auto getVector_3 = [](const Approximate_point_2& ap) { return Approximate_Vector_3(ap.dx(), ap.dy(), ap.dz()); };
auto vs = getVector_3(as);
auto vt = getVector_3(at);
auto vn = Approximate_Vector_3(CGAL::to_double(n.dx()), CGAL::to_double(n.dy()), CGAL::to_double(n.dz()));
normalize(vs);
normalize(vt);
normalize(vn);

// direction vectors for the coordinate system where we are going to make the approximation:
auto axisX = vs; // x-axis will coincide with the vector from the origin to the normalized SOURCE-vector
auto axisZ = vn; // this will make sure that the orientation of the curve during approximation is consistent with the curve
auto axisY = CGAL::cross_product(axisZ, axisX);
normalize(axisY);

// In this coordinate system the source has local coords (0,0), hence its initial angle with the X-axis is 0
// now compute the local coordinates and the angle it makes with the X-axis
Approximate_number_type theta;
if (xcv.is_full())
{
theta = 2.0 * CGAL_PI;
}
else
{
auto ltx = CGAL::scalar_product(axisX, vt);
auto lty = CGAL::scalar_product(axisY, vt);
theta = std::atan2(lty, ltx);
if (theta < 0)
theta += 2.0 * CGAL_PI;
}

// compute the number of divisions given the error
const Approximate_number_type R = 1.0; // radius is always 1
Approximate_number_type dtheta = 2.0 * acos(1 - error / R);
int N = ceil(theta / dtheta);
dtheta = theta / N;

// generate the points
vector<Approximate_Vector_3> pts;
pts.reserve(N + 1);
pts.push_back(vs);
for (int i = 1; i < N; i++)
{
const Approximate_number_type angle = i * dtheta;
auto p = cos(angle) * axisX + sin(angle) * axisY;
pts.push_back(p);
}
pts.push_back(vt);
auto s = xcv.source();
auto t = xcv.target();
auto n = xcv.normal();

// get the approximate points;
auto as = (*this)(s);
auto at = (*this)(t);

// convert the approximate points to vectors with approximate-kernel
auto vs = get_approximate_vector_3(as);
auto vt = get_approximate_vector_3(at);
auto vn = get_approximate_vector_3(n);

// normalize the vectors
auto normalize = [](auto& x) { x /= std::sqrt(x.squared_length()); };
normalize(vs);
normalize(vt);
normalize(vn);

// Define the spanning vectors of the coordinate system where we are
// going to make the approximation:
auto axisX = vs; // x-axis will coincide with the vector from the
// origin to the normalized SOURCE-vector
auto axisZ = vn; // this will make sure that the orientation of the
// approximated curve is consistent with the curve
auto axisY = CGAL::cross_product(axisZ, axisX);
normalize(axisY);

// In this coordinate system the source has local coords (0,0), hence its
// initial angle with the X-axis is 0 degrees (radians)
// Compute the local coordinates and the angle it makes with the X-axis
Approximate_number_type theta;
if (xcv.is_full())
{
theta = 2.0 * CGAL_PI;
}
else
{
auto ltx = CGAL::scalar_product(axisX, vt);
auto lty = CGAL::scalar_product(axisY, vt);
theta = std::atan2(lty, ltx);
if (theta < 0)
theta += 2.0 * CGAL_PI;
}

// compute the number of divisions given the requested error
const Approximate_number_type R = 1.0; // radius is always 1
Approximate_number_type dtheta = 2.0 * std::acos(1 - error / R);
int N = std::ceil(theta / dtheta);
dtheta = theta / N;

// generate the points approximating the curve
const auto loc = Approximate_point_2::NO_BOUNDARY_LOC;
*oi++ = get_approximate_point_2(vs, loc); // source vector
for (int i = 1; i < N; ++i)
{
const Approximate_number_type angle = i * dtheta;
auto p = std::cos(angle) * axisX + std::sin(angle) * axisY;
*oi++ = get_approximate_point_2(p, loc);
}
*oi++ = get_approximate_point_2(vt, loc); // target vector

for (const auto& p : pts)
{
Approximate_point_2 ap(Approximate_Direction_3(p.x(), p.y(), p.z()), Approximate_point_2::NO_BOUNDARY_LOC);
*oi++ = ap;
}
return oi;
}

return oi;
private:
Approximate_vector_3 get_approximate_vector_3(const Approximate_point_2& p)
const {
return Approximate_vector_3(p.dx(), p.dy(), p.dz());
};

Approximate_vector_3 get_approximate_vector_3(const Direction_3& d) const {
return Approximate_vector_3(CGAL::to_double(d.dx()),
CGAL::to_double(d.dy()), CGAL::to_double(d.dz()));
};

Approximate_point_2 get_approximate_point_2(const Approximate_vector_3& v,
const Approximate_point_2::Location_type loc) const {
Approximate_direction_3 d(v.x(), v.y(), v.z());
return Approximate_point_2(d, loc);
}
};

Expand Down

0 comments on commit d565320

Please sign in to comment.