From 3668ca16b219ef0d7904d0e329c45e31f2ac0f41 Mon Sep 17 00:00:00 2001
From: Laurent Rineau <laurent.rineau@cgal.org>
Date: Wed, 17 Aug 2022 17:21:05 +0200
Subject: [PATCH] Fix test_read_write_point_set

The LAS format, even binary, is not an exact representation.
The coordinates are encoded using and offset and a scaling, for each
coordinate.

In `write_LAS_with_properties()`, Simon has hard-coded:
```
  LASheader header;
  header.x_scale_factor = 1e-9 * (bbox.xmax() - bbox.xmin());
  header.y_scale_factor = 1e-9 * (bbox.ymax() - bbox.ymin());
  header.z_scale_factor = 1e-9 * (bbox.zmax() - bbox.zmin());
  header.x_offset = bbox.xmin();
  header.y_offset = bbox.ymin();
  header.z_offset = bbox.zmin();
```
So, the approximate comparison of coordinates, for the I/O test,
should compare using about the same precision.
---
 .../test_read_write_point_set.cpp             | 36 ++++++++++---------
 1 file changed, 19 insertions(+), 17 deletions(-)

diff --git a/Point_set_processing_3/test/Point_set_processing_3/test_read_write_point_set.cpp b/Point_set_processing_3/test/Point_set_processing_3/test_read_write_point_set.cpp
index abcbcf01747..0717221b8eb 100644
--- a/Point_set_processing_3/test/Point_set_processing_3/test_read_write_point_set.cpp
+++ b/Point_set_processing_3/test/Point_set_processing_3/test_read_write_point_set.cpp
@@ -19,24 +19,12 @@ typedef CGAL::Simple_cartesian<double> Kernel;
 typedef Kernel::Point_3 Point_3;
 typedef Kernel::Vector_3 Vector_3;
 
-const double epsilon = 5e-4;
+const double precision = 1e-8; // relates to 1e9 in write_LAS_with_properties
 
 template < typename Type >
-bool approx_equal_nt(const Type &t1, const Type &t2)
+bool approx_equal_nt(const Type &t1, const Type &t2, const double epsilon)
 {
-  if(t1 == t2)
-    return true;
-  if(CGAL::abs(t1 - t2) / (CGAL::max)(CGAL::abs(t1), CGAL::abs(t2)) < std::abs(t1)*epsilon)
-    return true;
-
-  std::cout << " Approximate comparison failed between : " << t1 << "  and  " << t2 << std::endl;
-  std::cout << "abs(t1 - t2) = " <<CGAL::abs(t1 - t2) << "n";
-  std::cout << "abs(t1) = " <<CGAL::abs(t1) << "n";
-  std::cout << "abs(t2) = " <<CGAL::abs(t2) << "n";
-  std::cout << "std::abs(t1)*epsilon = " <<std::abs(t1) * epsilon << "n";
-  std::cout << "CGAL::abs(t1 - t2) / (CGAL::max)(CGAL::abs(t1), CGAL::abs(t2) == "
-            << CGAL::abs(t1 - t2) / (CGAL::max)(CGAL::abs(t1), CGAL::abs(t2)) << std::endl;
-  return false;
+  return (std::abs(t1 - t2) < epsilon);
 }
 
 typedef std::pair<Point_3, Vector_3> PointVectorPair;
@@ -46,12 +34,19 @@ bool ps_are_equal(const CGAL::Point_set_3<Point_3, Vector_3>& ps,
   if(ps.size() != ps2.size())
     return false;
 
+  const auto bbox = bbox_3(ps.points().begin(), ps.points().end());
+  const auto dx = bbox.xmax() - bbox.xmin();
+  const auto dy = bbox.ymax() - bbox.ymin();
+  const auto dz = bbox.zmax() - bbox.zmin();
+
   typedef CGAL::Point_set_3<Point_3, Vector_3>::const_iterator Iterator;
   for(Iterator it1 = ps.begin(), it2 = ps2.begin(); it1!=ps.end() && it2!=ps2.end(); ++it1, ++it2)
   {
     const Point_3& p1 = ps.point(*it1);
     const Point_3& p2 = ps2.point(*it2);
-    if(!(approx_equal_nt(p1.x(), p2.x()) && approx_equal_nt(p1.y(), p2.y()) && approx_equal_nt(p1.z(), p2.z())))
+    if (!(approx_equal_nt(p1.x(), p2.x(), dx * precision) &&
+          approx_equal_nt(p1.y(), p2.y(), dy * precision) &&
+          approx_equal_nt(p1.z(), p2.z(), dz * precision)))
       return false;
   }
 
@@ -65,12 +60,19 @@ bool points_are_equal(const std::vector<Point_3>& ps,
   if(ps.size() != ps2.size())
     return false;
 
+  const auto bbox = bbox_3(ps.begin(), ps.end());
+  const auto dx = bbox.xmax() - bbox.xmin();
+  const auto dy = bbox.ymax() - bbox.ymin();
+  const auto dz = bbox.zmax() - bbox.zmin();
+
   typedef std::vector<Point_3>::const_iterator Iterator;
   for(Iterator it1 = ps.begin(), it2 = ps2.begin(); it1!=ps.end() && it2!=ps2.end(); ++it1, ++it2)
   {
     const Point_3& p1 = *it1;
     const Point_3& p2 = *it2;
-    if(!(approx_equal_nt(p1.x(), p2.x()) && approx_equal_nt(p1.y(), p2.y()) && approx_equal_nt(p1.z(), p2.z())))
+    if (!(approx_equal_nt(p1.x(), p2.x(), dx * precision) &&
+          approx_equal_nt(p1.y(), p2.y(), dy * precision) &&
+          approx_equal_nt(p1.z(), p2.z(), dz * precision)))
       return false;
   }