diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8bcbc5a..777e763 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -8,6 +8,11 @@ option(TNTN_DOWNLOAD_DEPS "download dependencies during cmake configure" ON)
 
 set(CMAKE_CXX_STANDARD 14)
 
+if(MSVC)
+    if(MSVC_VERSION LESS 1910)
+        message(FATAL_ERROR "Visual Studio 2017 version 15.0 or higher is required")
+    endif()
+endif()
 
 # get the current working branch
 execute_process(
@@ -31,14 +36,19 @@ find_package(GDAL REQUIRED)
 if(NOT GDAL_FOUND)
     message(FATAL_ERROR "GDAL not found, cannot proceed")
 endif()
-if(NOT GDAL_CONFIG)
-    message(FATAL_ERROR "gdal-config command not found (not in PATH?), cannot proceed")
-endif()
 
-execute_process(
-  COMMAND ${GDAL_CONFIG} --version
-  OUTPUT_VARIABLE SYSTEM_GDAL_VERSION
-)
+if(GDAL_VERSION)
+    set(SYSTEM_GDAL_VERSION ${GDAL_VERSION})
+else()
+    if(NOT GDAL_CONFIG)
+        message(FATAL_ERROR "gdal-config command not found (not in PATH?), cannot proceed")
+    endif()
+
+    execute_process(
+      COMMAND ${GDAL_CONFIG} --version
+      OUTPUT_VARIABLE SYSTEM_GDAL_VERSION
+    )
+endif()
 
 if(SYSTEM_GDAL_VERSION VERSION_LESS "2.2")
     message(FATAL_ERROR "GDAL version \"${SYSTEM_GDAL_VERSION}\" is too old, at least 2.2 is required")
diff --git a/include/tntn/File.h b/include/tntn/File.h
index 33151c8..b48a149 100644
--- a/include/tntn/File.h
+++ b/include/tntn/File.h
@@ -6,6 +6,10 @@
 #include <string>
 #include <vector>
 
+#ifdef _MSC_VER 
+#define ftello ftell
+#endif
+
 namespace tntn {
 
 class FileLike
diff --git a/include/tntn/FileFormat.h b/include/tntn/FileFormat.h
index 6253585..2767482 100644
--- a/include/tntn/FileFormat.h
+++ b/include/tntn/FileFormat.h
@@ -4,6 +4,11 @@
 #include <string.h>
 #include "tntn/MeshMode.h"
 
+#ifdef _MSC_VER 
+#define strncasecmp _strnicmp
+#define strcasecmp _stricmp
+#endif
+
 namespace tntn {
 
 class FileFormat
@@ -114,8 +119,8 @@ class FileFormat
         switch(m_value)
         {
             case OBJ: //fallthrough
-            case OFF: //fallthrough
-            case TERRAIN: return MeshMode::decomposed;
+            case OFF: return MeshMode::decomposed;
+            case TERRAIN: return MeshMode::triangles;
             default: return MeshMode::none;
         }
     }
diff --git a/include/tntn/Mesh.h b/include/tntn/Mesh.h
index 89e77dd..8f961d4 100644
--- a/include/tntn/Mesh.h
+++ b/include/tntn/Mesh.h
@@ -44,6 +44,7 @@ class Mesh
     SimpleRange<const Triangle*> triangles() const;
     SimpleRange<const Face*> faces() const;
     SimpleRange<const Vertex*> vertices() const;
+    SimpleRange<const Normal*> vertex_normals() const;
 
     void grab_triangles(std::vector<Triangle>& into);
     void grab_decomposed(std::vector<Vertex>& vertices, std::vector<Face>& faces);
@@ -62,6 +63,9 @@ class Mesh
     bool is_square() const;
     bool check_for_holes_in_square_mesh() const;
 
+    void compute_vertex_normals();
+    bool has_normals() const;
+
   private:
     bool semantic_equal_tri_tri(const Mesh& other) const;
     bool semantic_equal_dec_dec(const Mesh& other) const;
@@ -70,6 +74,7 @@ class Mesh
     std::vector<Vertex> m_vertices;
     std::vector<Face> m_faces;
     std::vector<Triangle> m_triangles;
+    std::vector<Normal> m_normals;
 
     void decompose_triangle(const Triangle& t);
 };
diff --git a/include/tntn/TerraMesh.h b/include/tntn/TerraMesh.h
index 819fe97..3210f85 100644
--- a/include/tntn/TerraMesh.h
+++ b/include/tntn/TerraMesh.h
@@ -29,7 +29,7 @@ class TerraMesh : public TerraBaseMesh
                             const double no_data_value);
 
   public:
