better walk grid routine. V hard to validate exactly

Author: thomas-lowe

raylib/raycuboid.cpp

@@ -12,12 +12,12 @@ Cuboid::Cuboid(const Eigen::Vector3d &min_bound, const Eigen::Vector3d &max_boun
   max_bound_ = max_bound;
 }
 
-bool Cuboid::clipRay(Eigen::Vector3d &start, Eigen::Vector3d &end) const
+bool Cuboid::clipRay(Eigen::Vector3d &start, Eigen::Vector3d &end, double eps) const
 {
   double max_near_d = 0;
   double min_far_d = 1.0;
   Eigen::Vector3d centre = (min_bound_ + max_bound_) / 2.0;
-  Eigen::Vector3d extent = (max_bound_ - min_bound_) / 2.0;
+  Eigen::Vector3d extent = (max_bound_ - min_bound_) / 2.0 - Eigen::Vector3d(eps, eps, eps);
   Eigen::Vector3d to_centre = centre - start;
   Eigen::Vector3d dir = end - start;
   for (int ax = 0; ax < 3; ax++)

raylib/raycuboid.h

@@ -37,7 +37,7 @@ class RAYLIB_EXPORT Cuboid
   }
 
   /// clip ray to cuboid. Return false if no ray left.
-  bool clipRay(Eigen::Vector3d &start, Eigen::Vector3d &end) const;
+  bool clipRay(Eigen::Vector3d &start, Eigen::Vector3d &end, double eps = 0.0) const;
 
   Eigen::Vector3d min_bound_, max_bound_;
 };

raylib/raydecimation.cpp

@@ -182,11 +182,12 @@ bool decimateRaysSpatial(const std::string &file_stub, double vox_width)
     subsampler.subsample.clear();
     for (int i = 0; i < (int)ends.size(); i++)
     {
+      subsampler.index = i;
       #define END_FIRST // Testing on building.ply it finds more rays, so deemed to be more successful at filling space
       #if defined END_FIRST
-      walkGrid(ends[i] / width, starts[i] / width, subsampler, i);
+      walkGrid(ends[i] / width, starts[i] / width, subsampler);
       #else
-      walkGrid(starts[i] / width, ends[i] / width, subsampler, i);
+      walkGrid(starts[i] / width, ends[i] / width, subsampler);
       #endif 
     }
     chunk.resize(subsampler.subsample.size());

raylib/raydecimation.h

@@ -40,58 +40,19 @@ bool RAYLIB_EXPORT decimateAngular(const std::string &file_stub, double radius_p
 
 struct Subsampler
 {
-  inline bool operator()(const Eigen::Vector3i &p, int i)
+  inline bool operator()(const Eigen::Vector3i &p, const Eigen::Vector3i &/*target*/, double /*in_length*/, double /*out_length*/, double /*max_length*/)
   {
     if (voxel_set.insert(p).second)
     {
-      subsample.push_back(i);
+      subsample.push_back(index);
       return true;
     }
     return false;
   }
   std::vector<int64_t> subsample;
   std::set<Eigen::Vector3i, ray::Vector3iLess> voxel_set;
+  int index;
 };
-
-inline int sign(double x)
-{
-  return (x > 0.0) - (x < 0.0);
-}
-
-// for similar appraoch see: https://github.com/StrandedKitty/tiles-intersect/blob/master/src/index.js
-template<class T> 
-void walkGrid(const Eigen::Vector3d &start, const Eigen::Vector3d &end, T &object, int i)
-{
-  Eigen::Vector3i p = Eigen::Vector3d(std::floor(start[0]), std::floor(start[1]), std::floor(start[2])).cast<int>();
-  
-  if (object(p, i))
-  {
-    return; // only adding to one cell
-  }
-  
-  Eigen::Vector3i step(sign(end[0] - start[0]), sign(end[1] - start[1]), sign(end[2] - start[2]));
-  Eigen::Vector3d tmax, tdelta;
-  for (int j = 0; j<3; j++)
-  {
-    step[j] = sign(end[j] - start[j]);
-    double to = std::abs(start[j] - p[j] - (double)std::max(0, step[j]));        
-    double dir = std::max(std::numeric_limits<double>::min(), std::abs(start[j] - end[j]));
-    tmax[j] = to / dir;
-    tdelta[j] = 1.0 / dir;
-  }
-  
-  Eigen::Vector3i target = Eigen::Vector3d(std::floor(end[0]), std::floor(end[1]), std::floor(end[2])).cast<int>();
-  while (p != target) 
-  {
-    int ax = tmax[0] < tmax[1] && tmax[0] < tmax[2] ? 0 : (tmax[1] < tmax[2] ? 1 : 2);
-    p[ax] += step[ax];
-    tmax[ax] += tdelta[ax];
-    if (object(p, i))
-    {
-      break; // only adding to one cell
-    }          
-  }     
-}
 }  // namespace ray
 
 #endif  // RAYLIB_RAYDECIMATION_H

raylib/rayrenderer.cpp

