OpenDroneMap · pierotofy · Dec 14, 2022 · Dec 1, 2022 · Dec 2, 2022 · Dec 2, 2022
diff --git a/SuperBuild/CMakeLists.txt b/SuperBuild/CMakeLists.txt
@@ -159,6 +159,7 @@ SETUP_EXTERNAL_PROJECT(Hexer 1.4 ON)
 set(custom_libs OpenSfM
 				LASzip
 				PDAL
+				PDALPython
 				Untwine
 				Entwine
 				MvsTexturing

diff --git a/SuperBuild/cmake/External-OpenSfM.cmake b/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-------------

diff --git a/SuperBuild/cmake/External-PDAL.cmake b/SuperBuild/cmake/External-PDAL.cmake
@@ -16,7 +16,7 @@ ExternalProject_Add(${_proj_name}
   STAMP_DIR         ${_SB_BINARY_DIR}/stamp
   #--Download step--------------
   DOWNLOAD_DIR      ${SB_DOWNLOAD_DIR}
-  URL               https://github.com/PDAL/PDAL/archive/refs/tags/2.3RC1.zip
+  URL               https://github.com/PDAL/PDAL/archive/refs/tags/2.4.3.zip
   #--Update/Patch step----------
   UPDATE_COMMAND    ""
   #--Configure step-------------
@@ -60,3 +60,4 @@ ExternalProject_Add(${_proj_name}
   LOG_CONFIGURE     OFF
   LOG_BUILD         OFF
 )
+
diff --git a/SuperBuild/cmake/External-PDALPython.cmake b/SuperBuild/cmake/External-PDALPython.cmake
@@ -0,0 +1,36 @@
+set(_proj_name pdal-python)
+set(_SB_BINARY_DIR "${SB_BINARY_DIR}/${_proj_name}")
+
+if (WIN32)
+  set(PP_EXTRA_ARGS -DPYTHON3_EXECUTABLE=${PYTHON_EXE_PATH}
+                    -DPython3_NumPy_INCLUDE_DIRS=${PYTHON_HOME}/lib/site-packages/numpy/core/include)
+endif()
+
+ExternalProject_Add(${_proj_name}
+  DEPENDS           pdal
+  PREFIX            ${_SB_BINARY_DIR}
+  TMP_DIR           ${_SB_BINARY_DIR}/tmp
+  STAMP_DIR         ${_SB_BINARY_DIR}/stamp
+  #--Download step--------------
+  DOWNLOAD_DIR      ${SB_DOWNLOAD_DIR}
+  GIT_REPOSITORY    https://github.com/OpenDroneMap/pdal-python
+  GIT_TAG           main
+  #--Update/Patch step----------
+  UPDATE_COMMAND    ""
+  #--Configure step-------------
+  SOURCE_DIR        ${SB_SOURCE_DIR}/${_proj_name}
+  CMAKE_ARGS
+    -DPDAL_DIR=${SB_INSTALL_DIR}/lib/cmake/PDAL
+    -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
+    -DCMAKE_INSTALL_PREFIX:PATH=${SB_INSTALL_DIR}/lib/python3.8/dist-packages
+    ${WIN32_CMAKE_ARGS}
+    ${PP_EXTRA_ARGS}
+  #--Build step-----------------
+  BINARY_DIR        ${_SB_BINARY_DIR}
+  #--Install step---------------
+  INSTALL_DIR       ${SB_INSTALL_DIR}
+  #--Output logging-------------
+  LOG_DOWNLOAD      OFF
+  LOG_CONFIGURE     OFF
+  LOG_BUILD         OFF
+)
diff --git a/VERSION b/VERSION
@@ -1 +1 @@
-3.0.1
+3.0.2
diff --git a/opendm/align.py b/opendm/align.py
@@ -0,0 +1,147 @@
+import os
+import shutil
+import json
+import codem
+import dataclasses
+import pdal
+import numpy as np
+import rasterio
+from rasterio.crs import CRS
+from opendm.utils import double_quote
+from opendm import log
+from opendm import io
+from opendm import system
+from opendm.concurrency import get_max_memory
+
+def get_point_cloud_crs(file):
+    pipeline = pdal.Pipeline(json.dumps([ file ]))
+    metadata = pipeline.quickinfo
+
+    reader_metadata = [val for key, val in metadata.items() if "readers" in key]
+    crs = CRS.from_string(reader_metadata[0]["srs"]["horizontal"])
+    return str(crs)
+
+def get_raster_crs(file):
+    with rasterio.open(file, 'r') as f:
+        return str(f.crs)
+
+def reproject_point_cloud(file, out_srs):
+    out_file = io.related_file_path(file, postfix="_reprojected_tmp")
+    pipeline = pdal.Pipeline(json.dumps([ file, {
+        "type": "filters.reprojection",
+        "out_srs": out_srs
+    }, out_file]))
+    pipeline.execute()
+    return out_file
+
+def reproject_raster(file, out_srs):
+    out_file = io.related_file_path(file, postfix="_reprojected_tmp")
+    kwargs = {
+        'input': double_quote(file),
+        'output': double_quote(out_file),
+        'out_srs': out_srs,
+        'max_memory': get_max_memory()
+    }
+    system.run('gdalwarp '
+                '-t_srs {out_srs} '
+                '{input} '
+                '{output} '
+                '--config GDAL_CACHEMAX {max_memory}% '.format(**kwargs))
+    return out_file
+
+def compute_alignment_matrix(input_laz, align_file, stats_dir):
+    if os.path.exists(stats_dir):
+        shutil.rmtree(stats_dir)
+    os.mkdir(stats_dir)
+
+    # Check if we need to reproject align file
+    input_crs = get_point_cloud_crs(input_laz)
+    log.ODM_INFO("Input CRS: %s" % input_crs)
+
+    _, ext = os.path.splitext(align_file)
+    repr_func = None
+
+    if ext.lower() in [".tif"]:
+        align_crs = get_raster_crs(align_file)
+        repr_func = reproject_raster
+    elif ext.lower() in [".las", ".laz"]:
+        align_crs = get_point_cloud_crs(align_file)
+        repr_func = reproject_point_cloud
+    else:
+        log.ODM_WARNING("Unsupported alignment file: %s" % align_file)
+        return
+
+    to_delete = []
+
+    try:
+        log.ODM_INFO("Align CRS: %s" % align_crs)
+        if input_crs != align_crs:
+            # Reprojection needed
+            log.ODM_INFO("Reprojecting %s to %s" % (align_file, input_crs))
+            align_file = repr_func(align_file, input_crs)
+            to_delete.append(align_file)
+
+        conf = dataclasses.asdict(codem.CodemRunConfig(align_file, input_laz, OUTPUT_DIR=stats_dir))
+        fnd_obj, aoi_obj = codem.preprocess(conf)
+        fnd_obj.prep()
+        aoi_obj.prep()
+        log.ODM_INFO("Aligning reconstruction to %s" % align_file)
+        log.ODM_INFO("Coarse registration...")
+        dsm_reg = codem.coarse_registration(fnd_obj, aoi_obj, conf)
+        log.ODM_INFO("Fine registration...")
+        icp_reg = codem.fine_registration(fnd_obj, aoi_obj, dsm_reg, conf)
+
+        app_reg = codem.registration.ApplyRegistration(
+            fnd_obj,
+            aoi_obj,
+            icp_reg.registration_parameters,
+            icp_reg.residual_vectors,
+            icp_reg.residual_origins,
+            conf,
+            None,
+        )
+
+        reg = app_reg.get_registration_transformation()
+
+        # Write JSON to stats folder
+        with open(os.path.join(stats_dir, "registration.json"), 'w') as f:
+            del dsm_reg.registration_parameters['matrix']
+            del icp_reg.registration_parameters['matrix']
+
+            f.write(json.dumps({
+                'coarse':  dsm_reg.registration_parameters,
+                'fine': icp_reg.registration_parameters,
+            }, indent=4))
+
+        matrix = np.fromstring(reg['matrix'], dtype=float, sep=' ').reshape((4, 4))
+        return matrix
+    finally:
+        for f in to_delete:
+            if os.path.isfile(f):
+                os.unlink(f)
+
+def transform_point_cloud(input_laz, a_matrix, output_laz):
+    pipe = [
+        input_laz,
+        {
+            'type': 'filters.transformation',
+            'matrix': " ".join(list(map(str, a_matrix.flatten()))),
+        },
+        output_laz,
+    ]
+    p = pdal.Pipeline(json.dumps(pipe))
+    p.execute()
+
+def transform_obj(input_obj, a_matrix, geo_offset, output_obj):
+    g_off = np.array([geo_offset[0], geo_offset[1], 0, 0])
+
+    with open(input_obj, 'r') as fin:
+        with open(output_obj, 'w') as fout:
+            lines = fin.readlines()
+            for line in lines:
+                if line.startswith("v "):
+                    v = np.fromstring(line.strip()[2:] + " 1",  sep=' ', dtype=float)
+                    vt = (a_matrix.dot((v + g_off)) - g_off)[:3]
+                    fout.write("v " + " ".join(map(str, list(vt))) + '\n')
+                else:
+                    fout.write(line)
diff --git a/opendm/config.py b/opendm/config.py
@@ -470,6 +470,14 @@ def config(argv=None, parser=None):
                               'EPSG:<code> or <+proj definition>\n'
                               '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',
+                    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. Experimental. '
+                            'Default: %(default)s'))
 
     parser.add_argument('--use-exif',
                         action=StoreTrue,

diff --git a/opendm/context.py b/opendm/context.py
@@ -14,7 +14,7 @@
     'install/lib/python3.9/dist-packages',
     'install/lib/python3.8/dist-packages',
     'install/lib/python3/dist-packages',
-    'install/bin/opensfm'
+    'install/bin/opensfm',
 ]]
 for p in python_packages_paths:
     sys.path.append(p)

