Improved memory usage, using vtkGreedyTerrainDecimator for 2.5D mesh generation

Former-commit-id: c3aedc3e07

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})

modules/odm_dem2mesh/CMakeLists.txt.user

@@ -1,6 +1,6 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <!DOCTYPE QtCreatorProject>
-<!-- Written by QtCreator 4.7.0, 2018-10-08T20:39:47. -->
+<!-- Written by QtCreator 4.7.0, 2018-10-11T10:46:34. -->
 <qtcreator>
  <data>
   <variable>EnvironmentId</variable>
@@ -452,7 +452,7 @@
      <value type="int">13</value>
      <value type="int">14</value>
     </valuelist>
-    <value type="QString" key="CMakeProjectManager.CMakeRunConfiguration.Arguments">-inputFile /data/drone/TBays/mesh_dsm.tif -outputFile /data/drone/TBays/mesh_dsm.ply -verbose</value>
+    <value type="QString" key="CMakeProjectManager.CMakeRunConfiguration.Arguments">-inputFile /data/drone/sheffield_cross/mesh_dsm.tif -outputFile /data/drone/sheffield_cross/mesh_dsm.ply -verbose</value>
     <value type="QString" key="CMakeProjectManager.CMakeRunConfiguration.UserWorkingDirectory"></value>
     <value type="QString" key="CMakeProjectManager.CMakeRunConfiguration.UserWorkingDirectory.default">/data/OpenDroneMap/modules/build-odm_dem2mesh-Desktop-Default</value>
     <value type="int" key="PE.EnvironmentAspect.Base">2</value>

modules/odm_dem2mesh/src/main.cpp

@@ -4,7 +4,13 @@
 #include <cmath>
 #include "CmdLineParser.h"
 #include "Logger.h"
-#include "Simplify.h"
+#include <vtkSmartPointer.h>
+#include <vtkGreedyTerrainDecimation.h>
+#include <vtkTransform.h>
+#include <vtkTransformFilter.h>
+#include <vtkPLYWriter.h>
+#include <vtkAdaptiveSubdivisionFilter.h>
+#include <vtkImageData.h>
 
 #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();
+        int width = dataset->GetRasterXSize();
-        arr_height = dataset->GetRasterYSize();
+        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;
+        double extentX = extent.max.x - extent.min.x;
-        float ext_height = extent.max.y - extent.min.y;
+        double extentY = extent.max.y - extent.min.y;
 
-        unsigned long long int vertex_count = static_cast<unsigned long long int>(arr_height - 2) *
+        vtkSmartPointer<vtkImageData> image =
-                                              static_cast<unsigned long long int>(arr_width - 2);
+          vtkSmartPointer<vtkImageData>::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,
-        if (band->RasterIO( GF_Read, 0, 0, arr_width, arr_height,
+                                rasterData, width, 1, GDT_Float32, 0, 0 ) == CE_Failure){
-                            rasterData, arr_width, arr_height, GDT_Float32, 0, 0 ) == CE_Failure){
                 std::cerr << "Cannot access raster data" << std::endl;
                 exit(1);
             }
 
-        // TODO:
+            for (int x = 0; x < width; x++){
-        // 2. Figure out bad_alloc in cmparks dataset
+                (static_cast<float*>(image->GetScalarPointer(x,y,0)))[0] = rasterData[x];
-
-        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<int>(std::ceil(((arr_height - 2) / static_cast<double>(stride))) *
-                                            std::ceil(((arr_width - 2) / static_cast<double>(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<float>(x) / static_cast<float>(arr_width)) * ext_width;
-                v.p.y = extent.max.y - (static_cast<float>(y) / static_cast<float>(arr_height)) * ext_height;
-                v.p.z = rasterData[y * arr_width + x];
-
-                Simplify::vertices.push_back(v);
             }
         }
 
-        unsigned int cols = static_cast<unsigned int>(std::ceil(((arr_width - 2) / static_cast<double>(stride))));
-        unsigned int rows = static_cast<unsigned int>(std::ceil(((arr_height - 2) / static_cast<double>(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<int>(Simplify::triangles.size()));
-
-        logWriter("Sampled %d faces, target is %d\n", static_cast<int>(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<float>(v.p.x);
-            p.y = static_cast<float>(v.p.y);
-            p.z = static_cast<float>(v.p.z);
-            f.write(reinterpret_cast<char *>(&p), psize);
-        }
-
-        uint8_t three = 3;
-        for(Simplify::Triangle &t : Simplify::triangles){
-            face.p1 = static_cast<uint32_t>(t.v[0]);
-            face.p2 = static_cast<uint32_t>(t.v[1]);
-            face.p3 = static_cast<uint32_t>(t.v[2]);
-
-            f.write(reinterpret_cast<char *>(&three), sizeof(three));
-            f.write(reinterpret_cast<char *>(&face), fsize);
-        }
-
-        logWriter(" done!\n");
-
-        f.close();
         GDALClose(dataset);
+        delete[] rasterData;
+
+        vtkSmartPointer<vtkGreedyTerrainDecimation> terrain =
+              vtkSmartPointer<vtkGreedyTerrainDecimation>::New();
+        terrain->SetErrorMeasureToNumberOfTriangles();
+        terrain->SetNumberOfTriangles(MaxVertexCount.value * 2); // Approximate
+        terrain->SetInputData(image);
+        terrain->BoundaryVertexDeletionOn();
+
+        vtkSmartPointer<vtkAdaptiveSubdivisionFilter> subdivision =
+                vtkSmartPointer<vtkAdaptiveSubdivisionFilter>::New();
+        subdivision->SetMaximumTriangleArea(50);
+        subdivision->SetInputConnection(terrain->GetOutputPort());
+
+        logWriter("OK\nTransform... ");
+        vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
+        transform->Scale(extentX / width, extentY / height, 1);
+        transform->Translate(extent.min.x * (width / extentX), extent.min.y * (height / extentY), 0);
+
+        vtkSmartPointer<vtkTransformFilter> transformFilter =
+          vtkSmartPointer<vtkTransformFilter>::New();
+        transformFilter->SetInputConnection(subdivision->GetOutputPort());
+        transformFilter->SetTransform(transform);
+
+        logWriter("OK\nSaving mesh to file... ");
+
+        vtkSmartPointer<vtkPLYWriter> plyWriter =
+              vtkSmartPointer<vtkPLYWriter>::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;
     }

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):
+    # if os.path.exists(tmp_directory):
-        shutil.rmtree(tmp_directory)
+    #     shutil.rmtree(tmp_directory)
 
     return mesh