-    void greedy_insert(double max_error);
+    void greedy_insert(double max_error, int max_iterations);
     void scan_triangle(dt_ptr t) override;
     std::unique_ptr<Mesh> convert_to_mesh();
 };
diff --git a/include/tntn/TileMaker.h b/include/tntn/TileMaker.h
index f007c2c..20e0ec5 100644
--- a/include/tntn/TileMaker.h
+++ b/include/tntn/TileMaker.h
@@ -10,10 +10,13 @@ namespace tntn {
 class TileMaker
 {
     std::unique_ptr<Mesh> m_mesh;
-
+	bool m_normals_enabled;
+	
   public:
-    TileMaker() : m_mesh(std::make_unique<Mesh>()) {}
+    TileMaker(bool normals_enabled=false) : m_mesh(std::make_unique<Mesh>()),
+											m_normals_enabled(normals_enabled) {}
 
+    bool normals_enabled() const;
     void setMeshWriter(MeshWriter* w);
     bool loadObj(const char* filename);
     void loadMesh(std::unique_ptr<Mesh> mesh);
diff --git a/include/tntn/dem2tintiles_workflow.h b/include/tntn/dem2tintiles_workflow.h
index b50f3c9..88eb0de 100644
--- a/include/tntn/dem2tintiles_workflow.h
+++ b/include/tntn/dem2tintiles_workflow.h
@@ -21,6 +21,7 @@ std::vector<Partition> create_partitions_for_zoom_level(const RasterDouble& dem,
 bool create_tiles_for_zoom_level(const RasterDouble& dem,
                                  const std::vector<Partition>& partitions,
                                  int zoom,
+                                 bool include_normals,
                                  const std::string& output_basedir,
                                  const double method_parameter,
                                  const std::string& meshing_method,
diff --git a/include/tntn/endianness.h b/include/tntn/endianness.h
index ab9c58a..8eadb14 100644
--- a/include/tntn/endianness.h
+++ b/include/tntn/endianness.h
@@ -8,7 +8,7 @@
 #    define TNTN_BIG_ENDIAN
 #elif defined(__BYTE_ORDER) && __BYTE_ORDER == __LITTLE_ENDIAN || defined(__LITTLE_ENDIAN__) || \
     defined(__ARMEL__) || defined(__THUMBEL__) || defined(__AARCH64EL__) || defined(_MIPSEL) || \
-    defined(__MIPSEL) || defined(__MIPSEL__)
+    defined(__MIPSEL) || defined(__MIPSEL__) || defined(_WIN32)
 #    define TNTN_LITTLE_ENDIAN
 #else
 #    error unknown architecture
diff --git a/include/tntn/geometrix.h b/include/tntn/geometrix.h
index 48322a9..6183fd6 100644
--- a/include/tntn/geometrix.h
+++ b/include/tntn/geometrix.h
@@ -18,6 +18,8 @@ typedef std::array<Vertex, 3> Triangle;
 typedef size_t VertexIndex; //0-based index into an array/vector of vertices
 typedef std::array<VertexIndex, 3> Face;
 
+typedef glm::dvec3 Normal;
+	
 struct Edge
 {
     Edge() = default;
diff --git a/include/tntn/terra_meshing.h b/include/tntn/terra_meshing.h
index bc4e983..c4b61b4 100644
--- a/include/tntn/terra_meshing.h
+++ b/include/tntn/terra_meshing.h
@@ -8,10 +8,10 @@
 
 namespace tntn {
 
-std::unique_ptr<Mesh> generate_tin_terra(std::unique_ptr<RasterDouble> raster, double max_error);
+std::unique_ptr<Mesh> generate_tin_terra(std::unique_ptr<RasterDouble> raster, double max_error, int max_iterations=0);
 
 std::unique_ptr<Mesh> generate_tin_terra(std::unique_ptr<SurfacePoints> surface_points,
-                                         double max_error);
-std::unique_ptr<Mesh> generate_tin_terra(const SurfacePoints& surface_points, double max_error);
+                                         double max_error, int max_iterations=0);
+std::unique_ptr<Mesh> generate_tin_terra(const SurfacePoints& surface_points, double max_error, int max_iterations=0);
 
 } //namespace tntn
diff --git a/src/Mesh.cpp b/src/Mesh.cpp
index fe478f5..f517165 100644
--- a/src/Mesh.cpp
+++ b/src/Mesh.cpp
@@ -5,11 +5,11 @@
 #include <limits>
 #include <unordered_map>
 #include <unordered_set>
+#include <glm/gtx/normal.hpp>
 
 #include "tntn/tntn_assert.h"
 #include "tntn/logging.h"
 #include "tntn/geometrix.h"
-
 namespace tntn {
 
 void Mesh::clear()
@@ -17,6 +17,7 @@ void Mesh::clear()
     m_triangles.clear();
     m_vertices.clear();
     m_faces.clear();
+    m_normals.clear();
 }
 
 void Mesh::clear_triangles()
@@ -36,6 +37,7 @@ Mesh Mesh::clone() const
     out.m_faces = m_faces;
     out.m_vertices = m_vertices;
     out.m_triangles = m_triangles;
+    out.m_normals = m_normals;
     return out;
 }
 
@@ -217,6 +219,15 @@ SimpleRange<const Vertex*> Mesh::vertices() const
     return {m_vertices.data(), m_vertices.data() + m_vertices.size()};
 }
 
+SimpleRange<const Normal*> Mesh::vertex_normals() const
+{
+    if(m_normals.empty())
+    {
+        return {nullptr, nullptr};
+    }
+    return {m_normals.data(), m_normals.data() + m_normals.size()};
+}
+
 void Mesh::grab_triangles(std::vector<Triangle>& into)
 {
     into.clear();
@@ -710,4 +721,42 @@ bool Mesh::check_tin_properties() const
     return true;
 }
 
+bool Mesh::has_normals() const
+{
+    return !m_normals.empty();
+}
+
+void Mesh::compute_vertex_normals()
+{
+    using namespace glm;
+
+    std::vector<dvec3> face_normals;
+    face_normals.reserve(m_faces.size());
+    m_normals.resize(m_vertices.size());
+
+    auto face_normal = [](const Triangle& t) {
+        return normalize(cross(t[0] - t[2], t[1] - t[2]));
+    };
+
+    std::transform(m_triangles.begin(), m_triangles.end(), face_normals.begin(), face_normal);
+
+    const int faces_count = m_faces.size();
+
+    for(int f = 0; f < faces_count; f++)
+    {
+        const auto& face = m_faces[f];
+        const auto& face_normal = face_normals[f];
+
+        for(int v = 0; v < 3; v++)
+        {
+            const int vertex_id = face[v];
+            m_normals[vertex_id] += face_normal;
+        }
+    }
+
+    std::transform(m_normals.begin(), m_normals.end(), m_normals.begin(), [](const auto& n) {
+        return normalize(n);
+    });
+}
+
 } //namespace tntn
diff --git a/src/QuantizedMeshIO.cpp b/src/QuantizedMeshIO.cpp
index a865878..881e1da 100644
--- a/src/QuantizedMeshIO.cpp
+++ b/src/QuantizedMeshIO.cpp
@@ -91,6 +91,12 @@ struct QuantizedMeshHeader
     // double HorizonOcclusionPointZ;
 };
 
