2018-06-12 21:36:24 +00:00
|
|
|
from __future__ import absolute_import
|
2018-07-07 04:15:15 +00:00
|
|
|
import os, shutil, sys, struct, random, math
|
2018-06-10 18:57:16 +00:00
|
|
|
from gippy import GeoImage
|
|
|
|
from opendm.dem import commands
|
|
|
|
from opendm import system
|
|
|
|
from opendm import log
|
2018-07-03 15:24:47 +00:00
|
|
|
from opendm import context
|
2018-06-18 13:54:09 +00:00
|
|
|
from scipy import signal, ndimage
|
|
|
|
import numpy as np
|
2018-06-10 18:57:16 +00:00
|
|
|
|
2018-07-03 17:01:18 +00:00
|
|
|
def create_25dmesh(inPointCloud, outMesh, dsm_resolution=0.05, depth=8, samples=1, maxVertexCount=100000, verbose=False, max_workers=None):
|
2018-06-10 18:57:16 +00:00
|
|
|
# Create DSM from point cloud
|
|
|
|
|
|
|
|
# Create temporary directory
|
|
|
|
mesh_directory = os.path.dirname(outMesh)
|
|
|
|
tmp_directory = os.path.join(mesh_directory, 'tmp')
|
|
|
|
if os.path.exists(tmp_directory):
|
|
|
|
shutil.rmtree(tmp_directory)
|
|
|
|
os.mkdir(tmp_directory)
|
|
|
|
log.ODM_INFO('Created temporary directory: %s' % tmp_directory)
|
|
|
|
|
2018-07-07 04:15:15 +00:00
|
|
|
radius_steps = [dsm_resolution * math.sqrt(2)]
|
2018-06-10 18:57:16 +00:00
|
|
|
|
|
|
|
log.ODM_INFO('Creating DSM for 2.5D mesh')
|
|
|
|
|
|
|
|
commands.create_dems(
|
2018-06-12 21:36:24 +00:00
|
|
|
[inPointCloud],
|
2018-06-10 18:57:16 +00:00
|
|
|
'mesh_dsm',
|
|
|
|
radius=map(str, radius_steps),
|
|
|
|
gapfill=True,
|
|
|
|
outdir=tmp_directory,
|
|
|
|
resolution=dsm_resolution,
|
2018-06-18 13:54:09 +00:00
|
|
|
products=['idw'],
|
2018-07-03 17:01:18 +00:00
|
|
|
verbose=verbose,
|
|
|
|
max_workers=max_workers
|
2018-06-10 18:57:16 +00:00
|
|
|
)
|
|
|
|
|
2018-07-06 20:41:56 +00:00
|
|
|
dsm_points = dem_to_points(os.path.join(tmp_directory, 'mesh_dsm.tif'), os.path.join(tmp_directory, 'dsm_points.ply'), verbose)
|
2018-07-06 14:29:55 +00:00
|
|
|
mesh = screened_poisson_reconstruction(dsm_points, outMesh, depth=depth,
|
|
|
|
samples=samples,
|
|
|
|
maxVertexCount=maxVertexCount,
|
|
|
|
threads=max_workers,
|
|
|
|
verbose=verbose)
|
2018-07-06 14:01:33 +00:00
|
|
|
|
2018-06-10 18:57:16 +00:00
|
|
|
# Cleanup tmp
|
2018-07-02 19:21:30 +00:00
|
|
|
if os.path.exists(tmp_directory):
|
|
|
|
shutil.rmtree(tmp_directory)
|
2018-06-10 18:57:16 +00:00
|
|
|
|
|
|
|
return mesh
|
|
|
|
|
2018-07-06 20:41:56 +00:00
|
|
|
def dem_to_points(inGeotiff, outPointCloud, verbose=False):
|
2018-06-10 18:57:16 +00:00
|
|
|
log.ODM_INFO('Sampling points from DSM: %s' % inGeotiff)
|
|
|
|
|
2018-07-06 20:41:56 +00:00
|
|
|
kwargs = {
|
|
|
|
'bin': context.odm_modules_path,
|
|
|
|
'outfile': outPointCloud,
|
|
|
|
'infile': inGeotiff,
|
|
|
|
'verbose': '-verbose' if verbose else ''
|
|
|
|
}
|
2018-06-10 18:57:16 +00:00
|
|
|
|
2018-07-06 20:41:56 +00:00
|
|
|
system.run('{bin}/odm_dem2points -inputFile {infile} '
|
|
|
|
'-outputFile {outfile} '
|
|
|
|
'-skirtHeightThreshold 1.5 '
|
|
|
|
'-skirtIncrements 0.2 '
|
|
|
|
'-skirtHeightCap 100 '
|
|
|
|
' {verbose} '.format(**kwargs))
|
2018-06-10 18:57:16 +00:00
|
|
|
|
|
|
|
return outPointCloud
|
|
|
|
|
2018-07-06 20:41:56 +00:00
|
|
|
# Old Python implementation of dem_to_points
|
|
|
|
# def dem_to_points(inGeotiff, outPointCloud):
|
|
|
|
# log.ODM_INFO('Sampling points from DSM: %s' % inGeotiff)
|
|
|
|
|
|
|
|
# image = GeoImage.open([inGeotiff], bandnames=['z'], nodata=-9999)
|
|
|
|
# arr = image['z'].read_raw()
|
|
|
|
|
|
|
|
# mem_file = BytesIO()
|
|
|
|
|
|
|
|
# xmin, xmax, ymin, ymax = image.extent().x0(), image.extent().x1(), image.extent().y0(), image.extent().y1()
|
|
|
|
# ext_width, ext_height = xmax - xmin, ymax - ymin
|
|
|
|
# arr_height, arr_width = arr.shape
|
|
|
|
# vertex_count = (arr_height - 4) * (arr_width - 4)
|
|
|
|
# skirt_points = 0
|
|
|
|
|
|
|
|
# skirt_height_threshold = 1 # meter
|
|
|
|
# skirt_increments = 0.1
|
|
|
|
|
|
|
|
# for x in range(2, arr_width - 2):
|
|
|
|
# for y in range(2, arr_height - 2):
|
|
|
|
# z = arr[y][x]
|
|
|
|
# tx = xmin + (float(x) / float(arr_width)) * ext_width
|
|
|
|
# ty = ymax - (float(y) / float(arr_height)) * ext_height
|
|
|
|
# mem_file.write(struct.pack('ffffff', tx, ty, z, 0, 0, 1))
|
|
|
|
|
|
|
|
# # Skirting
|
|
|
|
# for (nx, ny) in ((x, y + 1), (x, y - 1), (x + 1, y), (x - 1, y)):
|
|
|
|
# current_z = z
|
|
|
|
# neighbor_z = arr[ny][nx]
|
|
|
|
|
|
|
|
# if current_z - neighbor_z > skirt_height_threshold:
|
|
|
|
# while current_z > neighbor_z:
|
|
|
|
# current_z -= skirt_increments
|
|
|
|
# mem_file.write(struct.pack('ffffff', tx, ty, current_z, 0, 0, 1))
|
|
|
|
# skirt_points += 1
|
|
|
|
|
|
|
|
# mem_file.write(struct.pack('ffffff', tx, ty, neighbor_z, 0, 0, 1))
|
|
|
|
# skirt_points += 1
|
|
|
|
|
|
|
|
# log.ODM_INFO("Points count: %s (%s samples, %s skirts)", vertex_count + skirt_points, vertex_count, skirt_points)
|
|
|
|
# log.ODM_INFO('Writing points...')
|
|
|
|
|
|
|
|
# mem_file.seek(0)
|
|
|
|
# with open(outPointCloud, "wb") as f:
|
|
|
|
# f.write("ply\n")
|
|
|
|
# f.write("format binary_%s_endian 1.0\n" % sys.byteorder)
|
|
|
|
# f.write("element vertex %s\n" % (vertex_count + skirt_points))
|
|
|
|
# f.write("property float x\n")
|
|
|
|
# f.write("property float y\n")
|
|
|
|
# f.write("property float z\n")
|
|
|
|
# f.write("property float nx\n")
|
|
|
|
# f.write("property float ny\n")
|
|
|
|
# f.write("property float nz\n")
|
|
|
|
# f.write("end_header\n")
|
|
|
|
# shutil.copyfileobj(mem_file, f)
|
|
|
|
|
|
|
|
# mem_file.close()
|
|
|
|
|
|
|
|
# log.ODM_INFO('Wrote points to: %s' % outPointCloud)
|
|
|
|
|
|
|
|
# return outPointCloud
|
|
|
|
|
2018-06-10 18:57:16 +00:00
|
|
|
|
2018-07-06 14:01:33 +00:00
|
|
|
def screened_poisson_reconstruction(inPointCloud, outMesh, depth = 8, samples = 1, maxVertexCount=100000, pointWeight=4, threads=context.num_cores, verbose=False):
|
2018-07-03 15:24:47 +00:00
|
|
|
|
|
|
|
mesh_path, mesh_filename = os.path.split(outMesh)
|
|
|
|
# mesh_path = path/to
|
|
|
|
# mesh_filename = odm_mesh.ply
|
|
|
|
|
|
|
|
basename, ext = os.path.splitext(mesh_filename)
|
|
|
|
# basename = odm_mesh
|
|
|
|
# ext = .ply
|
|
|
|
|
|
|
|
outMeshDirty = os.path.join(mesh_path, "{}.dirty{}".format(basename, ext))
|
|
|
|
|
|
|
|
poissonReconArgs = {
|
2018-07-06 14:01:33 +00:00
|
|
|
'bin': context.poisson_recon_path,
|
2018-07-03 15:24:47 +00:00
|
|
|
'outfile': outMeshDirty,
|
2018-06-10 18:57:16 +00:00
|
|
|
'infile': inPointCloud,
|
|
|
|
'depth': depth,
|
|
|
|
'samples': samples,
|
2018-07-06 14:01:33 +00:00
|
|
|
'pointWeight': pointWeight,
|
2018-07-06 14:29:55 +00:00
|
|
|
'threads': threads,
|
2018-06-10 18:57:16 +00:00
|
|
|
'verbose': '--verbose' if verbose else ''
|
|
|
|
}
|
|
|
|
|
2018-07-06 14:29:55 +00:00
|
|
|
# Run PoissonRecon
|
2018-07-06 14:01:33 +00:00
|
|
|
system.run('{bin} --in {infile} '
|
|
|
|
'--out {outfile} '
|
|
|
|
'--depth {depth} '
|
2018-07-06 14:29:55 +00:00
|
|
|
'--pointWeight {pointWeight} '
|
|
|
|
'--samplesPerNode {samples} '
|
|
|
|
'--threads {threads} '
|
2018-07-06 14:01:33 +00:00
|
|
|
'--linearFit '
|
|
|
|
'{verbose}'.format(**poissonReconArgs))
|
2018-07-03 15:24:47 +00:00
|
|
|
|
|
|
|
# Cleanup and reduce vertex count if necessary
|
|
|
|
cleanupArgs = {
|
|
|
|
'bin': context.odm_modules_path,
|
|
|
|
'outfile': outMesh,
|
|
|
|
'infile': outMeshDirty,
|
|
|
|
'max_vertex': maxVertexCount,
|
|
|
|
'verbose': '-verbose' if verbose else ''
|
|
|
|
}
|
|
|
|
|
|
|
|
system.run('{bin}/odm_cleanmesh -inputFile {infile} '
|
|
|
|
'-outputFile {outfile} '
|
|
|
|
'-removeIslands '
|
|
|
|
'-decimateMesh {max_vertex} {verbose} '.format(**cleanupArgs))
|
|
|
|
|
|
|
|
# Delete intermediate results
|
|
|
|
os.remove(outMeshDirty)
|
2018-06-10 18:57:16 +00:00
|
|
|
|
|
|
|
return outMesh
|