Skip to content

Commit

Permalink
Fixed artifact by adiing uniformly distributed generated points in th…
Browse files Browse the repository at this point in the history 
…e sphere to the triangulations
  • Loading branch information
@efifogel
efifogel committed Apr 15, 2024
1 parent 92a2139 commit eb330d6
Show file tree
Hide file tree
Showing 4 changed files with 101 additions and 60 deletions.
145 changes: 89 additions & 56 deletions Arrangement_on_surface_2/demo/Arrangement_on_surface_2_earth/Aos_triangulator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,11 @@

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/draw_triangulation_2.h>
// #include <CGAL/draw_triangulation_2.h>
#include <CGAL/mark_domain_in_triangulation.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Arr_batched_point_location.h>

#include "Aos_defs.h"
#include "Aos_triangulator.h"
Expand Down Expand Up @@ -64,8 +66,7 @@ std::vector<QVector3D> Aos_triangulator::get_all(Aos::Arr_handle arrh) {
// loop on all faces of the arrangement
for (auto fh : arr.face_handles()) {
// skip any face with no OUTER-CCB
if (0 == fh->number_of_outer_ccbs())
continue;
if (0 == fh->number_of_outer_ccbs()) continue;

// COMPUTE THE CENTROID OF ALL FACE-POINTS
std::vector<QVector3D> face_points;
Expand All @@ -84,13 +85,13 @@ std::vector<QVector3D> Aos_triangulator::get_all(Aos::Arr_handle arrh) {
}

// RESULTING TRIANGLE POINTS (every 3 point => triangle)
std::vector<QVector3D> triangles;
std::vector<QVector3D> triangles;

std::cout << "triangulating individual faces\n";

// loop on all approximated faces
for (auto& face_points : all_faces) {
//std::cout << "num face points = " << face_points.size() << std::endl;
std::cout << "num face points = " << face_points.size() << std::endl;
// no need to triangulate if the number of points is 3
if (face_points.size() == 3) {
triangles.insert(triangles.end(), face_points.begin(), face_points.end());
Expand Down Expand Up @@ -141,7 +142,7 @@ std::vector<QVector3D> Aos_triangulator::get_all(Aos::Arr_handle arrh) {
// loop on all the triangles ("faces" in triangulation doc)
for (Face_handle f : cdt.finite_face_handles()) {
// if the current triangles is not inside the polygon -> skip it
if (false == get(in_domain, f)) continue;
if (! get(in_domain, f)) continue;

for (int i = 0; i < 3; ++i) {
auto tp = f->vertex(i)->point();
Expand All @@ -157,48 +158,77 @@ std::vector<QVector3D> Aos_triangulator::get_all(Aos::Arr_handle arrh) {
}

//! \brief
Country_triangles_map Aos_triangulator::get_by_country(Aos::Arr_handle arrh,
float error) {
Country_triangles_map
Aos_triangulator::get_by_country(Aos::Arr_handle arrh, float error,
std::size_t num_uniform_points) {

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Pnt_2 = typename K::Point_2;
using Pnt_3 = typename K::Point_3;
using Vb = CGAL::Triangulation_vertex_base_2<K>;
using Fb = CGAL::Constrained_triangulation_face_base_2<K>;
using TDS = CGAL::Triangulation_data_structure_2<Vb, Fb>;
using Itag = CGAL::Exact_predicates_tag;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<K, TDS, Itag>;
using Face_handle = CDT::Face_handle;
// using Point = CDT::Point;
using Polygon_2 = CGAL::Polygon_2<K>;
using Approximate_point_2 = Geom_traits::Approximate_point_2;
using Aos = Countries_arr;
using Aos_pnt = Aos::Point_2;

auto& arr = *reinterpret_cast<Countries_arr*>(arrh.get());
auto& arr = *reinterpret_cast<Aos*>(arrh.get());
auto approx = s_traits.approximate_2_object();

// group the faces by their country name
using Face_ = Countries_arr::Face_handle::value_type;
std::map<std::string, std::vector<Face_*>> country_faces_map;
using Face_const_handle = Countries_arr::Face_const_handle;
std::map<std::string, std::vector<Face_const_handle>> country_faces_map;
for (auto fh : arr.face_handles()) {
auto& face = *fh;
const auto& country_name = fh->data();
// skipping spherical-face
if (country_name.empty()) continue;
country_faces_map[country_name].push_back(&face);
country_faces_map[country_name].push_back(fh);
}

std::cout << "triangulating individual faces\n";
// Generate random point uniformly distributed
std::vector<Pnt_3> rps;
rps.reserve(num_uniform_points);
CGAL::Random_points_in_sphere_3<Pnt_3> g(1.0);
std::copy_n(g, rps.capacity(), std::back_inserter(rps));
std::map<Face_const_handle, std::vector<Pnt_3>> face_random_points;
using Point_location_result = CGAL::Arr_point_location_result<Aos>;
using Query_result = std::pair<Aos_pnt, Point_location_result::Type>;
std::vector<Aos_pnt> rqs(rps.size());
auto ctr_p = arr.geometry_traits()->construct_point_2_object();
std::transform(rps.begin(), rps.end(), rqs.begin(),
[&](const Pnt_3& rp) -> Aos_pnt {
return ctr_p(rp.x(), rp.y(), rp.z());
});
std::list<Query_result> results;
CGAL::locate(arr, rqs.begin(), rqs.end(), std::back_inserter(results));
for (auto& [q, res] : results) {
const Aos::Face_const_handle* f;
if ((f = std::get_if<Face_const_handle>(&res))) {
auto x = CGAL::to_double(q.dx());
auto y = CGAL::to_double(q.dy());
auto z = CGAL::to_double(q.dz());
face_random_points[*f].push_back(Pnt_3(x, y, z));
}
}

// Triangulate the faces
Country_triangles_map result;
for (auto& [country_name, faces] : country_faces_map) {
for (auto& [country_name, fhs] : country_faces_map) {
// std::cout << "processing country " << country_name << std::endl;
auto& triangles = result[country_name];
// CONVERT the face-points to QVector3D
using Face_points = std::vector<QVector3D>;
using Faces_ = std::vector<Face_points>;
Faces_ all_faces_of_current_country;
for (auto* face : faces) {
for (auto fh : fhs) {
// skip any face with no OUTER-CCB
if (0 == face->number_of_outer_ccbs()) continue;
if (0 == fh->number_of_outer_ccbs()) continue;

std::vector<QVector3D> face_points;
// loop on the egdes of the current outer-ccb
auto first = face->outer_ccb();
// Loop on the egdes of the current outer-ccb
auto first = fh->outer_ccb();
auto curr = first;
do {
std::vector<Approximate_point_2> apx_points;
Expand All @@ -212,65 +242,68 @@ Country_triangles_map Aos_triangulator::get_by_country(Aos::Arr_handle arrh,
}
} while (++curr != first);

all_faces_of_current_country.push_back(std::move(face_points));
}

// RESULTING TRIANGLE POINTS (every 3 point => triangle)
auto& triangles = result[country_name];
// loop on all approximated faces
for (auto& face_points : all_faces_of_current_country) {
// no need to triangulate if the number of points is 3
if (face_points.size() == 3) {
triangles.insert(triangles.end(), face_points.begin(),
face_points.end());
continue;
}

// find the centroid of all face-points
// Calculate the centroid of all face-points
QVector3D centroid(0, 0, 0);
for (const auto& fp : face_points) centroid += fp;
centroid /= face_points.size();
centroid.normalize();
auto normal = centroid;

K::Point_3 plane_origin(centroid.x(), centroid.y(), centroid.z());
Pnt_3 plane_origin(centroid.x(), centroid.y(), centroid.z());
K::Vector_3 plane_normal(normal.x(), normal.y(), normal.z());
K::Plane_3 plane(plane_origin, plane_normal);

Polygon_2 polygon;

// project all points onto the plane
K::Point_3 origin(0, 0, 0);
for (const auto& fp : face_points) {
// define a ray through the origin and the current point
K::Point_3 current_point(fp.x(), fp.y(), fp.z());
K::Ray_3 ray(origin, current_point);

auto intersection = CGAL::intersection(plane, ray);
if (! intersection.has_value())
std::cout << "INTERSECTION ASSERTION ERROR!!!\n";
auto ip = std::get<K::Point_3>(intersection.value());
auto ip2 = plane.to_2d(ip);
// Project all points onto the plane
std::vector<Pnt_2> planar_points(face_points.size());
std::transform(face_points.begin(), face_points.end(),
planar_points.begin(),
[&](const QVector3D& fp) -> Pnt_2 {
// define a ray through the origin and the current point
Pnt_3 p3(fp.x(), fp.y(), fp.z());
K::Ray_3 ray(CGAL::ORIGIN, p3);
auto intersection = CGAL::intersection(plane, ray);
if (! intersection.has_value())
std::cout << "INTERSECTION ASSERTION ERROR!!!\n";
auto ip = std::get<Pnt_3>(intersection.value());
return plane.to_2d(ip);
});
CDT cdt;

// add this to the polygon constraint
polygon.push_back(ip2);
// Insert points uniformly distributed into the triangulation
auto it = face_random_points.find(fh);
if (it != face_random_points.end()) {
const auto& points = it->second;
for (const auto& p3 : points) {
K::Ray_3 ray(CGAL::ORIGIN, p3);
auto intersection = CGAL::intersection(plane, ray);
if (! intersection.has_value())
std::cout << "INTERSECTION ASSERTION ERROR!!!\n";
auto ip = std::get<Pnt_3>(intersection.value());
auto p2 = plane.to_2d(ip);
cdt.insert(p2);
}
}

CDT cdt;
cdt.insert_constraint(polygon.vertices_begin(), polygon.vertices_end(),
true);
// Insert the constraints into the triangulation
cdt.insert_constraint(planar_points.begin(), planar_points.end(), true);

// Mark facets that are inside the domain bounded by the polygon
std::unordered_map<Face_handle, bool> in_domain_map;
boost::associative_property_map< std::unordered_map<Face_handle, bool> >
boost::associative_property_map< std::unordered_map<Face_handle, bool>>
in_domain(in_domain_map);

// Mark facets that are inside the domain bounded by the polygon
CGAL::mark_domain_in_triangulation(cdt, in_domain);

// loop on all the triangles ("faces" in triangulation doc)
// Loop on all the triangles ("faces" in triangulation doc)
for (Face_handle f : cdt.finite_face_handles()) {
// if the current triangles is not inside the polygon -> skip it
if (false == get(in_domain, f)) continue;
// If the current triangles is not inside the polygon -> skip it
if (! get(in_domain, f)) continue;

for (int i = 0; i < 3; ++i) {
auto tp = f->vertex(i)->point();
Expand Down
5 changes: 3 additions & 2 deletions Arrangement_on_surface_2/demo/Arrangement_on_surface_2_earth/Aos_triangulator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,9 @@ class Aos_triangulator {
public:
static std::vector<QVector3D> get_all(Aos::Arr_handle arrh);

static Country_triangles_map get_by_country(Aos::Arr_handle arrh, float error);
static Country_triangles_map get_by_country(Aos::Arr_handle arrh,
float error,
std::size_t num_uniform_points);
};


#endif
10 changes: 8 additions & 2 deletions Arrangement_on_surface_2/demo/Arrangement_on_surface_2_earth/Main_widget.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,7 @@ void Main_widget::initializeGL() {
init_shader_programs();

m_current_approx_error = 0.001f;
m_num_uniform_points = 4096;

qDebug() << "loading arrangement..";
m_arrh = Aos::load_arr(m_file_name.toStdString());
Expand All @@ -130,7 +131,8 @@ void Main_widget::initializeGL() {
//auto triangle_points = Aos_triangulator::get_all(arrh);
//auto country_triangles_map = Aos::get_triangles_by_country(m_arrh);
auto country_triangles_map =
Aos_triangulator::get_by_country(m_arrh, m_current_approx_error);
Aos_triangulator::get_by_country(m_arrh, m_current_approx_error,
m_num_uniform_points);
//auto color_map = Aos::get_color_mapping(m_arrh);
//qDebug() << "color map size = " << color_map.size();
qDebug() << "num countries = " << country_triangles_map.size();
Expand Down Expand Up @@ -267,11 +269,15 @@ void draw_safe(T& ptr) { if (ptr) ptr->draw(); }
void Main_widget::paintGL() {
if (m_update_approx_error) {
const auto error = compute_backprojected_error(0.5);
qDebug() << "new approx error = " << error;
qDebug() << "current error = " << m_current_approx_error;
qDebug() << "new approx error = " << error;
if (error < m_current_approx_error) {
init_country_borders(error);
m_current_approx_error = error;
auto ratio = static_cast<double>(m_current_approx_error) / error;
m_num_uniform_points *= ratio*ratio;
std::cout << "num uniform points = " << m_num_uniform_points
<< std::endl;
}
m_update_approx_error = false;
}
Expand Down
1 change: 1 addition & 0 deletions Arrangement_on_surface_2/demo/Arrangement_on_surface_2_earth/Main_widget.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,7 @@ class Main_widget : public QOpenGLWidget, protected OpenGLFunctionsBase {
// INSIDE the paintGL (or whereever the OpenGL context is active)!
bool m_update_approx_error = false;
float m_current_approx_error;
std::size_t m_num_uniform_points;

// Timer for continuous screen-updates
// QBasicTimer m_timer;
Expand Down

0 comments on commit eb330d6

Please sign in to comment.