CGAL users discussion list

[costantinos <>]

I have a requirement to merge two Polyhedron_3 objects, one is closed
geometry and the other is a surface. There are common vertices and edges
between them. I have written a routine based on
Polyhedron_incremental_builder_3, I am able to successfully combine the
Polyhedron_3 objects (as checked by outputting results into OFF file),
however assertion errors due to shared halfedges unfortunately inhibit the
method. I am not even sure whether the approach I have taken works, I would
be very grateful for any help/advice or a solution to my problem? 

Ideally, I would get to a point where I am able to combine the geometries
by:

CGAL::poly_merge<Polyhedron, HalfedgeDS> modifier(P1,P2);
P3.delegate(modifier);

or,

CGAL::poly_merge<Polyhedron, HalfedgeDS> modifier(P1,P2);
 P2.delegate(modifier);

(I have checked orientation of surface but no effect + please note the code
below has worked well so far for combining geometry but due to the sorting
method there are exceptions!)



int compare_Point_3(const Point_3* p, const Point_3* q)
{
    if ((fabs((*p)[0] - (*q)[0]) < GEOM_TOL) && (fabs((*p)[1] - (*q)[1]) <
GEOM_TOL) && (fabs((*p)[2] - (*q)[2]) < GEOM_TOL))
    {
        return 0;
    } else {
        double delta = ((*p)[0] * (*p)[0] + (*p)[1] * (*p)[1] + (*p)[2] *
(*p)[2]) - ((*q)[0] * (*q)[0] + (*q)[1] * (*q)[1] + (*q)[2] * (*q)[2]);
        if (delta < 0) {
            return -1;
        } else {
            return 1;
        }
    }
}

bool compare_lessthan(const Point_3* p, const Point_3* q)
{
    return compare_Point_3(p, q) == -1;
}



namespace CGAL {
    
    template <class Poly, class HDS>
    class poly_merge : public Modifier_base<HDS> {
    public:
        typedef Poly Polyhedron;
        typedef HDS  Halfedge_data_structure;
        
        poly_merge(const Polyhedron& poly_1, const Polyhedron& poly_2) :
input_1(poly_1), input_2(poly_2) {}

        void operator()(HDS& output);

    protected:
        
