Fix Isosurfacing_3's Image_3 IO: ref to 'grid' in 'values' must be stable

Isosurfacing_3/include/CGAL/Isosurfacing_3/Cartesian_grid_3.h

@@ -89,7 +89,39 @@ private:
 
   Geom_traits m_gt;
 
+private:
+  void compute_spacing()
+  {
+    auto x_coord = m_gt.compute_x_3_object();
+    auto y_coord = m_gt.compute_y_3_object();
+    auto z_coord = m_gt.compute_z_3_object();
+    auto vertex = m_gt.construct_vertex_3_object();
+
+    // calculate grid spacing
+    const Point_3& min_p = vertex(m_bbox, 0);
+    const Point_3& max_p = vertex(m_bbox, 7);
+    const FT x_span = x_coord(max_p) - x_coord(min_p);
+    const FT y_span = y_coord(max_p) - y_coord(min_p);
+    const FT z_span = z_coord(max_p) - z_coord(min_p);
+
+    const FT d_x = x_span / (m_sizes[0] - 1);
+    const FT d_y = y_span / (m_sizes[1] - 1);
+    const FT d_z = z_span / (m_sizes[2] - 1);
+
+    m_spacing = make_array(d_x, d_y, d_z);
+  }
+
 public:
+  /*!
+   * \brief Default constructor
+   */
+  Cartesian_grid_3()
+    : m_bbox{Point_3{0, 0, 0}, Point_3{0, 0, 0}},
+      m_sizes{2, 2, 2},
+      m_spacing{0, 0, 0},
+      m_gt{Geom_traits()}
+  { }
+
   /**
    * \brief creates a %Cartesian grid with `xdim * ydim * zdim` grid vertices.
    *
@@ -112,22 +144,7 @@ public:
       m_sizes{xdim, ydim, zdim},
       m_gt{gt}
   {
-    auto x_coord = m_gt.compute_x_3_object();
-    auto y_coord = m_gt.compute_y_3_object();
-    auto z_coord = m_gt.compute_z_3_object();
-    auto vertex = m_gt.construct_vertex_3_object();
-
-    // calculate grid spacing
-    const Point_3& min_p = vertex(bbox, 0);
-    const Point_3& max_p = vertex(bbox, 7);
-    const FT x_span = x_coord(max_p) - x_coord(min_p);
-    const FT y_span = y_coord(max_p) - y_coord(min_p);
-    const FT z_span = z_coord(max_p) - z_coord(min_p);
-
-    const FT d_x = x_span / (xdim - 1);
-    const FT d_y = y_span / (ydim - 1);
-    const FT d_z = z_span / (zdim - 1);
-    m_spacing = make_array(d_x, d_y, d_z);
+    compute_spacing();
   }
 
   /**
@@ -221,6 +238,15 @@ public:
    */
   const Iso_cuboid_3& bbox() const { return m_bbox; }
 
+  /**
+   * sets the bounding box of the %Cartesian grid and recomputes the spacing.
+   */
+  void set_bbox(const Iso_cuboid_3& bbox)
+  {
+    m_bbox = bbox;
+    compute_spacing();
+  }
+
   /**
    * \return the spacing of the %Cartesian grid, that is a vector whose coordinates are
    *         the grid steps in the `x`, `y`, and `z` directions, respectively
@@ -242,6 +268,16 @@ public:
    */
   std::size_t zdim() const { return m_sizes[2]; }
 
+  /**
+   * sets the number of grid vertices in the `x`, `y`, and `z` directions, respectively,
+   * and recomputes the spacing.
+   */
+  void set_sizes(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim)
+  {
+    m_sizes = {xdim, ydim, zdim};
+    compute_spacing();
+  }
+
 public:
   /**
    * \brief gets the canonical index of a grid cell given its indices.
@@ -251,8 +287,6 @@ public:
                            const std::size_t k) const
   {
     CGAL_precondition(i < m_sizes[0] && j < m_sizes[1] && k < m_sizes[2]);
-
-    // convert (i, j, k) into a linear index, e.g. to access the scalar values / gradient vectors
     return (k * m_sizes[1] + j) * m_sizes[0] + i;
   }
 

Isosurfacing_3/include/CGAL/Isosurfacing_3/IO/Image_3.h

@@ -32,59 +32,49 @@ namespace IO {
  * The dimensions and bounding box are read from the image. The values stored
  * in the image must be of type `Geom_traits::FT` or implicitly convertible to it.
  *
- * \tparam K must be a model of `IsosurfacingTraits_3`
+ * \tparam Grid must be `CGAL::Cartesian_grid_3<GeomTraits>` whose `GeomTraits` is a model of `IsosurfacingTraits_3`
+ * \tparam Values must be `CGAL::Interpolated_discrete_values_3<Grid>`
  *
  * \param image the image providing the data
- * \param k the traits
+ * \param grid the grid
+ * \tparam values the values
  */