+enum QuantizedMeshExtensionList
+{
+    QMExtensionNormals = 1,
+    QMExtensionWaterMask = 2
+};
+
 namespace detail {
 
 uint16_t zig_zag_encode(int16_t i)
@@ -317,6 +323,37 @@ bool point_to_ecef(Vertex& p)
     return (bool)tr->Transform(1, &p.x, &p.y, &p.z);
 }
 
+void oct_encode(const Normal& v, uint8_t out[2])
+{
+    const glm::dvec2 p = v.xy * (1.0 / (abs(v.x) + abs(v.y) + abs(v.z)));
+    const auto signNotZero = glm::dvec2(v.x >= 0 ? 1.0 : -1.0, v.y > 0 ? 1.0 : -1.0);
+    const glm::dvec2 abs_v = p.yx;
+    glm::dvec2 ov = (v.z <= 0.0) ? ((1.0 - glm::abs(abs_v)) * signNotZero) : p;
+
+    ov = (ov + 1.0) / 2.0;
+
+    out[0] = static_cast<uint8_t>(ov.x * 255);
+    out[1] = static_cast<uint8_t>(ov.y * 255);
+}
+
+void write_normals(BinaryIO& bio,
+                   BinaryIOErrorTracker& e,
+                   const SimpleRange<const Normal*>& normals)
+{
+    // Write extension header
+    bio.write_byte(QMExtensionNormals, e);
+    bio.write_uint32(normals.distance() * 2, e);
+
+    uint8_t oct_normal[2];
+
+    for(auto n = normals.begin; n != normals.end; n++)
+    {
+        oct_encode(*n, oct_normal);
+        bio.write_byte(oct_normal[0], e);
+        bio.write_byte(oct_normal[1], e);
+    }
+}
+
 bool write_mesh_as_qm(const std::shared_ptr<FileLike>& f,
                       const Mesh& m,
                       const BBox3D& bbox,
@@ -496,6 +533,11 @@ bool write_mesh_as_qm(const std::shared_ptr<FileLike>& f,
         write_indices<uint32_t>(bio, e, northlings);
     }
 
+    if(m.has_normals())
+    {
+        write_normals(bio, e, m.vertex_normals());
+    }
+
     if(e.has_error())
     {
         TNTN_LOG_ERROR("{} in file {}", e.to_string(), f->name());
@@ -503,6 +545,7 @@ bool write_mesh_as_qm(const std::shared_ptr<FileLike>& f,
     }
 
     TNTN_LOG_INFO("writer log: {}", log.to_string());
+
     return true;
 }
 
diff --git a/src/RasterOverviews.cpp b/src/RasterOverviews.cpp
index 5d07c71..ccc0bbf 100644
--- a/src/RasterOverviews.cpp
+++ b/src/RasterOverviews.cpp
@@ -22,35 +22,34 @@ RasterOverviews::RasterOverviews(UniqueRasterPointer input_raster, int min_zoom,
 // Guesses (numerically) maximal zoom level from a raster resolution
 int RasterOverviews::guess_max_zoom_level(double resolution)
 {
-	// pixel_size_z0 is a magic number that is number of meters per pixel on zoom level 0, given tile size is 256 pixels.
-	// This number is approximate and does not account for latitude, i.e. uses pixel size on equator.
-	// The formula is: earth circumference * 2 * pi / 256
-	// "Real" number can be off up to 30% depending on the latitude. More details here: https://msdn.microsoft.com/en-us/library/aa940990.aspx
-	const double pixel_size_z0 = 156543.04; 
+    // pixel_size_z0 is a magic number that is number of meters per pixel on zoom level 0, given tile size is 256 pixels.
+    // This number is approximate and does not account for latitude, i.e. uses pixel size on equator.
+    // The formula is: earth circumference * 2 * pi / 256
+    // "Real" number can be off up to 30% depending on the latitude. More details here: https://msdn.microsoft.com/en-us/library/aa940990.aspx
+    const double pixel_size_z0 = 156543.04;
     return static_cast<int>(round(log2(pixel_size_z0 / resolution)));
 }
 
 // Guesses (numerically) minimal zoom level from a raster resolution and it's size
 int RasterOverviews::guess_min_zoom_level(int max_zoom_level)
 {
-	// This constant is an arbitrary number representing some minimal size to which the raster can be downsized when 'zooming out' a map.
-	// Math is simple:
-	// 2**max_zoom = raster_size; 2**min_zoom = 128
-	// Solve it in regards to min_zoom and there you have it.
-	const int MINIMAL_RASTER_SIZE = 128;
+    // This constant is an arbitrary number representing some minimal size to which the raster can be downsized when 'zooming out' a map.
+    // Math is simple:
+    // 2**max_zoom = raster_size; 2**min_zoom = 128
+    // Solve it in regards to min_zoom and there you have it.
+    const int MINIMAL_RASTER_SIZE = 128;
     const double quotient = MINIMAL_RASTER_SIZE * (1 << max_zoom_level);
 
     int raster_width = m_base_raster->get_width();
     int raster_height = m_base_raster->get_height();
 
-    TNTN_LOG_DEBUG("guess_min_zoom_level: raster_width: {}, raster_height: {}",
-                       raster_width,
-                       raster_height);
+    TNTN_LOG_DEBUG(
+        "guess_min_zoom_level: raster_width: {}, raster_height: {}", raster_width, raster_height);
 
     int zoom_x = static_cast<int>(floor(log2(quotient / raster_width)));
     int zoom_y = static_cast<int>(floor(log2(quotient / raster_height)));
 
-	// Cap the resulting value so it won't go below 0 and we wouldn't end up with negative zoom levels
+    // Cap the resulting value so it won't go below 0 and we wouldn't end up with negative zoom levels
     return std::max(0, std::min(zoom_x, zoom_y));
 }
 
@@ -60,12 +59,21 @@ void RasterOverviews::compute_zoom_levels()
     m_estimated_min_zoom = guess_min_zoom_level(m_estimated_max_zoom);
 
     m_min_zoom = std::max(m_min_zoom, m_estimated_min_zoom);
-    m_max_zoom = std::min(std::max(0, m_max_zoom), m_estimated_max_zoom);
+
+    if(m_max_zoom < 0 || m_max_zoom > m_estimated_max_zoom)
+    {
+        m_max_zoom = m_estimated_max_zoom;
+    }
 
     if(m_max_zoom < m_min_zoom)
     {
         std::swap(m_min_zoom, m_max_zoom);
     }
+
+    TNTN_LOG_INFO(
+        "After checking with data, tiles will be generated in a range between {} and {}",
+        m_min_zoom,
+        m_max_zoom);
 }
 
 bool RasterOverviews::next(RasterOverview& overview)
@@ -88,11 +96,12 @@ bool RasterOverviews::next(RasterOverview& overview)
     overview.resolution = output_raster->get_cell_size();
     overview.raster = std::move(output_raster);
 
-    TNTN_LOG_DEBUG("Generated next overview at zoom {}, window size {}, min zoom level {}, max zoom level {}",
-                   m_current_zoom + 1,
-                   window_size,
-                   m_min_zoom,
-                   m_max_zoom);
+    TNTN_LOG_DEBUG(
+        "Generated next overview at zoom {}, window size {}, min zoom level {}, max zoom level {}",
+        m_current_zoom + 1,
+        window_size,
+        m_min_zoom,
+        m_max_zoom);
 
     return true;
 }
