diff --git a/README.md b/README.md index 75980d3..ba87e70 100644 --- a/README.md +++ b/README.md @@ -2,9 +2,13 @@ Fast generation of 2.5D meshes from elevation models. +![image](https://user-images.githubusercontent.com/1951843/47350997-15d3da00-d685-11e8-8d9f-e394fc17859e.png) + +![image](https://user-images.githubusercontent.com/1951843/47351205-7e22bb80-d685-11e8-87c5-33b21ae05b75.png) + ## Dependencies -GDAL is the only requirement. To install it run: +GDAL is the only dependency. To install it run: ``` sudo apt-get install -y libgdal-dev @@ -24,6 +28,8 @@ make ./dem2mesh -rtc -verbose -inputFile dem.tif -outputFile mesh.ply ``` +:warning: The DEM should use a coordinate reference system (CRS) that has matching horizontal and vertical units (for example, UTM). If you are getting stretched or deformed results, double check your CRS and use `gdalwarp` to reproject the raster DEM. https://www.gdal.org/gdalwarp.html + ## Options | Flag | Description | Required | @@ -32,5 +38,6 @@ make | -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`. | | +| -bandNum | Raster band # to use. Defaults to `1`. | | | -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/main.cpp b/src/main.cpp index db967d1..0731779 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -71,13 +71,15 @@ cmdLineParameter< char* > OutputFile( "outputFile" ); cmdLineParameter< int > MaxVertexCount( "maxVertexCount" ) , - MaxTileLength( "maxTileLength" ); + MaxTileLength( "maxTileLength" ) , + BandNum ( "bandNum" ); cmdLineReadable Rtc ( "rtc" ), Verbose( "verbose" ); cmdLineReadable* params[] = { - &InputFile , &OutputFile , &MaxVertexCount , &MaxTileLength, + &InputFile , &OutputFile , + &MaxVertexCount , &MaxTileLength, &BandNum, &Rtc, &Verbose , NULL }; @@ -88,6 +90,7 @@ void help(char *ex){ << "\t -" << OutputFile.name << " " << std::endl << "\t [-" << MaxVertexCount.name << " (Default: 100000)]" << std::endl << "\t [-" << MaxTileLength.name << " (Default: 1000)]" << std::endl + << "\t [-" << BandNum.name << " (Default: 1)]" << std::endl << "\t [-" << Rtc.name << "]" << std::endl << "\t [-" << Verbose.name << "]" << std::endl; exit(EXIT_FAILURE); @@ -280,12 +283,17 @@ void readBin(const std::string &filename, int blockX, int blockY){ void simplify(int target_count){ unsigned long start_size = Simplify::triangles.size(); + if (target_count >= static_cast(start_size)){ + logWriter("No simplification needed\n"); + return; + } + 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); + exit(EXIT_FAILURE); } } @@ -314,6 +322,7 @@ int main(int argc, char **argv) { if( !InputFile.set || !OutputFile.set ) help(argv[0]); if ( !MaxVertexCount.set ) MaxVertexCount.value = 100000; if ( !MaxTileLength.set ) MaxTileLength.value = 1000; + if ( !BandNum.set ) BandNum.value = 1; logWriter.verbose = Verbose.set; logWriter.outputFile = "dem2mesh.txt"; @@ -332,14 +341,22 @@ int main(int argc, char **argv) { logWriter("Extent is (%f, %f), (%f, %f)\n", extent.min.x, extent.max.x, extent.min.y, extent.max.y); + GDALRasterBand *band = dataset->GetRasterBand(BandNum.value); + + int hasNoData = FALSE; + double nodata = band->GetNoDataValue(&hasNoData); + + if (hasNoData){ + logWriter("NoData value: %.18g\n", nodata); + } + logWriter("Description: %s\n", band->GetDescription()); + 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){ @@ -378,7 +395,7 @@ int main(int argc, char **argv) { 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); + exit(EXIT_FAILURE); } for (int x = 0; x < blockSizeX + blockXPad; x++){ @@ -403,7 +420,12 @@ int main(int argc, char **argv) { if (y == 0 || x == 0 || y == rows - 2 || x == cols - 2) t1.deleted = -1; // freeze else t1.deleted = 0; - Simplify::triangles.push_back(t1); + if (!hasNoData || + (Simplify::vertices[t1.v[0]].p.z != nodata && + Simplify::vertices[t1.v[1]].p.z != nodata && + Simplify::vertices[t1.v[2]].p.z != nodata)){ + Simplify::triangles.push_back(t1); + } Simplify::Triangle t2; t2.v[0] = cols * (y + 1) + x; @@ -412,7 +434,13 @@ int main(int argc, char **argv) { if (y == 0 || x == 0 || y == rows - 2 || x == cols - 2) t2.deleted = -1; // freeze else t2.deleted = 0; - Simplify::triangles.push_back(t2); + if (!hasNoData || + (Simplify::vertices[t2.v[0]].p.z != nodata && + Simplify::vertices[t2.v[1]].p.z != nodata && + Simplify::vertices[t2.v[2]].p.z != nodata)){ + Simplify::triangles.push_back(t2); + } + } } @@ -428,7 +456,6 @@ int main(int argc, char **argv) { logWriter("Sampled %d faces, target is %d\n", static_cast(Simplify::triangles.size()), target_count); logWriter("Simplifying...\n"); - simplify(target_count); if (qtreeLevels == 0){