Merge pull request #42 from AutoModality/develop

Updates of band alignment, DEM generation, and texturing projection

Author: usplm

contrib/pc2dem/pc2dem.py

@@ -47,6 +47,7 @@
                     args.type,
                     output_type='idw' if args.type == 'dtm' else 'max',
                     radiuses=list(map(str, radius_steps)),
+                    power=1,
                     gapfill=args.gapfill_steps > 0,
                     outdir=outdir,
                     resolution=args.resolution,

opendm/config.py

@@ -475,6 +475,12 @@ def config(argv=None, parser=None):
         help='Simple Morphological Filter window radius parameter (meters). '
                 'Default: %(default)s')
 
+    parser.add_argument('--texturing-use-dtm',
+                        action=StoreTrue,
+                        nargs=0,
+                        default=False,
+                        help=('Choose to use 2.5D DSM or DTM to project textures. Default:  %(default)s'))
+
     parser.add_argument('--texturing-data-term',
                         metavar='<string>',
                         action=StoreValue,

opendm/dem/commands.py

@@ -66,7 +66,7 @@ def rectify(lasFile, reclassify_threshold=5, min_area=750, min_points=500):
 
 error = None
 
-def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56'], gapfill=True,
+def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56'], power=1, gapfill=True,
                 outdir='', resolution=0.1, max_workers=1, max_tile_size=4096,
                 decimation=None, keep_unfilled_copy=False,
                 apply_smoothing=True):
@@ -156,9 +156,9 @@ def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56']
     def process_tile(q):
         log.ODM_INFO("Generating %s (%s, radius: %s, resolution: %s)" % (q['filename'], output_type, q['radius'], resolution))
 
-        d = pdal.json_gdal_base(q['filename'], output_type, q['radius'], resolution, q['bounds'])
+        d = pdal.json_gdal_base(q['filename'], output_type, q['radius'], power, resolution, q['bounds'])
 
-        if dem_type == 'dtm':
+        if dem_type == 'dtm' or dem_type == 'mesh_dtm':
             d = pdal.json_add_classification_filter(d, 2)
 
         if decimation is not None:

opendm/dem/pdal.py

@@ -49,14 +49,15 @@ def json_base():
     return {'pipeline': []}
 
 
