Transform shots.geojson, bump opensfm, bump version

SuperBuild/cmake/External-OpenSfM.cmake

@@ -25,7 +25,7 @@ ExternalProject_Add(${_proj_name}
   #--Download step--------------
   DOWNLOAD_DIR      ${SB_DOWNLOAD_DIR}
   GIT_REPOSITORY    https://github.com/OpenDroneMap/OpenSfM/
-  GIT_TAG           292
+  GIT_TAG           302
   #--Update/Patch step----------
   UPDATE_COMMAND    git submodule update --init --recursive
   #--Configure step-------------

VERSION

@@ -1 +1 @@
-3.0.1
+3.0.2

opendm/config.py

@@ -471,12 +471,12 @@ def config(argv=None, parser=None):
                               'image_name geo_x geo_y geo_z [omega (degrees)] [phi (degrees)] [kappa (degrees)] [horz accuracy (meters)] [vert accuracy (meters)]\n'
                               'Default: %(default)s'))
     
-    parser.add_argument('--align-to',
+    parser.add_argument('--align',
                     metavar='<path string>',
                     action=StoreValue,
                     default=None,
                     help=('Path to a GeoTIFF DEM or a LAS/LAZ point cloud '
-                            'that the reconstruction outputs should be automatically aligned to. '
+                            'that the reconstruction outputs should be automatically aligned to. Experimental. '
                             'Default: %(default)s'))
 
     parser.add_argument('--use-exif',

opendm/shots.py

@@ -23,7 +23,7 @@ def get_origin(shot):
     """The origin of the pose in world coordinates."""
     return -get_rotation_matrix(np.array(shot['rotation'])).T.dot(np.array(shot['translation']))
 
-def get_geojson_shots_from_opensfm(reconstruction_file, utm_srs=None, utm_offset=None, pseudo_geotiff=None):
+def get_geojson_shots_from_opensfm(reconstruction_file, utm_srs=None, utm_offset=None, pseudo_geotiff=None, a_matrix=None):
     """
     Extract shots from OpenSfM's reconstruction.json
     """
@@ -92,6 +92,11 @@ def get_geojson_shots_from_opensfm(reconstruction_file, utm_srs=None, utm_offset
                         utm_coords = [origin[0] + utm_offset[0],
                                        origin[1] + utm_offset[1],
                                        origin[2]]
+
+                        if a_matrix is not None:
+                            rotation = list(np.array(rotation).dot(a_matrix[:3,:3]))
+                            utm_coords = list(a_matrix.dot(np.hstack((np.array(utm_coords), 1)))[:-1])
+                        
                         translation = utm_coords
                         trans_coords = crstrans.TransformPoint(utm_coords[0], utm_coords[1], utm_coords[2])
 

opendm/types.py

@@ -307,6 +307,9 @@ class ODM_Tree(object):
             self.odm_georeferencing, 'odm_georeferenced_model.laz')
         self.odm_georeferencing_model_las = os.path.join(
             self.odm_georeferencing, 'odm_georeferenced_model.las')
+        self.odm_georeferencing_alignment_matrix = os.path.join(
+            self.odm_georeferencing, 'alignment_matrix.json'
+        )
 
         # odm_orthophoto
         self.odm_orthophoto_render = os.path.join(self.odm_orthophoto, 'odm_orthophoto_render.tif')

opendm/utils.py

@@ -1,9 +1,18 @@
 import os, shutil
+import numpy as np
+import json
 from opendm import log
 from opendm.photo import find_largest_photo_dims
 from osgeo import gdal
 from opendm.loghelpers import double_quote
 
+class NumpyEncoder(json.JSONEncoder):
+    def default(self, obj):
+        if isinstance(obj, np.ndarray):
+            return obj.tolist()
+        return json.JSONEncoder.default(self, obj)
+
+
 def get_depthmap_resolution(args, photos):
     max_dims = find_largest_photo_dims(photos)
     min_dim = 320 # Never go lower than this
@@ -98,4 +107,10 @@ def rm_r(path):
         elif os.path.exists(path):
             os.remove(path)
     except:
-        log.ODM_WARNING("Cannot remove %s" % path)
+        log.ODM_WARNING("Cannot remove %s" % path)
+
+def np_to_json(arr):
+    return json.dumps(arr, cls=NumpyEncoder)
+
+def np_from_json(json_dump):
+    return np.asarray(json.loads(json_dump))

stages/dataset.py

@@ -46,7 +46,7 @@ def load_images_database(database_file):
 class ODMLoadDatasetStage(types.ODM_Stage):
     def process(self, args, outputs):
         outputs['start_time'] = system.now_raw()
-        tree = types.ODM_Tree(args.project_path, args.gcp, args.geo, args.align_to)
+        tree = types.ODM_Tree(args.project_path, args.gcp, args.geo, args.align)
         outputs['tree'] = tree
 
         if io.file_exists(tree.benchmarking):

stages/odm_georeferencing.py

@@ -1,4 +1,5 @@
 import os
+import shutil
 import struct
 import pipes
 import fiona
@@ -19,6 +20,7 @@ from opendm.multispectral import get_primary_band_name
 from opendm.osfm import OSFMContext
 from opendm.boundary import as_polygon, export_to_bounds_files
 from opendm.align import compute_alignment_matrix, transform_point_cloud, transform_obj
+from opendm.utils import np_to_json
 
 class ODMGeoreferencingStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -172,11 +174,17 @@ class ODMGeoreferencingStage(types.ODM_Stage):
                 system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))
 
 
