CGAL difference on Nef_polyhedron produces shape with duplicate vertices

classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

CGAL difference on Nef_polyhedron produces shape with duplicate vertices

crobar
Hi,

I'm working on a computational solid geometry package for the Matlab/Octave programming language. This interface is based on the project pyPolyCSG, which makes use of CGAL to perform it's operations. I have forked and modified this project to make a more generally useful library that is a wrapper for CGAL. In testing, I have found that CGAL appears to make meshes with duplicate nodes under certain conditions which then lead to errors in subsequent operations.

The particular shape that causes the problem for me is constructed from the difference between a sphere and a cube.
Problem shape

This was constructed by subtracting this sphere:



from this cube:



The resulting shape has a duplicate node at the point (273.4741e-003, 0.0000e+000, 418.5832e-003) i.e. on a point on the z-y plane. Neither shape has duplicate nodels before the difference operation. And when I try to the difference this object with another cube, I get:

Error !!!!!!!!!!!!!!!!!!!!!!!
Error !!!!!!!!!!!!!!!!!!!!!!!

from CGAL

However, the problem goes away if I first rotate the sphere by 90 degrees around the y-axis. I suspect this is because one of the lines making up the sphere (which has 20 x 20 segments in these examples) is no longer lying in the same plane as the edge of the cube.

is there a possibility of a bug here?

To reproduce the issue, you can find my modified pyPolyCSG project here:

https://github.com/crobarcro/pyPolyCSG

however, you don't need the whole project, only the "libpolyhcsg" subproject, which does not use python or matlab. It is a cmake project and can be built and installed with:

cmake .
make
sudo make install

The main files of relevance are polyhedron.h/cpp, ane polyhedron_binary_op.h/cpp

The actual issue is reproduced in a test program, which you can find in the "test/cornerrd" subdirectory here:

https://github.com/crobarcro/pyPolyCSG/tree/master/libpolyhcsg/test/cornerrad

This is not a cmake project, but is only a single file, main.cpp.  A code::blocks project is provided.

If I run this I get the following exception thrown when trying to difference the polyhedron with the duplicate vertices with a cube:

Debugger name and version: GNU gdb (Ubuntu 7.7-0ubuntu3.1) 7.7
In __cxa_throw () (/usr/lib/x86_64-linux-gnu/libstdc++.so.6)
#1  0x00007ffff799dea8 in CGAL::Uncertain<bool>::make_certain (this=0x7fffffffd8f0) at /usr/local/include/CGAL/Uncertain.h:125
/usr/local/include/CGAL/Uncertain.h:125:2852:beg:0x7ffff799dea8
At /usr/local/include/CGAL/Uncertain.h:125

The backtrace:

#0  0x00007ffff71a3a30 in __cxa_throw () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#1  0x00007ffff799dea8 in CGAL::Uncertain<bool>::make_certain (this=0x7fffffffd8f0) at /usr/local/include/CGAL/Uncertain.h:125
#2  0x00007ffff7979636 in CGAL::Uncertain<bool>::operator bool (this=0x7fffffffd8f0) at /usr/local/include/CGAL/Uncertain.h:133
#3  0x00007ffff799d927 in CGAL::operator==<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > > (v=..., w=...) at /usr/local/include/CGAL/Cartesian/Vector_3.h:143
#4  0x00007ffff7978e09 in CGAL::CommonKernelFunctors::Equal_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >::operator() (this=0x7fffffffda11, v1=..., v2=...) at /usr/local/include/CGAL/Kernel/function_objects.h:2229
#5  0x00007ffff795aa09 in CGAL::Filtered_predicate<CGAL::CommonKernelFunctors::Equal_3<CGAL::Simple_cartesian<CGAL::Gmpq> >, CGAL::CommonKernelFunctors::Equal_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >, CGAL::Exact_converter<CGAL::Epeck, CGAL::Simple_cartesian<CGAL::Gmpq> >, CGAL::Approx_converter<CGAL::Epeck, CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >, true>::operator()<CGAL::Vector_3<CGAL::Epeck>, CGAL::Vector_3<CGAL::Epeck> > (this=0x7fffffffda10, a1=..., a2=...) at /usr/local/include/CGAL/Filtered_predicate.h:214
#6  0x00007ffff7941af0 in CGAL::operator==<CGAL::Epeck> (p=..., q=...) at /usr/local/include/CGAL/Kernel/global_functions_3.h:784
#7  0x00007ffff792e540 in CGAL::Facet_plane_3::operator()<CGAL::HalfedgeDS_in_place_list_face<CGAL::I_Polyhedron_facet<CGAL::HalfedgeDS_face_base<CGAL::HalfedgeDS_list_types<CGAL::Epeck, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> > > > > (this=0x7fffffffdb60, f=...) at /usr/local/include/CGAL/Nef_3/polyhedron_3_to_nef_3.h:105
#8  0x00007ffff791fd2c in std::transform<CGAL::internal::In_place_list_iterator<CGAL::HalfedgeDS_in_place_list_face<CGAL::I_Polyhedron_facet<CGAL::HalfedgeDS_face_base<CGAL::HalfedgeDS_list_types<CGAL::Epeck, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> > > >, std::allocator<CGAL::HalfedgeDS_in_place_list_face<CGAL::I_Polyhedron_facet<CGAL::HalfedgeDS_face_base<CGAL::HalfedgeDS_list_types<CGAL::Epeck, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> > > > > >, CGAL::Iterator_project<CGAL::internal::In_place_list_iterator<CGAL::HalfedgeDS_in_place_list_face<CGAL::I_Polyhedron_facet<CGAL::HalfedgeDS_face_base<CGAL::HalfedgeDS_list_types<CGAL::Epeck, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> > > >, std::allocator<CGAL::HalfedgeDS_in_place_list_face<CGAL::I_Polyhedron_facet<CGAL::HalfedgeDS_face_base<CGAL::HalfedgeDS_list_types<CGAL::Epeck, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> > > > > >, CGAL::Project_plane<CGAL::HalfedgeDS_in_place_list_face<CGAL::I_Polyhedron_facet<CGAL::HalfedgeDS_face_base<CGAL::HalfedgeDS_list_types<CGAL::Epeck, CGAL::I_Polyhedron_derived_items_3<CGAL::Polyhedron_items_3>, std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> > > > >, int, int, int, int>, CGAL::Facet_plane_3> (__first=..., __last=..., __result=..., __unary_op=...) at /usr/include/c++/4.8/bits/stl_algo.h:4926
#9  0x00007ffff791a764 in CGAL::polyhedron_3_to_nef_3<CGAL::Polyhedron_3<CGAL::Epeck, CGAL::Polyhedron_items_3, CGAL::HalfedgeDS_default, std::allocator<int> >, CGAL::SNC_structure<CGAL::Epeck, CGAL::SNC_indexed_items, bool> > (P=..., S=...) at /usr/local/include/CGAL/Nef_3/polyhedron_3_to_nef_3.h:198
#10 0x00007ffff7917d2c in CGAL::Nef_polyhedron_3<CGAL::Epeck, CGAL::SNC_indexed_items, bool>::Nef_polyhedron_3<CGAL::Epeck, CGAL::Polyhedron_items_3, CGAL::HalfedgeDS_default, std::allocator<int> > (this=0x7fffffffe270, P=...) at /usr/local/include/CGAL/Nef_polyhedron_3.h:566
#11 0x00007ffff79112b8 in polyhcsg::polyhedron_to_cgal (p=...) at /home/rcrozier/src/pyPolyCSG-crobarcro/libpolyhcsg/source/polyhedron_binary_op.cpp:86
#12 0x00007ffff7911bb4 in polyhcsg::polyhedron_difference::operator() (this=0x7fffffffe300, A=..., B=...) at /home/rcrozier/src/pyPolyCSG-crobarcro/libpolyhcsg/source/polyhedron_binary_op.cpp:140
#13 0x00007ffff7acbe29 in polyhcsg::polyhedron::operator- (this=0x7fffffffe340, in=...) at /home/rcrozier/src/pyPolyCSG-crobarcro/libpolyhcsg/source/polyhedron.cpp:758
#14 0x00000000004013e6 in main (argc=1, argv=0x7fffffffe5c8) at /home/rcrozier/src/pyPolyCSG-crobarcro/libpolyhcsg/test/cornerrad/main.cpp:16


