diff --git a/Dockerfile b/Dockerfile index a480be8f..b19af620 100644 --- a/Dockerfile +++ b/Dockerfile @@ -14,9 +14,9 @@ RUN apt-get update -y # All packages (Will install much faster) RUN apt-get install --no-install-recommends -y git cmake python-pip build-essential software-properties-common python-software-properties libgdal-dev gdal-bin libgeotiff-dev \ libgtk2.0-dev libavcodec-dev libavformat-dev libswscale-dev python-dev python-numpy libtbb2 libtbb-dev libjpeg-dev libpng-dev libtiff-dev libjasper-dev libflann-dev \ -libproj-dev libxext-dev liblapack-dev libeigen3-dev libvtk5-dev python-networkx libgoogle-glog-dev libsuitesparse-dev libboost-filesystem-dev libboost-iostreams-dev \ +libproj-dev libxext-dev liblapack-dev libeigen3-dev libvtk6-dev python-networkx libgoogle-glog-dev libsuitesparse-dev libboost-filesystem-dev libboost-iostreams-dev \ libboost-regex-dev libboost-python-dev libboost-date-time-dev libboost-thread-dev python-pyproj python-empy python-nose python-pyside python-pyexiv2 python-scipy \ -libexiv2-dev liblas-bin python-matplotlib libatlas-base-dev libgmp-dev libmpfr-dev swig2.0 python-wheel libboost-log-dev libjsoncpp-dev +libexiv2-dev liblas-bin python-matplotlib libatlas-base-dev swig2.0 python-wheel libboost-log-dev libjsoncpp-dev python-gdal RUN apt-get remove libdc1394-22-dev RUN pip install --upgrade pip @@ -53,15 +53,14 @@ COPY VERSION /code/VERSION RUN cd SuperBuild && mkdir build && cd build && cmake .. && make -j$(nproc) && cd ../.. && mkdir build && cd build && cmake .. && make -j$(nproc) RUN apt-get -y remove libgl1-mesa-dri git cmake python-pip build-essential -RUN apt-get install -y libvtk5-dev +RUN apt-get install -y libvtk6-dev # Cleanup APT RUN apt-get clean && rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* # Clean Superbuild -RUN rm -rf /code/SuperBuild/download -RUN rm -rf /code/SuperBuild/src/opencv/samples /code/SuperBuild/src/pcl/test /code/SuperBuild/src/pcl/doc /code/SuperBuild/src/pdal/test /code/SuperBuild/src/pdal/doc +RUN rm -rf /code/SuperBuild/download /code/SuperBuild/src/vtk7 /code/SuperBuild/src/opencv /code/SuperBuild/src/pcl /code/SuperBuild/src/pdal /code/SuperBuild/src/opengv /code/SuperBuild/src/mvstexturing /code/SuperBuild/src/ceres /code/SuperBuild/build/vtk7 /code/SuperBuild/build/opencv # Entry point ENTRYPOINT ["python", "/code/run.py", "code"] diff --git a/README.md b/README.md index 7dbabc6e..caf3472c 100644 --- a/README.md +++ b/README.md @@ -147,8 +147,8 @@ You can also view the orthophoto GeoTIFF in [QGIS](http://www.qgis.org/) or othe ## Build and Run Using Docker -(Instructions below apply to Ubuntu 14.04, but the Docker image workflow -has equivalent procedures for Mac OS X and Windows. See [docs.docker.com](docs.docker.com)) +(Instructions below apply to Ubuntu 14.04, but the Docker image workflow +has equivalent procedures for Mac OS X and Windows. See [docs.docker.com](https://docs.docker.com/)) OpenDroneMap is Dockerized, meaning you can use containerization to build and run it without tampering with the configuration of libraries and packages already installed on your machine. Docker software is free to install and use in this context. If you don't have it installed, diff --git a/SuperBuild/CMakeLists.txt b/SuperBuild/CMakeLists.txt index 453c3fb4..ec40a608 100644 --- a/SuperBuild/CMakeLists.txt +++ b/SuperBuild/CMakeLists.txt @@ -91,12 +91,15 @@ SETUP_EXTERNAL_PROJECT(Ceres ${ODM_Ceres_Version} ${ODM_BUILD_Ceres}) # --------------------------------------------------------------------------------------------- -# CGAL -# -set(ODM_CGAL_Version 4.9) -option(ODM_BUILD_CGAL "Force to build CGAL library" OFF) +# VTK7 +# We need to build VTK from sources because Debian packages +# are built with DVTK_SMP_IMPLEMENTATION_TYPE set to +# "Sequential" which means no multithread support. -SETUP_EXTERNAL_PROJECT(CGAL ${ODM_CGAL_Version} ${ODM_BUILD_CGAL}) +set(ODM_VTK7_Version 7.1.1) +option(ODM_BUILD_VTK7 "Force to build VTK7 library" OFF) + +SETUP_EXTERNAL_PROJECT(VTK7 ${ODM_VTK7_Version} ${ODM_BUILD_VTK7}) # --------------------------------------------------------------------------------------------- # Hexer diff --git a/SuperBuild/cmake/External-CGAL.cmake b/SuperBuild/cmake/External-VTK7.cmake similarity index 77% rename from SuperBuild/cmake/External-CGAL.cmake rename to SuperBuild/cmake/External-VTK7.cmake index 51e09c4b..02a00efd 100644 --- a/SuperBuild/cmake/External-CGAL.cmake +++ b/SuperBuild/cmake/External-VTK7.cmake @@ -1,4 +1,4 @@ -set(_proj_name cgal) +set(_proj_name vtk7) set(_SB_BINARY_DIR "${SB_BINARY_DIR}/${_proj_name}") ExternalProject_Add(${_proj_name} @@ -7,14 +7,17 @@ ExternalProject_Add(${_proj_name} STAMP_DIR ${_SB_BINARY_DIR}/stamp #--Download step-------------- DOWNLOAD_DIR ${SB_DOWNLOAD_DIR}/${_proj_name} - URL https://github.com/CGAL/cgal/releases/download/releases%2FCGAL-4.9/CGAL-4.9.zip - URL_MD5 31c08d762a72fda785df194c89b833df + URL https://github.com/Kitware/VTK/archive/v7.1.1.zip #--Update/Patch step---------- UPDATE_COMMAND "" #--Configure step------------- SOURCE_DIR ${SB_SOURCE_DIR}/${_proj_name} CMAKE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${SB_INSTALL_DIR} + -DVTK_SMP_IMPLEMENTATION_TYPE=TBB + -DCMAKE_BUILD_TYPE=Release + -DVTK_Group_Rendering=OFF + -DBUILD_TESTING=OFF #--Build step----------------- BINARY_DIR ${_SB_BINARY_DIR} #--Install step--------------- diff --git a/configure.sh b/configure.sh index b06a1ec8..101e6f0e 100755 --- a/configure.sh +++ b/configure.sh @@ -28,7 +28,8 @@ install() { gdal-bin \ libgeotiff-dev \ pkg-config \ - libjsoncpp-dev + libjsoncpp-dev \ + python-gdal echo "Getting CMake 3.1 for MVS-Texturing" apt-get install -y software-properties-common python-software-properties @@ -79,9 +80,6 @@ install() { appsettings \ loky - echo "Installing CGAL dependencies" - apt-get install -y -qq libgmp-dev libmpfr-dev - echo "Installing Ecto Dependencies" pip install -U catkin-pkg apt-get install -y -qq python-empy \ diff --git a/core2.Dockerfile b/core2.Dockerfile index e300fbce..f129ca44 100644 --- a/core2.Dockerfile +++ b/core2.Dockerfile @@ -12,9 +12,9 @@ RUN apt-get update -y # All packages (Will install much faster) RUN apt-get install --no-install-recommends -y git cmake python-pip build-essential software-properties-common python-software-properties libgdal-dev gdal-bin libgeotiff-dev \ libgtk2.0-dev libavcodec-dev libavformat-dev libswscale-dev python-dev python-numpy libtbb2 libtbb-dev libjpeg-dev libpng-dev libtiff-dev libjasper-dev libflann-dev \ -libproj-dev libxext-dev liblapack-dev libeigen3-dev libvtk5-dev python-networkx libgoogle-glog-dev libsuitesparse-dev libboost-filesystem-dev libboost-iostreams-dev \ +libproj-dev libxext-dev liblapack-dev libeigen3-dev libvtk6-dev python-networkx libgoogle-glog-dev libsuitesparse-dev libboost-filesystem-dev libboost-iostreams-dev \ libboost-regex-dev libboost-python-dev libboost-date-time-dev libboost-thread-dev python-pyproj python-empy python-nose python-pyside python-pyexiv2 python-scipy \ -libexiv2-dev liblas-bin python-matplotlib libatlas-base-dev libgmp-dev libmpfr-dev swig2.0 python-wheel libboost-log-dev libjsoncpp-dev +libexiv2-dev liblas-bin python-matplotlib libatlas-base-dev swig2.0 python-wheel libboost-log-dev libjsoncpp-dev python-gdal RUN apt-get remove libdc1394-22-dev RUN pip install --upgrade pip @@ -54,15 +54,14 @@ RUN mv -v /usr/bin/gcc /usr/bin/gcc_real && mv -v /usr/bin/g++ /usr/bin/g++_real RUN cd SuperBuild && mkdir build && cd build && cmake .. && make -j$(nproc) && cd ../.. && mkdir build && cd build && cmake .. && make -j$(nproc) RUN apt-get -y remove libgl1-mesa-dri git cmake python-pip build-essential -RUN apt-get install -y libvtk5-dev +RUN apt-get install -y libvtk6-dev # Cleanup APT RUN apt-get clean && rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* # Clean Superbuild -RUN rm -rf /code/SuperBuild/download -RUN rm -rf /code/SuperBuild/src/opencv/samples /code/SuperBuild/src/pcl/test /code/SuperBuild/src/pcl/doc /code/SuperBuild/src/pdal/test /code/SuperBuild/src/pdal/doc +RUN rm -rf /code/SuperBuild/download /code/SuperBuild/src/vtk7 /code/SuperBuild/src/opencv /code/SuperBuild/src/pcl /code/SuperBuild/src/pdal /code/SuperBuild/src/opengv /code/SuperBuild/src/mvstexturing /code/SuperBuild/src/ceres /code/SuperBuild/build/vtk7 /code/SuperBuild/build/opencv # Entry point ENTRYPOINT ["python", "/code/run.py", "code"] diff --git a/modules/odm_25dmeshing/CMakeLists.txt b/modules/odm_25dmeshing/CMakeLists.txt index 4c19bbbe..708555e9 100644 --- a/modules/odm_25dmeshing/CMakeLists.txt +++ b/modules/odm_25dmeshing/CMakeLists.txt @@ -3,45 +3,21 @@ cmake_minimum_required(VERSION 2.8) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}) -# Set pcl dir to the input spedified with option -DCGAL_DIR="path" -set(CGAL_DIR "CGAL_DIR-NOTFOUND" CACHE "CGAL_DIR" "Path to the CGAL installation directory") +set (CMAKE_CXX_STANDARD 11) +#set(VTK_DIR "/data/packages/VTK-7.1.1-build") +set(VTK_SMP_IMPLEMENTATION_TYPE TBB) +find_package(VTK 7.1.1 REQUIRED) +include(${VTK_USE_FILE}) -# Link -find_package(CGAL COMPONENTS Core HINTS "${CGAL_DIR}") +# Add compiler options. +#add_definitions(-Wall -Wextra -O0 -g3) +add_definitions(-Wall -Wextra) -if ( CGAL_FOUND ) - - find_package( Eigen3 ) - - if( EIGEN3_FOUND ) - include_directories(${EIGEN3_INCLUDE_DIR}) - add_definitions(-DCGAL_EIGEN3_ENABLED) - else() - message(FATAL_ERROR "This program requires the Eigen3 library, and will not be compiled.") - endif() - - include( ${CGAL_USE_FILE} ) - - find_package( TBB ) - if( TBB_FOUND ) - include(${TBB_USE_FILE}) - else() - message(WARNING "TBB not found, parallel processing will be disabled.") - endif() - - # Add compiler options. - add_definitions(-Wall -Wextra -std=c++11) +# Add source directory +aux_source_directory("./src" SRC_LIST) - # Add source directory - aux_source_directory("./src" SRC_LIST) +# Add exectuteable +add_executable(${PROJECT_NAME} ${SRC_LIST}) - # Add exectuteable - add_executable(${PROJECT_NAME} ${SRC_LIST}) - - target_link_libraries(odm_25dmeshing ${CGAL_LIBRARIES} ${TBB_LIBRARIES}) -else() - - message(FATAL_ERROR "This program requires the CGAL library, and will not be compiled.") - -endif() +target_link_libraries(odm_25dmeshing ${VTK_LIBRARIES}) diff --git a/modules/odm_25dmeshing/FindEigen3.cmake b/modules/odm_25dmeshing/FindEigen3.cmake deleted file mode 100644 index e56005eb..00000000 --- a/modules/odm_25dmeshing/FindEigen3.cmake +++ /dev/null @@ -1,84 +0,0 @@ -# - Try to find Eigen3 lib -# -# This module supports requiring a minimum version, e.g. you can do -# find_package(Eigen3 3.1.2) -# to require version 3.1.2 or newer of Eigen3. -# -# Once done this will define -# -# EIGEN3_FOUND - system has eigen lib with correct version -# EIGEN3_INCLUDE_DIR - the eigen include directory -# EIGEN3_VERSION - eigen version - -# Copyright (c) 2006, 2007 Montel Laurent, -# Copyright (c) 2008, 2009 Gael Guennebaud, -# Copyright (c) 2009 Benoit Jacob -# Redistribution and use is allowed according to the terms of the 2-clause BSD license. - -include(FindPackageHandleStandardArgs) - -if(NOT Eigen3_FIND_VERSION) - if(NOT Eigen3_FIND_VERSION_MAJOR) - set(Eigen3_FIND_VERSION_MAJOR 2) - endif(NOT Eigen3_FIND_VERSION_MAJOR) - if(NOT Eigen3_FIND_VERSION_MINOR) - set(Eigen3_FIND_VERSION_MINOR 91) - endif(NOT Eigen3_FIND_VERSION_MINOR) - if(NOT Eigen3_FIND_VERSION_PATCH) - set(Eigen3_FIND_VERSION_PATCH 0) - endif(NOT Eigen3_FIND_VERSION_PATCH) - - set(Eigen3_FIND_VERSION "${Eigen3_FIND_VERSION_MAJOR}.${Eigen3_FIND_VERSION_MINOR}.${Eigen3_FIND_VERSION_PATCH}") -endif(NOT Eigen3_FIND_VERSION) - -macro(_eigen3_get_version) - file(READ "${EIGEN3_INCLUDE_DIR}/Eigen/src/Core/util/Macros.h" _eigen3_version_header) - - string(REGEX MATCH "define[ \t]+EIGEN_WORLD_VERSION[ \t]+([0-9]+)" _eigen3_world_version_match "${_eigen3_version_header}") - set(EIGEN3_WORLD_VERSION "${CMAKE_MATCH_1}") - string(REGEX MATCH "define[ \t]+EIGEN_MAJOR_VERSION[ \t]+([0-9]+)" _eigen3_major_version_match "${_eigen3_version_header}") - set(EIGEN3_MAJOR_VERSION "${CMAKE_MATCH_1}") - string(REGEX MATCH "define[ \t]+EIGEN_MINOR_VERSION[ \t]+([0-9]+)" _eigen3_minor_version_match "${_eigen3_version_header}") - set(EIGEN3_MINOR_VERSION "${CMAKE_MATCH_1}") - - set(EIGEN3_VERSION ${EIGEN3_WORLD_VERSION}.${EIGEN3_MAJOR_VERSION}.${EIGEN3_MINOR_VERSION}) -endmacro(_eigen3_get_version) - -if (EIGEN3_INCLUDE_DIR) - - if (EXISTS ${EIGEN3_INCLUDE_DIR}/signature_of_eigen3_matrix_library) - # in cache already and valid - _eigen3_get_version() - set(EIGEN3_FOUND ${EIGEN3_VERSION_OK}) - - find_package_handle_standard_args(Eigen3 - REQUIRED_VARS EIGEN3_INCLUDE_DIR - VERSION_VAR EIGEN3_VERSION) - - else() - message(STATUS "Eigen3 path specified in cmake variable EIGEN3_INCLUDE_DIR is " - "set to ${EIGEN3_INCLUDE_DIR}, but that path does not contains the file " - "signature_of_eigen3_matrix_library and is considered as invalid.") - endif() - - - -else (EIGEN3_INCLUDE_DIR) - - find_path(EIGEN3_INCLUDE_DIR NAMES signature_of_eigen3_matrix_library - HINTS ENV EIGEN3_INC_DIR - ENV EIGEN3_DIR - PATHS ${KDE4_INCLUDE_DIR} - PATH_SUFFIXES include eigen3 eigen - DOC "Directory containing the Eigen3 header files" - ) - - if(EIGEN3_INCLUDE_DIR) - _eigen3_get_version() - endif(EIGEN3_INCLUDE_DIR) - - find_package_handle_standard_args(Eigen3 - REQUIRED_VARS EIGEN3_INCLUDE_DIR - VERSION_VAR EIGEN3_VERSION) - -endif(EIGEN3_INCLUDE_DIR) diff --git a/modules/odm_25dmeshing/src/CGAL.hpp b/modules/odm_25dmeshing/src/CGAL.hpp deleted file mode 100644 index b203cae9..00000000 --- a/modules/odm_25dmeshing/src/CGAL.hpp +++ /dev/null @@ -1,38 +0,0 @@ -#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 9d1c1788..26816336 100644 --- a/modules/odm_25dmeshing/src/Odm25dMeshing.cpp +++ b/modules/odm_25dmeshing/src/Odm25dMeshing.cpp @@ -37,6 +37,290 @@ int Odm25dMeshing::run(int argc, char **argv) { return EXIT_SUCCESS; } +void Odm25dMeshing::loadPointCloud() { + log << "Loading point cloud... "; + + try{ + std::ifstream ss(inputFile, std::ios::binary); + if (ss.fail()) throw Odm25dMeshingException("Failed to open " + inputFile); + PlyFile file; + + file.parse_header(ss); + + std::shared_ptr vertices = file.request_properties_from_element("vertex", { "x", "y", "z" }); + file.read(ss); + + const size_t numVerticesBytes = vertices->buffer.size_bytes(); + struct float3 { float x, y, z; }; + struct double3 { double x, y, z; }; + + if (vertices->t == tinyply::Type::FLOAT32) { + std::vector verts(vertices->count); + std::memcpy(verts.data(), vertices->buffer.get(), numVerticesBytes); + for (float3 &v : verts){ + points->InsertNextPoint(v.x, v.y, v.z); + } + }else if (vertices->t == tinyply::Type::FLOAT64) { + std::vector verts(vertices->count); + std::memcpy(verts.data(), vertices->buffer.get(), numVerticesBytes); + for (double3 &v : verts){ + points->InsertNextPoint(v.x, v.y, v.z); + } + }else{ + throw Odm25dMeshingException("Invalid data type (only float32 and float64 are supported): " + std::to_string((int)vertices->t)); + } + } + catch (const std::exception & e) + { + throw Odm25dMeshingException("Error while loading point cloud: " + std::string(e.what())); + } + + log << "loaded " << points->GetNumberOfPoints() << " points\n"; +} + +void Odm25dMeshing::buildMesh(){ + vtkThreadedImageAlgorithm::SetGlobalDefaultEnableSMP(true); + + log << "Remove outliers... "; + + vtkSmartPointer polyPoints = + vtkSmartPointer::New(); + polyPoints->SetPoints(points); + + vtkSmartPointer locator = vtkSmartPointer::New(); + + vtkSmartPointer radiusRemoval = + vtkSmartPointer::New(); + radiusRemoval->SetInputData(polyPoints); + radiusRemoval->SetLocator(locator); + radiusRemoval->SetRadius(20); // 20 meters + radiusRemoval->SetNumberOfNeighbors(2); + + vtkSmartPointer statsRemoval = + vtkSmartPointer::New(); + statsRemoval->SetInputConnection(radiusRemoval->GetOutputPort()); + statsRemoval->SetLocator(locator); + statsRemoval->SetSampleSize(neighbors); + statsRemoval->SetStandardDeviationFactor(1.5); + statsRemoval->GenerateOutliersOff(); + statsRemoval->Update(); + + log << (radiusRemoval->GetNumberOfPointsRemoved() + statsRemoval->GetNumberOfPointsRemoved()) << " points removed\n"; + vtkSmartPointer cleanedPoints = statsRemoval->GetOutput()->GetPoints(); + + statsRemoval = nullptr; + radiusRemoval = nullptr; + polyPoints = nullptr; + + log << "Squash point cloud to plane... "; + + vtkSmartPointer elevation = vtkSmartPointer::New(); + elevation->SetName("elevation"); + elevation->SetNumberOfComponents(1); + double p[2]; + + for (vtkIdType i = 0; i < cleanedPoints->GetNumberOfPoints(); i++){ + cleanedPoints->GetPoint(i, p); + elevation->InsertNextValue(p[2]); + p[2] = 0.0; + cleanedPoints->SetPoint(i, p); + } + + log << "OK\n"; + + vtkSmartPointer polydataToProcess = + vtkSmartPointer::New(); + polydataToProcess->SetPoints(cleanedPoints); + polydataToProcess->GetPointData()->SetScalars(elevation); + + const float NODATA = -9999; + + double *bounds = polydataToProcess->GetBounds(); + + double centerX = polydataToProcess->GetCenter()[0]; + double centerY = polydataToProcess->GetCenter()[1]; + double centerZ = polydataToProcess->GetCenter()[2]; + + double extentX = bounds[1] - bounds[0]; + double extentY = bounds[3] - bounds[2]; + + if (resolution == 0.0){ + resolution = (double)maxVertexCount / (sqrt(extentX * extentY) * 75.0); + log << "Automatically set resolution to " << std::fixed << resolution << "\n"; + } + + int width = ceil(extentX * resolution); + int height = ceil(extentY * resolution); + + log << "Plane extentX: " << extentX << + ", extentY: " << extentY << "\n"; + + double planeCenter[3]; + planeCenter[0] = centerX; + planeCenter[1] = centerY; + planeCenter[2] = centerZ; + + vtkSmartPointer plane = + vtkSmartPointer::New(); + plane->SetResolution(width, height); + plane->SetOrigin(0.0, 0.0, 0.0); + plane->SetPoint1(extentX, 0.0, 0.0); + plane->SetPoint2(0.0, extentY, 0); + plane->SetCenter(planeCenter); + plane->SetNormal(0.0, 0.0, 1.0); + + vtkSmartPointer shepardKernel = + vtkSmartPointer::New(); + shepardKernel->SetPowerParameter(2.0); + shepardKernel->SetKernelFootprintToNClosest(); + shepardKernel->SetNumberOfPoints(neighbors); + + vtkSmartPointer image = + vtkSmartPointer::New(); + image->SetDimensions(width, height, 1); + log << "DSM size is " << width << "x" << height << " (" << ceil(width * height * sizeof(float) * 1e-6) << " MB) \n"; + image->AllocateScalars(VTK_FLOAT, 1); + + log << "Point interpolation using shepard's kernel... "; + + vtkSmartPointer interpolator = + vtkSmartPointer::New(); + interpolator->SetInputConnection(plane->GetOutputPort()); + interpolator->SetSourceData(polydataToProcess); + interpolator->SetKernel(shepardKernel); + interpolator->SetLocator(locator); + interpolator->SetNullValue(NODATA); + interpolator->Update(); + + vtkSmartPointer interpolatedPoly = + interpolator->GetPolyDataOutput(); + + log << "OK\nTransfering interpolation results to DSM... "; + + interpolator = nullptr; + polydataToProcess = nullptr; + elevation = nullptr; + cleanedPoints = nullptr; + plane = nullptr; + shepardKernel = nullptr; + locator = nullptr; + + vtkSmartPointer interpolatedElevation = + vtkFloatArray::SafeDownCast(interpolatedPoly->GetPointData()->GetArray("elevation")); + + for (int i = 0; i < width; i++){ + for (int j = 0; j < height; j++){ + float* pixel = static_cast(image->GetScalarPointer(i,j,0)); + vtkIdType cellId = interpolatedPoly->GetCell(j * width + i)->GetPointId(0); + pixel[0] = interpolatedElevation->GetValue(cellId); + } + } + + log << "OK\nMedian filter..."; + + vtkSmartPointer medianFilter = + vtkSmartPointer::New(); + medianFilter->SetInputData(image); + medianFilter->SetKernelSize( + std::max(1.0, resolution), + std::max(1.0, resolution), + 1); + medianFilter->Update(); + + log << "OK\n"; + +// double diffuseIterations = std::max(1.0, resolution / 2.0); +// vtkSmartPointer diffuse1 = +// vtkSmartPointer::New(); +// diffuse1->SetInputConnection(medianFilter->GetOutputPort()); +// diffuse1->FacesOn(); +// diffuse1->EdgesOn(); +// diffuse1->CornersOn(); +// diffuse1->SetDiffusionFactor(1); // Full strength +// diffuse1->GradientMagnitudeThresholdOn(); +// diffuse1->SetDiffusionThreshold(0.2); // Don't smooth jumps in elevation > than 0.20m +// diffuse1->SetNumberOfIterations(diffuseIterations); +// diffuse1->Update(); + + if (outputDsmFile != ""){ + log << "Saving DSM to file... "; + vtkSmartPointer tiffWriter = + vtkSmartPointer::New(); + tiffWriter->SetFileName(outputDsmFile.c_str()); + tiffWriter->SetInputData(medianFilter->GetOutput()); + tiffWriter->Write(); + log << "OK\n"; + } + + log << "Triangulate... "; + + vtkSmartPointer terrain = + vtkSmartPointer::New(); + terrain->SetErrorMeasureToNumberOfTriangles(); + terrain->SetNumberOfTriangles(maxVertexCount * 2); // Approximate + terrain->SetInputData(medianFilter->GetOutput()); + terrain->BoundaryVertexDeletionOn(); + + log << "OK\nTransform... "; + vtkSmartPointer transform = + vtkSmartPointer::New(); + transform->Translate(-extentX / 2.0 + centerX, + -extentY / 2.0 + centerY, 0); + transform->Scale(extentX / width, extentY / height, 1); + + vtkSmartPointer transformFilter = + vtkSmartPointer::New(); + transformFilter->SetInputConnection(terrain->GetOutputPort()); + transformFilter->SetTransform(transform); + + log << "OK\n"; + + log << "Saving mesh to file... "; + + vtkSmartPointer plyWriter = + vtkSmartPointer::New(); + plyWriter->SetFileName(outputFile.c_str()); + plyWriter->SetInputConnection(transformFilter->GetOutputPort()); + plyWriter->SetFileTypeToASCII(); + plyWriter->Write(); + + log << "OK\n"; + +#ifdef SUPPORTDEBUGWINDOW + if (showDebugWindow){ + vtkSmartPointer mapper = + vtkSmartPointer::New(); + mapper->SetInputConnection(transformFilter->GetOutputPort()); + mapper->SetScalarRange(150, 170); + + // vtkSmartPointer mapper = + // vtkSmartPointer::New(); + // mapper->SetInputData(image); + // mapper->SetScalarRange(150, 170); + + vtkSmartPointer actor = + vtkSmartPointer::New(); + actor->SetMapper(mapper); + actor->GetProperty()->SetPointSize(5); + + vtkSmartPointer renderer = + vtkSmartPointer::New(); + vtkSmartPointer renderWindow = + vtkSmartPointer::New(); + renderWindow->AddRenderer(renderer); + vtkSmartPointer renderWindowInteractor = + vtkSmartPointer::New(); + renderWindowInteractor->SetRenderWindow(renderWindow); + + renderer->AddActor(actor); + renderer->SetBackground(0.1804,0.5451,0.3412); // Sea green + + renderWindow->Render(); + renderWindowInteractor->Start(); + } +#endif +} + void Odm25dMeshing::parseArguments(int argc, char **argv) { for (int argIndex = 1; argIndex < argc; ++argIndex) { // The argument to be parsed. @@ -55,24 +339,23 @@ void Odm25dMeshing::parseArguments(int argc, char **argv) { if (ss.bad()) throw Odm25dMeshingException("Argument '" + argument + "' has a bad value (wrong type)."); maxVertexCount = std::max(maxVertexCount, 0); log << "Vertex count was manually set to: " << maxVertexCount << "\n"; - } else if (argument == "-outliersRemovalPercentage" && argIndex < argc) { + } else if (argument == "-resolution" && argIndex < argc) { ++argIndex; if (argIndex >= argc) throw Odm25dMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided."); std::stringstream ss(argv[argIndex]); - ss >> outliersRemovalPercentage; + ss >> resolution; if (ss.bad()) throw Odm25dMeshingException("Argument '" + argument + "' has a bad value (wrong type)."); - outliersRemovalPercentage = std::min(99.99, std::max(outliersRemovalPercentage, 0)); - log << "Outliers removal was manually set to: " << outliersRemovalPercentage << "\n"; - } else if (argument == "-wlopIterations" && argIndex < argc) { + resolution = std::min(100000, std::max(resolution, 0)); + log << "Resolution was manually set to: " << resolution << "\n"; + } else if (argument == "-neighbors" && argIndex < argc) { ++argIndex; if (argIndex >= argc) throw Odm25dMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided."); std::stringstream ss(argv[argIndex]); - ss >> wlopIterations; + ss >> neighbors; if (ss.bad()) throw Odm25dMeshingException("Argument '" + argument + "' has a bad value (wrong type)."); - - wlopIterations = std::min(1000, std::max(wlopIterations, 1)); - log << "WLOP iterations was manually set to: " << wlopIterations << "\n"; + neighbors = std::min(1000, std::max(neighbors, 1)); + log << "Neighbors was manually set to: " << neighbors << "\n"; } else if (argument == "-inputFile" && argIndex < argc) { ++argIndex; if (argIndex >= argc) { @@ -102,6 +385,23 @@ void Odm25dMeshing::parseArguments(int argc, char **argv) { } testFile.close(); log << "Writing output to: " << outputFile << "\n"; + }else if (argument == "-outputDsmFile" && argIndex < argc) { + ++argIndex; + if (argIndex >= argc) { + throw Odm25dMeshingException( + "Argument '" + argument + + "' expects 1 more input following it, but no more inputs were provided."); + } + outputDsmFile = std::string(argv[argIndex]); + std::ofstream testFile(outputDsmFile.c_str()); + if (!testFile.is_open()) { + throw Odm25dMeshingException( + "Argument '" + argument + "' has a bad value. (file not accessible)"); + } + testFile.close(); + log << "Saving DSM output to: " << outputDsmFile << "\n"; + } else if (argument == "-showDebugWindow") { + showDebugWindow = true; } else if (argument == "-logFile" && argIndex < argc) { ++argIndex; if (argIndex >= argc) { @@ -124,302 +424,22 @@ void Odm25dMeshing::parseArguments(int argc, char **argv) { } } -void Odm25dMeshing::loadPointCloud(){ - PlyInterpreter interpreter(points); - - std::ifstream in(inputFile); - if (!in || !CGAL::read_ply_custom_points (in, interpreter, Kernel())){ - throw Odm25dMeshingException( - "Error when reading points and normals from:\n" + inputFile + "\n"); - } - - flipFaces = interpreter.flip_faces(); - - log << "Successfully loaded " << points.size() << " points from file\n"; -} - -void Odm25dMeshing::buildMesh(){ - const unsigned int NEIGHBORS = 24; - - size_t pointCount = points.size(); - size_t pointCountBeforeOutlierRemoval = pointCount; - - log << "Removing outliers... "; - - points.erase(CGAL::remove_outliers(points.begin(), points.end(), - NEIGHBORS, - outliersRemovalPercentage), - points.end()); - std::vector(points).swap(points); - pointCount = points.size(); - - log << "removed " << pointCountBeforeOutlierRemoval - pointCount << " points\n"; - - log << "Computing average spacing... "; - - FT avgSpacing = CGAL::compute_average_spacing( - points.begin(), - points.end(), - NEIGHBORS); - - log << avgSpacing << "\n"; - - log << "Grid Z sampling... "; - - size_t pointCountBeforeGridSampling = pointCount; - - double gridStep = avgSpacing / 2; - Kernel::Iso_cuboid_3 bbox = CGAL::bounding_box(points.begin(), points.end()); - Vector3 boxDiag = bbox.max() - bbox.min(); - - int gridWidth = 1 + static_cast(boxDiag.x() / gridStep + 0.5); - int gridHeight = 1 + static_cast(boxDiag.y() / gridStep + 0.5); - - #define KEY(i, j) (i * gridWidth + j) - - std::unordered_map grid; - - for (size_t c = 0; c < pointCount; c++){ - const Point3 &p = points[c]; - Vector3 relativePos = p - bbox.min(); - int i = static_cast((relativePos.x() / gridStep + 0.5)); - int j = static_cast((relativePos.y() / gridStep + 0.5)); - - if ((i >= 0 && i < gridWidth) && (j >= 0 && j < gridHeight)){ - int key = KEY(i, j); - - if (grid.find(key) == grid.end()){ - grid[key] = p; - }else if ((!flipFaces && p.z() > grid[key].z()) || (flipFaces && p.z() < grid[key].z())){ - grid[key] = p; - } - } - } - - std::vector gridPoints; - for ( auto it = grid.begin(); it != grid.end(); ++it ){ - gridPoints.push_back(it->second); - } - - pointCount = gridPoints.size(); - log << "sampled " << (pointCountBeforeGridSampling - pointCount) << " points\n"; - - const double RETAIN_PERCENTAGE = std::min(80., 100. * static_cast(maxVertexCount) / 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"; - - CGAL::wlop_simplify_and_regularize_point_set( - gridPoints.begin(), - gridPoints.end(), - std::back_inserter(simplifiedPoints), - RETAIN_PERCENTAGE, - 8 * avgSpacing, - wlopIterations, - true); - - pointCount = simplifiedPoints.size(); - - if (pointCount < 3){ - throw Odm25dMeshingException("Not enough points"); - } - - log << "Vertex count is " << pointCount << "\n"; - - typedef CDT::Point cgalPoint; - typedef CDT::Vertex_circulator Vertex_circulator; - - std::vector< std::pair > pts; - try{ - pts.reserve(pointCount); - } catch (const std::bad_alloc&){ - throw Odm25dMeshingException("Not enough memory"); - } - - for (size_t i = 0; i < pointCount; ++i){ - pts.push_back(std::make_pair(cgalPoint(simplifiedPoints[i].x(), simplifiedPoints[i].y()), i)); - } - - log << "Computing delaunay triangulation... "; - - CDT cdt; - cdt.insert(pts.begin(), pts.end()); - - unsigned int numberOfTriangles = static_cast(cdt.number_of_faces()); - unsigned int triIndexes = cdt.number_of_faces()*3; - - if (numberOfTriangles == 0) throw Odm25dMeshingException("No triangles in resulting mesh"); - - log << numberOfTriangles << " triangles\n"; - - std::vector vertices; - std::vector vertexIndices; - - try{ - vertices.reserve(pointCount); - vertexIndices.reserve(triIndexes); - } catch (const std::bad_alloc&){ - throw Odm25dMeshingException("Not enough memory"); - } - - - for (size_t i = 0; i < pointCount; ++i){ - vertices.push_back(simplifiedPoints[i].x()); - vertices.push_back(simplifiedPoints[i].y()); - vertices.push_back(simplifiedPoints[i].z()); - } - - for (CDT::Face_iterator face = cdt.faces_begin(); face != cdt.faces_end(); ++face) { - if (flipFaces){ - vertexIndices.push_back(face->vertex(2)->info()); - vertexIndices.push_back(face->vertex(1)->info()); - vertexIndices.push_back(face->vertex(0)->info()); - }else{ - vertexIndices.push_back(face->vertex(0)->info()); - vertexIndices.push_back(face->vertex(1)->info()); - vertexIndices.push_back(face->vertex(2)->info()); - } - } - - log << "Removing spikes... "; - - const float THRESHOLD = avgSpacing; - std::vector heights; - unsigned int spikesRemoved = 0; - - for (CDT::Vertex_iterator vertex = cdt.vertices_begin(); vertex != cdt.vertices_end(); ++vertex){ - // Check if the height between this vertex and its - // incident vertices is greater than THRESHOLD - Vertex_circulator vc = cdt.incident_vertices(vertex), done(vc); - - if (vc != 0){ - float height = vertices[vertex->info() * 3 + 2]; - int threshold_over_count = 0; - int vertexCount = 0; - - do{ - if (cdt.is_infinite(vc)) continue; - - float ivHeight = vertices[vc->info() * 3 + 2]; - - if (fabs(height - ivHeight) > THRESHOLD){ - threshold_over_count++; - heights.push_back(ivHeight); - } - - vertexCount++; - }while(++vc != done); - - if (vertexCount == threshold_over_count){ - // Replace the height of the vertex by the median height - // of its incident vertices - std::sort(heights.begin(), heights.end()); - - vertices[vertex->info() * 3 + 2] = heights[heights.size() / 2]; - - spikesRemoved++; - } - - heights.clear(); - } - } - - log << "removed " << spikesRemoved << " spikes\n"; - - log << "Building polyhedron... "; - - Polyhedron poly; - PolyhedronBuilder builder(vertices, vertexIndices); - poly.delegate( builder ); - - log << "done\n"; - - log << "Refining... "; - - typedef Polyhedron::Vertex_handle Vertex_handle; - std::vector new_facets; - std::vector new_vertices; - CGAL::Polygon_mesh_processing::refine(poly, - faces(poly), - std::back_inserter(new_facets), - std::back_inserter(new_vertices), - CGAL::Polygon_mesh_processing::parameters::density_control_factor(2.)); - - log << "added " << new_vertices.size() << " new vertices\n"; - -// log << "Edge collapsing... "; -// -// SMS::Count_stop_predicate stop(maxVertexCount * 3); -// int redgesRemoved = SMS::edge_collapse(poly, stop, -// CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index, poly)) -// .halfedge_index_map (get(CGAL::halfedge_external_index, poly)) -// .get_cost (SMS::Edge_length_cost ()) -// .get_placement(SMS::Midpoint_placement()) -// ); -// -// log << redgesRemoved << " edges removed.\n"; - - log << "Final vertex count is " << poly.size_of_vertices() << "\n"; - - log << "Saving mesh to file.\n"; - - typedef typename Polyhedron::Vertex_const_iterator VCI; - typedef typename Polyhedron::Facet_const_iterator FCI; - typedef typename Polyhedron::Halfedge_around_facet_const_circulator HFCC; - - std::filebuf fb; - fb.open(outputFile, std::ios::out); - std::ostream os(&fb); - - os << "ply\n" - << "format ascii 1.0\n" - << "element vertex " << poly.size_of_vertices() << "\n" - << "property float x\n" - << "property float y\n" - << "property float z\n" - << "element face " << poly.size_of_facets() << "\n" - << "property list uchar int vertex_index\n" - << "end_header\n"; - - for (auto it = poly.vertices_begin(); it != poly.vertices_end(); it++){ - os << it->point().x() << " " << it->point().y() << " " << it->point().z() << std::endl; - } - - typedef CGAL::Inverse_index Index; - Index index(poly.vertices_begin(), poly.vertices_end()); - - for( FCI fi = poly.facets_begin(); fi != poly.facets_end(); ++fi) { - HFCC hc = fi->facet_begin(); - HFCC hc_end = hc; - - os << circulator_size(hc) << " "; - do { - os << index[VCI(hc->vertex())] << " "; - ++hc; - } while( hc != hc_end); - - os << "\n"; - } - - fb.close(); - - log << "Successfully wrote mesh to: " << outputFile << "\n"; -} - void Odm25dMeshing::printHelp() { bool printInCoutPop = log.isPrintingInCout(); log.setIsPrintingInCout(true); log << "Usage: odm_25dmeshing -inputFile [plyFile] [optional-parameters]\n"; - log << "Create a 2.5D mesh from an oriented point cloud (points with normals) using a constrained delaunay triangulation. " + log << "Create a 2.5D mesh from a point cloud. " << "The program requires a path to an input PLY point cloud file, all other input parameters are optional.\n\n"; log << " -inputFile to PLY point cloud\n" << " -outputFile where the output PLY 2.5D mesh should be saved (default: " << outputFile << ")\n" + << " -outputDsmFile Optionally output the Digital Surface Model (DSM) computed for generating the mesh. (default: " << outputDsmFile << ")\n" << " -logFile log file path (default: " << logFilePath << ")\n" << " -verbose whether to print verbose output (default: " << (printInCoutPop ? "true" : "false") << ")\n" << " -maxVertexCount <0 - N> Maximum number of vertices in the output mesh. The mesh might have fewer vertices, but will not exceed this limit. (default: " << maxVertexCount << ")\n" - << " -wlopIterations <1 - 1000> Iterations of the Weighted Locally Optimal Projection (WLOP) simplification algorithm. Higher values take longer but produce a smoother mesh. (default: " << wlopIterations << ")\n" + << " -neighbors <1 - 1000> Number of nearest neighbors to consider when doing shepard's interpolation and outlier removal. Higher values lead to smoother meshes but take longer to process. (default: " << neighbors << ")\n" + << " -resolution <0 - N> Size of the interpolated digital surface model (DSM) used for deriving the 2.5D mesh, expressed in pixels per meter unit. When set to zero, the program automatically attempts to find a good value based on the point cloud extent and target vertex count. (default: " << resolution << ")\n" << "\n"; diff --git a/modules/odm_25dmeshing/src/Odm25dMeshing.hpp b/modules/odm_25dmeshing/src/Odm25dMeshing.hpp index 4ebf9df3..bc1ddd6c 100644 --- a/modules/odm_25dmeshing/src/Odm25dMeshing.hpp +++ b/modules/odm_25dmeshing/src/Odm25dMeshing.hpp @@ -1,25 +1,48 @@ #pragma once -// STL -#include -#include -#include -#include +//#define SUPPORTDEBUGWINDOW 1 -#include -#include -#include -#include -#include -#include -#include -#include -#include +#ifdef SUPPORTDEBUGWINDOW + +#include +#include +#include +#include +#include +#include +#include +#include + +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// For compatibility with new VTK generic data arrays +#ifdef vtkGenericDataArray_h +#define InsertNextTupleValue InsertNextTypedTuple +#endif + +#include +#include "tinyply.h" +using namespace tinyply; -#include "CGAL.hpp" #include "Logger.hpp" -#include "PlyInterpreter.hpp" -#include "PolyhedronBuilder.hpp" class Odm25dMeshing { public: @@ -46,31 +69,26 @@ private: */ void parseArguments(int argc, char** argv); - /*! - * \brief loadPointCloud Loads a PLY file with points and normals from file. - */ - void loadPointCloud(); - - /*! - * \brief loadPointCloud Builds a 2.5D mesh from loaded points - */ - void buildMesh(); - /*! * \brief printHelp Prints help, explaining usage. Can be shown by calling the program with argument: "-help". */ void printHelp(); + void loadPointCloud(); + void buildMesh(); + Logger log; std::string inputFile = ""; std::string outputFile = "odm_25dmesh.ply"; std::string logFilePath = "odm_25dmeshing_log.txt"; - unsigned int maxVertexCount = 100000; - double outliersRemovalPercentage = 2; - unsigned int wlopIterations = 35; - std::vector points; - bool flipFaces = false; + int maxVertexCount = 100000; + double resolution = 0; + unsigned int neighbors = 24; + std::string outputDsmFile = ""; + bool showDebugWindow = false; + + vtkSmartPointer points = vtkSmartPointer::New(); }; class Odm25dMeshingException: public std::exception { diff --git a/modules/odm_25dmeshing/src/PlyInterpreter.cpp b/modules/odm_25dmeshing/src/PlyInterpreter.cpp deleted file mode 100644 index f8a07488..00000000 --- a/modules/odm_25dmeshing/src/PlyInterpreter.cpp +++ /dev/null @@ -1,40 +0,0 @@ -#include "PlyInterpreter.hpp" - -// Init and test if input file contains the right properties -bool PlyInterpreter::is_applicable(CGAL::Ply_reader& reader) { - return reader.does_tag_exist ("x") - && reader.does_tag_exist ("y") - && reader.does_tag_exist ("z") - && reader.does_tag_exist ("nx") - && reader.does_tag_exist ("ny") - && reader.does_tag_exist ("nz"); -} - -// Describes how to process one line (= one point object) -void PlyInterpreter::process_line(CGAL::Ply_reader& reader) { - FT x = (FT)0., y = (FT)0., z = (FT)0., - nx = (FT)0., ny = (FT)0., nz = (FT)0.; - - reader.assign (x, "x"); - reader.assign (y, "y"); - reader.assign (z, "z"); - reader.assign (nx, "nx"); - reader.assign (ny, "ny"); - reader.assign (nz, "nz"); - - Point3 p(x, y, z); -// Vector3 n(nx, ny, nz); - - if (nz >= 0 && zNormalsDirectionCount < std::numeric_limits::max()){ - zNormalsDirectionCount++; - }else if (nz < 0 && zNormalsDirectionCount > std::numeric_limits::min()){ - zNormalsDirectionCount--; - } - -// points.push_back(std::make_pair(p, n)); - points.push_back(p); -} - -bool PlyInterpreter::flip_faces(){ - return zNormalsDirectionCount < 0; -} diff --git a/modules/odm_25dmeshing/src/PlyInterpreter.hpp b/modules/odm_25dmeshing/src/PlyInterpreter.hpp deleted file mode 100644 index eca6bcb7..00000000 --- a/modules/odm_25dmeshing/src/PlyInterpreter.hpp +++ /dev/null @@ -1,27 +0,0 @@ -#pragma once - -#include -#include -#include -#include - -#include -#include - -#include "CGAL.hpp" - -// points, normals -//typedef std::pair Pwn; - -class PlyInterpreter { - std::vector& points; - long zNormalsDirectionCount; - - public: - PlyInterpreter (std::vector& points) - : points (points), zNormalsDirectionCount(0) - { } - bool is_applicable (CGAL::Ply_reader& reader); - void process_line (CGAL::Ply_reader& reader); - bool flip_faces(); -}; diff --git a/modules/odm_25dmeshing/src/PolyhedronBuilder.cpp b/modules/odm_25dmeshing/src/PolyhedronBuilder.cpp deleted file mode 100644 index 69b06d44..00000000 --- a/modules/odm_25dmeshing/src/PolyhedronBuilder.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#include "PolyhedronBuilder.hpp" - diff --git a/modules/odm_25dmeshing/src/PolyhedronBuilder.hpp b/modules/odm_25dmeshing/src/PolyhedronBuilder.hpp deleted file mode 100644 index f68dbe0c..00000000 --- a/modules/odm_25dmeshing/src/PolyhedronBuilder.hpp +++ /dev/null @@ -1,43 +0,0 @@ -#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/modules/odm_25dmeshing/src/tinyply.cpp b/modules/odm_25dmeshing/src/tinyply.cpp new file mode 100644 index 00000000..8f0d78e3 --- /dev/null +++ b/modules/odm_25dmeshing/src/tinyply.cpp @@ -0,0 +1,621 @@ +// This software is in the public domain. Where that dedication is not +// recognized, you are granted a perpetual, irrevocable license to copy, +// distribute, and modify this file as you see fit. +// Authored in 2015 by Dimitri Diakopoulos (http://www.dimitridiakopoulos.com) +// https://github.com/ddiakopoulos/tinyply +// Version 2.0 + +#include "tinyply.h" +#include +#include +#include +#include +#include + +using namespace tinyply; +using namespace std; + +////////////////// +// Endian Utils // +////////////////// + +template T endian_swap(const T & v) { return v; } +template<> inline uint16_t endian_swap(const uint16_t & v) { return (v << 8) | (v >> 8); } +template<> inline uint32_t endian_swap(const uint32_t & v) { return (v << 24) | ((v << 8) & 0x00ff0000) | ((v >> 8) & 0x0000ff00) | (v >> 24); } +template<> inline uint64_t endian_swap(const uint64_t & v) +{ + return (((v & 0x00000000000000ffLL) << 56) | + ((v & 0x000000000000ff00LL) << 40) | + ((v & 0x0000000000ff0000LL) << 24) | + ((v & 0x00000000ff000000LL) << 8) | + ((v & 0x000000ff00000000LL) >> 8) | + ((v & 0x0000ff0000000000LL) >> 24) | + ((v & 0x00ff000000000000LL) >> 40) | + ((v & 0xff00000000000000LL) >> 56)); +} +template<> inline int16_t endian_swap(const int16_t & v) { uint16_t r = endian_swap(*(uint16_t*)&v); return *(int16_t*)&r; } +template<> inline int32_t endian_swap(const int32_t & v) { uint32_t r = endian_swap(*(uint32_t*)&v); return *(int32_t*)&r; } +template<> inline int64_t endian_swap(const int64_t & v) { uint64_t r = endian_swap(*(uint64_t*)&v); return *(int64_t*)&r; } +inline float endian_swap_float(const uint32_t & v) { union { float f; uint32_t i; }; i = endian_swap(v); return f; } +inline double endian_swap_double(const uint64_t & v) { union { double d; uint64_t i; }; i = endian_swap(v); return d; } + +///////////////////////////// +// Internal Implementation // +///////////////////////////// + +inline Type property_type_from_string(const std::string & t) +{ + if (t == "int8" || t == "char") return Type::INT8; + else if (t == "uint8" || t == "uchar") return Type::UINT8; + else if (t == "int16" || t == "short") return Type::INT16; + else if (t == "uint16" || t == "ushort") return Type::UINT16; + else if (t == "int32" || t == "int") return Type::INT32; + else if (t == "uint32" || t == "uint") return Type::UINT32; + else if (t == "float32" || t == "float") return Type::FLOAT32; + else if (t == "float64" || t == "double") return Type::FLOAT64; + return Type::INVALID; +} + +struct PlyFile::PlyFileImpl +{ + struct PlyCursor + { + size_t byteOffset; + size_t totalSizeBytes; + }; + + struct ParsingHelper + { + std::shared_ptr data; + std::shared_ptr cursor; + }; + + std::map userData; + + bool isBinary = false; + bool isBigEndian = false; + std::vector elements; + std::vector comments; + std::vector objInfo; + + void read(std::istream & is); + void write(std::ostream & os, bool isBinary); + + std::shared_ptr request_properties_from_element(const std::string & elementKey, const std::initializer_list propertyKeys); + void add_properties_to_element(const std::string & elementKey, const std::initializer_list propertyKeys, const Type type, const size_t count, uint8_t * data, const Type listType, const size_t listCount); + + size_t read_property_binary(const Type t, void * dest, size_t & destOffset, std::istream & is); + size_t read_property_ascii(const Type t, void * dest, size_t & destOffset, std::istream & is); + size_t skip_property_binary(const PlyProperty & property, std::istream & is); + size_t skip_property_ascii(const PlyProperty & property, std::istream & is); + + bool parse_header(std::istream & is); + void parse_data(std::istream & is, bool firstPass); + void read_header_format(std::istream & is); + void read_header_element(std::istream & is); + void read_header_property(std::istream & is); + void read_header_text(std::string line, std::vector & place, int erase = 0); + + void write_header(std::ostream & os); + void write_ascii_internal(std::ostream & os); + void write_binary_internal(std::ostream & os); + void write_property_ascii(Type t, std::ostream & os, uint8_t * src, size_t & srcOffset); + void write_property_binary(Type t, std::ostream & os, uint8_t * src, size_t & srcOffset); +}; + +////////////////// +// PLY Property // +////////////////// + +PlyProperty::PlyProperty(std::istream & is) : isList(false) +{ + std::string type; + is >> type; + if (type == "list") + { + std::string countType; + is >> countType >> type; + listType = property_type_from_string(countType); + isList = true; + } + propertyType = property_type_from_string(type); + is >> name; +} + +///////////////// +// PLY Element // +///////////////// + +PlyElement::PlyElement(std::istream & is) +{ + is >> name >> size; +} + +/////////// +// Utils // +/////////// + +std::string make_key(const std::string & a, const std::string & b) +{ + return (a + "-" + b); +} + +template void ply_cast(void * dest, const char * src, bool be) +{ + *(static_cast(dest)) = (be) ? endian_swap(*(reinterpret_cast(src))) : *(reinterpret_cast(src)); +} + +template void ply_cast_float(void * dest, const char * src, bool be) +{ + *(static_cast(dest)) = (be) ? endian_swap_float(*(reinterpret_cast(src))) : *(reinterpret_cast(src)); +} + +template void ply_cast_double(void * dest, const char * src, bool be) +{ + *(static_cast(dest)) = (be) ? endian_swap_double(*(reinterpret_cast(src))) : *(reinterpret_cast(src)); +} + +template T ply_read_ascii(std::istream & is) +{ + T data; + is >> data; + return data; +} + +template void ply_cast_ascii(void * dest, std::istream & is) +{ + *(static_cast(dest)) = ply_read_ascii(is); +} + +size_t find_element(const std::string & key, const std::vector & list) +{ + for (size_t i = 0; i < list.size(); i++) if (list[i].name == key) return i; + return -1; +} + +size_t find_property(const std::string & key, const std::vector & list) +{ + for (size_t i = 0; i < list.size(); ++i) if (list[i].name == key) return i; + return -1; +} + +////////////// +// PLY File // +////////////// + +bool PlyFile::PlyFileImpl::parse_header(std::istream & is) +{ + std::string line; + while (std::getline(is, line)) + { + std::istringstream ls(line); + std::string token; + ls >> token; + if (token == "ply" || token == "PLY" || token == "") continue; + else if (token == "comment") read_header_text(line, comments, 8); + else if (token == "format") read_header_format(ls); + else if (token == "element") read_header_element(ls); + else if (token == "property") read_header_property(ls); + else if (token == "obj_info") read_header_text(line, objInfo, 9); + else if (token == "end_header") break; + else return false; + } + return true; +} + +void PlyFile::PlyFileImpl::read_header_text(std::string line, std::vector& place, int erase) +{ + place.push_back((erase > 0) ? line.erase(0, erase) : line); +} + +void PlyFile::PlyFileImpl::read_header_format(std::istream & is) +{ + std::string s; + (is >> s); + if (s == "binary_little_endian") isBinary = true; + else if (s == "binary_big_endian") isBinary = isBigEndian = true; +} + +void PlyFile::PlyFileImpl::read_header_element(std::istream & is) +{ + elements.emplace_back(is); +} + +void PlyFile::PlyFileImpl::read_header_property(std::istream & is) +{ + elements.back().properties.emplace_back(is); +} + +size_t PlyFile::PlyFileImpl::skip_property_binary(const PlyProperty & p, std::istream & is) +{ + static std::vector skip(PropertyTable[p.propertyType].stride); + if (p.isList) + { + size_t listSize = 0; + size_t dummyCount = 0; + read_property_binary(p.listType, &listSize, dummyCount, is); + for (size_t i = 0; i < listSize; ++i) is.read(skip.data(), PropertyTable[p.propertyType].stride); + return listSize * PropertyTable[p.propertyType].stride; // in bytes + } + else + { + is.read(skip.data(), PropertyTable[p.propertyType].stride); + return PropertyTable[p.propertyType].stride; + } +} + +size_t PlyFile::PlyFileImpl::skip_property_ascii(const PlyProperty & p, std::istream & is) +{ + std::string skip; + if (p.isList) + { + size_t listSize = 0; + size_t dummyCount = 0; + read_property_ascii(p.listType, &listSize, dummyCount, is); + for (size_t i = 0; i < listSize; ++i) is >> skip; + return listSize * PropertyTable[p.propertyType].stride; // in bytes + } + else + { + is >> skip; + return PropertyTable[p.propertyType].stride; + } +} + +size_t PlyFile::PlyFileImpl::read_property_binary(const Type t, void * dest, size_t & destOffset, std::istream & is) +{ + destOffset += PropertyTable[t].stride; + + std::vector src(PropertyTable[t].stride); + is.read(src.data(), PropertyTable[t].stride); + + switch (t) + { + case Type::INT8: ply_cast(dest, src.data(), isBigEndian); break; + case Type::UINT8: ply_cast(dest, src.data(), isBigEndian); break; + case Type::INT16: ply_cast(dest, src.data(), isBigEndian); break; + case Type::UINT16: ply_cast(dest, src.data(), isBigEndian); break; + case Type::INT32: ply_cast(dest, src.data(), isBigEndian); break; + case Type::UINT32: ply_cast(dest, src.data(), isBigEndian); break; + case Type::FLOAT32: ply_cast_float(dest, src.data(), isBigEndian); break; + case Type::FLOAT64: ply_cast_double(dest, src.data(), isBigEndian); break; + case Type::INVALID: throw std::invalid_argument("invalid ply property"); + } + + return PropertyTable[t].stride; +} + +size_t PlyFile::PlyFileImpl::read_property_ascii(const Type t, void * dest, size_t & destOffset, std::istream & is) +{ + destOffset += PropertyTable[t].stride; + + switch (t) + { + case Type::INT8: *((int8_t *)dest) = ply_read_ascii(is); break; + case Type::UINT8: *((uint8_t *)dest) = ply_read_ascii(is); break; + case Type::INT16: ply_cast_ascii(dest, is); break; + case Type::UINT16: ply_cast_ascii(dest, is); break; + case Type::INT32: ply_cast_ascii(dest, is); break; + case Type::UINT32: ply_cast_ascii(dest, is); break; + case Type::FLOAT32: ply_cast_ascii(dest, is); break; + case Type::FLOAT64: ply_cast_ascii(dest, is); break; + case Type::INVALID: throw std::invalid_argument("invalid ply property"); + } + return PropertyTable[t].stride; +} + +void PlyFile::PlyFileImpl::write_property_ascii(Type t, std::ostream & os, uint8_t * src, size_t & srcOffset) +{ + switch (t) + { + case Type::INT8: os << static_cast(*reinterpret_cast(src)); break; + case Type::UINT8: os << static_cast(*reinterpret_cast(src)); break; + case Type::INT16: os << *reinterpret_cast(src); break; + case Type::UINT16: os << *reinterpret_cast(src); break; + case Type::INT32: os << *reinterpret_cast(src); break; + case Type::UINT32: os << *reinterpret_cast(src); break; + case Type::FLOAT32: os << *reinterpret_cast(src); break; + case Type::FLOAT64: os << *reinterpret_cast(src); break; + case Type::INVALID: throw std::invalid_argument("invalid ply property"); + } + os << " "; + srcOffset += PropertyTable[t].stride; +} + +void PlyFile::PlyFileImpl::write_property_binary(Type t, std::ostream & os, uint8_t * src, size_t & srcOffset) +{ + os.write((char *)src, PropertyTable[t].stride); + srcOffset += PropertyTable[t].stride; +} + +void PlyFile::PlyFileImpl::read(std::istream & is) +{ + // Parse but only get the data size + parse_data(is, true); + + std::vector> buffers; + for (auto & entry : userData) buffers.push_back(entry.second.data); + + // Since group-requested properties share the same cursor, we need to find unique cursors so we only allocate once + std::sort(buffers.begin(), buffers.end()); + buffers.erase(std::unique(buffers.begin(), buffers.end()), buffers.end()); + + // Not great, but since we sorted by ptrs on PlyData, need to remap back onto its cursor in the userData table + for (auto & b : buffers) + { + for (auto & entry : userData) + { + if (entry.second.data == b && b->buffer.get() == nullptr) + { + b->buffer = Buffer(entry.second.cursor->totalSizeBytes); + } + } + } + + // Populate the data + parse_data(is, false); +} + +void PlyFile::PlyFileImpl::write(std::ostream & os, bool _isBinary) +{ + if (_isBinary) write_binary_internal(os); + else write_ascii_internal(os); +} + +void PlyFile::PlyFileImpl::write_binary_internal(std::ostream & os) +{ + isBinary = true; + write_header(os); + + for (auto & e : elements) + { + for (size_t i = 0; i < e.size; ++i) + { + for (auto & p : e.properties) + { + auto & helper = userData[make_key(e.name, p.name)]; + if (p.isList) + { + uint8_t listSize[4] = {0, 0, 0, 0}; + std::memcpy(listSize, &p.listCount, sizeof(uint32_t)); + size_t dummyCount = 0; + write_property_binary(p.listType, os, listSize, dummyCount); + for (int j = 0; j < p.listCount; ++j) + { + write_property_binary(p.propertyType, os, (helper.data->buffer.get() + helper.cursor->byteOffset), helper.cursor->byteOffset); + } + } + else + { + write_property_binary(p.propertyType, os, (helper.data->buffer.get() + helper.cursor->byteOffset), helper.cursor->byteOffset); + } + } + } + } +} + +void PlyFile::PlyFileImpl::write_ascii_internal(std::ostream & os) +{ + write_header(os); + + for (auto & e : elements) + { + for (size_t i = 0; i < e.size; ++i) + { + for (auto & p : e.properties) + { + auto & helper = userData[make_key(e.name, p.name)]; + if (p.isList) + { + os << p.listCount << " "; + for (int j = 0; j < p.listCount; ++j) + { + write_property_ascii(p.propertyType, os, (helper.data->buffer.get() + helper.cursor->byteOffset), helper.cursor->byteOffset); + } + } + else + { + write_property_ascii(p.propertyType, os, (helper.data->buffer.get() + helper.cursor->byteOffset), helper.cursor->byteOffset); + } + } + os << "\n"; + } + } +} + +void PlyFile::PlyFileImpl::write_header(std::ostream & os) +{ + const std::locale & fixLoc = std::locale("C"); + os.imbue(fixLoc); + + os << "ply\n"; + if (isBinary) os << ((isBigEndian) ? "format binary_big_endian 1.0" : "format binary_little_endian 1.0") << "\n"; + else os << "format ascii 1.0\n"; + + for (const auto & comment : comments) os << "comment " << comment << "\n"; + + for (auto & e : elements) + { + os << "element " << e.name << " " << e.size << "\n"; + for (const auto & p : e.properties) + { + if (p.isList) + { + os << "property list " << PropertyTable[p.listType].str << " " + << PropertyTable[p.propertyType].str << " " << p.name << "\n"; + } + else + { + os << "property " << PropertyTable[p.propertyType].str << " " << p.name << "\n"; + } + } + } + os << "end_header\n"; +} + +// Returns the size (in bytes) +std::shared_ptr PlyFile::PlyFileImpl::request_properties_from_element(const std::string & elementKey, const std::initializer_list propertyKeys) +{ + // All requested properties in the userDataTable share the same cursor (thrown into the same flat array) + ParsingHelper helper; + helper.data = std::make_shared(); + helper.data->count = 0; + helper.data->t = Type::INVALID; + helper.cursor = std::make_shared(); + helper.cursor->byteOffset = 0; + helper.cursor->totalSizeBytes = 0; + + if (elements.size() == 0) throw std::runtime_error("parsed header had no elements defined. malformed file?"); + if (!propertyKeys.size()) throw std::invalid_argument("`propertyKeys` argument is empty"); + if (elementKey.size() == 0) throw std::invalid_argument("`elementKey` argument is empty"); + + const int elementIndex = find_element(elementKey, elements); + + // Sanity check if the user requested element is in the pre-parsed header + if (elementIndex >= 0) + { + // We found the element + const PlyElement & element = elements[elementIndex]; + + helper.data->count = element.size; + + // Find each of the keys + for (auto key : propertyKeys) + { + const int propertyIndex = find_property(key, element.properties); + + if (propertyIndex >= 0) + { + // We found the property + const PlyProperty & property = element.properties[propertyIndex]; + + helper.data->t = property.propertyType; // hmm.... + + auto result = userData.insert(std::pair(make_key(element.name, property.name), helper)); + if (result.second == false) throw std::invalid_argument("element-property key has already been requested: " + make_key(element.name, property.name)); + } + else throw std::invalid_argument("one of the property keys was not found in the header: " + key); + } + } + else throw std::invalid_argument("the element key was not found in the header: " + elementKey); + + return helper.data; +} + +void PlyFile::PlyFileImpl::add_properties_to_element(const std::string & elementKey, const std::initializer_list propertyKeys, const Type type, const size_t count, uint8_t * data, const Type listType, const size_t listCount) +{ + ParsingHelper helper; + helper.data = std::make_shared(); + helper.data->count = count; + helper.data->t = type; + helper.data->buffer = Buffer(data); + helper.cursor = std::make_shared(); + helper.cursor->byteOffset = 0; + helper.cursor->totalSizeBytes = 0; + + auto create_property_on_element = [&](PlyElement & e) + { + for (auto key : propertyKeys) + { + PlyProperty newProp = (listType == Type::INVALID) ? PlyProperty(type, key) : PlyProperty(listType, type, key, listCount); + /* auto result = */userData.insert(std::pair(make_key(elementKey, key), helper)); + e.properties.push_back(newProp); + } + }; + + int idx = find_element(elementKey, elements); + if (idx >= 0) + { + PlyElement & e = elements[idx]; + create_property_on_element(e); + } + else + { + PlyElement newElement = (listType == Type::INVALID) ? PlyElement(elementKey, count / propertyKeys.size()) : PlyElement(elementKey, count / listCount); + create_property_on_element(newElement); + elements.push_back(newElement); + } +} + +void PlyFile::PlyFileImpl::parse_data(std::istream & is, bool firstPass) +{ + std::function read; + std::function skip; + + const auto start = is.tellg(); + + if (isBinary) + { + read = [&](const Type t, void * dest, size_t & destOffset, std::istream & _is) { return read_property_binary(t, dest, destOffset, _is); }; + skip = [&](const PlyProperty & p, std::istream & _is) { return skip_property_binary(p, _is); }; + } + else + { + read = [&](const Type t, void * dest, size_t & destOffset, std::istream & _is) { return read_property_ascii(t, dest, destOffset, _is); }; + skip = [&](const PlyProperty & p, std::istream & _is) { return skip_property_ascii(p, _is); }; + } + + for (auto & element : elements) + { + for (size_t count = 0; count < element.size; ++count) + { + for (auto & property : element.properties) + { + auto cursorIt = userData.find(make_key(element.name, property.name)); + if (cursorIt != userData.end()) + { + auto & helper = cursorIt->second; + if (!firstPass) + { + if (property.isList) + { + size_t listSize = 0; + size_t dummyCount = 0; + read(property.listType, &listSize, dummyCount, is); + for (size_t i = 0; i < listSize; ++i) + { + read(property.propertyType, (helper.data->buffer.get() + helper.cursor->byteOffset), helper.cursor->byteOffset, is); + } + } + else + { + read(property.propertyType, (helper.data->buffer.get() + helper.cursor->byteOffset), helper.cursor->byteOffset, is); + } + } + else + { + helper.cursor->totalSizeBytes += skip(property, is); + } + } + else + { + skip(property, is); + } + } + } + } + + // Reset istream reader to the beginning + if (firstPass) is.seekg(start, is.beg); +} + +/////////////////////////////////// +// Pass-Through Public Interface // +/////////////////////////////////// + +PlyFile::PlyFile() { impl.reset(new PlyFileImpl()); }; +PlyFile::~PlyFile() { }; +bool PlyFile::parse_header(std::istream & is) { return impl->parse_header(is); } +void PlyFile::read(std::istream & is) { return impl->read(is); } +void PlyFile::write(std::ostream & os, bool isBinary) { return impl->write(os, isBinary); } +std::vector PlyFile::get_elements() const { return impl->elements; } +std::vector & PlyFile::get_comments() { return impl->comments; } +std::vector PlyFile::get_info() const { return impl->objInfo; } +std::shared_ptr PlyFile::request_properties_from_element(const std::string & elementKey, const std::initializer_list propertyKeys) +{ + return impl->request_properties_from_element(elementKey, propertyKeys); +} +void PlyFile::add_properties_to_element(const std::string & elementKey, const std::initializer_list propertyKeys, const Type type, const size_t count, uint8_t * data, const Type listType, const size_t listCount) +{ + return impl->add_properties_to_element(elementKey, propertyKeys, type, count, data, listType, listCount); +} diff --git a/modules/odm_25dmeshing/src/tinyply.h b/modules/odm_25dmeshing/src/tinyply.h new file mode 100644 index 00000000..08fedf6f --- /dev/null +++ b/modules/odm_25dmeshing/src/tinyply.h @@ -0,0 +1,119 @@ +// This software is in the public domain. Where that dedication is not +// recognized, you are granted a perpetual, irrevocable license to copy, +// distribute, and modify this file as you see fit. +// Authored in 2015 by Dimitri Diakopoulos (http://www.dimitridiakopoulos.com) +// https://github.com/ddiakopoulos/tinyply +// Version 2.0 + +#ifndef tinyply_h +#define tinyply_h + +#include +#include +#include +#include +#include +#include + +namespace tinyply +{ + + enum class Type : uint8_t + { + INVALID, + INT8, + UINT8, + INT16, + UINT16, + INT32, + UINT32, + FLOAT32, + FLOAT64 + }; + + struct PropertyInfo + { + int stride; + std::string str; + }; + + static std::map PropertyTable + { + { Type::INT8,{ 1, "char" } }, + { Type::UINT8,{ 1, "uchar" } }, + { Type::INT16,{ 2, "short" } }, + { Type::UINT16,{ 2, "ushort" } }, + { Type::INT32,{ 4, "int" } }, + { Type::UINT32,{ 4, "uint" } }, + { Type::FLOAT32,{ 4, "float" } }, + { Type::FLOAT64,{ 8, "double" } }, + { Type::INVALID,{ 0, "INVALID" } } + }; + + class Buffer + { + uint8_t * alias{ nullptr }; + struct delete_array { void operator()(uint8_t * p) { delete[] p; } }; + std::unique_ptr data; + size_t size; + public: + Buffer() {}; + Buffer(const size_t size) : data(new uint8_t[size], delete_array()), size(size) { alias = data.get(); } // allocating + Buffer(uint8_t * ptr) { alias = ptr; } // non-allocating, fixme: set size? + uint8_t * get() { return alias; } + size_t size_bytes() const { return size; } + }; + + struct PlyData + { + Type t; + size_t count; + Buffer buffer; + }; + + struct PlyProperty + { + PlyProperty(std::istream & is); + PlyProperty(Type type, std::string & _name) : name(_name), propertyType(type) {} + PlyProperty(Type list_type, Type prop_type, std::string & _name, int list_count) : name(_name), propertyType(prop_type), isList(true), listType(list_type), listCount(list_count) {} + std::string name; + Type propertyType; + bool isList{ false }; + Type listType{ Type::INVALID }; + int listCount{ 0 }; + }; + + struct PlyElement + { + PlyElement(std::istream & istream); + PlyElement(const std::string & _name, size_t count) : name(_name), size(count) {} + std::string name; + size_t size; + std::vector properties; + }; + + struct PlyFile + { + struct PlyFileImpl; + std::unique_ptr impl; + + PlyFile(); + ~PlyFile(); + + bool parse_header(std::istream & is); + + void read(std::istream & is); + + void write(std::ostream & os, bool isBinary); + + std::vector get_elements() const; + std::vector & get_comments(); + std::vector get_info() const; + + std::shared_ptr request_properties_from_element(const std::string & elementKey, const std::initializer_list propertyKeys); + void add_properties_to_element(const std::string & elementKey, const std::initializer_list propertyKeys, const Type type, const size_t count, uint8_t * data, const Type listType, const size_t listCount); + }; + +} // namesapce tinyply + +#endif // tinyply_h diff --git a/opendm/config.py b/opendm/config.py index 5e7e91a3..846278a8 100644 --- a/opendm/config.py +++ b/opendm/config.py @@ -104,7 +104,7 @@ def config(): parser.add_argument('--min-num-features', metavar='', - default=4000, + default=8000, type=int, help=('Minimum number of features to extract per image. ' 'More features leads to better results but slower ' @@ -248,23 +248,43 @@ def config(): 'times slightly but helps reduce memory usage. ' 'Default: %(default)s')) - parser.add_argument('--mesh-remove-outliers', - metavar='', - default=2, - type=float, - help=('Percentage of outliers to remove from the point set. Set to 0 to disable. ' + parser.add_argument('--mesh-neighbors', + metavar='', + default=24, + type=int, + help=('Number of neighbors to select when estimating the surface model used to compute the mesh and for statistical outlier removal. Higher values lead to smoother meshes but take longer to process. ' 'Applies to 2.5D mesh only. ' 'Default: %(default)s')) - parser.add_argument('--mesh-wlop-iterations', - metavar='', - default=35, - type=int, - help=('Iterations of the Weighted Locally Optimal Projection (WLOP) simplification algorithm. ' - 'Higher values take longer but produce a smoother mesh. ' + parser.add_argument('--mesh-resolution', + metavar='', + default=0, + type=float, + help=('Size of the interpolated surface model used for deriving the 2.5D mesh, expressed in pixels per meter. ' + 'Higher values work better for complex or urban terrains. ' + 'Lower values work better on flat areas. ' + 'Resolution has no effect on the number of vertices, but high values can severely impact runtime speed and memory usage. ' + 'When set to zero, the program automatically attempts to find a good value based on the point cloud extent and target vertex count. ' 'Applies to 2.5D mesh only. ' 'Default: %(default)s')) + parser.add_argument('--fast-orthophoto', + action='store_true', + default=False, + help='Skips dense reconstruction and 3D model generation. ' + 'It generates an orthophoto directly from the sparse reconstruction. ' + 'If you just need an orthophoto and do not need a full 3D model, turn on this option. ' + 'Experimental.') + + parser.add_argument('--crop', + metavar='', + default=3, + type=float, + help=('Automatically crop image outputs by creating a smooth buffer ' + 'around the dataset boundaries, shrinked by N meters. ' + 'Use 0 to disable cropping. ' + 'Default: %(default)s')) + parser.add_argument('--texturing-data-term', metavar='', default='gmi', @@ -490,4 +510,13 @@ def config(): 'run.py --help` for more information. ') sys.exit(1) + if args.fast_orthophoto: + log.ODM_INFO('Fast orthophoto is turned on, automatically setting --use-25dmesh') + args.use_25dmesh = True + + # Cannot use pmvs + if args.use_pmvs: + log.ODM_INFO('Fast orthophoto is turned on, cannot use pmvs (removing --use-pmvs)') + args.use_pmvs = False + return args diff --git a/opendm/cropper.py b/opendm/cropper.py new file mode 100644 index 00000000..f41b6ccd --- /dev/null +++ b/opendm/cropper.py @@ -0,0 +1,200 @@ +from opendm import context +from opendm import system +from opendm import log +from osgeo import ogr +import json, os + +def run(command): + env_paths = [context.superbuild_bin_path] + return system.run(command, env_paths) + +class Cropper: + def __init__(self, storage_dir, files_prefix = "crop"): + self.storage_dir = storage_dir + self.files_prefix = files_prefix + + def path(self, suffix): + """ + @return a path relative to storage_dir and prefixed with files_prefix + """ + return os.path.join(self.storage_dir, '{}.{}'.format(self.files_prefix, suffix)) + + @staticmethod + def crop(shapefile_path, geotiff_path, gdal_options): + if not os.path.exists(shapefile_path) or not os.path.exists(geotiff_path): + log.ODM_WARNING("Either {} or {} does not exist, will skip cropping.".format(shapefile_path, geotiff_path)) + return geotiff_path + + # Rename original file + # path/to/odm_orthophoto.tif --> path/to/odm_orthophoto.original.tif + + path, filename = os.path.split(geotiff_path) + # path = path/to + # filename = odm_orthophoto.tif + + basename, ext = os.path.splitext(filename) + # basename = odm_orthophoto + # ext = .tif + + original_geotiff = os.path.join(path, "{}.original{}".format(basename, ext)) + os.rename(geotiff_path, original_geotiff) + + try: + kwargs = { + 'shapefile_path': shapefile_path, + 'geotiffInput': original_geotiff, + 'geotiffOutput': geotiff_path, + 'options': ' '.join(map(lambda k: '-co {}={}'.format(k, gdal_options[k]), gdal_options)) + } + + run('gdalwarp -cutline {shapefile_path} ' + '-crop_to_cutline ' + '{options} ' + '{geotiffInput} ' + '{geotiffOutput} '.format(**kwargs)) + + except Exception as e: + log.ODM_WARNING('Something went wrong while cropping: {}'.format(e.message)) + + # Revert rename + os.rename(original_geotiff, geotiff_path) + + return geotiff_path + + def create_bounds_geojson(self, pointcloud_path, buffer_distance = 0): + """ + Compute a buffered polygon around the data extents (not just a bounding box) + of the given point cloud. + + @return filename to GeoJSON containing the polygon + """ + if not os.path.exists(pointcloud_path): + log.ODM_WARNING('Point cloud does not exist, cannot generate shapefile bounds {}'.format(pointcloud_path)) + return '' + + # Do basic outlier removal prior to extracting boundary information + filtered_pointcloud_path = self.path('filtered.las') + + run("pdal translate -i \"{}\" " + "-o \"{}\" " + "decimation outlier range " + "--filters.decimation.step=40 " + "--filters.outlier.method=radius " + "--filters.outlier.radius=20 " + "--filters.outlier.min_k=2 " + "--filters.range.limits='Classification![7:7]'".format(pointcloud_path, filtered_pointcloud_path)) + + if not os.path.exists(filtered_pointcloud_path): + log.ODM_WARNING('Could not filter point cloud, cannot generate shapefile bounds {}'.format(filtered_pointcloud_path)) + return '' + + # Use PDAL to dump boundary information + # then read the information back + + boundary_file_path = self.path('boundary.json') + + run('pdal info --boundary --filters.hexbin.edge_length=1 --filters.hexbin.threshold=0 {0} > {1}'.format(filtered_pointcloud_path, boundary_file_path)) + + pc_geojson_boundary_feature = None + + with open(boundary_file_path, 'r') as f: + json_f = json.loads(f.read()) + pc_geojson_boundary_feature = json_f['boundary']['boundary_json'] + + if pc_geojson_boundary_feature is None: raise RuntimeError("Could not determine point cloud boundaries") + + # Write bounds to GeoJSON + bounds_geojson_path = self.path('bounds.geojson') + with open(bounds_geojson_path, "w") as f: + f.write(json.dumps({ + "type": "FeatureCollection", + "features": [{ + "type": "Feature", + "geometry": pc_geojson_boundary_feature + }] + })) + + # Create a convex hull around the boundary + # as to encompass the entire area (no holes) + driver = ogr.GetDriverByName('GeoJSON') + ds = driver.Open(bounds_geojson_path, 0) # ready-only + layer = ds.GetLayer() + + # Collect all Geometry + geomcol = ogr.Geometry(ogr.wkbGeometryCollection) + for feature in layer: + geomcol.AddGeometry(feature.GetGeometryRef()) + + # Calculate convex hull + convexhull = geomcol.ConvexHull() + + # If buffer distance is specified + # Create two buffers, one shrinked by + # N + 3 and then that buffer expanded by 3 + # so that we get smooth corners. \m/ + BUFFER_SMOOTH_DISTANCE = 3 + + if buffer_distance > 0: + convexhull = convexhull.Buffer(-(buffer_distance + BUFFER_SMOOTH_DISTANCE)) + convexhull = convexhull.Buffer(BUFFER_SMOOTH_DISTANCE) + + # Save to a new file + bounds_geojson_path = self.path('bounds.geojson') + if os.path.exists(bounds_geojson_path): + driver.DeleteDataSource(bounds_geojson_path) + + out_ds = driver.CreateDataSource(bounds_geojson_path) + layer = out_ds.CreateLayer("convexhull", geom_type=ogr.wkbPolygon) + + feature_def = layer.GetLayerDefn() + feature = ogr.Feature(feature_def) + feature.SetGeometry(convexhull) + layer.CreateFeature(feature) + feature = None + + # Save and close data sources + out_ds = ds = None + + # Remove filtered point cloud + if os.path.exists(filtered_pointcloud_path): + os.remove(filtered_pointcloud_path) + + return bounds_geojson_path + + + def create_bounds_shapefile(self, pointcloud_path, buffer_distance = 0): + """ + Compute a buffered polygon around the data extents (not just a bounding box) + of the given point cloud. + + @return filename to Shapefile containing the polygon + """ + if not os.path.exists(pointcloud_path): + log.ODM_WARNING('Point cloud does not exist, cannot generate shapefile bounds {}'.format(pointcloud_path)) + return '' + + bounds_geojson_path = self.create_bounds_geojson(pointcloud_path, buffer_distance) + + summary_file_path = os.path.join(self.storage_dir, '{}.summary.json'.format(self.files_prefix)) + run('pdal info --summary {0} > {1}'.format(pointcloud_path, summary_file_path)) + + pc_proj4 = None + with open(summary_file_path, 'r') as f: + json_f = json.loads(f.read()) + pc_proj4 = json_f['summary']['srs']['proj4'] + + if pc_proj4 is None: raise RuntimeError("Could not determine point cloud proj4 declaration") + + bounds_shapefile_path = os.path.join(self.storage_dir, '{}.bounds.shp'.format(self.files_prefix)) + + # Convert bounds to Shapefile + kwargs = { + 'input': bounds_geojson_path, + 'output': bounds_shapefile_path, + 'proj4': pc_proj4 + } + + run('ogr2ogr -overwrite -a_srs "{proj4}" {output} {input}'.format(**kwargs)) + + return bounds_shapefile_path + diff --git a/opendm/types.py b/opendm/types.py index a981c93d..5ee2728d 100644 --- a/opendm/types.py +++ b/opendm/types.py @@ -221,82 +221,6 @@ class ODM_GeoRef(object): system.run('{bin}/pdal pipeline -i {json} --readers.ply.filename={f_in} ' '--writers.las.filename={f_out}'.format(**kwargs)) - def utm_to_latlon(self, _file, _photo, idx): - - gcp = self.gcps[idx] - - kwargs = {'epsg': self.epsg, - 'file': _file, - 'x': gcp.x + self.utm_east_offset, - 'y': gcp.y + self.utm_north_offset, - 'z': gcp.z} - - latlon = system.run_and_return('echo {x} {y} {z} '.format(**kwargs), - 'gdaltransform -s_srs \"EPSG:{epsg}\" ' - '-t_srs \"EPSG:4326\"'.format(**kwargs)).split() - - # Example: 83d18'16.285"W - # Example: 41d2'11.789"N - # Example: 0.998 - - if len(latlon) == 3: - lon_str, lat_str, alt_str = latlon - elif len(latlon) == 2: - lon_str, lat_str = latlon - alt_str = '' - else: - log.ODM_ERROR('Something went wrong %s' % latlon) - - lat_frac = self.coord_to_fractions(latlon[1], ['N', 'S']) - lon_frac = self.coord_to_fractions(latlon[0], ['E', 'W']) - - # read image metadata - metadata = pyexiv2.ImageMetadata(_photo.path_file) - metadata.read() - - # #set values - # - # # GPS latitude - # key = 'Exif.GPSInfo.GPSLatitude' - # value = lat_frac[0].split(' ') - # log.ODM_DEBUG('lat_frac: %s %s %s' % (value[0], value[1], value[2])) - # metadata[key] = pyexiv2.ExifTag(key, - # [Fraction(value[0]), - # Fraction(value[1]), - # Fraction(value[2])]) - # - # key = 'Exif.GPSInfo.GPSLatitudeRef' - # value = lat_frac[1] - # metadata[key] = pyexiv2.ExifTag(key, value) - # - # # GPS longitude - # key = 'Exif.GPSInfo.GPSLongitude' - # value = lon_frac[0].split(' ') - # metadata[key] = pyexiv2.ExifTag(key, - # [Fraction(value[0]), - # Fraction(value[1]), - # Fraction(value[2])]) - # - # key = 'Exif.GPSInfo.GPSLongitudeRef' - # value = lon_frac[1] - # metadata[key] = pyexiv2.ExifTag(key, value) - # - # # GPS altitude - # altitude = abs(int(float(latlon[2]) * 100)) - # key = 'Exif.GPSInfo.GPSAltitude' - # value = Fraction(altitude, 1) - # metadata[key] = pyexiv2.ExifTag(key, value) - # - # if latlon[2] >= 0: - # altref = '0' - # else: - # altref = '1' - # key = 'Exif.GPSInfo.GPSAltitudeRef' - # metadata[key] = pyexiv2.ExifTag(key, altref) - # - # # write values - # metadata.write() - def parse_coordinate_system(self, _file): """Write attributes to jobOptions from coord file""" # check for coordinate file existence diff --git a/scripts/mvstex.py b/scripts/mvstex.py index a5a59733..47adea24 100644 --- a/scripts/mvstex.py +++ b/scripts/mvstex.py @@ -50,13 +50,24 @@ class ODMMvsTexCell(ecto.Cell): runs = [{ 'out_dir': tree.odm_texturing, - 'model': tree.odm_mesh + 'model': tree.odm_mesh, + 'force_skip_vis_test': False }] + if args.fast_orthophoto: + runs = [] + if args.use_25dmesh: runs += [{ 'out_dir': tree.odm_25dtexturing, - 'model': tree.odm_25dmesh + 'model': tree.odm_25dmesh, + + # We always skip the visibility test when using the 2.5D mesh + # because many faces end up being narrow, and almost perpendicular + # to the ground plane. The visibility test improperly classifies + # them as "not seen" since the test is done on a single triangle vertex, + # and while one vertex might be occluded, the other two might not. + 'force_skip_vis_test': True }] for r in runs: @@ -74,7 +85,7 @@ class ODMMvsTexCell(ecto.Cell): skipHoleFilling = "" keepUnseenFaces = "" - if (self.params.skip_vis_test): + if (self.params.skip_vis_test or r['force_skip_vis_test']): skipGeometricVisibilityTest = "--skip_geometric_visibility_test" if (self.params.skip_glob_seam_leveling): skipGlobalSeamLeveling = "--skip_global_seam_leveling" diff --git a/scripts/odm_app.py b/scripts/odm_app.py index 1bc5aa24..943e40a4 100644 --- a/scripts/odm_app.py +++ b/scripts/odm_app.py @@ -59,8 +59,6 @@ class ODMApp(ecto.BlackBox): oct_tree=p.args.mesh_octree_depth, samples=p.args.mesh_samples, solver=p.args.mesh_solver_divide, - remove_outliers=p.args.mesh_remove_outliers, - wlop_iterations=p.args.mesh_wlop_iterations, verbose=p.args.verbose), 'texturing': ODMMvsTexCell(data_term=p.args.texturing_data_term, outlier_rem_type=p.args.texturing_outlier_removal_type, diff --git a/scripts/odm_dem.py b/scripts/odm_dem.py index 57d8870f..eb30c07a 100644 --- a/scripts/odm_dem.py +++ b/scripts/odm_dem.py @@ -61,46 +61,7 @@ class ODMDEMCell(ecto.Cell): if (args.dtm and not io.file_exists(dtm_output_filename)) or \ (args.dsm and not io.file_exists(dsm_output_filename)) or \ rerun_cell: - - # Extract boundaries and srs of point cloud - summary_file_path = os.path.join(odm_dem_root, 'odm_georeferenced_model.summary.json') - boundary_file_path = os.path.join(odm_dem_root, 'odm_georeferenced_model.boundary.json') - - system.run('pdal info --summary {0} > {1}'.format(tree.odm_georeferencing_model_las, summary_file_path), env_paths) - system.run('pdal info --boundary {0} > {1}'.format(tree.odm_georeferencing_model_las, boundary_file_path), env_paths) - - pc_proj4 = "" - pc_geojson_bounds_feature = None - - with open(summary_file_path, 'r') as f: - json_f = json.loads(f.read()) - pc_proj4 = json_f['summary']['srs']['proj4'] - - with open(boundary_file_path, 'r') as f: - json_f = json.loads(f.read()) - pc_geojson_boundary_feature = json_f['boundary']['boundary_json'] - - # Write bounds to GeoJSON - bounds_geojson_path = os.path.join(odm_dem_root, 'odm_georeferenced_model.bounds.geojson') - with open(bounds_geojson_path, "w") as f: - f.write(json.dumps({ - "type": "FeatureCollection", - "features": [{ - "type": "Feature", - "geometry": pc_geojson_boundary_feature - }] - })) - - bounds_shapefile_path = os.path.join(odm_dem_root, 'bounds.shp') - - # Convert bounds to Shapefile - kwargs = { - 'input': bounds_geojson_path, - 'output': bounds_shapefile_path, - 'proj4': pc_proj4 - } - system.run('ogr2ogr -overwrite -a_srs "{proj4}" {output} {input}'.format(**kwargs)) - + # Process with lidar2dems terrain_params_map = { 'flatnonforest': (1, 3), @@ -115,12 +76,17 @@ class ODMDEMCell(ecto.Cell): 'slope': terrain_params[0], 'cellsize': terrain_params[1], 'outdir': odm_dem_root, - 'site': bounds_shapefile_path + 'site': '' } + if args.crop > 0: + bounds_shapefile_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.shp') + if os.path.exists(bounds_shapefile_path): + kwargs['site'] = '-s {}'.format(bounds_shapefile_path) + l2d_params = '--slope {slope} --cellsize {cellsize} ' \ '{verbose} ' \ - '-o -s {site} ' \ + '-o {site} ' \ '--outdir {outdir}'.format(**kwargs) approximate = '--approximate' if args.dem_approximate else '' @@ -135,7 +101,8 @@ class ODMDEMCell(ecto.Cell): args.dem_initial_distance, tree.odm_georeferencing), env_paths) else: log.ODM_INFO("Will skip classification, only DSM is needed") - copyfile(tree.odm_georeferencing_model_las, os.path.join(odm_dem_root, 'bounds-0_l2d_s{slope}c{cellsize}.las'.format(**kwargs))) + l2d_classified_pattern = 'odm_georeferenced_model.bounds-0_l2d_s{slope}c{cellsize}.las' if args.crop > 0 else 'l2d_s{slope}c{cellsize}.las' + copyfile(tree.odm_georeferencing_model_las, os.path.join(odm_dem_root, l2d_classified_pattern.format(**kwargs))) products = [] if args.dsm: products.append('dsm') @@ -168,9 +135,11 @@ class ODMDEMCell(ecto.Cell): # Rename final output if product == 'dsm': - os.rename(os.path.join(odm_dem_root, 'bounds-0_dsm.idw.tif'), dsm_output_filename) + dsm_pattern = 'odm_georeferenced_model.bounds-0_dsm.idw.tif' if args.crop > 0 else 'dsm.idw.tif' + os.rename(os.path.join(odm_dem_root, dsm_pattern), dsm_output_filename) elif product == 'dtm': - os.rename(os.path.join(odm_dem_root, 'bounds-0_dtm.idw.tif'), dtm_output_filename) + dtm_pattern = 'odm_georeferenced_model.bounds-0_dsm.idw.tif' if args.crop > 0 else 'dtm.idw.tif' + os.rename(os.path.join(odm_dem_root, dtm_pattern), dtm_output_filename) else: log.ODM_WARNING('Found existing outputs in: %s' % odm_dem_root) diff --git a/scripts/odm_georeferencing.py b/scripts/odm_georeferencing.py index 0ca6e167..a7e2cf9b 100644 --- a/scripts/odm_georeferencing.py +++ b/scripts/odm_georeferencing.py @@ -7,6 +7,7 @@ from opendm import log from opendm import types from opendm import system from opendm import context +from opendm.cropper import Cropper class ODMGeoreferencingCell(ecto.Cell): @@ -42,7 +43,7 @@ class ODMGeoreferencingCell(ecto.Cell): tree = self.inputs.tree gcpfile = io.join_paths(tree.root_path, self.params.gcp_file) \ if self.params.gcp_file else find('gcp_list.txt', tree.root_path) - geocreated = True + doPointCloudGeo = True verbose = '-verbose' if self.params.verbose else '' # define paths and create working directories @@ -92,6 +93,10 @@ class ODMGeoreferencingCell(ecto.Cell): 'texturing_dir': tree.odm_texturing, 'model': os.path.join(tree.odm_texturing, tree.odm_textured_model_obj) }] + + if args.fast_orthophoto: + runs = [] + if args.use_25dmesh: runs += [{ 'georeferencing_dir': tree.odm_25dgeoreferencing, @@ -126,7 +131,10 @@ class ODMGeoreferencingCell(ecto.Cell): } if not args.use_pmvs: - kwargs['pc'] = tree.opensfm_model + if args.fast_orthophoto: + kwargs['pc'] = os.path.join(tree.opensfm, 'reconstruction.ply') + else: + kwargs['pc'] = tree.opensfm_model else: kwargs['pc'] = tree.pmvs_model @@ -153,17 +161,12 @@ class ODMGeoreferencingCell(ecto.Cell): log.ODM_WARNING('Georeferencing failed. Make sure your ' 'photos have geotags in the EXIF or you have ' 'provided a GCP file. ') - geocreated = False # skip the rest of the georeferencing + doPointCloudGeo = False # skip the rest of the georeferencing - odm_georeferencing_model_ply_geo = os.path.join(tree.odm_georeferencing, tree.odm_georeferencing_model_ply_geo) - if geocreated: - # update images metadata + if doPointCloudGeo: geo_ref = types.ODM_GeoRef() geo_ref.parse_coordinate_system(tree.odm_georeferencing_coords) - for idx, photo in enumerate(self.inputs.photos): - geo_ref.utm_to_latlon(tree.odm_georeferencing_latlon, photo, idx) - # convert ply model to LAS reference system geo_ref.convert_to_las(odm_georeferencing_model_ply_geo, tree.odm_georeferencing_model_las, @@ -186,6 +189,16 @@ class ODMGeoreferencingCell(ecto.Cell): reachedpoints = True csvfile.close() + if args.crop > 0: + log.ODM_INFO("Calculating cropping area and generating bounds shapefile from point cloud") + cropper = Cropper(tree.odm_georeferencing, 'odm_georeferenced_model') + cropper.create_bounds_shapefile(tree.odm_georeferencing_model_las, args.crop) + + # Do not execute a second time, since + # We might be doing georeferencing for + # multiple models (3D, 2.5D, ...) + doPointCloudGeo = False + else: log.ODM_WARNING('Found a valid georeferenced model in: %s' % odm_georeferencing_model_ply_geo) diff --git a/scripts/odm_meshing.py b/scripts/odm_meshing.py index 8a05552c..e46af468 100644 --- a/scripts/odm_meshing.py +++ b/scripts/odm_meshing.py @@ -1,4 +1,4 @@ -import ecto +import ecto, os from opendm import log from opendm import io @@ -19,13 +19,6 @@ class ODMeshingCell(ecto.Cell): 'is solved in the surface reconstruction step. ' 'Increasing this value increases computation ' 'times slightly but helps reduce memory usage.', 9) - - params.declare("remove_outliers", 'Percentage of outliers to remove from the point set. Set to 0 to disable. ' - 'Applies to 2.5D mesh only.', 2) - params.declare("wlop_iterations", 'Iterations of the Weighted Locally Optimal Projection (WLOP) simplification algorithm. ' - 'Higher values take longer but produce a smoother mesh. ' - 'Applies to 2.5D mesh only. ', 70) - params.declare("verbose", 'print additional messages to console', False) def declare_io(self, params, inputs, outputs): @@ -59,32 +52,37 @@ class ODMeshingCell(ecto.Cell): infile = tree.opensfm_model if args.use_pmvs: infile = tree.pmvs_model + elif args.fast_orthophoto: + infile = os.path.join(tree.opensfm, 'reconstruction.ply') - if not io.file_exists(tree.odm_mesh) or rerun_cell: - log.ODM_DEBUG('Writing ODM Mesh file in: %s' % tree.odm_mesh) + # Do not create full 3D model with fast_orthophoto + if not args.fast_orthophoto: + if not io.file_exists(tree.odm_mesh) or rerun_cell: + log.ODM_DEBUG('Writing ODM Mesh file in: %s' % tree.odm_mesh) - kwargs = { - 'bin': context.odm_modules_path, - 'outfile': tree.odm_mesh, - 'infile': infile, - 'log': tree.odm_meshing_log, - 'max_vertex': self.params.max_vertex, - 'oct_tree': self.params.oct_tree, - 'samples': self.params.samples, - 'solver': self.params.solver, - 'verbose': verbose - } + kwargs = { + 'bin': context.odm_modules_path, + 'outfile': tree.odm_mesh, + 'infile': infile, + 'log': tree.odm_meshing_log, + 'max_vertex': self.params.max_vertex, + 'oct_tree': self.params.oct_tree, + 'samples': self.params.samples, + 'solver': self.params.solver, + 'verbose': verbose + } - # run meshing binary - system.run('{bin}/odm_meshing -inputFile {infile} ' - '-outputFile {outfile} -logFile {log} ' - '-maxVertexCount {max_vertex} -octreeDepth {oct_tree} {verbose} ' - '-samplesPerNode {samples} -solverDivide {solver}'.format(**kwargs)) - else: - log.ODM_WARNING('Found a valid ODM Mesh file in: %s' % - tree.odm_mesh) + # run meshing binary + system.run('{bin}/odm_meshing -inputFile {infile} ' + '-outputFile {outfile} -logFile {log} ' + '-maxVertexCount {max_vertex} -octreeDepth {oct_tree} {verbose} ' + '-samplesPerNode {samples} -solverDivide {solver}'.format(**kwargs)) + else: + log.ODM_WARNING('Found a valid ODM Mesh file in: %s' % + tree.odm_mesh) # Do we need to generate a 2.5D mesh also? + # This is always set if fast_orthophoto is set 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) @@ -96,15 +94,15 @@ class ODMeshingCell(ecto.Cell): 'log': tree.odm_25dmeshing_log, 'verbose': verbose, 'max_vertex': self.params.max_vertex, - 'remove_outliers': self.params.remove_outliers, - 'wlop_iterations': self.params.wlop_iterations + 'neighbors': args.mesh_neighbors, + 'resolution': args.mesh_resolution } # run 2.5D meshing binary system.run('{bin}/odm_25dmeshing -inputFile {infile} ' '-outputFile {outfile} -logFile {log} ' - '-maxVertexCount {max_vertex} -outliersRemovalPercentage {remove_outliers} ' - '-wlopIterations {wlop_iterations} {verbose}'.format(**kwargs)) + '-maxVertexCount {max_vertex} -neighbors {neighbors} ' + '-resolution {resolution} {verbose}'.format(**kwargs)) else: log.ODM_WARNING('Found a valid ODM 2.5D Mesh file in: %s' % tree.odm_25dmesh) diff --git a/scripts/odm_orthophoto.py b/scripts/odm_orthophoto.py index 4174af42..b1470e96 100644 --- a/scripts/odm_orthophoto.py +++ b/scripts/odm_orthophoto.py @@ -5,6 +5,7 @@ from opendm import log from opendm import system from opendm import context from opendm import types +from opendm.cropper import Cropper class ODMOrthoPhotoCell(ecto.Cell): @@ -128,6 +129,18 @@ class ODMOrthoPhotoCell(ecto.Cell): '-a_srs \"EPSG:{epsg}\" ' '{png} {tiff} > {log}'.format(**kwargs)) + if args.crop > 0: + shapefile_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.shp') + Cropper.crop(shapefile_path, tree.odm_orthophoto_tif, { + 'TILED': 'NO' if self.params.no_tiled else 'YES', + 'COMPRESS': self.params.compress, + 'PREDICTOR': '2' if self.params.compress in ['LZW', 'DEFLATE'] else '1', + 'BIGTIFF': self.params.bigtiff, + 'BLOCKXSIZE': 512, + 'BLOCKYSIZE': 512, + 'NUM_THREADS': 'ALL_CPUS' + }) + if self.params.build_overviews: log.ODM_DEBUG("Building Overviews") kwargs = { diff --git a/scripts/opensfm.py b/scripts/opensfm.py index d61f36b5..c59d0250 100644 --- a/scripts/opensfm.py +++ b/scripts/opensfm.py @@ -52,6 +52,8 @@ class ODMOpenSfMCell(ecto.Cell): if not args.use_pmvs: output_file = tree.opensfm_model + if args.fast_orthophoto: + output_file = io.join_paths(tree.opensfm, 'reconstruction.ply') else: output_file = tree.opensfm_reconstruction @@ -141,8 +143,16 @@ class ODMOpenSfMCell(ecto.Cell): system.run('PYTHONPATH=%s %s/bin/opensfm undistort %s' % (context.pyopencv_path, context.opensfm_path, tree.opensfm)) - system.run('PYTHONPATH=%s %s/bin/opensfm compute_depthmaps %s' % + + # Skip dense reconstruction if necessary and export + # sparse reconstruction instead + if args.fast_orthophoto: + system.run('PYTHONPATH=%s %s/bin/opensfm export_ply --no-cameras %s' % (context.pyopencv_path, context.opensfm_path, tree.opensfm)) + else: + system.run('PYTHONPATH=%s %s/bin/opensfm compute_depthmaps %s' % + (context.pyopencv_path, context.opensfm_path, tree.opensfm)) + else: log.ODM_WARNING('Found a valid OpenSfM reconstruction file in: %s' % tree.opensfm_reconstruction)