I've implemented some code that uses filtered predicate of cgal. I am defining my own filtered predicate to use Exact_predicates_inexact_constructions_kernel at first and then Exact_predicates_exact_constructions_kernel at second if the first call fails. I thougt that if there occurs some ambiguity with the computations inside filtered predicate, then it understands and turns to use exact constructions and reaches the exact results in the end. However, I saw this is not the case. Inside my filtered predicate, I am computing a point on a plane and checking the orientation of this point with respect to this plane. My expectation was that it would return CGAL::COPLANAR. However, it returns CGAL::NEGATIVE. If there occurred some ambiguity with inexact constructions, then does not filtered kernel understand it? How does filtered kernel work? Does not it throw an exception and repeat the operations with the second template argument? Or are the data types that I use a wrong choice?
Edit: By the way, when I put the point into the plane equation it prints -1.05501e-10. Should I change EPSILON value used by CGAL. If yes, how should I do it?
You can find my code below:
template <typename K>
struct ClosestPlaneToPoint_1D {
typedef typename K::RT RT;
typedef typename K::Point_3 Point_3;
typedef typename K::Plane_3 Plane_3;
typedef typename K::Line_3 Line_3;
typedef typename CGAL::Comparison_result result_type;
result_type operator()( Point_3 plane1_p1, Point_3 plane1_p2, Point_3 plane1_p3,
Point_3 plane2_p1, Point_3 plane2_p2, Point_3 plane2_p3,
Point_3 candidate_plane_p1, Point_3 candidate_plane_p2, Point_3 candidate_plane_p3,
Point_3 front_plane_p1, Point_3 front_plane_p2, Point_3 front_plane_p3,
Point_3 back_plane_p1, Point_3 back_plane_p2, Point_3 back_plane_p3) const
{
Plane_3 plane1(plane1_p1, plane1_p2, plane1_p3);
Plane_3 plane2(plane2_p1, plane2_p2, plane2_p3);
Plane_3 candidate_plane(candidate_plane_p1, candidate_plane_p2, candidate_plane_p3);
if (CGAL::parallel(candidate_plane, plane1) || CGAL::parallel(candidate_plane, plane2))
return CGAL::SMALLER; // not valid
auto auto_line = CGAL::intersection(plane1, plane2);
Line_3* line = boost::get<Line_3>(&*auto_line);
auto auto_candidate_point = CGAL::intersection(candidate_plane, *line);
Point_3* candidate_point = boost::get<Point_3>(&*auto_candidate_point);
if (candidate_point == NULL)
return CGAL::SMALLER;
CGAL::Sign sign = CGAL::orientation(back_plane_p1, back_plane_p2, back_plane_p3, *candidate_point);
cout << "SIGN: " << sign << endl;
Plane_3 back_plane2(back_plane_p1, back_plane_p2, back_plane_p3);
cout << "PLACEMENT: " << back_plane2.a() * candidate_point->x() + back_plane2.b() * candidate_point->y() + back_plane2.c() * candidate_point->z() + back_plane2.d() << endl;
if (sign == CGAL::POSITIVE)
return CGAL::SMALLER;
if (sign == CGAL::COPLANAR) {
Plane_3 back_plane(back_plane_p1, back_plane_p2, back_plane_p3);
if (CGAL::parallel(back_plane, candidate_plane) &&
CGAL::scalar_product(back_plane.orthogonal_vector(), candidate_plane.orthogonal_vector()) > 0)
return CGAL::SMALLER;
}
Plane_3 front_plane(front_plane_p1, front_plane_p2, front_plane_p3);
auto auto_front_point = CGAL::intersection(front_plane, *line);
Point_3* front_point = boost::get<Point_3>(&*auto_front_point);
Plane_3 back_plane(back_plane_p1, back_plane_p2, back_plane_p3);
auto auto_back_point = CGAL::intersection(back_plane, *line);
Point_3* back_point = boost::get<Point_3>(&*auto_back_point);
RT distance1 = CGAL::squared_distance(*back_point, *front_point);
RT distance2 = CGAL::squared_distance(*back_point, *candidate_point);
return CGAL::compare(distance1, distance2);
}
};
typedef CGAL::Filtered_predicate<ClosestPlaneToPoint_1D<ExK>,
ClosestPlaneToPoint_1D<InK>, C2ExK, C2InK> CP_1D;
The above code computes the intersection point of plane1, plane2 and candidate_plane whose points are given in the arguments. Next the point is checked whether it is in front of the back_plane or not. The input that I give produces an intersection point locating on back_plane. But when it checks the point position with respect to back_plane, it says that it is on the positive side of the back_plane.
Something like the following should be working:
of course it will crash at runtime but it's because all points are default constructed.