If you are interested in the matlab.Octave interface, it can be found here:

https://github.com/crobarcro/mpolycsg

It is quite functional, but not yet documented.
Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

Sebastien Loriot (GeometryFactory)
What kernel are you using internally to do the computation?
I guess you export the result into double/float which result in rounding
of the exact output, which result in points having the same coordinates
in your case.

Sebastien.


On 08/31/2014 02:51 AM, crobar wrote:

> Hi,
>
> I'm working on a computational solid geometry package for the Matlab/Octave
> programming language. This interface is based on the project pyPolyCSG,
> which makes use of CGAL to perform it's operations. I have forked and
> modified this project to make a more generally useful library that is a
> wrapper for CGAL. In testing, I have found that CGAL appears to make meshes
> with duplicate nodes under certain conditions which then lead to errors in
> subsequent operations.
>
> The particular shape that causes the problem for me is constructed from the
> difference between a sphere and a cube.
> <http://cgal-discuss.949826.n4.nabble.com/file/n4659760/fillet.png>
>
> This was constructed by subtracting this sphere:
>
> <http://cgal-discuss.949826.n4.nabble.com/file/n4659760/sphere_unrotated.png>
>
> from this cube:
>
> <http://cgal-discuss.949826.n4.nabble.com/file/n4659760/cube.png>
>
> The resulting shape has a duplicate node at the point (273.4741e-003,
> 0.0000e+000, 418.5832e-003) i.e. on a point on the z-y plane. Neither shape
> has duplicate nodels before the difference operation. And when I try to the
> difference this object with another cube, I get:
>
> Error !!!!!!!!!!!!!!!!!!!!!!!
> Error !!!!!!!!!!!!!!!!!!!!!!!
>
> from CGAL
>
> However, the problem goes away if I first rotate the sphere by 90 degrees
> around the y-axis. I suspect this is because one of the lines making up the
> sphere (which has 20 x 20 segments in these examples) is no longer lying in
> the same plane as the edge of the cube.
>
> is there a possibility of a bug here?
>
> To reproduce the issue, you can find my modified pyPolyCSG project here:
>
> https://github.com/crobarcro/pyPolyCSG
>
> however, you don't need the whole project, only the "libpolyhcsg"
> subproject, which does not use python or matlab. It is a cmake project and
> can be built and installed with:
>
> cmake .
> make
> sudo make install
>
> The main files of relevance are polyhedron.h/cpp, ane
> polyhedron_binary_op.h/cpp
>
> The actual issue is reproduced in a test program, which you can find in the
> "test/cornerrd" subdirectory here:
>
> https://github.com/crobarcro/pyPolyCSG/tree/master/libpolyhcsg/test/cornerrad
>
> This is not a cmake project, but is only a single file, main.cpp.  A
> code::blocks project is provided.
>
> If I run this I get the following exception thrown when trying to difference
> the polyhedron with the duplicate vertices with a cube:
>
> Debugger name and version: GNU gdb (Ubuntu 7.7-0ubuntu3.1) 7.7
> In __cxa_throw () (/usr/lib/x86_64-linux-gnu/libstdc++.so.6)
> #1  0x00007ffff799dea8 in CGAL::Uncertain<bool>::make_certain
> (this=0x7fffffffd8f0) at /usr/local/include/CGAL/Uncertain.h:125
> /usr/local/include/CGAL/Uncertain.h:125:2852:beg:0x7ffff799dea8
> At /usr/local/include/CGAL/Uncertain.h:125
>
> The backtrace:
>
> #0  0x00007ffff71a3a30 in __cxa_throw () from
> /usr/lib/x86_64-linux-gnu/libstdc++.so.6
> #1  0x00007ffff799dea8 in CGAL::Uncertain<bool>::make_certain
> (this=0x7fffffffd8f0) at /usr/local/include/CGAL/Uncertain.h:125
> #2  0x00007ffff7979636 in CGAL::Uncertain<bool>::operator bool
> (this=0x7fffffffd8f0) at /usr/local/include/CGAL/Uncertain.h:133
> #3  0x00007ffff799d927 in
> CGAL::operator==<CGAL::Simple_cartesian&lt;CGAL::Interval_nt&lt;false> > >
> (v=..., w=...) at /usr/local/include/CGAL/Cartesian/Vector_3.h:143
> #4  0x00007ffff7978e09 in
> CGAL::CommonKernelFunctors::Equal_3<CGAL::Simple_cartesian&lt;CGAL::Interval_nt&lt;false>
>>> ::operator() (this=0x7fffffffda11, v1=..., v2=...) at
> /usr/local/include/CGAL/Kernel/function_objects.h:2229
> #5  0x00007ffff795aa09 in
> CGAL::Filtered_predicate<CGAL::CommonKernelFunctors::Equal_3&lt;CGAL::Simple_cartesian&lt;CGAL::Gmpq>
>> ,
> CGAL::CommonKernelFunctors::Equal_3<CGAL::Simple_cartesian&lt;CGAL::Interval_nt&lt;false>
>>> , CGAL::Exact_converter<CGAL::Epeck,
> CGAL::Simple_cartesian&lt;CGAL::Gmpq> >, CGAL::Approx_converter<CGAL::Epeck,
> CGAL::Simple_cartesian&lt;CGAL::Interval_nt&lt;false> > >,
> true>::operator()<CGAL::Vector_3&lt;CGAL::Epeck>,
> CGAL::Vector_3<CGAL::Epeck> > (this=0x7fffffffda10, a1=..., a2=...) at
> /usr/local/include/CGAL/Filtered_predicate.h:214
> #6  0x00007ffff7941af0 in CGAL::operator==<CGAL::Epeck> (p=..., q=...) at
> /usr/local/include/CGAL/Kernel/global_functions_3.h:784
> #7  0x00007ffff792e540 in
> CGAL::Facet_plane_3::operator()<CGAL::HalfedgeDS_in_place_list_face&lt;CGAL::I_Polyhedron_facet&lt;CGAL::HalfedgeDS_face_base&lt;CGAL::HalfedgeDS_list_types&lt;CGAL::Epeck,
> CGAL::I_Polyhedron_derived_items_3&lt;CGAL::Polyhedron_items_3>,
> std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> >
>>>> (this=0x7fffffffdb60, f=...) at
> /usr/local/include/CGAL/Nef_3/polyhedron_3_to_nef_3.h:105
> #8  0x00007ffff791fd2c in
> std::transform<CGAL::internal::In_place_list_iterator&lt;CGAL::HalfedgeDS_in_place_list_face&lt;CGAL::I_Polyhedron_facet&lt;CGAL::HalfedgeDS_face_base&lt;CGAL::HalfedgeDS_list_types&lt;CGAL::Epeck,
> CGAL::I_Polyhedron_derived_items_3&lt;CGAL::Polyhedron_items_3>,
> std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> >
>>> ,
> std::allocator<CGAL::HalfedgeDS_in_place_list_face&lt;CGAL::I_Polyhedron_facet&lt;CGAL::HalfedgeDS_face_base&lt;CGAL::HalfedgeDS_list_types&lt;CGAL::Epeck,
> CGAL::I_Polyhedron_derived_items_3&lt;CGAL::Polyhedron_items_3>,
> std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> >
>>>>> ,
> CGAL::Iterator_project<CGAL::internal::In_place_list_iterator&lt;CGAL::HalfedgeDS_in_place_list_face&lt;CGAL::I_Polyhedron_facet&lt;CGAL::HalfedgeDS_face_base&lt;CGAL::HalfedgeDS_list_types&lt;CGAL::Epeck,
> CGAL::I_Polyhedron_derived_items_3&lt;CGAL::Polyhedron_items_3>,
> std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> >
>>> ,
> std::allocator<CGAL::HalfedgeDS_in_place_list_face&lt;CGAL::I_Polyhedron_facet&lt;CGAL::HalfedgeDS_face_base&lt;CGAL::HalfedgeDS_list_types&lt;CGAL::Epeck,
> CGAL::I_Polyhedron_derived_items_3&lt;CGAL::Polyhedron_items_3>,
> std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> >
>>>>> ,
> CGAL::Project_plane<CGAL::HalfedgeDS_in_place_list_face&lt;CGAL::I_Polyhedron_facet&lt;CGAL::HalfedgeDS_face_base&lt;CGAL::HalfedgeDS_list_types&lt;CGAL::Epeck,
> CGAL::I_Polyhedron_derived_items_3&lt;CGAL::Polyhedron_items_3>,
> std::allocator<int> >, CGAL::Boolean_tag<true>, CGAL::Plane_3<CGAL::Epeck> >
>>>> , int, int, int, int>, CGAL::Facet_plane_3> (__first=..., __last=...,
> __result=..., __unary_op=...) at /usr/include/c++/4.8/bits/stl_algo.h:4926
> #9  0x00007ffff791a764 in
> CGAL::polyhedron_3_to_nef_3<CGAL::Polyhedron_3&lt;CGAL::Epeck,
> CGAL::Polyhedron_items_3, CGAL::HalfedgeDS_default, std::allocator&lt;int>
>> , CGAL::SNC_structure<CGAL::Epeck, CGAL::SNC_indexed_items, bool> > (P=...,
> S=...) at /usr/local/include/CGAL/Nef_3/polyhedron_3_to_nef_3.h:198
> #10 0x00007ffff7917d2c in CGAL::Nef_polyhedron_3<CGAL::Epeck,
> CGAL::SNC_indexed_items, bool>::Nef_polyhedron_3<CGAL::Epeck,
> CGAL::Polyhedron_items_3, CGAL::HalfedgeDS_default, std::allocator&lt;int> >
> (this=0x7fffffffe270, P=...) at
> /usr/local/include/CGAL/Nef_polyhedron_3.h:566
> #11 0x00007ffff79112b8 in polyhcsg::polyhedron_to_cgal (p=...) at
> /home/rcrozier/src/pyPolyCSG-crobarcro/libpolyhcsg/source/polyhedron_binary_op.cpp:86
> #12 0x00007ffff7911bb4 in polyhcsg::polyhedron_difference::operator()
> (this=0x7fffffffe300, A=..., B=...) at
> /home/rcrozier/src/pyPolyCSG-crobarcro/libpolyhcsg/source/polyhedron_binary_op.cpp:140
> #13 0x00007ffff7acbe29 in polyhcsg::polyhedron::operator-
> (this=0x7fffffffe340, in=...) at
> /home/rcrozier/src/pyPolyCSG-crobarcro/libpolyhcsg/source/polyhedron.cpp:758
> #14 0x00000000004013e6 in main (argc=1, argv=0x7fffffffe5c8) at
> /home/rcrozier/src/pyPolyCSG-crobarcro/libpolyhcsg/test/cornerrad/main.cpp:16
>
>
> If you are interested in the matlab.Octave interface, it can be found here:
>
> https://github.com/crobarcro/mpolycsg
>
> It is quite functional, but not yet documented.
>
>
>
> --
> View this message in context: http://cgal-discuss.949826.n4.nabble.com/CGAL-difference-on-Nef-polyhedron-produces-shape-with-duplicate-vertices-tp4659760.html
> Sent from the cgal-discuss mailing list archive at Nabble.com.
>


--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

crobar
Sebastien Loriot (GeometryFactory) wrote
What kernel are you using internally to do the computation?
I guess you export the result into double/float which result in rounding
of the exact output, which result in points having the same coordinates
in your case.

Sebastien.
Here is the bulk of the CGAL wrapper:

#include <CGAL/Polyhedron_items_with_id_3.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include<CGAL/Polyhedron_incremental_builder_3.h>
#include<CGAL/Polyhedron_3.h>
#include<CGAL/Nef_polyhedron_3.h>
#include<CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/Iterator_project.h>
#include <CGAL/function_objects.h>


typedef CGAL::Exact_predicates_exact_constructions_kernel     Kernel;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;
typedef Polyhedron::HalfedgeDS             HalfedgeDS;
typedef CGAL::Nef_polyhedron_3<Kernel>     Nef_polyhedron;

// A modifier creating a triangle with the incremental builder.
template<class HDS>
class polyhedron_builder : public CGAL::Modifier_base<HDS> {
public:
    polyhedron t;

    polyhedron_builder( const polyhedron &p ){
        t = p;
    }
    void operator()( HDS& hds) {
        typedef typename HDS::Vertex   Vertex;
        typedef typename Vertex::Point Point;


        // create a cgal incremental builder
        CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
        B.begin_surface( t.num_vertices(), t.num_faces() );

        // add the polyhedron vertices
        for( int i=0; i<t.num_vertices(); i++ ){
            double x, y, z;
            t.get_vertex( i, x, y, z );
            B.add_vertex( Point( x, y, z ) );
        }

        std::vector<double> coords;
        std::vector<int> faces;
        t.output_store_in_mesh( coords, faces );

        int nverts, tmpi = 0;
        while( tmpi < (int)faces.size() ){
            nverts = faces[tmpi++];
            if( nverts != 3 )
                std::cout << "face has " << nverts << " vertices" << std::endl;
            B.begin_facet();
            for( int i=0; i<nverts; i++ ){
                B.add_vertex_to_facet( faces[tmpi++] );
            }
            B.end_facet();
        }

        // finish up the surface
        B.end_surface();
    }
};

Nef_polyhedron polyhedron_to_cgal( const polyhedron &p ){
    polyhedron tmp = p.triangulate();
    Polyhedron P;
    polyhedron_builder<HalfedgeDS> builder( tmp );
    P.delegate( builder );
    if( P.is_closed() )
        return Nef_polyhedron( P );
    else
        std::cout << "input polyhedron is not closed!" << std::endl;

    return Nef_polyhedron();
}

polyhedron cgal_to_polyhedron( const Nef_polyhedron &NP ){
    Polyhedron P;
    polyhedron ret;

    if( NP.is_simple() ){
        NP.convert_to_polyhedron(P);
        std::vector<double> coords;
        std::vector<int> tris;
        int next_id = 0;
        std::map< Polyhedron::Vertex*, int > vid;
        for( Polyhedron::Vertex_iterator iter=P.vertices_begin(); iter!=P.vertices_end(); iter++ ){
            coords.push_back( CGAL::to_double( (*iter).point().x() ) );
            coords.push_back( CGAL::to_double( (*iter).point().y() ) );
            coords.push_back( CGAL::to_double( (*iter).point().z() ) );
            vid[ &(*iter) ] = next_id++;
        }

        for( Polyhedron::Facet_iterator iter=P.facets_begin(); iter!=P.facets_end(); iter++ ){
            Polyhedron::Halfedge_around_facet_circulator j = iter->facet_begin();
            tris.push_back( CGAL::circulator_size(j) );
            do {
                tris.push_back( std::distance(P.vertices_begin(), j->vertex()) );
            } while ( ++j != iter->facet_begin());
        }

        ret.initialize_load_from_mesh( coords, tris );
    } else {
        std::cout << "resulting polyhedron is not simple!" << std::endl;
    }
    return ret;
}

polyhedron polyhedron_union::operator()( const polyhedron &A, const polyhedron &B ){
    Nef_polyhedron a, b, c;
    try {
        a = polyhedron_to_cgal( A );
        b = polyhedron_to_cgal( B );
        c = (a + b).interior().closure();
        return cgal_to_polyhedron( c );
    } catch( std::exception &e ){
        return A;
    }
}

polyhedron polyhedron_difference::operator()( const polyhedron &A, const polyhedron &B ){
    Nef_polyhedron a, b, c;
    try {
        a = polyhedron_to_cgal( A );
        b = polyhedron_to_cgal( B );
        c = (a - b).interior().closure();
        return cgal_to_polyhedron( c );
    } catch( std::exception &e ){
        return A;
    }
}

polyhedron polyhedron_symmetric_difference::operator()( const polyhedron &A, const polyhedron &B ){
    Nef_polyhedron a, b, c;
    try {
        a = polyhedron_to_cgal( A );
        b = polyhedron_to_cgal( B );
        c = (a ^ b).interior().closure();
        return cgal_to_polyhedron( c );
    } catch( std::exception &e ){
        return A;
    }
}

polyhedron polyhedron_intersection::operator()( const polyhedron &A, const polyhedron &B ){
    Nef_polyhedron a, b, c;
    try {
        a = polyhedron_to_cgal( A );
        b = polyhedron_to_cgal( B );
        c = (a * b).interior().closure();
        return cgal_to_polyhedron( c );
    } catch( std::exception &e ){
        return A;
    }
}



I'm going to level with you that I did not write this code, and am just making use of it. I once looked at using CGAL myself, but it was just too complex and daunting, I could never have justified the time needed to learn it. Luckily someone else has done the hard work for me in this project.

Therefore I don't know all the details  You can see all the CGAL code in this file:

https://github.com/crobarcro/pyPolyCSG/blob/master/libpolyhcsg/source/polyhedron_binary_op.cpp

As an incentive, consider that having such a simple wrapper available, and also available in higher level languages like python and Matlab/Octave might increase the CGAL user base. :-)

Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

crobar
In reply to this post by Sebastien Loriot (GeometryFactory)
A simpler question, is there a way to get CGAL to merge vertices to a given tolerance prior to extracting the data?

Sebastien Loriot (GeometryFactory) wrote
What kernel are you using internally to do the computation?
I guess you export the result into double/float which result in rounding
of the exact output, which result in points having the same coordinates
in your case.

Sebastien.

Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

Sebastien Loriot (GeometryFactory)
On 09/01/2014 06:23 PM, crobar wrote:
> A simpler question, is there a way to get CGAL to merge vertices to a given
> tolerance prior to extracting the data?
>

You can try using the Surface mesh simplification package and remove
short edges:

http://doc.cgal.org/latest/Surface_mesh_simplification/index.html

Sebastien.

>
> Sebastien Loriot (GeometryFactory) wrote
>> What kernel are you using internally to do the computation?
>> I guess you export the result into double/float which result in rounding
>> of the exact output, which result in points having the same coordinates
>> in your case.
>>
>> Sebastien.
>
>
>
>
>
>
> --
> View this message in context: http://cgal-discuss.949826.n4.nabble.com/CGAL-difference-on-Nef-polyhedron-produces-shape-with-duplicate-vertices-tp4659760p4659764.html
> Sent from the cgal-discuss mailing list archive at Nabble.com.
>


--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

crobar
Sebastien Loriot (GeometryFactory) wrote
You can try using the Surface mesh simplification package and remove
short edges:

http://doc.cgal.org/latest/Surface_mesh_simplification/index.html

Sebastien.

