CGAL найти точки интерьера в сетке
Я пытаюсь использовать CGAL, чтобы найти точки, которые находятся внутри треугольной сетки. (т.е. набор точек, которые не находятся на его границе. Здесь есть аналогичный пример для трехмерной сетки, которая использует функцию CGAL::Side_of_triangle_mesh<>. Кто-нибудь может помочь / посоветовать, как изменить это для 2D триангуляции?
Мой тестовый код просто создает два квадрата, один внутри другого, затем помещает начальное число в начало координат (чтобы сделать отверстие) и затем выполняет 2D триангуляцию Делоне. Когда я вызываю класс side_of_triangle_mesh<> с помощью:
Point_2 p = points2D[i];
CGAL::Bounded_side res = inside2D(p);
Я получаю ошибку:
/usr/local/include/CGAL/Side_of_triangle_mesh.h:164:16: Candidate function not viable: no known conversion from 'Point_2' (aka 'Point_2<CGAL::Epick>') to 'const Point' (aka 'const Point_3<CGAL::Epick>') for 1st argument
Означает ли это, что side_of_triangle_mesh работает только для сетки 3D Polyhedron_3? Если да, то может ли кто-нибудь предложить способ сделать это для 2D-сетки?
Мой тестовый код ниже: Спасибо
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_mesh_vertex_base_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <CGAL/Side_of_triangle_mesh.h>
#include <vector>
#include <CGAL/Delaunay_mesher_2.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef CGAL::Triangle_2<K> triangle;
typedef CGAL::Delaunay_mesh_vertex_base_2<K> Vb;
typedef CGAL::Delaunay_mesh_face_base_2<K> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT;
typedef CDT::Vertex_handle Vertex_handle;
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;
int main(int argc, char* argv[])
{
// Create a vector of the points
//
std::vector<Point_2> points2D ;
points2D.push_back(Point_2(-5.0, -5.0)); // ----------
points2D.push_back(Point_2( 5.0, -5.0)); // | OUTER
points2D.push_back(Point_2( 5.0, 5.0)); // | SQUARE
points2D.push_back(Point_2(-5.0, 5.0)); // ----------
points2D.push_back(Point_2(-2.5, -2.5)); // ----------
points2D.push_back(Point_2( 2.5, -2.5)); // | INNER
points2D.push_back(Point_2( 2.5, 2.5)); // | SQUARE
points2D.push_back(Point_2(-2.5, 2.5)); // ----------
size_t numTestPoints = points2D.size();
// Create a constrained delaunay triangulation and add the points
//
CDT cdt;
std::vector<Vertex_handle> vhs;
for (unsigned int i=0; i<numTestPoints; ++i){
vhs.push_back(cdt.insert(points2D[i]));
}
// Creare constraints of the sides of the mesh
//
cdt.insert_constraint(vhs[0],vhs[1]);
cdt.insert_constraint(vhs[1],vhs[2]);
cdt.insert_constraint(vhs[2],vhs[3]);
cdt.insert_constraint(vhs[3],vhs[0]);
cdt.insert_constraint(vhs[4],vhs[5]);
cdt.insert_constraint(vhs[5],vhs[6]);
cdt.insert_constraint(vhs[6],vhs[7]);
cdt.insert_constraint(vhs[7],vhs[4]);
// Create a seed to make sure the inner square is a hole
//
std::list<Point_2> list_of_seeds;
list_of_seeds.push_back(Point_2(0, 0));
// Refine the mesh into triangles of max side length '1' and ensuring seeds are 'holes'
//
CGAL::refine_Delaunay_mesh_2(cdt, list_of_seeds.begin(), list_of_seeds.end(),
Criteria(0.125, 1),false);
// Call side_of_triangle_mesh
//
CGAL::Side_of_triangle_mesh<CDT, K> inside2D(cdt) ;
int n_inside = 0;
int n_boundary = 0;
for(unsigned int i=0; i<numTestPoints; ++i)
{
Point_2 p = points2D[i];
CGAL::Bounded_side res = inside2D(p);
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
// NO MATCHING FUNCTION Call
if (res == CGAL::ON_BOUNDED_SIDE) { ++n_inside; }
if (res == CGAL::ON_BOUNDARY) { ++n_boundary; }
}
std::cerr << "2D results for query size: " << cdt.number_of_vertices() << std::endl;
std::cerr << " " << n_inside << " points inside " << std::endl;
std::cerr << " " << n_boundary << " points on boundary " << std::endl;
return 0 ;
}
2 ответа
Решение
Благодаря @sloriot я дошел до сути. Суть в том, что вы не можете использовать CGAL::Side_of_triangle_mesh<> для двумерных сеток Делоне. Вместо этого вам нужно пройти по всем вершинам сетки и проверить каждую из них, просматривая все соседние грани вокруг точки вершины. если какая-либо из граней находится за пределами области, то вершина или точка находится на краю сетки.
Вот исправленный код:
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_mesh_vertex_base_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <vector>
#include <CGAL/Delaunay_mesher_2.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef CGAL::Delaunay_mesh_vertex_base_2<K> Vb;
typedef CGAL::Delaunay_mesh_face_base_2<K> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT;
typedef CDT::Vertex_handle Vertex_handle;
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;
typedef CDT::Vertex_iterator Vertex_iterator;
typedef CDT::Vertex Vertex;
typedef CDT::Face Face ;
typedef Face::Face_handle Face_handle ;
typedef CDT::Face_circulator Face_circulator ;
int main(int argc, char* argv[])
{
// Create a vector of the points
//
std::vector<Point_2> points2D ;
points2D.push_back(Point_2(-5.0, -5.0)); // ----------
points2D.push_back(Point_2( 5.0, -5.0)); // | OUTER
points2D.push_back(Point_2( 5.0, 5.0)); // | SQUARE
points2D.push_back(Point_2(-5.0, 5.0)); // ----------
points2D.push_back(Point_2(-2.5, -2.5)); // ----------
points2D.push_back(Point_2( 2.5, -2.5)); // | INNER
points2D.push_back(Point_2( 2.5, 2.5)); // | SQUARE
points2D.push_back(Point_2(-2.5, 2.5)); // ----------
size_t numTestPoints = points2D.size();
// Create a constrained delaunay triangulation and add the points
//
CDT cdt;
std::vector<Vertex_handle> vhs;
for (unsigned int i=0; i<numTestPoints; ++i){
vhs.push_back(cdt.insert(points2D[i]));
}
// Creare constraints of the sides of the mesh
//
cdt.insert_constraint(vhs[0],vhs[1]);
cdt.insert_constraint(vhs[1],vhs[2]);
cdt.insert_constraint(vhs[2],vhs[3]);
cdt.insert_constraint(vhs[3],vhs[0]);
cdt.insert_constraint(vhs[4],vhs[5]);
cdt.insert_constraint(vhs[5],vhs[6]);
cdt.insert_constraint(vhs[6],vhs[7]);
cdt.insert_constraint(vhs[7],vhs[4]);
// Create a seed to make sure the inner square is a hole
//
std::list<Point_2> list_of_seeds;
list_of_seeds.push_back(Point_2(0, 0));
// Refine the mesh into triangles of max side length '1' and ensuring seeds are 'holes'
//
CGAL::refine_Delaunay_mesh_2(cdt, list_of_seeds.begin(), list_of_seeds.end(),
Criteria(0.125, 1.5),false);
// The mesh is now created. The next bit swings around each vertex point checking that
// all faces around it are in the domain. If any are not then the vertex is on the
// edge of the mesh.
// thanks to @sloriot for this
//
Vertex v ;
std::vector<Point_2> interior_points ;
std::vector<Point_2> boundary_points ;
CDT::Locate_type loc = CDT::Locate_type::VERTEX ;
int li;
Vertex_iterator vit = cdt.vertices_begin(), beyond = cdt.vertices_end() ;
while (vit != beyond) {
v = *vit ;
Point_2 query = vit->point();
Face_handle f = cdt.locate(query, loc, li);
bool current_face_in_domain ;
Face_handle start = f ;
do {
current_face_in_domain = f->is_in_domain() ;
Vertex_handle vh = f->vertex(li);
f = f->neighbor(cdt.ccw(li));
li = f->index(vh) ;
}
while(current_face_in_domain && f != start) ;
if (f == start) {
interior_points.push_back(query) ;
}else{
boundary_points.push_back(query) ;
}
++vit ;
}
printf("Boundary points:\n");
for(auto p = boundary_points.begin(); p != boundary_points.end(); ++p){
printf("%f,%f\n",p->x(), p->y()) ;
}
printf("interior points:\n");
for(auto p = interior_points.begin(); p != interior_points.end(); ++p){
printf("%f,%f\n",p->x(), p->y()) ;
}
return 0 ;
}
Когда вы запустите его, вы должны получить что-то вроде:
Вы должны использовать locate функция. Он вернет дескриптор лица, содержащий вашу точку запроса. Затем вам нужно проверить, находится ли лицо в домене или нет, используя is_in_domain() функция-член.
Обратите внимание, что если у вас есть несколько запросов, вы должны сначала поместить все точки запроса в контейнер, отсортировать их, используя hilbert_sort и использовать местоположение предыдущей точки в качестве второго параметра функции locate.
Вот пример кода того, как получить грани инцидента в граничных случаях:
CDT_2 mesh;
CGAL::Locate_type loc;
int li;
Point_2 query;
Face_handle f = mesh.locate(query, loc, li);
switch(loc)
{
case FACE:
f->is_in_domain();
break;
case EDGE:
{
bool face_1_in_domain = f->is_in_domain(); // first incident face status
bool face_2_in_domain = f->neighbor(li)->is_in_domain(); // second incident face status
break;
}
case VERTEX:
{
// turn around f->vertex(li)
Face_handle start = f;
do{
bool current_face_in_domain = f->is_in_domain();
Vertex_handle v = f->vertex(mesh.cw(li));
f = f->neighbor(mesh.ccw(li));
li = f->index(v);
}
while(f!=current);
break;
}
default:
// outside of the domain.
}