        const Poly& input_1;
        const Poly& input_2;
    };
    
    
    template <class Poly, class HDS>
    void poly_merge<Poly,HDS>:: operator()(HDS& output) {
        
        typedef typename Poly::Vertex_const_iterator Vertex_const_iterator;
        typedef typename Poly::Facet_const_iterator Facet_const_iterator;
        typedef typename HDS::Vertex::Point Point;
        typedef typename Poly::Point_3 Point_3;
        typedef typename std::vector<const Point_3*> pPoints_3;
        typename std::vector<const Point_3*>::iterator p_it;
        typedef typename std::map<const Point_3*, int> Point_map;
        
        pPoints_3 vert_input_1;
        for (Vertex_const_iterator vi = input_1.vertices_begin(); vi !=
input_1.vertices_end(); ++vi) vert_input_1.push_back(&(vi->point()));
        
        std::sort(vert_input_1.begin(), vert_input_1.end(),
compare_lessthan);
        
        
        pPoints_3 vert_input_2;
        for (Vertex_const_iterator vi = input_2.vertices_begin(); vi !=
input_2.vertices_end(); ++vi) vert_input_2.push_back(&(vi->point()));
        
        std::sort(vert_input_2.begin(), vert_input_2.end(),
compare_lessthan);

        
        Point_map pm;
        pPoints_3 vert_input_2_diff, input_intersection;
       
        int i = 0, j = 0;
        while (i < vert_input_1.size() && j < vert_input_2.size()) {
            std::cout << "i: " << i << "; j: " << j << std::endl;
            int c = compare_Point_3(vert_input_1[i], vert_input_2[j]);
            std::cout << "c: " << c << std::endl;
            if (c == -1) {
                ++i;
            } else if (c == 1) {
                vert_input_2_diff.push_back(vert_input_2[j]);
                ++j;
            } else if (c == 0) {
                input_intersection.push_back(vert_input_1[i]);
                pm[vert_input_2[j]] = i;
               // std::cout << "Intersection: " << i << std::endl; 
                ++i; ++j;
            }
        }
        
        if (vert_input_2_diff.size() != vert_input_2.size())
vert_input_2_diff.push_back(vert_input_2[vert_input_2.size() - 1]);

        
        /*
        std::cout << "Input 1 verts: " << vert_input_1.size() << "\n";
        for (p_it = vert_input_1.begin(); p_it != vert_input_1.end();
++p_it) std::cout << *(*p_it) << std::endl;
        std::cout << std::endl << std::endl;
        std::cout << "Input 2 verts: " << vert_input_2.size() << "\n";
        for (p_it = vert_input_2.begin(); p_it != vert_input_2.end();
++p_it) std::cout << *(*p_it) << std::endl;
        std::cout << std::endl << std::endl;
        std::cout << "Intersect points: " << input_intersection.size() <<
"\n";
       for (p_it = input_intersection.begin(); p_it !=
input_intersection.end(); ++p_it) std::cout << *(*p_it) << std::endl;
        std::cout << std::endl << std::endl;
        std::cout << "Input 2 diff: " << vert_input_2_diff.size() << "\n";
        for (p_it = vert_input_2_diff.begin(); p_it !=
vert_input_2_diff.end(); ++p_it) std::cout << *(*p_it) << std::endl; 
        std::cout << std::endl << std::endl;
        std::cout << "intersection map\n"
        "Printing point and position from original vector through pointer"
<< std::endl;
        for (auto it = pm.begin(); it != pm.end(); ++it) {
            std::cout << "Position: " << it->second << std::endl;
            const Point_3* p = it->first;
            std::cout << *p << std::endl;
        }
        */

        Polyhedron_incremental_builder_3<HDS> B(output, true);
        B.begin_surface(vert_input_1.size() + vert_input_2_diff.size(),
input_1.size_of_facets() + input_2.size_of_facets());
        //std::cout << "OFF" << std::endl;
        //std::cout << (vert_input_1.size() + vert_input_2_diff.size()) << "
" << (input_1.size_of_facets() + input_2.size_of_facets())  << " 0" <<
std::endl;
        std::cout << std::endl;
        
        for (i = 0; i < vert_input_1.size(); ++i)
        {
            B.add_vertex(Point_3(*(vert_input_1[i])));
            pm[vert_input_1[i]] = i;
            
            //std::cout << *(vert_input_1[i]) << std::endl;
            //std::cout << "Index 1: " << pm[vert_input_1[i]] << std::endl;
        }
        std::cout << std::endl;
        
        for (j = 0; j < vert_input_2_diff.size(); ++j, ++i)
        {
            B.add_vertex(Point_3(*(vert_input_2_diff[j])));
            pm[vert_input_2_diff[j]] = i;
            
            //std::cout << *(vert_input_2_diff[j]) << std::endl;
            //std::cout << "Index 1: " << pm[vert_input_2_diff[j]] <<
std::endl;
        }

        for (Facet_const_iterator fi = input_1.facets_begin(); fi !=
input_1.facets_end(); ++fi) 
        {
            //std::cout << "3 ";
            B.begin_facet();
            typedef typename Poly::Halfedge_around_facet_const_circulator
            Halfedge_around_facet_const_circulator;
            Halfedge_around_facet_const_circulator hc = fi->facet_begin();
            Halfedge_around_facet_const_circulator hc_end = hc;
            CGAL_assertion( hc != NULL);
            do {
                B.add_vertex_to_facet(pm[&(hc->vertex()->point())]);
                //std::cout << "Point around face 1: " <<
hc->vertex()->point() << std::endl;
                //std::cout << pm[&(hc->vertex()->point())] << " ";
                ++hc;
            } while( hc != hc_end);
            B.end_facet();
            
            //std::cout << std::endl;
        }

        for (Facet_const_iterator fi = input_2.facets_begin(); fi !=
input_2.facets_end(); ++fi) 
        {
            //std::cout << "3 ";
            B.begin_facet();
            typedef typename Poly::Halfedge_around_facet_const_circulator
            Halfedge_around_facet_const_circulator;
            Halfedge_around_facet_const_circulator hc = fi->facet_begin();
            Halfedge_around_facet_const_circulator hc_end = hc;
            CGAL_assertion( hc != NULL);
            do {
                B.add_vertex_to_facet(pm[&(hc->vertex()->point())]);
                //std::cout << "Point around face 2: " <<
hc->vertex()->point() << std::endl;
                //std::cout << pm[&(hc->vertex()->point())] << " ";
                ++hc;
            } while( hc != hc_end);
            B.end_facet();
            //std::cout << std::endl;
        }
        
        B.end_surface();
        output.normalize_border();
    }
    
} 



--
View this message in context: 
http://cgal-discuss.949826.n4.nabble.com/Merging-Combining-two-Polyhedron-3-with-shared-vertices-and-edges-tp4658412.html
Sent from the cgal-discuss mailing list archive at Nabble.com.