CGAL users discussion list

[Danyang Su <>]

Dear All,

I want to check if a point is inside a polyhedron by the example Polygon_mesh_processing/point_inside_example.cpp

http://doc.cgal.org/latest/Polygon_mesh_processing/Polygon_mesh_processing_2point_inside_example_8cpp-example.html

The actual code section used in my program is as follows

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point3d;
...
    for (subdomainIterator=subdomains.begin();
            subdomainIterator!=subdomains.end();++subdomainIterator){
        CGAL::Side_of_triangle_mesh<Polyhedron, K> inside(subdomainIterator->polyHedron());

        std::cout << "Number of vertices in polyhedron " << subdomainIterator->polyHedron().size_of_vertices() << endl;
        std::cout << "Number of facets in polyhedron " << subdomainIterator->polyHedron().size_of_facets() << endl;
        std::cout << "Number of halfedges in polyhedron " << subdomainIterator->polyHedron().size_of_halfedges() << endl;
        std::cout << "is polyhedron closed " << subdomainIterator->polyHedron().is_closed() << endl;
        std::cout << "is polyhedron pure triangle " << subdomainIterator->polyHedron().is_pure_triangle() << endl;
        std::cout << "is polyhedron pure quadrilaterals " << subdomainIterator->polyHedron().is_pure_quad() << endl;

        for (int iz=0; iz<zdim; iz++){
            for (int iy=0; iy<ydim; iy++){
                for (int ix=0; ix<xdim; ix++){
                    int idx = iz*xdim*ydim + iy*xdim + ix;
                    double x = ix*dxmin + xcoord;
                    double y = iy*dymin + ycoord;
                    double z = iz*dzmin + zcoord;
                    Point3d pt(x,y,z);
                    CGAL::Bounded_side res = inside(pt);
                    if (res == CGAL::ON_BOUNDED_SIDE){
                        std::cout << "idx " << idx << std::endl;
                        data[idx] = subdomainIterator->DomainLabel();
                    }
                }
            }
        }
    }

When run the program, I get Segmentation fault at the line CGAL::Bounded_side res = inside(pt). The output of the polyhedron is as follows? Would anyone tell me what's wrong in my code?

Number of vertices in polyhedron 7
Number of facets in polyhedron 10
Number of halfedges in polyhedron 30
is polyhedron closed 1
is polyhedron pure triangle 1
is polyhedron pure quadrilaterals 0
Segmentation fault (From CGAL::Bounded_side res = inside(pt))

Thanks and Regards,

Danyang