@@ -279,11 +279,16 @@ void DensityGrid::calculateDensities(const std::string &file_name)
     {
       Eigen::Vector3d start = starts[i];
       Eigen::Vector3d end = ends[i];
-      if (!bounds_.clipRay(start, end))
+      if (!bounds_.clipRay(start, end, 1e-10))
       {
         continue; // ray is outside of bounds
       }
 
+      #define NEW_METHOD
+      #if defined NEW_METHOD
+      bounded_ = colours[i].alpha > 0;
+      walkGrid((start - bounds_.min_bound_) / voxel_width_, (end - bounds_.min_bound_) / voxel_width_, *this);
+      #else
       // now walk the voxels
       Eigen::Vector3d dir = end - start;
       const Eigen::Vector3d source = (start - bounds_.min_bound_) / voxel_width_;
@@ -352,6 +357,7 @@ void DensityGrid::calculateDensities(const std::string &file_name)
           voxels_[index].addMissRay(static_cast<float>(minL * voxel_width_));
         }
       } while (depth <= maxDist);
+      #endif
     }
   };
   Cloud::read(file_name, calculate);

raylib/rayrenderer.h

@@ -107,12 +107,15 @@ struct RAYLIB_EXPORT DensityGrid
   inline const std::vector<Voxel> &voxels() const { return voxels_; }
   inline Eigen::Vector3i dimensions(){ return voxel_dims_; }
   inline Cuboid bounds(){ return bounds_; }
-
+  inline double voxelWidth() const { return voxel_width_; }
+  // used in walking grid only
+  inline bool operator()(const Eigen::Vector3i &p, const Eigen::Vector3i &target, double in_length, double out_length, double max_length);
 private:
   Cuboid bounds_;
   std::vector<Voxel> voxels_;
   double voxel_width_;
   Eigen::Vector3i voxel_dims_;
+  bool bounded_;
 };
 
 // inline functions
@@ -168,6 +171,23 @@ int DensityGrid::getIndexFromPos(const Eigen::Vector3d &pos) const
   Eigen::Vector3d gridspace = (pos - bounds_.min_bound_) / voxel_width_;
   return getIndex(gridspace.cast<int>());
 }
-
+inline bool DensityGrid::operator()(const Eigen::Vector3i &p, const Eigen::Vector3i &target, double in_length, double out_length, double max_length)
+{
+  int index = getIndex(p);
+  if (p == target && bounded_)
+  {
+    double length_in_voxel = std::min(out_length, max_length) - in_length;
+    if (length_in_voxel > 1.733 || length_in_voxel < -1e-7)
+      std::cout << "what's going on?" << out_length << ", " << max_length << ", " << in_length << std::endl;
+    voxels_[index].addHitRay(static_cast<float>(length_in_voxel * voxel_width_));
+  }
+  else
+  {
+    if ((out_length - in_length) > 1.733 || (out_length - in_length) < -1e-7)
+      std::cout << "what on earth's going on?" << (out_length - in_length) << std::endl;
+    voxels_[index].addMissRay(static_cast<float>((out_length - in_length) * voxel_width_));
+  }
+  return false;
+}
 }  // namespace ray
 #endif  // RAYLIB_RAYRENDERER_H

raylib/rayutils.h

@@ -361,6 +361,49 @@ inline std::ostream &logDuration(std::ostream &out, const D &duration)
   }
   return out;
 }
+
+inline int sign(double x)
+{
+  return (x > 0.0) - (x < 0.0);
+}
+
+// for similar appraoch see: https://github.com/StrandedKitty/tiles-intersect/blob/master/src/index.js
+template<class T> 
+void walkGrid(const Eigen::Vector3d &start, const Eigen::Vector3d &end, T &object)
+{
+  Eigen::Vector3d direction = end - start;
+  double max_length = direction.norm();
+  Eigen::Vector3i p = Eigen::Vector3d(std::floor(start[0]), std::floor(start[1]), std::floor(start[2])).cast<int>();
+  const Eigen::Vector3i target = Eigen::Vector3d(std::floor(end[0]), std::floor(end[1]), std::floor(end[2])).cast<int>();
+  
+  const Eigen::Vector3i step(sign(direction[0]), sign(direction[1]), sign(direction[2]));
+  direction /= max_length;
+  Eigen::Vector3d lengths, length_delta;
+  for (int j = 0; j<3; j++)
+  {
+    const double to = std::abs(start[j] - p[j] - (double)std::max(0, step[j]));        
+    const double dir = std::max(std::numeric_limits<double>::epsilon(), std::abs(direction[j]));
+    lengths[j] = to / dir;
+    length_delta[j] = 1.0 / dir;
+  }
+  int ax = lengths[0] < lengths[1] && lengths[0] < lengths[2] ? 0 : (lengths[1] < lengths[2] ? 1 : 2);
+  if (object(p, target, 0.0, lengths[ax], max_length))
+  {      
+    return; // only adding to one cell
+  }
+  
+  while (p != target) 
+  {
+    p[ax] += step[ax];
+    const double in_length = lengths[ax];
+    lengths[ax] += length_delta[ax];
+    ax = lengths[0] < lengths[1] && lengths[0] < lengths[2] ? 0 : (lengths[1] < lengths[2] ? 1 : 2);
+    if (object(p, target, in_length, lengths[ax], max_length))
+    {
+      break; // only adding to one cell
+    }          
+  }     
+}
 }  // namespace ray
 
 #endif  // RAYLIB_RAYUTILS_H