-            if tree.odm_align_file is not None:
-                alignment_file = os.path.join(tree.odm_georeferencing, 'alignment_done.txt')
+            stats_dir = tree.path("opensfm", "stats", "codem")
+            if os.path.exists(stats_dir) and self.rerun():
+                shutil.rmtree(stats_dir)
+
+            if tree.odm_align_file is not None:
+                alignment_file_exists = io.file_exists(tree.odm_georeferencing_alignment_matrix)
+
+                if not alignment_file_exists or self.rerun():
+                    if alignment_file_exists:
+                        os.unlink(tree.odm_georeferencing_alignment_matrix)
 
-                if not io.file_exists(alignment_file) or self.rerun():
-                    stats_dir = tree.path("opensfm", "stats", "codem")
                     a_matrix = compute_alignment_matrix(tree.odm_georeferencing_model_laz, tree.odm_align_file, stats_dir)
                     if a_matrix is not None:
                         log.ODM_INFO("Alignment matrix: %s" % a_matrix)
@@ -207,12 +215,15 @@ class ODMGeoreferencingStage(types.ODM_Stage):
                                 except Exception as e:
                                     log.ODM_WARNING("Cannot transform textured model: %s" % str(e))
                                     os.rename(unaligned_obj, obj)
+                        
+                        with open(tree.odm_georeferencing_alignment_matrix, "w") as f:
+                            f.write(np_to_json(a_matrix))
                     else:
                         log.ODM_WARNING("Alignment to %s will be skipped." % tree.odm_align_file)
-                    
-                    io.touch(alignment_file)
                 else:
                     log.ODM_WARNING("Already computed alignment")
+            elif io.file_exists(tree.odm_georeferencing_alignment_matrix):
+                os.unlink(tree.odm_georeferencing_alignment_matrix)
 
             point_cloud.post_point_cloud_steps(args, tree, self.rerun())
         else:

stages/odm_report.py

@@ -14,7 +14,7 @@ from opendm.point_cloud import export_info_json
 from opendm.cropper import Cropper
 from opendm.orthophoto import get_orthophoto_vars, get_max_memory, generate_png
 from opendm.tiles.tiler import generate_colored_hillshade
-from opendm.utils import get_raster_stats
+from opendm.utils import get_raster_stats, np_from_json
 
 def hms(seconds):
     h = seconds // 3600
@@ -49,7 +49,14 @@ class ODMReport(types.ODM_Stage):
         if not io.file_exists(shots_geojson) or self.rerun():
             # Extract geographical camera shots
             if reconstruction.is_georeferenced():
-                shots = get_geojson_shots_from_opensfm(tree.opensfm_reconstruction, utm_srs=reconstruction.get_proj_srs(), utm_offset=reconstruction.georef.utm_offset())    
+                # Check if alignment has been performed (we need to transform our shots if so)
+                a_matrix = None
+                if io.file_exists(tree.odm_georeferencing_alignment_matrix):
+                    with open(tree.odm_georeferencing_alignment_matrix, 'r') as f:
+                        a_matrix = np_from_json(f.read())
+                        log.ODM_INFO("Aligning shots to %s" % a_matrix)
+
+                shots = get_geojson_shots_from_opensfm(tree.opensfm_reconstruction, utm_srs=reconstruction.get_proj_srs(), utm_offset=reconstruction.georef.utm_offset(), a_matrix=a_matrix)
             else:
                 # Pseudo geo
                 shots = get_geojson_shots_from_opensfm(tree.opensfm_reconstruction, pseudo_geotiff=tree.odm_orthophoto_tif)
@@ -73,6 +80,7 @@ class ODMReport(types.ODM_Stage):
         odm_stats_json = os.path.join(tree.odm_report, "stats.json")
         octx = OSFMContext(tree.opensfm)
         osfm_stats_json = octx.path("stats", "stats.json")
+        codem_stats_json = octx.path("stats", "codem", "registration.json")
         odm_stats = None
         point_cloud_file = None
         views_dimension = None
@@ -111,6 +119,11 @@ class ODMReport(types.ODM_Stage):
                     'average_gsd': gsd.opensfm_reconstruction_average_gsd(octx.recon_file(), use_all_shots=reconstruction.has_gcp()),
                 }
 
+                # Add CODEM stats
+                if os.path.exists(codem_stats_json):
+                    with open(codem_stats_json, 'r') as f:
+                        odm_stats['align'] = json.loads(f.read())
+
                 with open(odm_stats_json, 'w') as f:
                     f.write(json.dumps(odm_stats))
             else: