diff options
author | skachano <skachano@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-12-07 14:45:43 +0000 |
---|---|---|
committer | skachano <skachano@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-12-07 14:45:43 +0000 |
commit | c4078affdbf6fac7150c10ade96fcb72270ac013 (patch) | |
tree | 1ad197bb90078a56036a49c6ee3766a032f85e63 /src/Witness_complex/example/protected_sets/protected_sets_paper.cpp | |
parent | f70e386fc98f1dbd8287d1cb7cc715710a8f751b (diff) | |
parent | 061e43a2a48525bc5a69482a1ea80f20ff505e55 (diff) |
Merged with trunk and removed unnecessary files
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@934 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: d0ec52d222d22c102e9fe57590882cd0024c82d5
Diffstat (limited to 'src/Witness_complex/example/protected_sets/protected_sets_paper.cpp')
-rw-r--r-- | src/Witness_complex/example/protected_sets/protected_sets_paper.cpp | 610 |
1 files changed, 0 insertions, 610 deletions
diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp b/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp deleted file mode 100644 index f3df3f1e..00000000 --- a/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp +++ /dev/null @@ -1,610 +0,0 @@ -#ifndef PROTECTED_SETS_H -#define PROTECTED_SETS_H - -#include <algorithm> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Euclidean_distance.h> -#include <CGAL/Kernel_d/Sphere_d.h> -#include <CGAL/Kernel_d/Hyperplane_d.h> -#include <CGAL/Kernel_d/Vector_d.h> - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -typedef K::Vector_d Vector_d; -typedef K::Oriented_side_d Oriented_side_d; -typedef K::Has_on_positive_side_d Has_on_positive_side_d; -typedef K::Sphere_d Sphere_d; -typedef K::Hyperplane_d Hyperplane_d; - -typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; -typedef Delaunay_triangulation::Facet Facet; -typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; -typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; - -typedef std::vector<Point_d> Point_Vector; -typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; - -FT _sfty = pow(10,-14); - -/////////////////////////////////////////////////////////////////////////////////////////////////////////// -// AUXILLARY FUNCTIONS -/////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well - * - * W is the initial point vector - * chosen_landmark is the index of the chosen point in W - * landmarks_ind is the vector of indices of already chosen points in W - * delaunay is the Delaunay triangulation - * landmark_count is the current number of chosen vertices - * torus is true iff you are working on a flat torus [-1,1]^d - * OUT: Vertex handle to the newly inserted point - */ -Delaunay_vertex insert_delaunay_landmark_with_copies(Point_d& p, Delaunay_triangulation& delaunay, int& landmark_count, bool torus) -{ - if (!torus) - { - Delaunay_vertex v =delaunay.insert(p); - landmark_count++; - return v; - } - else - { - int D = W[0].size(); - int nb_cells = pow(3, D); - Delaunay_vertex v; - for (int i = 0; i < nb_cells; ++i) - { - std::vector<FT> point; - int cell_i = i; - for (int l = 0; l < D; ++l) - { - point.push_back(p[l] + 2.0*(cell_i%3-1)); - cell_i /= 3; - } - v = delaunay.insert(point); - } - landmark_count++; - return v; - } -} - -/** Small check if the vertex v is in the full cell fc - */ - -bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc) -{ - for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it) - if (*v_it == v) - return true; - return false; -} - -/** Fill chosen point vector from indices with copies if you are working on a flat torus - * - * IN: W is the point vector - * OUT: landmarks is the output vector - * IN: landmarks_ind is the vector of indices - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus) -{ - if (!torus) - for (unsigned j = 0; j < landmarks_ind.size(); ++j) - landmarks.push_back(W[landmarks_ind[j]]); - else - { - int D = W[0].size(); - int nb_cells = pow(3, D); - int nbL = landmarks_ind.size(); - // Fill landmarks - for (int i = 0; i < nb_cells-1; ++i) - for (int j = 0; j < nbL; ++j) - { - int cell_i = i; - Point_d point; - for (int l = 0; l < D; ++l) - { - point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1)); - cell_i /= 3; - } - landmarks.push_back(point); - } - } -} - -/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices - * - * IN: t is the Delaunay triangulation - * OUT: full_cells is the output vector - */ - -void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) -{ - // Store vertex indices in a map - int ind = 0; //index of a vertex - std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (t.is_infinite(v_it)) - continue; - else - index_of_vertex[v_it] = ind++; - // Write full cells as vectors in full_cells - for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) - { - if (t.is_infinite(fc_it)) - continue; - Point_Vector vertices; - for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) - vertices.push_back((*fc_v_it)->point()); - Sphere_d cs( vertices.begin(), vertices.end()); - Point_d csc = cs.center(); - bool in_cube = true; - for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) - if (*xi > 1.0 || *xi < -1.0) - { - in_cube = false; break; - } - if (!in_cube) - continue; - std::vector<int> cell; - for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) - cell.push_back(index_of_vertex[*v_it]); - full_cells.push_back(cell); - } -} - -//////////////////////////////////////////////////////////////////////////////////////////////////////////// -// IS VIOLATED TEST -//////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Check if a newly created cell is protected from old vertices - * - * t is the Delaunay triangulation - * vertices is the vector containing the point to insert and a facet f in t - * v1 is the vertex of t, such that f and v1 form a simplex - * v2 is the vertex of t, such that f and v2 form another simplex - * delta is the protection constant - * power_protection is true iff the delta-power protection is used - */ - -bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0) -{ - assert(vertices.size() == vertices[0].size() || - vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1 - assert(v1 != v2); - if (vertices.size() == vertices[0].size() + 1) - // FINITE CASE - { - Sphere_d cs(vertices.begin(), vertices.end()); - Point_d center_cs = cs.center(); - FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0])); - /* - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (!t.is_infinite(v_it)) - { - //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); - if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - } - */ - // Check if the simplex is thick enough - Hyperplane_d tau_h(vertices.begin()+1, vertices.end()); - Vector_d orth_tau = tau_h.orthogonal_vector(); - /* - p_s1 = Vector_d(*(vertices.begin()), *(vertices.begin()+1)); - */ - //std::cout << "||orth_tau|| = " << sqrt(orth_tau.squared_length()) << "\n"; - FT orth_length = sqrt(orth_tau.squared_length()); - K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; - // Compute the altitude - FT h = 0; - for (o_it = orth_tau.cartesian_begin(), - p_it = vertices.begin()->cartesian_begin(), - s_it = (vertices.begin()+1)->cartesian_begin(); - o_it != orth_tau.cartesian_end(); - ++o_it, ++p_it, ++s_it) - h += (*o_it)*(*p_it - *s_it)/orth_length; - h = fabs(h); - // Is the center inside the box? - bool inside_the_box = true; - for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it) - if (*c_it > 1.0 || *c_it < -1.0) - { - inside_the_box = false; break; - } - if (inside_the_box && h/r < theta0) - return true; - if (!t.is_infinite(v1)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - if (!t.is_infinite(v2)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - } - else - // INFINITE CASE - { - Delaunay_triangulation::Vertex_iterator v = t.vertices_begin(); - while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end()) - v++; - Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE); - Vector_d orth_v = facet_plane.orthogonal_vector(); - /* - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (!t.is_infinite(v_it)) - if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) - { - std::vector<FT> coords; - Point_d p = v_it->point(); - auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); - for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) - coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); - Point_d p_delta = Point_d(coords); - bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - if (!p_is_inside && p_delta_is_inside) - return true; - } - */ - if (!t.is_infinite(v1)) - { - std::vector<FT> coords; - Point_d p = v1->point(); - auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); - for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) - coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); - Point_d p_delta = Point_d(coords); - bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - if (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - if (!t.is_infinite(v2)) - { - std::vector<FT> coords; - Point_d p = v2->point(); - auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); - for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) - coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); - Point_d p_delta = Point_d(coords); - bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - if (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - } - return false; -} - -/** Auxillary recursive function to check if the point p violates the protection of the cell c and - * if there is a violation of an eventual new cell - * - * p is the point to insert - * t is the current triangulation - * c is the current cell (simplex) - * parent_cell is the parent cell (simplex) - * index is the index of the facet between c and parent_cell from parent_cell's point of view - * D is the dimension of the triangulation - * delta is the protection constant - * marked_cells is the vector of all visited cells containing p in their circumscribed ball - * power_protection is true iff you are working with delta-power protection - * - * OUT: true iff inserting p hasn't produced any violation so far - */ - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0) -{ - Euclidean_distance ed; - std::vector<Point_d> vertices; - if (!t.is_infinite(c)) - { - // if the cell is finite, we look if the protection is violated - for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it) - vertices.push_back((*v_it)->point()); - Sphere_d cs( vertices.begin(), vertices.end()); - Point_d center_cs = cs.center(); - FT r = sqrt(ed.transformed_distance(center_cs, vertices[0])); - FT dist2 = ed.transformed_distance(center_cs, p); - // if the new point is inside the protection ball of a non conflicting simplex - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - // if the new point is inside the circumscribing ball : continue violation searching on neighbours - //if (dist2 < r*r) - //if (dist2 < (5*r+delta)*(5*r+delta)) - if (dist2 < r*r) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - for (int i = 0; i < D+1; ++i) - { - Full_cell_handle next_c = c->neighbor(i); - if (next_c->tds_data().is_clear() && - is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0)) - return true; - } - } - // if the new point is outside the protection sphere - else - { - // facet f is on the border of the conflict zone : check protection of simplex {p,f} - // the new simplex is guaranteed to be finite - vertices.clear(); vertices.push_back(p); - for (int i = 0; i < D+1; ++i) - if (i != index) - vertices.push_back(parent_cell->vertex(i)->point()); - Delaunay_vertex vertex_to_check = t.infinite_vertex(); - for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) - if (!vertex_is_in_full_cell(*vh_it, parent_cell)) - { - vertex_to_check = *vh_it; break; - } - if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) - //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) - return true; - } - } - else - { - // Inside of the convex hull is + side. Outside is - side. - for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) - if (!t.is_infinite(*vh_it)) - vertices.push_back((*vh_it)->point()); - Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin(); - while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c)) - v_it++; - Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE); - //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point()); - Vector_d orth_v = facet_plane.orthogonal_vector(); - std::vector<FT> coords; - auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); - for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) - coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); - Point_d p_delta = Point_d(coords); - bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - - // If we work with power protection, we just ignore any conflicts - if (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - //if the cell is infinite we look at the neighbours regardless - if (p_is_inside) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - for (int i = 0; i < D+1; ++i) - { - Full_cell_handle next_c = c->neighbor(i); - if (next_c->tds_data().is_clear() && - is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0)) - return true; - } - } - else - { - // facet f is on the border of the conflict zone : check protection of simplex {p,f} - // the new simplex is finite if the parent cell is finite - vertices.clear(); vertices.push_back(p); - for (int i = 0; i < D+1; ++i) - if (i != index) - if (!t.is_infinite(parent_cell->vertex(i))) - vertices.push_back(parent_cell->vertex(i)->point()); - Delaunay_vertex vertex_to_check = t.infinite_vertex(); - for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) - if (!vertex_is_in_full_cell(*vh_it, parent_cell)) - { - vertex_to_check = *vh_it; break; - } - if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) - //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) - return true; - } - } - //c->tds_data().clear_visited(); - //marked_cells.pop_back(); - return false; -} - -/** Checks if inserting the point p in t will make conflicts - * - * p is the point to insert - * t is the current triangulation - * D is the dimension of triangulation - * delta is the protection constant - * power_protection is true iff you are working with delta-power protection - * OUT: true iff inserting p produces a violation of delta-protection. - */ - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0) -{ - Euclidean_distance ed; - Delaunay_triangulation::Vertex_handle v; - Delaunay_triangulation::Face f(t.current_dimension()); - Delaunay_triangulation::Facet ft; - Delaunay_triangulation::Full_cell_handle c; - Delaunay_triangulation::Locate_type lt; - std::vector<Full_cell_handle> marked_cells; - c = t.locate(p, lt, f, ft, v); - bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells, power_protection, theta0); - for (Full_cell_handle fc : marked_cells) - fc->tds_data().clear(); - return violation_existing_cells; -} - -////////////////////////////////////////////////////////////////////// -// INITIALIZATION -////////////////////////////////////////////////////////////////////// - -void initialize(Search_Tree& W, Delaunay& t, int D, int width, bool torus) -{ - if (!torus) - std::cout << "Non-toric case is not supported\n"; - else - { - if (D == 2) - { - FT stepx = 2.0/width; - FT stepy = sqrt(3)/width; - for (int i = 0; i < width; ++i) - for (int j = 0; j < floor(2*width/sqrt(3)); ++j) - { - insert_delaunay_landmark_with_copies(Point_d(step*i,)) - } - } - else (D == 3) - { - - } - else std::cout << "T^d with d>3 not supported"; - } -} - -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// -//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!// -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// - -/////////////////////////////////////////////////////////////////////// -// LANDMARK CHOICE PROCEDURE AS IN PAPER -/////////////////////////////////////////////////////////////////////// - -/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call. - * - * IN: W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d - * IN: nbP is the size of W - * OUT: landmarks is the output vector for the points - * OUT: landmarks_ind is the output vector for the indices of the selected points in W - * IN: delta is the constant of protection - * OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -template<class Search_Tree> -void protected_delaunay_refinement(Search_Tree& W, int nbP, Point_Vector& landmarks, FT delta, bool torus, bool power_protection, FT theta0) -{ - bool return_ = true; - unsigned D = W[0].size(); - Torus_distance td; - Euclidean_distance ed; - Delaunay_triangulation t(D); - CGAL::Random rand; - int landmark_count = 0; - //std::list<int> index_list; - // shuffle the list of indexes (via a vector) - // { - // std::vector<int> temp_vector; - // for (int i = 0; i < nbP; ++i) - // temp_vector.push_back(i); - // unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - // std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); - // //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end()); - // for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) - // index_list.push_front(*it); - // } - if (torus) - if (D == 2) - // \T^2 - { - for (int i = 0; i < 4; ++i) - for (int j = 0; j < 2; ++j) - { - W[index_list.front()] = Point_d(std::vector<FT>{i*0.5, j*1.0}); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); - index_list.pop_front(); - W[index_list.front()] = Point_d(std::vector<FT>{0.25+i*0.5, 0.5+j}); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); - index_list.pop_front(); - } - } - else if (D == 3) - { - - } - //std::cout << "No torus starter available for dim>2\n"; - std::list<int>::iterator list_it = index_list.begin(); - while (list_it != index_list.end()) - { - if (!is_violating_protection(W[*list_it], t, D, delta, power_protection, theta0)) - { - // If no conflicts then insert in every copy of T^3 - - insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count, torus); - if (return_) - { - index_list.erase(list_it); - list_it = index_list.begin(); - } - else - index_list.erase(list_it++); - /* - // PIECE OF CODE FOR DEBUGGING PURPOSES - - Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); - if (triangulation_is_protected(t, delta)) - { - index_list.erase(list_it); - list_it = index_list.begin(); - } - else - { //THAT'S WHERE SOMETHING'S WRONG - t.remove(inserted_v); - landmarks_ind.pop_back(); - landmark_count--; - write_delaunay_mesh(t, W[*list_it], is2d); - is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore - } - */ - //std::cout << "index_list_size() = " << index_list.size() << "\n"; - } - else - { - list_it++; - //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; - } - //if (list_it != index_list.end()) - // write_delaunay_mesh(t, W[*list_it], is2d); - } - fill_landmarks(W, landmarks, landmarks_ind, torus); - fill_full_cell_vector(t, full_cells); - /* - if (triangulation_is_protected(t, delta)) - std::cout << "Triangulation is ok\n"; - else - { - std::cout << "Triangulation is BAD!! T_T しくしく!\n"; - } - */ - //write_delaunay_mesh(t, W[0], is2d); - //std::cout << t << std::endl; -} - -#endif |