Fix duplicate rows in ordered pointsets (#25)

* Fix bugs in segmentation algorithms which created duplicate rows

* Add `vc_repair_pointsets` for fixing existing pointsets

segmentation/include/vc/segmentation/OpticalFlowSegmentation.hpp

@@ -102,7 +102,10 @@ class OpticalFlowSegmentation : public ChainSegmentationAlgorithm
     [[nodiscard]] auto progressIterations() const -> size_t override;
 
 private:
-    /** @brief Compute the segmentation 1 Line */
+    /**
+     * @brief Compute the curve for z + 1 given a curve on z using the optical
+     * flow between the two slices
+     */
     auto compute_curve_(const FittedCurve& currentCurve, int zIndex)
         -> std::vector<Voxel>;
 

segmentation/src/LocalResliceParticleSim.cpp

@@ -76,8 +76,9 @@ LocalResliceSegmentation::PointSet LocalResliceSegmentation::compute()
     points.push_back(currentVs);
 
     // Iterate over z-slices
+    auto stepSize = static_cast<int>(stepSize_);
     size_t iteration{0};
-    for (int zIndex = startIndex; zIndex <= endIndex_; zIndex += stepSize_) {
+    for (int zIndex = startIndex; zIndex < endIndex_; zIndex += stepSize) {
         // Update progress
         progressUpdated(iteration++);
 
@@ -115,7 +116,7 @@ LocalResliceSegmentation::PointSet LocalResliceSegmentation::compute()
         // XXX DEBUG
         for (int i = 0; i < int(currentCurve.size()); ++i) {
             // Estimate normal and reslice along it
-            const cv::Vec3d normal = estimate_normal_at_index_(currentCurve, i);
+            const auto normal = estimate_normal_at_index_(currentCurve, i);
             const auto reslice = vol_->reslice(
                 currentCurve(i), normal, {0, 0, 1}, resliceSize_, resliceSize_);
             reslices.push_back(reslice);

segmentation/src/OpticalFlowSegmentation.cpp

@@ -300,7 +300,7 @@ auto OpticalFlowSegmentation::compute() -> PointSet
     std::size_t iteration{0};
     auto stepSize = static_cast<int>(stepSize_);
     const int padding = vol_->numSlices();
-    for (int zIndex = startIndex; zIndex <= endIndex_; zIndex += stepSize) {
+    for (int zIndex = startIndex; zIndex < endIndex_; zIndex += stepSize) {
         // Update progress
         progressUpdated(iteration++);
 
@@ -393,7 +393,7 @@ auto OpticalFlowSegmentation::compute() -> PointSet
         }
 
         // Generate nextVs by evenly spacing points in the stitched curve
-        FittedCurve stitchedFittedCurve(stitched, zIndex);
+        FittedCurve stitchedFittedCurve(stitched, zIndex + 1);
         auto nextVs = stitchedFittedCurve.evenlySpacePoints();
 
         // Check if any points in nextVs are outside volume boundaries. If so,

utils/CMakeLists.txt

@@ -164,6 +164,11 @@ add_executable(vc_visualize_graph src/VisualizeRenderGraph.cpp)
 target_link_libraries(vc_visualize_graph VC::core VC::graph Boost::program_options smgl::smgl)
 list(APPEND utils_install_list vc_visualize_graph)
 
+# vc_fix_pointset_duplicates
+add_executable(vc_repair_pointsets src/RepairPointSets.cpp)
+target_link_libraries(vc_repair_pointsets VC::core Boost::program_options)
+list(APPEND utils_install_list vc_repair_pointsets)
+
 # Install targets
 if(VC_INSTALL_UTILS)
 install(

utils/src/ConvertPointSet.cpp

@@ -103,17 +103,15 @@ void PointSetToMesh(const fs::path& inputPath, const fs::path& outputPath)
             [](const auto& l, const auto& r) { return l[2] < r[2]; });
 
         // Get color info
-        for (const auto& pt :
-             vc::ProgressWrap(inputCloud, "Getting point intensities")) {
+        for (const auto& pt : inputCloud) {
             intensities.emplace_back(volume->interpolateAt(pt));
         }
     }
 
     // Convert to ITKMesh
     vc::ITKPoint tmpPt;
     auto mesh = vc::ITKMesh::New();
-    for (const auto it :
-         vc::ProgressWrap(enumerate(inputCloud), "Converting points")) {
+    for (const auto it : enumerate(inputCloud)) {
         tmpPt[0] = it.second[0];
         tmpPt[1] = it.second[1];
         tmpPt[2] = it.second[2];

utils/src/RepairPointSets.cpp

@@ -0,0 +1,106 @@
+#include <boost/program_options.hpp>
+
+#include "vc/core/filesystem.hpp"
+#include "vc/core/io/PointSetIO.hpp"
+#include "vc/core/types/OrderedPointSet.hpp"
+#include "vc/core/util/Iteration.hpp"
+#include "vc/core/util/Logging.hpp"
+
+namespace fs = volcart::filesystem;
+namespace po = boost::program_options;
+namespace vc = volcart;
+using vc::enumerate;
+using vc::range;
+
+using psio = vc::PointSetIO<cv::Vec3d>;
+using PointSet = vc::OrderedPointSet<cv::Vec3d>;
+
+auto main(int argc, char* argv[]) -> int
+{
+    ///// Parse the command line options /////
+    // All command line options
+    // clang-format off
+    po::options_description all("Usage");
+    all.add_options()
+        ("help,h", "Show this message")
+        ("input,i", po::value<std::string>()->required(),
+             "Path to the input PointSet")
+        ("output,o", po::value<std::string>(), "Path to the output PointSet")
+        ("verbose,v", "Print verbose information");
+    // clang-format on
+
+    // parsed will hold the values of all parsed options as a Map
+    po::variables_map args;
+    po::store(po::command_line_parser(argc, argv).options(all).run(), args);
+
+    // Show the help message
+    if (args.count("help") > 0 or argc < 3) {
+        std::cout << all << std::endl;
+        return EXIT_SUCCESS;
+    }
+
+    // Warn of missing options
+    try {
+        po::notify(args);
+    } catch (po::error& e) {
+        vc::Logger()->error(e.what());
+        return EXIT_FAILURE;
+    }
+    // Logging
+    auto verbose = args.count("verbose") > 0;
+
+    // Load the pointset
+    const fs::path inputPath = args["input"].as<std::string>();
+    auto cloud = psio::ReadOrderedPointSet(inputPath);
+
+    // Report the point set properties
+    auto cols = cloud.width();
+    auto rows = cloud.height();
+    auto startZ = static_cast<std::size_t>(cloud[0][2]);
+    auto endZ = static_cast<std::size_t>(cloud.getRow(rows - 1)[0][2]);
+    vc::Logger()->info(
+        "Original pointset :: Shape: ({}, {}), Z-Range: [{}, {}]", rows, cols,
+        startZ, endZ);
+
+    // Remap all the rows
+    int repaired{0};
+    vc::Logger()->info("Repairing pointset...");
+    for (const auto row : range(rows)) {
+        // get the z-values
+        auto newZ = row + startZ;
+        auto oldZ = static_cast<std::size_t>(cloud(row, 0)[2]);
+
+        // if the z-values don't match...
+        if (oldZ != newZ) {
+            // report this row if verbose
+            if (verbose) {
+                std::cout << "  [" << row << "] ";
+                std::cout << oldZ << " -> " << newZ;
+                std::cout << "\n";
+            }
+            // update the row
+            for (const auto col : range(cols)) {
+                cloud(row, col)[2] = static_cast<double>(newZ);
+            }
+            // track num repairs
+            repaired++;
+        }
+    }
+    vc::Logger()->info("Repaired {} rows", repaired);
+
+    // Report the new point set properties
+    cols = cloud.width();
+    rows = cloud.height();
+    startZ = static_cast<std::size_t>(cloud[0][2]);
+    endZ = static_cast<std::size_t>(cloud.getRow(rows - 1)[0][2]);
+    vc::Logger()->info(
+        "New pointset :: Shape: ({}, {}), Z-Range: [{}, {}]", rows, cols,
+        startZ, endZ);
+
+    // (optional) Save the new point cloud
+    if (args.count("output") > 0) {
+        vc::Logger()->info("Writing pointset...");
+        const fs::path outputPath = args["output"].as<std::string>();
+        psio::WriteOrderedPointSet(outputPath, cloud);
+    }
+}