diff --git a/opendm/io.py b/opendm/io.py
@@ -85,3 +85,7 @@ def path_or_json_string_to_dict(string):
             raise ValueError("{0} is not a valid JSON file.".format(string))
     else:
         raise ValueError("{0} is not a valid JSON file or string.".format(string))
+
+def touch(file):
+    with open(file, 'w') as fout:
+        fout.write("Done!\n")
diff --git a/opendm/shots.py b/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])
 

diff --git a/opendm/types.py b/opendm/types.py
@@ -241,7 +241,7 @@ def utm_offset(self):
         return (self.utm_east_offset, self.utm_north_offset)
 
 class ODM_Tree(object):
-    def __init__(self, root_path, gcp_file = None, geo_file = None):
+    def __init__(self, root_path, gcp_file = None, geo_file = None, align_file = None):
         # root path to the project
         self.root_path = io.absolute_path_file(root_path)
         self.input_images = os.path.join(self.root_path, 'images')
@@ -296,6 +296,7 @@ def __init__(self, root_path, gcp_file = None, geo_file = None):
         self.odm_georeferencing_gcp = gcp_file or io.find('gcp_list.txt', self.root_path)
         self.odm_georeferencing_gcp_utm = os.path.join(self.odm_georeferencing, 'gcp_list_utm.txt')
         self.odm_geo_file = geo_file or io.find('geo.txt', self.root_path)
