CGAL users discussion list

[Panagiotis Cheilaris <>]

On Fri, May 03, 2013 at 11:48:19AM -0700, Pablo Miranda Carranza wrote:
> After some work I have found a better example of the bug I am chasing in
> getting the circumcenters of the segment delaunay (I attach another sample
> program):
> 
> main.cpp <http://cgal-discuss.949826.n4.nabble.com/file/n4657350/main.cpp>  
>  
> 
> I have tried a number of different kernels and number types, all of them
> seem to fail (see sample code). In the output text result of running the
> included program one can see that the centre for "circle 14" is:  (
> 1.61455e+17, 11.9295 ). 
> 
> I attach snapshots of the program I am using to draw the results, obviously
> the program is just to large to include here, but the data is the same. It
> only draws the interior part of the segment voronoi (the medial axis), but
> the error is clearly visible in the largely off position of one of the
> circumcenters.
> 
> <http://cgal-discuss.949826.n4.nabble.com/file/n4657350/Screen_Shot_2013-05-03_at_8.22.26_PM.jpg>
>  
> 
> Do you have any suggestions that may work (any kernels, traits, number
> types...)?  
> Thank you for your time!

I am not an expert in number types, but CORE seems to work fine
for the exact number type, instead of CGAL::Quotient<CGAL::MP_Float>:

$ diff main.cpp maincore.cpp
37c37,39
< typedef CGAL::Quotient<CGAL::MP_Float>     ENT;
---
> #include <CGAL/CORE_Expr.h>
>
> typedef CORE::Expr     ENT;

Remember to add support for CORE, e.g. with:

$ cgal_create_CMakeLists -c Core

Here is the correct circumcenter:

$ ./maincore | grep "circle 14"
circle 14:   -0.104718 11.9295

Best,
Panos

Here is maincore.cpp:

#include <iostream>

//OTHER TYPES THAT ALSO FAIL:
/*
#include <CGAL/MP_Float.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Segment_Delaunay_graph_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_2.h>
 
typedef CGAL::Quotient<CGAL::MP_Float>     ENT;
typedef CGAL::Simple_cartesian<ENT >       CK;
typedef 
CGAL::Segment_Delaunay_graph_traits_without_intersections_2<CK,CGAL::Integral_domain_without_division_tag>
  Gt;
typedef CGAL::Segment_Delaunay_graph_2<Gt>             SDG_2;
*/

/*
#include <CGAL/Filtered_kernel.h>
#include <CGAL/Quotient.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Segment_Delaunay_graph_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_2.h>
 
typedef CGAL::Lazy_exact_nt<CGAL::Quotient<CGAL::MP_Float> >        LENT;
typedef CGAL::Simple_cartesian<LENT>                                 CK;
typedef 
CGAL::Segment_Delaunay_graph_traits_without_intersections_2<CK,CGAL::Integral_domain_without_division_tag>
  Gt;
typedef CGAL::Segment_Delaunay_graph_2<Gt>             SDG_2;
 */


#include <CGAL/Quotient.h>
#include <CGAL/MP_Float.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Segment_Delaunay_graph_filtered_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_2.h>

#include <CGAL/CORE_Expr.h>

typedef CORE::Expr     ENT;
typedef CGAL::Simple_cartesian<double>     CK;
typedef CGAL::Simple_cartesian<ENT>        EK;
typedef 
CGAL::Segment_Delaunay_graph_filtered_traits_2<CK,CGAL::Field_with_sqrt_tag, 
EK, CGAL::Field_tag>  Gt;
typedef CGAL::Segment_Delaunay_graph_2<Gt>             SDG_2;


typedef SDG_2::Point_2                                Point2;

int main(int argc, char *argv[]){
        
    //THESE ARE FORE TEST  OF CIRCUMCENTER
    const int testSize=10;
        Point2 ptData[testSize]={
    Point2(-0.7204357120380354 ,13.89277201765595),
    Point2(-0.7204357120380394 ,12.40288551592062),
    Point2(-3.471103184719311 ,12.40288551592062),
    Point2(-3.471103184719311 ,11.4560803911127),
    Point2(-0.7204357120380394 ,11.4560803911127),
    Point2(-0.7204357120380394 ,11.1683175492763),
    Point2(-1.27464817265971 ,11.1683175492763),
    Point2(-1.27464817265971 ,8.270305214503029),
    Point2(0.6719538274731851 ,8.270305214503029),
    Point2(0.6719538274731853 ,13.89277201765595),
    };
    
        SDG_2   sdg;
        
        //create segment sites for the closed polygon with the points as 
vertices
        std::cout<< "Calculating Segment Delaunay:"<<std::endl;
    
        for (int i = 0; i < testSize; ++i) {
                std::cout<< "processing "<< i <<" of "<<testSize << std::endl;
                int next = (i + 1) % testSize;
                
                const SDG_2::Site_2 
&site=SDG_2::Site_2::construct_site_2(ptData[i], ptData[next]);
                sdg.insert(site);
        }
    
    std::cout<<std::endl;
    std::cout<< "Testing Circumcenter (Primal):"<<std::endl;
    
    int circNum=0;
    for(auto fit=sdg.finite_faces_begin(); fit!= sdg.finite_faces_end(); 
++fit){
        Point2  center = sdg.primal(fit);
        std::cout<<"circle "<< circNum++ <<":   "<<center<<std::endl;
    }
    
    std::cout<<std::endl;
    std::cout<< "Done " << std::endl;
        
}