diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Shape_base.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Shape_base.h index 4175d064ca1..59d30096295 100644 --- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Shape_base.h +++ b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Shape_base.h @@ -163,7 +163,6 @@ namespace CGAL { std::vector map; map.reserve(64); map.resize(2); - unsigned int label = 2; for (std::size_t y = 0;y 3) + post_wrap(bitmap, u_extent, v_extent, map); // Update labels for (std::size_t y = 0;y &labels, unsigned int i, unsigned int &new_value) const { + void inline update_label(std::vector &labels, unsigned int i, + unsigned int &new_value) const { if (labels[i] != i) update_label(labels, labels[i], new_value); @@ -466,6 +467,7 @@ namespace CGAL { // Avoid compiler warnings about unused parameters. (void)indices; (void)parameter_space; + (void)cluster_epsilon; (void)min; (void)max; } @@ -522,7 +524,8 @@ namespace CGAL { arg != -std::numeric_limits::infinity(); } - void compute_bound(const std::size_t num_evaluated_points, const std::size_t num_available_points) { + void compute_bound(const std::size_t num_evaluated_points, + const std::size_t num_available_points) { hypergeometrical_dist(-2 - num_evaluated_points, -2 - num_available_points, -1 - signed(m_indices.size()), @@ -557,7 +560,8 @@ namespace CGAL { // /// \cond SKIP_IN_MANUAL /*! - Contains indices of the inliers of the candidate, access to the point and normal data is provided via property maps. + Contains indices of the inliers of the candidate, access + to the point and normal data is provided via property maps. */ std::vector m_indices; diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Sphere.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Sphere.h index e5cbf8c13ac..c3e8b44a117 100644 --- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Sphere.h +++ b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Sphere.h @@ -24,6 +24,17 @@ #include #include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif +#ifndef M_PI_2 +#define M_PI_2 1.57079632679489661923 +#endif +#ifndef M_PI_4 +#define M_PI_4 0.785398163397448309616 +#endif /*! \file Sphere.h @@ -220,12 +231,304 @@ namespace CGAL { return 3; } - virtual bool supports_connected_component() const { - return false; + // Maps to the range [-1,1]^2 + static void concentric_mapping(FT phi, FT proj, FT rad, FT &x, FT &y) { + phi = (phi < FT(-M_PI_4)) ? phi + FT(2 * M_PI) : phi; + FT r = FT(acos(double(CGAL::abs(proj)))) / FT(M_PI_2); + + FT a = 0, b = 0; + if (phi < FT(M_PI_4)) { + a = r; + b = phi * r / FT(M_PI_4); + } + else if (phi < FT(3.0 * M_PI_4)) { + a = -FT(phi - M_PI_2) * r / FT(M_PI_4); + b = r; + } + else if (phi < FT(5.0 * M_PI_4)) { + a = -r; + b = (phi - FT(M_PI)) * (-r) / FT(M_PI_4); + } + else { + a = (phi - 3 * FT(M_PI_2)) * r / FT(M_PI_4); + b = -r; + } + + x = a; + y = b; + + // Map into hemisphere + if (proj >= 0) + y += 1; + else + y = -1 - y; + + // Scale to surface distance + x = FT(x * M_PI_2 * rad); + y = FT(y * M_PI_2 * rad); } + virtual void parameters(const std::vector &indices, + std::vector > ¶meterSpace, + FT &cluster_epsilon, + FT min[2], + FT max[2]) const { + Vector_3 axis; + FT rad = radius(); + // Take average normal as axis + for (std::size_t i = 0;inormal(indices[i]); + axis = axis / (CGAL::sqrt(axis.squared_length())); + + // create basis x,y + Vector_3 d1 = Vector_3((FT) 0, (FT) 0, (FT) 1); + Vector_3 d2 = CGAL::cross_product(axis, d1); + FT l = d2.squared_length(); + if (l < (FT)0.0001) { + d1 = Vector_3((FT) 1, (FT) 0, (FT) 0); + d2 = CGAL::cross_product(axis, d1); + l = d2.squared_length(); + } + d2 = d2 / CGAL::sqrt(l); + + d1 = CGAL::cross_product(axis, d2); + l = CGAL::sqrt(d1.squared_length()); + if (l == 0) + return; + + d1 = d1 * (FT)1.0 / l; + + // Process first point separately to initialize min/max + + Vector_3 vec = this->point(indices[0]) - m_sphere.center(); + FT proj = axis * vec; // sign indicates northern or southern hemisphere + FT phi = atan2(vec * d2, vec * d1); + FT x, y; + concentric_mapping(phi, proj, rad, x, y); + + min[0] = max[0] = x; + min[1] = max[1] = y; + + parameterSpace[0] = std::pair(x, y); + + for (std::size_t i = 1;ipoint(indices[i]) - m_sphere.center(); + proj = axis * vec; // sign indicates northern or southern hemisphere + phi = atan2(vec * d2, vec * d1); + + concentric_mapping(phi, proj, rad, x, y); + + min[0] = (std::min)(min[0], x); + max[0] = (std::max)(max[0], x); + + min[1] = (std::min)(min[1], y); + max[1] = (std::max)(max[1], y); + + parameterSpace[i] = std::pair(x, y); + } + + // Is close to wrapping around? Check all three directions separately + m_wrap_right = abs(max[0] - M_PI_2 * rad) < (cluster_epsilon * 0.5); + m_wrap_left = abs(min[0] + M_PI_2 * rad) < (cluster_epsilon * 0.5); + + FT diff_top = CGAL::abs(-FT(M_PI * rad) - min[1]) + + FT(M_PI * rad) - max[1]; + m_wrap_top = diff_top < cluster_epsilon; + + if (m_wrap_top || m_wrap_left || m_wrap_right) { + cluster_epsilon = FT(M_PI * rad) + / FT(floor((M_PI * rad) / cluster_epsilon)); + + // center bitmap at equator + FT required_space = ceil( + (std::max)(CGAL::abs(min[1]), max[1]) / cluster_epsilon) + * cluster_epsilon; + min[1] = -required_space; + max[1] = required_space; + } + + m_equator = std::size_t((abs(min[1])) / cluster_epsilon - 0.5); + } + + virtual void post_wrap(const std::vector &bitmap, + const std::size_t &u_extent, + const std::size_t &v_extent, + std::vector &labels) const { + unsigned int l; + unsigned int nw, n, ne; + if (m_wrap_top && v_extent > 2) { + // Handle first index separately. + l = bitmap[0]; + if (l) { + n = bitmap[(v_extent - 1) * u_extent]; + + if (u_extent == 1) { + if (n && l != n) { + update_label(labels, (std::max)(n, l), + l = (std::min)(n, l)); + return; + } + } + + ne = bitmap[(v_extent - 1) * u_extent + 1]; + + // diagonal wrapping + if (m_wrap_left && v_extent > 2) { + nw = bitmap[v_extent * u_extent - 1]; + if (nw && nw != l) + update_label(labels, (std::max)(nw, l), + l = (std::min)(nw, l)); + } + + if (n && n != l) + update_label(labels, (std::max)(n, l), + l = (std::min)(n, l)); + else if (ne && ne != l) + update_label(labels, (std::max)(ne, l), + l = (std::min)(ne, l)); + } + + for (std::size_t i = 1;i)(nw, l), + l = (std::min)(nw, l)); + if (n && n != l) + update_label(labels, (std::max)(n, l), + l = (std::min)(n, l)); + else if (ne && ne != l) + update_label(labels, (std::max)(ne, l), + l = (std::min)(ne, l)); + } + + // Handle last index separately + l = bitmap[u_extent - 1]; + if (l) { + n = bitmap[u_extent * v_extent - 1]; + nw = bitmap[u_extent * v_extent - 2]; + + // diagonal wrapping + if (m_wrap_right && v_extent > 2) { + ne = bitmap[(v_extent - 1) * u_extent]; + if (ne && ne != l) + update_label(labels, (std::max)(ne, l), + l = (std::min)(ne, l)); + } + + if (n && n != l) + update_label(labels, (std::max)(n, l), + l = (std::min)(n, l)); + else if (nw && nw != l) + update_label(labels, (std::max)(nw, l), + l = (std::min)(nw, l)); + } + } + + // Walk upwards on the right side in the northern hemisphere + if (m_wrap_right && v_extent > 2) { + // First index + l = bitmap[(m_equator + 1) * u_extent - 1]; + unsigned int ws = bitmap[(m_equator + 3) * u_extent - 1]; + + if (l && ws && l != ws) + update_label(labels, (std::max)(ws, l), + l = (std::min)(ws, l)); + + for (std::size_t i = 1;i<(v_extent>>1) - 1;i++) { + l = bitmap[(m_equator - i + 1) * u_extent - 1]; + if (!l) + continue; + + unsigned int wn = bitmap[(m_equator + i) * u_extent - 1]; + unsigned int w = bitmap[(m_equator + i + 1) * u_extent - 1]; + ws = bitmap[(m_equator + i + 2) * u_extent - 1]; + + if (wn && wn != l) + update_label(labels, (std::max)(wn, l), + l = (std::min)(wn, l)); + if (w && w != l) + update_label(labels, (std::max)(w, l), + l = (std::min)(w, l)); + else if (ws && ws != l) + update_label(labels, (std::max)(ws, l), + l = (std::min)(ws, l)); + } + + // Last index + l = bitmap[u_extent - 1]; + if (l) { + unsigned int w = bitmap[u_extent * v_extent - 1]; + unsigned int wn = bitmap[(v_extent - 1) * u_extent - 1]; + + if (w && w != l) + update_label(labels, (std::max)(w, l), + l = (std::min)(w, l)); + else if (wn && wn != l) + update_label(labels, (std::max)(wn, l), + l = (std::min)(wn, l)); + } + } + + if (m_wrap_left && v_extent > 2) { + // First index + l = bitmap[(m_equator) * u_extent]; + unsigned int es = bitmap[(m_equator + 2) * u_extent]; + + if (l && l != es) + update_label(labels, (std::max)(es, l), + l = (std::min)(es, l)); + + for (std::size_t i = 1;i<(v_extent>>1) - 1;i++) { + l = bitmap[(m_equator - i) * u_extent]; + if (!l) + continue; + + unsigned int en = bitmap[(m_equator + i) * u_extent]; + unsigned int e = bitmap[(m_equator + i + 1) * u_extent]; + es = bitmap[(m_equator + i + 2) * u_extent]; + + if (en && en != l) + update_label(labels, (std::max)(en, l), + l = (std::min)(en, l)); + if (e && e != l) + update_label(labels, (std::max)(e, l), + l = (std::min)(e, l)); + else if (es && es != l) + update_label(labels, (std::max)(es, l), + l = (std::min)(es, l)); + } + + // Last index + l = bitmap[0]; + if (l) { + unsigned int w = bitmap[(v_extent - 1) * u_extent]; + unsigned int wn = bitmap[(v_extent - 2) * u_extent]; + + if (w && w != l) + update_label(labels, (std::max)(w, l), + l = (std::min)(w, l)); + else if (wn && wn != l) + update_label(labels, (std::max)(wn, l), + l = (std::min)(wn, l)); + } + } + } + + virtual bool supports_connected_component() const { + return true; + } + private: Sphere_3 m_sphere; + mutable bool m_wrap_right, m_wrap_top, m_wrap_left; + mutable std::size_t m_equator; /// \endcond }; }