I've had a go using this, but I'm not getting the results I'm expecting, or rather I'm failing to prune the short edges. Here is what I tried (largely based on the example shown in http://doc.cgal.org/latest/Surface_mesh_simplification/Surface_mesh_simplification_2edge_collapse_enriched_polyhedron_8cpp-example.html ):


#include <CGAL/Polyhedron_items_with_id_3.h>
#include<CGAL/Nef_polyhedron_3.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include<CGAL/Polyhedron_incremental_builder_3.h>
#include<CGAL/Polyhedron_3.h>
#include<CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/Iterator_project.h>
#include <CGAL/function_objects.h>

// Adaptor for Polyhedron_3
#include <CGAL/Surface_mesh_simplification/HalfedgeGraph_Polyhedron_3.h>
// Simplification function
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
// Visitor base
#include <CGAL/Surface_mesh_simplification/Edge_collapse_visitor_base.h>
// Stop-condition policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_ratio_stop_predicate.h>
// Non-default cost and placement policies
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_and_length.h>
//#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>


typedef CGAL::Exact_predicates_exact_constructions_kernel     Kernel;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;
typedef Polyhedron::HalfedgeDS             HalfedgeDS;
typedef CGAL::Nef_polyhedron_3<Kernel>     Nef_polyhedron;

// the following is used for checking for degenerate meshes
typedef Kernel::Point_3 Point ;
typedef CGAL::Polyhedron_3<Kernel,CGAL::Polyhedron_items_with_id_3> Polyhedron_id;
typedef Polyhedron_id::Halfedge_handle Halfedge_handle ;
typedef Polyhedron_id::Vertex_handle Vertex_handle ;
typedef CGAL::Surface_mesh_simplification::Edge_profile<Polyhedron_id> Profile ;
typedef  boost::graph_traits<Polyhedron_id>::edges_size_type size_type ;

//*******************************************************************************************************************
// -= stopping condition predicate =-
//
// Determines whether the simplification has finished.
// The arguments are (current_cost,vertex,vertex,is_edge,initial_pair_count,current_pair_count,surface) and the result is bool
//
//*******************************************************************************************************************
//
// Stops when the cost is below some tolerance
//
template<class ECM_>
        class length_stop_predicate
{
        public:
            typedef ECM_ ECM ;
            typedef CGAL::Surface_mesh_simplification::Edge_profile<ECM> Profile ;
            typedef typename boost::graph_traits<ECM>::edge_descriptor edge_descriptor ;
            typedef typename boost::graph_traits<ECM>::edges_size_type size_type ;
            typedef typename CGAL::halfedge_graph_traits<ECM>::Point Point ;
            typedef typename CGAL::Kernel_traits<Point>::Kernel Kernel ;
            typedef typename Kernel::FT FT ;
        public :
            length_stop_predicate( double tolerance ) : _tol(tolerance) {}
            bool operator()( FT const& aCurrentCost // aCurrentCost
                    , Profile const& //aEdgeProfile
                    , size_type //aInitialCount
                    , size_type //aCurrentCount
                    ) const
            {
                std::cout << "sqrt aCurrentCost Is " << sqrt(CGAL::to_double(aCurrentCost)) << " _tol is " << _tol << std::endl;
                return ( sqrt(CGAL::to_double(aCurrentCost)) > _tol );
            }
        private:
            double _tol;
};

 // The following is a Visitor that keeps track of the simplification process.
 // In this example the progress is printed real-time and a few statistics are
 // recorded (and printed in the end).
 //
 struct Stats
 {
     Stats()
     : collected(0)
     , processed(0)
     , collapsed(0)
     , non_collapsable(0)
     , cost_uncomputable(0)
     , placement_uncomputable(0)
     {}
     std::size_t collected ;
     std::size_t processed ;
     std::size_t collapsed ;
     std::size_t non_collapsable ;
     std::size_t cost_uncomputable ;
     std::size_t placement_uncomputable ;
 } ;

 struct My_visitor : CGAL::Surface_mesh_simplification::Edge_collapse_visitor_base<Polyhedron_id>
 {
     My_visitor( Stats* s) : stats(s){}
     // Called during the collecting phase for each edge collected.
     void OnCollected( Profile const&, boost::optional<Kernel::FT> const& )
     {
         ++ stats->collected ;
         std::cerr << "Edges collected: " << stats->collected << std::endl << std::flush ;
     }

     // Called during the processing phase for each edge selected.
     // If cost is absent the edge won't be collapsed.
     void OnSelected(Profile const&
             ,boost::optional<Kernel::FT> cost
             ,std::size_t initial
             ,std::size_t current
             )
     {
         ++ stats->processed ;
         if ( !cost )
             ++ stats->cost_uncomputable ;
         if ( current == initial )
             std::cerr << "\n" << std::flush ;
         std::cerr << "\r" << current << std::flush ;
     }

     // Called during the processing phase for each edge being collapsed.
     // If placement is absent the edge is left uncollapsed.
     void OnCollapsing(Profile const&
             ,boost::optional<Point> placement
             )
     {
         if ( !placement )
             ++ stats->placement_uncomputable ;
     }

     // Called for each edge which failed the so called link-condition,
     // that is, which cannot be collapsed because doing so would
     // turn the surface mesh into a non-manifold.
     void OnNonCollapsable( Profile const& )
     {
         ++ stats->non_collapsable;
     }

     // Called AFTER each edge has been collapsed
     void OnCollapsed( Profile const&, Vertex_handle )
     {
         ++ stats->collapsed;
     }
     Stats* stats ;
 } ;


Then, where I convert from CGAL to my own polyhedron construct, I first do the mesh simplification like so:


polyhedron cgal_to_polyhedron( const Nef_polyhedron &NP ){
    Polyhedron_id P;
    polyhedron ret;

    if( NP.is_simple() )
    {
        NP.convert_to_polyhedron(P);

        // ensure there are no degenerate (too small) edges

        // The items in this polyhedron have an "id()" field
        // which the default index maps used in the algorithm
        // need to get the index of a vertex/edge.
        // However, the Polyhedron_3 class doesn't assign any value to
        // this id(), so we must do it here:
        int index = 0 ;
        for( Polyhedron_id::Halfedge_iterator eb = P.halfedges_begin()
            , ee = P.halfedges_end()
            ; eb != ee
            ; ++ eb
           )
        {
            eb->id() = index++;
        }

        index = 0 ;
        for( Polyhedron_id::Vertex_iterator vb = P.vertices_begin()
            , ve = P.vertices_end()
            ; vb != ve
            ; ++ vb
           )
        {
            vb->id() = index++;
        }

        length_stop_predicate<Polyhedron_id> stop (0.001);
        //CGAL::Surface_mesh_simplification::Count_ratio_stop_predicate<Polyhedron_id> stop(0.99);

        Stats stats ;

        My_visitor vis(&stats) ;

        // The index maps are not explicitetly passed as in the previous
        // example because the surface mesh items have a proper id() field.
        // On the other hand, we pass here explicit cost and placement
        // function which differ from the default policies, ommited in
        // the previous example.
        int r = CGAL::Surface_mesh_simplification::edge_collapse (
                                    P
                                    ,stop
                                    ,CGAL::get_cost (CGAL::Surface_mesh_simplification::Edge_length_cost <Polyhedron_id>())
                                    .get_placement(CGAL::Surface_mesh_simplification::LindstromTurk_placement<Polyhedron_id>())
                                    .visitor (vis)
                                                                 );

        std::cout << "\nEdges collected: " << stats.collected
                  << "\nEdges processed: " << stats.processed
                  << "\nEdges collapsed: " << stats.collapsed
                  << std::endl
                  << "\nEdges not collapsed due to topological constriants: " << stats.non_collapsable
                  << "\nEdge not collapsed due to cost computation constriants: " << stats.cost_uncomputable
                  << "\nEdge not collapsed due to placement computation constriants: " << stats.placement_uncomputable
                  << std::endl ;

        std::cout << "\nFinished...\n" << r << " edges removed.\n"
                  << (P.size_of_halfedges()/2) << " final edges.\n" ;


        std::vector<double> coords;
        std::vector<int> tris;
        int next_id = 0;
        std::map< Polyhedron_id::Vertex*, int > vid;
        for( Polyhedron_id::Vertex_iterator iter = P.vertices_begin(); iter != P.vertices_end(); iter++ )
        {
            coords.push_back( CGAL::to_double( (*iter).point().x() ) );
            coords.push_back( CGAL::to_double( (*iter).point().y() ) );
            coords.push_back( CGAL::to_double( (*iter).point().z() ) );
            vid[ &(*iter) ] = next_id++;
        }

        for( Polyhedron_id::Facet_iterator iter = P.facets_begin(); iter != P.facets_end(); iter++ )
        {
            Polyhedron_id::Halfedge_around_facet_circulator j = iter->facet_begin();
            tris.push_back( CGAL::circulator_size(j) );
            do
            {
                tris.push_back( std::distance(P.vertices_begin(), j->vertex()) );
            } while ( ++j != iter->facet_begin());
        }

        ret.initialize_load_from_mesh( coords, tris );
    }
    else
    {
        std::cout << "resulting polyhedron is not simple!" << std::endl;
    }
    return ret;
}




This time I triangulated my polygons before passing them into CGAL. Now at the point where the problem shape is converted I get the following output:


252sqrt aCurrentCost Is 5.10207e-18 _tol is 0.001
252sqrt aCurrentCost Is 8.29476e-18 _tol is 0.001
252sqrt aCurrentCost Is 9.05531e-18 _tol is 0.001
252sqrt aCurrentCost Is 1.0065e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.11633e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.26797e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.36607e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.39187e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.43805e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.47533e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.49587e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.54581e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.59169e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.81504e-17 _tol is 0.001
252sqrt aCurrentCost Is 1.90392e-17 _tol is 0.001
252sqrt aCurrentCost Is 2.21188e-17 _tol is 0.001
252sqrt aCurrentCost Is 2.28058e-17 _tol is 0.001
252sqrt aCurrentCost Is 2.57659e-17 _tol is 0.001
252sqrt aCurrentCost Is 2.59503e-17 _tol is 0.001
252sqrt aCurrentCost Is 2.8387e-17 _tol is 0.001
252sqrt aCurrentCost Is 2.8972e-17 _tol is 0.001
252sqrt aCurrentCost Is 3.00493e-17 _tol is 0.001
252sqrt aCurrentCost Is 3.08919e-17 _tol is 0.001
252sqrt aCurrentCost Is 3.08919e-17 _tol is 0.001
252sqrt aCurrentCost Is 3.16258e-17 _tol is 0.001
252sqrt aCurrentCost Is 3.36404e-17 _tol is 0.001
252sqrt aCurrentCost Is 0.0257483 _tol is 0.001

Edges collected: 252
Edges processed: 27
Edges collapsed: 0

Edges not collapsed due to topological constriants: 5
Edge not collapsed due to cost computation constriants: 0
Edge not collapsed due to placement computation constriants: 0

Finished...
0 edges removed.
252 final edges.


Could you possibly give me some tips here? Why can't CGAL collapse the edges, and what strategy should I take (or what am I doing wrong/misunderstanding).

For info, below is the problem polyhedron again with the offending point marked.

Problem shape with duplicate node location.


Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

crobar
crobar wrote
        int r = CGAL::Surface_mesh_simplification::edge_collapse (
                                    P
                                    ,stop
                                    ,CGAL::get_cost (CGAL::Surface_mesh_simplification::Edge_length_cost <Polyhedron_id>())
                                    .get_placement(CGAL::Surface_mesh_simplification::LindstromTurk_placement<Polyhedron_id>())
                                    .visitor (vis)
                                                                 );

actually this wasn't write, I actually have:

        int r = CGAL::Surface_mesh_simplification::edge_collapse (
                                    P
                                    ,stop
                                    ,CGAL::get_cost (CGAL::Surface_mesh_simplification::Edge_length_cost <Polyhedron_id>())
                                    .get_placement(CGAL::Surface_mesh_simplification::Midpoint_placement<Polyhedron_id>())
                                    .visitor (vis)
                                                                 );


I was mid changing it back as I wrote this after trying out the Lingstom cost.
Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

Sebastien Loriot (GeometryFactory)
In reply to this post by crobar
If I understand correctly, you're left with 5 edges that cannot be
collapse because of topological constraints.
You should have a look at where these are in your data set because
I'm rather surprised you have such edges considering the type of model
you have.

Sebastien.

On 09/04/2014 05:15 PM, crobar wrote:

> Sebastien Loriot (GeometryFactory) wrote
>> You can try using the Surface mesh simplification package and remove
>> short edges:
>>
>> http://doc.cgal.org/latest/Surface_mesh_simplification/index.html
>>
>> Sebastien.
>
>
> I've had a go using this, but I'm not getting the results I'm expecting, or
> rather I'm failing to prune the short edges. Here is what I tried (largely
> based on the example shown in
> http://doc.cgal.org/latest/Surface_mesh_simplification/Surface_mesh_simplification_2edge_collapse_enriched_polyhedron_8cpp-example.html
> ):
>
>
> #include <CGAL/Polyhedron_items_with_id_3.h>
> #include<CGAL/Nef_polyhedron_3.h>
> #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
> #include<CGAL/Polyhedron_incremental_builder_3.h>
> #include<CGAL/Polyhedron_3.h>
> #include<CGAL/IO/Polyhedron_iostream.h>
> #include <CGAL/Iterator_project.h>
> #include <CGAL/function_objects.h>
>
> // Adaptor for Polyhedron_3
> #include <CGAL/Surface_mesh_simplification/HalfedgeGraph_Polyhedron_3.h>
> // Simplification function
> #include <CGAL/Surface_mesh_simplification/edge_collapse.h>
> // Visitor base
> #include <CGAL/Surface_mesh_simplification/Edge_collapse_visitor_base.h>
> // Stop-condition policy
> #include
> <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_ratio_stop_predicate.h>
> // Non-default cost and placement policies
> #include <CGAL/Surface_mesh_simplification/Detail/Common.h>
> #include
> <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_and_length.h>
> //#include
> <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h>
> #include
> <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
>
>
> typedef CGAL::Exact_predicates_exact_constructions_kernel     Kernel;
> typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;
> typedef Polyhedron::HalfedgeDS             HalfedgeDS;
> typedef CGAL::Nef_polyhedron_3<Kernel>     Nef_polyhedron;
>
> // the following is used for checking for degenerate meshes
> typedef Kernel::Point_3 Point ;
> typedef CGAL::Polyhedron_3<Kernel,CGAL::Polyhedron_items_with_id_3>
> Polyhedron_id;
> typedef Polyhedron_id::Halfedge_handle Halfedge_handle ;
> typedef Polyhedron_id::Vertex_handle Vertex_handle ;
> typedef CGAL::Surface_mesh_simplification::Edge_profile<Polyhedron_id>
> Profile ;
> typedef  boost::graph_traits<Polyhedron_id>::edges_size_type size_type ;
>
> //*******************************************************************************************************************
> // -= stopping condition predicate =-
> //
> // Determines whether the simplification has finished.
> // The arguments are
> (current_cost,vertex,vertex,is_edge,initial_pair_count,current_pair_count,surface)
> and the result is bool
> //
> //*******************************************************************************************************************
> //
> // Stops when the cost is below some tolerance
> //
> template<class ECM_>
>          class length_stop_predicate
> {
>          public:
>              typedef ECM_ ECM ;
>              typedef CGAL::Surface_mesh_simplification::Edge_profile<ECM>
> Profile ;
>              typedef typename boost::graph_traits<ECM>::edge_descriptor
> edge_descriptor ;
>              typedef typename boost::graph_traits<ECM>::edges_size_type
> size_type ;
>              typedef typename CGAL::halfedge_graph_traits<ECM>::Point Point ;
>              typedef typename CGAL::Kernel_traits<Point>::Kernel Kernel ;
>              typedef typename Kernel::FT FT ;
>          public :
>              length_stop_predicate( double tolerance ) : _tol(tolerance) {}
>              bool operator()( FT const& aCurrentCost // aCurrentCost
>                      , Profile const& //aEdgeProfile
>                      , size_type //aInitialCount
>                      , size_type //aCurrentCount
>                      ) const
>              {
>                  std::cout << "sqrt aCurrentCost Is " <<
> sqrt(CGAL::to_double(aCurrentCost)) << " _tol is " << _tol << std::endl;
>                  return ( sqrt(CGAL::to_double(aCurrentCost)) > _tol );
>              }
>          private:
>              double _tol;
> };
>
>   // The following is a Visitor that keeps track of the simplification
> process.
>   // In this example the progress is printed real-time and a few statistics
> are
>   // recorded (and printed in the end).
>   //
>   struct Stats
>   {
>       Stats()
>       : collected(0)
>       , processed(0)
>       , collapsed(0)
>       , non_collapsable(0)
>       , cost_uncomputable(0)
>       , placement_uncomputable(0)
>       {}
>       std::size_t collected ;
>       std::size_t processed ;
>       std::size_t collapsed ;
>       std::size_t non_collapsable ;
>       std::size_t cost_uncomputable ;
>       std::size_t placement_uncomputable ;
>   } ;
>
>   struct My_visitor :
> CGAL::Surface_mesh_simplification::Edge_collapse_visitor_base<Polyhedron_id>
>   {
>       My_visitor( Stats* s) : stats(s){}
>       // Called during the collecting phase for each edge collected.
>       void OnCollected( Profile const&, boost::optional<Kernel::FT> const& )
>       {
>           ++ stats->collected ;
>           std::cerr << "Edges collected: " << stats->collected << std::endl
> << std::flush ;
>       }
>
>       // Called during the processing phase for each edge selected.
>       // If cost is absent the edge won't be collapsed.
>       void OnSelected(Profile const&
>               ,boost::optional<Kernel::FT> cost
>               ,std::size_t initial
>               ,std::size_t current
>               )
>       {
>           ++ stats->processed ;
>           if ( !cost )
>               ++ stats->cost_uncomputable ;
>           if ( current == initial )
>               std::cerr << "\n" << std::flush ;
>           std::cerr << "\r" << current << std::flush ;
>       }
>
>       // Called during the processing phase for each edge being collapsed.
>       // If placement is absent the edge is left uncollapsed.
>       void OnCollapsing(Profile const&
>               ,boost::optional<Point> placement
>               )
>       {
>           if ( !placement )
>               ++ stats->placement_uncomputable ;
>       }
>
>       // Called for each edge which failed the so called link-condition,
>       // that is, which cannot be collapsed because doing so would
>       // turn the surface mesh into a non-manifold.
>       void OnNonCollapsable( Profile const& )
>       {
>           ++ stats->non_collapsable;
>       }
>
>       // Called AFTER each edge has been collapsed
>       void OnCollapsed( Profile const&, Vertex_handle )
>       {
>           ++ stats->collapsed;
>       }
>       Stats* stats ;
>   } ;
>
>
> Then, where I convert from CGAL to my own polyhedron construct, I first do
> the mesh simplification like so:
>
>
> polyhedron cgal_to_polyhedron( const Nef_polyhedron &NP ){
>      Polyhedron_id P;
>      polyhedron ret;
>
>      if( NP.is_simple() )
>      {
>          NP.convert_to_polyhedron(P);
>
>          // ensure there are no degenerate (too small) edges
>
>          // The items in this polyhedron have an "id()" field
>          // which the default index maps used in the algorithm
>          // need to get the index of a vertex/edge.
>          // However, the Polyhedron_3 class doesn't assign any value to
>          // this id(), so we must do it here:
>          int index = 0 ;
>          for( Polyhedron_id::Halfedge_iterator eb = P.halfedges_begin()
>              , ee = P.halfedges_end()
>              ; eb != ee
>              ; ++ eb
>             )
>          {
>              eb->id() = index++;
>          }
>
>          index = 0 ;
>          for( Polyhedron_id::Vertex_iterator vb = P.vertices_begin()
>              , ve = P.vertices_end()
>              ; vb != ve
>              ; ++ vb
>             )
>          {
>              vb->id() = index++;
>          }
>
>          length_stop_predicate<Polyhedron_id> stop (0.001);
>
> //CGAL::Surface_mesh_simplification::Count_ratio_stop_predicate<Polyhedron_id>
> stop(0.99);
>
>          Stats stats ;
>
>          My_visitor vis(&stats) ;
>
>          // The index maps are not explicitetly passed as in the previous
>          // example because the surface mesh items have a proper id() field.
>          // On the other hand, we pass here explicit cost and placement
>          // function which differ from the default policies, ommited in
>          // the previous example.
>          int r = CGAL::Surface_mesh_simplification::edge_collapse (
>                                      P
>                                      ,stop
>                                      ,CGAL::get_cost
> (CGAL::Surface_mesh_simplification::Edge_length_cost <Polyhedron_id>())
>
> .get_placement(CGAL::Surface_mesh_simplification::LindstromTurk_placement<Polyhedron_id>())
>                                      .visitor (vis)
>                                                                   );
>
>          std::cout << "\nEdges collected: " << stats.collected
>                    << "\nEdges processed: " << stats.processed
>                    << "\nEdges collapsed: " << stats.collapsed
>                    << std::endl
>                    << "\nEdges not collapsed due to topological constriants:
> " << stats.non_collapsable
>                    << "\nEdge not collapsed due to cost computation
> constriants: " << stats.cost_uncomputable
>                    << "\nEdge not collapsed due to placement computation
> constriants: " << stats.placement_uncomputable
>                    << std::endl ;
>
>          std::cout << "\nFinished...\n" << r << " edges removed.\n"
>                    << (P.size_of_halfedges()/2) << " final edges.\n" ;
>
>
>          std::vector<double> coords;
>          std::vector<int> tris;
>          int next_id = 0;
>          std::map< Polyhedron_id::Vertex*, int > vid;
>          for( Polyhedron_id::Vertex_iterator iter = P.vertices_begin(); iter
> != P.vertices_end(); iter++ )
>          {
>              coords.push_back( CGAL::to_double( (*iter).point().x() ) );
>              coords.push_back( CGAL::to_double( (*iter).point().y() ) );
>              coords.push_back( CGAL::to_double( (*iter).point().z() ) );
>              vid[ &(*iter) ] = next_id++;
>          }
>
>          for( Polyhedron_id::Facet_iterator iter = P.facets_begin(); iter !=
> P.facets_end(); iter++ )
>          {
>              Polyhedron_id::Halfedge_around_facet_circulator j =
> iter->facet_begin();
>              tris.push_back( CGAL::circulator_size(j) );
>              do
>              {
>                  tris.push_back( std::distance(P.vertices_begin(),
> j->vertex()) );
>              } while ( ++j != iter->facet_begin());
>          }
>
>          ret.initialize_load_from_mesh( coords, tris );
>      }
>      else
>      {
>          std::cout << "resulting polyhedron is not simple!" << std::endl;
>      }
>      return ret;
> }
>
>
>
>
> This time I triangulated my polygons before passing them into CGAL. Now at
> the point where the problem shape is converted I get the following output:
>
>
> 252sqrt aCurrentCost Is 5.10207e-18 _tol is 0.001
> 252sqrt aCurrentCost Is 8.29476e-18 _tol is 0.001
> 252sqrt aCurrentCost Is 9.05531e-18 _tol is 0.001
> 252sqrt aCurrentCost Is 1.0065e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.11633e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.26797e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.36607e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.39187e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.43805e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.47533e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.49587e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.54581e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.59169e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.81504e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 1.90392e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 2.21188e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 2.28058e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 2.57659e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 2.59503e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 2.8387e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 2.8972e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 3.00493e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 3.08919e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 3.08919e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 3.16258e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 3.36404e-17 _tol is 0.001
> 252sqrt aCurrentCost Is 0.0257483 _tol is 0.001
>
> Edges collected: 252
> Edges processed: 27
> Edges collapsed: 0
>
> Edges not collapsed due to topological constriants: 5
> Edge not collapsed due to cost computation constriants: 0
> Edge not collapsed due to placement computation constriants: 0
>
> Finished...
> 0 edges removed.
> 252 final edges.
>
>
> Could you possibly give me some tips here? Why can't CGAL collapse the
> edges, and what strategy should I take (or what am I doing
> wrong/misunderstanding).
>
> For info, below is the problem polyhedron again with the offending point
> marked.
>
> <http://cgal-discuss.949826.n4.nabble.com/file/n4659781/cgalprob.png>
>
>
>
>
>
>
> --
> View this message in context: http://cgal-discuss.949826.n4.nabble.com/CGAL-difference-on-Nef-polyhedron-produces-shape-with-duplicate-vertices-tp4659760p4659781.html
> Sent from the cgal-discuss mailing list archive at Nabble.com.
>


