Bug fixes, updated README

README.md

@@ -8,10 +8,10 @@ Fast generation of 2.5D meshes from elevation models.
 
 ## Dependencies
 
-GDAL is the only dependency. To install it run:
+GDAL and OpenMP are the only dependencies. To install it run:
 
 ```
-sudo apt-get install -y libgdal-dev 
+sudo apt-get install -y libgdal-dev libomp-dev
 ```
 
 ## Building
@@ -40,6 +40,7 @@ make
 | -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`. | |
 | -aggressiveness | Value between `1` and `10` that specifies how "aggressive" the mesh simplification process should be at each iteration. Higher values simplify the mesh more aggressively but can decrease the fidelity of the mesh. Defaults to `5`. | |
 | -bandNum | Raster band # to use. Defaults to `1`. | |
+| -maxConcurrency | Maximum number of threads to use. Defaults to all CPUs available. | |
 | -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. | |
 

src/Simplify.h

@@ -390,12 +390,12 @@ namespace Simplify
 					// not flipped, so remove edge
 					v0.p=p;
 					v0.q=v1.q+v0.q;
-					int tstart=refs.size();
+                    int tstart=refs[thread]->size();
 
                     update_triangles(i0,v0,deleted0,deleted_triangles,thread);
                     update_triangles(i0,v1,deleted1,deleted_triangles,thread);
 
-					int tcount=refs.size()-tstart;
+                    int tcount=refs[thread]->size()-tstart;
 
 					if(tcount<=v0.tcount)
 					{
@@ -513,7 +513,7 @@ namespace Simplify
                 loopj(0,3) (*vertices[thread])[t.v[j]].q =
                     (*vertices[thread])[t.v[j]].q+SymetricMatrix(n.x,n.y,n.z,-n.dot(p[0]));
 			}
-            loopi(0,(*triangles[thread]).size())
+            loopi(0,triangles[thread]->size())
 			{
 				// Calc Edge Error
                 Triangle &t=(*triangles[thread])[i];vec3f p;

src/main.cpp

@@ -382,12 +382,17 @@ int main(int argc, char **argv) {
         logWriter("Splitting area in %d\n", numBlocks);
         logWriter("Block size is %d, %d\n", blockSizeX, blockSizeY);
 
-        float *rasterData = new float[blockSizeX + 1];
-
         logWriter("Concurrency set to %d\n", MaxConcurrency.value);
         Simplify::allocate(MaxConcurrency.value);
         omp_set_num_threads(MaxConcurrency.value);
 
+        int rasterDataBlocks = std::min(MaxConcurrency.value, numBlocks);
+        logWriter("Allocating %d raster data blocks of %d bytes\n", rasterDataBlocks, sizeof(float) * (blockSizeX + 1));
+        float *rasterData = new float[rasterDataBlocks * (blockSizeX + 1)];
+
+        omp_lock_t readLock;
+        omp_init_lock(&readLock);
+
         #pragma omp parallel for collapse(2)
         for (int blockX = 0; blockX < subdivisions; blockX++){
             for (int blockY = 0; blockY < subdivisions; blockY++){
@@ -404,19 +409,19 @@ int main(int argc, char **argv) {
 
                 for (int y = 0; y < blockSizeY + blockYPad; y++){
 
-                    // TODO: needs lock
+                    omp_set_lock(&readLock);
                     if (band->RasterIO( GF_Read, xOffset, yOffset + y, blockSizeX + blockXPad, 1,
-                                        rasterData, blockSizeX + blockXPad, 1, GDT_Float32, 0, 0 ) == CE_Failure){
+                                        rasterData + t * (blockSizeX + 1), blockSizeX + blockXPad, 1, GDT_Float32, 0, 0 ) == CE_Failure){
                         std::cerr << "Cannot access raster data" << std::endl;
                         exit(EXIT_FAILURE);
                     }
-                    // TODO: END LOCK
+                    omp_unset_lock(&readLock);
 
                     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];
+                        v.p.z = (rasterData + t * (blockSizeX + 1))[x];
 
                         Simplify::vertices[t]->push_back(v);
                     }
@@ -519,7 +524,7 @@ int main(int argc, char **argv) {
             vertexToVertexMap.clear();
 
             logWriter("Simplifying final mesh...\n");
-            int target_count = std::min(MaxVertexCount.value * 2, static_cast<int>(Simplify::triangles.size()));
+            int target_count = std::min(MaxVertexCount.value * 2, static_cast<int>(Simplify::triangles[0]->size()));
             simplify(target_count, 0);
             transform(extent, 0);
             logWriter("Writing to file... ");