pull/1851/merge
s.zhong 2025-04-09 11:08:31 +07:00 zatwierdzone przez GitHub
commit 9659bd702b
Nie znaleziono w bazie danych klucza dla tego podpisu
ID klucza GPG: B5690EEEBB952194
8 zmienionych plików z 186 dodań i 51 usunięć

56
opendm/georef.py 100644
Wyświetl plik

@ -0,0 +1,56 @@
import numpy as np
from numpy import ndarray
from typing import Tuple
from pyproj import Proj
from opensfm.geo import TopocentricConverter
def topocentric_to_georef(
reflat: float,
reflon: float,
refalt: float,
a_srs: str,
x: ndarray,
y: ndarray,
z: ndarray,
x_offset: float = 0,
y_offset: float = 0,
) -> Tuple[ndarray, ndarray, ndarray]:
reference = TopocentricConverter(reflat, reflon, refalt)
projection = Proj(a_srs)
lat, lon, alt = reference.to_lla(x, y, z)
easting, northing = projection(lon, lat)
return easting - x_offset, northing - y_offset, alt
class TopocentricToProj:
def __init__(self, reflat:float, reflon:float, refalt:float, a_srs:str):
self.reference = TopocentricConverter(reflat, reflon, refalt)
self.projection = Proj(a_srs)
def convert_array(self, arr:ndarray, offset_x:float=0, offset_y:float=0):
x, y, z = arr['X'], arr['Y'], arr['Z']
easting, northing, alt = self.convert_points(x, y, z, offset_x, offset_y)
arr['X'] = easting
arr['Y'] = northing
arr['Z'] = alt
return arr
def convert_points(self, x:ndarray, y:ndarray, z:ndarray, offset_x:float=0, offset_y:float=0):
lat, lon, alt = self.reference.to_lla(x, y, z)
easting, northing = self.projection(lon, lat)
return easting - offset_x, northing - offset_y, alt
def convert_obj(self, input_obj:str, output_obj:str, offset_x:float=0, offset_y:float=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 = self.convert_points(v[0], v[1], v[2], offset_x, offset_y)
fout.write("v " + " ".join(map(str, list(vt))) + '\n')
else:
fout.write(line)

Wyświetl plik

@ -629,9 +629,12 @@ class OSFMContext:
return result
def reference(self):
ds = DataSet(self.opensfm_project_path)
return ds.load_reference()
def name(self):
return os.path.basename(os.path.abspath(self.path("..")))
return os.path.basename(os.path.abspath(self.path("..")))
def get_submodel_argv(args, submodels_path = None, submodel_name = None):
"""

Wyświetl plik

@ -361,16 +361,20 @@ class ODM_Tree(object):
self.openmvs_model = os.path.join(self.openmvs, 'scene_dense_dense_filtered.ply')
# filter points
self.filtered_point_cloud_topo = os.path.join(self.odm_filterpoints, "point_cloud_topocentric.ply")
self.filtered_point_cloud = os.path.join(self.odm_filterpoints, "point_cloud.ply")
self.filtered_point_cloud_stats = os.path.join(self.odm_filterpoints, "point_cloud_stats.json")
# odm_meshing
self.odm_mesh_topo = os.path.join(self.odm_meshing, 'odm_mesh_topocentric.ply')
self.odm_mesh = os.path.join(self.odm_meshing, 'odm_mesh.ply')
self.odm_meshing_log = os.path.join(self.odm_meshing, 'odm_meshing_log.txt')
self.odm_25dmesh_topo = os.path.join(self.odm_meshing, 'odm_25dmesh_topocentric.ply')
self.odm_25dmesh = os.path.join(self.odm_meshing, 'odm_25dmesh.ply')
self.odm_25dmeshing_log = os.path.join(self.odm_meshing, 'odm_25dmeshing_log.txt')
# texturing
self.odm_textured_model_obj_topo = 'odm_textured_model_topocentric.obj'
self.odm_textured_model_obj = 'odm_textured_model_geo.obj'
self.odm_textured_model_glb = 'odm_textured_model_geo.glb'

Wyświetl plik

@ -33,7 +33,7 @@ class ODMMvsTexStage(types.ODM_Stage):
if not args.skip_3dmodel and (primary or args.use_3dmesh):
nonloc.runs += [{
'out_dir': os.path.join(tree.odm_texturing, subdir),
'model': tree.odm_mesh,
'model': tree.odm_mesh_topo,
'nadir': False,
'primary': primary,
'nvm_file': nvm_file,
@ -43,7 +43,7 @@ class ODMMvsTexStage(types.ODM_Stage):
if not args.use_3dmesh:
nonloc.runs += [{
'out_dir': os.path.join(tree.odm_25dtexturing, subdir),
'model': tree.odm_25dmesh,
'model': tree.odm_25dmesh_topo,
'nadir': True,
'primary': primary,
'nvm_file': nvm_file,
@ -69,9 +69,8 @@ class ODMMvsTexStage(types.ODM_Stage):
if not io.dir_exists(r['out_dir']):
system.mkdir_p(r['out_dir'])
odm_textured_model_obj = os.path.join(r['out_dir'], tree.odm_textured_model_obj)
odm_textured_model_obj = os.path.join(r['out_dir'], tree.odm_textured_model_obj_topo)
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)
@ -91,10 +90,12 @@ class ODMMvsTexStage(types.ODM_Stage):
if (r['nadir']):
nadir = '--nadir_mode'
# mvstex definitions
# mtl and texture files would be the same between topo and proj geomodel, so create with the final name
kwargs = {
'bin': context.mvstex_path,
'out_dir': os.path.join(r['out_dir'], "odm_textured_model_geo"),
'out_dir': os.path.join(r['out_dir'], 'odm_textured_model_geo'),
'model': r['model'],
'dataTerm': 'gmi',
'outlierRemovalType': 'gauss_clamping',
@ -125,6 +126,9 @@ class ODMMvsTexStage(types.ODM_Stage):
'{nadirMode} '
'{labelingFile} '
'{maxTextureSize} '.format(**kwargs))
# update the obj file name to topo for further conversion
shutil.move(os.path.join(r['out_dir'], tree.odm_textured_model_obj), odm_textured_model_obj)
if r['primary'] and (not r['nadir'] or args.skip_3dmodel):
# GlTF?
@ -147,7 +151,7 @@ class ODMMvsTexStage(types.ODM_Stage):
shutil.rmtree(packed_dir)
try:
obj_pack(os.path.join(r['out_dir'], tree.odm_textured_model_obj), packed_dir, _info=log.ODM_INFO)
obj_pack(os.path.join(r['out_dir'], tree.odm_textured_model_obj_topo), packed_dir, _info=log.ODM_INFO)
# Move packed/* into texturing folder
system.delete_files(r['out_dir'], (".vec", ))

Wyświetl plik

@ -19,7 +19,7 @@ class ODMFilterPoints(types.ODM_Stage):
inputPointCloud = ""
# check if reconstruction was done before
if not io.file_exists(tree.filtered_point_cloud) or self.rerun():
if not io.file_exists(tree.filtered_point_cloud_topo) or self.rerun():
if args.fast_orthophoto:
inputPointCloud = os.path.join(tree.opensfm, 'reconstruction.ply')
else:
@ -49,14 +49,14 @@ class ODMFilterPoints(types.ODM_Stage):
else:
log.ODM_WARNING("Not a georeferenced reconstruction, will ignore --auto-boundary")
point_cloud.filter(inputPointCloud, tree.filtered_point_cloud, tree.filtered_point_cloud_stats,
point_cloud.filter(inputPointCloud, tree.filtered_point_cloud_topo, tree.filtered_point_cloud_stats,
standard_deviation=args.pc_filter,
sample_radius=args.pc_sample,
boundary=boundary_offset(outputs.get('boundary'), reconstruction.get_proj_offset()),
max_concurrency=args.max_concurrency)
# Quick check
info = point_cloud.ply_info(tree.filtered_point_cloud)
info = point_cloud.ply_info(tree.filtered_point_cloud_topo)
if info["vertex_count"] == 0:
extra_msg = ''
if 'boundary' in outputs:
@ -64,7 +64,7 @@ class ODMFilterPoints(types.ODM_Stage):
raise system.ExitException("Uh oh! We ended up with an empty point cloud. This means that the reconstruction did not succeed. Have you followed best practices for data acquisition? See https://docs.opendronemap.org/flying/%s" % extra_msg)
else:
log.ODM_WARNING('Found a valid point cloud file in: %s' %
tree.filtered_point_cloud)
tree.filtered_point_cloud_topo)
if args.optimize_disk_space and inputPointCloud:
if os.path.isfile(inputPointCloud):

Wyświetl plik

@ -9,6 +9,7 @@ import zipfile
import math
from collections import OrderedDict
from pyproj import CRS
import pdal
from opendm import io
from opendm import log
@ -23,6 +24,7 @@ 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
from opendm.georef import TopocentricToProj
class ODMGeoreferencingStage(types.ODM_Stage):
def process(self, args, outputs):
@ -113,19 +115,72 @@ class ODMGeoreferencingStage(types.ODM_Stage):
else:
log.ODM_WARNING("GCPs could not be loaded for writing to %s" % gcp_export_file)
if reconstruction.is_georeferenced():
# prepare pipeline stage for topocentric to georeferenced conversion
octx = OSFMContext(tree.opensfm)
reference = octx.reference()
converter = TopocentricToProj(reference.lat, reference.lon, reference.alt, reconstruction.georef.proj4())
if not io.file_exists(tree.filtered_point_cloud) or self.rerun():
log.ODM_INFO("Georeferecing filtered point cloud")
if reconstruction.is_georeferenced():
pipeline = pdal.Reader.ply(tree.filtered_point_cloud_topo).pipeline()
pipeline.execute()
arr = pipeline.arrays[0]
arr = converter.convert_array(
arr,
reconstruction.georef.utm_east_offset,
reconstruction.georef.utm_north_offset
)
pipeline = pdal.Writer.ply(
filename = tree.filtered_point_cloud,
storage_mode = "little endian",
).pipeline(arr)
pipeline.execute()
else:
shutil.copy(tree.filtered_point_cloud_topo, tree.filtered_point_cloud)
def georefernce_textured_model(obj_in, obj_out):
log.ODM_INFO("Georeferecing textured model %s" % obj_in)
if not io.file_exists(obj_out) or self.rerun():
if reconstruction.is_georeferenced():
converter.convert_obj(
obj_in,
obj_out,
reconstruction.georef.utm_east_offset,
reconstruction.georef.utm_north_offset
)
else:
shutil.copy(obj_in, obj_out)
#TODO: maybe parallelize this
#TODO: gltf export? Should problably move the exporting process after this
for texturing in [tree.odm_texturing, tree.odm_25dtexturing]:
if reconstruction.multi_camera:
primary = get_primary_band_name(reconstruction.multi_camera, args.primary_band)
for band in reconstruction.multi_camera:
subdir = "" if band['name'] == primary else band['name'].lower()
obj_in = os.path.join(texturing, subdir, tree.odm_textured_model_obj_topo)
obj_out = os.path.join(texturing, subdir, tree.odm_textured_model_obj)
georefernce_textured_model(obj_in, obj_out)
else:
obj_in = os.path.join(texturing, tree.odm_textured_model_obj_topo)
obj_out = os.path.join(texturing, tree.odm_textured_model_obj)
transform_textured_model(obj_in, obj_out)
if not io.file_exists(tree.odm_georeferencing_model_laz) or self.rerun():
cmd = f'pdal translate -i "{tree.filtered_point_cloud}" -o \"{tree.odm_georeferencing_model_laz}\"'
stages = ["ferry"]
params = [
'--filters.ferry.dimensions="views => UserData"'
]
pipeline = pdal.Pipeline()
pipeline |= pdal.Reader.ply(tree.filtered_point_cloud)
pipeline |= pdal.Filter.ferry(dimensions="views => UserData")
if reconstruction.is_georeferenced():
log.ODM_INFO("Georeferencing point cloud")
stages.append("transformation")
utmoffset = reconstruction.georef.utm_offset()
pipeline |= pdal.Filter.transformation(
matrix=f"1 0 0 {utmoffset[0]} 0 1 0 {utmoffset[1]} 0 0 1 0 0 0 0 1"
)
# Establish appropriate las scale for export
las_scale = 0.001
@ -148,27 +203,37 @@ class ODMGeoreferencingStage(types.ODM_Stage):
else:
log.ODM_INFO("No point_cloud_stats.json found. Using default las scale: %s" % las_scale)
params += [
f'--filters.transformation.matrix="1 0 0 {utmoffset[0]} 0 1 0 {utmoffset[1]} 0 0 1 0 0 0 0 1"',
f'--writers.las.offset_x={reconstruction.georef.utm_east_offset}' ,
f'--writers.las.offset_y={reconstruction.georef.utm_north_offset}',
f'--writers.las.scale_x={las_scale}',
f'--writers.las.scale_y={las_scale}',
f'--writers.las.scale_z={las_scale}',
'--writers.las.offset_z=0',
f'--writers.las.a_srs="{reconstruction.georef.proj4()}"' # HOBU this should maybe be WKT
]
las_writer_def = {
"filename": tree.odm_georeferencing_model_laz,
"a_srs": reconstruction.georef.proj4(),
"offset_x": utmoffset[0],
"offset_y": utmoffset[1],
"offset_z": 0,
"scale_x": las_scale,
"scale_y": las_scale,
"scale_z": las_scale,
}
if reconstruction.has_gcp() and io.file_exists(gcp_geojson_zip_export_file):
if os.path.getsize(gcp_geojson_zip_export_file) <= 65535:
log.ODM_INFO("Embedding GCP info in point cloud")
params += [
'--writers.las.vlrs="{\\\"filename\\\": \\\"%s\\\", \\\"user_id\\\": \\\"ODM\\\", \\\"record_id\\\": 2, \\\"description\\\": \\\"Ground Control Points (zip)\\\"}"' % gcp_geojson_zip_export_file.replace(os.sep, "/")
]
las_writer_def["vlrs"] = json.dumps(
{
"filename": gcp_geojson_zip_export_file.replace(os.sep, "/"),
"user_id": "ODM",
"record_id": 2,
"description": "Ground Control Points (zip)"
}
)
else:
log.ODM_WARNING("Cannot embed GCP info in point cloud, %s is too large" % gcp_geojson_zip_export_file)
system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))
pipeline |= pdal.Writer.las(
**las_writer_def
)
pipeline.execute()
self.update_progress(50)
@ -202,7 +267,10 @@ class ODMGeoreferencingStage(types.ODM_Stage):
export_to_bounds_files(outputs['boundary'], reconstruction.get_proj_srs(), bounds_json, bounds_gpkg)
else:
log.ODM_INFO("Converting point cloud (non-georeferenced)")
system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))
pipeline |= pdal.Writer.las(
tree.odm_georeferencing_model_laz
)
pipeline.execute()
stats_dir = tree.path("opensfm", "stats", "codem")
@ -250,7 +318,7 @@ class ODMGeoreferencingStage(types.ODM_Stage):
except Exception as e:
log.ODM_WARNING("Cannot transform textured model: %s" % str(e))
os.rename(unaligned_obj, obj)
#TODO: seems gltf file is not converted in alignment?
for texturing in [tree.odm_texturing, tree.odm_25dtexturing]:
if reconstruction.multi_camera:
primary = get_primary_band_name(reconstruction.multi_camera, args.primary_band)

Wyświetl plik

@ -19,11 +19,11 @@ class ODMeshingStage(types.ODM_Stage):
# Create full 3D model unless --skip-3dmodel is set
if not args.skip_3dmodel:
if not io.file_exists(tree.odm_mesh) or self.rerun():
log.ODM_INFO('Writing ODM Mesh file in: %s' % tree.odm_mesh)
if not io.file_exists(tree.odm_mesh_topo) or self.rerun():
log.ODM_INFO('Writing ODM Mesh file in: %s' % tree.odm_mesh_topo)
mesh.screened_poisson_reconstruction(tree.filtered_point_cloud,
tree.odm_mesh,
mesh.screened_poisson_reconstruction(tree.filtered_point_cloud_topo,
tree.odm_mesh_topo,
depth=self.params.get('oct_tree'),
samples=self.params.get('samples'),
maxVertexCount=self.params.get('max_vertex'),
@ -31,16 +31,16 @@ class ODMeshingStage(types.ODM_Stage):
threads=max(1, self.params.get('max_concurrency') - 1)), # poissonrecon can get stuck on some machines if --threads == all cores
else:
log.ODM_WARNING('Found a valid ODM Mesh file in: %s' %
tree.odm_mesh)
tree.odm_mesh_topo)
self.update_progress(50)
# Always generate a 2.5D mesh
# unless --use-3dmesh is set.
if not args.use_3dmesh:
if not io.file_exists(tree.odm_25dmesh) or self.rerun():
if not io.file_exists(tree.odm_25dmesh_topo) or self.rerun():
log.ODM_INFO('Writing ODM 2.5D Mesh file in: %s' % tree.odm_25dmesh)
log.ODM_INFO('Writing ODM 2.5D Mesh file in: %s' % tree.odm_25dmesh_topo)
multiplier = math.pi / 2.0
radius_steps = commands.get_dem_radius_steps(tree.filtered_point_cloud_stats, 3, args.orthophoto_resolution, multiplier=multiplier)
@ -51,7 +51,7 @@ class ODMeshingStage(types.ODM_Stage):
if args.fast_orthophoto:
dsm_resolution *= 8.0
mesh.create_25dmesh(tree.filtered_point_cloud, tree.odm_25dmesh,
mesh.create_25dmesh(tree.filtered_point_cloud_topo, tree.odm_25dmesh_topo,
radius_steps,
dsm_resolution=dsm_resolution,
depth=self.params.get('oct_tree'),
@ -63,5 +63,5 @@ class ODMeshingStage(types.ODM_Stage):
max_tiles=None if reconstruction.has_geotagged_photos() else math.ceil(len(reconstruction.photos) / 2))
else:
log.ODM_WARNING('Found a valid ODM 2.5D Mesh file in: %s' %
tree.odm_25dmesh)
tree.odm_25dmesh_topo)

Wyświetl plik

@ -69,13 +69,13 @@ class ODMOpenSfMStage(types.ODM_Stage):
self.update_progress(75)
# We now switch to a geographic CRS
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))
shutil.move(tree.opensfm_reconstruction, tree.opensfm_topocentric_reconstruction)
shutil.move(tree.opensfm_geocoords_reconstruction, tree.opensfm_reconstruction)
else:
log.ODM_WARNING("Will skip exporting %s" % tree.opensfm_geocoords_reconstruction)
# 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))
# shutil.move(tree.opensfm_reconstruction, tree.opensfm_topocentric_reconstruction)
# shutil.move(tree.opensfm_geocoords_reconstruction, tree.opensfm_reconstruction)
# else:
# log.ODM_WARNING("Will skip exporting %s" % tree.opensfm_geocoords_reconstruction)
self.update_progress(80)