--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

crobar
To be honest, it also surprises me. I have attached two .OFF files which I wrote out from CGAL, just prior to performing the boolean operation. These are a cube and a sphere.

I subtract the sphere from the box. When I do this, and then run the resulting shape through the mesh simplification discussed previously I get the following results:

252sqrt aCurrentCost Is 5.10207e-18 _tol is 1e-08


252sqrt aCurrentCost Is 1.0065e-17 _tol is 1e-08


252sqrt aCurrentCost Is 1.47533e-17 _tol is 1e-08


252sqrt aCurrentCost Is 1.49587e-17 _tol is 1e-08


252sqrt aCurrentCost Is 1.68221e-17 _tol is 1e-08


252sqrt aCurrentCost Is 1.73215e-17 _tol is 1e-08


252sqrt aCurrentCost Is 1.86195e-17 _tol is 1e-08


252sqrt aCurrentCost Is 1.90392e-17 _tol is 1e-08


252sqrt aCurrentCost Is 1.96004e-17 _tol is 1e-08


252sqrt aCurrentCost Is 2.11169e-17 _tol is 1e-08


252sqrt aCurrentCost Is 2.28058e-17 _tol is 1e-08


252sqrt aCurrentCost Is 2.39854e-17 _tol is 1e-08


252sqrt aCurrentCost Is 2.59503e-17 _tol is 1e-08


