Merge pull request #4198 from afabri/Distance-robustify-GF

Make 2d segment segment distance more robust for kernel with inexact constructions
This commit is contained in:
Sebastien Loriot 2020-04-02 15:01:46 +02:00 committed by GitHub
commit 560243e72a
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
10 changed files with 273 additions and 30 deletions

View File

@ -526,6 +526,39 @@ namespace CartesianKernelFunctors {
}
};
template <typename K>
class Compare_signed_distance_to_line_2
{
typedef typename K::Point_2 Point_2;
typedef typename K::Line_2 Line_2;
typedef typename K::Equal_2 Equal_2;
typedef typename K::Less_signed_distance_to_line_2 Less_signed_distance_to_line_2;
public:
typedef typename K::Comparison_result result_type;
result_type
operator()(const Point_2& a, const Point_2& b,
const Point_2& c, const Point_2& d) const
{
CGAL_kernel_precondition_code(Equal_2 equal;)
CGAL_kernel_precondition(! equal(a,b));
return cmp_signed_dist_to_lineC2( a.x(), a.y(),
b.x(), b.y(),
c.x(), c.y(),
d.x(), d.y());
}
result_type
operator()(const Line_2& l, const Point_2& p, const Point_2& q) const
{
Less_signed_distance_to_line_2 less = K().less_signed_distance_to_line_2_object();
if (less(l, p, q)) return SMALLER;
if (less(l, q, p)) return LARGER;
return EQUAL;
}
};
template <typename K>
class Compare_squared_radius_3
{
@ -3884,7 +3917,7 @@ namespace CartesianKernelFunctors {
}
};
// TODO ...
template <typename K>
class Less_signed_distance_to_line_2
{

View File

@ -188,7 +188,108 @@ namespace internal {
typename K::FT
squared_distance(const typename K::Segment_2 &seg1,
const typename K::Segment_2 &seg2,
const K& k)
const K& k,
const Cartesian_tag&)
{
typedef typename K::RT RT;
typedef typename K::FT FT;
bool crossing1, crossing2;
RT c1s, c1e, c2s, c2e;
if (seg1.source() == seg1.target())
return internal::squared_distance(seg1.source(), seg2, k);
if (seg2.source() == seg2.target())
return internal::squared_distance(seg2.source(), seg1, k);
Orientation o1s = orientation(seg2.source(), seg2.target(), seg1.source());
Orientation o1e = orientation(seg2.source(), seg2.target(), seg1.target());
if (o1s == RIGHT_TURN) {
crossing1 = (o1e != RIGHT_TURN);
} else {
if (o1e != LEFT_TURN) {
if (o1s == COLLINEAR && o1e == COLLINEAR)
return internal::squared_distance_parallel(seg1, seg2, k);
crossing1 = true;
} else {
crossing1 = (o1s == COLLINEAR);
}
}
Orientation o2s = orientation(seg1.source(), seg1.target(), seg2.source());
Orientation o2e = orientation(seg1.source(), seg1.target(), seg2.target());
if (o2s == RIGHT_TURN) {
crossing2 = (o2e != RIGHT_TURN);
} else {
if (o2e != LEFT_TURN) {
if (o2s == COLLINEAR && o2e == COLLINEAR)
return internal::squared_distance_parallel(seg1, seg2, k);
crossing2 = true;
} else {
crossing2 = (o2s == COLLINEAR);
}
}
if (crossing1) {
if (crossing2)
return (FT)0;
c2s = CGAL::abs(wcross(seg1.source(), seg1.target(), seg2.source(), k));
c2e = CGAL::abs(wcross(seg1.source(), seg1.target(), seg2.target(), k));
Comparison_result dm = compare(c2s,c2e);
if (dm == SMALLER) {
return internal::squared_distance(seg2.source(), seg1, k);
} else {
if (dm == LARGER) {
return internal::squared_distance(seg2.target(), seg1, k);
} else {
// parallel, should not happen (no crossing)
return internal::squared_distance_parallel(seg1, seg2, k);
}
}
} else {
c1s = CGAL::abs(wcross(seg2.source(), seg2.target(), seg1.source(), k));
c1e = CGAL::abs(wcross(seg2.source(), seg2.target(), seg1.target(), k));
Comparison_result dm = compare(c1s,c1e);
if (crossing2) {
if (dm == SMALLER) {
return internal::squared_distance(seg1.source(), seg2, k);
} else {
if (dm == LARGER) {
return internal::squared_distance(seg1.target(), seg2, k);
} else {
// parallel, should not happen (no crossing)
return internal::squared_distance_parallel(seg1, seg2, k);
}
}
} else {
FT min1, min2;
if (dm == EQUAL)
return internal::squared_distance_parallel(seg1, seg2, k);
min1 = (dm == SMALLER) ?
internal::squared_distance(seg1.source(), seg2, k):
internal::squared_distance(seg1.target(), seg2, k);
c2s = CGAL::abs(wcross(seg1.source(), seg1.target(), seg2.source(), k));
c2e = CGAL::abs(wcross(seg1.source(), seg1.target(), seg2.target(), k));
dm = compare(c2s,c2e);
if (dm == EQUAL) // should not happen.
return internal::squared_distance_parallel(seg1, seg2, k);
min2 = (dm == SMALLER) ?
internal::squared_distance(seg2.source(), seg1, k):
internal::squared_distance(seg2.target(), seg1, k);
return (min1 < min2) ? min1 : min2;
}
}
}
template <class K>
typename K::FT
squared_distance(const typename K::Segment_2 &seg1,
const typename K::Segment_2 &seg2,
const K& k,
const Homogeneous_tag&)
{
typedef typename K::RT RT;
typedef typename K::FT FT;
@ -277,6 +378,7 @@ namespace internal {
}
}
template <class K>
inline typename K::RT
_distance_measure_sub(const typename K::RT &startwcross,
@ -690,7 +792,8 @@ template <class K>
inline typename K::FT
squared_distance(const Segment_2<K> &seg1, const Segment_2<K> &seg2)
{
return internal::squared_distance(seg1, seg2, K());
typedef typename K::Kernel_tag Tag;
return internal::squared_distance(seg1, seg2, K(), Tag());
}
template <class K>

View File

@ -0,0 +1,34 @@
#include <iostream>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_rational.h>
#include <CGAL/Simple_cartesian.h>
typedef CGAL::Simple_cartesian<double> SC;
typedef CGAL::Simple_cartesian<CGAL::Exact_rational> EC;
template <typename Kernel>
double fct() {
typedef typename Kernel::Segment_2 Segment_2;
const Segment_2 segi = {
{ -4.0380854964382, -1.9947196614192 },
{ 10.43442091460618, -0.5886833953492263 } };
const Segment_2 segj = {
{ -11.5138934277993, -2.721011070186227 },
{ -8.822747585009402, -2.459560251317805 } };
const auto dist = CGAL::squared_distance(segi, segj);
std::cout << "#dist: " << dist << std::endl;
return CGAL::to_double(dist);
}
int main()
{
auto approx_dist = fct<SC>();
fct<CGAL::Epick>();
auto exact_dist = fct<EC>();
assert(CGAL::abs(approx_dist - exact_dist) < 0.05 * CGAL::abs(exact_dist));
return 0;
}

View File

@ -761,6 +761,52 @@ namespace HomogeneousKernelFunctors {
}
};
template <typename K>
class Compare_signed_distance_to_line_2
{
typedef typename K::Point_2 Point_2;
typedef typename K::Line_2 Line_2;
typedef typename K::Less_signed_distance_to_line_2 Less_signed_distance_to_line_2;
public:
typedef Comparison_result result_type;
result_type
operator()(const Point_2& p, const Point_2& q,
const Point_2& r, const Point_2& s) const
{
typedef typename K::RT RT;
const RT & phx = p.hx();
const RT & phy = p.hy();
const RT & phw = p.hw();
const RT & qhx = q.hx();
const RT & qhy = q.hy();
const RT & qhw = q.hw();
const RT & rhx = r.hx();
const RT & rhy = r.hy();
const RT & rhw = r.hw();
const RT & shx = s.hx();
const RT & shy = s.hy();
const RT & shw = s.hw();
RT scaled_dist_r_minus_scaled_dist_s =
( rhx*shw - shx*rhw ) * (phy*qhw - qhy*phw)
- ( rhy*shw - shy*rhw ) * (phx*qhw - qhx*phw);
return compare(scaled_dist_r_minus_scaled_dist_s, 0);
}
result_type
operator()(const Line_2& l, const Point_2& p, const Point_2& q) const
{
Less_signed_distance_to_line_2 less = K().less_signed_distance_to_line_2_object();
if (less(l, p, q)) return SMALLER;
if (less(l, q, p)) return LARGER;
return EQUAL;
}
};
template <typename K>
class Compare_slope_2
{
@ -4109,26 +4155,7 @@ namespace HomogeneousKernelFunctors {
operator()(const Point_2& p, const Point_2& q,
const Point_2& r, const Point_2& s) const
{
typedef typename K::RT RT;
const RT & phx= p.hx();
const RT & phy= p.hy();
const RT & phw= p.hw();
const RT & qhx= q.hx();
const RT & qhy= q.hy();
const RT & qhw= q.hw();
const RT & rhx= r.hx();
const RT & rhy= r.hy();
const RT & rhw= r.hw();
const RT & shx= s.hx();
const RT & shy= s.hy();
const RT & shw= s.hw();
RT scaled_dist_r_minus_scaled_dist_s =
( rhx*shw - shx*rhw ) * (phy*qhw - qhy*phw)
- ( rhy*shw - shy*rhw ) * (phx*qhw - qhx*phw);
return scaled_dist_r_minus_scaled_dist_s < 0;
return Compare_signed_distance_to_line_2<K>().operator()(p, q, r, s) == SMALLER;
}
result_type

View File

@ -2,6 +2,14 @@ Release History
===============
[Release 5.1] (https://github.com/CGAL/cgal/releases/tag/releases%2FCGAL-5.1)
### 2D and 3D Linear Geometry Kernel
- Add `CompareSignedDistanceToLine_2` in the 2D/3D Kernel concept to compare
the signed distance of two points to a line, or the line passing through two given points.
Corresponding functors in the model (`Compare_signed_distance_to_line_2`) are also added.
Release 5.0
-----------
Release date: June 2020

View File

@ -976,6 +976,40 @@ public:
/// @}
}; /* end Kernel::ComparePowerDistance_3 */
/*!
\ingroup PkgKernel23ConceptsFunctionObjects
\cgalConcept
\cgalRefines `AdaptableFunctor` (with four arguments)
*/
class CompareSignedDistanceToLine_2 {
public:
/// \name Operations
/// A model of this concept must provide:
/// @{
/*!
compares the signed distance of `r` and `s` to the directed line through `p` and `q`.
*/
Comparison_result operator()(const Kernel::Point_2& p,
const Kernel::Point_2& q,
const Kernel::Point_2& r,
const Kernel::Point_2& s);
/*!
compares the signed distance of `r` and `s` to the directed line `l`.
*/
Comparison_result operator()(const Kernel::Line_2& l,
const Kernel::Point_2& r,
const Kernel::Point_2& s);
/// @}
}; /* end Kernel::CompareSignedDistanceToLine_2 */
/*!
\ingroup PkgKernel23ConceptsFunctionObjects
\cgalConcept
@ -8757,7 +8791,7 @@ public:
\ingroup PkgKernel23ConceptsFunctionObjects
\cgalConcept
\cgalRefines `AdaptableFunctor`` (with four arguments)
\cgalRefines `AdaptableFunctor` (with four arguments)
\sa `has_smaller_signed_distance_to_line_grp`

View File

@ -647,6 +647,11 @@ public:
*/
typedef unspecified_type Compare_angle_with_x_axis_2;
/*!
a model of `Kernel::CompareSignedDistanceToLine_2`
*/
typedef unspecified_type Compare_signed_distance_to_line_2;
/*!
a model of `Kernel::CompareSlope_2`
*/

View File

@ -281,6 +281,7 @@
- `Kernel::CompareSlope_3`
- `Kernel::ComparePowerDistance_2`
- `Kernel::ComparePowerDistance_3`
- `Kernel::CompareSignedDistanceToLine_2`
- `Kernel::CompareSquaredDistance_2`
- `Kernel::CompareSquaredDistance_3`
- `Kernel::CompareSquaredRadius_3`

View File

@ -341,9 +341,7 @@ compare_signed_distance_to_line(const typename K::Point_2& p,
const typename K::Point_2& s,
const K& k)
{
if (k.less_signed_distance_to_line_2_object()(p, q, r, s)) return SMALLER;
if (k.less_signed_distance_to_line_2_object()(p, q, s, r)) return LARGER;
return EQUAL;
return k.compare_signed_distance_to_line_2_object()(p, q, r, s);
}
template <class K>
@ -354,9 +352,7 @@ compare_signed_distance_to_line(const typename K::Line_2& l,
const typename K::Point_2& q,
const K& k)
{
if (k.less_signed_distance_to_line_2_object()(l, p, q)) return SMALLER;
if (k.less_signed_distance_to_line_2_object()(l, q, p)) return LARGER;
return EQUAL;
return k.compare_signed_distance_to_line_2_object()(l, p, q);
}
template < class K >

View File

@ -116,6 +116,8 @@ CGAL_Kernel_pred_RT(Compare_power_distance_2,
compare_power_distance_2_object)
CGAL_Kernel_pred_RT(Compare_power_distance_3,
compare_power_distance_3_object)
CGAL_Kernel_pred(Compare_signed_distance_to_line_2,
compare_signed_distance_to_line_2_object)
CGAL_Kernel_pred(Compare_slope_2,
compare_slope_2_object)
CGAL_Kernel_pred(Compare_slope_3,