Use OpenSfM GCP stats, keep copy of topocentric rec, export GCPs to GeoJSON

pull/1351/head
Piero Toffanin 2021-09-24 15:06:40 +00:00
rodzic 0589483b9b
commit 5259491165
6 zmienionych plików z 58 dodań i 50 usunięć

Wyświetl plik

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

Wyświetl plik

@ -1 +1 @@
2.6.3
2.6.4

Wyświetl plik

@ -23,7 +23,7 @@ from opensfm import report
from opendm.multispectral import get_photos_by_band
from opendm.gpu import has_gpus
from opensfm import multiview
from opensfm.actions.export_geocoords import _get_transformation
from opensfm.actions.export_geocoords import _transform
class OSFMContext:
def __init__(self, opensfm_project_path):
@ -458,47 +458,33 @@ class OSFMContext:
"""
Load ground control point information.
"""
ds = DataSet(self.opensfm_project_path)
gcps = ds.load_ground_control_points()
if not gcps:
gcp_stats_file = self.path("stats", "ground_control_points.json")
if not io.file_exists(gcp_stats_file):
return []
reconstructions = ds.load_reconstruction()
reference = ds.load_reference()
gcps_stats = {}
try:
with open(gcp_stats_file) as f:
gcps_stats = json.loads(f.read())
except:
log.ODM_INFO("Cannot parse %s" % gcp_stats_file)
if not gcps_stats:
return []
ds = DataSet(self.opensfm_project_path)
reference = ds.load_reference()
projection = pyproj.Proj(proj4)
t = _get_transformation(reference, projection, (0, 0))
A, b = t[:3, :3], t[:3, 3]
result = []
for gcp in gcps:
if not gcp.coordinates.has_value:
continue
triangulated = None
for rec in reconstructions:
triangulated = multiview.triangulate_gcp(gcp, rec.shots, 1.0, 0.1)
if triangulated is None:
continue
else:
break
if triangulated is None:
continue
triangulated_topocentric = np.dot(A.T, triangulated)
coordinates_topocentric = np.array(gcp.coordinates.value)
coordinates = np.dot(A, coordinates_topocentric) + b
triangulated = triangulated + b
for gcp in gcps_stats:
geocoords = _transform(gcp['coordinates'], reference, projection)
result.append({
'id': gcp.id,
'observations': [obs.shot_id for obs in gcp.observations],
'triangulated': triangulated,
'coordinates': coordinates,
'error': np.abs(triangulated_topocentric - coordinates_topocentric)
'id': gcp['id'],
'observations': gcp['observations'],
'coordinates': geocoords,
'error': gcp['error']
})
return result

Wyświetl plik

@ -235,6 +235,7 @@ class ODM_Tree(object):
self.opensfm_reconstruction = os.path.join(self.opensfm, 'reconstruction.json')
self.opensfm_reconstruction_nvm = os.path.join(self.opensfm, 'undistorted/reconstruction.nvm')
self.opensfm_geocoords_reconstruction = os.path.join(self.opensfm, 'reconstruction.geocoords.json')
self.opensfm_topocentric_reconstruction = os.path.join(self.opensfm, 'reconstruction.topocentric.json')
# OpenMVS
self.openmvs_model = os.path.join(self.openmvs, 'scene_dense_dense_filtered.ply')

Wyświetl plik