diff --git a/src/TerraMesh.cpp b/src/TerraMesh.cpp
index c11a2bf..6645f3b 100644
--- a/src/TerraMesh.cpp
+++ b/src/TerraMesh.cpp
@@ -13,7 +13,7 @@
 namespace tntn {
 namespace terra {
 
-void TerraMesh::greedy_insert(double max_error)
+void TerraMesh::greedy_insert(double max_error, int max_iterations)
 {
     m_max_error = max_error;
     m_counter = 0;
@@ -56,12 +56,14 @@ void TerraMesh::greedy_insert(double max_error)
         t = t->getLink();
     }
 
+    int iterations_count = 0;
     // Iterate until the error threshold is met
     while(!m_candidates.empty())
     {
         Candidate candidate = m_candidates.grab_greatest();
 
         if(candidate.importance < m_max_error) continue;
+        if(max_iterations && iterations_count >= max_iterations) continue;
 
         // Skip if the candidate is not the latest
         if(m_token.value(candidate.y, candidate.x) != candidate.token) continue;
@@ -70,6 +72,7 @@ void TerraMesh::greedy_insert(double max_error)
 
         //TNTN_LOG_DEBUG("inserting point: ({}, {}, {})", candidate.x, candidate.y, candidate.z);
         this->insert(glm::dvec2(candidate.x, candidate.y), candidate.triangle);
+        iterations_count++;
     }
 
     TNTN_LOG_INFO("finished greedy insertion");
diff --git a/src/TileMaker.cpp b/src/TileMaker.cpp
index 8615d1a..5e26961 100644
--- a/src/TileMaker.cpp
+++ b/src/TileMaker.cpp
@@ -41,6 +41,11 @@ void TileMaker::loadMesh(std::unique_ptr<Mesh> mesh)
     m_mesh = std::move(mesh);
 }
 
+bool TileMaker::normals_enabled() const
+{
+    return m_normals_enabled;
+}
+
 // Dump a tile into an terrain tile in format determined by a MeshWriter
 bool TileMaker::dumpTile(int tx, int ty, int zoom, const char* filename, MeshWriter& mesh_writer)
 {
@@ -82,10 +87,8 @@ bool TileMaker::dumpTile(int tx, int ty, int zoom, const char* filename, MeshWri
     {
         for(int i = 0; i < 3; i++)
         {
-            if(triangle[i].z < tileSpaceBbox.min.z)
-                tileSpaceBbox.min.z = triangle[i].z;
-            if(triangle[i].z > tileSpaceBbox.max.z)
-                tileSpaceBbox.max.z = triangle[i].z;
+            if(triangle[i].z < tileSpaceBbox.min.z) tileSpaceBbox.min.z = triangle[i].z;
+            if(triangle[i].z > tileSpaceBbox.max.z) tileSpaceBbox.max.z = triangle[i].z;
         }
     }
 
@@ -120,6 +123,11 @@ bool TileMaker::dumpTile(int tx, int ty, int zoom, const char* filename, MeshWri
     tileMesh.from_triangles(std::move(trianglesInTile));
     tileMesh.generate_decomposed();
 
+    if(normals_enabled())
+    {
+        tileMesh.compute_vertex_normals();
+    }
+
     return mesh_writer.write_mesh_to_file(filename, tileMesh, tileSpaceBbox);
 }
 
diff --git a/src/benchmark_workflow.cpp b/src/benchmark_workflow.cpp
index ae34870..7a5539f 100644
--- a/src/benchmark_workflow.cpp
+++ b/src/benchmark_workflow.cpp
@@ -158,7 +158,7 @@ static bool prepare_output_directory(const fs::path& output_dir, const bool resu
         if(e)
         {
             TNTN_LOG_ERROR("unable to create output directory {} with message: {}",
-                           output_dir.c_str(),
+                           output_dir.string().c_str(),
                            e.message());
             return false;
         }
@@ -167,13 +167,13 @@ static bool prepare_output_directory(const fs::path& output_dir, const bool resu
     //make sure it is a directory
     if(!fs::is_directory(output_dir, e))
     {
-        TNTN_LOG_ERROR("output directory {} is not a directory", output_dir.c_str());
+        TNTN_LOG_ERROR("output directory {} is not a directory", output_dir.string().c_str());
         return false;
     }
     if(e)
     {
         TNTN_LOG_ERROR("filesystem error while checking output directory {} with message: {}",
-                       output_dir.c_str(),
+                       output_dir.string().c_str(),
                        e.message());
         return false;
     }
@@ -183,13 +183,13 @@ static bool prepare_output_directory(const fs::path& output_dir, const bool resu
         //make sure it is empty
         if(!fs::is_empty(output_dir, e))
         {
-            TNTN_LOG_ERROR("output directory {} is not empty", output_dir.c_str());
+            TNTN_LOG_ERROR("output directory {} is not empty", output_dir.string().c_str());
             return false;
         }
         if(e)
         {
             TNTN_LOG_ERROR("filesystem error while checking output directory {} with message: {}",
-                           output_dir.c_str(),
+                           output_dir.string().c_str(),
                            e.message());
             return false;
         }
@@ -662,13 +662,13 @@ static bool write_mesh_as_obj_and_off(const fs::path& parametrization_subdir,
 
     bool rc = true;
     File f;
-    if(f.open(obj_filename.c_str(), File::OM_RWC))
+    if(f.open(obj_filename.string().c_str(), File::OM_RWC))
     {
         rc = rc && write_mesh_as_obj(f, m);
         f.close();
     }
 
-    if(f.open(off_filename.c_str(), File::OM_RWC))
+    if(f.open(off_filename.string().c_str(), File::OM_RWC))
     {
         rc = rc && write_mesh_as_off(f, m);
         f.close();
@@ -688,7 +688,7 @@ static bool write_raster_as_asc_with_prefix(const fs::path& parametrization_subd
     raster_filename = parametrization_subdir / raster_filename;
 
     File f;
-    if(!f.open(raster_filename.c_str(), File::OM_RWC))
+    if(!f.open(raster_filename.string().c_str(), File::OM_RWC))
     {
         return false;
     }
@@ -941,7 +941,7 @@ static bool run_all_dem2tin_method_benchmarks_on_single_file(
 
             //create is_done_file to signal to later resume runs
             File f;
-            f.open(is_done_file.c_str(), File::OM_RWC);
+            f.open(is_done_file.string().c_str(), File::OM_RWC);
         }
     }
     return true;
@@ -1006,7 +1006,7 @@ bool run_dem2tin_method_benchmarks(const std::string& output_dir,
     const auto benchmark_csv_filename = output_dir_p / "tin_terrain_benchmarks.csv";
     if(resume && fs::exists(benchmark_csv_filename))
     {
-        if(!csv_file->open(benchmark_csv_filename.c_str(), File::OM_RW))
+        if(!csv_file->open(benchmark_csv_filename.string().c_str(), File::OM_RW))
         {
             TNTN_LOG_ERROR("unable to open CSV output file {} for resume",
                            benchmark_csv_filename.string());
@@ -1015,7 +1015,7 @@ bool run_dem2tin_method_benchmarks(const std::string& output_dir,
     }
     else
     {
-        if(!csv_file->open(benchmark_csv_filename.c_str(), File::OM_RWC))
+        if(!csv_file->open(benchmark_csv_filename.string().c_str(), File::OM_RWC))
         {
             TNTN_LOG_ERROR("unable to create CSV output file {} for writing",
                            benchmark_csv_filename.string());
diff --git a/src/cmd.cpp b/src/cmd.cpp
index c303e12..84ec665 100644
--- a/src/cmd.cpp
+++ b/src/cmd.cpp
@@ -55,6 +55,7 @@ static int subcommand_dem2tintiles(bool need_help,
         ("max-error", po::value<double>(), "max error parameter when using terra or zemlya method")
         ("step", po::value<int>()->default_value(1), "grid spacing in pixels when using dense method")
         ("output-format", po::value<std::string>()->default_value("terrain"), "output tiles in terrain (quantized mesh) or obj")
+		("normals", po::bool_switch()->default_value(false), "generate normals for terrain")
 #if defined(TNTN_USE_ADDONS) && TNTN_USE_ADDONS
         ("method", po::value<std::string>()->default_value("terra"), "meshing algorithm. one of: terra, zemlya, curvature or dense")
         ("threshold", po::value<double>(), "threshold when using curvature method");
@@ -176,6 +177,8 @@ static int subcommand_dem2tintiles(bool need_help,
         throw po::error(std::string("unknown method ") + meshing_method);
     }
 
+    bool write_normals = local_varmap["normals"].as<bool>();
+
     RasterOverviews overviews(std::move(input_raster), min_zoom, max_zoom);
 
     RasterOverview overview;
@@ -214,6 +217,7 @@ static int subcommand_dem2tintiles(bool need_help,
         if(!create_tiles_for_zoom_level(*overview.raster,
                                         partitions,
                                         zoom_level,
+                                        write_normals,
                                         output_basedir,
                                         max_error,
                                         meshing_method,
@@ -283,6 +287,7 @@ static int subcommand_dem2tin(bool need_help,
         ("output", po::value<std::string>(), "output filename")
         ("output-format", po::value<std::string>()->default_value("auto"), "output file format, can be any of: auto, obj, off, terrain (quantized mesh), json/geojson")
         ("max-error", po::value<double>(), "max error parameter when using terra or zemlya method")
+        ("max-iterations", po::value<int>()->default_value(0), "max add point iterations when using terra or zemlya method")
         ("step", po::value<int>(), "grid spacing in pixels when using dense method")
 #if defined(TNTN_USE_ADDONS) && TNTN_USE_ADDONS
         ("threshold", po::value<double>(), "threshold when using curvature method")
@@ -359,6 +364,11 @@ static int subcommand_dem2tin(bool need_help,
 
     const auto t_start = std::chrono::high_resolution_clock::now();
 
+    if(method != "terra" && local_varmap.count("max-iterations"))
+    {
+        throw po::error("parameter --max-iterations implemented for terra method only");
+    }
+
     if(method == "terra" || method == "zemlya")
     {
         double max_error = raster->get_cell_size();
@@ -370,7 +380,7 @@ static int subcommand_dem2tin(bool need_help,
         if("terra" == method)
         {
             TNTN_LOG_INFO("performing terra meshing...");
-            mesh = generate_tin_terra(std::move(raster), max_error);
+            mesh = generate_tin_terra(std::move(raster), max_error, local_varmap["max-iterations"].as<int>());
         }
         else if("zemlya" == method)
         {
diff --git a/src/dem2tintiles_workflow.cpp b/src/dem2tintiles_workflow.cpp
index 4310145..42fa5a9 100644
--- a/src/dem2tintiles_workflow.cpp
+++ b/src/dem2tintiles_workflow.cpp
@@ -66,6 +66,7 @@ std::vector<Partition> create_partitions_for_zoom_level(const RasterDouble& dem,
 bool create_tiles_for_zoom_level(const RasterDouble& dem,
                                  const std::vector<Partition>& partitions,
                                  int zoom,
+                                 bool include_normals,
                                  const std::string& output_basedir,
                                  const double method_parameter,
                                  const std::string& meshing_method,
@@ -128,7 +129,7 @@ bool create_tiles_for_zoom_level(const RasterDouble& dem,
         }
 
         // Cut the TIN into tiles
-        TileMaker tm;
+        TileMaker tm(include_normals);
         tm.loadMesh(std::move(mesh));
 
         fs::create_directory(fs::path(output_basedir));
@@ -147,7 +148,7 @@ bool create_tiles_for_zoom_level(const RasterDouble& dem,
                 auto file_path =
                     tile_dir / (std::to_string(ty) + "." + mesh_writer.file_extension());
 
-                if(!tm.dumpTile(tx, ty, zoom, file_path.c_str(), mesh_writer))
+                if(!tm.dumpTile(tx, ty, zoom, file_path.string().c_str(), mesh_writer))
                 {
                     TNTN_LOG_ERROR("error dumping tile z:{} x:{} y:{}", zoom, tx, ty);
                     return false;
diff --git a/src/raster_tools.cpp b/src/raster_tools.cpp
index 547f886..dd8b80c 100644
--- a/src/raster_tools.cpp
+++ b/src/raster_tools.cpp
@@ -324,22 +324,6 @@ static double subsample_raster_3x3(const RasterDouble& src,
 
 static constexpr int MAX_AVERAGING_SAMPLES = 64;
 
-static inline double average(const std::array<double, MAX_AVERAGING_SAMPLES>& to_average,
-                             const int avg_count)
-{
-    if(avg_count == 0)
-    {
-        return NAN;
-    }
-    double sum = 0;
-    for(int i = 0; i < MAX_AVERAGING_SAMPLES; i++)
-    {
-        sum += to_average[i];
-    }
-    const double avg = sum / avg_count;
-    return avg;
-}
-
 double raster_tools::sample_nearest_valid_avg(const RasterDouble& src,
                                               const unsigned int _row,
                                               const unsigned int _column,
@@ -365,19 +349,18 @@ double raster_tools::sample_nearest_valid_avg(const RasterDouble& src,
         return z;
     }
 
-    std::array<double, MAX_AVERAGING_SAMPLES> to_average;
-    std::fill(to_average.begin(), to_average.end(), 0.0);
+    double avg = 0.0;
     int avg_count = 0;
 
-    auto putpixel = [row, column, w, h, no_data_value, &src, &to_average, &avg_count](int x,
+    auto putpixel = [row, column, w, h, no_data_value, &src, &avg, &avg_count](int x,
                                                                                       int y) {
         const int64_t dest_r = row + y;
         const int64_t dest_c = column + x;
         double z = subsample_raster_3x3(src, no_data_value, w, h, dest_r, dest_c);
         if(!is_no_data(z, no_data_value))
         {
-            to_average[avg_count] = z;
-            avg_count++;
+          avg_count++;
+          avg = avg + ( z - avg ) / (double) avg_count;
         }
     };
 
@@ -417,13 +400,7 @@ double raster_tools::sample_nearest_valid_avg(const RasterDouble& src,
         }
     }
 
-    //short path for only one sample
-    if(avg_count == 1)
-    {
-        return to_average[0];
-    }
-
-    return average(to_average, avg_count);
+    return ( avg_count == 0 ? NAN : avg );
 }
 
 } // namespace tntn
diff --git a/src/terra_meshing.cpp b/src/terra_meshing.cpp
index 18cdd4f..9d8f11d 100644
--- a/src/terra_meshing.cpp
+++ b/src/terra_meshing.cpp
@@ -5,34 +5,34 @@
 
 namespace tntn {
 
-std::unique_ptr<Mesh> generate_tin_terra(std::unique_ptr<RasterDouble> raster, double max_error)
+std::unique_ptr<Mesh> generate_tin_terra(std::unique_ptr<RasterDouble> raster, double max_error, int max_iterations)
 {
     TNTN_ASSERT(raster != nullptr);
     terra::TerraMesh g;
     g.load_raster(std::move(raster));
-    g.greedy_insert(max_error);
+    g.greedy_insert(max_error, max_iterations);
     return g.convert_to_mesh();
 }
 
 std::unique_ptr<Mesh> generate_tin_terra(std::unique_ptr<SurfacePoints> surface_points,
-                                         double max_error)
+                                         double max_error, int max_iterations)
 {
     auto raster = surface_points->to_raster();
     surface_points.reset();
 
     terra::TerraMesh g;
     g.load_raster(std::move(raster));
-    g.greedy_insert(max_error);
+    g.greedy_insert(max_error, max_iterations);
     return g.convert_to_mesh();
 }
 
-std::unique_ptr<Mesh> generate_tin_terra(const SurfacePoints& surface_points, double max_error)
+std::unique_ptr<Mesh> generate_tin_terra(const SurfacePoints& surface_points, double max_error, int max_iterations)
 {
     auto raster = surface_points.to_raster();
 
     terra::TerraMesh g;
     g.load_raster(std::move(raster));
-    g.greedy_insert(max_error);
+    g.greedy_insert(max_error, max_iterations);
     return g.convert_to_mesh();
 }
 
diff --git a/test/src/Mesh_tests.cpp b/test/src/Mesh_tests.cpp
index 2da3c01..eb8e3e6 100644
--- a/test/src/Mesh_tests.cpp
+++ b/test/src/Mesh_tests.cpp
@@ -322,5 +322,34 @@ TEST_CASE("Mesh::check_tin_properties detects unreferenced vertices")
     CHECK(!m.check_tin_properties());
 }
 
+TEST_CASE("Mesh::compute_normals normals")
+{
+    Mesh m;
+
+    m.from_decomposed(
+        {
+            {0, 0, 1},
+            {1, 0, 2},
+            {1, 1, 3},
+            {0, 1, 4},
+            {0.5, -1, 5},
+        },
+        {
+            {{0, 1, 2}}, //ccw
+            {{0, 2, 3}}, //ccw
+            {{0, 4, 1}}, //ccw
+        });
+    m.generate_triangles();
+    m.compute_vertex_normals();
+
+    const auto vnormals = m.vertex_normals();
+    CHECK(std::distance(vnormals.begin, vnormals.end) == 5);
+
+    for(auto n = vnormals.begin; n != vnormals.end; ++n)
+    {
+        CHECK((glm::length(*n) - 1.0) < 0.00001);
+    }
+}
+
 } //namespace unittests
 } //namespace tntn