test run of half space_intersection_3

test.cpp

#include <stdio.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Point_3.h>
#include <CGAL/Plane_3.h>
#include <CGAL/Vector_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Segment_3.h>
#include <CGAL/squared_distance_3.h>
#include <CGAL/Nef_polyhedron_3.h>
#include <CGAL/Point_2.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/Nef_polyhedron_2.h>
#include <CGAL/Extended_cartesian.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_incremental_builder_3.h>
#include <CGAL/Convex_hull_3/dual/halfspace_intersection_3.h>


typedef CGAL::Exact_predicates_exact_constructions_kernel                           K;
typedef CGAL::Nef_polyhedron_3<K>                                                  Nef;
typedef K::Plane_3                                                                 Plane;
typedef CGAL::Polyhedron_3<K>                                                      Polyhdedron;
typedef K::Point_3                                                                 Point;
typedef Polyhdedron::HalfedgeDS                                                    HalfedgeDS;


template <class HDS>
class Build_box : public CGAL::Modifier_base<HDS> {
public:
    Build_box() {}
    void operator()( HDS& hds) {
        // Postcondition: hds is a valid polyhedral surface.
        CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
        B.begin_surface( 3, 1, 6);
        typedef typename HDS::Vertex   Vertex;
        typedef typename Vertex::Point Point;
        
        B.add_vertex(Point(-74.3247,-207.04535,-200.456825));  //0
        B.add_vertex(Point(74.3247,-207.04535,-200.456825));   //1
        B.add_vertex(Point(74.3247,158.60275,-200.456825));    //2
        B.add_vertex(Point(-74.3247,158.60275,-200.456825));   //3
        
        B.add_vertex(Point(-74.3247,-207.04535,201.188725));   //4
        B.add_vertex(Point(74.3247,-207.04535,201.188725));    //5
        B.add_vertex(Point(74.3247,158.60275,201.188725));     //6
        B.add_vertex(Point(-74.3247,158.60275,201.188725));    //7

        
        B.begin_facet();
        B.add_vertex_to_facet( 6);
        B.add_vertex_to_facet( 7);
        B.add_vertex_to_facet( 3);
        B.add_vertex_to_facet( 2);
        B.end_facet();
        
        B.begin_facet();
        B.add_vertex_to_facet( 0);
        B.add_vertex_to_facet( 4);
        B.add_vertex_to_facet( 5);
        B.add_vertex_to_facet( 1);
        B.end_facet();
        
        B.begin_facet();
        B.add_vertex_to_facet( 0);
        B.add_vertex_to_facet( 1);
        B.add_vertex_to_facet( 2);
        B.add_vertex_to_facet( 3);
        B.end_facet();
        
        B.begin_facet();
        B.add_vertex_to_facet( 0);
        B.add_vertex_to_facet( 3);
        B.add_vertex_to_facet( 7);
        B.add_vertex_to_facet( 4);
        B.end_facet();
        
        B.begin_facet();
        B.add_vertex_to_facet( 5);
        B.add_vertex_to_facet( 6);
        B.add_vertex_to_facet( 2);
        B.add_vertex_to_facet( 1);
        B.end_facet();
        
        B.begin_facet();
        B.add_vertex_to_facet( 4);
        B.add_vertex_to_facet( 7);
        B.add_vertex_to_facet( 6);
        B.add_vertex_to_facet( 5);
        B.end_facet();
        
        B.end_surface();
    }
};

struct Plane_equations {	// only for convex facets
    template <class Facet>
    typename Facet::Plane_3 operator()( Facet& f) {
        typename Facet::Halfedge_handle h = f.halfedge();
        typedef typename Facet::Plane_3  Plane;
        return Plane( h->vertex()->point(),
                     h->next()->next()->vertex()->point(),
                     h->next()->vertex()->point()
                     //order reversed so that polyhderon is an intersection of the halfspaces of its face planes
                     );
    }
};


int main(){
    Polyhdedron P;
    Build_box<HalfedgeDS> triangle;
    P.delegate( triangle);
    CGAL_assertion( P.is_pure_quad());

    std::transform( P.facets_begin(), P.facets_end(),P.planes_begin(),Plane_equations());

    //sample testplane to cut polyhderon made by Build_box
    Plane testplane(97.22703654, -0.5915931443, -0.2859889977, -198.4034144);  // almost parallel to yz plane, at x=-2
    
    //list contaiing planes for 2 halfspace intersections, for both sides of the plane
    std::list<Plane> listplanes;
    for (Polyhdedron::Plane_iterator pol=P.planes_begin(); pol!=P.planes_end(); pol++) {
        listplanes.push_back(*pol);
    }
    
    Polyhdedron tempPoly;
    
    listplanes.push_back(testplane.opposite());
    CGAL::halfspace_intersection_3(listplanes.begin(), listplanes.end(), tempPoly,boost::none);
    listplanes.pop_back();
    
    listplanes.push_back(testplane);
    CGAL::halfspace_intersection_3(listplanes.begin(), listplanes.end(), tempPoly,boost::none); // gives error
    listplanes.pop_back();
    
    return 0;
}

74.2757 -206.925 -200.324 (passes test)
-74.2731 158.476 201.049 (gives error)

CGAL error: assertion violation!
Expression : boost::is_floating_point::value || Convex_hull_3::internal::point_inside_convex_polyhedron(P, *origin)
File : /usr/local/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h
Line : 292
Explanation: halfspace_intersection_3: origin not in the polyhedron
Refer to the bug-reporting instructions at http://www.cgal.org/bug_report.html
libc++abi.dylib: terminating with uncaught exception of type CGAL::Assertion_exception: CGAL ERROR: assertion violation!
Expr: boost::is_floating_point::value || Convex_hull_3::internal::point_inside_convex_polyhedron(P, *origin)
File: /usr/local/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h
Line: 292
Explanation: halfspace_intersection_3: origin not in the polyhedron