252sqrt aCurrentCost Is 2.8387e-17 _tol is 1e-08


252sqrt aCurrentCost Is 2.8972e-17 _tol is 1e-08


252sqrt aCurrentCost Is 2.99057e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.00493e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.02278e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.08919e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.08919e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.17361e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.22802e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.3361e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.52399e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.54382e-17 _tol is 1e-08


252sqrt aCurrentCost Is 3.55536e-17 _tol is 1e-08


252sqrt aCurrentCost Is 0.0257483 _tol is 1e-08

Edges collected: 252
Edges processed: 27
Edges collapsed: 0

Edges not collapsed due to topological constraints: 2
Edge not collapsed due to cost computation constraints: 0
Edge not collapsed due to placement computation constraints: 0

Finished...
0 edges removed.
252 final edges.


So even this shape has edges that cannot be collapsed. I can import these two OFF files into FreeCAD (based on OpenCASCADE) and perform the boolean operation. Exporting the cube and sphere from FreeCAD results in shapes with the same number of faces and vertices as the originals, implying to me that there is nothing unusual or problematic about these shapes which FreeCAD is repairing somehow.

I wrote out the OFF files using the following function and the union is performed immediately after (just to show everything happens in CGAL, and does not seem to be introduced by my code):

void off_write (Polyhedron P)
{
    // Write polyhedron in Object File Format (OFF).
    CGAL::set_ascii_mode( std::cout);
    std::cout << "OFF" << std::endl
              << P.size_of_vertices()
              << ' '
              << P.size_of_facets()
              << " 0" << std::endl;

    std::copy( P.points_begin(), P.points_end(),
               std::ostream_iterator<Point>( std::cout, "\n"));

    for ( Facet_iterator i = P.facets_begin(); i != P.facets_end(); ++i)
    {
        Halfedge_facet_circulator j = i->facet_begin();
        // Facets in polyhedral surfaces are at least triangles.
        CGAL_assertion( CGAL::circulator_size(j) >= 3);
        std::cout << CGAL::circulator_size(j) << ' ';
        do
        {
            std::cout << ' ' << std::distance(P.vertices_begin(), j->vertex());
        }
        while ( ++j != i->facet_begin());
        std::cout << std::endl;
    }

}

