From 7ec33ffb097101a92a0263ee696089d2cbc4c724 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Fri, 6 Jul 2018 16:41:56 -0400 Subject: [PATCH] Added odm_dem2points module, better skirts algo --- modules/CMakeLists.txt | 1 + modules/odm_cleanmesh/src/main.cpp | 2 +- modules/odm_dem2points/CMakeLists.txt | 19 ++ modules/odm_dem2points/src/CmdLineParser.h | 106 +++++++ modules/odm_dem2points/src/CmdLineParser.inl | 300 +++++++++++++++++++ modules/odm_dem2points/src/Logger.h | 33 ++ modules/odm_dem2points/src/main.cpp | 234 +++++++++++++++ opendm/mesh.py | 134 +++++---- 8 files changed, 770 insertions(+), 59 deletions(-) create mode 100644 modules/odm_dem2points/CMakeLists.txt create mode 100644 modules/odm_dem2points/src/CmdLineParser.h create mode 100644 modules/odm_dem2points/src/CmdLineParser.inl create mode 100644 modules/odm_dem2points/src/Logger.h create mode 100644 modules/odm_dem2points/src/main.cpp diff --git a/modules/CMakeLists.txt b/modules/CMakeLists.txt index 7f26197a..637a6711 100644 --- a/modules/CMakeLists.txt +++ b/modules/CMakeLists.txt @@ -8,6 +8,7 @@ add_subdirectory(odm_extract_utm) add_subdirectory(odm_georef) add_subdirectory(odm_orthophoto) add_subdirectory(odm_cleanmesh) +add_subdirectory(odm_dem2points) if (ODM_BUILD_SLAM) add_subdirectory(odm_slam) endif () diff --git a/modules/odm_cleanmesh/src/main.cpp b/modules/odm_cleanmesh/src/main.cpp index 81a9d429..c5251c0f 100644 --- a/modules/odm_cleanmesh/src/main.cpp +++ b/modules/odm_cleanmesh/src/main.cpp @@ -58,7 +58,7 @@ int main(int argc, char **argv) { logWriter.verbose = Verbose.set; - logWriter.outputFile = "odm_cleanmesh_log.txt"; + // logWriter.outputFile = "odm_cleanmesh_log.txt"; logArgs(params, logWriter); vtkSmartPointer reader = diff --git a/modules/odm_dem2points/CMakeLists.txt b/modules/odm_dem2points/CMakeLists.txt new file mode 100644 index 00000000..ff727a2e --- /dev/null +++ b/modules/odm_dem2points/CMakeLists.txt @@ -0,0 +1,19 @@ +project(odm_dem2points) +cmake_minimum_required(VERSION 2.8) + +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}) + +set (CMAKE_CXX_STANDARD 11) +find_package(GDAL REQUIRED) +include_directories(${GDAL_INCLUDE_DIR}) + +# Add compiler options. +add_definitions(-Wall -Wextra) + +# Add source directory +aux_source_directory("./src" SRC_LIST) + +# Add exectuteable +add_executable(${PROJECT_NAME} ${SRC_LIST}) + +target_link_libraries(${PROJECT_NAME} ${GDAL_LIBRARY}) diff --git a/modules/odm_dem2points/src/CmdLineParser.h b/modules/odm_dem2points/src/CmdLineParser.h new file mode 100644 index 00000000..a1fbd0e3 --- /dev/null +++ b/modules/odm_dem2points/src/CmdLineParser.h @@ -0,0 +1,106 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. Redistributions in binary form must reproduce +the above copyright notice, this list of conditions and the following disclaimer +in the documentation and/or other materials provided with the distribution. + +Neither the name of the Johns Hopkins University nor the names of its contributors +may be used to endorse or promote products derived from this software without specific +prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#ifndef CMD_LINE_PARSER_INCLUDED +#define CMD_LINE_PARSER_INCLUDED + +#include +#include +#include +#include +#include + +#ifdef WIN32 +int strcasecmp( const char* c1 , const char* c2 ); +#endif // WIN32 + +class cmdLineReadable +{ +public: + bool set; + char *name; + cmdLineReadable( const char *name ); + virtual ~cmdLineReadable( void ); + virtual int read( char** argv , int argc ); + virtual void writeValue( char* str ) const; +}; + +template< class Type > void cmdLineWriteValue( Type t , char* str ); +template< class Type > void cmdLineCleanUp( Type* t ); +template< class Type > Type cmdLineInitialize( void ); +template< class Type > Type cmdLineCopy( Type t ); +template< class Type > Type cmdLineStringToType( const char* str ); + +template< class Type > +class cmdLineParameter : public cmdLineReadable +{ +public: + Type value; + cmdLineParameter( const char *name ); + cmdLineParameter( const char *name , Type v ); + ~cmdLineParameter( void ); + int read( char** argv , int argc ); + void writeValue( char* str ) const; + bool expectsArg( void ) const { return true; } +}; + +template< class Type , int Dim > +class cmdLineParameterArray : public cmdLineReadable +{ +public: + Type values[Dim]; + cmdLineParameterArray( const char *name, const Type* v=NULL ); + ~cmdLineParameterArray( void ); + int read( char** argv , int argc ); + void writeValue( char* str ) const; + bool expectsArg( void ) const { return true; } +}; + +template< class Type > +class cmdLineParameters : public cmdLineReadable +{ +public: + int count; + Type *values; + cmdLineParameters( const char* name ); + ~cmdLineParameters( void ); + int read( char** argv , int argc ); + void writeValue( char* str ) const; + bool expectsArg( void ) const { return true; } +}; + +void cmdLineParse( int argc , char **argv, cmdLineReadable** params ); +char* FileExtension( char* fileName ); +char* LocalFileName( char* fileName ); +char* DirectoryName( char* fileName ); +char* GetFileExtension( const char* fileName ); +char* GetLocalFileName( const char* fileName ); +char** ReadWords( const char* fileName , int& cnt ); + +#include "CmdLineParser.inl" +#endif // CMD_LINE_PARSER_INCLUDED diff --git a/modules/odm_dem2points/src/CmdLineParser.inl b/modules/odm_dem2points/src/CmdLineParser.inl new file mode 100644 index 00000000..84574a57 --- /dev/null +++ b/modules/odm_dem2points/src/CmdLineParser.inl @@ -0,0 +1,300 @@ +/* -*- C++ -*- +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. Redistributions in binary form must reproduce +the above copyright notice, this list of conditions and the following disclaimer +in the documentation and/or other materials provided with the distribution. + +Neither the name of the Johns Hopkins University nor the names of its contributors +may be used to endorse or promote products derived from this software without specific +prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#include +#include + +#if defined( WIN32 ) || defined( _WIN64 ) +inline int strcasecmp( const char* c1 , const char* c2 ){ return _stricmp( c1 , c2 ); } +#endif // WIN32 || _WIN64 + +template< > void cmdLineCleanUp< int >( int* t ){ *t = 0; } +template< > void cmdLineCleanUp< float >( float* t ){ *t = 0; } +template< > void cmdLineCleanUp< double >( double* t ){ *t = 0; } +template< > void cmdLineCleanUp< char* >( char** t ){ if( *t ) free( *t ) ; *t = NULL; } +template< > int cmdLineInitialize< int >( void ){ return 0; } +template< > float cmdLineInitialize< float >( void ){ return 0.f; } +template< > double cmdLineInitialize< double >( void ){ return 0.; } +template< > char* cmdLineInitialize< char* >( void ){ return NULL; } +template< > void cmdLineWriteValue< int >( int t , char* str ){ sprintf( str , "%d" , t ); } +template< > void cmdLineWriteValue< float >( float t , char* str ){ sprintf( str , "%f" , t ); } +template< > void cmdLineWriteValue< double >( double t , char* str ){ sprintf( str , "%f" , t ); } +template< > void cmdLineWriteValue< char* >( char* t , char* str ){ if( t ) sprintf( str , "%s" , t ) ; else str[0]=0; } +template< > int cmdLineCopy( int t ){ return t; } +template< > float cmdLineCopy( float t ){ return t; } +template< > double cmdLineCopy( double t ){ return t; } +#if defined( WIN32 ) || defined( _WIN64 ) +template< > char* cmdLineCopy( char* t ){ return _strdup( t ); } +#else // !WIN32 && !_WIN64 +template< > char* cmdLineCopy( char* t ){ return strdup( t ); } +#endif // WIN32 || _WIN64 +template< > int cmdLineStringToType( const char* str ){ return atoi( str ); } +template< > float cmdLineStringToType( const char* str ){ return float( atof( str ) ); } +template< > double cmdLineStringToType( const char* str ){ return double( atof( str ) ); } +#if defined( WIN32 ) || defined( _WIN64 ) +template< > char* cmdLineStringToType( const char* str ){ return _strdup( str ); } +#else // !WIN32 && !_WIN64 +template< > char* cmdLineStringToType( const char* str ){ return strdup( str ); } +#endif // WIN32 || _WIN64 + + +///////////////////// +// cmdLineReadable // +///////////////////// +#if defined( WIN32 ) || defined( _WIN64 ) +inline cmdLineReadable::cmdLineReadable( const char *name ) : set(false) { this->name = _strdup( name ); } +#else // !WIN32 && !_WIN64 +inline cmdLineReadable::cmdLineReadable( const char *name ) : set(false) { this->name = strdup( name ); } +#endif // WIN32 || _WIN64 + +inline cmdLineReadable::~cmdLineReadable( void ){ if( name ) free( name ) ; name = NULL; } +inline int cmdLineReadable::read( char** , int ){ set = true ; return 0; } +inline void cmdLineReadable::writeValue( char* str ) const { str[0] = 0; } + +////////////////////// +// cmdLineParameter // +////////////////////// +template< class Type > cmdLineParameter< Type >::~cmdLineParameter( void ) { cmdLineCleanUp( &value ); } +template< class Type > cmdLineParameter< Type >::cmdLineParameter( const char *name ) : cmdLineReadable( name ){ value = cmdLineInitialize< Type >(); } +template< class Type > cmdLineParameter< Type >::cmdLineParameter( const char *name , Type v ) : cmdLineReadable( name ){ value = cmdLineCopy< Type >( v ); } +template< class Type > +int cmdLineParameter< Type >::read( char** argv , int argc ) +{ + if( argc>0 ) + { + cmdLineCleanUp< Type >( &value ) , value = cmdLineStringToType< Type >( argv[0] ); + set = true; + return 1; + } + else return 0; +} +template< class Type > +void cmdLineParameter< Type >::writeValue( char* str ) const { cmdLineWriteValue< Type >( value , str ); } + + +/////////////////////////// +// cmdLineParameterArray // +/////////////////////////// +template< class Type , int Dim > +cmdLineParameterArray< Type , Dim >::cmdLineParameterArray( const char *name , const Type* v ) : cmdLineReadable( name ) +{ + if( v ) for( int i=0 ; i( v[i] ); + else for( int i=0 ; i(); +} +template< class Type , int Dim > +cmdLineParameterArray< Type , Dim >::~cmdLineParameterArray( void ){ for( int i=0 ; i( values+i ); } +template< class Type , int Dim > +int cmdLineParameterArray< Type , Dim >::read( char** argv , int argc ) +{ + if( argc>=Dim ) + { + for( int i=0 ; i( values+i ) , values[i] = cmdLineStringToType< Type >( argv[i] ); + set = true; + return Dim; + } + else return 0; +} +template< class Type , int Dim > +void cmdLineParameterArray< Type , Dim >::writeValue( char* str ) const +{ + char* temp=str; + for( int i=0 ; i( values[i] , temp ); + temp = str+strlen( str ); + } +} +/////////////////////// +// cmdLineParameters // +/////////////////////// +template< class Type > +cmdLineParameters< Type >::cmdLineParameters( const char* name ) : cmdLineReadable( name ) , values(NULL) , count(0) { } +template< class Type > +cmdLineParameters< Type >::~cmdLineParameters( void ) +{ + if( values ) delete[] values; + values = NULL; + count = 0; +} +template< class Type > +int cmdLineParameters< Type >::read( char** argv , int argc ) +{ + if( values ) delete[] values; + values = NULL; + + if( argc>0 ) + { + count = atoi(argv[0]); + if( count <= 0 || argc <= count ) return 1; + values = new Type[count]; + if( !values ) return 0; + for( int i=0 ; i( argv[i+1] ); + set = true; + return count+1; + } + else return 0; +} +template< class Type > +void cmdLineParameters< Type >::writeValue( char* str ) const +{ + char* temp=str; + for( int i=0 ; i( values[i] , temp ); + temp = str+strlen( str ); + } +} + + +inline char* FileExtension( char* fileName ) +{ + char* temp = fileName; + for( unsigned int i=0 ; i=0 ; i-- ) + if( fileName[i] =='\\' ) + { + fileName[i] = 0; + return fileName; + } + fileName[0] = 0; + return fileName; +} + +inline void cmdLineParse( int argc , char **argv , cmdLineReadable** params ) +{ + while( argc>0 ) + { + if( argv[0][0]=='-' ) + { + cmdLineReadable* readable=NULL; + for( int i=0 ; params[i]!=NULL && readable==NULL ; i++ ) if( !strcasecmp( params[i]->name , argv[0]+1 ) ) readable = params[i]; + if( readable ) + { + int j = readable->read( argv+1 , argc-1 ); + argv += j , argc -= j; + } + else + { + fprintf( stderr , "[WARNING] Invalid option: %s\n" , argv[0] ); + for( int i=0 ; params[i]!=NULL ; i++ ) printf( "\t-%s\n" , params[i]->name ); + } + } + else fprintf( stderr , "[WARNING] Parameter name should be of the form -: %s\n" , argv[0] ); + ++argv , --argc; + } +} + +inline char** ReadWords(const char* fileName,int& cnt) +{ + char** names; + char temp[500]; + FILE* fp; + + fp=fopen(fileName,"r"); + if(!fp){return NULL;} + cnt=0; + while(fscanf(fp," %s ",temp)==1){cnt++;} + fclose(fp); + + names=new char*[cnt]; + if(!names){return NULL;} + + fp=fopen(fileName,"r"); + if(!fp){ + delete[] names; + cnt=0; + return NULL; + } + cnt=0; + while(fscanf(fp," %s ",temp)==1){ + names[cnt]=new char[strlen(temp)+1]; + if(!names){ + for(int j=0;j +#include +#include "CmdLineParser.h" + +struct Logger{ + bool verbose; + const char* outputFile; + + Logger(){ + this->verbose = false; + this->outputFile = NULL; + } + + void operator() ( const char* format , ... ) + { + if( outputFile ) + { + FILE* fp = fopen( outputFile , "a" ); + va_list args; + va_start( args , format ); + vfprintf( fp , format , args ); + fclose( fp ); + va_end( args ); + } + if( verbose ) + { + va_list args; + va_start( args , format ); + vprintf( format , args ); + va_end( args ); + } + } +}; diff --git a/modules/odm_dem2points/src/main.cpp b/modules/odm_dem2points/src/main.cpp new file mode 100644 index 00000000..227a59d8 --- /dev/null +++ b/modules/odm_dem2points/src/main.cpp @@ -0,0 +1,234 @@ +#include +#include +#include +#include "CmdLineParser.h" +#include "Logger.h" + +#include "gdal_priv.h" +#include "cpl_conv.h" // for CPLMalloc() + +Logger logWriter; + +typedef struct Point2D{ + double x; + double y; + + Point2D(): x(0.0), y(0.0){} + Point2D(double x, double y): x(x), y(y){} +} Point2D; + +typedef struct BoundingBox{ + Point2D max; + Point2D min; + + BoundingBox(): max(Point2D()), min(Point2D()){} + BoundingBox(Point2D min, Point2D max): max(max), min(min){} +} BoundingBox; + +struct PlyPoint{ + float x; + float y; + float z; + float nx; + float ny; + float nz; +} p; +size_t psize = sizeof(float) * 6; + +float *rasterData; +int arr_width, arr_height; + +#define IS_BIG_ENDIAN (*(uint16_t *)"\0\xff" < 0x100) + +cmdLineParameter< char* > + InputFile( "inputFile" ) , + OutputFile( "outputFile" ); +cmdLineParameter< float > + SkirtHeightThreshold( "skirtHeightThreshold" ), + SkirtIncrements("skirtIncrements"), + SkirtHeightCap( "skirtHeightCap"); +cmdLineReadable + Verbose( "verbose" ); + +cmdLineReadable* params[] = { + &InputFile , &OutputFile , &SkirtHeightThreshold, &SkirtIncrements, &SkirtHeightCap, &Verbose , + NULL +}; + +void help(char *ex){ + std::cout << "Usage: " << ex << std::endl + << "\t -" << InputFile.name << " " << std::endl + << "\t -" << OutputFile.name << " " << std::endl + << "\t [-" << SkirtHeightThreshold.name << " ]" << std::endl + << "\t [-" << SkirtIncrements.name << " ]" << std::endl + << "\t [-" << SkirtHeightCap.name << " ]" << std::endl + << "\t [-" << Verbose.name << "]" << std::endl; + exit(EXIT_FAILURE); +} + + +void logArgs(cmdLineReadable* params[], Logger& logWriter){ + logWriter("Running with parameters:\n"); + char str[1024]; + for( int i=0 ; params[i] ; i++ ){ + if( params[i]->set ){ + params[i]->writeValue( str ); + if( strlen( str ) ) logWriter( "\t--%s %s\n" , params[i]->name , str ); + else logWriter( "\t--%s\n" , params[i]->name ); + } + } +} + +Point2D geoLoc(float xloc, float yloc, double *affine){ + return Point2D(affine[0] + xloc*affine[1] + yloc*affine[2], affine[3] + xloc*affine[4] + yloc*affine[5]); +} + +BoundingBox getExtent(GDALDataset *dataset){ + double affine[6]; + dataset->GetGeoTransform(affine); + return BoundingBox(geoLoc(0, dataset->GetRasterYSize(), affine), geoLoc(dataset->GetRasterXSize(), 0, affine)); +} + +int countSkirts(int nx, int ny){ + float neighbor_z = rasterData[ny * arr_width + nx]; + float current_z = p.z; + int result = 0; + + if (current_z - neighbor_z > SkirtHeightThreshold.value && current_z - neighbor_z < SkirtHeightCap.value){ + while (current_z > neighbor_z){ + current_z -= SkirtIncrements.value; + result++; + } + result++; + } + + return result; +} + +void writeSkirts(std::ofstream &f, int nx, int ny){ + float neighbor_z = rasterData[ny * arr_width + nx]; + float z_val = p.z; + + if (p.z - neighbor_z > SkirtHeightThreshold.value && p.z - neighbor_z < SkirtHeightCap.value){ + while (p.z > neighbor_z){ + p.z -= SkirtIncrements.value; + f.write(reinterpret_cast(&p), psize); + } + p.z = neighbor_z; + f.write(reinterpret_cast(&p), psize); + } + + // Restore original p.z value + p.z = z_val; +} + +int main(int argc, char **argv) { + cmdLineParse( argc-1 , &argv[1] , params ); + if( !InputFile.set || !OutputFile.set ) help(argv[0]); + + if (!SkirtHeightThreshold.set) SkirtHeightThreshold.value = 1.5f; + if (!SkirtIncrements.set) SkirtIncrements.value = 0.1f; + if (!SkirtHeightCap.set) SkirtHeightCap.value = 100.0f; + + logWriter.verbose = Verbose.set; + // logWriter.outputFile = "odm_dem2points_log.txt"; + logArgs(params, logWriter); + + + GDALDataset *dataset; + GDALAllRegister(); + dataset = (GDALDataset *) GDALOpen( InputFile.value, GA_ReadOnly ); + if( dataset != NULL ) + { + arr_width = dataset->GetRasterXSize(); + arr_height = dataset->GetRasterYSize(); + + logWriter("Raster Size is %dx%d\n", arr_width, arr_height); + BoundingBox extent = getExtent(dataset); + + logWriter("Extent is (%f, %f), (%f, %f)\n", extent.min.x, extent.max.x, extent.min.y, extent.max.y); + + float ext_width = extent.max.x - extent.min.x; + float ext_height = extent.max.y - extent.min.y; + + int vertex_count = (arr_height - 4) * (arr_width - 4); + int skirt_points = 0; + + GDALRasterBand *band = dataset->GetRasterBand(1); + + rasterData = new float[arr_width * arr_height]; + + if (band->RasterIO( GF_Read, 0, 0, arr_width, arr_height, + rasterData, arr_width, arr_height, GDT_Float32, 0, 0 ) == CE_Failure){ + std::cerr << "Cannot access raster data" << std::endl; + exit(1); + } + + // On the first pass we calculate the number of skirts + // on the second pass we sample and write the points + + logWriter("Calculating skirts... "); + + for (int y = 2; y < arr_height - 2; y++){ + for (int x = 2; x < arr_width - 2; x++){ + p.z = rasterData[y * arr_width + x]; + + skirt_points += countSkirts(x, y + 1); + skirt_points += countSkirts(x, y - 1); + skirt_points += countSkirts(x + 1, y); + skirt_points += countSkirts(x - 1, y); + } + } + + logWriter("%d vertices will be added\n", skirt_points); + logWriter("Total vertices: %d\n", (skirt_points + vertex_count)); + logWriter("Sampling and writing to file..."); + + // Starting writing ply file + std::ofstream f (OutputFile.value); + f << "ply" << std::endl; + + if (IS_BIG_ENDIAN){ + f << "format binary_big_endian 1.0" << std::endl; + }else{ + f << "format binary_little_endian 1.0" << std::endl; + } + + f << "element vertex " << (vertex_count + skirt_points) << std::endl + << "property float x" << std::endl + << "property float y" << std::endl + << "property float z" << std::endl + << "property float nx" << std::endl + << "property float ny" << std::endl + << "property float nz" << std::endl + << "end_header" << std::endl; + + // Normals are always pointing up + p.nx = 0; + p.ny = 0; + p.nz = 1; + + for (int y = 2; y < arr_height - 2; y++){ + for (int x = 2; x < arr_width - 2; x++){ + p.z = rasterData[y * arr_width + x]; + p.x = extent.min.x + (static_cast(x) / static_cast(arr_width)) * ext_width; + p.y = extent.max.y - (static_cast(y) / static_cast(arr_height)) * ext_height; + + f.write(reinterpret_cast(&p), psize); + + writeSkirts(f, x, y + 1); + writeSkirts(f, x, y - 1); + writeSkirts(f, x + 1, y); + writeSkirts(f, x - 1, y); + } + } + + logWriter(" done!\n"); + + f.close(); + GDALClose(dataset); + + }else{ + std::cerr << "Cannot open " << InputFile.value << std::endl; + } +} diff --git a/opendm/mesh.py b/opendm/mesh.py index 4c9e08b6..c0f6ddbd 100644 --- a/opendm/mesh.py +++ b/opendm/mesh.py @@ -1,6 +1,5 @@ from __future__ import absolute_import import os, shutil, sys, struct, random -from io import BytesIO from gippy import GeoImage from opendm.dem import commands from opendm import system @@ -36,7 +35,7 @@ def create_25dmesh(inPointCloud, outMesh, dsm_resolution=0.05, depth=8, samples= max_workers=max_workers ) - dsm_points = dem_to_points(os.path.join(tmp_directory, 'mesh_dsm.tif'), os.path.join(tmp_directory, 'dsm_points.ply')) + dsm_points = dem_to_points(os.path.join(tmp_directory, 'mesh_dsm.tif'), os.path.join(tmp_directory, 'dsm_points.ply'), verbose) mesh = screened_poisson_reconstruction(dsm_points, outMesh, depth=depth, samples=samples, maxVertexCount=maxVertexCount, @@ -49,68 +48,87 @@ def create_25dmesh(inPointCloud, outMesh, dsm_resolution=0.05, depth=8, samples= return mesh -def dem_to_points(inGeotiff, outPointCloud): +def dem_to_points(inGeotiff, outPointCloud, verbose=False): log.ODM_INFO('Sampling points from DSM: %s' % inGeotiff) - image = GeoImage.open([inGeotiff], bandnames=['z'], nodata=-9999) - arr = image['z'].read_raw() - resolution = max(abs(image.resolution().x()), abs(image.resolution().y())) + kwargs = { + 'bin': context.odm_modules_path, + 'outfile': outPointCloud, + 'infile': inGeotiff, + 'verbose': '-verbose' if verbose else '' + } - log.ODM_INFO('Writing points...') - - mem_file = BytesIO() - - xmin, xmax, ymin, ymax = image.extent().x0(), image.extent().x1(), image.extent().y0(), image.extent().y1() - ext_width, ext_height = xmax - xmin, ymax - ymin - arr_height, arr_width = arr.shape - vertex_count = (arr_height - 4) * (arr_width - 4) - skirt_points = 0 - - skirt_height_threshold = 1 # meter - skirt_increments = 0.1 - - for x in range(2, arr_width - 2): - for y in range(2, arr_height - 2): - z = arr[y][x] - tx = xmin + (float(x) / float(arr_width)) * ext_width - ty = ymax - (float(y) / float(arr_height)) * ext_height - mem_file.write(struct.pack('ffffff', tx, ty, z, 0, 0, 1)) - - # Skirting - for (nx, ny) in ((x, y + 1), (x, y - 1), (x + 1, y), (x - 1, y)): - current_z = z - neighbor_z = arr[ny][nx] - - if current_z - neighbor_z > skirt_height_threshold: - while current_z > neighbor_z: - current_z -= skirt_increments - mem_file.write(struct.pack('ffffff', tx, ty, current_z, 0, 0, 1)) - skirt_points += 1 - - mem_file.write(struct.pack('ffffff', tx, ty, neighbor_z, 0, 0, 1)) - skirt_points += 1 - - - with open(outPointCloud, "wb") as f: - f.write("ply\n") - f.write("format binary_%s_endian 1.0\n" % sys.byteorder) - f.write("element vertex %s\n" % (vertex_count + skirt_points)) - f.write("property float x\n") - f.write("property float y\n") - f.write("property float z\n") - f.write("property float nx\n") - f.write("property float ny\n") - f.write("property float nz\n") - f.write("end_header\n") - f.write(mem_file.getvalue()) - - mem_file.close() - - log.ODM_INFO("Points count: %s (%s samples, %s skirts)", vertex_count + skirt_points, vertex_count, skirt_points) - log.ODM_INFO('Wrote points to: %s' % outPointCloud) + system.run('{bin}/odm_dem2points -inputFile {infile} ' + '-outputFile {outfile} ' + '-skirtHeightThreshold 1.5 ' + '-skirtIncrements 0.2 ' + '-skirtHeightCap 100 ' + ' {verbose} '.format(**kwargs)) return outPointCloud +# Old Python implementation of dem_to_points +# def dem_to_points(inGeotiff, outPointCloud): +# log.ODM_INFO('Sampling points from DSM: %s' % inGeotiff) + +# image = GeoImage.open([inGeotiff], bandnames=['z'], nodata=-9999) +# arr = image['z'].read_raw() + +# mem_file = BytesIO() + +# xmin, xmax, ymin, ymax = image.extent().x0(), image.extent().x1(), image.extent().y0(), image.extent().y1() +# ext_width, ext_height = xmax - xmin, ymax - ymin +# arr_height, arr_width = arr.shape +# vertex_count = (arr_height - 4) * (arr_width - 4) +# skirt_points = 0 + +# skirt_height_threshold = 1 # meter +# skirt_increments = 0.1 + +# for x in range(2, arr_width - 2): +# for y in range(2, arr_height - 2): +# z = arr[y][x] +# tx = xmin + (float(x) / float(arr_width)) * ext_width +# ty = ymax - (float(y) / float(arr_height)) * ext_height +# mem_file.write(struct.pack('ffffff', tx, ty, z, 0, 0, 1)) + +# # Skirting +# for (nx, ny) in ((x, y + 1), (x, y - 1), (x + 1, y), (x - 1, y)): +# current_z = z +# neighbor_z = arr[ny][nx] + +# if current_z - neighbor_z > skirt_height_threshold: +# while current_z > neighbor_z: +# current_z -= skirt_increments +# mem_file.write(struct.pack('ffffff', tx, ty, current_z, 0, 0, 1)) +# skirt_points += 1 + +# mem_file.write(struct.pack('ffffff', tx, ty, neighbor_z, 0, 0, 1)) +# skirt_points += 1 + +# log.ODM_INFO("Points count: %s (%s samples, %s skirts)", vertex_count + skirt_points, vertex_count, skirt_points) +# log.ODM_INFO('Writing points...') + +# mem_file.seek(0) +# with open(outPointCloud, "wb") as f: +# f.write("ply\n") +# f.write("format binary_%s_endian 1.0\n" % sys.byteorder) +# f.write("element vertex %s\n" % (vertex_count + skirt_points)) +# f.write("property float x\n") +# f.write("property float y\n") +# f.write("property float z\n") +# f.write("property float nx\n") +# f.write("property float ny\n") +# f.write("property float nz\n") +# f.write("end_header\n") +# shutil.copyfileobj(mem_file, f) + +# mem_file.close() + +# log.ODM_INFO('Wrote points to: %s' % outPointCloud) + +# return outPointCloud + def screened_poisson_reconstruction(inPointCloud, outMesh, depth = 8, samples = 1, maxVertexCount=100000, pointWeight=4, threads=context.num_cores, verbose=False):