Added confidence filtering in odm_filterpoints

modules/odm_filterpoints/src/main.cpp

@@ -1,4 +1,5 @@
 #include <iostream>
+#include <algorithm>
 #include <pdal/filters/OutlierFilter.hpp>
 #include <pdal/filters/RangeFilter.hpp>
 #include "CmdLineParser.h"
@@ -13,12 +14,13 @@ cmdLineParameter< char* >
     OutputFile( "outputFile" );
 cmdLineParameter< float >
     StandardDeviation( "sd" ) ,
-    MeanK ( "meank" );
+    MeanK ( "meank" ) ,
+    Confidence ( "confidence" );
 cmdLineReadable
 	Verbose( "verbose" );
 
 cmdLineReadable* params[] = {
-    &InputFile , &OutputFile , &StandardDeviation, &MeanK, &Verbose ,
+    &InputFile , &OutputFile , &StandardDeviation, &MeanK, &Confidence, &Verbose ,
     NULL
 };
 
@@ -28,6 +30,7 @@ void help(char *ex){
               << "\t -" << OutputFile.name << " <output PLY point cloud>" << std::endl
               << "\t [-" << StandardDeviation.name << " <standard deviation threshold>]" << std::endl
               << "\t [-" << MeanK.name << " <mean number of neighbors >]" << std::endl
+              << "\t [-" << Confidence.name << " <lower bound filter for confidence property>]" << std::endl
 
               << "\t [-" << Verbose.name << "]" << std::endl;
     exit(EXIT_FAILURE);
@@ -65,13 +68,31 @@ int main(int argc, char **argv) {
     pdal::FloatPlyReader plyReader;
     plyReader.setOptions(inPlyOpts);
 
+    pdal::RangeFilter confidenceFilter;
+    if (Confidence.set){
+        pdal::Options confidenceFilterOpts;
+        float confidenceValue = std::min(1.0f, std::max(Confidence.value, 0.0f));
+        std::ostringstream confidenceLimit;
+        confidenceLimit << "confidence[" << confidenceValue << ":1]";
+        confidenceFilterOpts.add("limits", confidenceLimit.str());
+
+        confidenceFilter.setInput(plyReader);
+        confidenceFilter.setOptions(confidenceFilterOpts);
+    }
+
     pdal::Options outlierOpts;
     outlierOpts.add("method", "statistical");
     outlierOpts.add("mean_k", MeanK.value);
     outlierOpts.add("multiplier", StandardDeviation.value);
 
     pdal::OutlierFilter outlierFilter;
-    outlierFilter.setInput(plyReader);
+    if (Confidence.set){
+        logWriter("Filtering confidence\n");
+        outlierFilter.setInput(confidenceFilter);
+    }else{
+        outlierFilter.setInput(plyReader);
+
+    }
     outlierFilter.setOptions(outlierOpts);
 
     pdal::Options rangeOpts;

opendm/config.py

@@ -171,6 +171,15 @@ def config():
                         help='Run local bundle adjustment for every image added to the reconstruction and a global '
                              'adjustment every 100 images. Speeds up reconstruction for very large datasets.')
 
+    parser.add_argument('--mve-confidence',
+                        metavar='<float: 0 <= x <= 1>',
+                        type=float,
+                        default=0.1,
+                        help=('Discard points that have less than a certain confidence threshold. '
+                              'This only affects dense reconstructions performed with MVE. '
+                              'Higher values discard more points. '
+                              'Default: %(default)s'))
+
     parser.add_argument('--use-3dmesh',
                     action='store_true',
                     default=False,

opendm/point_cloud.py

@@ -3,7 +3,7 @@ from opendm import system
 from opendm import log
 from opendm import context
 
-def filter(inputPointCloud, outputPointCloud, standard_deviation=2.5, meank=16, verbose=False):
+def filter(inputPointCloud, outputPointCloud, standard_deviation=2.5, meank=16, confidence=None, verbose=False):
     """
     Filters a point cloud
     """
@@ -12,6 +12,8 @@ def filter(inputPointCloud, outputPointCloud, standard_deviation=2.5, meank=16,
         return
 
     log.ODM_INFO("Filtering point cloud (statistical, meanK {}, standard deviation {})".format(meank, standard_deviation))
+    if confidence:
+        log.ODM_INFO("Keeping only points with > %s confidence" % confidence)
 
     if not os.path.exists(inputPointCloud):
         log.ODM_ERROR("{} does not exist, cannot filter point cloud. The program will now exit.".format(inputPointCloud))
@@ -29,13 +31,14 @@ def filter(inputPointCloud, outputPointCloud, standard_deviation=2.5, meank=16,
       'outputFile': outputPointCloud,
       'sd': standard_deviation,
       'meank': meank,
-      'verbose': '--verbose' if verbose else '',
+      'verbose': '-verbose' if verbose else '',
+      'confidence': '-confidence %s' % confidence if confidence else '',
     }
 
     system.run('{bin} -inputFile {inputFile} '
          '-outputFile {outputFile} '
          '-sd {sd} '
-         '-meank {meank} {verbose} '.format(**filterArgs))
+         '-meank {meank} {confidence} {verbose} '.format(**filterArgs))
 
     # Remove input file, swap temp file
     if not os.path.exists(outputPointCloud):

scripts/mve.py

@@ -59,29 +59,43 @@ class ODMMveCell(ecto.Cell):
             if not io.dir_exists(tree.mve_views):
                 system.run('%s %s %s' % (context.makescene_path, tree.mve_path, tree.mve))
 
+            # Compute mve output scale based on depthmap_resolution
+            max_width = 0
+            max_height = 0
+            for photo in photos:
+                max_width = max(photo.width, max_width)
+                max_height = max(photo.height, max_height)
+
+            min_depthmap_pixels = args.depthmap_resolution * args.depthmap_resolution
+            image_pixels = max_width * max_height
+
+            mve_output_scale = 0
+            while True:
+                image_pixels = int(image_pixels / 4)
+                if image_pixels <= min_depthmap_pixels: break
+                mve_output_scale += 1
+
+            # Can we get closer to the requested resolution by dividing one more time?
+            if abs(int(image_pixels / 4) - min_depthmap_pixels) < abs(image_pixels - min_depthmap_pixels):
+                mve_output_scale += 1
+
             dmrecon_config = [
-                "--max-pixels=%s" % int(args.depthmap_resolution * args.depthmap_resolution),
-                # "-s%s" % args.mve_output_scale,
+                "-s%s" % mve_output_scale,
 	            "--progress=silent",
                 "--force",
             ]
 
             # Run MVE's dmrecon
+            log.ODM_INFO("Running dense reconstruction. This might take a while. Please be patient, the process is not dead or hang.")
             system.run('%s %s %s' % (context.dmrecon_path, ' '.join(dmrecon_config), tree.mve))
 
             scene2pset_config = [
-                "-F1",
-                "--with-conf"
+                "-F%s" % mve_output_scale
             ]
 
             # run scene2pset
             system.run('%s %s "%s" "%s"' % (context.scene2pset_path, ' '.join(scene2pset_config), tree.mve, tree.mve_model))
 
-            # run pdal
-            # system.run('pdal translate -i {} '
-            #            ' -o {} '
-            #            ' -f range --filters.range.limits="confidence[0.25:1]'.format(tree.mve_model, tree.mve_model_output))
-        
             # find and rename the output file for simplicity
             mve_files = glob.glob(os.path.join(tree.mve, 'mve-*'))
             mve_files.sort(key=os.path.getmtime) # sort by last modified date

scripts/odm_filterpoints.py

@@ -41,8 +41,11 @@ class ODMFilterPoints(ecto.Cell):
             else:
                 inputPointCloud = tree.mve_model
 
-            # TODO add MVE filter
-            point_cloud.filter(inputPointCloud, tree.filtered_point_cloud, standard_deviation=args.pc_filter, verbose=args.verbose)
+            confidence = None
+            if not args.use_opensfm_dense and not args.fast_orthophoto:
+                confidence = args.mve_confidence
+
+            point_cloud.filter(inputPointCloud, tree.filtered_point_cloud, standard_deviation=args.pc_filter, confidence=confidence, verbose=args.verbose)
         else:
             log.ODM_WARNING('Found a valid point cloud file in: %s' %
                             tree.filtered_point_cloud)