I'm using CGAL 4.4 in case it is relevant.

cgal_box_triagulated_by_cgal.off

sphere_triangulated_by_cgal.off

Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

Sebastien Loriot (GeometryFactory)
Looking at your OFF files, I can see that they are not valid polyhedral
surface (some facet orientations are incorrect). You can fix that by
loading the off in the polyhedron demo and run the plugin
"Orient polygon soup".

Looking at the operation you want to do, you have a face of the cube
that pass really close to a circle used to sample the points on the
sphere.

The small edges/swap triangles you have when going back to double comes
from there I guess. I don't have a magic solution now but we are working
on that.

Sebastien.


On 09/15/2014 02:58 PM, crobar wrote:

> To be honest, it also surprises me. I have attached two .OFF files which I
> wrote out from CGAL, just prior to performing the boolean operation. These
> are a cube and a sphere.
>
> I subtract the sphere from the box. When I do this, and then run the
> resulting shape through the mesh simplification discussed previously I get
> the following results:
>
> 252sqrt aCurrentCost Is 5.10207e-18 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 1.0065e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 1.47533e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 1.49587e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 1.68221e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 1.73215e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 1.86195e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 1.90392e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 1.96004e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 2.11169e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 2.28058e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 2.39854e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 2.59503e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 2.8387e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 2.8972e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 2.99057e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.00493e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.02278e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.08919e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.08919e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.17361e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.22802e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.3361e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.52399e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.54382e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 3.55536e-17 _tol is 1e-08
>
>
> 252sqrt aCurrentCost Is 0.0257483 _tol is 1e-08
>
> Edges collected: 252
> Edges processed: 27
> Edges collapsed: 0
>
> Edges not collapsed due to topological constraints: 2
> Edge not collapsed due to cost computation constraints: 0
> Edge not collapsed due to placement computation constraints: 0
>
> Finished...
> 0 edges removed.
> 252 final edges.
>
>
> So even this shape has edges that cannot be collapsed. I can import these
> two OFF files into FreeCAD (based on OpenCASCADE) and perform the boolean
> operation. Exporting the cube and sphere from FreeCAD results in shapes with
> the same number of faces and vertices as the originals, implying to me that
> there is nothing unusual or problematic about these shapes which FreeCAD is
> repairing somehow.
>
> I wrote out the OFF files using the following function and the union is
> performed immediately after (just to show everything happens in CGAL, and
> does not seem to be introduced by my code):
>
> void off_write (Polyhedron P)
> {
>      // Write polyhedron in Object File Format (OFF).
>      CGAL::set_ascii_mode( std::cout);
>      std::cout << "OFF" << std::endl
>                << P.size_of_vertices()
>                << ' '
>                << P.size_of_facets()
>                << " 0" << std::endl;
>
>      std::copy( P.points_begin(), P.points_end(),
>                 std::ostream_iterator<Point>( std::cout, "\n"));
>
>      for ( Facet_iterator i = P.facets_begin(); i != P.facets_end(); ++i)
>      {
>          Halfedge_facet_circulator j = i->facet_begin();
>          // Facets in polyhedral surfaces are at least triangles.
>          CGAL_assertion( CGAL::circulator_size(j) >= 3);
>          std::cout << CGAL::circulator_size(j) << ' ';
>          do
>          {
>              std::cout << ' ' << std::distance(P.vertices_begin(),
> j->vertex());
>          }
>          while ( ++j != i->facet_begin());
>          std::cout << std::endl;
>      }
>
> }
>
> I'm using CGAL 4.4 in case it is relevant.
>
> cgal_box_triagulated_by_cgal.off
> <http://cgal-discuss.949826.n4.nabble.com/file/n4659841/cgal_box_triagulated_by_cgal.off>
>
> sphere_triangulated_by_cgal.off
> <http://cgal-discuss.949826.n4.nabble.com/file/n4659841/sphere_triangulated_by_cgal.off>
>
>
>
>
>
> --
> View this message in context: http://cgal-discuss.949826.n4.nabble.com/CGAL-difference-on-Nef-polyhedron-produces-shape-with-duplicate-vertices-tp4659760p4659841.html
> Sent from the cgal-discuss mailing list archive at Nabble.com.
>


--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: CGAL difference on Nef_polyhedron produces shape with duplicate vertices

kintel
Sebastien Loriot wrote:
> The small edges/swap triangles you have when going back to double comes
> from there I guess. I don't have a magic solution now but we are working
> on that.

I've been investigating some similar issues I've had (see separate thread from earlier this week), and it looks like the source of my non-simple polygons could indeed be conversion to double. This can cause an otherwise valid polygon to become self-intersecting, thus triggering the infinite recursion in CDT.

You mention working on a "magic solution". Any news on that?
What I had in mind would be to perform some rough cleaning of the Nef polyhedron (e.g. edge collapse of very short edges). I'm not sure how I would implement that using the Nef API though.
Doing naive vertex clustering to merge vertices won't work as that just pushes the problem ahead (and may distort topology).

Ideas are welcome :)

 -Marius