diff --git a/README.md b/README.md index 1f2926d8..d40f7c8c 100644 --- a/README.md +++ b/README.md @@ -14,22 +14,22 @@ In a word, OpenDroneMap is a toolchain for processing raw civilian UAS imagery t 2. Digital Surface Models 3. Textured Digital Surface Models 4. Orthorectified Imagery -5. Classified Point Clouds +5. Classified Point Clouds (coming soon) 6. Digital Elevation Models 7. etc. -So far, it does Point Clouds, Digital Surface Models, Textured Digital Surface Models, and Orthorectified Imagery. Open Drone Map now includes state-of-the-art 3D reconstruction work by Michael Waechter, Nils Moehrle, and Michael Goesele. See their publication at http://www.gcc.tu-darmstadt.de/media/gcc/papers/Waechter-2014-LTB.pdf. +Open Drone Map now includes state-of-the-art 3D reconstruction work by Michael Waechter, Nils Moehrle, and Michael Goesele. See their publication at http://www.gcc.tu-darmstadt.de/media/gcc/papers/Waechter-2014-LTB.pdf. ## QUICKSTART -OpenDroneMap can run natively on Ubuntu 14.04 or later, see [Build and Run Using Docker](#build-and-run-using-docker) for running on Windows / MacOS. A Vagrant VM is also available: https://github.com/OpenDroneMap/odm_vagrant. +OpenDroneMap can run natively on Ubuntu 14.04 or later, see [Build and Run Using Docker](#build-and-run-using-docker) for running on Windows / MacOS. A [Vagrant VM](https://github.com/OpenDroneMap/odm_vagrant) is also available. *Support for Ubuntu 12.04 is currently BROKEN with the addition of OpenSfM and Ceres-Solver. It is likely to remain broken unless a champion is found to fix it.* **[Download the latest release here](https://github.com/OpenDroneMap/OpenDroneMap/releases)** -Current version: 0.2 (this software is in beta) +Current version: 0.3 (this software is in beta) 1. Extract and enter the OpenDroneMap directory 2. Run `bash configure.sh install` @@ -74,7 +74,7 @@ Note that using `run.sh` sets these temporarily in the shell. First you need a set of images, taken from a drone or otherwise. Example data can be obtained from https://github.com/OpenDroneMap/odm_data -Next, you need to copy over the settings file `default.settings.yaml` and edit it. The only setting you must edit is the `project-path` key. Set this to an empty directory within projects will be saved. There are many options for tuning your project. See the [wiki](https://github.com/OpenDroneMap/OpenDroneMap/wiki/Run-Time-Parameters) or run `python run.py -h` +Next, you need to edit the `settings.yaml` file. The only setting you must edit is the `project-path` key. Set this to an empty directory within projects will be saved. There are many options for tuning your project. See the [wiki](https://github.com/OpenDroneMap/OpenDroneMap/wiki/Run-Time-Parameters) or run `python run.py -h` Then run: @@ -146,7 +146,6 @@ instructions through "Create a Docker group". Once Docker is installed, the fast If you want to build your own Docker image from sources, type: - docker build -t packages -f packages.Dockerfile . docker build -t my_odm_image . docker run -it --rm -v $(pwd)/images:/code/images -v $(pwd)/odm_orthophoto:/code/odm_orthophoto -v $(pwd)/odm_texturing:/code/odm_texturing my_odm_image diff --git a/VERSION b/VERSION index f55f8fd9..1d71ef97 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.3-development \ No newline at end of file +0.3 \ No newline at end of file diff --git a/docker.settings.yaml b/docker.settings.yaml index dd9dfb4d..8109ea52 100644 --- a/docker.settings.yaml +++ b/docker.settings.yaml @@ -42,6 +42,10 @@ project_path: '/' #DO NOT CHANGE THIS OR DOCKER WILL NOT WORK. It should be '/' #texturing_keep_unseen_faces: False #texturing_tone_mapping: 'none' #gcp: !!null # YAML tag for None +#dem: False +#dem_sample_radius: 1.0 +#dem_resolution: 2 +#dem_radius: 0.5 #use_exif: False # Set to True if you have a GCP file (it auto-detects) and want to use EXIF #orthophoto_resolution: 20.0 # Pixels/meter #orthophoto_target_srs: !!null # Currently does nothing diff --git a/modules/odm_25dmeshing/src/CGAL.hpp b/modules/odm_25dmeshing/src/CGAL.hpp new file mode 100644 index 00000000..b203cae9 --- /dev/null +++ b/modules/odm_25dmeshing/src/CGAL.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::FT FT; +typedef Kernel::Point_3 Point3; +typedef Kernel::Vector_3 Vector3; + +//We define a vertex_base with info. The "info" (size_t) allow us to keep track of the original point index. +typedef CGAL::Triangulation_vertex_base_with_info_2 Vb; +typedef CGAL::Constrained_triangulation_face_base_2 Fb; +typedef CGAL::Triangulation_data_structure_2 Tds; +typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; + +typedef CGAL::Polyhedron_3 Polyhedron; +typedef Polyhedron::HalfedgeDS HalfedgeDS; + +namespace SMS = CGAL::Surface_mesh_simplification; + +// Concurrency +#ifdef CGAL_LINKED_WITH_TBB +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + +//typedef CGAL::First_of_pair_property_map Point_map; +//typedef CGAL::Second_of_pair_property_map Normal_map; + diff --git a/modules/odm_25dmeshing/src/Odm25dMeshing.cpp b/modules/odm_25dmeshing/src/Odm25dMeshing.cpp index f036a3a5..35307e26 100644 --- a/modules/odm_25dmeshing/src/Odm25dMeshing.cpp +++ b/modules/odm_25dmeshing/src/Odm25dMeshing.cpp @@ -204,7 +204,7 @@ void Odm25dMeshing::buildMesh(){ pointCount = gridPoints.size(); log << "sampled " << (pointCountBeforeGridSampling - pointCount) << " points\n"; - const double RETAIN_PERCENTAGE = std::min(100., 100. * static_cast(maxVertexCount) / static_cast(pointCount)); // percentage of points to retain. + const double RETAIN_PERCENTAGE = std::min(100., 100. * static_cast(maxVertexCount / 4) / static_cast(pointCount)); // percentage of points to retain. std::vector simplifiedPoints; log << "Performing weighted locally optimal projection simplification and regularization (retain: " << RETAIN_PERCENTAGE << "%, iterate: " << wlopIterations << ")" << "\n"; diff --git a/modules/odm_25dmeshing/src/PolyhedronBuilder.cpp b/modules/odm_25dmeshing/src/PolyhedronBuilder.cpp new file mode 100644 index 00000000..69b06d44 --- /dev/null +++ b/modules/odm_25dmeshing/src/PolyhedronBuilder.cpp @@ -0,0 +1,2 @@ +#include "PolyhedronBuilder.hpp" + diff --git a/modules/odm_25dmeshing/src/PolyhedronBuilder.hpp b/modules/odm_25dmeshing/src/PolyhedronBuilder.hpp new file mode 100644 index 00000000..f68dbe0c --- /dev/null +++ b/modules/odm_25dmeshing/src/PolyhedronBuilder.hpp @@ -0,0 +1,43 @@ +#include +#include + +#include +#include + +#include "CGAL.hpp" + +// A modifier creating a triangle with the incremental builder. +template +class PolyhedronBuilder : public CGAL::Modifier_base { +public: + std::vector &vertices; + std::vector &vertexIndices; + + PolyhedronBuilder( std::vector &vertices, std::vector &vertexIndices ) + : vertices(vertices), vertexIndices(vertexIndices) {} + + void operator()( HDS& hds) { + typedef typename HDS::Vertex Vertex; + typedef typename Vertex::Point Point; + + CGAL::Polyhedron_incremental_builder_3 builder( hds, true); + builder.begin_surface( vertices.size() / 3, vertexIndices.size() / 3 ); + + for(size_t i = 0; i < vertices.size(); i+=3 ){ + builder.add_vertex(Point(vertices[i+0], vertices[i+1], vertices[i+2])); + } + + for(size_t i = 0; i < vertexIndices.size(); i+=3){ + builder.begin_facet(); + builder.add_vertex_to_facet(vertexIndices[i+0]); + builder.add_vertex_to_facet(vertexIndices[i+1]); + builder.add_vertex_to_facet(vertexIndices[i+2]); + builder.end_facet(); + } + + // finish up the surface + builder.end_surface(); + } +}; + + diff --git a/opendm/config.py b/opendm/config.py index fe880d6a..70e19293 100644 --- a/opendm/config.py +++ b/opendm/config.py @@ -153,10 +153,10 @@ def config(): help=('The maximum number of processes to use in dense ' 'reconstruction. Default: %(default)s')) - parser.add_argument('--skip-25dmesh', + parser.add_argument('--use-25dmesh', action='store_true', default=False, - help='Do not build a 2.5D mesh and use the poisson mesh to compute the orthophoto') + help='Use a 2.5D mesh to compute the orthophoto. This option tends to provide better results for planar surfaces. Experimental.') parser.add_argument('--use-pmvs', action='store_true', @@ -336,6 +336,34 @@ def config(): help=('Use this tag if you have a gcp_list.txt but ' 'want to use the exif geotags instead')) + parser.add_argument('--dem', + action='store_true', + default=False, + help='Use this tag to build a DEM using a progressive ' + 'morphological filter in PDAL.') + + parser.add_argument('--dem-sample-radius', + metavar='', + default=1.0, + type=float, + help='Minimum distance between samples for DEM ' + 'generation.\nDefault=%(default)s') + + parser.add_argument('--dem-resolution', + metavar='', + type=float, + default=2, + help='Length of raster cell edges in X/Y units.' + '\nDefault: %(default)s') + + parser.add_argument('--dem-radius', + metavar='', + type=float, + default=0.5, + help='Radius about cell center bounding points to ' + 'use to calculate a cell value.\nDefault: ' + '%(default)s') + parser.add_argument('--orthophoto-resolution', metavar=' 0.0>', default=20.0, diff --git a/opendm/types.py b/opendm/types.py index 9e4799b8..b874bcec 100644 --- a/opendm/types.py +++ b/opendm/types.py @@ -143,13 +143,13 @@ class ODM_GeoRef(object): def coord_to_fractions(self, coord, refs): deg_dec = abs(float(coord)) deg = int(deg_dec) - minute_dec = (deg_dec-deg)*60 + minute_dec = (deg_dec - deg) * 60 minute = int(minute_dec) - sec_dec = (minute_dec-minute)*60 - sec_dec = round(sec_dec,3) + sec_dec = (minute_dec - minute) * 60 + sec_dec = round(sec_dec, 3) sec_denominator = 1000 - sec_numerator = int(sec_dec*sec_denominator) + sec_numerator = int(sec_dec * sec_denominator) if float(coord) >= 0: latRef = refs[0] else: @@ -158,7 +158,7 @@ class ODM_GeoRef(object): output = str(deg) + '/1 ' + str(minute) + '/1 ' + str(sec_numerator) + '/' + str(sec_denominator) return output, latRef - def convert_to_las(self, _file, pdalXML): + def convert_to_las(self, _file, _file_out, json_file): if not self.epsg: log.ODM_ERROR('Empty EPSG: Could not convert to LAS') @@ -166,50 +166,82 @@ class ODM_GeoRef(object): kwargs = {'bin': context.pdal_path, 'f_in': _file, - 'f_out': _file + '.las', + 'f_out': _file_out, 'east': self.utm_east_offset, 'north': self.utm_north_offset, 'epsg': self.epsg, - 'xml': pdalXML} + 'json': json_file} - # call txt2las - # system.run('{bin}/txt2las -i {f_in} -o {f_out} -skip 30 -parse xyzRGBssss ' \ - # '-set_scale 0.01 0.01 0.01 -set_offset {east} {north} 0 ' \ - # '-translate_xyz 0 -epsg {epsg}'.format(**kwargs)) - # # create pipeline file transform.xml to enable transformation - pipelineXml = '' - pipelineXml += '' - pipelineXml += ' ' - pipelineXml += ' ' - pipelineXml += ' ' - pipelineXml += ' ' - pipelineXml += ' ' - pipelineXml += ' ' - pipelineXml += ' ' - pipelineXml += ' ' - pipelineXml += ' ' - pipelineXml += ' ' - pipelineXml += '' + pipeline = '{{' \ + ' "pipeline":[' \ + ' "untransformed.ply",' \ + ' {{' \ + ' "type":"filters.transformation",' \ + ' "matrix":"1 0 0 {east} 0 1 0 {north} 0 0 1 0 0 0 0 1"' \ + ' }},' \ + ' {{' \ + ' "a_srs":"EPSG:{epsg}",' \ + ' "forward":"scale",' \ + ' "filename":"transformed.las"' \ + ' }}' \ + ' ]' \ + '}}'.format(**kwargs) - with open(pdalXML, 'w') as f: - f.write(pipelineXml) + with open(json_file, 'w') as f: + f.write(pipeline) # call pdal - system.run('{bin}/pdal pipeline -i {xml} --readers.ply.filename={f_in} ' + system.run('{bin}/pdal pipeline -i {json} --readers.ply.filename={f_in} ' '--writers.las.filename={f_out}'.format(**kwargs)) + def convert_to_dem(self, _file, _file_out, pdalJSON, sample_radius, gdal_res, gdal_radius): + # Check if exists f_in + if not io.file_exists(_file): + log.ODM_ERROR('LAS file does not exist') + return False + + kwargs = { + 'bin': context.pdal_path, + 'f_in': _file, + 'sample_radius': sample_radius, + 'gdal_res': gdal_res, + 'gdal_radius': gdal_radius, + 'f_out': _file_out, + 'json': pdalJSON + } + + pipelineJSON = '{{' \ + ' "pipeline":[' \ + ' "input.las",' \ + ' {{' \ + ' "type":"filters.sample",' \ + ' "radius":"{sample_radius}"' \ + ' }},' \ + ' {{' \ + ' "type":"filters.pmf",' \ + ' "extract":"true"' \ + ' }},' \ + ' {{' \ + ' "resolution": {gdal_res},' \ + ' "radius": {gdal_radius},' \ + ' "output_type":"idw",' \ + ' "filename":"outputfile.tif"' \ + ' }}' \ + ' ]' \ + '}}'.format(**kwargs) + + with open(pdalJSON, 'w') as f: + f.write(pipelineJSON) + + system.run('{bin}/pdal pipeline {json} --readers.las.filename={f_in} ' + '--writers.gdal.filename={f_out}'.format(**kwargs)) + + if io.file_exists(kwargs['f_out']): + return True + else: + return False + def utm_to_latlon(self, _file, _photo, idx): gcp = self.gcps[idx] @@ -271,7 +303,7 @@ class ODM_GeoRef(object): metadata[key] = pyexiv2.ExifTag(key, value) # GPS altitude - altitude = abs(int(float(latlon[2])*100)) + altitude = abs(int(float(latlon[2]) * 100)) key = 'Exif.GPSInfo.GPSAltitude' value = Fraction(altitude, 1) metadata[key] = pyexiv2.ExifTag(key, value) @@ -301,7 +333,7 @@ class ODM_GeoRef(object): log.ODM_DEBUG('Line: %s' % line) ref = line.split(' ') # match_wgs_utm = re.search('WGS84 UTM (\d{1,2})(N|S)', line, re.I) - if ref[0] == 'WGS84' and ref[1] == 'UTM': # match_wgs_utm: + if ref[0] == 'WGS84' and ref[1] == 'UTM': # match_wgs_utm: self.datum = ref[0] self.utm_pole = ref[2][len(ref) - 1] self.utm_zone = int(ref[2][:len(ref) - 1]) @@ -331,6 +363,7 @@ class ODM_GeoRef(object): x, y = xyz[:2] z = 0 self.gcps.append(ODM_GCPoint(float(x), float(y), float(z))) + # Write to json file class ODM_Tree(object): @@ -410,8 +443,14 @@ class ODM_Tree(object): self.odm_georeferencing_model_obj_geo = 'odm_textured_model_geo.obj' self.odm_georeferencing_xyz_file = io.join_paths( self.odm_georeferencing, 'odm_georeferenced_model.csv') - self.odm_georeferencing_pdal = io.join_paths( - self.odm_georeferencing, 'pipeline.xml') + self.odm_georeferencing_las_json = io.join_paths( + self.odm_georeferencing, 'las.json') + self.odm_georeferencing_model_las = io.join_paths( + self.odm_georeferencing, 'odm_georeferenced_model.las') + self.odm_georeferencing_dem = io.join_paths( + self.odm_georeferencing, 'odm_georeferencing_model_dem.tif') + self.odm_georeferencing_dem_json = io.join_paths( + self.odm_georeferencing, 'dem.json') # odm_orthophoto self.odm_orthophoto_file = io.join_paths(self.odm_orthophoto, 'odm_orthophoto.png') diff --git a/scripts/mvstex.py b/scripts/mvstex.py index 401b9ddd..c4e0c8ee 100644 --- a/scripts/mvstex.py +++ b/scripts/mvstex.py @@ -39,7 +39,7 @@ class ODMMvsTexCell(ecto.Cell): # define paths and create working directories system.mkdir_p(tree.odm_texturing) - if not args.skip_25dmesh: system.mkdir_p(tree.odm_25dtexturing) + if args.use_25dmesh: system.mkdir_p(tree.odm_25dtexturing) # check if we rerun cell or not rerun_cell = (args.rerun is not None and @@ -53,7 +53,7 @@ class ODMMvsTexCell(ecto.Cell): 'model': tree.odm_mesh }] - if not args.skip_25dmesh: + if args.use_25dmesh: runs += [{ 'out_dir': tree.odm_25dtexturing, 'model': tree.odm_25dmesh diff --git a/scripts/odm_app.py b/scripts/odm_app.py index 88632a73..2e9abdb6 100644 --- a/scripts/odm_app.py +++ b/scripts/odm_app.py @@ -72,6 +72,10 @@ class ODMApp(ecto.BlackBox): 'georeferencing': ODMGeoreferencingCell(img_size=p.args.resize_to, gcp_file=p.args.gcp, use_exif=p.args.use_exif, + dem=p.args.dem, + sample_radius=p.args.dem_sample_radius, + gdal_res=p.args.dem_resolution, + gdal_radius=p.args.dem_radius, verbose=p.args.verbose), 'orthophoto': ODMOrthoPhotoCell(resolution=p.args.orthophoto_resolution, t_srs=p.args.orthophoto_target_srs, diff --git a/scripts/odm_georeferencing.py b/scripts/odm_georeferencing.py index c735db4f..66395fab 100644 --- a/scripts/odm_georeferencing.py +++ b/scripts/odm_georeferencing.py @@ -17,6 +17,10 @@ class ODMGeoreferencingCell(ecto.Cell): 'northing height pixelrow pixelcol imagename', 'gcp_list.txt') params.declare("img_size", 'image size used in calibration', 2400) params.declare("use_exif", 'use exif', False) + params.declare("dem", 'Generate a dem', False) + params.declare("sample_radius", "Minimum distance between samples for DEM gen", 3) + params.declare("gdal_res", "Length of raster cell edges in X/Y units ", 2) + params.declare("gdal_radius", "Radius about cell center bounding points to use to calculate a cell value", 0.5) params.declare("verbose", 'print additional messages to console', False) def declare_io(self, params, inputs, outputs): @@ -48,7 +52,7 @@ class ODMGeoreferencingCell(ecto.Cell): # define paths and create working directories system.mkdir_p(tree.odm_georeferencing) - if not args.skip_25dmesh: system.mkdir_p(tree.odm_25dgeoreferencing) + if args.use_25dmesh: system.mkdir_p(tree.odm_25dgeoreferencing) # in case a gcp file it's not provided, let's try to generate it using # images metadata. Internally calls jhead. @@ -93,7 +97,7 @@ class ODMGeoreferencingCell(ecto.Cell): 'texturing_dir': tree.odm_texturing, 'model': os.path.join(tree.odm_texturing, tree.odm_textured_model_obj) }] - if not args.skip_25dmesh: + if args.use_25dmesh: runs += [{ 'georeferencing_dir': tree.odm_25dgeoreferencing, 'texturing_dir': tree.odm_25dtexturing, @@ -166,10 +170,24 @@ class ODMGeoreferencingCell(ecto.Cell): # convert ply model to LAS reference system geo_ref.convert_to_las(odm_georeferencing_model_ply_geo, - tree.odm_georeferencing_pdal) + tree.odm_georeferencing_model_las, + tree.odm_georeferencing_las_json) + + # If --dem, create a DEM + if args.dem: + demcreated = geo_ref.convert_to_dem(tree.odm_georeferencing_model_las, + tree.odm_georeferencing_dem, + tree.odm_georeferencing_dem_json, + self.params.sample_radius, + self.params.gdal_res, + self.params.gdal_radius) + if not demcreated: + log.ODM_WARNING('Something went wrong. Check the logs in odm_georeferencing.') + else: + log.ODM_INFO('DEM created at {0}'.format(tree.odm_georeferencing_dem)) # XYZ point cloud output - log.ODM_INFO("Creating geo-referenced CSV file (XYZ format, can be used with GRASS to create DEM)") + log.ODM_INFO("Creating geo-referenced CSV file (XYZ format)") with open(tree.odm_georeferencing_xyz_file, "wb") as csvfile: csvfile_writer = csv.writer(csvfile, delimiter=",") reachedpoints = False diff --git a/scripts/odm_meshing.py b/scripts/odm_meshing.py index 2b764c6e..8a05552c 100644 --- a/scripts/odm_meshing.py +++ b/scripts/odm_meshing.py @@ -85,7 +85,7 @@ class ODMeshingCell(ecto.Cell): tree.odm_mesh) # Do we need to generate a 2.5D mesh also? - if not args.skip_25dmesh: + if args.use_25dmesh: if not io.file_exists(tree.odm_25dmesh) or rerun_cell: log.ODM_DEBUG('Writing ODM 2.5D Mesh file in: %s' % tree.odm_25dmesh) diff --git a/scripts/odm_orthophoto.py b/scripts/odm_orthophoto.py index 7ecb204b..4174af42 100644 --- a/scripts/odm_orthophoto.py +++ b/scripts/odm_orthophoto.py @@ -58,12 +58,12 @@ class ODMOrthoPhotoCell(ecto.Cell): # Have geo coordinates? if io.file_exists(tree.odm_georeferencing_coords): - if not args.skip_25dmesh: + if args.use_25dmesh: kwargs['model_geo'] = os.path.join(tree.odm_25dtexturing, tree.odm_georeferencing_model_obj_geo) else: kwargs['model_geo'] = os.path.join(tree.odm_texturing, tree.odm_georeferencing_model_obj_geo) else: - if not args.skip_25dmesh: + if args.use_25dmesh: kwargs['model_geo'] = os.path.join(tree.odm_25dtexturing, tree.odm_textured_model_obj) else: kwargs['model_geo'] = os.path.join(tree.odm_texturing, tree.odm_textured_model_obj) diff --git a/settings.yaml b/settings.yaml index 7768477d..7eac080f 100644 --- a/settings.yaml +++ b/settings.yaml @@ -43,6 +43,10 @@ project_path: '' # Example: '/home/user/ODMProjects #texturing_tone_mapping: 'none' #gcp: !!null # YAML tag for None #use_exif: False # Set to True if you have a GCP file (it auto-detects) and want to use EXIF +#dem: False +#dem_sample_radius: 1.0 +#dem_resolution: 2 +#dem_radius: 0.5 #orthophoto_resolution: 20.0 # Pixels/meter #orthophoto_target_srs: !!null # Currently does nothing #orthophoto_no_tiled: False