Merge pull request #1 from OpenDroneMap/master

Pull in latest master

contrib/grass/odm_grass.py

@@ -0,0 +1,172 @@
+#!/usr/bin/env python
+
+# To run, set the following env variables:
+# PYTHONHOME location of Python
+# PYTHONPATH location of GRASS Python libs
+# PATH include GRASS bin and lib
+# GISBASE location of GRASS
+
+import os
+import sys
+import grass.script as gscript
+import grass.script.core
+import grass.script.setup
+
+vlidarName = 'odm_vlidar'
+rsurfName = 'odm_rsurf'
+contourName = 'odm_contour'
+orthophotoName = 'odm_orthophoto'
+reliefName = 'odm_relief'
+shadedReliefName = reliefName + '_shaded'
+
+overwrite = True
+
+
+def main():
+    if len(sys.argv) < 2:
+        sys.exit('Please provide the ODM project path.')
+
+    projectHome = sys.argv[1]
+
+    gisdb = projectHome+'/grassdata'
+    location = 'odm'
+    gisrc = gscript.setup.init(os.environ['GISBASE'], gisdb, location)
+
+    # get srs and initial extents
+    with open(projectHome+'/odm_georeferencing/coords.txt') as f:
+        srs = f.readline().split()
+        mean = f.readline().split()
+        meanX = float(mean[0])
+        meanY = float(mean[1])
+        minX = float('inf')
+        maxX = float('-inf')
+        minY = float('inf')
+        maxY = float('-inf')
+        for line in f:
+            xy = line.split()
+            x = float(xy[0])
+            y = float(xy[1])
+            minX = min(x, minX)
+            maxX = max(x, maxX)
+            minY = min(y, minY)
+            maxY = max(y, maxY)
+
+    datum = srs[0]
+    proj = srs[1]
+    zone = srs[2]
+    gscript.core.create_location(gisdb, location, datum=datum,
+                                 proj4='+proj='+proj+' +zone='+zone,
+                                 overwrite=overwrite)
+
+    n = meanY + maxY
+    s = meanY + minY
+    e = meanX + maxX
+    w = meanX + minX
+    gscript.run_command('g.region', flags='s', n=n, s=s, e=e, w=w, res=0.01,
+                        res3=0.01, overwrite=overwrite)
+
+    dem(projectHome)
+    contour(projectHome)
+    relief(projectHome)
+
+    os.remove(gisrc)
+
+
+def dem(projectHome):
+    """
+    Creates a DEM in GeoTIFF format.
+    NB: this is a data raster, not an RGBA raster
+    and so is normally only readable by GIS and not image software.
+    """
+    print 'Creating DEM'
+
+    step = 0.5
+    gscript.run_command('v.in.lidar', flags='beo',
+                        input=projectHome +
+                        '/odm_georeferencing/odm_georeferenced_model.ply.las',
+                        output=vlidarName, overwrite=overwrite)
+
+    gscript.run_command('v.surf.bspline', input=vlidarName,
+                        raster_output=rsurfName,
+                        ew_step=step, ns_step=step, method='bicubic',
+                        memory=4096, overwrite=overwrite)
+
+    gscript.run_command('r.out.gdal', flags='cfm', input=rsurfName,
+                        output=projectHome+'/odm_georeferencing/odm_dem.tif',
+                        format='GTiff', type='Float32',
+                        createopt='TILED=yes,COMPRESS=DEFLATE,PREDICTOR=2,' +
+                        'BLOCKXSIZE=512,BLOCKYSIZE=512', nodata=0,
+                        overwrite=overwrite)
+
+
+def contour(projectHome):
+    """
+    Creates a contour map.
+    """
+    print 'Creating contour map'
+
+    step = 0.25
+    gscript.run_command('r.contour', input=rsurfName, output=contourName,
+                        step=step, overwrite=overwrite)
+
+    gscript.run_command('v.out.ogr', input=contourName,
+                        output=projectHome +
+                        '/odm_georeferencing/odm_contour.shp',
+                        overwrite=overwrite)
+
+
+def relief(projectHome):
+    """
+    Creates a textured relief map in GeoTIFF format.
+    NB: this is an RGBA raster and so is readable by image software.
+    """
+    print 'Creating relief map'
+
+    gscript.run_command('r.in.gdal', flags='e',
+                        input=projectHome+'/odm_orthophoto/odm_orthophoto.tif',
+                        output=orthophotoName, memory=2047,
+                        overwrite=overwrite)
+
+    gscript.run_command('r.composite', red=orthophotoName+'.red',
+                        green=orthophotoName+'.green',
+                        blue=orthophotoName+'.blue',
+                        output=orthophotoName+'.rgb',
+                        overwrite=overwrite)
+
+    gscript.run_command('r.relief', input=rsurfName, output=reliefName,
+                        overwrite=overwrite)
+
+    gscript.run_command('r.shade', shade=reliefName,
+                        color=orthophotoName+'.rgb', output=shadedReliefName,
+                        overwrite=overwrite)
+
+    calc = ';'.join([
+        '$shadedRelief.red = ' +
+        'if(isnull($orthophoto.red), 0, r#$shadedRelief)',
+        '$shadedRelief.green = ' +
+        'if(isnull($orthophoto.green), 0, g#$shadedRelief)',
+        '$shadedRelief.blue = ' +
+        'if(isnull($orthophoto.blue), 0, b#$shadedRelief)',
+        '$shadedRelief.alpha = ' +
+        'if(isnull($orthophoto.alpha), 0, 255)'
+    ])
+    gscript.mapcalc(calc, shadedRelief=shadedReliefName,
+                    orthophoto=orthophotoName, overwrite=overwrite)
+
+    gscript.run_command('i.group', group=shadedReliefName+'.group',
+                        input=shadedReliefName+'.red,' +
+                        shadedReliefName+'.green,' +
+                        shadedReliefName+'.blue,' +
+                        shadedReliefName+'.alpha')
+
+    gscript.run_command('r.out.gdal', flags='cm',
+                        input=shadedReliefName+'.group',
+                        output=projectHome+'/odm_orthophoto/odm_relief.tif',
+                        format='GTiff', type='Byte',
+                        createopt='TILED=yes,COMPRESS=DEFLATE,PREDICTOR=2,' +
+                        'BLOCKXSIZE=512,BLOCKYSIZE=512',
+                        nodata=0, overwrite=overwrite)
+
+
+if __name__ == '__main__':
+    main()