@ -3,18 +3,22 @@ import struct
import pipes
import fiona
import fiona.crs
import json
from collections import OrderedDict
from pyproj import CRS
from opendm import io
from opendm import log
from opendm import types
from opendm import system
from opendm import context
from opendm import location
from opendm.cropper import Cropper
from opendm import point_cloud
from opendm.multispectral import get_primary_band_name
from opendm.osfm import OSFMContext
class ODMGeoreferencingStage(types.ODM_Stage):
def process(self, args, outputs):
tree = outputs['tree']
@ -24,6 +28,7 @@ class ODMGeoreferencingStage(types.ODM_Stage):
gcp_export_file = tree.path("odm_georeferencing", "ground_control_points.gpkg")
gcp_gml_export_file = tree.path("odm_georeferencing", "ground_control_points.gml")
gcp_geojson_export_file = tree.path("odm_georeferencing", "ground_control_points.geojson")
if reconstruction.has_gcp() and (not io.file_exists(gcp_export_file) or self.rerun()):
octx = OSFMContext(tree.opensfm)
@ -36,9 +41,6 @@ class ODMGeoreferencingStage(types.ODM_Stage):
('id', 'str'),
('observations_count', 'int'),
('observations_list', 'str'),
('triangulated_x', 'float'),
('triangulated_y', 'float'),
('triangulated_z', 'float'),
('error_x', 'float'),
('error_y', 'float'),
('error_z', 'float'),
@ -58,10 +60,7 @@ class ODMGeoreferencingStage(types.ODM_Stage):
'properties': OrderedDict([
('id', gcp['id']),
('observations_count', len(gcp['observations'])),
('observations_list', ",".join(gcp['observations'])),
('triangulated_x', gcp['triangulated'][0]),
('triangulated_y', gcp['triangulated'][1]),
('triangulated_z', gcp['triangulated'][2]),
('observations_list', ",".join([obs['shot_id'] for obs in gcp['observations']])),
('error_x', gcp['error'][0]),
('error_y', gcp['error'][1]),
('error_z', gcp['error'][2]),
@ -73,10 +72,35 @@ class ODMGeoreferencingStage(types.ODM_Stage):
system.run('ogr2ogr -of GML "{}" "{}"'.format(gcp_gml_export_file, gcp_export_file))
except Exception as e:
log.ODM_WARNING("Cannot generate ground control points GML file: %s" % str(e))
# Write GeoJSON
geojson = {
'type': 'FeatureCollection',
'features': []
}
from_srs = CRS.from_proj4(reconstruction.georef.proj4())
to_srs = CRS.from_epsg(4326)
transformer = location.transformer(from_srs, to_srs)
for gcp in gcps:
properties = gcp.copy()
del properties['coordinates']
geojson['features'].append({
'geometry': {
'type': 'Point',
'coordinates': transformer.TransformPoint(*gcp['coordinates']),
},
'properties': properties
})
with open(gcp_geojson_export_file, 'w') as f:
f.write(json.dumps(geojson, indent=4))
else:
log.ODM_WARNING("GCPs could not be loaded for writing to %s" % gcp_export_file)
if not io.file_exists(tree.odm_georeferencing_model_laz) or self.rerun():
cmd = ('pdal translate -i "%s" -o \"%s\"' % (tree.filtered_point_cloud, tree.odm_georeferencing_model_laz))
stages = ["ferry"]

Wyświetl plik

@ -68,14 +68,11 @@ class ODMOpenSfMStage(types.ODM_Stage):
self.update_progress(75)
# We now switch to a geographic CRS
geocoords_flag_file = octx.path("exported_geocoords.txt")
if reconstruction.is_georeferenced() and (not io.file_exists(geocoords_flag_file) or self.rerun()):
if reconstruction.is_georeferenced() and (not io.file_exists(tree.opensfm_topocentric_reconstruction) or self.rerun()):
octx.run('export_geocoords --reconstruction --proj "%s" --offset-x %s --offset-y %s' %
(reconstruction.georef.proj4(), reconstruction.georef.utm_east_offset, reconstruction.georef.utm_north_offset))
# Destructive
shutil.move(tree.opensfm_reconstruction, tree.opensfm_topocentric_reconstruction)
shutil.move(tree.opensfm_geocoords_reconstruction, tree.opensfm_reconstruction)
octx.touch(geocoords_flag_file)
else:
log.ODM_WARNING("Will skip exporting %s" % tree.opensfm_geocoords_reconstruction)