Merge pull request #2141 from ThorstenHarter/master

Avoid huge index jumps in RandomSample

2d/CMakeLists.txt

@@ -1,6 +1,6 @@
 set(SUBSYS_NAME 2d)
 set(SUBSYS_DESC "Point cloud 2d")
-set(SUBSYS_DEPS common io filters)
+set(SUBSYS_DEPS common filters)
 
 set(build TRUE)
 PCL_SUBSYS_OPTION(build "${SUBSYS_NAME}" "${SUBSYS_DESC}" ON)

filters/include/pcl/filters/impl/random_sample.hpp

@@ -92,9 +92,8 @@ template<typename PointT>
 void
 pcl::RandomSample<PointT>::applyFilter (std::vector<int> &indices)
 {
-  unsigned N = static_cast<unsigned> (indices_->size ());
-
-  unsigned int sample_size = negative_ ? N - sample_ : sample_;
+  size_t N = indices_->size ();  
+  size_t sample_size = negative_ ? N - sample_ : sample_;
   // If sample size is 0 or if the sample size is greater then input cloud size
   //   then return all indices
   if (sample_size >= N)
@@ -105,48 +104,44 @@ pcl::RandomSample<PointT>::applyFilter (std::vector<int> &indices)
   else
   {
     // Resize output indices to sample size
-    indices.resize (static_cast<size_t> (sample_size));
+    indices.resize (sample_size);
     if (extract_removed_indices_)
-      removed_indices_->resize (static_cast<size_t> (N - sample_size));
+      removed_indices_->resize (N - sample_size);
 
     // Set random seed so derived indices are the same each time the filter runs
     std::srand (seed_);
 
-    // Algorithm A
-    unsigned top = N - sample_size;
-    unsigned i = 0;
-    unsigned index = 0;
+    // Algorithm S
+    size_t i = 0;
+    size_t index = 0;
     std::vector<bool> added;
     if (extract_removed_indices_)
       added.resize (indices_->size (), false);
-    for (size_t n = sample_size; n >= 2; n--)
+    size_t n = sample_size;
+    while (n > 0)
     {
-      float V = unifRand ();
-      unsigned S = 0;
-      float quot = static_cast<float> (top) / static_cast<float> (N);
-      while (quot > V)
+      // Step 1: [Generate U.] Generate a random variate U that is uniformly distributed between 0 and 1.
+      const float U = unifRand ();
+      // Step 2: [Test.] If N * U > n, go to Step 4. 
+      if ((N * U) <= n)
       {
-        S++;
-        top--;
-        N--;
-        quot = quot * static_cast<float> (top) / static_cast<float> (N);
+        // Step 3: [Select.] Select the next record in the file for the sample, and set n : = n - 1.
+        if (extract_removed_indices_)
+          added[index] = true;
+        indices[i++] = (*indices_)[index];
+        --n;
       }
-      index += S;
-      if (extract_removed_indices_)
-        added[index] = true;
-      indices[i++] = (*indices_)[index++];
-      N--;
+      // Step 4: [Don't select.] Skip over the next record (do not include it in the sample).      
+      // Set N : = N - 1.
+      --N;
+      ++index;
+      // If n > 0, then return to Step 1; otherwise, the sample is complete and the algorithm terminates.
     }
 
-    index += N * static_cast<unsigned> (unifRand ());
-    if (extract_removed_indices_)
-      added[index] = true;
-    indices[i++] = (*indices_)[index++];
-
     // Now populate removed_indices_ appropriately
     if (extract_removed_indices_)
     {
-      unsigned ri = 0;
+      size_t ri = 0;
       for (size_t i = 0; i < added.size (); i++)
       {
         if (!added[i])

filters/include/pcl/filters/random_sample.h

@@ -232,7 +232,6 @@ namespace pcl
       unifRand ()
       {
         return (static_cast<float> (rand () / double (RAND_MAX)));
-        //return (((214013 * seed_ + 2531011) >> 16) & 0x7FFF);
       }
    };
 }