modules/odm_georef/src/Georef.cpp

@@ -1500,6 +1500,7 @@ bool Georef::loadObjFile(std::string inputFile, pcl::TextureMesh &mesh)
     if (vt_idx != v_idx)
     {
         std::vector<Eigen::Vector2f> texcoordinates = std::vector<Eigen::Vector2f>(0);
+        texcoordinates.reserve(3*f_idx);
 
         for (size_t faceIndex = 0; faceIndex < f_idx; ++faceIndex)
         {

modules/odm_orthophoto/src/OdmOrthoPhoto.cpp

@@ -395,6 +395,12 @@ void OdmOrthoPhoto::createOrthoPhoto()
         std::vector<Eigen::Vector2f> textureCoordinates = std::vector<Eigen::Vector2f>(0);
 
         size_t vertexIndexCount = 0;
+        for(size_t t = 0; t < mesh.tex_polygons.size(); ++t)
+        {
+            vertexIndexCount += 3 * mesh.tex_polygons[t].size();
+        }
+        textureCoordinates.reserve(vertexIndexCount);
+
         for(size_t t = 0; t < mesh.tex_polygons.size(); ++t)
         {
 
@@ -419,17 +425,14 @@ void OdmOrthoPhoto::createOrthoPhoto()
                 meshCloudSplit->points.push_back(v1);
                 textureCoordinates.push_back(vt1);
                 mesh.tex_polygons[t][faceIndex].vertices[0] = vertexIndexCount;
-                ++vertexIndexCount;
 
                 meshCloudSplit->points.push_back(v2);
                 textureCoordinates.push_back(vt2);
                 mesh.tex_polygons[t][faceIndex].vertices[1] = vertexIndexCount;
-                ++vertexIndexCount;
 
                 meshCloudSplit->points.push_back(v3);
                 textureCoordinates.push_back(vt3);
                 mesh.tex_polygons[t][faceIndex].vertices[2] = vertexIndexCount;
-                ++vertexIndexCount;
             }
         }
 
@@ -450,6 +453,7 @@ void OdmOrthoPhoto::createOrthoPhoto()
 
     // Flatten texture coordinates.
     std::vector<Eigen::Vector2f> uvs;
+    uvs.reserve(mesh.tex_coordinates.size());
     for(size_t t = 0; t < mesh.tex_coordinates.size(); ++t)
     {
         uvs.insert(uvs.end(), mesh.tex_coordinates[t].begin(), mesh.tex_coordinates[t].end());
@@ -1240,6 +1244,7 @@ bool OdmOrthoPhoto::loadObjFile(std::string inputFile, pcl::TextureMesh &mesh)
     if (vt_idx != v_idx)
     {
         std::vector<Eigen::Vector2f> texcoordinates = std::vector<Eigen::Vector2f>(0);
+        texcoordinates.reserve(3*f_idx);
 
         for (size_t faceIndex = 0; faceIndex < f_idx; ++faceIndex)
         {

opendm/config.py

@@ -129,7 +129,7 @@ def config():
                         help=('The maximum number of processes to use in dense '
                               'reconstruction. Default: %(default)s'))
 
-    parser.add_argument('--use-opensfm-pointcloud',
+    parser.add_argument('--use-pmvs',
                         action='store_true',
                         default=False,
                         help='Use OpenSfM to compute the point cloud instead '
@@ -232,7 +232,7 @@ def config():
 
     parser.add_argument('--texturing-outlier-removal-type',
                         metavar='<string>',
-                        default='none',
+                        default='gauss_clamping',
                         help=('Type of photometric outlier removal method: ' 
                               '[none, gauss_damping, gauss_clamping]. Default: '  
                               '%(default)s'))

scripts/mvstex.py

@@ -85,7 +85,7 @@ class ODMMvsTexCell(ecto.Cell):
                 'keepUnseenFaces': keepUnseenFaces
             }
 
-            if args.use_opensfm_pointcloud:
+            if not args.use_pmvs:
                 kwargs['nvm_file'] = io.join_paths(tree.opensfm,
                                                    "reconstruction.nvm")
             else:

scripts/odm_app.py

@@ -112,7 +112,7 @@ class ODMApp(ecto.BlackBox):
                         self.args[:] >> self.opensfm['args'],
                         self.resize['photos'] >> self.opensfm['photos']]
 
-        if _p.args.use_opensfm_pointcloud:
+        if not _p.args.use_pmvs:
             # create odm mesh from opensfm point cloud
             connections += [self.tree[:] >> self.meshing['tree'],
                             self.args[:] >> self.meshing['args'],

scripts/odm_georeferencing.py

@@ -107,7 +107,7 @@ class ODMGeoreferencingCell(ecto.Cell):
                 'verbose': verbose
 
             }
-            if args.use_opensfm_pointcloud:
+            if not args.use_pmvs:
                 kwargs['pc'] = tree.opensfm_model
             else:
                 kwargs['pc'] = tree.pmvs_model

scripts/odm_meshing.py

@@ -62,7 +62,7 @@ class ODMeshingCell(ecto.Cell):
                 'solver': self.params.solver,
                 'verbose': verbose
             }
-            if args.use_opensfm_pointcloud:
+            if not args.use_pmvs:
                 kwargs['infile'] = tree.opensfm_model
             else:
                 kwargs['infile'] = tree.pmvs_model

scripts/opensfm.py

@@ -48,7 +48,7 @@ class ODMOpenSfMCell(ecto.Cell):
                      (args.rerun_from is not None and
                       'opensfm' in args.rerun_from)
 
-        if args.use_opensfm_pointcloud:
+        if not args.use_pmvs:
             output_file = tree.opensfm_model
         else:
             output_file = tree.opensfm_reconstruction
@@ -81,7 +81,7 @@ class ODMOpenSfMCell(ecto.Cell):
             # run OpenSfM reconstruction
             system.run('PYTHONPATH=%s %s/bin/run_all %s' %
                        (context.pyopencv_path, context.opensfm_path, tree.opensfm))
-            if args.use_opensfm_pointcloud:
+            if not args.use_pmvs:
                 system.run('PYTHONPATH=%s %s/bin/opensfm export_visualsfm %s' %
                            (context.pyopencv_path, context.opensfm_path, tree.opensfm))
                 system.run('PYTHONPATH=%s %s/bin/opensfm undistort %s' %
@@ -101,7 +101,7 @@ class ODMOpenSfMCell(ecto.Cell):
             log.ODM_WARNING('Found a valid Bundler file in: %s' %
                             tree.opensfm_reconstruction)
 
-        if not args.use_opensfm_pointcloud:
+        if args.use_pmvs:
             # check if reconstruction was exported to pmvs before
             if not io.file_exists(tree.pmvs_visdat) or rerun_cell:
                 # run PMVS converter