last before removing old method

Author: thomas-lowe

raylib/CMakeLists.txt

@@ -108,7 +108,7 @@ set(SOURCES
 )
 
 get_target_property(SIMPLE_FFT_INCLUDE_DIRS simple_fft INTERFACE_INCLUDE_DIRECTORIES)
-add_compile_options("-fPIC")
+add_compile_options("-fPIC" -O0)
 
 if(WITH_QHULL)
 set(QHULL_LIBS

raylib/raydecimation.cpp

@@ -164,6 +164,10 @@ bool decimateSpatioTemporal(const std::string &file_stub, double vox_width, int
   return true;
 }
 
+int sign(double x)
+{
+  return (x > 0.0) - (x < 0.0);
+}
 
 bool decimateRaysSpatial(const std::string &file_stub, double vox_width)
 {
@@ -193,7 +197,11 @@ bool decimateRaysSpatial(const std::string &file_stub, double vox_width)
       const Eigen::Vector3d source = starts[i] / width;
       const Eigen::Vector3d target = ends[i] / width;
       #endif
+
+      #define OLD_METHOD
+      #if defined OLD_METHOD
       const double length = dir.norm();
+      dir /= length; 
       for (int a = 0; a<3; a++)
       {
         if (dir[a] == 0.0)
@@ -236,16 +244,76 @@ bool decimateRaysSpatial(const std::string &file_stub, double vox_width)
                          (std::round(p[2] + offsets[2]) - p[2]) / dir[2] };
         int axis = (ls[0] < ls[1] && ls[0] < ls[2]) ? 0 : (ls[1] < ls[2] ? 1 : 2);
         inds[axis] += adds[axis];
-        double minL = ls[axis] * length;
-        depth += minL + eps;
-        p = source + dir * (depth / length);
+        depth += ls[axis] + eps;
+        p = source + dir * depth;
 
         if (voxel_set.insert(inds).second)
         {
           subsample.push_back(i);
           break; // only adding to one cell
         }
       } while (depth <= maxDist);
+      #else
+        Eigen::Vector3i p = Eigen::Vector3d(std::floor(source[0]), std::floor(source[1]), std::floor(source[2])).cast<int>();
+        Eigen::Vector3i end = Eigen::Vector3d(std::floor(target[0]), std::floor(target[1]), std::floor(target[2])).cast<int>();
+        
+        if (voxel_set.insert(p).second)
+        {
+          std::cout << "missed p: " << p.transpose() << " between " << source.transpose() << " and " << target.transpose() << std::endl;
+          subsample.push_back(i);
+          continue; // only adding to one cell
+        }
+        
+        Eigen::Vector3i step(sign(target[0] - source[0]), sign(target[1] - source[1]), sign(target[2] - source[2]));
+        Eigen::Vector3d tmax, tdelta;
+        for (int j = 0; j<3; j++)
+        {
+          step[j] = sign(target[j] - source[j]);
+          double to = std::abs(source[j] - p[j] - (double)std::max(0, step[j]));        
+          double dir = std::max(std::numeric_limits<double>::min(), std::abs(source[j] - target[j]));
+          tmax[j] = to / dir;
+          tdelta[j] = 1.0 / dir;
+        }
+        
+        while (p != end) 
+        {
+          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 (tmax[0] < tmax[1]) 
+          {
+            if (tmax[0] < tmax[2]) 
+            {
+              p[0] += step[0];
+              tmax[0] += tdelta[0];
+            } 
+            else 
+            {
+              p[2] += step[2];
+              tmax[2] += tdelta[2];
+            }
+          } 
+          else 
+          {
+            if (tmax[1] < tmax[2]) 
+            {
+              p[1] += step[1];
+              tmax[1] += tdelta[1];
+            } 
+            else 
+            {
+              p[2] += step[2];
+              tmax[2] += tdelta[2];
+            }
+          }*/
+          if (voxel_set.insert(p).second)
+          {
+            std::cout << "missed p: " << p.transpose() << " between " << source.transpose() << " and " << target.transpose() << std::endl;
+            subsample.push_back(i);
+            break; // only adding to one cell
+          }          
+        }      
+       #endif
     }
     chunk.resize(subsample.size());
     for (int64_t i = 0; i < (int64_t)subsample.size(); i++)