From b7fd10163e3597e5ebc2d4283b97883016613ccd Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Mon, 22 Oct 2018 17:27:40 -0400 Subject: [PATCH] First commit --- .gitignore | 2 + CMakeLists.txt | 17 + README.md | 34 ++ src/CmdLineParser.h | 106 ++++++ src/CmdLineParser.inl | 300 +++++++++++++++++ src/Logger.h | 33 ++ src/Simplify.h | 749 ++++++++++++++++++++++++++++++++++++++++++ src/main.cpp | 496 ++++++++++++++++++++++++++++ 8 files changed, 1737 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 src/CmdLineParser.h create mode 100644 src/CmdLineParser.inl create mode 100644 src/Logger.h create mode 100644 src/Simplify.h create mode 100644 src/main.cpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..dc84959 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +build/ + diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..f754184 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,17 @@ +project(dem2mesh) +cmake_minimum_required(VERSION 2.8) + +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/README.md b/README.md index 157f944..75980d3 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,36 @@ # dem2mesh + Fast generation of 2.5D meshes from elevation models. + +## Dependencies + +GDAL is the only requirement. To install it run: + +``` +sudo apt-get install -y libgdal-dev +``` + +## Building + +```bash +mkdir build && cd build +cmake .. +make +``` + +## Running + +```bash +./dem2mesh -rtc -verbose -inputFile dem.tif -outputFile mesh.ply +``` + +## Options + +| Flag | Description | Required | +| --- | --- | --- | +| -inputFile | Path to DEM raster datasource (GDAL compatible) | ✔ | +| -outputFile | Path to PLY mesh output | ✔ | +| -maxVertexCount | Target number of vertices for the output mesh. This number is not always guaranteed to be reached. Defaults to `100000` | | +| -maxTileLength | Max length of a tile. Smaller values take longer to process but reduce memory usage by splitting the meshing process into tiles. Defaults to `1000`. | | +| -rtc | Use Relative To Center (RTC) X/Y coordinates in the output PLY mesh. This can be useful since most 3D visualization software use floating coordinate precision to represent vertices and using absolute coordinates might lead to jittering artifacts. | | +| -verbose | Print verbose output. | | diff --git a/src/CmdLineParser.h b/src/CmdLineParser.h new file mode 100644 index 0000000..a1fbd0e --- /dev/null +++ b/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/src/CmdLineParser.inl b/src/CmdLineParser.inl new file mode 100644 index 0000000..84574a5 --- /dev/null +++ b/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 ); + vfprintf( stdout, format , args ); + va_end( args ); + } + } +}; diff --git a/src/Simplify.h b/src/Simplify.h new file mode 100644 index 0000000..406cac4 --- /dev/null +++ b/src/Simplify.h @@ -0,0 +1,749 @@ +///////////////////////////////////////////// +// +// Mesh Simplification Tutorial +// +// (C) by Sven Forstmann in 2014 +// +// License : MIT +// http://opensource.org/licenses/MIT +// +//https://github.com/sp4cerat/Fast-Quadric-Mesh-Simplification +// +// 5/2016: Chris Rorden created minimal version for OSX/Linux/Windows compile + +//#include +//#include +//#include +//#include +//#include +#include +//#include +//#include +#include +#include +#include +#include +#include +#include +#include //FLT_EPSILON, DBL_EPSILON + +#define loopi(start_l,end_l) for ( int i=start_l;i1) input=1; + return (double) acos ( input ); + } + + inline double angle2( const vec3f& v , const vec3f& w ) + { + vec3f a = v , b= *this; + double dot = a.x*b.x + a.y*b.y + a.z*b.z; + double len = a.length() * b.length(); + if(len==0)len=1; + + vec3f plane; plane.cross( b,w ); + + if ( plane.x * a.x + plane.y * a.y + plane.z * a.z > 0 ) + return (double) -acos ( dot / len ); + + return (double) acos ( dot / len ); + } + + inline vec3f rot_x( double a ) + { + double yy = cos ( a ) * y + sin ( a ) * z; + double zz = cos ( a ) * z - sin ( a ) * y; + y = yy; z = zz; + return *this; + } + inline vec3f rot_y( double a ) + { + double xx = cos ( -a ) * x + sin ( -a ) * z; + double zz = cos ( -a ) * z - sin ( -a ) * x; + x = xx; z = zz; + return *this; + } + inline void clamp( double min, double max ) + { + if (xmax) x=max; + if (y>max) y=max; + if (z>max) z=max; + } + inline vec3f rot_z( double a ) + { + double yy = cos ( a ) * y + sin ( a ) * x; + double xx = cos ( a ) * x - sin ( a ) * y; + y = yy; x = xx; + return *this; + } + inline vec3f invert() + { + x=-x;y=-y;z=-z;return *this; + } + inline vec3f frac() + { + return vec3f( + x-double(int(x)), + y-double(int(y)), + z-double(int(z)) + ); + } + + inline vec3f integer() + { + return vec3f( + double(int(x)), + double(int(y)), + double(int(z)) + ); + } + + inline double length() const + { + return (double)sqrt(x*x + y*y + z*z); + } + + inline vec3f normalize( double desired_length = 1 ) + { + double square = sqrt(x*x + y*y + z*z); + /* + if (square <= 0.00001f ) + { + x=1;y=0;z=0; + return *this; + }*/ + //double len = desired_length / square; + x/=square;y/=square;z/=square; + + return *this; + } + static vec3f normalize( vec3f a ); + + static void random_init(); + static double random_double(); + static vec3f random(); + + static int random_number; + + double random_double_01(double a){ + double rnf=a*14.434252+a*364.2343+a*4213.45352+a*2341.43255+a*254341.43535+a*223454341.3523534245+23453.423412; + int rni=((int)rnf)%100000; + return double(rni)/(100000.0f-1.0f); + } + + vec3f random01_fxyz(){ + x=(double)random_double_01(x); + y=(double)random_double_01(y); + z=(double)random_double_01(z); + return *this; + } + +}; + +double min(double v1, double v2) { + return fmin(v1,v2); +} + + +class SymetricMatrix { + + public: + + // Constructor + + SymetricMatrix(double c=0) { loopi(0,10) m[i] = c; } + + SymetricMatrix( double m11, double m12, double m13, double m14, + double m22, double m23, double m24, + double m33, double m34, + double m44) { + m[0] = m11; m[1] = m12; m[2] = m13; m[3] = m14; + m[4] = m22; m[5] = m23; m[6] = m24; + m[7] = m33; m[8] = m34; + m[9] = m44; + } + + // Make plane + + SymetricMatrix(double a,double b,double c,double d) + { + m[0] = a*a; m[1] = a*b; m[2] = a*c; m[3] = a*d; + m[4] = b*b; m[5] = b*c; m[6] = b*d; + m[7 ] =c*c; m[8 ] = c*d; + m[9 ] = d*d; + } + + double operator[](int c) const { return m[c]; } + + // Determinant + + double det( int a11, int a12, int a13, + int a21, int a22, int a23, + int a31, int a32, int a33) + { + double det = m[a11]*m[a22]*m[a33] + m[a13]*m[a21]*m[a32] + m[a12]*m[a23]*m[a31] + - m[a13]*m[a22]*m[a31] - m[a11]*m[a23]*m[a32]- m[a12]*m[a21]*m[a33]; + return det; + } + + const SymetricMatrix operator+(const SymetricMatrix& n) const + { + return SymetricMatrix( m[0]+n[0], m[1]+n[1], m[2]+n[2], m[3]+n[3], + m[4]+n[4], m[5]+n[5], m[6]+n[6], + m[ 7]+n[ 7], m[ 8]+n[8 ], + m[ 9]+n[9 ]); + } + + SymetricMatrix& operator+=(const SymetricMatrix& n) + { + m[0]+=n[0]; m[1]+=n[1]; m[2]+=n[2]; m[3]+=n[3]; + m[4]+=n[4]; m[5]+=n[5]; m[6]+=n[6]; m[7]+=n[7]; + m[8]+=n[8]; m[9]+=n[9]; + return *this; + } + + double m[10]; +}; +/////////////////////////////////////////// + +namespace Simplify +{ + // Global Variables & Strctures + struct Triangle { int v[3];double err[4];int8_t deleted,dirty;vec3f n; }; + struct Vertex { vec3f p;int tstart,tcount;SymetricMatrix q;int8_t border;}; + struct Ref { int tid,tvertex; }; + std::vector triangles; + std::vector vertices; + std::vector refs; + + // Helper functions + + double vertex_error(SymetricMatrix q, double x, double y, double z); + double calculate_error(int id_v1, int id_v2, vec3f &p_result); + bool flipped(vec3f p,int i0,int i1,Vertex &v0,Vertex &v1,std::vector &deleted); + void update_triangles(int i0,Vertex &v,std::vector &deleted,int &deleted_triangles); + void update_mesh(int iteration); + void compact_mesh(); + // + // Main simplification function + // + // target_count : target nr. of triangles + // agressiveness : sharpness to increase the threashold. + // 5..8 are good numbers + // more iterations yield higher quality + // + + void simplify_mesh(int target_count, double agressiveness=7, bool verbose=false) + { + // init +// loopi(0,triangles.size()) +// { +// triangles[i].deleted=0; +// } + + // main iteration loop + int deleted_triangles=0; + std::vector deleted0,deleted1; + int triangle_count=triangles.size(); + //int iteration = 0; + //loop(iteration,0,100) + for (int iteration = 0; iteration < 100; iteration ++) + { + if(triangle_count-deleted_triangles<=target_count)break; + // + // All triangles with edges below the threshold will be removed + // + // The following numbers works well for most models. + // If it does not, try to adjust the 3 parameters + // + double threshold = 0.000000001*pow(double(iteration+3),agressiveness); + + // update mesh once in a while + if(iteration%5==0) + { + update_mesh(iteration); + } + + // clear dirty flag + loopi(0,triangles.size()) triangles[i].dirty=0; + + // target number of triangles reached ? Then break + if ((verbose) && (iteration%5==0)) { + printf("iteration %d - triangles %d threshold %g\n",iteration,triangle_count-deleted_triangles, threshold); + } + + // remove vertices & mark deleted triangles + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + if(t.err[3]>threshold) continue; + if(t.deleted == 1) continue; + if(t.deleted == -1) continue; + if(t.dirty) continue; + + loopj(0,3)if(t.err[j] deleted0,deleted1; + int triangle_count=triangles.size(); + //int iteration = 0; + //loop(iteration,0,100) + for (int iteration = 0; iteration < 9999; iteration ++) + { + // update mesh constantly + update_mesh(iteration); + // clear dirty flag + loopi(0,triangles.size()) triangles[i].dirty=0; + // + // All triangles with edges below the threshold will be removed + // + // The following numbers works well for most models. + // If it does not, try to adjust the 3 parameters + // + if (verbose) { + printf("lossless iteration %d\n", iteration); + } + + // remove vertices & mark deleted triangles + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + if(t.err[3]>threshold) continue; + if(t.deleted == 1) continue; + if(t.deleted == -1) continue; + if(t.dirty) continue; + + loopj(0,3)if(t.err[j] &deleted) + { + + loopk(0,v0.tcount) + { + Triangle &t=triangles[refs[v0.tstart+k].tid]; + if(t.deleted == 1)continue; + + int s=refs[v0.tstart+k].tvertex; + int id1=t.v[(s+1)%3]; + int id2=t.v[(s+2)%3]; + + if(id1==i1 || id2==i1) // delete ? + { + + deleted[k]=1; + continue; + } + vec3f d1 = vertices[id1].p-p; d1.normalize(); + vec3f d2 = vertices[id2].p-p; d2.normalize(); + if(fabs(d1.dot(d2))>0.999) return true; + vec3f n; + n.cross(d1,d2); + n.normalize(); + deleted[k]=0; + if(n.dot(t.n)<0.2) return true; + } + return false; + } + + // Update triangle connections and edge error after a edge is collapsed + + void update_triangles(int i0,Vertex &v,std::vector &deleted,int &deleted_triangles) + { + vec3f p; + loopk(0,v.tcount) + { + Ref &r=refs[v.tstart+k]; + Triangle &t=triangles[r.tid]; + if(t.deleted == 1)continue; + + if(deleted[k]) + { + t.deleted=1; + deleted_triangles++; + continue; + } + t.v[r.tvertex]=i0; + t.dirty=1; + t.err[0]=calculate_error(t.v[0],t.v[1],p); + t.err[1]=calculate_error(t.v[1],t.v[2],p); + t.err[2]=calculate_error(t.v[2],t.v[0],p); + t.err[3]=min(t.err[0],min(t.err[1],t.err[2])); + refs.push_back(r); + } + } + + // compact triangles, compute edge error and build reference list + + void update_mesh(int iteration) + { + if(iteration>0) // compact triangles + { + int dst=0; + loopi(0,triangles.size()) + if(triangles[i].deleted == 0 || triangles[i].deleted == -1) + { + triangles[dst++]=triangles[i]; + } + triangles.resize(dst); + } + // + // Init Quadrics by Plane & Edge Errors + // + // required at the beginning ( iteration == 0 ) + // recomputing during the simplification is not required, + // but mostly improves the result for closed meshes + // + if( iteration == 0 ) + { + loopi(0,vertices.size()) + vertices[i].q=SymetricMatrix(0.0); + + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + vec3f n,p[3]; + loopj(0,3) p[j]=vertices[t.v[j]].p; + n.cross(p[1]-p[0],p[2]-p[0]); + n.normalize(); + t.n=n; + loopj(0,3) vertices[t.v[j]].q = + vertices[t.v[j]].q+SymetricMatrix(n.x,n.y,n.z,-n.dot(p[0])); + } + loopi(0,triangles.size()) + { + // Calc Edge Error + Triangle &t=triangles[i];vec3f p; + loopj(0,3) t.err[j]=calculate_error(t.v[j],t.v[(j+1)%3],p); + t.err[3]=min(t.err[0],min(t.err[1],t.err[2])); + } + } + + // Init Reference ID list + loopi(0,vertices.size()) + { + vertices[i].tstart=0; + vertices[i].tcount=0; + } + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + loopj(0,3) vertices[t.v[j]].tcount++; + } + int tstart=0; + loopi(0,vertices.size()) + { + Vertex &v=vertices[i]; + v.tstart=tstart; + tstart+=v.tcount; + v.tcount=0; + } + + // Write References + refs.resize(triangles.size()*3); + loopi(0,triangles.size()) + { + Triangle &t=triangles[i]; + loopj(0,3) + { + Vertex &v=vertices[t.v[j]]; + refs[v.tstart+v.tcount].tid=i; + refs[v.tstart+v.tcount].tvertex=j; + v.tcount++; + } + } + + // Identify boundary : vertices[].border=0,1 + if( iteration == 0 ) + { + std::vector vcount,vids; + + loopi(0,vertices.size()) + vertices[i].border=0; + + loopi(0,vertices.size()) + { + Vertex &v=vertices[i]; + vcount.clear(); + vids.clear(); + loopj(0,v.tcount) + { + int k=refs[v.tstart+j].tid; + Triangle &t=triangles[k]; + loopk(0,3) + { + int ofs=0,id=t.v[k]; + while(ofs try to find best result + vec3f p1=vertices[id_v1].p; + vec3f p2=vertices[id_v2].p; + vec3f p3=(p1+p2)/2; + double error1 = vertex_error(q, p1.x,p1.y,p1.z); + double error2 = vertex_error(q, p2.x,p2.y,p2.z); + double error3 = vertex_error(q, p3.x,p3.y,p3.z); + error = min(error1, min(error2, error3)); + if (error1 == error) p_result=p1; + if (error2 == error) p_result=p2; + if (error3 == error) p_result=p3; + } + return error; + } +}; +/////////////////////////////////////////// diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 0000000..db967d1 --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,496 @@ +/* +Quickly generate 2.5D meshes from elevation models. +Copyright (C) 2018 Piero Toffanin + +dem2mesh is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +dem2mesh is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with dem2mesh. If not, see . +*/ +#include +#include +#include +#include +#include +#include +#include "CmdLineParser.h" +#include "Logger.h" +#include "Simplify.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; +} p; +size_t psize = sizeof(float) * 3; + +struct PlyFace{ + uint32_t p1; + uint32_t p2; + uint32_t p3; +} face; +size_t fsize = sizeof(uint32_t) * 3; + +int arr_width, arr_height; +int subdivisions; +int blockSizeX, blockSizeY; + +#define IS_BIG_ENDIAN (*(uint16_t *)"\0\xff" < 0x100) + +cmdLineParameter< char* > + InputFile( "inputFile" ) , + OutputFile( "outputFile" ); +cmdLineParameter< int > + MaxVertexCount( "maxVertexCount" ) , + MaxTileLength( "maxTileLength" ); +cmdLineReadable + Rtc ( "rtc" ), + Verbose( "verbose" ); + +cmdLineReadable* params[] = { + &InputFile , &OutputFile , &MaxVertexCount , &MaxTileLength, + &Rtc, &Verbose , + NULL +}; + +void help(char *ex){ + std::cout << "Usage: " << ex << std::endl + << "\t -" << InputFile.name << " " << std::endl + << "\t -" << OutputFile.name << " " << std::endl + << "\t [-" << MaxVertexCount.name << " (Default: 100000)]" << std::endl + << "\t [-" << MaxTileLength.name << " (Default: 1000)]" << std::endl + << "\t [-" << Rtc.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)); +} + +void writePly(const std::string &filename){ + // Start writing ply file + std::ofstream f (filename); + 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 " << Simplify::vertices.size() << std::endl + << "property float x" << std::endl + << "property float y" << std::endl + << "property float z" << std::endl + << "element face " << Simplify::triangles.size() << std::endl + << "property list uint8 uint32 vertex_indices" << std::endl + << "end_header" << std::endl; + + for(Simplify::Vertex &v : Simplify::vertices){ + p.x = static_cast(v.p.x); + p.y = static_cast(v.p.y); + p.z = static_cast(v.p.z); + f.write(reinterpret_cast(&p), psize); + } + + uint8_t three = 3; + for(Simplify::Triangle &t : Simplify::triangles){ + face.p1 = static_cast(t.v[0]); + face.p2 = static_cast(t.v[1]); + face.p3 = static_cast(t.v[2]); + + f.write(reinterpret_cast(&three), sizeof(three)); + f.write(reinterpret_cast(&face), fsize); + } + + f.close(); +} + +void writeBin(const std::string &filename, int blockWidth, int blockHeight){ + std::ofstream f (filename, std::ios::binary); + + unsigned long vsize = Simplify::vertices.size(); + unsigned long tsize = Simplify::triangles.size(); + + f.write(reinterpret_cast(&blockWidth), sizeof(blockWidth)); + f.write(reinterpret_cast(&blockHeight), sizeof(blockHeight)); + f.write(reinterpret_cast(&vsize), sizeof(vsize)); + f.write(reinterpret_cast(&tsize), sizeof(tsize)); + + for(Simplify::Vertex &v : Simplify::vertices){ + p.x = static_cast(v.p.x); + p.y = static_cast(v.p.y); + p.z = static_cast(v.p.z); + f.write(reinterpret_cast(&p), psize); + } + + for(Simplify::Triangle &t : Simplify::triangles){ + face.p1 = static_cast(t.v[0]); + face.p2 = static_cast(t.v[1]); + face.p3 = static_cast(t.v[2]); + f.write(reinterpret_cast(&face), fsize); + } + + f.close(); +} + + +// Keep track of edge points's vertex IDs +// During the merge step we need to replace +// duplicate points with existing point's vertex IDs +// Half the array maps to X coords +// the other half maps to Y coords +// 1 2 +// __|__|__ 3 +// __|__|__ 4 +// | | +// +// 1---2---3---4 +long *pointToVertexIdMap = nullptr; + +// (current vertex index) --> (existing vertex index) for duplicate points during merge +std::unordered_map vertexToVertexMap; + +void readBin(const std::string &filename, int blockX, int blockY){ + std::ifstream f (filename, std::ios::binary); + + int blockWidth, blockHeight; + unsigned long vcount, tcount; + unsigned long voffset = Simplify::vertices.size(); + float vertices[3]; + uint32_t triangles[3]; + + f.read(reinterpret_cast(&blockWidth), sizeof(blockWidth)); + f.read(reinterpret_cast(&blockHeight), sizeof(blockHeight)); + f.read(reinterpret_cast(&vcount), sizeof(vcount)); + f.read(reinterpret_cast(&tcount), sizeof(tcount)); + + int blockXPad = blockWidth - blockSizeX; + int blockYPad = blockHeight - blockSizeY; + int xOffset = blockX * blockSizeX - blockXPad; + int yOffset = blockY * blockSizeY - blockYPad; + + int ptvidMapSize = arr_height * (subdivisions - 1) + arr_width * (subdivisions - 1) + 1; + int ptvidYOffset = arr_height * (subdivisions - 1) + 1; + if (pointToVertexIdMap == nullptr){ + pointToVertexIdMap = new long[ptvidMapSize]; + memset(pointToVertexIdMap, -1, ptvidMapSize*sizeof(*pointToVertexIdMap)); + } + + int x, y; + for (unsigned long i = 0; i < vcount; i++){ + f.read(reinterpret_cast(&vertices), sizeof(float) * 3); + + x = static_cast(vertices[0]); + y = static_cast(vertices[1]); + + // Detect edge points + if ((blockX > 0 && x == xOffset) || (blockX < subdivisions - 1 && x == xOffset + blockWidth - 1)){ + if (pointToVertexIdMap[y + (x / blockSizeX) * arr_height] == -1){ + pointToVertexIdMap[y + (x / blockSizeX) * arr_height] = static_cast(i + voffset); + }else{ + // Already added from another block + // Keep a reference to the point from the other block + vertexToVertexMap[i + voffset] = pointToVertexIdMap[y + (x / blockSizeX) * arr_height]; + } + } + + else if ((blockY > 0 && y == yOffset) || (blockY < subdivisions - 1 && y == yOffset + blockHeight - 1)){ + if (pointToVertexIdMap[ptvidYOffset + x + (y / blockSizeY) * arr_width] == -1){ + pointToVertexIdMap[ptvidYOffset + x + (y / blockSizeY) * arr_width] = i + voffset; + }else{ + // Already added from another block + // Keep a reference to the point from the other block + vertexToVertexMap[i + voffset] = pointToVertexIdMap[ptvidYOffset + x + (y / blockSizeY) * arr_width]; + } + } + + Simplify::Vertex v; + v.p.x = vertices[0]; + v.p.y = vertices[1]; + v.p.z = vertices[2]; + + Simplify::vertices.push_back(v); + } + + for (unsigned long i = 0; i < tcount; i++){ + f.read(reinterpret_cast(&triangles), sizeof(uint32_t) * 3); + Simplify::Triangle t; + loopk(0, 3) t.v[k] = triangles[k] + voffset; + + // Check vertices for substitutions + loopk(0, 3) if (vertexToVertexMap.find(t.v[k]) != vertexToVertexMap.end()) t.v[k] = vertexToVertexMap[t.v[k]]; + + // Skip degenerate triangles + if (t.v[0] != t.v[1] && t.v[0] != t.v[2] && t.v[1] != t.v[2]){ + t.deleted = 0; + Simplify::triangles.push_back(t); + } + } +} + +void simplify(int target_count){ + unsigned long start_size = Simplify::triangles.size(); + const double AGRESSIVENESS = 5.0; + + Simplify::simplify_mesh(target_count, AGRESSIVENESS, Verbose.set); + if ( Simplify::triangles.size() >= start_size) { + std::cerr << "Unable to reduce mesh.\n"; + exit(1); + } +} + +// From grid X/Y to world X/Y +void transform(const BoundingBox &extent){ + double ext_width = extent.max.x - extent.min.x; + double ext_height = extent.max.y - extent.min.y; + + double center_x = 0.0; + double center_y = 0.0; + + if (Rtc.set){ + center_x = (extent.min.x + extent.max.x) / 2.0; + center_y = (extent.min.y + extent.max.y) / 2.0; + logWriter("Relative to center (RTC) coordinate set to (%f, %f)\n", center_x, center_y); + } + + for(Simplify::Vertex &v : Simplify::vertices){ + v.p.x = center_x - extent.min.x + (static_cast(v.p.x) / static_cast(arr_width)) * ext_width; + v.p.y = center_y - extent.max.y - (static_cast(v.p.y) / static_cast(arr_height)) * ext_height; + } +} + +int main(int argc, char **argv) { + cmdLineParse( argc-1 , &argv[1] , params ); + if( !InputFile.set || !OutputFile.set ) help(argv[0]); + if ( !MaxVertexCount.set ) MaxVertexCount.value = 100000; + if ( !MaxTileLength.set ) MaxTileLength.value = 1000; + + logWriter.verbose = Verbose.set; + logWriter.outputFile = "dem2mesh.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); + + unsigned long long int vertex_count = static_cast(arr_height) * + static_cast(arr_width); + + logWriter("Reading raster...\n"); + logWriter("Total vertices before simplification: %llu\n", vertex_count); + + GDALRasterBand *band = dataset->GetRasterBand(1); + + int qtreeLevels = 0; + int numBlocks = 1; + while(true){ + subdivisions = (int)pow(2, qtreeLevels); + numBlocks = subdivisions * subdivisions; + blockSizeX = arr_width / subdivisions; + blockSizeY = arr_height / subdivisions; + if (blockSizeX > MaxTileLength.value || blockSizeY > MaxTileLength.value){ + qtreeLevels++; + }else{ + break; + } + } + + int blockXPad = 0; // Blocks > 0 need to re-add a column for seamless meshing + int blockYPad = 0; // Blocks > 0 need to re-add a row for seamless meshing + + logWriter("Blocks depth: %d\n", qtreeLevels); + logWriter("Splitting area in %d\n", numBlocks); + logWriter("Block size is %d, %d\n", blockSizeX, blockSizeY); + + float *rasterData = new float[blockSizeX + 1]; + + for (int blockX = 0; blockX < subdivisions; blockX++){ + int xOffset = blockX * blockSizeX - blockXPad; + + for (int blockY = 0; blockY < subdivisions; blockY++){ + int yOffset = blockY * blockSizeY - blockYPad; + + Simplify::vertices.clear(); + Simplify::triangles.clear(); + + logWriter("Processing block (%d,%d)\n", blockX, blockY); + + for (int y = 0; y < blockSizeY + blockYPad; y++){ + if (band->RasterIO( GF_Read, xOffset, yOffset + y, blockSizeX + blockXPad, 1, + rasterData, blockSizeX + blockXPad, 1, GDT_Float32, 0, 0 ) == CE_Failure){ + std::cerr << "Cannot access raster data" << std::endl; + exit(1); + } + + for (int x = 0; x < blockSizeX + blockXPad; x++){ + Simplify::Vertex v; + v.p.x = xOffset + x; + v.p.y = yOffset + y; + v.p.z = rasterData[x]; + + Simplify::vertices.push_back(v); + } + } + + unsigned int cols = blockSizeX + blockXPad; + unsigned int rows = blockSizeY + blockYPad; + + for (unsigned int y = 0; y < rows - 1; y++){ + for (unsigned int x = 0; x < cols - 1; x++){ + Simplify::Triangle t1; + t1.v[0] = cols * (y + 1) + x; + t1.v[1] = cols * y + x + 1; + t1.v[2] = cols * y + x; + if (y == 0 || x == 0 || y == rows - 2 || x == cols - 2) t1.deleted = -1; // freeze + else t1.deleted = 0; + + Simplify::triangles.push_back(t1); + + Simplify::Triangle t2; + t2.v[0] = cols * (y + 1) + x; + t2.v[1] = cols * (y + 1) + x + 1; + t2.v[2] = cols * y + x + 1; + if (y == 0 || x == 0 || y == rows - 2 || x == cols - 2) t2.deleted = -1; // freeze + else t2.deleted = 0; + + Simplify::triangles.push_back(t2); + } + } + + int trianglesPerBlock = (MaxVertexCount.value * 2) / numBlocks; + + // If we have a merge step, + // overshoot the triangle count requirement + // since we'll simplify the final mesh anyway. + // This leads to more uniform meshes. + if (qtreeLevels > 0) trianglesPerBlock = trianglesPerBlock * 3 / 2; + + int target_count = std::min(trianglesPerBlock, static_cast(Simplify::triangles.size())); + + logWriter("Sampled %d faces, target is %d\n", static_cast(Simplify::triangles.size()), target_count); + logWriter("Simplifying...\n"); + + simplify(target_count); + + if (qtreeLevels == 0){ + transform(extent); + logWriter("Single quad tree level, saving to PLY\n"); + logWriter("Writing to file..."); + writePly(OutputFile.value); + }else{ + logWriter("Writing to binary file..."); + std::stringstream ss; + ss << OutputFile.value << "." << blockX << "-" << blockY << ".bin"; + writeBin(ss.str(), blockSizeX + blockXPad, blockSizeY + blockYPad); + } + + logWriter(" done!\n"); + + blockYPad = 1; + } + + blockYPad = 0; + blockXPad = 1; + } + + delete[] rasterData; + GDALClose(dataset); + + if (qtreeLevels > 0){ + // Merge + logWriter("Merge step...\n"); + + Simplify::vertices.clear(); + Simplify::triangles.clear(); + + for (int blockX = 0; blockX < subdivisions; blockX++){ + for (int blockY = 0; blockY < subdivisions; blockY++){ + std::stringstream ss; + ss << OutputFile.value << "." << blockX << "-" << blockY << ".bin"; + logWriter("Reading %s\n", ss.str().c_str()); + readBin(ss.str(), blockX, blockY); + + if (std::remove(ss.str().c_str()) != 0){ + logWriter("Error while deleting intermediate file: %s\n", ss.str().c_str()); + } + } + } + + // Cleanup + if (pointToVertexIdMap != nullptr){ + delete[] pointToVertexIdMap; + pointToVertexIdMap = nullptr; + } + vertexToVertexMap.clear(); + + logWriter("Simplifying final mesh...\n"); + int target_count = std::min(MaxVertexCount.value * 2, static_cast(Simplify::triangles.size())); + simplify(target_count); + transform(extent); + logWriter("Writing to file... "); + writePly(OutputFile.value); + logWriter(" done!\n"); + } + }else{ + std::cerr << "Cannot open " << InputFile.value << std::endl; + } +}