-def json_gdal_base(filename, output_type, radius, resolution=1, bounds=None):
+def json_gdal_base(filename, output_type, radius, power=1, resolution=1, bounds=None):
     """ Create initial JSON for PDAL pipeline containing a Writer element """
     json = json_base()
 
     d = {
         'type': 'writers.gdal',
         'resolution': resolution,
         'radius': radius,
+        'power': power, # IDW distance exponent
         'filename': filename,
         'output_type': output_type,
         'data_type': 'float'

opendm/mesh.py

@@ -9,7 +9,8 @@
 from scipy import signal
 import numpy as np
 
-def create_25dmesh(inPointCloud, outMesh, radius_steps=["0.05"], dsm_resolution=0.05, depth=8, samples=1, maxVertexCount=100000, available_cores=None, method='gridded', smooth_dsm=True):
+def create_25dmesh(inPointCloud, outMesh, radius_steps=["0.05"], dsm_resolution=0.05, depth=8, samples=1, maxVertexCount=100000, 
+                   available_cores=None, method='gridded', smooth_dsm=True, use_dtm=False):
     # Create DSM from point cloud
 
     # Create temporary directory
@@ -22,11 +23,15 @@ def create_25dmesh(inPointCloud, outMesh, radius_steps=["0.05"], dsm_resolution=
 
     log.ODM_INFO('Creating DSM for 2.5D mesh')
 
+    mesh_dem = 'mesh_dtm' if use_dtm else 'mesh_dsm'
+    output_type = 'max' if use_dtm else 'idw'
+
     commands.create_dem(
             inPointCloud,
-            'mesh_dsm',
-            output_type='idw',
+            mesh_dem,
+            output_type=output_type,
             radiuses=radius_steps,
+            power=1,
             gapfill=True,
             outdir=tmp_directory,
             resolution=dsm_resolution,
@@ -35,9 +40,9 @@ def create_25dmesh(inPointCloud, outMesh, radius_steps=["0.05"], dsm_resolution=
         )
 
     if method == 'gridded':
-        mesh = dem_to_mesh_gridded(os.path.join(tmp_directory, 'mesh_dsm.tif'), outMesh, maxVertexCount, maxConcurrency=max(1, available_cores))
+        mesh = dem_to_mesh_gridded(os.path.join(tmp_directory, mesh_dem + '.tif'), outMesh, maxVertexCount, maxConcurrency=max(1, available_cores))
     elif method == 'poisson':
-        dsm_points = dem_to_points(os.path.join(tmp_directory, 'mesh_dsm.tif'), os.path.join(tmp_directory, 'dsm_points.ply'))
+        dsm_points = dem_to_points(os.path.join(tmp_directory, mesh_dem + '.tif'), os.path.join(tmp_directory, 'dsm_points.ply'))
         mesh = screened_poisson_reconstruction(dsm_points, outMesh, depth=depth, 
                                     samples=samples, 
                                     maxVertexCount=maxVertexCount, 

opendm/multispectral.py

@@ -552,16 +552,19 @@ def compute_using(algorithm):
         log.ODM_WARNING("Compute homography: %s" % str(e))
         return None, (None, None), None
 
-def find_ecc_homography(image_gray, align_image_gray, number_of_iterations=1000, termination_eps=1e-8, start_eps=1e-4, warp_matrix_init=None):
+def find_ecc_homography(image_gray, align_image_gray, number_of_iterations=2000, termination_eps=1e-8, start_eps=1e-4, warp_matrix_init=None):
     # Major props to Alexander Reynolds for his insight into the pyramided matching process found at
     # https://stackoverflow.com/questions/45997891/cv2-motion-euclidean-for-the-warp-mode-in-ecc-image-alignment-method
     pyramid_levels = 0
     h,w = image_gray.shape
     min_dim = min(h, w)
 
-    number_of_iterations = 1000 if min_dim > 300 else 5000
-    termination_eps = 1e-8 if min_dim > 300 else 1e-6
-    gaussian_filter_size = 9 # if min_dim > 300 else 5
+    if (min_dim <= 300):
+        number_of_iterations = 5000
+        termination_eps = 1e-6
+        gaussian_filter_size = 9 # a constant since there is only one pyramid level
+    else:
+        gaussian_filter_size = 5 # will be increased in each pyramid level iteration
 
     while min_dim > 300:
         min_dim /= 2.0
@@ -609,8 +612,9 @@ def find_ecc_homography(image_gray, align_image_gray, number_of_iterations=1000,
                 number_of_iterations, eps)
 
         try:
+            gaussian_filter_size = gaussian_filter_size + level * 2
             log.ODM_INFO("Computing ECC pyramid level %s using Gaussian filter size %s" % (level, gaussian_filter_size))
-            _, warp_matrix = cv2.findTransformECC(ig, aig, warp_matrix, cv2.MOTION_HOMOGRAPHY, criteria, inputMask=None, gaussFiltSize=gaussian_filter_size)
+            _, warp_matrix = cv2.findTransformECC(ig, aig, warp_matrix, cv2.MOTION_HOMOGRAPHY, criteria, inputMask=None, gaussFiltSize=gaussian_filter_size)            
         except Exception as e:
             if level != pyramid_levels:
                 log.ODM_INFO("Could not compute ECC warp_matrix at pyramid level %s, resetting matrix" % level)
@@ -625,7 +629,7 @@ def find_ecc_homography(image_gray, align_image_gray, number_of_iterations=1000,
 
     return warp_matrix
 
-def find_features_homography(image_gray, align_image_gray, feature_retention=0.7, min_match_count=10):
+def find_features_homography(image_gray, align_image_gray, feature_retention=0.8, min_match_count=4):
 
     # Detect SIFT features and compute descriptors.
     detector = cv2.SIFT_create() # edgeThreshold=10, contrastThreshold=0.1 (default 0.04)

opendm/point_cloud.py

@@ -89,11 +89,12 @@ def filter(input_point_cloud, output_point_cloud, output_stats, standard_deviati
         log.ODM_INFO("Sampling points around a %sm radius" % sample_radius)
         args.append('--radius %s' % sample_radius)
 
-    meank = 16
-    log.ODM_INFO("Filtering {} (statistical, meanK {}, standard deviation {})".format(input_point_cloud, meank, standard_deviation))
-    args.append('--meank %s' % meank)
-    args.append('--std %s' % standard_deviation)
-    args.append('--stats "%s"' % output_stats)
+    if standard_deviation > 0:
+        meank = 16
+        log.ODM_INFO("Filtering {} (statistical, meanK {}, standard deviation {})".format(input_point_cloud, meank, standard_deviation))
+        args.append('--meank %s' % meank)
+        args.append('--std %s' % standard_deviation)
+        args.append('--stats "%s"' % output_stats)
 
     if boundary is not None:
         log.ODM_INFO("Boundary {}".format(boundary))

opendm/types.py

@@ -290,6 +290,7 @@ def __init__(self, root_path, gcp_file = None, geo_file = None, align_file = Non
 
         # filter points
         self.filtered_point_cloud = os.path.join(self.odm_filterpoints, "point_cloud.ply")
+        self.filtered_point_cloud_classified = os.path.join(self.odm_filterpoints, "point_cloud_classified.ply")
         self.filtered_point_cloud_stats = os.path.join(self.odm_filterpoints, "point_cloud_stats.json")
 
         # odm_meshing

stages/mvstex.py

@@ -36,6 +36,7 @@ def add_run(nvm_file, primary=True, band=None):
                     'model': tree.odm_mesh,
                     'nadir': False,
                     'primary': primary,
+                    'band': band,
                     'nvm_file': nvm_file,
                     'labeling_file': os.path.join(tree.odm_texturing, "odm_textured_model_geo_labeling.vec") if subdir else None
                 }]
@@ -46,6 +47,7 @@ def add_run(nvm_file, primary=True, band=None):
                     'model': tree.odm_25dmesh,
                     'nadir': True,
                     'primary': primary,
+                    'band': band,
                     'nvm_file': nvm_file,
                     'labeling_file': os.path.join(tree.odm_25dtexturing, "odm_textured_model_geo_labeling.vec") if subdir else None
                 }]
@@ -80,6 +82,8 @@ def add_run(nvm_file, primary=True, band=None):
                     os.unlink(unaligned_obj)
 
                 # Format arguments to fit Mvs-Texturing app
+                dataTerm = args.texturing_data_term
+                outlierRemovalType = args.texturing_outlier_removal_type
                 skipGlobalSeamLeveling = ""
                 skipLocalSeamLeveling = ""
                 keepUnseenFaces = ""
@@ -94,13 +98,21 @@ def add_run(nvm_file, primary=True, band=None):
                 if (r['nadir']):
                     nadir = '--nadir_mode'
 
+                # thermal band
+                if (r['band'] == 'lwir'):
+                    dataTerm = "gmi"
+                    outlierRemovalType = "none"
+                    skipGlobalSeamLeveling = ""
+                    skipLocalSeamLeveling = ""
+                    keepUnseenFaces = ""
+
                 # mvstex definitions
                 kwargs = {
                     'bin': context.mvstex_path,
                     'out_dir': os.path.join(r['out_dir'], "odm_textured_model_geo"),
                     'model': r['model'],
-                    'dataTerm': args.texturing_data_term,
-                    'outlierRemovalType': args.texturing_outlier_removal_type,
+                    'dataTerm': dataTerm,
+                    'outlierRemovalType': outlierRemovalType,
                     'skipGlobalSeamLeveling': skipGlobalSeamLeveling,
                     'skipLocalSeamLeveling': skipLocalSeamLeveling,
                     'keepUnseenFaces': keepUnseenFaces,

stages/odm_dem.py

@@ -102,6 +102,7 @@ def process(self, args, outputs):
                             product,
                             output_type=args.dtm_interpolation if product == 'dtm' else args.dsm_interpolation,
                             radiuses=list(map(str, radius_steps)),
+                            power=2,
                             gapfill=args.dem_gapfill_steps > 0,
                             outdir=odm_dem_root,
                             resolution=resolution / 100.0,