Terrain normals.fester heremaps#55

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

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

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,

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;

src/Mesh.cpp

@@ -5,18 +5,19 @@
 #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()
 {
     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

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,13 +533,19 @@ 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());
         return false;
     }
 
     TNTN_LOG_INFO("writer log: {}", log.to_string());
+
     return true;
 }
 

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

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

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,

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