CGAL users discussion list

[Bernd Gärtner <>]

Dear Andrea,

I need to apologize for the bad documentation, but the Polytope_distance package requires exact computations to work properly. In the code below, you are using

 typedef CGAL::Polytope_distance_d_traits_3<K> Traits;

It is ok if K is in inexact kernel, but you need to provide some exact number types to the traits for the internal computations through an additional template parameter. The example/Polytope_distance_d directory of the CGAL distribution contains a file polytope_distance_d_fast_exact.cpp that shows how to do this. I attach the file for convenience.

Best,
Bernd.


On 2/7/13 6:21 PM, andrea.tagliasacchi wrote:

I have a convex hull in 3D and need to find all the points (in a set) that fall within a distance \epsilon from it, or that are in its interior (distance<0). This CGAL class seems highly related. To approach this I simply keep "P" fixed, and create a lot of dummy sets "Q" to query the distance. I have tried this solution, but it seems that I am having problems, this is a demo code:

#include "CGAL/point_generators_3.h"
#include <CGAL/Polytope_distance_d.h>
#include <CGAL/Polytope_distance_d_traits_3.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel KEPIC;
typedef CGAL::Cartesian<double>                             KCART;
//typedef KEPIC                                               KERNEL;
typedef KCART                                               KERNEL;
typedef CGAL::Creator_uniform_3<double, KERNEL::Point_3>    Creator;

void function(){
    /// Fill in source hull
    typedef CGAL::Polytope_distance_d_traits_3<K> Traits;
    CGAL::Polytope_distance_d<Traits> pd;
    CGAL::Random_points_in_sphere_3<Point_3, Creator> pgen(100);
    for (int i = 0; i < 100 ; i++, ++pgen)
        pd.insert_p(*pgen);
        
    /// Query point sets
    CGAL::Random_points_in_sphere_3<Point_3, Creator> gen(100);
    Scalar off = 10.0;
    for (int i = 0; i < 100 ; i++, ++gen){
        Point_3 p = *gen;
        Point_3 q(p.x()+off, p.y()+off, p.z()+off);
        std::vector<Point_3> qv;
        qv.push_back(q);
        pd.set_q(qv.begin(), qv.end()); // Can I use ptr arith with just 1 element? (&q, (&q)+1)?

        assert( pd.is_valid(true) );        
        
        if( pd.is_zero() )
            qDebug() << "is_zero";
        if( pd.is_finite() )
            qDebug("dist %f", pd.squared_distance());
    }
 }

for which I get a crash during set_q(...) giving the message:

basis-inverse check: failed ( row=2 | col=0 )
basis-inverse check: failed ( row=2 | col=1 )
basis-inverse check: failed ( row=3 | col=1 )
basis-inverse check: failed ( row=2 | col=2 )
basis-inverse check: failed ( row=3 | col=3 )
basis-inverse check: failed ( row=4 | col=4 )
CGAL error: assertion violation!
_expression_ : check_basis_inverse()
File       : ../../starlab-opttrans/external/CGAL-4.0.2/include/CGAL/QP_solver/QP_solver_impl.h
Line       : 1155
Explanation: 
Refer to the bug-reporting instructions at http://www.cgal.org/bug_report.html
libc++abi.dylib: terminate called throwing an exception

Note however that in a previous test the single query point was actually working. For a while I thought CGAL had problems when using a set Q which was completely inside the CHull(P). With my own data (not using generator) I was getting a fail at the is_valid. Overall... quite unstable... Thoughts? Thanks! :) ---
Andrea Tagliasacchi
Ph.D. Candidate
School of Computer Science
Simon Fraser University

---

View this message in context: Distance from a Convex_hull_3 to Point_3
Sent from the cgal-discuss mailing list archive at Nabble.com.

// computes the distance between two cubes in R^3 using double
// as input type and some internal EXACT floating point type;
// the fast type double is also safely used for many of the
// internal computations
#include <iostream>
#include <cassert>
#include <CGAL/Homogeneous.h>
#include <CGAL/Polytope_distance_d.h>
#include <CGAL/Polytope_distance_d_traits_3.h>

#ifdef CGAL_USE_GMP
#include <CGAL/Gmpzf.h>
typedef CGAL::Gmpzf ET;
#else
#include <CGAL/MP_Float.h>
typedef CGAL::MP_Float ET;
#endif

// use an inexact kernel...
typedef CGAL::Homogeneous<double>                         K;
typedef K::Point_3                                        Point;
// ... and the EXACT traits class based on the inexcat kernel
typedef CGAL::Polytope_distance_d_traits_3<K, ET, double> Traits;
typedef CGAL::Polytope_distance_d<Traits>                 Polytope_distance;


int main()
{
  // the cube [0,1]^3
  Point P[8] = { Point(0,0,0), Point(0,0,1), Point(0,1,0), Point(0,1,1),
                 Point(1,0,0), Point(1,0,1), Point(1,1,0), Point(1,1,1)};

  // the cube [2,3]^3
  Point Q[8] = { Point(2,2,2), Point(2,2,3), Point(2,3,2), Point(2,3,3),
                 Point(3,2,2), Point(3,2,3), Point(3,3,2), Point(3,3,3)};

  Polytope_distance pd(P, P+8, Q, Q+8);
  assert (pd.is_valid());

  // get squared distance (2,2,2)-(1,1,1))^2 = 3
  std::cout << "Squared distance: " <<
    CGAL::to_double (pd.squared_distance_numerator()) /
    CGAL::to_double (pd.squared_distance_denominator()) << std::endl;

  // get points that realize the distance
  Polytope_distance::Coordinate_iterator  coord_it;

  std::cout << "p:"; // homogeneous point from first cube, (1,1,1,1)
  for (coord_it = pd.realizing_point_p_coordinates_begin();
       coord_it != pd.realizing_point_p_coordinates_end();
       ++coord_it)
    std::cout << " " << *coord_it;
  std::cout << std::endl;

  std::cout << "q:"; // homogeneous point from second cube, (2,2,2,1)
  for (coord_it = pd.realizing_point_q_coordinates_begin();
       coord_it != pd.realizing_point_q_coordinates_end();
       ++coord_it)
    std::cout << " " << *coord_it;
  std::cout << std::endl;

  return 0;

}