songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

added optimal bbox, assertions + fixed a few edge cases

This commit is contained in:
Dmitry Anisimov 2020-12-22 15:18:23 +01:00
parent 5c5c7c58a9
commit 6333b0e1d5
14 changed files with 181 additions and 506 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/CMakeLists.txt
View File

@ -19,11 +19,11 @@ if(CGAL_FOUND)
if(Boost_FOUND)
message(STATUS "Found Boost")
# find_package(Eigen3 3.1.0 REQUIRED)
# if(Eigen3_FOUND)
# message(STATUS "Found Eigen")
find_package(Eigen3 3.1.0 REQUIRED)
if(Eigen3_FOUND)
message(STATUS "Found Eigen")
# include(CGAL_Eigen_support)
include(CGAL_Eigen_support)
set(targets
# kinetic_2d_example
@ -37,13 +37,13 @@ if(CGAL_FOUND)
foreach(target ${targets})
create_single_source_cgal_program("${target}.cpp")
if(TARGET ${target})
target_link_libraries(${target} PUBLIC ${project_linked_libraries})
target_link_libraries(${target} PUBLIC ${project_linked_libraries} CGAL::Eigen_support)
target_compile_definitions(${target} PUBLIC ${project_compilation_definitions})
endif()
endforeach()
# else()
# message(ERROR "This program requires the Eigen library, and will not be compiled.")
# endif()
else()
message(ERROR "This program requires the Eigen library, and will not be compiled.")
endif()
else()
message(ERROR "This program requires the Boost library, and will not be compiled.")
endif()

12
Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/data/edge-case-test/test-flat-bbox.off Normal file
View File

@ -0,0 +1,12 @@
OFF
8 2 0
0.0 0.0 0.0
1.0 0.0 0.0
1.0 1.0 0.0
0.0 1.0 0.0
2.0 0.0 0.0
3.0 0.0 0.0
3.0 1.0 0.0
2.0 1.0 0.0
4 0 1 2 3
4 4 5 6 7

5
Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes_example.cpp
View File

@ -70,14 +70,15 @@ int main(const int argc, const char** argv) {
const unsigned int k = (argc > 2 ? std::atoi(argv[2]) : 1);
std::cout << "* number of intersections k: " << k << std::endl;
const unsigned int n = 0; // number of subdivisions per bbox side
const unsigned int n = 0;
std::cout << "* number of subdivisions per bbox side: " << n << std::endl;
const unsigned int num_blocks = std::pow(n + 1, 3);
std::cout << "* number of blocks: " << num_blocks << std::endl;
const double enlarge_bbox_ratio = 1.1;
std::cout << "* enlarge bbox ratio: " << enlarge_bbox_ratio << std::endl;
const bool reorient = true;
const bool reorient = false;
std::cout << "* reorient: " << (reorient ? "true" : "false") << std::endl;
// Algorithm.

12
Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_random_shapes_example.cpp
View File

@ -380,11 +380,11 @@ int main(const int argc, const char** argv) {
std::cout << std::endl;
std::cout << "--- INPUT STATS: " << std::endl;
std::cout << "* input kernel: " << boost::typeindex::type_id<EPICK>().pretty_name() << std::endl;
std::cout << "* polygon kernel: " << boost::typeindex::type_id<Kernel>().pretty_name() << std::endl;
std::cout << "* expected number of polygons: " << n << std::endl;
std::cout << "* generated number of polygons: " << rnd_polygons.size() << std::endl;
std::cout << "* number of vertices in a polygon: " << p << std::endl;
std::cout << "* input kernel: " << boost::typeindex::type_id<EPICK>().pretty_name() << std::endl;
std::cout << "* polygon kernel: " << boost::typeindex::type_id<Kernel>().pretty_name() << std::endl;
std::cout << "* expected number of polygons: " << n << std::endl;
std::cout << "* generated number of polygons: " << rnd_polygons.size() << std::endl;
std::cout << "* number of vertices in a polygon: " << p << std::endl;
// exit(EXIT_SUCCESS);
IPolygon_3 input_polygon;
@ -403,7 +403,7 @@ int main(const int argc, const char** argv) {
assert(input_polygons.size() == rnd_polygons.size());
// Algorithm.
KSR ksr(true, true);
KSR ksr(false, false);
const IPolygon_3_map polygon_map;
const unsigned int k = (argc > 3 ? std::atoi(argv[3]) : 1);
std::cout << "* input k: " << k << std::endl;

11
Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h
View File

@ -2422,11 +2422,11 @@ public:
pvertex, prev, next, crossed[i], future_points[i], future_directions[i]);
if (is_parallel) {
if (is_intersecting_iedge(min_time, max_time, prev, crossed[i])) {
CGAL_assertion_msg(i == crossed.size() - 1, "TODO: FRONT, CAN WE HAVE OTHER I HERE?");
CGAL_assertion_msg(i == crossed.size() - 1, "TODO: OPEN, PREV, CAN WE HAVE OTHER I HERE?");
prev_iedge = crossed[i];
}
if (is_intersecting_iedge(min_time, max_time, next, crossed[i])) {
CGAL_assertion_msg(i == 0, "TODO: FRONT, CAN WE HAVE OTHER I HERE?");
CGAL_assertion_msg(i == 0, "TODO: OPEN, NEXT, CAN WE HAVE OTHER I HERE?");
next_iedge = crossed[i];
}
}
@ -2654,7 +2654,6 @@ public:
CGAL_assertion(this->k(pface_front) == 1);
CGAL_assertion(this->k(pface_back) == 1);
}
// CGAL_assertion_msg(false, "TODO: FRONT AND BACK, FINISH THIS CASE!");
} else if ((!is_occupied_edge_front && !is_occupied_edge_back)) {
@ -2662,9 +2661,7 @@ public:
CGAL_assertion(this->k(pface_front) == this->k(pface_back));
add_new_pfaces(this->k(pface_front), pvertex, crossed,
future_points, future_directions, new_pvertices);
if (m_verbose) std::cout << "- continue !front && !back" << std::endl;
// CGAL_assertion_msg(false, "TODO: !FRONT AND !BACK, FINISH THIS CASE!");
} else if (is_occupied_edge_front || is_occupied_edge_back) {
@ -2681,9 +2678,7 @@ public:
} else {
CGAL_assertion_msg(false, "ERROR: WRONG OPEN CASE!");
}
if (m_verbose) std::cout << "- continue front || back" << std::endl;
// CGAL_assertion_msg(false, "TODO: FRONT OR BACK, FINISH THIS CASE!");
} else {
CGAL_assertion_msg(false, "TODO: ADD NEW OPEN CASE! DO NOT FORGET TO UPDATE K!");
@ -3608,7 +3603,7 @@ public:
auto vec2 = Vector_3(centroid, other_centroid);
vec2 = KSR::normalize(vec2);
const FT d = FT(1) / FT(100000); // TODO: CAN WE AVOID THIS VALUE?
const FT d = KSR::tolerance<FT>(); // TODO: CAN WE AVOID THIS VALUE?
const FT dot_product = vec1 * vec2;
if (dot_product < FT(0)) {

2
Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event.h
View File

@ -28,7 +28,7 @@
namespace CGAL {
namespace KSR_3 {
// TODO: Can we avoid forward declaration?
// TODO: CAN WE AVOID FORWARD DECLARATION?
template<typename Data_structure>
class Event_queue;

2
Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event_queue.h
View File

@ -146,7 +146,7 @@ public:
}
queue_by_pvertex_idx().erase(pv.first, pv.second);
// Erase by pother. TODO: Why is pother here?
// Erase by pother.
const auto po = queue_by_pother_idx().equal_range(pvertex);
const auto po_range = CGAL::make_range(po);

8
Kinetic_shape_reconstruction/include/CGAL/KSR_3/Experimental.h
View File

@ -560,7 +560,7 @@ class Experimental {
// continue ..
// std::cout << "crossed size: " << crossed.size() << std::endl;
// std::cout << "all crossed size: " << all_crossed.size() << std::endl;
// CGAL_assertion_msg(false, "TODO: BACK CROSSED > LIMIT!");
CGAL_assertion_msg(false, "TODO: BACK CROSSED > LIMIT!");
}
}
else if (front_constrained) // Border case
@ -893,7 +893,7 @@ class Experimental {
if (is_ok) {
if (num_extra_faces < 3) {
// CGAL_assertion_msg(false, "TODO: FRONT, CROSSED < LIMIT, 1 or 2 FACES!");
CGAL_assertion_msg(false, "TODO: FRONT, CROSSED < LIMIT, 1 or 2 FACES!");
CGAL_assertion(future_points.size() == asize);
CGAL_assertion(future_directions.size() == asize);
@ -943,7 +943,7 @@ class Experimental {
connect(propagated, all_crossed[i]);
crossed.push_back(all_crossed[i]); // remove events from this one
// CGAL_assertion_msg(false, "TODO: FRONT, NULL PROPAGATED CASE!");
CGAL_assertion_msg(false, "TODO: FRONT, NULL PROPAGATED CASE!");
} else {
@ -1116,7 +1116,7 @@ class Experimental {
// continue ..
// std::cout << "crossed size: " << crossed.size() << std::endl;
// std::cout << "all crossed size: " << all_crossed.size() << std::endl;
// CGAL_assertion_msg(false, "TODO: FRONT, CROSSED > LIMIT!");
CGAL_assertion_msg(false, "TODO: FRONT, CROSSED > LIMIT!");
}
}
else // Open case

92
Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h
View File

@ -24,7 +24,8 @@
// #include <CGAL/license/Kinetic_shape_reconstruction.h>
// CGAL includes.
// #include <CGAL/optimal_bounding_box.h>
#include <CGAL/optimal_bounding_box.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
// Internal includes.
#include <CGAL/KSR/utils.h>
@ -51,10 +52,13 @@ private:
using Data_structure = KSR_3::Data_structure<Kernel>;
using Polygon_splitter = KSR_3::Polygon_splitter<Data_structure, Kernel>;
using IVertex = typename Data_structure::IVertex;
using IVertex = typename Data_structure::IVertex;
using IK = CGAL::Exact_predicates_inexact_constructions_kernel;
using IFT = typename IK::FT;
using IPoint_3 = typename IK::Point_3;
using Bbox_3 = CGAL::Bbox_3;
// using OBB_traits = CGAL::Oriented_bounding_box_traits_3<Kernel>;
using Bbox_3 = CGAL::Bbox_3;
using OBB_traits = CGAL::Oriented_bounding_box_traits_3<IK>;
public:
Initializer(
@ -155,10 +159,11 @@ private:
std::array<Point_3, 8>& bbox,
FT& time_step) const {
if (reorient)
if (reorient) {
initialize_optimal_box(input_range, polygon_map, bbox);
else
} else {
initialize_axis_aligned_box(input_range, polygon_map, bbox);
}
CGAL_assertion(bbox.size() == 8);
time_step = KSR::distance(bbox.front(), bbox.back());
@ -194,24 +199,55 @@ private:
}
// Set points.
std::vector<Point_3> points;
points.reserve(num_points);
std::vector<IPoint_3> ipoints;
ipoints.reserve(num_points);
for (const auto& item : input_range) {
const auto& polygon = get(polygon_map, item);
for (const auto& p : polygon) {
const Point_3 point(p.x(), p.y(), p.z());
points.push_back(point);
for (const auto& point : polygon) {
const IPoint_3 ipoint(
static_cast<IFT>(CGAL::to_double(point.x())),
static_cast<IFT>(CGAL::to_double(point.y())),
static_cast<IFT>(CGAL::to_double(point.z())));
ipoints.push_back(ipoint);
}
}
// Compute optimal bbox.
// The order of faces corresponds to the standard order from here:
// https://doc.cgal.org/latest/BGL/group__PkgBGLHelperFct.html#gad9df350e98780f0c213046d8a257358e
// const OBB_traits obb_traits;
// CGAL::oriented_bounding_box(
// points, bbox,
// CGAL::parameters::use_convex_hull(true).
// geom_traits(obb_traits));
const OBB_traits obb_traits;
std::array<IPoint_3, 8> ibbox;
CGAL::oriented_bounding_box(
ipoints, ibbox,
CGAL::parameters::use_convex_hull(true).
geom_traits(obb_traits));
for (std::size_t i = 0; i < 8; ++i) {
const auto& ipoint = ibbox[i];
const Point_3 point(
static_cast<FT>(ipoint.x()),
static_cast<FT>(ipoint.y()),
static_cast<FT>(ipoint.z()));
bbox[i] = point;
}
const FT sq_length = CGAL::squared_distance(bbox[0], bbox[1]);
const FT sq_width = CGAL::squared_distance(bbox[0], bbox[3]);
const FT sq_height = CGAL::squared_distance(bbox[0], bbox[5]);
CGAL_assertion(sq_length >= FT(0));
CGAL_assertion(sq_width >= FT(0));
CGAL_assertion(sq_height >= FT(0));
const FT tol = KSR::tolerance<FT>();
if (sq_length < tol || sq_width < tol || sq_height < tol) {
if (m_verbose) {
std::cout << "* warning: optimal bounding box is flat, reverting ..." << std::endl;
}
initialize_axis_aligned_box(input_range, polygon_map, bbox);
} else {
if (m_verbose) {
std::cout << "* using optimal bounding box" << std::endl;
}
}
}
template<
@ -239,16 +275,34 @@ private:
Point_3(box.xmin(), box.ymin(), box.zmax()),
Point_3(box.xmax(), box.ymin(), box.zmax()),
Point_3(box.xmax(), box.ymax(), box.zmax()) };
const FT sq_length = CGAL::squared_distance(bbox[0], bbox[1]);
const FT sq_width = CGAL::squared_distance(bbox[0], bbox[3]);
const FT sq_height = CGAL::squared_distance(bbox[0], bbox[5]);
CGAL_assertion(sq_length >= FT(0));
CGAL_assertion(sq_width >= FT(0));
CGAL_assertion(sq_height >= FT(0));
const FT tol = KSR::tolerance<FT>();
if (sq_length < tol || sq_width < tol || sq_height < tol) {
CGAL_assertion_msg(false, "TODO: HANDLE FLAT AXIS ALIGNED BBOX!");
} else {
if (m_verbose) {
std::cout << "* using axis aligned bounding box" << std::endl;
}
}
}
void enlarge_bounding_box(
const FT enlarge_bbox_ratio,
std::array<Point_3, 8>& bbox) const {
CGAL_assertion_msg(
enlarge_bbox_ratio > FT(1), "TODO: HANDLE THE CASE ENLARGE_BBOX_RATIO = FT(1)");
FT enlarge_ratio = enlarge_bbox_ratio;
if (enlarge_bbox_ratio == FT(1)) {
enlarge_ratio += KSR::tolerance<FT>();
}
const auto a = CGAL::centroid(bbox.begin(), bbox.end());
Transform_3 scale(CGAL::Scaling(), enlarge_bbox_ratio);
Transform_3 scale(CGAL::Scaling(), enlarge_ratio);
for (auto& point : bbox)
point = scale.transform(point);

437
Kinetic_shape_reconstruction/include/CGAL/KSR_3/Polygon_splitter.h
View File

@ -25,7 +25,6 @@
// CGAL includes.
#include <CGAL/convex_hull_2.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>
@ -98,11 +97,6 @@ private:
using Face_index = typename Mesh_3::Face_index;
using Uchar_map = typename Mesh_3::template Property_map<Face_index, unsigned char>;
using Regular_triangulation = CGAL::Regular_triangulation_2<Kernel>;
using Weighted_point = typename Regular_triangulation::Weighted_point;
using RVertex_handle = typename Regular_triangulation::Vertex_handle;
using REdge = typename Regular_triangulation::Edge;
enum struct Merge_type { CONVEX_HULL = 0, RECTANGLE = 1 };
Data_structure& m_data;
@ -147,9 +141,6 @@ public:
if (original_faces.size() > 1) {
CGAL_assertion_msg(false,
"ERROR: WE CANNOT HAVE MULTIPLE COPLANAR PFACES!");
// merge_coplanar_pfaces(
// support_plane_idx, original_input, original_faces);
// add_bisectors(original_faces, support_plane_idx);
}
}
@ -180,9 +171,7 @@ private:
m_data.to_3d(support_plane_idx, q) << std::endl;
}
}
add_merged_pface(support_plane_idx, merged);
// CGAL_assertion_msg(false, "TODO: MERGE COPLANAR PFACES!");
}
void collect_pface_points(
@ -234,14 +223,14 @@ private:
const KSR::size_t support_plane_idx,
const std::vector<Point_2>& merged) const {
std::vector<Weighted_point> wps;
add_weighted_bbox(support_plane_idx, wps);
CGAL_assertion(wps.size() == 4);
std::vector<Point_2> bbox;
create_bbox(support_plane_idx, bbox);
CGAL_assertion(bbox.size() == 4);
for (std::size_t i = 0; i < 4; ++i) {
const std::size_t ip = (i + 1) % 4;
const auto& pi = wps[i].point();
const auto& qi = wps[ip].point();
const auto& pi = bbox[i];
const auto& qi = bbox[ip];
const Segment_2 edge(pi, qi);
for (std::size_t j = 0; j < merged.size(); ++j) {
@ -295,30 +284,25 @@ private:
original_faces.clear();
original_input.clear();
// std::vector<Vertex_handle> vhs; // TODO: Why we cannot use vhs here?
std::vector<TRI> triangulations;
// std::vector<Vertex_handle> vhs;
std::vector<Point_2> original_face;
const auto all_pfaces = m_data.pfaces(support_plane_idx);
for (const auto pface : all_pfaces) {
const auto pvertices = m_data.pvertices_of_pface(pface);
// vhs.clear();
// TRI triangulation; // TODO: remove triangulations here!
original_face.clear();
for (const auto pvertex : pvertices) {
CGAL_assertion(vhs_map.find(pvertex) != vhs_map.end());
const auto vh = vhs_map.at(pvertex);
original_face.push_back(vh->point());
// triangulation.insert(vh->point());
// vhs.push_back(vh);
}
// triangulations.push_back(triangulation);
original_faces.push_back(original_face);
original_input.push_back(m_data.input(pface));
// TODO: Why we cannot use vhs directly here? That should be more precise!
// TODO: WHY WE CANNOT USE VHS DIRECTLY HERE? THAT SHOULD BE MORE PRECISE!
// vhs.push_back(vhs.front());
// const auto cid = m_cdt.insert_constraint(vhs.begin(), vhs.end());
original_face.push_back(original_face.front());
@ -343,21 +327,8 @@ private:
// Finally, add original labels to the cdt.
if (all_pfaces.size() > 1) {
CGAL_assertion_msg(false, "ERROR: WE CANNOT HAVE MULTIPLE COPLANAR PFACES!");
// for (auto fit = m_cdt.finite_faces_begin(); fit != m_cdt.finite_faces_end(); ++fit) {
// const auto centroid = CGAL::centroid(
// fit->vertex(0)->point(),
// fit->vertex(1)->point(),
// fit->vertex(2)->point());
// for (KSR::size_t i = 0; i < triangulations.size(); ++i) {
// const auto fh = triangulations[i].locate(centroid);
// if (fh == nullptr || triangulations[i].is_infinite(fh)) {
// continue;
// }
// fit->info().input = i;
// break;
// }
// }
CGAL_assertion_msg(false,
"ERROR: WE CANNOT HAVE MULTIPLE COPLANAR PFACES!");
}
}
@ -525,16 +496,10 @@ private:
if (original_faces.size() != 1) {
CGAL_assertion_msg(original_faces.size() <= 1,
"ERROR: WE CANNOT HAVE MULTIPLE COPLANAR PFACES!");
// dump_original_faces(support_plane_idx, original_faces, "original-faces");
// dump_current_pface(pface, "current-pface-" + m_data.str(pface));
// locate_pface_among_coplanar_pfaces(pface);
} else {
m_data.input(pface) = original_input[0];
}
}
// if (original_faces.size() > 1) {
// CGAL_assertion_msg(false, "TODO: DEBUG THIS SUPPORT PLANE!");
// }
}
void reconnect_pvertices_to_ivertices() {
@ -695,159 +660,9 @@ private:
return is_okay;
}
void locate_pface_among_coplanar_pfaces(
const PFace& pface) {
// TODO: Change fh->info().input to vector!
// std::cout << "locating " << m_data.str(pface) << std::endl;
const auto centroid = m_data.centroid_of_pface(pface);
const auto fh = m_cdt.locate(m_data.to_2d(pface.first, centroid));
CGAL_assertion(fh != nullptr && !m_cdt.is_infinite(fh));
CGAL_assertion(fh->info().input != KSR::uninitialized());
m_data.input(pface) = fh->info().input;
// std::cout << "found original input: " << m_data.input(pface) << std::endl;
}
void merge_coplanar_pfaces(
void create_bbox(
const KSR::size_t support_plane_idx,
const std::vector< std::vector<KSR::size_t> >& original_input,
const std::vector< std::vector<Point_2> >& original_faces) {
const bool is_debug = false;
CGAL_assertion(support_plane_idx >= 6);
if (is_debug) {
std::cout << std::endl << "support plane idx: " << support_plane_idx << std::endl;
std::cout << "dumping support plane ... ";
dump(true, 0, support_plane_idx, "support-plane-" +
std::to_string(support_plane_idx) + ".ply");
std::cout << "done" << std::endl;
}
if (original_faces.size() > 2) {
CGAL_assertion_msg(false, "TODO: MERGE > 2 COPLANAR PFACES!");
// This code should be very similar to the one with 2 coplanar pfaces!
}
const auto all_pfaces = m_data.pfaces(support_plane_idx);
std::vector<PFace> pfaces;
pfaces.reserve(all_pfaces.size());
for (const auto pface : all_pfaces) {
pfaces.push_back(pface);
}
CGAL_assertion(pfaces.size() >= 2);
CGAL_assertion(pfaces.size() == all_pfaces.size());
if (is_debug) std::cout << "num pfaces: " << pfaces.size() << std::endl;
std::vector< std::pair<PFace, PFace> > to_be_merged;
const auto& iedges = m_data.iedges(support_plane_idx);
for (std::size_t i = 0; i < pfaces.size(); ++i) {
const auto& pface_i = pfaces[i];
const auto centroid_i = m_data.centroid_of_pface(pface_i);
for (std::size_t j = i + 1; j < pfaces.size(); ++j) {
const auto& pface_j = pfaces[j];
const auto centroid_j = m_data.centroid_of_pface(pface_j);
const Segment_2 query_segment(
m_data.to_2d(support_plane_idx, centroid_i),
m_data.to_2d(support_plane_idx, centroid_j));
bool found_intersection = false;
for (const auto& iedge : iedges) {
const auto source = m_data.source(iedge);
const auto target = m_data.target(iedge);
const auto s = m_data.to_2d(support_plane_idx, source);
const auto t = m_data.to_2d(support_plane_idx, target);
const Segment_2 segment(s, t); Point_2 inter;
if (KSR::intersection(query_segment, segment, inter)) {
found_intersection = true;
break;
}
}
if (!found_intersection) {
to_be_merged.push_back(std::make_pair(pface_i, pface_j));
}
}
}
if (is_debug) std::cout << "pairs to be merged: " << to_be_merged.size() << std::endl;
CGAL_assertion(to_be_merged.size() == 1);
if (is_debug) {
std::cout << "merge: " <<
m_data.str(to_be_merged[0].first) << " + " <<
m_data.str(to_be_merged[0].second) << std::endl;
std::cout << "centroid i: " << m_data.centroid_of_pface(to_be_merged[0].first) << std::endl;
std::cout << "centroid j: " << m_data.centroid_of_pface(to_be_merged[0].second) << std::endl;
}
CGAL_assertion_msg(false,
"TODO: POLYGON SPLITTER, MERGE COPLANAR PFACES!");
}
void add_bisectors(
const std::vector< std::vector<Point_2> >& original_faces,
const KSR::size_t support_plane_idx) {
CGAL_assertion_msg(false,
"WARNING: WHEN ADDING BISECTORS, WE SHOULD ALSO CHANGE THE KINETIC CODE! NOT FINISHED!");
const bool is_debug = false;
CGAL_assertion(support_plane_idx >= 6);
if (is_debug) {
std::cout << std::endl << "support plane idx: " << support_plane_idx << std::endl;
std::cout << "dumping support plane ... ";
dump(true, 0, support_plane_idx, "support-plane-" +
std::to_string(support_plane_idx) + ".ply");
std::cout << "done" << std::endl;
}
std::vector<Weighted_point> wps;
create_weighted_points(original_faces, support_plane_idx, wps);
std::vector<RVertex_handle> vhs;
vhs.reserve(original_faces.size());
Regular_triangulation triangulation;
create_regular_triangulation(wps, vhs, triangulation);
CGAL_assertion(triangulation.is_valid());
CGAL_assertion(vhs.size() == original_faces.size());
if (is_debug) {
std::cout << "dumping regular triangulation / power diagram ... ";
dump(triangulation, support_plane_idx);
std::cout << "done" << std::endl;
for (std::size_t i = 0; i < vhs.size(); ++i) {
std::cout << "vhs " << std::to_string(i) << ": "
<< m_data.to_3d(support_plane_idx, vhs[i]->point().point()) << std::endl;
}
}
std::vector< std::pair<Segment_2, bool> > bisectors;
find_bisectors(is_debug, support_plane_idx, vhs, triangulation, bisectors);
CGAL_assertion(bisectors.size() > 0);
convert_bisectors_to_iedges(is_debug, support_plane_idx, wps, bisectors);
// CGAL_assertion_msg(false,
// "TODO: POLYGON SPLITTER, ADD BISECTORS!");
}
void create_weighted_points(
const std::vector< std::vector<Point_2> >& original_faces,
const KSR::size_t support_plane_idx,
std::vector<Weighted_point>& wps) const {
wps.clear();
wps.reserve(original_faces.size() + 4);
add_weighted_bbox(support_plane_idx, wps);
CGAL_assertion(wps.size() == 4);
add_weighted_input(original_faces, support_plane_idx, wps);
CGAL_assertion(wps.size() == original_faces.size() + 4);
}
void add_weighted_bbox(
const KSR::size_t support_plane_idx,
std::vector<Weighted_point>& wps) const {
std::vector<Point_2>& bbox) const {
CGAL_assertion(support_plane_idx >= 6);
const auto& iedges = m_data.iedges(support_plane_idx);
@ -863,176 +678,18 @@ private:
}
CGAL_assertion(points.size() == iedges.size() * 2);
const auto bbox = CGAL::bbox_2(points.begin(), points.end());
const Weighted_point wp1(Point_2(bbox.xmin(), bbox.ymin()), FT(0));
const Weighted_point wp2(Point_2(bbox.xmax(), bbox.ymin()), FT(0));
const Weighted_point wp3(Point_2(bbox.xmax(), bbox.ymax()), FT(0));
const Weighted_point wp4(Point_2(bbox.xmin(), bbox.ymax()), FT(0));
wps.push_back(wp1);
wps.push_back(wp2);
wps.push_back(wp3);
wps.push_back(wp4);
}
const auto box = CGAL::bbox_2(points.begin(), points.end());
const Point_2 p1(box.xmin(), box.ymin());
const Point_2 p2(box.xmax(), box.ymin());
const Point_2 p3(box.xmax(), box.ymax());
const Point_2 p4(box.xmin(), box.ymax());
void add_weighted_input(
const std::vector< std::vector<Point_2> >& original_faces,
const KSR::size_t support_plane_idx,
std::vector<Weighted_point>& wps) const {
for (const auto& original_face : original_faces) {
const auto centroid = CGAL::centroid(
original_face.begin(), original_face.end());
FT max_sq_dist = -FT(1);
for (const auto& p : original_face) {
const FT sq_dist = CGAL::squared_distance(p, centroid);
max_sq_dist = (CGAL::max)(sq_dist, max_sq_dist);
}
CGAL_assertion(max_sq_dist > FT(0));
const FT weight = static_cast<FT>(CGAL::sqrt(CGAL::to_double(max_sq_dist)));
const Weighted_point wp(centroid, weight);
wps.push_back(wp);
}
}
void create_regular_triangulation(
const std::vector<Weighted_point>& wps,
std::vector<RVertex_handle>& vhs,
Regular_triangulation& triangulation) const {
for (std::size_t i = 0; i < wps.size(); ++i) {
const auto& wp = wps[i];
if (i < 4) { triangulation.insert(wp); }
else {
const auto vh = triangulation.insert(wp);
vhs.push_back(vh);
}
}
}
void find_bisectors(
const bool is_debug,
const KSR::size_t support_plane_idx,
const std::vector<RVertex_handle>& vhs,
const Regular_triangulation& triangulation,
std::vector< std::pair<Segment_2, bool> >& bisectors) const {
if (vhs.size() > 2) {
CGAL_assertion_msg(false, "TODO: FIND MULTIPLE BISECTORS!");
// The way to find multiple bisectors is actually almost the same as for
// one bisector. We first find all pairs of unique vhs connections:
// like vhs[0] -> vhs[1], vhs[1] -> vhs[3], vhs[3] -> vhs[2], vhs[2] -> vhs[1], etc.
// then we extract the corresponding segments + possibly two rays, which
// may intersect the bbox boundaries. We transform these rays into segments.
// That is we are going to have several interior bisectors and two boundary bisectors.
}
bisectors.clear();
const auto edge = find_bisecting_edge(
is_debug, vhs[0], vhs[1], support_plane_idx, vhs, triangulation);
const auto object = triangulation.dual(edge);
Segment_2 bisector;
if (CGAL::assign(bisector, object)) {
if (is_debug) {
std::cout << "found bisector: 2 " <<
m_data.to_3d(support_plane_idx, bisector.source()) << " " <<
m_data.to_3d(support_plane_idx, bisector.target()) << std::endl;
}
} else {
CGAL_assertion_msg(false, "TODO: ADD OTHER TYPES: RAY/LINE!");
}
bisectors.push_back(std::make_pair(bisector, true));
}
const REdge find_bisecting_edge(
const bool is_debug,
const RVertex_handle& vh0,
const RVertex_handle& vh1,
const KSR::size_t support_plane_idx,
const std::vector<RVertex_handle>& vhs,
const Regular_triangulation& triangulation) const {
auto curr = triangulation.incident_edges(vh0);
const auto end = curr;
do {
const auto fh = curr->first;
const auto id = curr->second;
const auto ip = (id + 1) % 3;
const auto im = (id + 2) % 3;
const auto& p = fh->vertex(ip)->point().point();
const auto& q = fh->vertex(im)->point().point();
// std::cout << "ip, p: " << m_data.to_3d(support_plane_idx, p) << std::endl;
// std::cout << "im, q: " << m_data.to_3d(support_plane_idx, q) << std::endl;
CGAL_assertion(
fh->vertex(ip) == vh0 || fh->vertex(im) == vh0);
if (fh->vertex(ip) == vh1) {
if (is_debug) {
std::cout << "ip, found connecting edge: 2 " <<
m_data.to_3d(support_plane_idx, q) << " " <<
m_data.to_3d(support_plane_idx, p) << std::endl;
}
return *curr;
}
if (fh->vertex(im) == vh1) {
if (is_debug) {
std::cout << "im, found connecting edge: 2 " <<
m_data.to_3d(support_plane_idx, p) << " " <<
m_data.to_3d(support_plane_idx, q) << std::endl;
}
return *curr;
}
++curr;
} while (curr != end);
CGAL_assertion_msg(curr == end, "ERROR: THE BISECTING EDGE IS NOT FOUND!");
return *curr;
}
void convert_bisectors_to_iedges(
const bool is_debug,
const KSR::size_t support_plane_idx,
const std::vector<Weighted_point>& wps,
const std::vector< std::pair<Segment_2, bool> >& bisectors) {
CGAL_assertion(bisectors.size() > 0);
if (bisectors.size() > 1) {
CGAL_assertion_msg(false, "TODO: HANDLE MULTIPLE BISECTORS!");
// Should be more or less the same code as for one bisector!
}
const bool is_boundary_bisector = bisectors[0].second;
CGAL_assertion(is_boundary_bisector);
const auto& bisector = bisectors[0].first;
std::vector<Point_2> inter_points;
for (std::size_t i = 0; i < 4; ++i) {
const std::size_t ip = (i + 1) % 4;
const Segment_2 bbox_edge(wps[i].point(), wps[ip].point());
Point_2 inter_point;
if (KSR::intersection(bisector, bbox_edge, inter_point)) {
inter_points.push_back(inter_point);
}
}
CGAL_assertion(inter_points.size() == 2);
if (is_debug) {
std::cout << "found iedge: 2 " <<
m_data.to_3d(support_plane_idx, inter_points[0]) << " " <<
m_data.to_3d(support_plane_idx, inter_points[1]) << std::endl;
}
const auto iv0 = m_data.igraph().
add_vertex(m_data.to_3d(support_plane_idx, inter_points[0])).first;
const auto iv1 = m_data.igraph().
add_vertex(m_data.to_3d(support_plane_idx, inter_points[1])).first;
const auto pair = m_data.igraph().add_edge(iv0, iv1, support_plane_idx);
const auto& iedge = pair.first;
const bool is_inserted = pair.second;
if (is_inserted) {
m_data.igraph().set_line(iedge, m_data.igraph().add_line());
}
m_data.support_plane(support_plane_idx).iedges().insert(iedge);
bbox.clear();
bbox.reserve(4);
bbox.push_back(p1);
bbox.push_back(p2);
bbox.push_back(p3);
bbox.push_back(p4);
}
void dump(
@ -1119,56 +776,6 @@ private:
KSR_3::Saver<Kernel> saver;
saver.export_polygon_soup_3(polygons, file_name);
}
void dump(
const Regular_triangulation& triangulation,
const KSR::size_t support_plane_idx) {
Mesh_3 mesh;
std::map<RVertex_handle, Vertex_index> map_v2i;
for (auto vit = triangulation.finite_vertices_begin();
vit != triangulation.finite_vertices_end(); ++vit) {
const auto wp = vit->point();
const Point_2 point(wp.x(), wp.y());
map_v2i.insert(std::make_pair(
vit, mesh.add_vertex(m_data.support_plane(support_plane_idx).to_3d(point))));
}
for (auto fit = triangulation.finite_faces_begin();
fit != triangulation.finite_faces_end(); ++fit) {
std::array<Vertex_index, 3> vertices;
for (std::size_t i = 0; i < 3; ++i) {
vertices[i] = map_v2i[fit->vertex(i)];
}
mesh.add_face(vertices);
}
std::ofstream out_cdt("power-cdt.ply");
out_cdt.precision(20);
CGAL::write_ply(out_cdt, mesh);
out_cdt.close();
std::ofstream out_diagram("power-diagram.polylines.txt");
out_diagram.precision(20);
for(auto eit = triangulation.finite_edges_begin();
eit != triangulation.finite_edges_end(); ++eit) {
const auto object = triangulation.dual(eit);
// typename Kernel::Ray_2 ray;
// typename Kernel::Line_2 line;
typename Kernel::Segment_2 segment;
if (CGAL::assign(segment, object)) {
out_diagram << "2 " <<
m_data.to_3d(support_plane_idx, segment.source()) << " " <<
m_data.to_3d(support_plane_idx, segment.target()) << std::endl;
}
// if (CGAL::assign(ray, object)) out_diagram << "2 " << r << std::endl;
// if (CGAL::assign(line, object)) out_diagram << "2 " << l << std::endl;
}
out_diagram.close();
}
};
} // namespace KSR_3

