Disk offloading, max_threshold flag

Former-commit-id: 98c340ad90

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-12T10:52:19. -->
+<!-- Written by QtCreator 4.7.0, 2018-10-15T12:36:52. -->
 <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/decapado3/mesh_dsm.tif -outputFile /data/drone/decapado3/mesh_dsm.ply -verbose</value>
+    <value type="QString" key="CMakeProjectManager.CMakeRunConfiguration.Arguments">-inputFile /data/drone/cmparks/mesh_dsm.tif -outputFile /data/drone/cmparks/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/Simplify.h

@@ -307,7 +307,7 @@ namespace Simplify
 	//                 more iterations yield higher quality
 	//
 
-	void simplify_mesh(int target_count, double agressiveness=7, bool verbose=false)
+    void simplify_mesh(int target_count, double max_threshold, double agressiveness=7, bool verbose=false)
 	{
 		// init
 //		loopi(0,triangles.size())
@@ -324,24 +324,24 @@ namespace Simplify
         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);
+            if (threshold > max_threshold) break;
 
 			// update mesh once in a while
-			if(iteration%5==0)
+            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);
@@ -404,7 +404,7 @@ namespace Simplify
 		compact_mesh();
 	} //simplify_mesh()
 
-	void simplify_mesh_lossless(bool verbose=false)
+    void simplify_mesh_lossless(double threshold, bool verbose=false)
 	{
 		// init
 //		loopi(0,triangles.size()) triangles[i].deleted=0;
@@ -427,7 +427,6 @@ namespace Simplify
 			// 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);
 			}

modules/odm_dem2mesh/src/main.cpp

@@ -92,6 +92,114 @@ BoundingBox getExtent(GDALDataset *dataset){
     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<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);
+    }
+
+    f.close();
+}
+
+void writeBin(const std::string &filename){
+    std::ofstream f (filename, std::ios::binary);
+
+    unsigned long vsize = Simplify::vertices.size();
+    unsigned long tsize = Simplify::triangles.size();
+
+    f.write(reinterpret_cast<char *>(&vsize), sizeof(vsize));
+    f.write(reinterpret_cast<char *>(&tsize), sizeof(tsize));
+
+    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);
+    }
+
+    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 *>(&face), fsize);
+    }
+
+    f.close();
+}
+
+void readBin(const std::string &filename){
+    std::ifstream f (filename, std::ios::binary);
+
+    unsigned long vcount, tcount;
+    unsigned long voffset = Simplify::vertices.size();
+    float vertices[3];
+    uint32_t triangles[3];
+
+    f.read(reinterpret_cast<char *>(&vcount), sizeof(vcount));
+    f.read(reinterpret_cast<char *>(&tcount), sizeof(tcount));
+
+    for (unsigned long i = 0; i < vcount; i++){
+        f.read(reinterpret_cast<char *>(&vertices), sizeof(float) * 3);
+        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<char *>(&triangles), sizeof(uint32_t) * 3);
+        Simplify::Triangle t;
+        t.v[0] = triangles[0] + voffset;
+        t.v[1] = triangles[1] + voffset;
+        t.v[2] = triangles[2] + voffset;
+        Simplify::triangles.push_back(t);
+    }
+}
+
+void simplify(int target_count){
+    unsigned long start_size = Simplify::triangles.size();
+    const double AGRESSIVENESS = 5.0;
+    const double MAX_THRESHOLD = 0.10;
+
+    Simplify::simplify_mesh(target_count, MAX_THRESHOLD, AGRESSIVENESS, Verbose.set);
+    if ( Simplify::triangles.size() >= start_size) {
+        std::cerr << "Unable to reduce mesh.\n";
+        exit(1);
+    }
+}
+
 int main(int argc, char **argv) {
     cmdLineParse( argc-1 , &argv[1] , params );
     if( !InputFile.set || !OutputFile.set ) help(argv[0]);
@@ -119,16 +227,13 @@ int main(int argc, char **argv) {
         unsigned long long int vertex_count = static_cast<unsigned long long int>(arr_height) *
                                               static_cast<unsigned long long int>(arr_width);
 
-        // TODO:
-        // 2. Figure out bad_alloc in cmparks dataset
-
         logWriter("Reading raster...\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;
+        int stride = 8;
         while (vertex_count > INT_MAX){
             stride *= 2;
             vertex_count = static_cast<int>(std::ceil((arr_height / static_cast<double>(stride))) *
@@ -142,8 +247,9 @@ int main(int argc, char **argv) {
 
         GDALRasterBand *band = dataset->GetRasterBand(1);
 
-        int subdivisions = 2;
+        int qtreeLevels = 0;
-        int numBlocks = (int)pow(2, subdivisions);
+        int subdivisions = (int)pow(2, qtreeLevels);
+        int numBlocks = subdivisions * subdivisions;
         int blockSizeX = arr_width / subdivisions;
         int blockSizeY = arr_height / subdivisions;
         int blockXPad = 0; // Blocks > 0 need to re-add a column for seamless meshing
@@ -191,8 +297,8 @@ int main(int argc, char **argv) {
                         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 - 1 || x == cols -1) t1.deleted = -1; // freeze
+                        if (y == 0 || x == 0 || y == rows - 2 || x == cols - 2) t1.deleted = -1; // freeze
-                        else */t1.deleted = 0;
+                        else t1.deleted = 0;
 
                         Simplify::triangles.push_back(t1);
 
@@ -200,78 +306,74 @@ int main(int argc, char **argv) {
                         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 - 1 || x == cols -1) t2.deleted = -1; // freeze
+                        if (y == 0 || x == 0 || y == rows - 2 || x == cols - 2) t2.deleted = -1; // freeze
-                        else */t1.deleted = 0;
+                        else t2.deleted = 0;
 
                         Simplify::triangles.push_back(t2);
                     }
                 }
 
-                double agressiveness = 7.0;
+                int trianglesPerBlock = (MaxVertexCount.value * 2) / numBlocks;
-                int target_count = std::min((MaxVertexCount.value * 2) / numBlocks, static_cast<int>(Simplify::triangles.size()));
+
+                // 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<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(target_count);
-                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...");
+                if (qtreeLevels == 0){
-
+                    logWriter("Single quad tree level, saving to PLY\n");
-                // Start writing ply file
+                    logWriter("Writing to file...");
-                std::stringstream ss;
+                    writePly(OutputFile.value);
-                ss << OutputFile.value << "." << blockX << "-" << blockY << ".ply";
-                std::ofstream f (ss.str());
-                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;
+                    logWriter("Writing to binary file...");
-                }
+                    std::stringstream ss;
-
+                    ss << OutputFile.value << "." << blockX << "-" << blockY << ".bin";
-                f   << "element vertex " << Simplify::vertices.size() << std::endl
+                    writeBin(ss.str());
-                    << "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();
-
                 blockYPad = stride;
             }
 
             blockYPad = 0;
-            blockXPad = stride;
+            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());
+                }
+            }
+
+            logWriter("Simplifying final mesh...\n");
+            int target_count = std::min(MaxVertexCount.value * 2, static_cast<int>(Simplify::triangles.size()));
+            simplify(target_count);
+            logWriter("Writing to file... ");
+            writePly(OutputFile.value);
+            logWriter(" done!\n");
+        }
     }else{
         std::cerr << "Cannot open " << InputFile.value << std::endl;
     }