Merge pull request #978 from pierotofy/clean25dmesh

Clean 2.5D mesh, keep vertex count in check

opendm/dem/commands.py

@@ -183,6 +183,7 @@ def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56']
     vrt_path = os.path.abspath(os.path.join(outdir, "merged.vrt"))
     run('gdalbuildvrt "%s" "%s"' % (vrt_path, '" "'.join(map(lambda t: t['filename'], tiles))))
 
+    geotiff_tmp_path = os.path.abspath(os.path.join(outdir, 'merged.tmp.tif'))
     geotiff_path = os.path.abspath(os.path.join(outdir, 'merged.tif'))
 
     # Build GeoTIFF
@@ -190,16 +191,25 @@ def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56']
         'max_memory': get_max_memory(),
         'threads': max_workers if max_workers else 'ALL_CPUS',
         'vrt': vrt_path,
-        'geotiff': geotiff_path
+        'geotiff': geotiff_path,
+        'geotiff_tmp': geotiff_tmp_path
     }
 
     if gapfill:
+        # Sometimes, for some reason gdal_fillnodata.py
+        # behaves strangely when reading data directly from a .VRT
+        # so we need to convert to GeoTIFF first.
+        run('gdal_translate '
+                '-co NUM_THREADS={threads} '
+                '--config GDAL_CACHEMAX {max_memory}% '
+                '{vrt} {geotiff_tmp}'.format(**kwargs))
+
         run('gdal_fillnodata.py '
             '-co NUM_THREADS={threads} '
             '--config GDAL_CACHEMAX {max_memory}% '
             '-b 1 '
             '-of GTiff '
-            '{vrt} {geotiff}'.format(**kwargs))
+            '{geotiff_tmp} {geotiff}'.format(**kwargs))
     else:
         run('gdal_translate '
                 '-co NUM_THREADS={threads} '
@@ -209,6 +219,7 @@ def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56']
     post_process(geotiff_path, output_path)
     os.remove(geotiff_path)
 
+    if os.path.exists(geotiff_tmp_path): os.remove(geotiff_tmp_path)
     if os.path.exists(vrt_path): os.remove(vrt_path)
     for t in tiles:
         if os.path.exists(t['filename']): os.remove(t['filename'])

opendm/mesh.py

@@ -74,12 +74,22 @@ def dem_to_points(inGeotiff, outPointCloud, verbose=False):
     return outPointCloud
 
 
-def dem_to_mesh_gridded(inGeotiff, outPointCloud, maxVertexCount, verbose=False):
+def dem_to_mesh_gridded(inGeotiff, outMesh, maxVertexCount, verbose=False):
     log.ODM_INFO('Creating mesh from DSM: %s' % inGeotiff)
 
+    mesh_path, mesh_filename = os.path.split(outMesh)
+    # mesh_path = path/to
+    # mesh_filename = odm_mesh.ply
+
+    basename, ext = os.path.splitext(mesh_filename)
+    # basename = odm_mesh
+    # ext = .ply
+
+    outMeshDirty = os.path.join(mesh_path, "{}.dirty{}".format(basename, ext))
+
     kwargs = {
         'bin': context.dem2mesh_path,
-        'outfile': outPointCloud,
+        'outfile': outMeshDirty,
         'infile': inGeotiff,
         'maxVertexCount': maxVertexCount,
         'verbose': '-verbose' if verbose else ''
@@ -90,7 +100,25 @@ def dem_to_mesh_gridded(inGeotiff, outPointCloud, maxVertexCount, verbose=False)
          '-maxVertexCount {maxVertexCount} '
          ' {verbose} '.format(**kwargs))
 
-    return outPointCloud
+    # Cleanup and reduce vertex count if necessary 
+    # (as dem2mesh cannot guarantee that we'll have the target vertex count)
+    cleanupArgs = {
+        'bin': context.odm_modules_path,
+        'outfile': outMesh,
+        'infile': outMeshDirty,
+        'max_vertex': maxVertexCount,
+        'verbose': '-verbose' if verbose else ''
+    }
+
+    system.run('{bin}/odm_cleanmesh -inputFile {infile} '
+         '-outputFile {outfile} '
+         '-removeIslands '
+         '-decimateMesh {max_vertex} {verbose} '.format(**cleanupArgs))
+
+    # Delete intermediate results
+    os.remove(outMeshDirty)
+
+    return outMesh
 
 
 def screened_poisson_reconstruction(inPointCloud, outMesh, depth = 8, samples = 1, maxVertexCount=100000, pointWeight=4, threads=context.num_cores, verbose=False):