-template <typename K>
-std::pair<Cartesian_grid_3<K>,
-          Interpolated_discrete_values_3<Cartesian_grid_3<K> > >
-read_Image_3(const CGAL::Image_3& image,
-             const K& k = K())
+template <typename Grid, typename Values>
+bool read_Image_3(const CGAL::Image_3& image,
+                  Grid& grid,
+                  Values& values)
 {
-  using Grid = Cartesian_grid_3<K>;
-  using Values = Interpolated_discrete_values_3<Grid>;
+  using Geom_traits = typename Grid::Geom_traits;
+  using FT = typename Geom_traits::FT;
+  using Point_3 = typename Geom_traits::Point_3;
+  using Iso_cuboid_3 = typename Geom_traits::Iso_cuboid_3;
 
-  using FT = typename K::FT;
-  using Point_3 = typename K::Point_3;
-  using Iso_cuboid_3 = typename K::Iso_cuboid_3;
-
-  auto point = k.construct_point_3_object();
-  auto iso_cuboid = k.construct_iso_cuboid_3_object();
-
-  Iso_cuboid_3 bbox;
+  auto point = grid.geom_traits().construct_point_3_object();
+  auto iso_cuboid = grid.geom_traits().construct_iso_cuboid_3_object();
 
   // compute bounding box
   const FT max_x = image.tx() + (image.xdim() - 1) * image.vx();
   const FT max_y = image.ty() + (image.ydim() - 1) * image.vy();
   const FT max_z = image.tz() + (image.zdim() - 1) * image.vz();
-  bbox = iso_cuboid(point(image.tx(), image.ty(), image.tz()),
-                    point(max_x, max_y, max_z));
+  Iso_cuboid_3 bbox = iso_cuboid(point(image.tx(), image.ty(), image.tz()),
+                                 point(max_x, max_y, max_z));
 
   // get spacing
   // std::array<FT, 3> spacing = make_array(image.vx(), image.vy(), image.vz());
 
   // get sizes
-  std::array<std::size_t, 3> sizes;
-  sizes[0] = image.xdim();
-  sizes[1] = image.ydim();
-  sizes[2] = image.zdim();
-
-  Grid grid { bbox, sizes[0], sizes[1], sizes[2], k };
-
-  Values values { grid };
+  grid.set_bbox(bbox);
+  grid.set_sizes(image.xdim(), image.ydim(), image.zdim());
 
   // copy values
-  for(std::size_t x=0; x<sizes[0]; ++x)
-    for(std::size_t y=0; y<sizes[1]; ++y)
-      for(std::size_t z=0; z<sizes[2]; ++z)
+  for(std::size_t x=0; x<image.xdim(); ++x)
+    for(std::size_t y=0; y<image.ydim(); ++y)
+      for(std::size_t z=0; z<image.zdim(); ++z)
         values(x, y, z) = image.value(x, y, z);
 
-  return { grid, values };
+  return true;
 }
 
-
 /**
  * \ingroup IS_IO_functions_grp
  *

Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h

@@ -99,6 +99,7 @@ public:
                 const std::size_t j,
                 const std::size_t k) const
   {
+    CGAL_precondition(i < m_grid.xdim() && j < m_grid.ydim() && k < m_grid.zdim());
     return m_values[m_grid.linear_index(i, j, k)];
   }
 

Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h

@@ -224,7 +224,6 @@ private:
   Point_index add_point(const Point_3& p, const Edge_index& e)
   {
     const Point_index i = m_point_counter++;
-    std::cout << "i =" << i << std::endl;
 
 #ifdef CGAL_LINKED_WITH_TBB
     typename Edge_point_map::accessor acc;

Isosurfacing_3/include/CGAL/Isosurfacing_3/interpolation_schemes_3.h

@@ -146,7 +146,7 @@ public:
     std::size_t j = static_cast<std::size_t>((y_coord(p) - y_coord(min_p)) / spacing[1]);
     std::size_t k = static_cast<std::size_t>((z_coord(p) - z_coord(min_p)) / spacing[2]);
 
-    if(i == g.xdim() - 1)
+    if(i == g.xdim() - 1) // dim is the point number
       --i;
     if(j == g.ydim() - 1)
       --j;