2
Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h
View File

@ -122,7 +122,7 @@ public:
const FT z = normal.z() + (pa.x() - pb.x()) * (pa.y() + pb.y());
normal = Vector_3(x, y, z);
}
CGAL_assertion_msg(normal != CGAL::NULL_VECTOR, "ERROR: POLYGON HAS FLAT BBOX!");
CGAL_assertion_msg(normal != CGAL::NULL_VECTOR, "ERROR: BBOX IS FLAT!");
m_data->k = 0;
m_data->plane = Plane_3(points[0], KSR::normalize(normal));

48
Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h
View File

@ -82,7 +82,7 @@ private:
public:
Kinetic_shape_reconstruction_3(
const bool verbose = true,
const bool debug = false) :
const bool debug = false) :
m_debug(debug),
m_verbose(verbose),
m_queue(m_debug),
@ -92,51 +92,50 @@ public:
m_data(m_debug)
{ }
// TODO: Use named parameters here!
// TODO: USE NAMED PARAMETERS!
template<
typename InputRange,
typename PolygonMap>
const bool partition(
const InputRange& input_range,
const PolygonMap polygon_map,
const unsigned int k = 1,
const unsigned int n = 0,
const double enlarge_bbox_ratio = 1.1,
const bool reorient = false) {
unsigned int k = 1,
unsigned int n = 0,
double enlarge_bbox_ratio = 1.1,
bool reorient = false) {
std::cout.precision(20);
if (input_range.size() == 0) {
CGAL_warning_msg(input_range.size() != 0,
"WARNING: YOUR INPUT IS EMPTY. RETURN WITH NO CHANGE!");
return false;
}
if (k == 0) {
CGAL_warning_msg(k != 0,
"WARNING: YOU SET K TO 0. THE VALID VALUES ARE {1,2,...}. RETURN WITH NO CHANGE!");
CGAL_warning_msg(input_range.size() > 0,
"WARNING: YOUR INPUT IS EMPTY! RETURN WITH NO CHANGE!");
return false;
}
if (n != 0) {
CGAL_assertion_msg(n == 0,
"TODO: IMPLEMENT KINETIC SUBDIVISION!");
CGAL_assertion_msg(false, "TODO: IMPLEMENT KINETIC SUBDIVISION!");
if (n > 3) {
CGAL_warning_msg(n <= 3,
"WARNING: DOES IT MAKE SENSE TO HAVE MORE THAN 64 INPUT BLOCKS? RETURN WITH NO CHANGE!");
return false;
"WARNING: DOES IT MAKE SENSE TO HAVE MORE THAN 64 INPUT BLOCKS? SETTING N TO 3!");
n = 3;
}
}
if (enlarge_bbox_ratio < 1.0) {
CGAL_warning_msg(enlarge_bbox_ratio >= 1.0,
"WARNING: YOU SET ENLARGE_BBOX_RATIO < 1.0. THE VALID RANGE IS [1.0, +INF). RETURN WITH NO CHANGE!");
return false;
"WARNING: YOU SET ENLARGE_BBOX_RATIO < 1.0! THE VALID RANGE IS [1.0, +INF). SETTING TO 1.0!");
enlarge_bbox_ratio = 1.0;
}
const FT time_step = static_cast<FT>(m_initializer.initialize(
input_range, polygon_map, k, enlarge_bbox_ratio, reorient));
m_initializer.convert(m_data);
m_data.check_integrity();
if (k == 0) {
CGAL_warning_msg(k > 0,
"WARNING: YOU SET K TO 0! THAT MEANS NO PROPAGATION! THE VALID VALUES ARE {1,2,...}. INTERSECT AND RETURN!");
return false;
}
// if (m_verbose) {
// std::cout << std::endl << "POLYGON SPLITTER SUCCESS!" << std::endl << std::endl;
@ -803,7 +802,8 @@ private:
const bool is_event_found = false;
return is_event_found;
CGAL_assertion_msg(false, "TODO: ADD THIS EXTRA TYPE OF EVENT!");
CGAL_assertion_msg(false,
"TODO: ADD PVERTEX TO IVERTEX UNCONSTRAINED EVENT!");
}
const std::size_t run(
@ -904,7 +904,7 @@ private:
const Event& /* event */) {
CGAL_assertion_msg(false,
"TODO: ADD CASE TWO CONSTRAINED PVERTICES MEET!");
"TODO: IMPLEMENT TWO CONSTRAINED PVERTICES MEET EVENT!");
}
void apply_event_two_unconstrained_pvertices_meet(
@ -1192,7 +1192,7 @@ private:
const Event& /* event */) {
CGAL_assertion_msg(false,
"TODO: ADD UNCONSTRAINED PVERTEX MEETS IVERTEX EVENT!");
"TODO: IMPLEMENT UNCONSTRAINED PVERTEX MEETS IVERTEX EVENT!");
}
void remove_events(const IEdge& iedge, const KSR::size_t support_plane_idx) {

35
Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/CMakeLists.txt
View File

@ -19,21 +19,30 @@ if(CGAL_FOUND)
if(Boost_FOUND)
message(STATUS "Found Boost")
set(targets
# kinetic_2d_stress_test
kinetic_3d_test_all
)
find_package(Eigen3 3.1.0 REQUIRED)
if(Eigen3_FOUND)
message(STATUS "Found Eigen")
set(project_linked_libraries)
set(project_compilation_definitions)
include(CGAL_Eigen_support)
foreach(target ${targets})
create_single_source_cgal_program("${target}.cpp")
if(TARGET ${target})
target_link_libraries(${target} PUBLIC ${project_linked_libraries})
target_compile_definitions(${target} PUBLIC ${project_compilation_definitions})
endif()
endforeach()
set(targets
# kinetic_2d_stress_test
kinetic_3d_test_all
)
set(project_linked_libraries)
set(project_compilation_definitions)
foreach(target ${targets})
create_single_source_cgal_program("${target}.cpp")
if(TARGET ${target})
target_link_libraries(${target} PUBLIC ${project_linked_libraries} CGAL::Eigen_support)
target_compile_definitions(${target} PUBLIC ${project_compilation_definitions})
endif()
endforeach()
else()
message(ERROR "This program requires the Eigen library, and will not be compiled.")
endif()
else()
message(ERROR "This program requires the Boost library, and will not be compiled.")
endif()

5
Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_3d_test_all.cpp
View File

@ -139,11 +139,8 @@ int main (const int argc, const char** argv) {
assert(run_test("data/edge-case-test/test-2-polygons.off" , ks, num_iters, num_tests)); // edge touch
assert(run_test("data/edge-case-test/test-4-polygons.off" , ks, num_iters, num_tests)); // edge touch / 2 coplanar
assert(run_test("data/edge-case-test/test-5-polygons.off" , ks, num_iters, num_tests)); // edge touch / vertex touch / 2 coplanar
// std::vector<unsigned int> ts;
// ts.push_back(1); ts.push_back(2); ts.push_back(4);
// ts.push_back(5); ts.push_back(6); ts.push_back(100);
assert(run_test("data/edge-case-test/test-20-polygons.off", ks, num_iters, num_tests)); // 2 overlap and coplanar
assert(run_test("data/edge-case-test/test-flat-bbox.off" , ks, num_iters, num_tests)); // flat bbox / 2 coplanar
std::cout << std::endl << "--OUTPUT STATS:" << std::endl;
std::cout << "* number of iterations per test: " << num_iters << std::endl;