The following code shows the intersection of two simple spherical polygons. The first is a ring around the z-axis defining a sphere truncated from the bottom. The other is a ring defining a half-sphere. When the latter is rotated about y by 90 degrees, you get a null result where you would expect it to cut the first polygon in half. A 270 degree rotation works as expected, as does a 89 and 91 degree rotation.
#include <array>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <cmath>
#include <iostream>
#include <vector>
namespace bg = boost::geometry;
using SphericalPoint =
bg::model::d2::point_xy<double, bg::cs::spherical_equatorial<bg::degree>>;
using SphericalPolygon =
boost::geometry::model::polygon<SphericalPoint,
false, // counter-clockwise
true, // closed
std::vector, std::vector, std::allocator,
std::allocator>;
auto rad(double deg) -> double { return M_PI * deg / 180.0; }
auto deg(double rad) -> double { return 180.0 * rad / M_PI; }
auto create_truncated_sphere_unit_vectors(double cone_half_angle_deg)
-> std::vector<std::array<double, 3>> {
auto const ring_radius = sin(rad(cone_half_angle_deg));
auto const ring_height = cos(rad(cone_half_angle_deg));
std::vector<std::array<double, 3>> result;
auto const num_ring_points = 361;
result.reserve(num_ring_points);
for (size_t i = 0; i < num_ring_points - 1; ++i) {
auto const angle_deg = (360.0 * i) / num_ring_points;
result.push_back({ring_radius * cos(rad(angle_deg)),
ring_radius * sin(rad(angle_deg)), ring_height});
}
result.push_back(result.front());
return result;
}
auto rotate_about_y(std::array<double, 3> const& vec, double rot_angle_deg)
-> std::array<double, 3> {
auto const rot_angle_rad = rad(rot_angle_deg);
auto const sin_angle = sin(rot_angle_rad);
auto const cos_angle = cos(rot_angle_rad);
return {cos_angle * vec[0] + sin_angle * vec[2], vec[1],
-sin_angle * vec[0] + cos_angle * vec[2]};
}
auto rotate_ring_about_y(std::vector<std::array<double, 3>> const& ring,
double rot_angle_deg)
-> std::vector<std::array<double, 3>> {
std::vector<std::array<double, 3>> result;
result.reserve(ring.size());
for (auto const& o : ring)
result.push_back(rotate_about_y(o, rot_angle_deg));
return result;
}
auto unit_vector_to_spherical_point(std::array<double, 3> const& unit_vector)
-> SphericalPoint {
auto calculate_theta = [](auto const& vec) {
if (vec[0] == 0.0 && vec[1] == 0.0) return 0.0;
return deg(atan2(vec[1], vec[0]));
};
return SphericalPoint(calculate_theta(unit_vector),
deg(asin(unit_vector[2])));
}
auto unit_vectors_to_spherical_polygon(
std::vector<std::array<double, 3>> const& unit_vectors)
-> SphericalPolygon {
SphericalPolygon result;
auto& outer = result.outer();
outer.resize(unit_vectors.size());
for (size_t i = 0; i < unit_vectors.size(); ++i)
outer.at(i) = unit_vector_to_spherical_point(unit_vectors[i]);
return result;
}
auto operator<<(std::ostream& o, SphericalPolygon const& sp) -> std::ostream& {
for (size_t i = 0; i < sp.outer().size(); ++i)
o << sp.outer().at(i).x() << ", " << sp.outer().at(i).y() << std::endl;
return o;
}
auto main() -> int {
auto const topper_unit_vectors = create_truncated_sphere_unit_vectors(30.0);
auto const topper_spherical_polygon =
unit_vectors_to_spherical_polygon(topper_unit_vectors);
auto const half_sphere_unit_vector =
create_truncated_sphere_unit_vectors(90.0);
auto const half_sphere_plus_90_unit_vector =
rotate_ring_about_y(half_sphere_unit_vector, 90.0);
auto const half_sphere_plus_270_unit_vector =
rotate_ring_about_y(half_sphere_unit_vector, 270.0);
auto const half_sphere_spherical_polygon =
unit_vectors_to_spherical_polygon(half_sphere_unit_vector);
auto const half_sphere_plus_90_spherical_polygon =
unit_vectors_to_spherical_polygon(half_sphere_plus_90_unit_vector);
auto const half_sphere_plus_270_spherical_polygon =
unit_vectors_to_spherical_polygon(half_sphere_plus_270_unit_vector);
// nominal
{
std::vector<SphericalPolygon> result;
bg::intersection(topper_spherical_polygon,
half_sphere_spherical_polygon, result);
std::cout << "nominal: " << result.size() << std::endl;
}
// plus 90
{
std::vector<SphericalPolygon> result;
bg::intersection(topper_spherical_polygon,
half_sphere_plus_90_spherical_polygon, result);
std::cout << "plus 90: " << result.size() << std::endl;
}
// plus 270
{
std::vector<SphericalPolygon> result;
bg::intersection(topper_spherical_polygon,
half_sphere_plus_270_spherical_polygon, result);
std::cout << "plus 270: " << result.size() << std::endl;
}
}
The following code shows the intersection of two simple spherical polygons. The first is a ring around the z-axis defining a sphere truncated from the bottom. The other is a ring defining a half-sphere. When the latter is rotated about y by 90 degrees, you get a null result where you would expect it to cut the first polygon in half. A 270 degree rotation works as expected, as does a 89 and 91 degree rotation.
I'm not sure, but the problem appears to be in add_rings.hpp where the area of the result comes out negative and is therefore rejected. Maybe this comment is relevant also.
https://godbolt.org/z/7h5Y9nxev