From c046309b7947c086d7c17480a8cd80e9a682a763 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 11 Oct 2018 11:31:49 -0400 Subject: [PATCH] Improved memory usage, using vtkGreedyTerrainDecimator for 2.5D mesh generation Former-commit-id: c3aedc3e07d57a369e0134f128b74cdd5b4c3017 --- modules/odm_dem2mesh/CMakeLists.txt | 4 +- modules/odm_dem2mesh/CMakeLists.txt.user | 4 +- modules/odm_dem2mesh/src/Simplify.h | 1028 ---------------------- modules/odm_dem2mesh/src/main.cpp | 201 ++--- opendm/mesh.py | 4 +- 5 files changed, 67 insertions(+), 1174 deletions(-) delete mode 100644 modules/odm_dem2mesh/src/Simplify.h diff --git a/modules/odm_dem2mesh/CMakeLists.txt b/modules/odm_dem2mesh/CMakeLists.txt index 85cf3561..d8498702 100644 --- a/modules/odm_dem2mesh/CMakeLists.txt +++ b/modules/odm_dem2mesh/CMakeLists.txt @@ -5,7 +5,9 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}) set (CMAKE_CXX_STANDARD 11) find_package(GDAL REQUIRED) +find_package(VTK REQUIRED) include_directories(${GDAL_INCLUDE_DIR}) +include(${VTK_USE_FILE}) # Add compiler options. add_definitions(-Wall -Wextra) @@ -16,4 +18,4 @@ aux_source_directory("./src" SRC_LIST) # Add exectuteable add_executable(${PROJECT_NAME} ${SRC_LIST}) -target_link_libraries(${PROJECT_NAME} ${GDAL_LIBRARY}) +target_link_libraries(${PROJECT_NAME} ${GDAL_LIBRARY} ${VTK_LIBRARIES}) diff --git a/modules/odm_dem2mesh/CMakeLists.txt.user b/modules/odm_dem2mesh/CMakeLists.txt.user index 70290431..dab68ab7 100644 --- a/modules/odm_dem2mesh/CMakeLists.txt.user +++ b/modules/odm_dem2mesh/CMakeLists.txt.user @@ -1,6 +1,6 @@ - + EnvironmentId @@ -452,7 +452,7 @@ 13 14 - -inputFile /data/drone/TBays/mesh_dsm.tif -outputFile /data/drone/TBays/mesh_dsm.ply -verbose + -inputFile /data/drone/sheffield_cross/mesh_dsm.tif -outputFile /data/drone/sheffield_cross/mesh_dsm.ply -verbose /data/OpenDroneMap/modules/build-odm_dem2mesh-Desktop-Default 2 diff --git a/modules/odm_dem2mesh/src/Simplify.h b/modules/odm_dem2mesh/src/Simplify.h deleted file mode 100644 index 63d8fb5d..00000000 --- a/modules/odm_dem2mesh/src/Simplify.h +++ /dev/null @@ -1,1028 +0,0 @@ -///////////////////////////////////////////// -// -// 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; - } - -}; - -vec3f barycentric(const vec3f &p, const vec3f &a, const vec3f &b, const vec3f &c){ - vec3f v0 = b-a; - vec3f v1 = c-a; - vec3f v2 = p-a; - double d00 = v0.dot(v0); - double d01 = v0.dot(v1); - double d11 = v1.dot(v1); - double d20 = v2.dot(v0); - double d21 = v2.dot(v1); - double denom = d00*d11-d01*d01; - double v = (d11 * d20 - d01 * d21) / denom; - double w = (d00 * d21 - d01 * d20) / denom; - double u = 1.0 - v - w; - return vec3f(u,v,w); -} - -vec3f interpolate(const vec3f &p, const vec3f &a, const vec3f &b, const vec3f &c, const vec3f attrs[3]) -{ - vec3f bary = barycentric(p,a,b,c); - vec3f out = vec3f(0,0,0); - out = out + attrs[0] * bary.x; - out = out + attrs[1] * bary.y; - out = out + attrs[2] * bary.z; - return out; -} - -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 - enum Attributes { - NONE, - NORMAL = 2, - TEXCOORD = 4, - COLOR = 8 - }; - struct Triangle { int v[3];double err[4];int deleted,dirty,attr;vec3f n;vec3f uvs[3];int material; }; - struct Vertex { vec3f p;int tstart,tcount;SymetricMatrix q;int border;}; - struct Ref { int tid,tvertex; }; - std::vector triangles; - std::vector vertices; - std::vector refs; - std::string mtllib; - std::vector materials; - - // 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_uvs(int i0,const Vertex &v,const vec3f &p,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; - - // update mesh once in a while - if(iteration%5==0) - { - 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 - // - double threshold = 0.000000001*pow(double(iteration+3),agressiveness); - - // 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) 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 - // - double threshold = DBL_EPSILON; //1.0E-3 EPS; - 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) 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)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_uvs - - void update_uvs(int i0,const Vertex &v,const vec3f &p,std::vector &deleted) - { - loopk(0,v.tcount) - { - Ref &r=refs[v.tstart+k]; - Triangle &t=triangles[r.tid]; - if(t.deleted)continue; - if(deleted[k])continue; - vec3f p1=vertices[t.v[0]].p; - vec3f p2=vertices[t.v[1]].p; - vec3f p3=vertices[t.v[2]].p; - t.uvs[r.tvertex] = interpolate(p,p1,p2,p3,t.uvs); - } - } - - // 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)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) - { - 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; - } - - char *trimwhitespace(char *str) - { - char *end; - - // Trim leading space - while(isspace((unsigned char)*str)) str++; - - if(*str == 0) // All spaces? - return str; - - // Trim trailing space - end = str + strlen(str) - 1; - while(end > str && isspace((unsigned char)*end)) end--; - - // Write new null terminator - *(end+1) = 0; - - return str; - } - - //Option : Load OBJ - void load_obj(const char* filename, bool process_uv=false){ - vertices.clear(); - triangles.clear(); - //printf ( "Loading Objects %s ... \n",filename); - FILE* fn; - if(filename==NULL) return ; - if((char)filename[0]==0) return ; - if ((fn = fopen(filename, "rb")) == NULL) - { - printf ( "File %s not found!\n" ,filename ); - return; - } - char line[1000]; - memset ( line,0,1000 ); - int vertex_cnt = 0; - int material = -1; - std::map material_map; - std::vector uvs; - std::vector > uvMap; - - while(fgets( line, 1000, fn ) != NULL) - { - Vertex v; - vec3f uv; - - if (strncmp(line, "mtllib", 6) == 0) - { - mtllib = trimwhitespace(&line[7]); - } - if (strncmp(line, "usemtl", 6) == 0) - { - std::string usemtl = trimwhitespace(&line[7]); - if (material_map.find(usemtl) == material_map.end()) - { - material_map[usemtl] = materials.size(); - materials.push_back(usemtl); - } - material = material_map[usemtl]; - } - - if ( line[0] == 'v' && line[1] == 't' ) - { - if ( line[2] == ' ' ) - if(sscanf(line,"vt %lf %lf", - &uv.x,&uv.y)==2) - { - uv.z = 0; - uvs.push_back(uv); - } else - if(sscanf(line,"vt %lf %lf %lf", - &uv.x,&uv.y,&uv.z)==3) - { - uvs.push_back(uv); - } - } - else if ( line[0] == 'v' ) - { - if ( line[1] == ' ' ) - if(sscanf(line,"v %lf %lf %lf", - &v.p.x, &v.p.y, &v.p.z)==3) - { - vertices.push_back(v); - } - } - int integers[9]; - if ( line[0] == 'f' ) - { - Triangle t; - bool tri_ok = false; - bool has_uv = false; - - if(sscanf(line,"f %d %d %d", - &integers[0],&integers[1],&integers[2])==3) - { - tri_ok = true; - }else - if(sscanf(line,"f %d// %d// %d//", - &integers[0],&integers[1],&integers[2])==3) - { - tri_ok = true; - }else - if(sscanf(line,"f %d//%d %d//%d %d//%d", - &integers[0],&integers[3], - &integers[1],&integers[4], - &integers[2],&integers[5])==6) - { - tri_ok = true; - }else - if(sscanf(line,"f %d/%d/%d %d/%d/%d %d/%d/%d", - &integers[0],&integers[6],&integers[3], - &integers[1],&integers[7],&integers[4], - &integers[2],&integers[8],&integers[5])==9) - { - tri_ok = true; - has_uv = true; - } - else - { - printf("unrecognized sequence\n"); - printf("%s\n",line); - while(1); - } - if ( tri_ok ) - { - t.v[0] = integers[0]-1-vertex_cnt; - t.v[1] = integers[1]-1-vertex_cnt; - t.v[2] = integers[2]-1-vertex_cnt; - t.attr = 0; - - if ( process_uv && has_uv ) - { - std::vector indices; - indices.push_back(integers[6]-1-vertex_cnt); - indices.push_back(integers[7]-1-vertex_cnt); - indices.push_back(integers[8]-1-vertex_cnt); - uvMap.push_back(indices); - t.attr |= TEXCOORD; - } - - t.material = material; - //geo.triangles.push_back ( tri ); - triangles.push_back(t); - //state_before = state; - //state ='f'; - } - } - } - - if ( process_uv && uvs.size() ) - { - loopi(0,triangles.size()) - { - loopj(0,3) - triangles[i].uvs[j] = uvs[uvMap[i][j]]; - } - } - - fclose(fn); - - //printf("load_obj: vertices = %lu, triangles = %lu, uvs = %lu\n", vertices.size(), triangles.size(), uvs.size() ); - } // load_obj() - - // Optional : Store as OBJ - - void write_obj(const char* filename) - { - FILE *file=fopen(filename, "w"); - int cur_material = -1; - bool has_uv = (triangles.size() && (triangles[0].attr & TEXCOORD) == TEXCOORD); - - if (!file) - { - printf("write_obj: can't write data file \"%s\".\n", filename); - exit(0); - } - if (!mtllib.empty()) - { - fprintf(file, "mtllib %s\n", mtllib.c_str()); - } - loopi(0,vertices.size()) - { - //fprintf(file, "v %lf %lf %lf\n", vertices[i].p.x,vertices[i].p.y,vertices[i].p.z); - fprintf(file, "v %g %g %g\n", vertices[i].p.x,vertices[i].p.y,vertices[i].p.z); //more compact: remove trailing zeros - } - if (has_uv) - { - loopi(0,triangles.size()) if(!triangles[i].deleted) - { - fprintf(file, "vt %g %g\n", triangles[i].uvs[0].x, triangles[i].uvs[0].y); - fprintf(file, "vt %g %g\n", triangles[i].uvs[1].x, triangles[i].uvs[1].y); - fprintf(file, "vt %g %g\n", triangles[i].uvs[2].x, triangles[i].uvs[2].y); - } - } - int uv = 1; - loopi(0,triangles.size()) if(!triangles[i].deleted) - { - if (triangles[i].material != cur_material) - { - cur_material = triangles[i].material; - fprintf(file, "usemtl %s\n", materials[triangles[i].material].c_str()); - } - if (has_uv) - { - fprintf(file, "f %d/%d %d/%d %d/%d\n", triangles[i].v[0]+1, uv, triangles[i].v[1]+1, uv+1, triangles[i].v[2]+1, uv+2); - uv += 3; - } - else - { - fprintf(file, "f %d %d %d\n", triangles[i].v[0]+1, triangles[i].v[1]+1, triangles[i].v[2]+1); - } - //fprintf(file, "f %d// %d// %d//\n", triangles[i].v[0]+1, triangles[i].v[1]+1, triangles[i].v[2]+1); //more compact: remove trailing zeros - } - fclose(file); - } -}; -/////////////////////////////////////////// diff --git a/modules/odm_dem2mesh/src/main.cpp b/modules/odm_dem2mesh/src/main.cpp index 1a91b7ad..3e63acb7 100644 --- a/modules/odm_dem2mesh/src/main.cpp +++ b/modules/odm_dem2mesh/src/main.cpp @@ -4,7 +4,13 @@ #include #include "CmdLineParser.h" #include "Logger.h" -#include "Simplify.h" +#include +#include +#include +#include +#include +#include +#include #include "gdal_priv.h" #include "cpl_conv.h" // for CPLMalloc() @@ -27,25 +33,6 @@ typedef struct BoundingBox{ 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; - -float *rasterData; -int arr_width, arr_height; - -#define IS_BIG_ENDIAN (*(uint16_t *)"\0\xff" < 0x100) - cmdLineParameter< char* > InputFile( "inputFile" ) , OutputFile( "outputFile" ); @@ -99,146 +86,78 @@ int main(int argc, char **argv) { logWriter.verbose = Verbose.set; logArgs(params, logWriter); + GDALDataset *dataset; GDALAllRegister(); dataset = (GDALDataset *) GDALOpen( InputFile.value, GA_ReadOnly ); if( dataset != NULL ) { - arr_width = dataset->GetRasterXSize(); - arr_height = dataset->GetRasterYSize(); + int width = dataset->GetRasterXSize(); + int height = dataset->GetRasterYSize(); - logWriter("Raster Size is %dx%d\n", arr_width, arr_height); + logWriter("Raster Size is %dx%d\n", width, 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; + double extentX = extent.max.x - extent.min.x; + double extentY = extent.max.y - extent.min.y; - unsigned long long int vertex_count = static_cast(arr_height - 2) * - static_cast(arr_width - 2); + vtkSmartPointer image = + vtkSmartPointer::New(); + image->SetDimensions(width, height, 1); + image->AllocateScalars(VTK_FLOAT, 1); GDALRasterBand *band = dataset->GetRasterBand(1); + float *rasterData = new float[width]; - rasterData = new float[arr_width * arr_height]; + for (int y = 0; y < height; y++){ + if (band->RasterIO( GF_Read, 0, height - y - 1, width, 1, + rasterData, width, 1, GDT_Float32, 0, 0 ) == CE_Failure){ + std::cerr << "Cannot access raster data" << std::endl; + exit(1); + } - 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); - } - - // TODO: - // 2. Figure out bad_alloc in cmparks dataset - - logWriter("Sampling...\n"); - - // If the DSM is really large, we only sample a subset of it - // to remain within INT_MAX vertices. This does not happen often, - // but it's a safeguard to make sure we'll get an output and not - // overflow. - int stride = 1; - while (vertex_count > INT_MAX){ - stride += 1; - vertex_count = static_cast(std::ceil(((arr_height - 2) / static_cast(stride))) * - std::ceil(((arr_width - 2) / static_cast(stride)))); - } - - if (stride != 1){ - logWriter("Warning: DSM is large, will sample using stride value of %d\n", stride); - } - logWriter("Total vertices before simplification: %llu\n", vertex_count); - - for (int y = 1; y < arr_height - 1; y += stride){ - for (int x = 1; x < arr_width - 1; x += stride){ - Simplify::Vertex v; - v.p.x = extent.min.x + (static_cast(x) / static_cast(arr_width)) * ext_width; - v.p.y = extent.max.y - (static_cast(y) / static_cast(arr_height)) * ext_height; - v.p.z = rasterData[y * arr_width + x]; - - Simplify::vertices.push_back(v); + for (int x = 0; x < width; x++){ + (static_cast(image->GetScalarPointer(x,y,0)))[0] = rasterData[x]; } } - unsigned int cols = static_cast(std::ceil(((arr_width - 2) / static_cast(stride)))); - unsigned int rows = static_cast(std::ceil(((arr_height - 2) / static_cast(stride)))); - - 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; - t1.attr = 0; - t1.material = -1; - - 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; - t2.attr = 0; - t2.material = -1; - - Simplify::triangles.push_back(t2); - } - } - - double agressiveness = 7.0; - int target_count = std::min(MaxVertexCount.value * 2, static_cast(Simplify::triangles.size())); - - logWriter("Sampled %d faces, target is %d\n", static_cast(Simplify::triangles.size()), target_count); - logWriter("Simplifying...\n"); - - unsigned long start_size = Simplify::triangles.size(); - Simplify::simplify_mesh(target_count, agressiveness, true); - if ( Simplify::triangles.size() >= start_size) { - std::cerr << "Unable to reduce mesh.\n"; - exit(1); - } - - logWriter("Writing to file..."); - - // Start 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 " << 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); - } - - logWriter(" done!\n"); - - f.close(); GDALClose(dataset); + delete[] rasterData; + + vtkSmartPointer terrain = + vtkSmartPointer::New(); + terrain->SetErrorMeasureToNumberOfTriangles(); + terrain->SetNumberOfTriangles(MaxVertexCount.value * 2); // Approximate + terrain->SetInputData(image); + terrain->BoundaryVertexDeletionOn(); + + vtkSmartPointer subdivision = + vtkSmartPointer::New(); + subdivision->SetMaximumTriangleArea(50); + subdivision->SetInputConnection(terrain->GetOutputPort()); + + logWriter("OK\nTransform... "); + vtkSmartPointer transform = vtkSmartPointer::New(); + transform->Scale(extentX / width, extentY / height, 1); + transform->Translate(extent.min.x * (width / extentX), extent.min.y * (height / extentY), 0); + + vtkSmartPointer transformFilter = + vtkSmartPointer::New(); + transformFilter->SetInputConnection(subdivision->GetOutputPort()); + transformFilter->SetTransform(transform); + + logWriter("OK\nSaving mesh to file... "); + + vtkSmartPointer plyWriter = + vtkSmartPointer::New(); + plyWriter->SetFileName(OutputFile.value); + plyWriter->SetInputConnection(transformFilter->GetOutputPort()); + plyWriter->SetFileTypeToBinary(); + plyWriter->Write(); + + logWriter("OK\n"); }else{ std::cerr << "Cannot open " << InputFile.value << std::endl; } diff --git a/opendm/mesh.py b/opendm/mesh.py index b48bfcf1..fabec9ba 100644 --- a/opendm/mesh.py +++ b/opendm/mesh.py @@ -38,8 +38,8 @@ def create_25dmesh(inPointCloud, outMesh, dsm_resolution=0.05, depth=8, samples= mesh = dem_to_mesh(os.path.join(tmp_directory, 'mesh_dsm.tif'), outMesh, maxVertexCount, verbose) # Cleanup tmp - if os.path.exists(tmp_directory): - shutil.rmtree(tmp_directory) + # if os.path.exists(tmp_directory): + # shutil.rmtree(tmp_directory) return mesh