+        self.odm_align_file = align_file or io.find('align.laz', self.root_path) or io.find('align.las', self.root_path) or io.find('align.tif', self.root_path)
 
         self.odm_georeferencing_proj = 'proj.txt'
         self.odm_georeferencing_model_txt_geo = os.path.join(
@@ -306,6 +307,9 @@ def __init__(self, root_path, gcp_file = None, geo_file = None):
             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')

diff --git a/opendm/utils.py b/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))
diff --git a/requirements.txt b/requirements.txt
@@ -31,3 +31,6 @@ xmltodict==0.12.0
 fpdf2==2.4.6
 Shapely==1.7.1
 >
+codem==0.24.0
+trimesh==3.17.1
+pandas==1.5.2
diff --git a/stages/dataset.py b/stages/dataset.py
@@ -46,7 +46,7 @@ class Empty:
 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)
+        tree = types.ODM_Tree(args.project_path, args.gcp, args.geo, args.align)
         outputs['tree'] = tree
 
         if io.file_exists(tree.benchmarking):

diff --git a/stages/mvstex.py b/stages/mvstex.py
@@ -68,11 +68,15 @@ def add_run(nvm_file, primary=True, band=None):
                 system.mkdir_p(r['out_dir'])
 
             odm_textured_model_obj = os.path.join(r['out_dir'], tree.odm_textured_model_obj)
+            unaligned_obj = io.related_file_path(odm_textured_model_obj, postfix="_unaligned")
 
             if not io.file_exists(odm_textured_model_obj) or self.rerun():
                 log.ODM_INFO('Writing MVS Textured file in: %s'
                               % odm_textured_model_obj)
 
+                if os.path.isfile(unaligned_obj):
+                    os.unlink(unaligned_obj)
+
                 # Format arguments to fit Mvs-Texturing app
                 skipGlobalSeamLeveling = ""
                 skipLocalSeamLeveling = ""
-Original file line number
+Diff line change
@@ Expand Up / @@ -159,6 +159,7 @@ SETUP_EXTERNAL_PROJECT(Hexer 1.4 ON) @@
     set(custom_libs OpenSfM
     				LASzip
     				PDAL
+    				PDALPython
     				Untwine
     				Entwine
     				MvsTexturing
@@ Expand Down @@