From 987b6665e925cfed86f6ce02c832df798f242720 Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Wed, 28 Mar 2018 01:41:07 +0200 Subject: [PATCH 01/26] Quantized-mesh output support Provides two new features: - New quantized-mesh format ouput - New metadata file output --- src/BoundingSphere.hpp | 184 +++++++++++++++++ src/CMakeLists.txt | 9 + src/Coordinate3D.hpp | 154 ++++++++++++++ src/HeightFieldChunker.hpp | 412 +++++++++++++++++++++++++++++++++++++ src/Mesh.hpp | 82 ++++++++ src/MeshIterator.hpp | 67 ++++++ src/MeshTile.cpp | 375 +++++++++++++++++++++++++++++++++ src/MeshTile.hpp | 127 ++++++++++++ src/MeshTiler.cpp | 208 +++++++++++++++++++ src/MeshTiler.hpp | 80 +++++++ src/types.hpp | 4 +- tools/ctb-tile.cpp | 348 +++++++++++++++++++++++++++++-- 12 files changed, 2037 insertions(+), 13 deletions(-) create mode 100644 src/BoundingSphere.hpp create mode 100644 src/Coordinate3D.hpp create mode 100644 src/HeightFieldChunker.hpp create mode 100644 src/Mesh.hpp create mode 100644 src/MeshIterator.hpp create mode 100644 src/MeshTile.cpp create mode 100644 src/MeshTile.hpp create mode 100644 src/MeshTiler.cpp create mode 100644 src/MeshTiler.hpp diff --git a/src/BoundingSphere.hpp b/src/BoundingSphere.hpp new file mode 100644 index 0000000..a082f91 --- /dev/null +++ b/src/BoundingSphere.hpp @@ -0,0 +1,184 @@ +#ifndef BBSPHERE_HPP +#define BBSPHERE_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file BoundingSphere.hpp + * @brief This declares and defines the `BoundingSphere` class + * @author Alvaro Huarte + */ + +#include +#include + +#include "Coordinate3D.hpp" +#include "types.hpp" + +namespace ctb { + template class BoundingSphere; + template class BoundingBox; +} + +/// A spherical bounding region which is defined by a center point and a radius +template +class ctb::BoundingSphere { +public: + Coordinate3D center; ///< The center of the BoundingSphere + double radius; ///< The radius of the BoundingSphere + + /// Create an empty BoundingSphere + BoundingSphere() { + } + /// Create a BoundingSphere from the specified point stream + BoundingSphere(const std::vector> &points) { + fromPoints(points); + } + + /// Calculate the center and radius from the specified point stream + /// Based on Ritter's algorithm + void fromPoints(const std::vector> &points) { + const T MAX = std::numeric_limits::infinity(); + const T MIN = -std::numeric_limits::infinity(); + + Coordinate3D minPointX(MAX, MAX, MAX); + Coordinate3D minPointY(MAX, MAX, MAX); + Coordinate3D minPointZ(MAX, MAX, MAX); + Coordinate3D maxPointX(MIN, MIN, MIN); + Coordinate3D maxPointY(MIN, MIN, MIN); + Coordinate3D maxPointZ(MIN, MIN, MIN); + + // Store the points containing the smallest and largest component + // Used for the naive approach + for (int i = 0, icount = points.size(); i < icount; i++) { + const Coordinate3D &point = points[i]; + + if (point.x < minPointX.x) minPointX = point; + if (point.y < minPointY.y) minPointY = point; + if (point.z < minPointZ.z) minPointZ = point; + if (point.x > maxPointX.x) maxPointX = point; + if (point.y > maxPointY.y) maxPointY = point; + if (point.z > maxPointZ.z) maxPointZ = point; + } + + // Squared distance between each component min and max + T xSpan = (maxPointX - minPointX).magnitudeSquared(); + T ySpan = (maxPointY - minPointY).magnitudeSquared(); + T zSpan = (maxPointZ - minPointZ).magnitudeSquared(); + + Coordinate3D diameter1 = minPointX; + Coordinate3D diameter2 = maxPointX; + T maxSpan = xSpan; + if (ySpan > maxSpan) { + diameter1 = minPointY; + diameter2 = maxPointY; + maxSpan = ySpan; + } + if (zSpan > maxSpan) { + diameter1 = minPointZ; + diameter2 = maxPointZ; + maxSpan = zSpan; + } + + Coordinate3D ritterCenter = Coordinate3D( + (diameter1.x + diameter2.x) * 0.5, + (diameter1.y + diameter2.y) * 0.5, + (diameter1.z + diameter2.z) * 0.5 + ); + T radiusSquared = (diameter2 - ritterCenter).magnitudeSquared(); + T ritterRadius = std::sqrt(radiusSquared); + + // Initial center and radius (naive) get min and max box + Coordinate3D minBoxPt(minPointX.x, minPointY.y, minPointZ.z); + Coordinate3D maxBoxPt(maxPointX.x, maxPointY.y, maxPointZ.z); + Coordinate3D naiveCenter = (minBoxPt + maxBoxPt) * 0.5; + T naiveRadius = 0; + + for (int i = 0, icount = points.size(); i < icount; i++) { + const Coordinate3D &point = points[i]; + + // Find the furthest point from the naive center to calculate the naive radius. + T r = (point - naiveCenter).magnitude(); + if (r > naiveRadius) naiveRadius = r; + + // Make adjustments to the Ritter Sphere to include all points. + T oldCenterToPointSquared = (point - ritterCenter).magnitudeSquared(); + + if (oldCenterToPointSquared > radiusSquared) { + T oldCenterToPoint = std::sqrt(oldCenterToPointSquared); + ritterRadius = (ritterRadius + oldCenterToPoint) * 0.5; + + // Calculate center of new Ritter sphere + T oldToNew = oldCenterToPoint - ritterRadius; + ritterCenter.x = (ritterRadius * ritterCenter.x + oldToNew * point.x) / oldCenterToPoint; + ritterCenter.y = (ritterRadius * ritterCenter.y + oldToNew * point.y) / oldCenterToPoint; + ritterCenter.z = (ritterRadius * ritterCenter.z + oldToNew * point.z) / oldCenterToPoint; + } + } + + // Keep the naive sphere if smaller + if (naiveRadius < ritterRadius) { + center = ritterCenter; + radius = ritterRadius; + } + else { + center = naiveCenter; + radius = naiveRadius; + } + } +}; + +/// A bounding box which is defined by a pair of minimum and maximum coordinates +template +class ctb::BoundingBox { +public: + Coordinate3D min; ///< The min coordinate of the BoundingBox + Coordinate3D max; ///< The max coordinate of the BoundingBox + + /// Create an empty BoundingBox + BoundingBox() { + } + /// Create a BoundingBox from the specified point stream + BoundingBox(const std::vector> &points) { + fromPoints(points); + } + + /// Calculate the BBOX from the specified point stream + void fromPoints(const std::vector> &points) { + const T MAX = std::numeric_limits::infinity(); + const T MIN = -std::numeric_limits::infinity(); + min.x = MAX; + min.y = MAX; + min.z = MAX; + max.x = MIN; + max.y = MIN; + max.z = MIN; + + for (int i = 0, icount = points.size(); i < icount; i++) { + const Coordinate3D &point = points[i]; + + if (point.x < min.x) min.x = point.x; + if (point.y < min.y) min.y = point.y; + if (point.z < min.z) min.z = point.z; + if (point.x > max.x) max.x = point.x; + if (point.y > max.y) max.y = point.y; + if (point.z > max.z) max.z = point.z; + } + } +}; + +#endif /* BBSPHERE_HPP */ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8fca685..dd18c3b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -22,6 +22,8 @@ add_library(ctb SHARED GDALTiler.cpp TerrainTiler.cpp TerrainTile.cpp + MeshTiler.cpp + MeshTile.cpp GlobalMercator.cpp GlobalGeodetic.cpp) target_link_libraries(ctb ${GDAL_LIBRARIES} ${ZLIB_LIBRARIES}) @@ -29,13 +31,20 @@ target_link_libraries(ctb ${GDAL_LIBRARIES} ${ZLIB_LIBRARIES}) # Install libctb set(HEADERS Bounds.hpp + BoundingSphere.hpp Coordinate.hpp + Coordinate3D.hpp GDALTile.hpp GDALTiler.hpp GlobalGeodetic.hpp GlobalMercator.hpp Grid.hpp GridIterator.hpp + HeightFieldChunker.hpp + Mesh.hpp + MeshIterator.hpp + MeshTile.hpp + MeshTiler.hpp RasterIterator.hpp RasterTiler.hpp CTBException.hpp diff --git a/src/Coordinate3D.hpp b/src/Coordinate3D.hpp new file mode 100644 index 0000000..d4f6579 --- /dev/null +++ b/src/Coordinate3D.hpp @@ -0,0 +1,154 @@ +#ifndef COORDINATE3D_HPP +#define COORDINATE3D_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +#include +#include + +/** + * @file Coordinate3D.hpp + * @brief This declares and defines the `Coordinate3D` class + * @author Alvaro Huarte + */ + +namespace ctb { + template class Coordinate3D; +} + +/// A representation of a 3-dimensional point coordinate +template +class ctb::Coordinate3D { +public: + T x, y, z; ///< The x, y and z coordinate members + + /// Create an empty coordinate + Coordinate3D(): + x(0), + y(0), + z(0) + {} + + /// The const copy constructor + Coordinate3D(const Coordinate3D &other): + x(other.x), + y(other.y), + z(other.z) + {} + + /// Instantiate a coordinate from an x, y and z value + Coordinate3D(T x, T y, T z): + x(x), + y(y), + z(z) + {} + + /// Overload the equality operator + virtual bool + operator==(const Coordinate3D &other) const { + return x == other.x + && y == other.y + && z == other.z; + } + + /// Overload the assignment operator + virtual void + operator=(const Coordinate3D &other) { + x = other.x; + y = other.y; + z = other.z; + } + + /// Gets a read-only index-ordinate of the coordinate + inline virtual T operator[](const int index) const { + return (index == 0) ? x : (index == 1 ? y : z); + } + + /// Add operator + inline virtual Coordinate3D operator+(const Coordinate3D& other) const { + return Coordinate3D(x + other.x, y + other.y, z + other.z); + } + /// Subtract operator + inline virtual Coordinate3D operator-(const Coordinate3D& other) const { + return Coordinate3D(x - other.x, y - other.y, z - other.z); + } + /// Multiply operator + inline virtual Coordinate3D operator*(const Coordinate3D& other) const { + return Coordinate3D(x * other.x, y * other.y, z * other.z); + } + /// Divide operator + inline virtual Coordinate3D operator/(const Coordinate3D& other) const { + return Coordinate3D(x / other.x, y / other.y, z / other.z); + } + + /// AddByScalar operator + inline virtual Coordinate3D operator+(const T scalar) const { + return Coordinate3D(x + scalar, y + scalar, z + scalar); + } + /// SubtractByScalar operator + inline virtual Coordinate3D operator-(const T scalar) const { + return Coordinate3D(x - scalar, y - scalar, z - scalar); + } + /// MultiplyByScalar operator + inline virtual Coordinate3D operator*(const T scalar) const { + return Coordinate3D(x * scalar, y * scalar, z * scalar); + } + /// DivideByScalar operator + inline virtual Coordinate3D operator/(const T scalar) const { + return Coordinate3D(x / scalar, y / scalar, z / scalar); + } + + /// Cross product + inline Coordinate3D cross(const Coordinate3D &other) const { + return Coordinate3D((y * other.z) - (other.y * z), + (z * other.x) - (other.z * x), + (x * other.y) - (other.x * y)); + } + /// Dot product + inline double dot(const Coordinate3D &other) const { + return (x * other.x) + (y * other.y) + (z * other.z); + } + + // Cartesian3d methods + inline T magnitudeSquared(void) const { + return (x * x) + (y * y) + (z * z); + } + inline T magnitude(void) const { + return std::sqrt(magnitudeSquared()); + } + inline static Coordinate3D add(const Coordinate3D &p1, const Coordinate3D &p2) { + return p1 + p2; + } + inline static Coordinate3D subtract(const Coordinate3D &p1, const Coordinate3D &p2) { + return p1 - p2; + } + inline static T distanceSquared(const Coordinate3D &p1, const Coordinate3D &p2) { + T xdiff = p1.x - p2.x; + T ydiff = p1.y - p2.y; + T zdiff = p1.z - p2.z; + return (xdiff * xdiff) + (ydiff * ydiff) + (zdiff * zdiff); + } + inline static T distance(const Coordinate3D &p1, const Coordinate3D &p2) { + return std::sqrt(distanceSquared(p1, p2)); + } + inline Coordinate3D normalize(void) const { + T mgn = magnitude(); + return Coordinate3D(x / mgn, y / mgn, z / mgn); + } +}; + +#endif /* COORDINATE3D_HPP */ diff --git a/src/HeightFieldChunker.hpp b/src/HeightFieldChunker.hpp new file mode 100644 index 0000000..d8f1bdd --- /dev/null +++ b/src/HeightFieldChunker.hpp @@ -0,0 +1,412 @@ +#ifndef HEIGHTFIELDCHUNKER_HPP +#define HEIGHTFIELDCHUNKER_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file HeightFieldChucker.hpp + * @brief This declares and defines the `mesh` and `heightfield` classes + * @author Alvaro Huarte + */ + +#include + +#include "cpl_config.h" +#include "cpl_string.h" + +/** + * Helper classes to fill an irregular mesh of triangles from a heightmap tile. + * They are a refactored version from 'heightfield_chunker.cpp' from + * http://tulrich.com/geekstuff/chunklod.html + * + * It applies the Chunked LOD strategy by 'Thatcher Ulrich' + * preserving the input geometric error. + */ +namespace ctb { namespace chunk { + struct gen_state; + class mesh; + class heightfield; +} } + +/// Helper struct with state info for chunking a HeightField. +struct ctb::chunk::gen_state { + int my_buffer[2][2]; // x,y coords of the last two vertices emitted by the generate_ functions. + int activation_level; // for determining whether a vertex is enabled in the block we're working on + int ptr; // indexes my_buffer. + int previous_level; // for keeping track of level changes during recursion. + + /// Returns true if the specified vertex is in my_buffer. + bool in_my_buffer(int x, int y) const + { + return ((x == my_buffer[0][0]) && (y == my_buffer[0][1])) + || ((x == my_buffer[1][0]) && (y == my_buffer[1][1])); + } + + /// Sets the current my_buffer entry to (x,y) + void set_my_buffer(int x, int y) + { + my_buffer[ptr][0] = x; + my_buffer[ptr][1] = y; + } +}; + +/// An irregular mesh of triangles target of the HeightField chunker process. +class ctb::chunk::mesh { +public: + /// Clear all data. + virtual void clear() = 0; + + /// New vertex (Call this in strip order). + virtual void emit_vertex(const heightfield &heightfield, int x, int y) = 0; +}; + +/// Defines a regular grid of heigths or HeightField. +class ctb::chunk::heightfield { +public: + /// Constructor + heightfield(float *tileHeights, int tileSize) { + int tileCellSize = tileSize * tileSize; + + m_heights = tileHeights; + m_size = tileSize; + m_log_size = (int)(log2((float)m_size - 1) + 0.5); + + // Initialize level array. + m_levels = (int*)CPLMalloc(tileCellSize * sizeof(int)); + for (int i = 0; i < tileCellSize; i++) m_levels[i] = 255; + } + ~heightfield() { + clear(); + } + + /// Apply the specified maximum geometric error to fill the level info of the grid. + void applyGeometricError(double maximumGeometricError, bool smoothSmallZooms = false) { + int tileCellSize = m_size * m_size; + + // Initialize level array. + for (int i = 0; i < tileCellSize; i++) m_levels[i] = 255; + + // Run a view-independent L-K style BTT update on the heightfield, + // to generate error and activation_level values for each element. + update(maximumGeometricError, 0, m_size - 1, m_size - 1, m_size - 1, 0, 0); // sw half of the square + update(maximumGeometricError, m_size - 1, 0, 0, 0, m_size - 1, m_size - 1); // ne half of the square + + // Make sure our corner verts are activated. + int size = (m_size - 1); + activate(size, 0, 0); + activate(0, 0, 0); + activate(0, size, 0); + activate(size, size, 0); + + // Activate some vertices to smooth the shape of the Globe for small zooms. + if (smoothSmallZooms) { + int step = size / 16; + + for (int x = 0; x <= size; x += step) { + for (int y = 0; y <= size; y += step) { + if (get_level(x, y) == -1) activate(x, y, 0); + } + } + } + + // Propagate the activation_level values of verts to their parent verts, + // quadtree LOD style. Gives same result as L-K. + for (int i = 0; i < m_log_size; i++) { + propagate_activation_level(m_size >> 1, m_size >> 1, m_log_size - 1, i); + propagate_activation_level(m_size >> 1, m_size >> 1, m_log_size - 1, i); + } + } + + /// Clear all object data + void clear() { + m_heights = NULL; + m_size = 0; + m_log_size = 0; + + if (m_levels) { + CPLFree(m_levels); + m_levels = NULL; + } + } + + /// Return the array-index of specified coordinate, row order by default. + virtual int indexOfGridCoordinate(int x, int y) const { + return (y * m_size) + x; + } + /// Return the height of specified coordinate. + virtual float height(int x, int y) const { + int index = indexOfGridCoordinate(x, y); + return m_heights[index]; + } + + /// Generates the mesh using verts which are active at the given level. + void generateMesh(ctb::chunk::mesh &mesh, int level) { + int x0 = 0; + int y0 = 0; + + int size = (1 << m_log_size); + int half_size = size >> 1; + int cx = x0 + half_size; + int cy = y0 + half_size; + + // Start making the mesh. + mesh.clear(); + + // !!! This needs to be done in propagate, or something (too late now) !!! + // Make sure our corner verts are activated on this level. + activate(x0 + size, y0, level); + activate(x0, y0, level); + activate(x0, y0 + size, level); + activate(x0 + size, y0 + size, level); + + // Generate the mesh. + const heightfield &hf = *this; + generate_block(hf, mesh, level, m_log_size, x0 + half_size, y0 + half_size); + } + +private: + int m_size; // Number of cols and rows of this Heightmap + int m_log_size; // size == (1 << log_size) + 1 + float *m_heights; // grid of heights + int *m_levels; // grid of activation levels + + /// Return the activation level at (x, y) + int get_level(int x, int y) const + { + int index = indexOfGridCoordinate(x, y); + int level = m_levels[index]; + + if (x & 1) { + level = level >> 4; + } + level &= 0x0F; + if (level == 0x0F) return -1; + else return level; + } + /// Set the activation level at (x, y) + void set_level(int x, int y, int newlevel) + { + newlevel &= 0x0F; + int index = indexOfGridCoordinate(x, y); + int level = m_levels[index]; + + if (x & 1) { + level = (level & 0x0F) | (newlevel << 4); + } + else { + level = (level & 0xF0) | (newlevel); + } + m_levels[index] = level; + } + /// Sets the activation_level to the given level. + /// if it's greater than the vert's current activation level. + void activate(int x, int y, int level) + { + int current_level = get_level(x, y); + if (level > current_level) set_level(x, y, level); + } + + /// Given the triangle, computes an error value and activation level + /// for its base vertex, and recurses to child triangles. + bool update(double base_max_error, int ax, int ay, int rx, int ry, int lx, int ly) + { + bool res = false; + + // Compute the coordinates of this triangle's base vertex. + int dx = lx - rx; + int dy = ly - ry; + + if (std::abs(dx) <= 1 && std::abs(dy) <= 1) { + // We've reached the base level. There's no base + // vertex to update, and no child triangles to + // recurse to. + + return false; + } + + // base vert is midway between left and right verts. + int bx = rx + (dx >> 1); + int by = ry + (dy >> 1); + + float heightB = height(bx, by); + float heightL = height(lx, ly); + float heightR = height(rx, ry); + float error_B = std::abs(heightB - 0.5 * (heightL + heightR)); + + if (error_B >= base_max_error) { + // Compute the mesh level above which this vertex + // needs to be included in LOD meshes. + int activation_level = (int)std::floor(log2(error_B / base_max_error) + 0.5); + + // Force the base vert to at least this activation level. + activate(bx, by, activation_level); + res = true; + } + + // Recurse to child triangles. + update(base_max_error, bx, by, ax, ay, rx, ry); // base, apex, right + update(base_max_error, bx, by, lx, ly, ax, ay); // base, left, apex + + return res; + } + + /// Does a quadtree descent through the heightfield, in the square with + /// center at (cx, cz) and size of (2 ^ (level + 1) + 1). Descends + /// until level == target_level, and then propagates this square's + /// child center verts to the corresponding edge vert, and the edge + /// verts to the center. Essentially the quadtree meshing update + /// dependency graph as in my Gamasutra article. Must call this with + /// successively increasing target_level to get correct propagation. + void propagate_activation_level(int cx, int cy, int level, int target_level) + { + int half_size = 1 << level; + int quarter_size = half_size >> 1; + + if (level > target_level) { + // Recurse to children. + for (int j = 0; j < 2; j++) { + for (int i = 0; i < 2; i++) { + propagate_activation_level( + cx - quarter_size + half_size * i, + cy - quarter_size + half_size * j, + level - 1, target_level); + } + } + return; + } + + // We're at the target level. Do the propagation on this square. + if (level > 0) { + int lev = 0; + + // Propagate child verts to edge verts. + lev = get_level(cx + quarter_size, cy - quarter_size); // ne. + activate(cx + half_size, cy, lev); + activate(cx, cy - half_size, lev); + + lev = get_level(cx - quarter_size, cy - quarter_size); // nw. + activate(cx, cy - half_size, lev); + activate(cx - half_size, cy, lev); + + lev = get_level(cx - quarter_size, cy + quarter_size); // sw. + activate(cx - half_size, cy, lev); + activate(cx, cy + half_size, lev); + + lev = get_level(cx + quarter_size, cy + quarter_size); // se. + activate(cx, cy + half_size, lev); + activate(cx + half_size, cy, lev); + } + + // Propagate edge verts to center. + activate(cx, cy, get_level(cx + half_size, cy)); + activate(cx, cy, get_level(cx, cy - half_size)); + activate(cx, cy, get_level(cx, cy + half_size)); + activate(cx, cy, get_level(cx - half_size, cy)); + } + + /// Auxiliary function for generate_block(). + /// Generates a mesh from a triangular quadrant of a square heightfield block. + /// Paraphrased directly out of Lindstrom et al, SIGGRAPH '96. + void generate_quadrant(const heightfield &hf, mesh &mesh, gen_state* state, int lx, int ly, int tx, int ty, int rx, int ry, int recursion_level) const { + if (recursion_level <= 0) return; + + if (hf.get_level(tx, ty) >= state->activation_level) { + // Find base vertex. + int bx = (lx + rx) >> 1; + int by = (ly + ry) >> 1; + + generate_quadrant(hf, mesh, state, lx, ly, bx, by, tx, ty, recursion_level - 1); // left half of quadrant + + if (state->in_my_buffer(tx, ty) == false) { + if ((recursion_level + state->previous_level) & 1) { + state->ptr ^= 1; + } + else { + int x = state->my_buffer[1 - state->ptr][0]; + int y = state->my_buffer[1 - state->ptr][1]; + mesh.emit_vertex(hf, x, y); // or, emit vertex(last - 1); + } + mesh.emit_vertex(hf, tx, ty); + state->set_my_buffer(tx, ty); + state->previous_level = recursion_level; + } + generate_quadrant(hf, mesh, state, tx, ty, bx, by, rx, ry, recursion_level - 1); + } + } + /// Generate the mesh for the specified square with the given center. + /// This is paraphrased directly out of Lindstrom et al, SIGGRAPH '96. + /// It generates a square mesh by walking counterclockwise around four + /// triangular quadrants. + /// The resulting mesh is composed of a single continuous triangle strip, + /// with a few corners turned via degenerate tris where necessary. + void generate_block(const heightfield &hf, mesh &mesh, int activation_level, int log_size, int cx, int cy) const { + int hs = 1 << (log_size - 1); + + // quadrant corner coordinates. + int q[4][2] = { + { cx + hs, cy + hs }, // se + { cx + hs, cy - hs }, // ne + { cx - hs, cy - hs }, // nw + { cx - hs, cy + hs }, // sw + }; + + // Init state for generating mesh. + gen_state state; + state.ptr = 0; + state.previous_level = 0; + state.activation_level = activation_level; + for (int i = 0; i < 4; i++) { + state.my_buffer[i >> 1][i & 1] = -1; + } + + mesh.emit_vertex(hf,q[0][0], q[0][1]); + state.set_my_buffer(q[0][0], q[0][1]); + + {for (int i = 0; i < 4; i++) { + if ((state.previous_level & 1) == 0) { + // tulrich: turn a corner? + state.ptr ^= 1; + } + else { + // tulrich: jump via degenerate? + int x = state.my_buffer[1 - state.ptr][0]; + int y = state.my_buffer[1 - state.ptr][1]; + + mesh.emit_vertex(hf, x, y); // or, emit vertex(last - 1); + } + + // Initial vertex of quadrant. + mesh.emit_vertex(hf,q[i][0], q[i][1]); + state.set_my_buffer(q[i][0], q[i][1]); + state.previous_level = 2 * log_size + 1; + + generate_quadrant(hf, mesh, + &state, + q[i][0], q[i][1], // q[i][l] + cx, cy, // q[i][t] + q[(i + 1) & 3][0], q[(i + 1) & 3][1], // q[i][r] + 2 * log_size + ); + }} + if (state.in_my_buffer(q[0][0], q[0][1]) == false) { + // finish off the strip. @@ may not be necessary? + mesh.emit_vertex(hf, q[0][0], q[0][1]); + } + } +}; + +#endif /* HEIGHTFIELDCHUNKER_HPP */ diff --git a/src/Mesh.hpp b/src/Mesh.hpp new file mode 100644 index 0000000..2e3559a --- /dev/null +++ b/src/Mesh.hpp @@ -0,0 +1,82 @@ +#ifndef CTBMESH_HPP +#define CTBMESH_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file Mesh.hpp + * @brief This declares the `Mesh` class + * @author Alvaro Huarte + */ + +#include +#include +#include +#include "types.hpp" +#include "CTBException.hpp" + +namespace ctb { + class Mesh; +} + +/** + * @brief An abstract base class for a mesh of triangles + */ +class ctb::Mesh +{ +public: + + /// Create an empty mesh + Mesh() + {} + +public: + + /// The array of shared vertices of a mesh + std::vector vertices; + + /// The index collection for each triangle in the mesh (3 for each triangle) + std::vector indices; + + /// Write mesh data to a WKT file + void writeWktFile(const char *fileName) const { + FILE *fp = fopen(fileName, "w"); + + if (fp == NULL) { + throw CTBException("Failed to open file"); + } + + char wktText[512]; + memset(wktText, 0, sizeof(wktText)); + + for (int i = 0, icount = indices.size(); i < icount; i += 3) { + CRSVertex v0 = vertices[indices[i]]; + CRSVertex v1 = vertices[indices[i+1]]; + CRSVertex v2 = vertices[indices[i+2]]; + + sprintf(wktText, "(%.8f %.8f %f, %.8f %.8f %f, %.8f %.8f %f, %.8f %.8f %f)", + v0.x, v0.y, v0.z, + v1.x, v1.y, v1.z, + v2.x, v2.y, v2.z, + v0.x, v0.y, v0.z); + fprintf(fp, "POLYGON Z(%s)\n", wktText); + } + fclose(fp); + }; +}; + +#endif /* CTBMESH_HPP */ diff --git a/src/MeshIterator.hpp b/src/MeshIterator.hpp new file mode 100644 index 0000000..fa083eb --- /dev/null +++ b/src/MeshIterator.hpp @@ -0,0 +1,67 @@ +#ifndef MESHITERATOR_HPP +#define MESHITERATOR_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file MeshIterator.hpp + * @brief This declares the `MeshIterator` class + * @author Alvaro Huarte + */ + +#include "MeshTiler.hpp" +#include "GridIterator.hpp" + +namespace ctb { + class MeshIterator; +} + +/** + * @brief This forward iterates over all `MeshTile`s in a `MeshTiler` + * + * Instances of this class take a `MeshTiler` in the constructor and are used + * to forward iterate over all tiles in the tiler, returning a `MeshTile *` + * when dereferenced. It is the caller's responsibility to call `delete` on the + * tile. + */ +class ctb::MeshIterator : + public GridIterator +{ +public: + + /// Instantiate an iterator with a tiler + MeshIterator(const MeshTiler &tiler) : + MeshIterator(tiler, tiler.maxZoomLevel(), 0) + {} + + MeshIterator(const MeshTiler &tiler, i_zoom startZoom, i_zoom endZoom = 0) : + GridIterator(tiler.grid(), tiler.bounds(), startZoom, endZoom), + tiler(tiler) + {} + + /// Override the dereference operator to return a Tile + virtual MeshTile * + operator*() const { + return tiler.createMesh(*(GridIterator::operator*())); + } + +protected: + + const MeshTiler &tiler; ///< The tiler we are iterating over +}; + +#endif /* MESHITERATOR_HPP */ diff --git a/src/MeshTile.cpp b/src/MeshTile.cpp new file mode 100644 index 0000000..f53d5bc --- /dev/null +++ b/src/MeshTile.cpp @@ -0,0 +1,375 @@ +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file MeshTile.cpp + * @brief This defines the `MeshTile` class + * @author Alvaro Huarte + */ + +#include +#include +#include +#include "zlib.h" + +#include "CTBException.hpp" +#include "MeshTile.hpp" +#include "BoundingSphere.hpp" + +using namespace ctb; + +//////////////////////////////////////////////////////////////////////////////// +// Utility functions + +// Constants taken from http://cesiumjs.org/2013/04/25/Horizon-culling +double llh_ecef_radiusX = 6378137.0; +double llh_ecef_radiusY = 6378137.0; +double llh_ecef_radiusZ = 6356752.3142451793; + +double llh_ecef_rX = 1.0 / llh_ecef_radiusX; +double llh_ecef_rY = 1.0 / llh_ecef_radiusY; +double llh_ecef_rZ = 1.0 / llh_ecef_radiusZ; + +// Stolen from https://github.com/bistromath/gr-air-modes/blob/master/python/mlat.py +// WGS84 reference ellipsoid constants +// http://en.wikipedia.org/wiki/Geodetic_datum#Conversion_calculations +// http://en.wikipedia.org/wiki/File%3aECEF.png +// +double llh_ecef_wgs84_a = llh_ecef_radiusX; // Semi - major axis +double llh_ecef_wgs84_b = llh_ecef_radiusZ; // Semi - minor axis +double llh_ecef_wgs84_e2 = 0.0066943799901975848; // First eccentricity squared + +// LLH2ECEF +static inline double llh_ecef_n(double x) { + double snx = std::sin(x); + return llh_ecef_wgs84_a / std::sqrt(1.0 - llh_ecef_wgs84_e2 * (snx * snx)); +} +static inline CRSVertex LLH2ECEF(const CRSVertex& coordinate) { + double lon = coordinate.x * (M_PI / 180.0); + double lat = coordinate.y * (M_PI / 180.0); + double alt = coordinate.z; + + double x = (llh_ecef_n(lat) + alt) * std::cos(lat) * std::cos(lon); + double y = (llh_ecef_n(lat) + alt) * std::cos(lat) * std::sin(lon); + double z = (llh_ecef_n(lat) * (1.0 - llh_ecef_wgs84_e2) + alt) * std::sin(lat); + + return CRSVertex(x, y, z); +} + +// HORIZON OCCLUSION POINT +// https://cesiumjs.org/2013/05/09/Computing-the-horizon-occlusion-point +// +static inline double ocp_computeMagnitude(const CRSVertex &position, const CRSVertex &sphereCenter) { + double magnitudeSquared = position.magnitudeSquared(); + double magnitude = std::sqrt(magnitudeSquared); + CRSVertex direction = position * (1.0 / magnitude); + + // For the purpose of this computation, points below the ellipsoid + // are considered to be on it instead. + magnitudeSquared = std::fmax(1.0, magnitudeSquared); + magnitude = std::fmax(1.0, magnitude); + + double cosAlpha = direction.dot(sphereCenter); + double sinAlpha = direction.cross(sphereCenter).magnitude(); + double cosBeta = 1.0 / magnitude; + double sinBeta = std::sqrt(magnitudeSquared - 1.0) * cosBeta; + + return 1.0 / (cosAlpha * cosBeta - sinAlpha * sinBeta); +} +static inline CRSVertex ocp_fromPoints(const std::vector &points, const BoundingSphere &boundingSphere) { + const double MIN = -std::numeric_limits::infinity(); + double max_magnitude = MIN; + + // Bring coordinates to ellipsoid scaled coordinates + const CRSVertex ¢er = boundingSphere.center; + CRSVertex scaledCenter = CRSVertex(center.x * llh_ecef_rX, center.y * llh_ecef_rY, center.z * llh_ecef_rZ); + + for (int i = 0, icount = points.size(); i < icount; i++) { + const CRSVertex &point = points[i]; + CRSVertex scaledPoint(point.x * llh_ecef_rX, point.y * llh_ecef_rY, point.z * llh_ecef_rZ); + + double magnitude = ocp_computeMagnitude(scaledPoint, scaledCenter); + if (magnitude > max_magnitude) max_magnitude = magnitude; + } + return scaledCenter * max_magnitude; +} + +// PACKAGE IO +const double SHORT_MAX = 32767.0; +const int BYTESPLIT = 65636; + +static inline int quantizeIndices(const double &origin, const double &factor, const double &value) { + return int(std::round((value - origin) * factor)); +} + +// Write the edge indices of the mesh +template int writeEdgeIndices(gzFile terrainFile, const Mesh &mesh, double edgeCoord, int componentIndex) { + std::vector indices; + std::map ihash; + + for (size_t i = 0, icount = mesh.indices.size(); i < icount; i++) { + uint32_t indice = mesh.indices[i]; + double val = mesh.vertices[indice][componentIndex]; + + if (val == edgeCoord) { + std::map::iterator it = ihash.find(indice); + + if (it == ihash.end()) { + ihash.insert(std::make_pair(indice, i)); + indices.push_back(indice); + } + } + } + + int edgeCount = indices.size(); + gzwrite(terrainFile, &edgeCount, sizeof(int)); + + for (size_t i = 0; i < edgeCount; i++) { + T indice = (T)indices[i]; + gzwrite(terrainFile, &indice, sizeof(T)); + } + return indices.size(); +} + +// ZigZag-Encodes a number (-1 = 1, -2 = 3, 0 = 0, 1 = 2, 2 = 4) +static inline uint16_t zigZagEncode(int n) { + return (n << 1) ^ (n >> 31); +} + +//////////////////////////////////////////////////////////////////////////////// + +MeshTile::MeshTile(): + Tile(), + mChildren(0) +{} + +MeshTile::MeshTile(const TileCoordinate &coord): + Tile(coord), + mChildren(0) +{} + +/** + * @details This writes gzipped terrain data to a file. + */ +void +MeshTile::writeFile(const char *fileName) const { + gzFile terrainFile = gzopen(fileName, "wb"); + + if (terrainFile == NULL) { + throw CTBException("Failed to open file"); + } + + // Calculate main header mesh data + std::vector cartesianVertices; + BoundingSphere cartesianBoundingSphere; + BoundingBox cartesianBounds; + BoundingBox bounds; + + cartesianVertices.resize(mMesh.vertices.size()); + for (size_t i = 0, icount = mMesh.vertices.size(); i < icount; i++) { + const CRSVertex &vertex = mMesh.vertices[i]; + cartesianVertices[i] = LLH2ECEF(vertex); + } + cartesianBoundingSphere.fromPoints(cartesianVertices); + cartesianBounds.fromPoints(cartesianVertices); + bounds.fromPoints(mMesh.vertices); + + + // # Write the mesh header data: + // # https://github.com/AnalyticalGraphicsInc/quantized-mesh + // + // The center of the tile in Earth-centered Fixed coordinates. + double centerX = cartesianBounds.min.x + 0.5 * (cartesianBounds.max.x - cartesianBounds.min.x); + double centerY = cartesianBounds.min.y + 0.5 * (cartesianBounds.max.y - cartesianBounds.min.y); + double centerZ = cartesianBounds.min.z + 0.5 * (cartesianBounds.max.z - cartesianBounds.min.z); + gzwrite(terrainFile, ¢erX, sizeof(double)); + gzwrite(terrainFile, ¢erY, sizeof(double)); + gzwrite(terrainFile, ¢erZ, sizeof(double)); + // + // The minimum and maximum heights in the area covered by this tile. + float minimumHeight = (float)bounds.min.z; + float maximumHeight = (float)bounds.max.z; + gzwrite(terrainFile, &minimumHeight, sizeof(float)); + gzwrite(terrainFile, &maximumHeight, sizeof(float)); + // + // The tile’s bounding sphere. The X,Y,Z coordinates are again expressed + // in Earth-centered Fixed coordinates, and the radius is in meters. + double boundingSphereCenterX = cartesianBoundingSphere.center.x; + double boundingSphereCenterY = cartesianBoundingSphere.center.y; + double boundingSphereCenterZ = cartesianBoundingSphere.center.z; + double boundingSphereRadius = cartesianBoundingSphere.radius; + gzwrite(terrainFile, &boundingSphereCenterX, sizeof(double)); + gzwrite(terrainFile, &boundingSphereCenterY, sizeof(double)); + gzwrite(terrainFile, &boundingSphereCenterZ, sizeof(double)); + gzwrite(terrainFile, &boundingSphereRadius , sizeof(double)); + // + // The horizon occlusion point, expressed in the ellipsoid-scaled Earth-centered Fixed frame. + CRSVertex horizonOcclusionPoint = ocp_fromPoints(cartesianVertices, cartesianBoundingSphere); + gzwrite(terrainFile, &horizonOcclusionPoint.x, sizeof(double)); + gzwrite(terrainFile, &horizonOcclusionPoint.y, sizeof(double)); + gzwrite(terrainFile, &horizonOcclusionPoint.z, sizeof(double)); + + + // # Write mesh vertices (X Y Z components of each vertex): + int vertexCount = mMesh.vertices.size(); + gzwrite(terrainFile, &vertexCount, sizeof(int)); + for (int c = 0; c < 3; c++) { + double origin = bounds.min[c]; + double factor = 0; + if (bounds.max[c] > bounds.min[c]) factor = SHORT_MAX / (bounds.max[c] - bounds.min[c]); + + // Move the initial value + int u0 = quantizeIndices(origin, factor, mMesh.vertices[0][c]), u1, ud; + uint16_t sval = zigZagEncode(u0); + gzwrite(terrainFile, &sval, sizeof(uint16_t)); + + for (size_t i = 1, icount = mMesh.vertices.size(); i < icount; i++) { + u1 = quantizeIndices(origin, factor, mMesh.vertices[i][c]); + ud = u1 - u0; + sval = zigZagEncode(ud); + gzwrite(terrainFile, &sval, sizeof(uint16_t)); + u0 = u1; + } + } + + // # Write mesh indices: + int triangleCount = mMesh.indices.size() / 3; + gzwrite(terrainFile, &triangleCount, sizeof(int)); + if (vertexCount > BYTESPLIT) { + uint32_t highest = 0; + uint32_t code; + + // Write main indices + for (size_t i = 0, icount = mMesh.indices.size(); i < icount; i++) { + code = highest - mMesh.indices[i]; + gzwrite(terrainFile, &code, sizeof(uint32_t)); + if (code == 0) highest++; + } + + // Write all vertices on the edge of the tile (W, S, E, N) + writeEdgeIndices(terrainFile, mMesh, bounds.min.x, 0); + writeEdgeIndices(terrainFile, mMesh, bounds.min.y, 1); + writeEdgeIndices(terrainFile, mMesh, bounds.max.x, 0); + writeEdgeIndices(terrainFile, mMesh, bounds.max.y, 1); + } + else { + uint16_t highest = 0; + uint16_t code; + + // Write main indices + for (size_t i = 0, icount = mMesh.indices.size(); i < icount; i++) { + code = highest - mMesh.indices[i]; + gzwrite(terrainFile, &code, sizeof(uint16_t)); + if (code == 0) highest++; + } + + // Write all vertices on the edge of the tile (W, S, E, N) + writeEdgeIndices(terrainFile, mMesh, bounds.min.x, 0); + writeEdgeIndices(terrainFile, mMesh, bounds.min.y, 1); + writeEdgeIndices(terrainFile, mMesh, bounds.max.x, 0); + writeEdgeIndices(terrainFile, mMesh, bounds.max.y, 1); + } + + // Try and close the file + switch (gzclose(terrainFile)) { + case Z_OK: + break; + case Z_STREAM_ERROR: + case Z_ERRNO: + case Z_MEM_ERROR: + case Z_BUF_ERROR: + default: + throw CTBException("Failed to close file"); + } +} + +bool +MeshTile::hasChildren() const { + return mChildren; +} + +bool +MeshTile::hasChildSW() const { + return ((mChildren & TERRAIN_CHILD_SW) == TERRAIN_CHILD_SW); +} + +bool +MeshTile::hasChildSE() const { + return ((mChildren & TERRAIN_CHILD_SE) == TERRAIN_CHILD_SE); +} + +bool +MeshTile::hasChildNW() const { + return ((mChildren & TERRAIN_CHILD_NW) == TERRAIN_CHILD_NW); +} + +bool +MeshTile::hasChildNE() const { + return ((mChildren & TERRAIN_CHILD_NE) == TERRAIN_CHILD_NE); +} + +void +MeshTile::setChildSW(bool on) { + if (on) { + mChildren |= TERRAIN_CHILD_SW; + } else { + mChildren &= ~TERRAIN_CHILD_SW; + } +} + +void +MeshTile::setChildSE(bool on) { + if (on) { + mChildren |= TERRAIN_CHILD_SE; + } else { + mChildren &= ~TERRAIN_CHILD_SE; + } +} + +void +MeshTile::setChildNW(bool on) { + if (on) { + mChildren |= TERRAIN_CHILD_NW; + } else { + mChildren &= ~TERRAIN_CHILD_NW; + } +} + +void +MeshTile::setChildNE(bool on) { + if (on) { + mChildren |= TERRAIN_CHILD_NE; + } else { + mChildren &= ~TERRAIN_CHILD_NE; + } +} + +void +MeshTile::setAllChildren(bool on) { + if (on) { + mChildren = TERRAIN_CHILD_SW | TERRAIN_CHILD_SE | TERRAIN_CHILD_NW | TERRAIN_CHILD_NE; + } else { + mChildren = 0; + } +} + +const Mesh & MeshTile::getMesh() const { + return mMesh; +} + +Mesh & MeshTile::getMesh() { + return mMesh; +} diff --git a/src/MeshTile.hpp b/src/MeshTile.hpp new file mode 100644 index 0000000..df7ae76 --- /dev/null +++ b/src/MeshTile.hpp @@ -0,0 +1,127 @@ +#ifndef MESHTILE_HPP +#define MESHTILE_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file MeshTile.hpp + * @brief This declares the `MeshTile` class + * @author Alvaro Huarte + */ + +#include "config.hpp" +#include "Mesh.hpp" +#include "TileCoordinate.hpp" +#include "Tile.hpp" + +namespace ctb { + class MeshTile; +} + +/** + * @brief `Terrain` data associated with a `Mesh` + * + * This aims to implement the Cesium [quantized-mesh-1.0 terrain + * format](https://github.com/AnalyticalGraphicsInc/quantized-mesh). + */ +class CTB_DLL ctb::MeshTile : + public Tile +{ + friend class MeshTiler; + +public: + + /// Create an empty mesh tile object + MeshTile(); + + /// Create a mesh tile from a tile coordinate + MeshTile(const TileCoordinate &coord); + + /// Write terrain data to the filesystem + void + writeFile(const char *fileName) const; + + /// Does the terrain tile have child tiles? + bool + hasChildren() const; + + /// Does the terrain tile have a south west child tile? + bool + hasChildSW() const; + + /// Does the terrain tile have a south east child tile? + bool + hasChildSE() const; + + /// Does the terrain tile have a north west child tile? + bool + hasChildNW() const; + + /// Does the terrain tile have a north east child tile? + bool + hasChildNE() const; + + /// Specify that there is a south west child tile + void + setChildSW(bool on = true); + + /// Specify that there is a south east child tile + void + setChildSE(bool on = true); + + /// Specify that there is a north west child tile + void + setChildNW(bool on = true); + + /// Specify that there is a north east child tile + void + setChildNE(bool on = true); + + /// Specify that all child tiles are present + void + setAllChildren(bool on = true); + + /// Get the mesh data as a const object + const ctb::Mesh & getMesh() const; + + /// Get the mesh data + ctb::Mesh & getMesh(); + +protected: + + /// The terrain mesh data + ctb::Mesh mMesh; + +private: + + char mChildren; ///< The child flags + + /** + * @brief Bit flags defining child tile existence + * + * There is a good discussion on bitflags + * [here](http://www.dylanleigh.net/notes/c-cpp-tricks.html#Using_"Bitflags"). + */ + enum Children { + TERRAIN_CHILD_SW = 1, // 2^0, bit 0 + TERRAIN_CHILD_SE = 2, // 2^1, bit 1 + TERRAIN_CHILD_NW = 4, // 2^2, bit 2 + TERRAIN_CHILD_NE = 8 // 2^3, bit 3 + }; +}; + +#endif /* MESHTILE_HPP */ diff --git a/src/MeshTiler.cpp b/src/MeshTiler.cpp new file mode 100644 index 0000000..b07cd2f --- /dev/null +++ b/src/MeshTiler.cpp @@ -0,0 +1,208 @@ +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file MeshTiler.cpp + * @brief This defines the `MeshTiler` class + * @author Alvaro Huarte + */ + +#include "CTBException.hpp" +#include "MeshTiler.hpp" +#include "HeightFieldChunker.hpp" + +using namespace ctb; + +//////////////////////////////////////////////////////////////////////////////// + +/** + * Implementation of ctb::chunk::mesh for ctb::Mesh class. + */ +class WrapperMesh : public ctb::chunk::mesh { +private: + CRSBounds &mBounds; + Mesh &mMesh; + double mCellSizeX; + double mCellSizeY; + + std::map mIndicesMap; + Coordinate mTriangles[3]; + bool mTriOddOrder; + int mTriIndex; + +public: + WrapperMesh(CRSBounds &bounds, Mesh &mesh, i_tile tileSizeX, i_tile tileSizeY): + mMesh(mesh), + mBounds(bounds), + mTriOddOrder(false), + mTriIndex(0) { + mCellSizeX = (bounds.getMaxX() - bounds.getMinX()) / (double)(tileSizeX - 1); + mCellSizeY = (bounds.getMaxY() - bounds.getMinY()) / (double)(tileSizeY - 1); + } + + virtual void clear() { + mMesh.vertices.clear(); + mMesh.indices.clear(); + mIndicesMap.clear(); + mTriOddOrder = false; + mTriIndex = 0; + } + virtual void emit_vertex(const ctb::chunk::heightfield &heightfield, int x, int y) { + mTriangles[mTriIndex].x = x; + mTriangles[mTriIndex].y = y; + mTriIndex++; + + if (mTriIndex == 3) { + mTriOddOrder = !mTriOddOrder; + + if (mTriOddOrder) { + appendVertex(heightfield, mTriangles[0].x, mTriangles[0].y); + appendVertex(heightfield, mTriangles[1].x, mTriangles[1].y); + appendVertex(heightfield, mTriangles[2].x, mTriangles[2].y); + } + else { + appendVertex(heightfield, mTriangles[1].x, mTriangles[1].y); + appendVertex(heightfield, mTriangles[0].x, mTriangles[0].y); + appendVertex(heightfield, mTriangles[2].x, mTriangles[2].y); + } + mTriangles[0].x = mTriangles[1].x; + mTriangles[0].y = mTriangles[1].y; + mTriangles[1].x = mTriangles[2].x; + mTriangles[1].y = mTriangles[2].y; + mTriIndex--; + } + } + void appendVertex(const ctb::chunk::heightfield &heightfield, int x, int y) { + int iv; + int index = heightfield.indexOfGridCoordinate(x, y); + + std::map::iterator it = mIndicesMap.find(index); + + if (it == mIndicesMap.end()) { + iv = mMesh.vertices.size(); + + double xmin = mBounds.getMinX(); + double ymax = mBounds.getMaxY(); + double height = heightfield.height(x, y); + + mMesh.vertices.push_back(CRSVertex(xmin + (x * mCellSizeX), ymax - (y * mCellSizeY), height)); + mIndicesMap.insert(std::make_pair(index, iv)); + } + else { + iv = it->second; + } + mMesh.indices.push_back(iv); + } +}; + +//////////////////////////////////////////////////////////////////////////////// + +MeshTile * +ctb::MeshTiler::createMesh(const TileCoordinate &coord) const { + // Get a mesh tile represented by the tile coordinate + MeshTile *terrainTile = new MeshTile(coord); + GDALTile *rasterTile = createRasterTile(coord); // the raster associated with this tile coordinate + GDALRasterBand *heightsBand = rasterTile->dataset->GetRasterBand(1); + + // Get the initial width and height of tile data in a tile + const ctb::i_tile TILE_SIZE = mGrid.tileSize(); + const ctb::i_tile TILE_CELL_SIZE = TILE_SIZE * TILE_SIZE; + + // Copy the raster data into an array + float *rasterHeights = (float *)CPLCalloc(TILE_CELL_SIZE, sizeof(float)); + if (heightsBand->RasterIO(GF_Read, 0, 0, TILE_SIZE, TILE_SIZE, + (void *) rasterHeights, TILE_SIZE, TILE_SIZE, GDT_Float32, + 0, 0) != CE_None) { + CPLFree(rasterHeights); + throw CTBException("Could not read heights from raster"); + } + + delete rasterTile; + + // Number of tiles in the horizontal direction at tile level zero. + double resolutionAtLevelZero = mGrid.resolution(0); + int numberOfTilesAtLevelZero = (int)(mGrid.getExtent().getWidth() / (TILE_SIZE * resolutionAtLevelZero)); + // Default quality of terrain created from heightmaps (TerrainProvider.js). + double heightmapTerrainQuality = 0.25; + // Earth semi-major-axis in meters. + const double semiMajorAxis = 6378137.0; + // Appropriate geometric error estimate when the geometry comes from a heightmap (TerrainProvider.js). + double maximumGeometricError = MeshTiler::getEstimatedLevelZeroGeometricErrorForAHeightmap( + semiMajorAxis, + heightmapTerrainQuality * mMeshQualityFactor, + TILE_SIZE, + numberOfTilesAtLevelZero + ); + // Geometric error for current Level. + maximumGeometricError /= (double)(1 << coord.zoom); + + // Convert the raster grid into an irregular mesh applying the Chunked LOD strategy by 'Thatcher Ulrich'. + // http://tulrich.com/geekstuff/chunklod.html + // + ctb::chunk::heightfield heightfield(rasterHeights, TILE_SIZE); + heightfield.applyGeometricError(maximumGeometricError, coord.zoom <= 6); + // + ctb::CRSBounds mGridBounds = mGrid.tileBounds(coord); + Mesh &tileMesh = terrainTile->getMesh(); + WrapperMesh mesh(mGridBounds, tileMesh, TILE_SIZE, TILE_SIZE); + heightfield.generateMesh(mesh, 0); + heightfield.clear(); + + CPLFree(rasterHeights); + + // If we are not at the maximum zoom level we need to set child flags on the + // tile where child tiles overlap the dataset bounds. + if (coord.zoom != maxZoomLevel()) { + CRSBounds tileBounds = mGrid.tileBounds(coord); + + if (! (bounds().overlaps(tileBounds))) { + terrainTile->setAllChildren(false); + } else { + if (bounds().overlaps(tileBounds.getSW())) { + terrainTile->setChildSW(); + } + if (bounds().overlaps(tileBounds.getNW())) { + terrainTile->setChildNW(); + } + if (bounds().overlaps(tileBounds.getNE())) { + terrainTile->setChildNE(); + } + if (bounds().overlaps(tileBounds.getSE())) { + terrainTile->setChildSE(); + } + } + } + + return terrainTile; +} + +MeshTiler & +ctb::MeshTiler::operator=(const MeshTiler &other) { + TerrainTiler::operator=(other); + + return *this; +} + +double ctb::MeshTiler::getEstimatedLevelZeroGeometricErrorForAHeightmap( + double maximumRadius, + double heightmapTerrainQuality, + int tileWidth, + int numberOfTilesAtLevelZero) +{ + double error = maximumRadius * 2 * M_PI * heightmapTerrainQuality; + error /= (double)(tileWidth * numberOfTilesAtLevelZero); + return error; +} diff --git a/src/MeshTiler.hpp b/src/MeshTiler.hpp new file mode 100644 index 0000000..575f19b --- /dev/null +++ b/src/MeshTiler.hpp @@ -0,0 +1,80 @@ +#ifndef MESHTILER_HPP +#define MESHTILER_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file MeshTiler.hpp + * @brief This declares the `MeshTiler` class + * @author Alvaro Huarte + */ + +#include "MeshTile.hpp" +#include "TerrainTiler.hpp" + +namespace ctb { + class MeshTiler; +} + +/** + * @brief Create `MeshTile`s from a GDAL Dataset + * + * This class derives from `GDALTiler` and `TerrainTiler` enabling `MeshTile`s + * to be created for a specific `TileCoordinate`. + */ +class CTB_DLL ctb::MeshTiler : + public TerrainTiler +{ +public: + + /// Instantiate a tiler with all required arguments + MeshTiler(GDALDataset *poDataset, const Grid &grid, const TilerOptions &options, double meshQualityFactor = 1.0): + TerrainTiler(poDataset, grid, options), + mMeshQualityFactor(meshQualityFactor) {} + + /// Instantiate a tiler with an empty GDAL dataset + MeshTiler(double meshQualityFactor = 1.0): + TerrainTiler(), + mMeshQualityFactor(meshQualityFactor) {} + + /// Instantiate a tiler with a dataset and grid but no options + MeshTiler(GDALDataset *poDataset, const Grid &grid, double meshQualityFactor = 1.0): + TerrainTiler(poDataset, grid, TilerOptions()), + mMeshQualityFactor(meshQualityFactor) {} + + /// Overload the assignment operator + MeshTiler & + operator=(const MeshTiler &other); + + /// Create a mesh from a tile coordinate + MeshTile * + createMesh(const TileCoordinate &coord) const; + +protected: + + // Specifies the factor of the quality to convert terrain heightmaps to meshes. + double mMeshQualityFactor; + + // Determines an appropriate geometric error estimate when the geometry comes from a heightmap. + static double getEstimatedLevelZeroGeometricErrorForAHeightmap( + double maximumRadius, + double heightmapTerrainQuality, + int tileWidth, + int numberOfTilesAtLevelZero); +}; + +#endif /* MESHTILER_HPP */ diff --git a/src/types.hpp b/src/types.hpp index 9fc93f3..7c61ad2 100644 --- a/src/types.hpp +++ b/src/types.hpp @@ -25,6 +25,7 @@ #include // uint16_t #include "Bounds.hpp" +#include "Coordinate3D.hpp" /// All terrain related data types reside in this namespace namespace ctb { @@ -38,7 +39,8 @@ namespace ctb { // Complex types typedef Bounds TileBounds; ///< Tile extents in tile coordinates typedef Coordinate PixelPoint; ///< The location of a pixel - typedef Coordinate CRSPoint; ///< A Coordinate Reference System coordinate + typedef Coordinate CRSPoint; ///< A Coordinate Reference System coordinate + typedef Coordinate3D CRSVertex; ///< A 3D-Vertex of a mesh or tile in CRS coordinates typedef Bounds CRSBounds; ///< Extents in CRS coordinates typedef Coordinate TilePoint; ///< The location of a tile diff --git a/tools/ctb-tile.cpp b/tools/ctb-tile.cpp index 3b4b679..1788bbf 100644 --- a/tools/ctb-tile.cpp +++ b/tools/ctb-tile.cpp @@ -50,6 +50,7 @@ #include "GlobalMercator.hpp" #include "RasterIterator.hpp" #include "TerrainIterator.hpp" +#include "MeshIterator.hpp" using namespace std; using namespace ctb; @@ -73,7 +74,10 @@ class TerrainBuild : public Command { startZoom(-1), endZoom(-1), verbosity(1), - resume(false) + resume(false), + meshQualityFactor(1.0), + metadata(false), + cesiumFriendly(false) {} void @@ -198,6 +202,21 @@ class TerrainBuild : public Command { return (command->argc == 1) ? command->argv[0] : NULL; } + static void + setMeshQualityFactor(command_t *command) { + static_cast(Command::self(command))->meshQualityFactor = atof(command->arg); + } + + static void + setMetadata(command_t *command) { + static_cast(Command::self(command))->metadata = true; + } + + static void + setCesiumFriendly(command_t *command) { + static_cast(Command::self(command))->cesiumFriendly = true; + } + const char *outputDir, *outputFormat, *profile; @@ -212,6 +231,10 @@ class TerrainBuild : public Command { CPLStringList creationOptions; TilerOptions tilerOptions; + + double meshQualityFactor; + bool metadata; + bool cesiumFriendly; }; /** @@ -339,9 +362,178 @@ fileExists(const std::string& filename) { return VSIStatExL(filename.c_str(), &statbuf, VSI_STAT_EXISTS_FLAG) == 0; } +static bool +fileCopy(const std::string& sourceFilename, const std::string& targetFilename) { + FILE *fp_s = VSIFOpen(sourceFilename.c_str(), "rb"); + if (!fp_s) return false; + FILE *fp_t = VSIFOpen(targetFilename.c_str(), "wb"); + if (!fp_t) return false; + + VSIFSeek(fp_s, 0, SEEK_END); + long fileSize = VSIFTell(fp_s); + + if (fileSize > 0) + { + VSIFSeek(fp_s, 0, SEEK_SET); + + void* buffer = VSIMalloc(fileSize); + VSIFRead(buffer, 1, fileSize, fp_s); + VSIFWrite(buffer, 1, fileSize, fp_t); + VSIFree(buffer); + } + VSIFClose(fp_t); + VSIFClose(fp_s); + + return fileSize > 0; +} + +/// Handle the terrain metadata +class TerrainMetadata { +public: + TerrainMetadata() { + } + + // Defines the valid tile indexes of a level in a Tileset + struct LevelInfo { + public: + LevelInfo() { + startX = startY = std::numeric_limits::max(); + finalX = finalY = std::numeric_limits::min(); + } + int startX, startY; + int finalX, finalY; + + inline void add(const TileCoordinate *coordinate) { + startX = std::min(startX, (int)coordinate->x); + startY = std::min(startY, (int)coordinate->y); + finalX = std::max(finalX, (int)coordinate->x); + finalY = std::max(finalY, (int)coordinate->y); + } + inline void add(const LevelInfo &level) { + startX = std::min(startX, level.startX); + startY = std::min(startY, level.startY); + finalX = std::max(finalX, level.finalX); + finalY = std::max(finalY, level.finalY); + } + }; + std::vector levels; + + // Defines the bounding box covered by the Terrain + CRSBounds bounds; + + // Add metadata of the specified Coordinate + void add(const Grid &grid, const TileCoordinate *coordinate) { + CRSBounds tileBounds = grid.tileBounds(*coordinate); + i_zoom zoom = coordinate->zoom; + + if ((1 + zoom) > levels.size()) { + for (size_t i = 0; i <= zoom; i++) { + levels.push_back(LevelInfo()); + } + } + LevelInfo &level = levels[zoom]; + level.add(coordinate); + + if (bounds.getMaxX() == bounds.getMinX()) { + bounds = tileBounds; + } + else { + bounds.setMinX(std::min(bounds.getMinX(), tileBounds.getMinX())); + bounds.setMinY(std::min(bounds.getMinY(), tileBounds.getMinY())); + bounds.setMaxX(std::max(bounds.getMaxX(), tileBounds.getMaxX())); + bounds.setMaxY(std::max(bounds.getMaxY(), tileBounds.getMaxY())); + } + } + + // Add metadata info + void add(const TerrainMetadata &otherMetadata) { + if (otherMetadata.levels.size() > 0) { + const CRSBounds &otherBounds = otherMetadata.bounds; + + for (size_t i = 0, icount = (otherMetadata.levels.size() - levels.size()); i < icount; i++) { + levels.push_back(LevelInfo()); + } + for (size_t i = 0; i < levels.size(); i++) { + levels[i].add(otherMetadata.levels[i]); + } + + bounds.setMinX(std::min(bounds.getMinX(), otherBounds.getMinX())); + bounds.setMinY(std::min(bounds.getMinY(), otherBounds.getMinY())); + bounds.setMaxX(std::max(bounds.getMaxX(), otherBounds.getMaxX())); + bounds.setMaxY(std::max(bounds.getMaxY(), otherBounds.getMaxY())); + } + } + + /// Output the layer.json metadata file + /// http://help.agi.com/TerrainServer/RESTAPIGuide.html + /// Example: + /// https://assets.agi.com/stk-terrain/v1/tilesets/world/tiles/layer.json + void writeJsonFile(const std::string &filename, const std::string &datasetName, const std::string &outputFormat = "Terrain", const std::string &profile = "geodetic") const { + FILE *fp = fopen(filename.c_str(), "w"); + + if (fp == NULL) { + throw CTBException("Failed to open metadata file"); + } + + fprintf(fp, "{\n"); + fprintf(fp, " \"tilejson\": \"2.1.0\",\n"); + fprintf(fp, " \"name\": \"%s\",\n", datasetName.c_str()); + fprintf(fp, " \"description\": \"\",\n"); + fprintf(fp, " \"version\": \"1.1.0\",\n"); + + if (strcmp(outputFormat.c_str(), "Terrain") == 0) { + fprintf(fp, " \"format\": \"heightmap-1.0\",\n"); + } + else if (strcmp(outputFormat.c_str(), "Mesh") == 0) { + fprintf(fp, " \"format\": \"quantized-mesh-1.0\",\n"); + } + else { + fprintf(fp, " \"format\": \"GDAL\",\n"); + } + fprintf(fp, " \"attribution\": \"\",\n"); + fprintf(fp, " \"schema\": \"tms\",\n"); + fprintf(fp, " \"tiles\": [ \"{z}/{x}/{y}.terrain?v={version}\" ],\n"); + + if (strcmp(profile.c_str(), "geodetic") == 0) { + fprintf(fp, " \"projection\": \"EPSG:4326\",\n"); + } + else { + fprintf(fp, " \"projection\": \"EPSG:3857\",\n"); + } + fprintf(fp, " \"bounds\": [ %.2f, %.2f, %.2f, %.2f ],\n", + bounds.getMinX(), + bounds.getMinY(), + bounds.getMaxX(), + bounds.getMaxY()); + + fprintf(fp, " \"available\": [\n"); + for (size_t i = 0, icount = levels.size(); i < icount; i++) { + const LevelInfo &level = levels[i]; + + if (i > 0) + fprintf(fp, " ,[ "); + else + fprintf(fp, " [ "); + + if (level.finalX >= level.startX) { + fprintf(fp, "{ \"startX\": %li, \"startY\": %li, \"endX\": %li, \"endY\": %li }", + level.startX, + level.startY, + level.finalX, + level.finalY); + } + fprintf(fp, " ]\n"); + } + fprintf(fp, " ]\n"); + + fprintf(fp, "}\n"); + fclose(fp); + } +}; + /// Output GDAL tiles represented by a tiler to a directory static void -buildGDAL(const RasterTiler &tiler, TerrainBuild *command) { +buildGDAL(const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(command->outputFormat); if (poDriver == NULL) { @@ -365,8 +557,8 @@ buildGDAL(const RasterTiler &tiler, TerrainBuild *command) { const TileCoordinate *coordinate = iter.GridIterator::operator*(); GDALDataset *poDstDS; const string filename = getTileFilename(coordinate, dirname, extension); + if (metadata) metadata->add(tiler.grid(), coordinate); - if( !command->resume || !fileExists(filename) ) { GDALTile *tile = *iter; const string temp_filename = concat(filename, ".tmp"); @@ -393,7 +585,7 @@ buildGDAL(const RasterTiler &tiler, TerrainBuild *command) { /// Output terrain tiles represented by a tiler to a directory static void -buildTerrain(const TerrainTiler &tiler, TerrainBuild *command) { +buildTerrain(const TerrainTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { const string dirname = string(command->outputDir) + osDirSep; i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom, endZoom = (command->endZoom < 0) ? 0 : command->endZoom; @@ -405,6 +597,7 @@ buildTerrain(const TerrainTiler &tiler, TerrainBuild *command) { while (!iter.exhausted()) { const TileCoordinate *coordinate = iter.GridIterator::operator*(); const string filename = getTileFilename(coordinate, dirname, "terrain"); + if (metadata) metadata->add(tiler.grid(), coordinate); if( !command->resume || !fileExists(filename) ) { TerrainTile *tile = *iter; @@ -423,26 +616,110 @@ buildTerrain(const TerrainTiler &tiler, TerrainBuild *command) { } } +/// Output mesh tiles represented by a tiler to a directory +static void +buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { + const string dirname = string(command->outputDir) + osDirSep; + i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom, + endZoom = (command->endZoom < 0) ? 0 : command->endZoom; + + // DEBUG Chunker: + #if 0 + TileCoordinate coordinate(13, 8102, 6047); + MeshTile *tile = tiler.createMesh(coordinate); + // + const string txtname = getTileFilename(&coordinate, dirname, "wkt"); + const Mesh &mesh = tile->getMesh(); + mesh.writeWktFile(txtname.c_str()); + // + CRSBounds bounds = tiler.grid().tileBounds(coordinate); + double x = bounds.getMinX() + 0.5 * (bounds.getMaxX() - bounds.getMinX()); + double y = bounds.getMinY() + 0.5 * (bounds.getMaxY() - bounds.getMinY()); + CRSPoint point(x,y); + TileCoordinate c = tiler.grid().crsToTile(point, coordinate.zoom); + // + const string filename = getTileFilename(&coordinate, dirname, "terrain"); + tile->writeFile(filename.c_str()); + delete tile; + return; + #endif + + MeshIterator iter(tiler, startZoom, endZoom); + int currentIndex = incrementIterator(iter, 0); + setIteratorSize(iter); + + while (!iter.exhausted()) { + const TileCoordinate *coordinate = iter.GridIterator::operator*(); + const string filename = getTileFilename(coordinate, dirname, "terrain"); + if (metadata) metadata->add(tiler.grid(), coordinate); + + if( !command->resume || !fileExists(filename) ) { + MeshTile *tile = *iter; + const string temp_filename = concat(filename, ".tmp"); + + tile->writeFile(temp_filename.c_str()); + delete tile; + + if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) { + throw new CTBException("Could not rename temporary file"); + } + } + + currentIndex = incrementIterator(iter, currentIndex); + showProgress(currentIndex, filename); + } +} + +static void +buildMetadata(const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { + const string dirname = string(command->outputDir) + osDirSep; + i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom, + endZoom = (command->endZoom < 0) ? 0 : command->endZoom; + + const std::string filename = concat(dirname, "layer.json"); + + RasterIterator iter(tiler, startZoom, endZoom); + int currentIndex = incrementIterator(iter, 0); + setIteratorSize(iter); + + while (!iter.exhausted()) { + const TileCoordinate *coordinate = iter.GridIterator::operator*(); + if (metadata) metadata->add(tiler.grid(), coordinate); + + currentIndex = incrementIterator(iter, currentIndex); + showProgress(currentIndex, filename); + } +} + /** * Perform a tile building operation * * This function is designed to be run in a separate thread. */ static int -runTiler(TerrainBuild *command, Grid *grid) { +runTiler(TerrainBuild *command, Grid *grid, TerrainMetadata *metadata) { GDALDataset *poDataset = (GDALDataset *) GDALOpen(command->getInputFilename(), GA_ReadOnly); if (poDataset == NULL) { cerr << "Error: could not open GDAL dataset" << endl; return 1; } + // Metadata of only this thread, it will be joined to global later + TerrainMetadata *threadMetadata = metadata ? new TerrainMetadata() : NULL; + try { - if (strcmp(command->outputFormat, "Terrain") == 0) { + if (command->metadata) { + const RasterTiler tiler(poDataset, *grid, command->tilerOptions); + buildMetadata(tiler, command, threadMetadata); + } else if (strcmp(command->outputFormat, "Terrain") == 0) { const TerrainTiler tiler(poDataset, *grid); - buildTerrain(tiler, command); + buildTerrain(tiler, command, threadMetadata); + } else if (strcmp(command->outputFormat, "Mesh") == 0) { + const MeshTiler tiler(poDataset, *grid, command->tilerOptions, command->meshQualityFactor); + buildMesh(tiler, command, threadMetadata); } else { // it's a GDAL format const RasterTiler tiler(poDataset, *grid, command->tilerOptions); - buildGDAL(tiler, command); + buildGDAL(tiler, command, threadMetadata); } } catch (CTBException &e) { @@ -451,6 +728,14 @@ runTiler(TerrainBuild *command, Grid *grid) { GDALClose(poDataset); + // Pass metadata to global instance. + if (threadMetadata) { + static std::mutex mutex; + std::lock_guard lock(mutex); + + metadata->add(*threadMetadata); + delete threadMetadata; + } return 0; } @@ -460,7 +745,7 @@ main(int argc, char *argv[]) { TerrainBuild command = TerrainBuild(argv[0], version.cstr); command.setUsage("[options] GDAL_DATASOURCE"); command.option("-o", "--output-dir ", "specify the output directory for the tiles (defaults to working directory)", TerrainBuild::setOutputDir); - command.option("-f", "--output-format ", "specify the output format for the tiles. This is either `Terrain` (the default) or any format listed by `gdalinfo --formats`", TerrainBuild::setOutputFormat); + command.option("-f", "--output-format ", "specify the output format for the tiles. This is either `Terrain` (the default), `Mesh` (Chunked LOD mesh), or any format listed by `gdalinfo --formats`", TerrainBuild::setOutputFormat); command.option("-p", "--profile ", "specify the TMS profile for the tiles. This is either `geodetic` (the default) or `mercator`", TerrainBuild::setProfile); command.option("-c", "--thread-count ", "specify the number of threads to use for tile generation. On multicore machines this defaults to the number of CPUs", TerrainBuild::setThreadCount); command.option("-t", "--tile-size ", "specify the size of the tiles in pixels. This defaults to 65 for terrain tiles and 256 for other GDAL formats", TerrainBuild::setTileSize); @@ -471,6 +756,9 @@ main(int argc, char *argv[]) { command.option("-z", "--error-threshold ", "specify the error threshold in pixel units for transformation approximation. Larger values should mean faster transforms. Defaults to 0.125", TerrainBuild::setErrorThreshold); command.option("-m", "--warp-memory ", "The memory limit in bytes used for warp operations. Higher settings should be faster. Defaults to a conservative GDAL internal setting.", TerrainBuild::setWarpMemory); command.option("-R", "--resume", "Do not overwrite existing files", TerrainBuild::setResume); + command.option("-g", "--mesh-qfactor ", "specify the factor to multiply the estimated geometric error to convert heightmaps to irregular meshes. Larger values should mean minor quality. Defaults to 1.0", TerrainBuild::setMeshQualityFactor); + command.option("-l", "--layer", "only output the layer.json metadata file", TerrainBuild::setMetadata); + command.option("-C", "--cesium-friendly", "Force the creation of missing root tiles to be CesiumJS-friendly", TerrainBuild::setCesiumFriendly); command.option("-q", "--quiet", "only output errors", TerrainBuild::setQuiet); command.option("-v", "--verbose", "be more noisy", TerrainBuild::setVerbose); @@ -514,11 +802,16 @@ main(int argc, char *argv[]) { vector> tasks; int threadCount = (command.threadCount > 0) ? command.threadCount : CPLGetNumCPUs(); + // Calculate metadata? + const string dirname = string(command.outputDir) + osDirSep; + const std::string filename = concat(dirname, "layer.json"); + TerrainMetadata *metadata = command.metadata ? new TerrainMetadata() : NULL; + // Instantiate the threads using futures from a packaged_task for (int i = 0; i < threadCount ; ++i) { - packaged_task task(runTiler); // wrap the function + packaged_task task(runTiler); // wrap the function tasks.push_back(task.get_future()); // get a future - thread(move(task), &command, &grid).detach(); // launch on a thread + thread(move(task), &command, &grid, metadata).detach(); // launch on a thread } // Synchronise the completion of the threads @@ -531,8 +824,39 @@ main(int argc, char *argv[]) { int retval = task.get(); // return on the first encountered problem - if (retval) + if (retval) { + delete metadata; return retval; + } + } + + // Write Json metadata file? + if (metadata) { + std::string datasetName(command.getInputFilename()); + datasetName = datasetName.substr(datasetName.find_last_of("/\\") + 1); + const size_t rfindpos = datasetName.rfind('.'); + if (std::string::npos != rfindpos) datasetName = datasetName.erase(rfindpos); + + metadata->writeJsonFile(filename, datasetName, std::string(command.outputFormat), std::string(command.profile)); + delete metadata; + } + + // CesiumJS friendly? + if (command.cesiumFriendly && (strcmp(command.profile, "geodetic") == 0) && command.endZoom <= 0) { + std::string dirName0 = string(command.outputDir) + osDirSep + "0" + osDirSep + "0"; + std::string dirName1 = string(command.outputDir) + osDirSep + "0" + osDirSep + "1"; + std::string tileName0 = dirName0 + osDirSep + "0.terrain"; + std::string tileName1 = dirName1 + osDirSep + "0.terrain"; + + if (fileExists(tileName0) && !fileExists(tileName1)) { + VSIMkdir(dirName1.c_str(), 0755); + fileCopy(tileName0, tileName1); + } + else + if (!fileExists(tileName0) && fileExists(tileName1)) { + VSIMkdir(dirName0.c_str(), 0755); + fileCopy(tileName1, tileName0); + } } return 0; From 3a804a9aff88363c6cb1db29d52d34769c0d3273 Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Wed, 28 Mar 2018 01:41:28 +0200 Subject: [PATCH 02/26] Fix nodata values --- src/GDALTiler.cpp | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/GDALTiler.cpp b/src/GDALTiler.cpp index 70fc735..e4b414d 100644 --- a/src/GDALTiler.cpp +++ b/src/GDALTiler.cpp @@ -246,7 +246,11 @@ getOverviewDataset(GDALDatasetH hSrcDS, GDALTransformerFunc pfnTransformer, void if( iOvr >= 0 ) { //std::cout << "CTB WARPING: Selecting overview level " << iOvr << " for output dataset " << nPixels << "x" << nLines << std::endl; - poSrcOvrDS = GDALCreateOverviewDataset( poSrcDS, iOvr, FALSE); + #if ( GDAL_VERSION_MAJOR >= 2 && GDAL_VERSION_MINOR >= 2 ) + poSrcOvrDS = GDALCreateOverviewDataset( poSrcDS, iOvr, FALSE ); + #else + poSrcOvrDS = GDALCreateOverviewDataset( poSrcDS, iOvr, FALSE, FALSE ); + #endif } } } @@ -304,7 +308,25 @@ GDALTiler::createRasterTile(double (&adfGeoTransform)[6]) const { psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount ); + psWarpOptions->padfSrcNoDataReal = + (double *)CPLCalloc(psWarpOptions->nBandCount, sizeof(double)); + psWarpOptions->padfSrcNoDataImag = + (double *)CPLCalloc(psWarpOptions->nBandCount, sizeof(double)); + psWarpOptions->padfDstNoDataReal = + (double *)CPLCalloc(psWarpOptions->nBandCount, sizeof(double)); + psWarpOptions->padfDstNoDataImag = + (double *)CPLCalloc(psWarpOptions->nBandCount, sizeof(double)); + for (short unsigned int i = 0; i < psWarpOptions->nBandCount; ++i) { + int bGotNoData = FALSE; + double noDataValue = dataset()->GetRasterBand(i + 1)->GetNoDataValue(&bGotNoData); + if (!bGotNoData) noDataValue = -32768; + + psWarpOptions->padfSrcNoDataReal[i] = noDataValue; + psWarpOptions->padfSrcNoDataImag[i] = 0; + psWarpOptions->padfDstNoDataReal[i] = noDataValue; + psWarpOptions->padfDstNoDataImag[i] = 0; + psWarpOptions->panDstBands[i] = psWarpOptions->panSrcBands[i] = i + 1; } From 358888cf2365f8e441addff9a2b27b994dba14df Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Sat, 7 Apr 2018 09:23:29 +0200 Subject: [PATCH 03/26] Fix 'Integer Overflow' error when creating tiles of low levels --- src/CMakeLists.txt | 2 + src/GDALDatasetReader.cpp | 156 ++++++++++++++++++++++++++++++++++++++ src/GDALDatasetReader.hpp | 100 ++++++++++++++++++++++++ src/GDALTile.cpp | 15 ++++ src/GDALTile.hpp | 3 + src/GDALTiler.cpp | 12 +-- src/GDALTiler.hpp | 9 ++- src/MeshIterator.hpp | 7 +- src/MeshTiler.cpp | 55 ++++++++------ src/MeshTiler.hpp | 9 ++- src/RasterTiler.hpp | 4 +- src/TerrainIterator.hpp | 4 + src/TerrainTiler.cpp | 52 ++++++++----- src/TerrainTiler.hpp | 11 ++- src/TilerIterator.hpp | 2 +- tools/ctb-tile.cpp | 7 +- 16 files changed, 385 insertions(+), 63 deletions(-) create mode 100644 src/GDALDatasetReader.cpp create mode 100644 src/GDALDatasetReader.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index dd18c3b..49d39d0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -20,6 +20,7 @@ include_directories(${ZLIB_INCLUDE_DIRS}) add_library(ctb SHARED GDALTile.cpp GDALTiler.cpp + GDALDatasetReader.cpp TerrainTiler.cpp TerrainTile.cpp MeshTiler.cpp @@ -36,6 +37,7 @@ set(HEADERS Coordinate3D.hpp GDALTile.hpp GDALTiler.hpp + GDALDatasetReader.hpp GlobalGeodetic.hpp GlobalMercator.hpp Grid.hpp diff --git a/src/GDALDatasetReader.cpp b/src/GDALDatasetReader.cpp new file mode 100644 index 0000000..43f97ef --- /dev/null +++ b/src/GDALDatasetReader.cpp @@ -0,0 +1,156 @@ +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file GDALDatasetReader.cpp + * @brief This defines the `GDALDatasetReader` class + */ + +#include "gdal_priv.h" +#include "gdalwarper.h" + +#include "CTBException.hpp" +#include "GDALDatasetReader.hpp" +#include "TerrainTiler.hpp" + +using namespace ctb; + +/** + * @details + * Read a region of raster heights for the specified Dataset and Coordinate. + * This method uses `GDALRasterBand::RasterIO` function. + */ +float * +ctb::GDALDatasetReader::readRasterHeights(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) { + GDALTile *rasterTile = createRasterTile(tiler, dataset, coord); // the raster associated with this tile coordinate + + const ctb::i_tile TILE_CELL_SIZE = tileSizeX * tileSizeY; + float *rasterHeights = (float *)CPLCalloc(TILE_CELL_SIZE, sizeof(float)); + + GDALRasterBand *heightsBand = rasterTile->dataset->GetRasterBand(1); + + if (heightsBand->RasterIO(GF_Read, 0, 0, tileSizeX, tileSizeY, + (void *) rasterHeights, tileSizeX, tileSizeY, GDT_Float32, + 0, 0) != CE_None) { + delete rasterTile; + CPLFree(rasterHeights); + + throw CTBException("Could not read heights from raster"); + } + delete rasterTile; + return rasterHeights; +} + +/// Create a raster tile from a tile coordinate +GDALTile * +ctb::GDALDatasetReader::createRasterTile(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord) { + return tiler.createRasterTile(dataset, coord); +} + +/// Create a VTR raster overview from a GDALDataset +GDALDataset * +ctb::GDALDatasetReader::createOverview(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord, int overviewIndex) { + int nFactorScale = 2 << overviewIndex; + int nRasterXSize = dataset->GetRasterXSize() / nFactorScale; + int nRasterYSize = dataset->GetRasterYSize() / nFactorScale; + + GDALDataset *poOverview = NULL; + double adfGeoTransform[6]; + + // Should We create an overview of the Dataset? + if (nRasterXSize > 4 && nRasterYSize > 4 && dataset->GetGeoTransform(adfGeoTransform) == CE_None) { + adfGeoTransform[1] *= nFactorScale; + adfGeoTransform[5] *= nFactorScale; + + TerrainTiler tempTiler(tiler.dataset(), tiler.grid(), tiler.options); + tempTiler.crsWKT = ""; + GDALTile *rasterTile = createRasterTile(tempTiler, dataset, coord); + if (rasterTile) { + poOverview = rasterTile->detach(); + delete rasterTile; + } + } + return poOverview; +} + +/// The destructor +ctb::GDALDatasetReaderWithOverviews::~GDALDatasetReaderWithOverviews() { + reset(); +} + +/// Read a region of raster heights into an array for the specified Dataset and Coordinate +float * +ctb::GDALDatasetReaderWithOverviews::readRasterHeights(GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) { + GDALDataset *mainDataset = dataset; + + const ctb::i_tile TILE_CELL_SIZE = tileSizeX * tileSizeY; + float *rasterHeights = (float *)CPLCalloc(TILE_CELL_SIZE, sizeof(float)); + + // Replace GDAL Dataset by last valid Overview. + for (int i = mOverviews.size() - 1; i >= 0; --i) { + if (mOverviews[i]) { + dataset = mOverviews[i]; + break; + } + } + + // Extract the raster data, using overviews when necessary + bool rasterOk = false; + while (!rasterOk) { + GDALTile *rasterTile = createRasterTile(poTiler, dataset, coord); // the raster associated with this tile coordinate + + GDALRasterBand *heightsBand = rasterTile->dataset->GetRasterBand(1); + + if (heightsBand->RasterIO(GF_Read, 0, 0, tileSizeX, tileSizeY, + (void *) rasterHeights, tileSizeX, tileSizeY, GDT_Float32, + 0, 0) != CE_None) { + + GDALDataset *psOverview = createOverview(poTiler, mainDataset, coord, mOverviewIndex++); + if (psOverview) { + mOverviews.push_back(psOverview); + dataset = psOverview; + } + else { + delete rasterTile; + CPLFree(rasterHeights); + throw CTBException("Could not create an overview of current GDAL dataset"); + } + } + else { + rasterOk = true; + } + delete rasterTile; + } + + // Everything ok? + if (!rasterOk) { + CPLFree(rasterHeights); + throw CTBException("Could not read heights from raster"); + } + return rasterHeights; +} + +/// Releases all overviews +void +ctb::GDALDatasetReaderWithOverviews::reset() { + mOverviewIndex = 0; + + for (int i = mOverviews.size() - 1; i >= 0; --i) { + GDALDataset *poOverview = mOverviews[i]; + GDALClose(poOverview); + } + mOverviews.clear(); +} diff --git a/src/GDALDatasetReader.hpp b/src/GDALDatasetReader.hpp new file mode 100644 index 0000000..a3ef18c --- /dev/null +++ b/src/GDALDatasetReader.hpp @@ -0,0 +1,100 @@ +#ifndef GDALDATASETREADER_HPP +#define GDALDATASETREADER_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file GDALDatasetReader.hpp + * @brief This declares the `GDALDatasetReader` class + */ + +#include +#include +#include "gdalwarper.h" + +#include "TileCoordinate.hpp" +#include "GDALTiler.hpp" + +namespace ctb { + class GDALDatasetReader; + class GDALDatasetReaderWithOverviews; +} + +/** + * @brief Read raster tiles from a GDAL Dataset + * + * This abstract base class is associated with a GDAL dataset. + * It allows to read a region of the raster according to + * a region defined by a Tile Coordinate. + * + * We can define our own + */ +class CTB_DLL ctb::GDALDatasetReader { +public: + /// Read a region of raster heights into an array for the specified Dataset and Coordinate + static float * + readRasterHeights(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY); + + /// Read a region of raster heights into an array for the specified Dataset and Coordinate + virtual float * + readRasterHeights(GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) = 0; + +protected: + /// Create a raster tile from a tile coordinate + static GDALTile * + createRasterTile(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord); + + /// Create a VTR raster overview from a GDALDataset + static GDALDataset * + createOverview(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord, int overviewIndex); +}; + +/** + * @brief Implements a GDALDatasetReader that takes care of 'Integer overflow' errors. + * + * This class creates Overviews to avoid 'Integer overflow' errors when extracting + * raster data. + */ +class CTB_DLL ctb::GDALDatasetReaderWithOverviews : public ctb::GDALDatasetReader { +public: + + /// Instantiate a GDALDatasetReaderWithOverviews + GDALDatasetReaderWithOverviews(const GDALTiler &tiler): + poTiler(tiler), + mOverviewIndex(0) {} + + /// The destructor + ~GDALDatasetReaderWithOverviews(); + + /// Read a region of raster heights into an array for the specified Dataset and Coordinate + virtual float * + readRasterHeights(GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) override; + + /// Releases all overviews + void reset(); + +protected: + /// The tiler to use + const GDALTiler &poTiler; + + /// List of VRT Overviews of the underlying GDAL dataset + std::vector mOverviews; + /// Current VRT Overview + int mOverviewIndex; +}; + +#endif /* GDALDATASETREADER_HPP */ diff --git a/src/GDALTile.cpp b/src/GDALTile.cpp index ad7df74..c3a6dae 100644 --- a/src/GDALTile.cpp +++ b/src/GDALTile.cpp @@ -34,3 +34,18 @@ GDALTile::~GDALTile() { } } } + +/// Detach the underlying GDAL dataset +GDALDataset *GDALTile::detach() { + if (dataset != NULL) { + GDALDataset *poDataset = dataset; + dataset = NULL; + + if (transformer != NULL) { + GDALDestroyGenImgProjTransformer(transformer); + transformer = NULL; + } + return poDataset; + } + return NULL; +} diff --git a/src/GDALTile.hpp b/src/GDALTile.hpp index 86ba270..85914b7 100644 --- a/src/GDALTile.hpp +++ b/src/GDALTile.hpp @@ -57,6 +57,9 @@ class CTB_DLL ctb::GDALTile : GDALDataset *dataset; + /// Detach the underlying GDAL dataset + GDALDataset *detach(); + protected: friend class GDALTiler; diff --git a/src/GDALTiler.cpp b/src/GDALTiler.cpp index e4b414d..fda83a7 100644 --- a/src/GDALTiler.cpp +++ b/src/GDALTiler.cpp @@ -173,7 +173,7 @@ GDALTiler::~GDALTiler() { } GDALTile * -GDALTiler::createRasterTile(const TileCoordinate &coord) const { +GDALTiler::createRasterTile(GDALDataset *dataset, const TileCoordinate &coord) const { // Convert the tile bounds into a geo transform double adfGeoTransform[6], resolution = mGrid.resolution(coord.zoom); @@ -186,7 +186,7 @@ GDALTiler::createRasterTile(const TileCoordinate &coord) const { adfGeoTransform[4] = 0; adfGeoTransform[5] = -resolution; - GDALTile *tile = createRasterTile(adfGeoTransform); + GDALTile *tile = createRasterTile(dataset, adfGeoTransform); static_cast(*tile) = coord; // Set the shifted geo transform to the VRT @@ -271,13 +271,13 @@ getOverviewDataset(GDALDatasetH hSrcDS, GDALTransformerFunc pfnTransformer, void * dataset. */ GDALTile * -GDALTiler::createRasterTile(double (&adfGeoTransform)[6]) const { - if (poDataset == NULL) { +GDALTiler::createRasterTile(GDALDataset *dataset, double (&adfGeoTransform)[6]) const { + if (dataset == NULL) { throw CTBException("No GDAL dataset is set"); } // The source and sink datasets - GDALDatasetH hSrcDS = (GDALDatasetH) dataset(); + GDALDatasetH hSrcDS = (GDALDatasetH) dataset; GDALDatasetH hDstDS; // The transformation option list @@ -319,7 +319,7 @@ GDALTiler::createRasterTile(double (&adfGeoTransform)[6]) const { for (short unsigned int i = 0; i < psWarpOptions->nBandCount; ++i) { int bGotNoData = FALSE; - double noDataValue = dataset()->GetRasterBand(i + 1)->GetNoDataValue(&bGotNoData); + double noDataValue = poDataset->GetRasterBand(i + 1)->GetNoDataValue(&bGotNoData); if (!bGotNoData) noDataValue = -32768; psWarpOptions->padfSrcNoDataReal[i] = noDataValue; diff --git a/src/GDALTiler.hpp b/src/GDALTiler.hpp index 1f3a3ee..87477db 100644 --- a/src/GDALTiler.hpp +++ b/src/GDALTiler.hpp @@ -33,6 +33,7 @@ namespace ctb { struct TilerOptions; class GDALTiler; + class GDALDatasetReader; // forward declaration } /// Options passed to a `GDALTiler` @@ -91,7 +92,7 @@ class CTB_DLL ctb::GDALTiler { /// Create a tile from a tile coordinate virtual Tile * - createTile(const TileCoordinate &coord) const = 0; + createTile(GDALDataset *dataset, const TileCoordinate &coord) const = 0; /// Get the maximum zoom level for the dataset inline i_zoom @@ -151,16 +152,18 @@ class CTB_DLL ctb::GDALTiler { } protected: + friend class GDALDatasetReader; + /// Close the underlying dataset void closeDataset(); /// Create a raster tile from a tile coordinate virtual GDALTile * - createRasterTile(const TileCoordinate &coord) const; + createRasterTile(GDALDataset *dataset, const TileCoordinate &coord) const; /// Create a raster tile from a geo transform virtual GDALTile * - createRasterTile(double (&adfGeoTransform)[6]) const; + createRasterTile(GDALDataset *dataset, double (&adfGeoTransform)[6]) const; /// The grid used for generating tiles Grid mGrid; diff --git a/src/MeshIterator.hpp b/src/MeshIterator.hpp index fa083eb..3383b67 100644 --- a/src/MeshIterator.hpp +++ b/src/MeshIterator.hpp @@ -56,7 +56,12 @@ class ctb::MeshIterator : /// Override the dereference operator to return a Tile virtual MeshTile * operator*() const { - return tiler.createMesh(*(GridIterator::operator*())); + return tiler.createMesh(tiler.dataset(), *(GridIterator::operator*())); + } + + virtual MeshTile * + operator*(ctb::GDALDatasetReader *reader) const { + return tiler.createMesh(tiler.dataset(), *(GridIterator::operator*()), reader); } protected: diff --git a/src/MeshTiler.cpp b/src/MeshTiler.cpp index b07cd2f..185d4bf 100644 --- a/src/MeshTiler.cpp +++ b/src/MeshTiler.cpp @@ -23,6 +23,7 @@ #include "CTBException.hpp" #include "MeshTiler.hpp" #include "HeightFieldChunker.hpp" +#include "GDALDatasetReader.hpp" using namespace ctb; @@ -110,31 +111,13 @@ class WrapperMesh : public ctb::chunk::mesh { //////////////////////////////////////////////////////////////////////////////// -MeshTile * -ctb::MeshTiler::createMesh(const TileCoordinate &coord) const { - // Get a mesh tile represented by the tile coordinate - MeshTile *terrainTile = new MeshTile(coord); - GDALTile *rasterTile = createRasterTile(coord); // the raster associated with this tile coordinate - GDALRasterBand *heightsBand = rasterTile->dataset->GetRasterBand(1); - - // Get the initial width and height of tile data in a tile - const ctb::i_tile TILE_SIZE = mGrid.tileSize(); - const ctb::i_tile TILE_CELL_SIZE = TILE_SIZE * TILE_SIZE; - - // Copy the raster data into an array - float *rasterHeights = (float *)CPLCalloc(TILE_CELL_SIZE, sizeof(float)); - if (heightsBand->RasterIO(GF_Read, 0, 0, TILE_SIZE, TILE_SIZE, - (void *) rasterHeights, TILE_SIZE, TILE_SIZE, GDT_Float32, - 0, 0) != CE_None) { - CPLFree(rasterHeights); - throw CTBException("Could not read heights from raster"); - } - - delete rasterTile; +void +ctb::MeshTiler::prepareSettingsOfTile(MeshTile *terrainTile, const TileCoordinate &coord, float *rasterHeights, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) const { + const ctb::i_tile TILE_SIZE = tileSizeX; // Number of tiles in the horizontal direction at tile level zero. double resolutionAtLevelZero = mGrid.resolution(0); - int numberOfTilesAtLevelZero = (int)(mGrid.getExtent().getWidth() / (TILE_SIZE * resolutionAtLevelZero)); + int numberOfTilesAtLevelZero = (int)(mGrid.getExtent().getWidth() / (tileSizeX * resolutionAtLevelZero)); // Default quality of terrain created from heightmaps (TerrainProvider.js). double heightmapTerrainQuality = 0.25; // Earth semi-major-axis in meters. @@ -157,12 +140,10 @@ ctb::MeshTiler::createMesh(const TileCoordinate &coord) const { // ctb::CRSBounds mGridBounds = mGrid.tileBounds(coord); Mesh &tileMesh = terrainTile->getMesh(); - WrapperMesh mesh(mGridBounds, tileMesh, TILE_SIZE, TILE_SIZE); + WrapperMesh mesh(mGridBounds, tileMesh, tileSizeX, tileSizeY); heightfield.generateMesh(mesh, 0); heightfield.clear(); - CPLFree(rasterHeights); - // If we are not at the maximum zoom level we need to set child flags on the // tile where child tiles overlap the dataset bounds. if (coord.zoom != maxZoomLevel()) { @@ -185,6 +166,30 @@ ctb::MeshTiler::createMesh(const TileCoordinate &coord) const { } } } +} + +MeshTile * +ctb::MeshTiler::createMesh(GDALDataset *dataset, const TileCoordinate &coord) const { + // Copy the raster data into an array + float *rasterHeights = ctb::GDALDatasetReader::readRasterHeights(*this, dataset, coord, mGrid.tileSize(), mGrid.tileSize()); + + // Get a mesh tile represented by the tile coordinate + MeshTile *terrainTile = new MeshTile(coord); + prepareSettingsOfTile(terrainTile, coord, rasterHeights, mGrid.tileSize(), mGrid.tileSize()); + CPLFree(rasterHeights); + + return terrainTile; +} + +MeshTile * +ctb::MeshTiler::createMesh(GDALDataset *dataset, const TileCoordinate &coord, ctb::GDALDatasetReader *reader) const { + // Copy the raster data into an array + float *rasterHeights = reader->readRasterHeights(dataset, coord, mGrid.tileSize(), mGrid.tileSize()); + + // Get a mesh tile represented by the tile coordinate + MeshTile *terrainTile = new MeshTile(coord); + prepareSettingsOfTile(terrainTile, coord, rasterHeights, mGrid.tileSize(), mGrid.tileSize()); + CPLFree(rasterHeights); return terrainTile; } diff --git a/src/MeshTiler.hpp b/src/MeshTiler.hpp index 575f19b..b2f6991 100644 --- a/src/MeshTiler.hpp +++ b/src/MeshTiler.hpp @@ -62,7 +62,11 @@ class CTB_DLL ctb::MeshTiler : /// Create a mesh from a tile coordinate MeshTile * - createMesh(const TileCoordinate &coord) const; + createMesh(GDALDataset *dataset, const TileCoordinate &coord) const; + + /// Create a mesh from a tile coordinate + MeshTile * + createMesh(GDALDataset *dataset, const TileCoordinate &coord, GDALDatasetReader *reader) const; protected: @@ -75,6 +79,9 @@ class CTB_DLL ctb::MeshTiler : double heightmapTerrainQuality, int tileWidth, int numberOfTilesAtLevelZero); + + /// Assigns settings of Tile just to use. + void prepareSettingsOfTile(MeshTile *tile, const TileCoordinate &coord, float *rasterHeights, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) const; }; #endif /* MESHTILER_HPP */ diff --git a/src/RasterTiler.hpp b/src/RasterTiler.hpp index f673e33..2bab229 100644 --- a/src/RasterTiler.hpp +++ b/src/RasterTiler.hpp @@ -54,8 +54,8 @@ class ctb::RasterTiler : /// Override to return a covariant data type virtual GDALTile * - createTile(const TileCoordinate &coord) const override { - return createRasterTile(coord); + createTile(GDALDataset *dataset, const TileCoordinate &coord) const override { + return createRasterTile(dataset, coord); } }; diff --git a/src/TerrainIterator.hpp b/src/TerrainIterator.hpp index f009e00..8d8396b 100644 --- a/src/TerrainIterator.hpp +++ b/src/TerrainIterator.hpp @@ -56,6 +56,10 @@ class ctb::TerrainIterator : operator*() const override { return static_cast(TilerIterator::operator*()); } + virtual TerrainTile * + operator*(ctb::GDALDatasetReader *reader) const { + return (static_cast(tiler)).createTile(tiler.dataset(), *(GridIterator::operator*()), reader); + } }; #endif /* TERRAINITERATOR_HPP */ diff --git a/src/TerrainTiler.cpp b/src/TerrainTiler.cpp index 174abf0..77463aa 100644 --- a/src/TerrainTiler.cpp +++ b/src/TerrainTiler.cpp @@ -21,32 +21,20 @@ #include "CTBException.hpp" #include "TerrainTiler.hpp" +#include "GDALDatasetReader.hpp" using namespace ctb; -TerrainTile * -ctb::TerrainTiler::createTile(const TileCoordinate &coord) const { - // Get a terrain tile represented by the tile coordinate - TerrainTile *terrainTile = new TerrainTile(coord); - GDALTile *rasterTile = createRasterTile(coord); // the raster associated with this tile coordinate - GDALRasterBand *heightsBand = rasterTile->dataset->GetRasterBand(1); - - // Copy the raster data into an array - float rasterHeights[TerrainTile::TILE_CELL_SIZE]; - if (heightsBand->RasterIO(GF_Read, 0, 0, TILE_SIZE, TILE_SIZE, - (void *) rasterHeights, TILE_SIZE, TILE_SIZE, GDT_Float32, - 0, 0) != CE_None) { - throw CTBException("Could not read heights from raster"); - } - - delete rasterTile; +void +ctb::TerrainTiler::prepareSettingsOfTile(TerrainTile *terrainTile, const TileCoordinate &coord, float *rasterHeights, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) const { + const ctb::i_tile TILE_CELL_SIZE = tileSizeX * tileSizeY; // Convert the raster data into the terrain tile heights. This assumes the // input raster data represents meters above sea level. Each terrain height // value is the number of 1/5 meter units above -1000 meters. // TODO: try doing this using a VRT derived band: // (http://www.gdal.org/gdal_vrttut.html) - for (unsigned short int i = 0; i < TerrainTile::TILE_CELL_SIZE; i++) { + for (unsigned short int i = 0; i < TILE_CELL_SIZE; i++) { terrainTile->mHeights[i] = (i_terrain_height) ((rasterHeights[i] + 1000) * 5); } @@ -72,14 +60,38 @@ ctb::TerrainTiler::createTile(const TileCoordinate &coord) const { } } } +} + +TerrainTile * +ctb::TerrainTiler::createTile(GDALDataset *dataset, const TileCoordinate &coord) const { + // Copy the raster data into an array + float *rasterHeights = ctb::GDALDatasetReader::readRasterHeights(*this, dataset, coord, TILE_SIZE, TILE_SIZE); + + // Get a terrain tile represented by the tile coordinate + TerrainTile *terrainTile = new TerrainTile(coord); + prepareSettingsOfTile(terrainTile, coord, rasterHeights, TILE_SIZE, TILE_SIZE); + CPLFree(rasterHeights); + + return terrainTile; +} + +TerrainTile * +ctb::TerrainTiler::createTile(GDALDataset *dataset, const TileCoordinate &coord, ctb::GDALDatasetReader *reader) const { + // Copy the raster data into an array + float *rasterHeights = reader->readRasterHeights(dataset, coord, TILE_SIZE, TILE_SIZE); + + // Get a mesh tile represented by the tile coordinate + TerrainTile *terrainTile = new TerrainTile(coord); + prepareSettingsOfTile(terrainTile, coord, rasterHeights, TILE_SIZE, TILE_SIZE); + CPLFree(rasterHeights); return terrainTile; } GDALTile * -ctb::TerrainTiler::createRasterTile(const TileCoordinate &coord) const { +ctb::TerrainTiler::createRasterTile(GDALDataset *dataset, const TileCoordinate &coord) const { // Ensure we have some data from which to create a tile - if (poDataset && poDataset->GetRasterCount() < 1) { + if (dataset && dataset->GetRasterCount() < 1) { throw CTBException("At least one band must be present in the GDAL dataset"); } @@ -97,7 +109,7 @@ ctb::TerrainTiler::createRasterTile(const TileCoordinate &coord) const { adfGeoTransform[4] = 0; adfGeoTransform[5] = -resolution; - GDALTile *tile = GDALTiler::createRasterTile(adfGeoTransform); + GDALTile *tile = GDALTiler::createRasterTile(dataset, adfGeoTransform); // The previous geotransform represented the data with an overlap as required // by the terrain specification. This now needs to be overwritten so that diff --git a/src/TerrainTiler.hpp b/src/TerrainTiler.hpp index 00da46b..171ae3e 100644 --- a/src/TerrainTiler.hpp +++ b/src/TerrainTiler.hpp @@ -59,13 +59,17 @@ class CTB_DLL ctb::TerrainTiler : /// Override to return a covariant data type TerrainTile * - createTile(const TileCoordinate &coord) const override; + createTile(GDALDataset *dataset, const TileCoordinate &coord) const override; + + /// Create a tile from a tile coordinate + TerrainTile * + createTile(GDALDataset *dataset, const TileCoordinate &coord, GDALDatasetReader *reader) const; protected: /// Create a `GDALTile` representing the required terrain tile data virtual GDALTile * - createRasterTile(const TileCoordinate &coord) const override; + createRasterTile(GDALDataset *dataset, const TileCoordinate &coord) const override; /** * @brief Get terrain bounds shifted to introduce a pixel overlap @@ -97,6 +101,9 @@ class CTB_DLL ctb::TerrainTiler : return tile; } + + /// Assigns settings of Tile just to use. + void prepareSettingsOfTile(TerrainTile *tile, const TileCoordinate &coord, float *rasterHeights, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) const; }; #endif /* TERRAINTILER_HPP */ diff --git a/src/TilerIterator.hpp b/src/TilerIterator.hpp index fca6f3a..83ed2d9 100644 --- a/src/TilerIterator.hpp +++ b/src/TilerIterator.hpp @@ -55,7 +55,7 @@ class ctb::TilerIterator : /// Override the dereference operator to return a Tile virtual Tile * operator*() const { - return tiler.createTile(*(GridIterator::operator*())); + return tiler.createTile(tiler.dataset(), *(GridIterator::operator*())); } protected: diff --git a/tools/ctb-tile.cpp b/tools/ctb-tile.cpp index 1788bbf..8187f04 100644 --- a/tools/ctb-tile.cpp +++ b/tools/ctb-tile.cpp @@ -51,6 +51,7 @@ #include "RasterIterator.hpp" #include "TerrainIterator.hpp" #include "MeshIterator.hpp" +#include "GDALDatasetReader.hpp" using namespace std; using namespace ctb; @@ -593,6 +594,7 @@ buildTerrain(const TerrainTiler &tiler, TerrainBuild *command, TerrainMetadata * TerrainIterator iter(tiler, startZoom, endZoom); int currentIndex = incrementIterator(iter, 0); setIteratorSize(iter); + GDALDatasetReaderWithOverviews reader(tiler); while (!iter.exhausted()) { const TileCoordinate *coordinate = iter.GridIterator::operator*(); @@ -600,7 +602,7 @@ buildTerrain(const TerrainTiler &tiler, TerrainBuild *command, TerrainMetadata * if (metadata) metadata->add(tiler.grid(), coordinate); if( !command->resume || !fileExists(filename) ) { - TerrainTile *tile = *iter; + TerrainTile *tile = iter.operator*(&reader); const string temp_filename = concat(filename, ".tmp"); tile->writeFile(temp_filename.c_str()); @@ -647,6 +649,7 @@ buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metada MeshIterator iter(tiler, startZoom, endZoom); int currentIndex = incrementIterator(iter, 0); setIteratorSize(iter); + GDALDatasetReaderWithOverviews reader(tiler); while (!iter.exhausted()) { const TileCoordinate *coordinate = iter.GridIterator::operator*(); @@ -654,7 +657,7 @@ buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metada if (metadata) metadata->add(tiler.grid(), coordinate); if( !command->resume || !fileExists(filename) ) { - MeshTile *tile = *iter; + MeshTile *tile = iter.operator*(&reader); const string temp_filename = concat(filename, ".tmp"); tile->writeFile(temp_filename.c_str()); From c448bc557b1dd0367e2e23d7fcca106381c4690a Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Wed, 11 Apr 2018 09:29:00 +0200 Subject: [PATCH 04/26] Option to write 'Oct-Encoded Per-Vertex Normals' for Terrain Lighting --- src/MeshTile.cpp | 82 +++++++++++++++++++++++++++++++++++++++++++++- src/MeshTile.hpp | 2 +- tools/ctb-tile.cpp | 20 +++++++---- 3 files changed, 96 insertions(+), 8 deletions(-) diff --git a/src/MeshTile.cpp b/src/MeshTile.cpp index f53d5bc..11ad15f 100644 --- a/src/MeshTile.cpp +++ b/src/MeshTile.cpp @@ -23,6 +23,7 @@ #include #include #include +#include "cpl_conv.h" #include "zlib.h" #include "CTBException.hpp" @@ -149,6 +150,46 @@ static inline uint16_t zigZagEncode(int n) { return (n << 1) ^ (n >> 31); } +// Triangle area +static inline double triangleArea(const CRSVertex &a, const CRSVertex &b) { + double i = std::pow(a[1] * b[2] - a[2] * b[1], 2); + double j = std::pow(a[2] * b[0] - a[0] * b[2], 2); + double k = std::pow(a[0] * b[1] - a[1] * b[0], 2); + return 0.5 * sqrt(i + j + k); +} + +// Constraint a value to lie between two values +static inline double clamp_value(double value, double min, double max) { + return value < min ? min : value > max ? max : value; +} + +// Converts a scalar value in the range [-1.0, 1.0] to a SNORM in the range [0, rangeMax] +static inline unsigned char snorm_value(double value, double rangeMax = 255) { + return (unsigned char)int(std::round((clamp_value(value, -1.0, 1.0) * 0.5 + 0.5) * rangeMax)); +} + +/** + * Encodes a normalized vector into 2 SNORM values in the range of [0-rangeMax] following the 'oct' encoding. + * + * Oct encoding is a compact representation of unit length vectors. + * The 'oct' encoding is described in "A Survey of Efficient Representations of Independent Unit Vectors", + * Cigolle et al 2014: {@link http://jcgt.org/published/0003/02/01/} + */ +static inline Coordinate octEncode(const CRSVertex &vector, double rangeMax = 255) { + Coordinate temp; + double llnorm = std::abs(vector.x) + std::abs(vector.y) + std::abs(vector.z); + temp.x = vector.x / llnorm; + temp.y = vector.y / llnorm; + + if (vector.z < 0) { + double x = temp.x; + double y = temp.y; + temp.x = (1.0 - std::abs(y)) * (x < 0.0 ? -1.0 : 1.0); + temp.y = (1.0 - std::abs(x)) * (y < 0.0 ? -1.0 : 1.0); + } + return Coordinate(snorm_value(temp.x, rangeMax), snorm_value(temp.y, rangeMax)); +} + //////////////////////////////////////////////////////////////////////////////// MeshTile::MeshTile(): @@ -165,7 +206,7 @@ MeshTile::MeshTile(const TileCoordinate &coord): * @details This writes gzipped terrain data to a file. */ void -MeshTile::writeFile(const char *fileName) const { +MeshTile::writeFile(const char *fileName, bool writeVertexNormals) const { gzFile terrainFile = gzopen(fileName, "wb"); if (terrainFile == NULL) { @@ -283,6 +324,45 @@ MeshTile::writeFile(const char *fileName) const { writeEdgeIndices(terrainFile, mMesh, bounds.max.y, 1); } + // # Write 'Oct-Encoded Per-Vertex Normals' for Terrain Lighting: + if (writeVertexNormals && triangleCount > 0) { + unsigned char extensionId = 1; + gzwrite(terrainFile, &extensionId, sizeof(unsigned char)); + int extensionLength = 2 * vertexCount; + gzwrite(terrainFile, &extensionLength, sizeof(int)); + + std::vector normalsPerVertex(vertexCount); + std::vector normalsPerFace(triangleCount); + std::vector areasPerFace(triangleCount); + + for (size_t i = 0, icount = mMesh.indices.size(), j = 0; i < icount; i+=3, j++) { + const CRSVertex &v0 = cartesianVertices[ mMesh.indices[i ] ]; + const CRSVertex &v1 = cartesianVertices[ mMesh.indices[i+1] ]; + const CRSVertex &v2 = cartesianVertices[ mMesh.indices[i+2] ]; + + CRSVertex normal = (v1 - v0).cross(v2 - v0); + double area = triangleArea(v0, v1); + normalsPerFace[j] = normal; + areasPerFace[j] = area; + } + for (size_t i = 0, icount = mMesh.indices.size(), j = 0; i < icount; i+=3, j++) { + int indexV0 = mMesh.indices[i ]; + int indexV1 = mMesh.indices[i+1]; + int indexV2 = mMesh.indices[i+2]; + + CRSVertex weightedNormal = normalsPerFace[j] * areasPerFace[j]; + + normalsPerVertex[indexV0] = normalsPerVertex[indexV0] + weightedNormal; + normalsPerVertex[indexV1] = normalsPerVertex[indexV1] + weightedNormal; + normalsPerVertex[indexV2] = normalsPerVertex[indexV2] + weightedNormal; + } + for (size_t i = 0; i < vertexCount; i++) { + Coordinate xy = octEncode(normalsPerVertex[i].normalize()); + gzwrite(terrainFile, &xy.x, sizeof(unsigned char)); + gzwrite(terrainFile, &xy.y, sizeof(unsigned char)); + } + } + // Try and close the file switch (gzclose(terrainFile)) { case Z_OK: diff --git a/src/MeshTile.hpp b/src/MeshTile.hpp index df7ae76..5cd3fc9 100644 --- a/src/MeshTile.hpp +++ b/src/MeshTile.hpp @@ -53,7 +53,7 @@ class CTB_DLL ctb::MeshTile : /// Write terrain data to the filesystem void - writeFile(const char *fileName) const; + writeFile(const char *fileName, bool writeVertexNormals = false) const; /// Does the terrain tile have child tiles? bool diff --git a/tools/ctb-tile.cpp b/tools/ctb-tile.cpp index 8187f04..3851ba0 100644 --- a/tools/ctb-tile.cpp +++ b/tools/ctb-tile.cpp @@ -78,7 +78,8 @@ class TerrainBuild : public Command { resume(false), meshQualityFactor(1.0), metadata(false), - cesiumFriendly(false) + cesiumFriendly(false), + vertexNormals(false) {} void @@ -218,6 +219,11 @@ class TerrainBuild : public Command { static_cast(Command::self(command))->cesiumFriendly = true; } + static void + setVertexNormals(command_t *command) { + static_cast(Command::self(command))->vertexNormals = true; + } + const char *outputDir, *outputFormat, *profile; @@ -236,6 +242,7 @@ class TerrainBuild : public Command { double meshQualityFactor; bool metadata; bool cesiumFriendly; + bool vertexNormals; }; /** @@ -620,7 +627,7 @@ buildTerrain(const TerrainTiler &tiler, TerrainBuild *command, TerrainMetadata * /// Output mesh tiles represented by a tiler to a directory static void -buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { +buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata, bool writeVertexNormals = false) { const string dirname = string(command->outputDir) + osDirSep; i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom, endZoom = (command->endZoom < 0) ? 0 : command->endZoom; @@ -628,7 +635,7 @@ buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metada // DEBUG Chunker: #if 0 TileCoordinate coordinate(13, 8102, 6047); - MeshTile *tile = tiler.createMesh(coordinate); + MeshTile *tile = tiler.createMesh(tiler.dataset(), coordinate); // const string txtname = getTileFilename(&coordinate, dirname, "wkt"); const Mesh &mesh = tile->getMesh(); @@ -641,7 +648,7 @@ buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metada TileCoordinate c = tiler.grid().crsToTile(point, coordinate.zoom); // const string filename = getTileFilename(&coordinate, dirname, "terrain"); - tile->writeFile(filename.c_str()); + tile->writeFile(filename.c_str(), writeVertexNormals); delete tile; return; #endif @@ -660,7 +667,7 @@ buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metada MeshTile *tile = iter.operator*(&reader); const string temp_filename = concat(filename, ".tmp"); - tile->writeFile(temp_filename.c_str()); + tile->writeFile(temp_filename.c_str(), writeVertexNormals); delete tile; if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) { @@ -719,7 +726,7 @@ runTiler(TerrainBuild *command, Grid *grid, TerrainMetadata *metadata) { buildTerrain(tiler, command, threadMetadata); } else if (strcmp(command->outputFormat, "Mesh") == 0) { const MeshTiler tiler(poDataset, *grid, command->tilerOptions, command->meshQualityFactor); - buildMesh(tiler, command, threadMetadata); + buildMesh(tiler, command, threadMetadata, command->vertexNormals); } else { // it's a GDAL format const RasterTiler tiler(poDataset, *grid, command->tilerOptions); buildGDAL(tiler, command, threadMetadata); @@ -762,6 +769,7 @@ main(int argc, char *argv[]) { command.option("-g", "--mesh-qfactor ", "specify the factor to multiply the estimated geometric error to convert heightmaps to irregular meshes. Larger values should mean minor quality. Defaults to 1.0", TerrainBuild::setMeshQualityFactor); command.option("-l", "--layer", "only output the layer.json metadata file", TerrainBuild::setMetadata); command.option("-C", "--cesium-friendly", "Force the creation of missing root tiles to be CesiumJS-friendly", TerrainBuild::setCesiumFriendly); + command.option("-N", "--vertex-normals", "Write 'Oct-Encoded Per-Vertex Normals' for Terrain Lighting, only for `Mesh` format", TerrainBuild::setVertexNormals); command.option("-q", "--quiet", "only output errors", TerrainBuild::setQuiet); command.option("-v", "--verbose", "be more noisy", TerrainBuild::setVerbose); From 1f5c1fbedf82d0e15f6089adf7b181c15db5bb01 Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Fri, 11 May 2018 16:53:37 +0200 Subject: [PATCH 05/26] Refactoring of code to easily write new output formats It simplifies the code to write new output formats (MBTiles...). It defines a set of abstract classes that anyone could override to implement its own logic when serializing. --- src/CMakeLists.txt | 10 ++ src/CTBFileOutputStream.cpp | 43 +++++++++ src/CTBFileOutputStream.hpp | 60 ++++++++++++ src/CTBFileTileSerializer.cpp | 169 ++++++++++++++++++++++++++++++++++ src/CTBFileTileSerializer.hpp | 74 +++++++++++++++ src/CTBOutputStream.hpp | 40 ++++++++ src/CTBZOutputStream.cpp | 72 +++++++++++++++ src/CTBZOutputStream.hpp | 55 +++++++++++ src/GDALSerializer.hpp | 52 +++++++++++ src/MeshSerializer.hpp | 50 ++++++++++ src/MeshTile.cpp | 94 +++++++++---------- src/MeshTile.hpp | 5 + src/TerrainSerializer.hpp | 50 ++++++++++ src/TerrainTile.cpp | 40 +++----- src/TerrainTile.hpp | 5 + tools/ctb-tile.cpp | 120 ++++++------------------ 16 files changed, 770 insertions(+), 169 deletions(-) create mode 100644 src/CTBFileOutputStream.cpp create mode 100644 src/CTBFileOutputStream.hpp create mode 100644 src/CTBFileTileSerializer.cpp create mode 100644 src/CTBFileTileSerializer.hpp create mode 100644 src/CTBOutputStream.hpp create mode 100644 src/CTBZOutputStream.cpp create mode 100644 src/CTBZOutputStream.hpp create mode 100644 src/GDALSerializer.hpp create mode 100644 src/MeshSerializer.hpp create mode 100644 src/TerrainSerializer.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 49d39d0..30be95c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,6 +21,9 @@ add_library(ctb SHARED GDALTile.cpp GDALTiler.cpp GDALDatasetReader.cpp + CTBFileTileSerializer.cpp + CTBFileOutputStream.cpp + CTBZOutputStream.cpp TerrainTiler.cpp TerrainTile.cpp MeshTiler.cpp @@ -35,9 +38,14 @@ set(HEADERS BoundingSphere.hpp Coordinate.hpp Coordinate3D.hpp + GDALSerializer.hpp GDALTile.hpp GDALTiler.hpp GDALDatasetReader.hpp + CTBFileTileSerializer.hpp + CTBFileOutputStream.hpp + CTBOutputStream.hpp + CTBZOutputStream.hpp GlobalGeodetic.hpp GlobalMercator.hpp Grid.hpp @@ -45,12 +53,14 @@ set(HEADERS HeightFieldChunker.hpp Mesh.hpp MeshIterator.hpp + MeshSerializer.hpp MeshTile.hpp MeshTiler.hpp RasterIterator.hpp RasterTiler.hpp CTBException.hpp TerrainIterator.hpp + TerrainSerializer.hpp TerrainTile.hpp TerrainTiler.hpp Tile.hpp diff --git a/src/CTBFileOutputStream.cpp b/src/CTBFileOutputStream.cpp new file mode 100644 index 0000000..a0dc4a3 --- /dev/null +++ b/src/CTBFileOutputStream.cpp @@ -0,0 +1,43 @@ +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file CTBFileOutputStream.cpp + * @brief This defines the `CTBFileOutputStream` class + */ + +#include "CTBFileOutputStream.hpp" + +using namespace ctb; + +/** + * @details + * Writes a sequence of memory pointed by ptr into the FILE*. + */ +uint32_t +ctb::CTBFileOutputStream::write(const void *ptr, uint32_t size) { + return (uint32_t)fwrite(ptr, size, 1, fp); +} + +/** + * @details + * Writes a sequence of memory pointed by ptr into the ostream. + */ +uint32_t +ctb::CTBStdOutputStream::write(const void *ptr, uint32_t size) { + mstream.write((const char *)ptr, size); + return size; +} diff --git a/src/CTBFileOutputStream.hpp b/src/CTBFileOutputStream.hpp new file mode 100644 index 0000000..a1bc3bd --- /dev/null +++ b/src/CTBFileOutputStream.hpp @@ -0,0 +1,60 @@ +#ifndef CTBFILEOUTPUTSTREAM_HPP +#define CTBFILEOUTPUTSTREAM_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file CTBFileOutputStream.hpp + * @brief This declares and defines the `CTBFileOutputStream` and `CTBStdOutputStream` classes + */ + +#include +#include +#include "CTBOutputStream.hpp" + +namespace ctb { + class CTBFileOutputStream; + class CTBStdOutputStream; +} + +/// Implements CTBOutputStream for `FILE*` objects +class CTB_DLL ctb::CTBFileOutputStream : public ctb::CTBOutputStream { +public: + CTBFileOutputStream(FILE *fptr): fp(fptr) {} + + /// Writes a sequence of memory pointed by ptr into the stream + virtual uint32_t write(const void *ptr, uint32_t size); + +protected: + /// The underlying FILE* + FILE *fp; +}; + +/// Implements CTBOutputStream for `std::ostream` objects +class CTB_DLL ctb::CTBStdOutputStream : public ctb::CTBOutputStream { +public: + CTBStdOutputStream(std::ostream &stream): mstream(stream) {} + + /// Writes a sequence of memory pointed by ptr into the stream + virtual uint32_t write(const void *ptr, uint32_t size); + +protected: + /// The underlying std::ostream + std::ostream &mstream; +}; + +#endif /* CTBFILEOUTPUTSTREAM_HPP */ diff --git a/src/CTBFileTileSerializer.cpp b/src/CTBFileTileSerializer.cpp new file mode 100644 index 0000000..2f24fd3 --- /dev/null +++ b/src/CTBFileTileSerializer.cpp @@ -0,0 +1,169 @@ +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file CTBFileTileSerializer.cpp + * @brief This defines the `CTBFileTileSerializer` class + */ + +#include +#include +#include + +#include "../deps/concat.hpp" +#include "cpl_vsi.h" +#include "CTBException.hpp" +#include "CTBFileTileSerializer.hpp" + +#include "CTBFileOutputStream.hpp" +#include "CTBZOutputStream.hpp" + +using namespace std; +using namespace ctb; + +#ifdef _WIN32 +static const char *osDirSep = "\\"; +#else +static const char *osDirSep = "/"; +#endif + + +/// Create a filename for a tile coordinate +std::string +ctb::CTBFileTileSerializer::getTileFilename(const TileCoordinate *coord, const string dirname, const char *extension) { + static mutex mutex; + VSIStatBufL stat; + string filename = concat(dirname, coord->zoom, osDirSep, coord->x); + + lock_guard lock(mutex); + + // Check whether the `{zoom}/{x}` directory exists or not + if (VSIStatExL(filename.c_str(), &stat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG)) { + filename = concat(dirname, coord->zoom); + + // Check whether the `{zoom}` directory exists or not + if (VSIStatExL(filename.c_str(), &stat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG)) { + // Create the `{zoom}` directory + if (VSIMkdir(filename.c_str(), 0755)) + throw CTBException("Could not create the zoom level directory"); + + } else if (!VSI_ISDIR(stat.st_mode)) { + throw CTBException("Zoom level file path is not a directory"); + } + + // Create the `{zoom}/{x}` directory + filename += concat(osDirSep, coord->x); + if (VSIMkdir(filename.c_str(), 0755)) + throw CTBException("Could not create the x level directory"); + + } else if (!VSI_ISDIR(stat.st_mode)) { + throw CTBException("X level file path is not a directory"); + } + + // Create the filename itself, adding the extension if required + filename += concat(osDirSep, coord->y); + if (extension != NULL) { + filename += "."; + filename += extension; + } + + return filename; +} + +/// Check if file exists +static bool +fileExists(const std::string& filename) { + VSIStatBufL statbuf; + return VSIStatExL(filename.c_str(), &statbuf, VSI_STAT_EXISTS_FLAG) == 0; +} + + +/** + * @details + * Returns if the specified Tile Coordinate should be serialized + */ +bool ctb::CTBFileTileSerializer::mustSerializeCoordinate(const ctb::TileCoordinate *coordinate) { + if (!mresume) + return true; + + const string filename = getTileFilename(coordinate, moutputDir, "terrain"); + return !fileExists(filename); +} + +/** + * @details + * Serialize a GDALTile to the Directory store + */ +bool +ctb::CTBFileTileSerializer::serializeTile(const ctb::GDALTile *tile, GDALDriver *driver, const char *extension, CPLStringList &creationOptions) { + const TileCoordinate *coordinate = tile; + const string filename = getTileFilename(coordinate, moutputDir, extension); + const string temp_filename = concat(filename, ".tmp"); + + GDALDataset *poDstDS; + poDstDS = driver->CreateCopy(temp_filename.c_str(), tile->dataset, FALSE, creationOptions, NULL, NULL); + + // Close the datasets, flushing data to destination + if (poDstDS == NULL) { + throw CTBException("Could not create GDAL tile"); + } + GDALClose(poDstDS); + + if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) { + throw new CTBException("Could not rename temporary file"); + } + return true; +} + +/** + * @details + * Serialize a TerrainTile to the Directory store + */ +bool +ctb::CTBFileTileSerializer::serializeTile(const ctb::TerrainTile *tile) { + const TileCoordinate *coordinate = tile; + const string filename = getTileFilename(tile, moutputDir, "terrain"); + const string temp_filename = concat(filename, ".tmp"); + + CTBZFileOutputStream ostream(temp_filename.c_str()); + tile->writeFile(ostream); + ostream.close(); + + if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) { + throw new CTBException("Could not rename temporary file"); + } + return true; +} + +/** + * @details + * Serialize a MeshTile to the Directory store + */ +bool +ctb::CTBFileTileSerializer::serializeTile(const ctb::MeshTile *tile, bool writeVertexNormals) { + const TileCoordinate *coordinate = tile; + const string filename = getTileFilename(coordinate, moutputDir, "terrain"); + const string temp_filename = concat(filename, ".tmp"); + + CTBZFileOutputStream ostream(temp_filename.c_str()); + tile->writeFile(ostream, writeVertexNormals); + ostream.close(); + + if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) { + throw new CTBException("Could not rename temporary file"); + } + return true; +} diff --git a/src/CTBFileTileSerializer.hpp b/src/CTBFileTileSerializer.hpp new file mode 100644 index 0000000..b723428 --- /dev/null +++ b/src/CTBFileTileSerializer.hpp @@ -0,0 +1,74 @@ +#ifndef CTBFILETILESERIALIZER_HPP +#define CTBFILETILESERIALIZER_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file CTBFileTileSerializer.hpp + * @brief This declares and defines the `CTBFileTileSerializer` class + */ + +#include + +#include "TileCoordinate.hpp" +#include "GDALSerializer.hpp" +#include "TerrainSerializer.hpp" +#include "MeshSerializer.hpp" + +namespace ctb { + class CTBFileTileSerializer; +} + +/// Implements a serializer of `Tile`s based in a directory of files +class CTB_DLL ctb::CTBFileTileSerializer : + public ctb::GDALSerializer, + public ctb::TerrainSerializer, + public ctb::MeshSerializer { +public: + CTBFileTileSerializer(const std::string &outputDir, bool resume): + moutputDir(outputDir), + mresume(resume) {} + + /// Start a new serialization task + virtual void startSerialization() {}; + + /// Returns if the specified Tile Coordinate should be serialized + virtual bool mustSerializeCoordinate(const ctb::TileCoordinate *coordinate); + + /// Serialize a GDALTile to the store + virtual bool serializeTile(const ctb::GDALTile *tile, GDALDriver *driver, const char *extension, CPLStringList &creationOptions); + /// Serialize a TerrainTile to the store + virtual bool serializeTile(const ctb::TerrainTile *tile); + /// Serialize a MeshTile to the store + virtual bool serializeTile(const ctb::MeshTile *tile, bool writeVertexNormals = false); + + /// Serialization finished, releases any resources loaded + virtual void endSerialization() {}; + + + /// Create a filename for a tile coordinate + static std::string + getTileFilename(const TileCoordinate *coord, const std::string dirname, const char *extension); + +protected: + /// The target directory where serializing + std::string moutputDir; + /// Do not overwrite existing files + bool mresume; +}; + +#endif /* CTBFILETILESERIALIZER_HPP */ diff --git a/src/CTBOutputStream.hpp b/src/CTBOutputStream.hpp new file mode 100644 index 0000000..7d0d10f --- /dev/null +++ b/src/CTBOutputStream.hpp @@ -0,0 +1,40 @@ +#ifndef CTBOUTPUTSTREAM_HPP +#define CTBOUTPUTSTREAM_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file CTBOutputStream.hpp + * @brief This declares and defines the `CTBOutputStream` class + */ + +#include "config.hpp" +#include "types.hpp" + +namespace ctb { + class CTBOutputStream; +} + +/// This represents a generic CTB output stream to write raw data +class CTB_DLL ctb::CTBOutputStream { +public: + + /// Writes a sequence of memory pointed by ptr into the stream + virtual uint32_t write(const void *ptr, uint32_t size) = 0; +}; + +#endif /* CTBOUTPUTSTREAM_HPP */ diff --git a/src/CTBZOutputStream.cpp b/src/CTBZOutputStream.cpp new file mode 100644 index 0000000..5420b5a --- /dev/null +++ b/src/CTBZOutputStream.cpp @@ -0,0 +1,72 @@ +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file CTBZOutputStream.cpp + * @brief This defines the `CTBZOutputStream` and `CTBZFileOutputStream` classes + */ + +#include "CTBException.hpp" +#include "CTBZOutputStream.hpp" + +using namespace ctb; + +/** + * @details + * Writes a sequence of memory pointed by ptr into the GZFILE*. + */ +uint32_t +ctb::CTBZOutputStream::write(const void *ptr, uint32_t size) { + if (size == 1) { + int c = *((const char *)ptr); + return gzputc(fp, c) == -1 ? 0 : 1; + } + else { + return gzwrite(fp, ptr, size) == 0 ? 0 : size; + } +} + +ctb::CTBZFileOutputStream::CTBZFileOutputStream(const char *fileName) : CTBZOutputStream(NULL) { + gzFile file = gzopen(fileName, "wb"); + + if (file == NULL) { + throw CTBException("Failed to open file"); + } + fp = file; +} + +ctb::CTBZFileOutputStream::~CTBZFileOutputStream() { + close(); +} + +void +ctb::CTBZFileOutputStream::close() { + + // Try and close the file + if (fp) { + switch (gzclose(fp)) { + case Z_OK: + break; + case Z_STREAM_ERROR: + case Z_ERRNO: + case Z_MEM_ERROR: + case Z_BUF_ERROR: + default: + throw CTBException("Failed to close file"); + } + fp = NULL; + } +} diff --git a/src/CTBZOutputStream.hpp b/src/CTBZOutputStream.hpp new file mode 100644 index 0000000..b4890ba --- /dev/null +++ b/src/CTBZOutputStream.hpp @@ -0,0 +1,55 @@ +#ifndef CTBZOUTPUTSTREAM_HPP +#define CTBZOUTPUTSTREAM_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file CTBZOutputStream.hpp + * @brief This declares and defines the `CTBZOutputStream` class + */ + +#include "zlib.h" +#include "CTBOutputStream.hpp" + +namespace ctb { + class CTBZFileOutputStream; + class CTBZOutputStream; +} + +/// Implements CTBOutputStream for `GZFILE` object +class CTB_DLL ctb::CTBZOutputStream : public ctb::CTBOutputStream { +public: + CTBZOutputStream(gzFile gzptr): fp(gzptr) {} + + /// Writes a sequence of memory pointed by ptr into the stream + virtual uint32_t write(const void *ptr, uint32_t size); + +protected: + /// The underlying GZFILE* + gzFile fp; +}; + +/// Implements CTBOutputStream for gzipped files +class CTB_DLL ctb::CTBZFileOutputStream : public ctb::CTBZOutputStream { +public: + CTBZFileOutputStream(const char *fileName); + ~CTBZFileOutputStream(); + + void close(); +}; + +#endif /* CTBZOUTPUTSTREAM_HPP */ diff --git a/src/GDALSerializer.hpp b/src/GDALSerializer.hpp new file mode 100644 index 0000000..130d8c5 --- /dev/null +++ b/src/GDALSerializer.hpp @@ -0,0 +1,52 @@ +#ifndef GDALSERIALIZER_HPP +#define GDALSERIALIZER_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file GDALSerializer.hpp + * @brief This declares and defines the `GDALSerializer` class + */ + +#include "config.hpp" +#include "TileCoordinate.hpp" +#include "GDALTile.hpp" + +#include "cpl_string.h" + +namespace ctb { + class GDALSerializer; +} + +/// Store `GDALTile`s from a GDAL Dataset +class CTB_DLL ctb::GDALSerializer { +public: + + /// Start a new serialization task + virtual void startSerialization() = 0; + + /// Returns if the specified Tile Coordinate should be serialized + virtual bool mustSerializeCoordinate(const ctb::TileCoordinate *coordinate) = 0; + + /// Serialize a GDALTile to the store + virtual bool serializeTile(const ctb::GDALTile *tile, GDALDriver *driver, const char *extension, CPLStringList &creationOptions) = 0; + + /// Serialization finished, releases any resources loaded + virtual void endSerialization() = 0; +}; + +#endif /* GDALSERIALIZER_HPP */ diff --git a/src/MeshSerializer.hpp b/src/MeshSerializer.hpp new file mode 100644 index 0000000..98d6fa0 --- /dev/null +++ b/src/MeshSerializer.hpp @@ -0,0 +1,50 @@ +#ifndef MESHSERIALIZER_HPP +#define MESHSERIALIZER_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file MeshSerializer.hpp + * @brief This declares and defines the `MeshSerializer` class + */ + +#include "config.hpp" +#include "TileCoordinate.hpp" +#include "MeshTile.hpp" + +namespace ctb { + class MeshSerializer; +} + +/// Store `MeshTile`s from a GDAL Dataset +class CTB_DLL ctb::MeshSerializer { +public: + + /// Start a new serialization task + virtual void startSerialization() = 0; + + /// Returns if the specified Tile Coordinate should be serialized + virtual bool mustSerializeCoordinate(const ctb::TileCoordinate *coordinate) = 0; + + /// Serialize a MeshTile to the store + virtual bool serializeTile(const ctb::MeshTile *tile, bool writeVertexNormals = false) = 0; + + /// Serialization finished, releases any resources loaded + virtual void endSerialization() = 0; +}; + +#endif /* MESHSERIALIZER_HPP */ diff --git a/src/MeshTile.cpp b/src/MeshTile.cpp index 11ad15f..7173dfb 100644 --- a/src/MeshTile.cpp +++ b/src/MeshTile.cpp @@ -24,11 +24,11 @@ #include #include #include "cpl_conv.h" -#include "zlib.h" #include "CTBException.hpp" #include "MeshTile.hpp" #include "BoundingSphere.hpp" +#include "CTBZOutputStream.hpp" using namespace ctb; @@ -117,7 +117,7 @@ static inline int quantizeIndices(const double &origin, const double &factor, co } // Write the edge indices of the mesh -template int writeEdgeIndices(gzFile terrainFile, const Mesh &mesh, double edgeCoord, int componentIndex) { +template int writeEdgeIndices(CTBOutputStream &ostream, const Mesh &mesh, double edgeCoord, int componentIndex) { std::vector indices; std::map ihash; @@ -136,11 +136,11 @@ template int writeEdgeIndices(gzFile terrainFile, const Mesh &mesh, } int edgeCount = indices.size(); - gzwrite(terrainFile, &edgeCount, sizeof(int)); + ostream.write(&edgeCount, sizeof(int)); for (size_t i = 0; i < edgeCount; i++) { T indice = (T)indices[i]; - gzwrite(terrainFile, &indice, sizeof(T)); + ostream.write(&indice, sizeof(T)); } return indices.size(); } @@ -207,11 +207,15 @@ MeshTile::MeshTile(const TileCoordinate &coord): */ void MeshTile::writeFile(const char *fileName, bool writeVertexNormals) const { - gzFile terrainFile = gzopen(fileName, "wb"); + CTBZFileOutputStream ostream(fileName); + writeFile(ostream); +} - if (terrainFile == NULL) { - throw CTBException("Failed to open file"); - } +/** + * @details This writes raw terrain data to an output stream. + */ +void +MeshTile::writeFile(CTBOutputStream &ostream, bool writeVertexNormals) const { // Calculate main header mesh data std::vector cartesianVertices; @@ -236,37 +240,37 @@ MeshTile::writeFile(const char *fileName, bool writeVertexNormals) const { double centerX = cartesianBounds.min.x + 0.5 * (cartesianBounds.max.x - cartesianBounds.min.x); double centerY = cartesianBounds.min.y + 0.5 * (cartesianBounds.max.y - cartesianBounds.min.y); double centerZ = cartesianBounds.min.z + 0.5 * (cartesianBounds.max.z - cartesianBounds.min.z); - gzwrite(terrainFile, ¢erX, sizeof(double)); - gzwrite(terrainFile, ¢erY, sizeof(double)); - gzwrite(terrainFile, ¢erZ, sizeof(double)); + ostream.write(¢erX, sizeof(double)); + ostream.write(¢erY, sizeof(double)); + ostream.write(¢erZ, sizeof(double)); // // The minimum and maximum heights in the area covered by this tile. float minimumHeight = (float)bounds.min.z; float maximumHeight = (float)bounds.max.z; - gzwrite(terrainFile, &minimumHeight, sizeof(float)); - gzwrite(terrainFile, &maximumHeight, sizeof(float)); + ostream.write(&minimumHeight, sizeof(float)); + ostream.write(&maximumHeight, sizeof(float)); // - // The tile’s bounding sphere. The X,Y,Z coordinates are again expressed + // The tile's bounding sphere. The X,Y,Z coordinates are again expressed // in Earth-centered Fixed coordinates, and the radius is in meters. double boundingSphereCenterX = cartesianBoundingSphere.center.x; double boundingSphereCenterY = cartesianBoundingSphere.center.y; double boundingSphereCenterZ = cartesianBoundingSphere.center.z; double boundingSphereRadius = cartesianBoundingSphere.radius; - gzwrite(terrainFile, &boundingSphereCenterX, sizeof(double)); - gzwrite(terrainFile, &boundingSphereCenterY, sizeof(double)); - gzwrite(terrainFile, &boundingSphereCenterZ, sizeof(double)); - gzwrite(terrainFile, &boundingSphereRadius , sizeof(double)); + ostream.write(&boundingSphereCenterX, sizeof(double)); + ostream.write(&boundingSphereCenterY, sizeof(double)); + ostream.write(&boundingSphereCenterZ, sizeof(double)); + ostream.write(&boundingSphereRadius , sizeof(double)); // // The horizon occlusion point, expressed in the ellipsoid-scaled Earth-centered Fixed frame. CRSVertex horizonOcclusionPoint = ocp_fromPoints(cartesianVertices, cartesianBoundingSphere); - gzwrite(terrainFile, &horizonOcclusionPoint.x, sizeof(double)); - gzwrite(terrainFile, &horizonOcclusionPoint.y, sizeof(double)); - gzwrite(terrainFile, &horizonOcclusionPoint.z, sizeof(double)); + ostream.write(&horizonOcclusionPoint.x, sizeof(double)); + ostream.write(&horizonOcclusionPoint.y, sizeof(double)); + ostream.write(&horizonOcclusionPoint.z, sizeof(double)); // # Write mesh vertices (X Y Z components of each vertex): int vertexCount = mMesh.vertices.size(); - gzwrite(terrainFile, &vertexCount, sizeof(int)); + ostream.write(&vertexCount, sizeof(int)); for (int c = 0; c < 3; c++) { double origin = bounds.min[c]; double factor = 0; @@ -275,20 +279,20 @@ MeshTile::writeFile(const char *fileName, bool writeVertexNormals) const { // Move the initial value int u0 = quantizeIndices(origin, factor, mMesh.vertices[0][c]), u1, ud; uint16_t sval = zigZagEncode(u0); - gzwrite(terrainFile, &sval, sizeof(uint16_t)); + ostream.write(&sval, sizeof(uint16_t)); for (size_t i = 1, icount = mMesh.vertices.size(); i < icount; i++) { u1 = quantizeIndices(origin, factor, mMesh.vertices[i][c]); ud = u1 - u0; sval = zigZagEncode(ud); - gzwrite(terrainFile, &sval, sizeof(uint16_t)); + ostream.write(&sval, sizeof(uint16_t)); u0 = u1; } } // # Write mesh indices: int triangleCount = mMesh.indices.size() / 3; - gzwrite(terrainFile, &triangleCount, sizeof(int)); + ostream.write(&triangleCount, sizeof(int)); if (vertexCount > BYTESPLIT) { uint32_t highest = 0; uint32_t code; @@ -296,15 +300,15 @@ MeshTile::writeFile(const char *fileName, bool writeVertexNormals) const { // Write main indices for (size_t i = 0, icount = mMesh.indices.size(); i < icount; i++) { code = highest - mMesh.indices[i]; - gzwrite(terrainFile, &code, sizeof(uint32_t)); + ostream.write(&code, sizeof(uint32_t)); if (code == 0) highest++; } // Write all vertices on the edge of the tile (W, S, E, N) - writeEdgeIndices(terrainFile, mMesh, bounds.min.x, 0); - writeEdgeIndices(terrainFile, mMesh, bounds.min.y, 1); - writeEdgeIndices(terrainFile, mMesh, bounds.max.x, 0); - writeEdgeIndices(terrainFile, mMesh, bounds.max.y, 1); + writeEdgeIndices(ostream, mMesh, bounds.min.x, 0); + writeEdgeIndices(ostream, mMesh, bounds.min.y, 1); + writeEdgeIndices(ostream, mMesh, bounds.max.x, 0); + writeEdgeIndices(ostream, mMesh, bounds.max.y, 1); } else { uint16_t highest = 0; @@ -313,23 +317,23 @@ MeshTile::writeFile(const char *fileName, bool writeVertexNormals) const { // Write main indices for (size_t i = 0, icount = mMesh.indices.size(); i < icount; i++) { code = highest - mMesh.indices[i]; - gzwrite(terrainFile, &code, sizeof(uint16_t)); + ostream.write(&code, sizeof(uint16_t)); if (code == 0) highest++; } // Write all vertices on the edge of the tile (W, S, E, N) - writeEdgeIndices(terrainFile, mMesh, bounds.min.x, 0); - writeEdgeIndices(terrainFile, mMesh, bounds.min.y, 1); - writeEdgeIndices(terrainFile, mMesh, bounds.max.x, 0); - writeEdgeIndices(terrainFile, mMesh, bounds.max.y, 1); + writeEdgeIndices(ostream, mMesh, bounds.min.x, 0); + writeEdgeIndices(ostream, mMesh, bounds.min.y, 1); + writeEdgeIndices(ostream, mMesh, bounds.max.x, 0); + writeEdgeIndices(ostream, mMesh, bounds.max.y, 1); } // # Write 'Oct-Encoded Per-Vertex Normals' for Terrain Lighting: if (writeVertexNormals && triangleCount > 0) { unsigned char extensionId = 1; - gzwrite(terrainFile, &extensionId, sizeof(unsigned char)); + ostream.write(&extensionId, sizeof(unsigned char)); int extensionLength = 2 * vertexCount; - gzwrite(terrainFile, &extensionLength, sizeof(int)); + ostream.write(&extensionLength, sizeof(int)); std::vector normalsPerVertex(vertexCount); std::vector normalsPerFace(triangleCount); @@ -358,22 +362,10 @@ MeshTile::writeFile(const char *fileName, bool writeVertexNormals) const { } for (size_t i = 0; i < vertexCount; i++) { Coordinate xy = octEncode(normalsPerVertex[i].normalize()); - gzwrite(terrainFile, &xy.x, sizeof(unsigned char)); - gzwrite(terrainFile, &xy.y, sizeof(unsigned char)); + ostream.write(&xy.x, sizeof(unsigned char)); + ostream.write(&xy.y, sizeof(unsigned char)); } } - - // Try and close the file - switch (gzclose(terrainFile)) { - case Z_OK: - break; - case Z_STREAM_ERROR: - case Z_ERRNO: - case Z_MEM_ERROR: - case Z_BUF_ERROR: - default: - throw CTBException("Failed to close file"); - } } bool diff --git a/src/MeshTile.hpp b/src/MeshTile.hpp index 5cd3fc9..2300c37 100644 --- a/src/MeshTile.hpp +++ b/src/MeshTile.hpp @@ -27,6 +27,7 @@ #include "Mesh.hpp" #include "TileCoordinate.hpp" #include "Tile.hpp" +#include "CTBOutputStream.hpp" namespace ctb { class MeshTile; @@ -55,6 +56,10 @@ class CTB_DLL ctb::MeshTile : void writeFile(const char *fileName, bool writeVertexNormals = false) const; + /// Write terrain data to an output stream + void + writeFile(CTBOutputStream &ostream, bool writeVertexNormals = false) const; + /// Does the terrain tile have child tiles? bool hasChildren() const; diff --git a/src/TerrainSerializer.hpp b/src/TerrainSerializer.hpp new file mode 100644 index 0000000..521f6a3 --- /dev/null +++ b/src/TerrainSerializer.hpp @@ -0,0 +1,50 @@ +#ifndef TERRAINSERIALIZER_HPP +#define TERRAINSERIALIZER_HPP + +/******************************************************************************* + * Copyright 2018 GeoData + * + * Licensed under the Apache License, Version 2.0 (the "License"); you may not + * use this file except in compliance with the License. You may obtain a copy + * of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + *******************************************************************************/ + +/** + * @file TerrainSerializer.hpp + * @brief This declares and defines the `TerrainSerializer` class + */ + +#include "config.hpp" +#include "TileCoordinate.hpp" +#include "TerrainTile.hpp" + +namespace ctb { + class TerrainSerializer; +} + +/// Store `TerrainTile`s from a GDAL Dataset +class CTB_DLL ctb::TerrainSerializer { +public: + + /// Start a new serialization task + virtual void startSerialization() = 0; + + /// Returns if the specified Tile Coordinate should be serialized + virtual bool mustSerializeCoordinate(const ctb::TileCoordinate *coordinate) = 0; + + /// Serialize a TerrainTile to the store + virtual bool serializeTile(const ctb::TerrainTile *tile) = 0; + + /// Serialization finished, releases any resources loaded + virtual void endSerialization() = 0; +}; + +#endif /* TERRAINSERIALIZER_HPP */ diff --git a/src/TerrainTile.cpp b/src/TerrainTile.cpp index 7b78262..393765c 100644 --- a/src/TerrainTile.cpp +++ b/src/TerrainTile.cpp @@ -28,6 +28,8 @@ #include "TerrainTile.hpp" #include "GlobalGeodetic.hpp" #include "Bounds.hpp" +#include "CTBFileOutputStream.hpp" +#include "CTBZOutputStream.hpp" using namespace ctb; @@ -136,9 +138,8 @@ Terrain::readFile(const char *fileName) { */ void Terrain::writeFile(FILE *fp) const { - fwrite(mHeights.data(), TILE_CELL_SIZE * 2, 1, fp); - fwrite(&mChildren, 1, 1, fp); - fwrite(mMask, mMaskLength, 1, fp); + CTBFileOutputStream ostream(fp); + writeFile(ostream); } /** @@ -146,41 +147,30 @@ Terrain::writeFile(FILE *fp) const { */ void Terrain::writeFile(const char *fileName) const { - gzFile terrainFile = gzopen(fileName, "wb"); + CTBZFileOutputStream ostream(fileName); + writeFile(ostream); +} - if (terrainFile == NULL) { - throw CTBException("Failed to open file"); - } +/** + * @details This writes raw terrain data to an output stream. + */ +void +Terrain::writeFile(CTBOutputStream &ostream) const { // Write the height data - if (gzwrite(terrainFile, mHeights.data(), TILE_CELL_SIZE * 2) == 0) { - gzclose(terrainFile); + if (ostream.write((const void *)mHeights.data(), TILE_CELL_SIZE * 2) != TILE_CELL_SIZE * 2) { throw CTBException("Failed to write height data"); } // Write the child flags - if (gzputc(terrainFile, mChildren) == -1) { - gzclose(terrainFile); + if (ostream.write(&mChildren, 1) != 1) { throw CTBException("Failed to write child flags"); } // Write the water mask - if (gzwrite(terrainFile, mMask, mMaskLength) == 0) { - gzclose(terrainFile); + if (ostream.write(&mMask, mMaskLength) != mMaskLength) { throw CTBException("Failed to write water mask"); } - - // Try and close the file - switch (gzclose(terrainFile)) { - case Z_OK: - break; - case Z_STREAM_ERROR: - case Z_ERRNO: - case Z_MEM_ERROR: - case Z_BUF_ERROR: - default: - throw CTBException("Failed to close file"); - } } std::vector diff --git a/src/TerrainTile.hpp b/src/TerrainTile.hpp index 5941166..ce739ad 100644 --- a/src/TerrainTile.hpp +++ b/src/TerrainTile.hpp @@ -29,6 +29,7 @@ #include "config.hpp" #include "Tile.hpp" #include "TileCoordinate.hpp" +#include "CTBOutputStream.hpp" namespace ctb { class Terrain; @@ -65,6 +66,10 @@ class CTB_DLL ctb::Terrain { void writeFile(const char *fileName) const; + /// Write terrain data to an output stream + void + writeFile(CTBOutputStream &ostream) const; + /// Get the water mask as a boolean mask std::vector mask(); diff --git a/tools/ctb-tile.cpp b/tools/ctb-tile.cpp index 3851ba0..de3bcbd 100644 --- a/tools/ctb-tile.cpp +++ b/tools/ctb-tile.cpp @@ -52,6 +52,7 @@ #include "TerrainIterator.hpp" #include "MeshIterator.hpp" #include "GDALDatasetReader.hpp" +#include "CTBFileTileSerializer.hpp" using namespace std; using namespace ctb; @@ -245,52 +246,6 @@ class TerrainBuild : public Command { bool vertexNormals; }; -/** - * Create a filename for a tile coordinate - * - * This also creates the tile directory structure. - */ -static string -getTileFilename(const TileCoordinate *coord, const string dirname, const char *extension) { - static mutex mutex; - VSIStatBufL stat; - string filename = concat(dirname, coord->zoom, osDirSep, coord->x); - - lock_guard lock(mutex); - - // Check whether the `{zoom}/{x}` directory exists or not - if (VSIStatExL(filename.c_str(), &stat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG)) { - filename = concat(dirname, coord->zoom); - - // Check whether the `{zoom}` directory exists or not - if (VSIStatExL(filename.c_str(), &stat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG)) { - // Create the `{zoom}` directory - if (VSIMkdir(filename.c_str(), 0755)) - throw CTBException("Could not create the zoom level directory"); - - } else if (!VSI_ISDIR(stat.st_mode)) { - throw CTBException("Zoom level file path is not a directory"); - } - - // Create the `{zoom}/{x}` directory - filename += concat(osDirSep, coord->x); - if (VSIMkdir(filename.c_str(), 0755)) - throw CTBException("Could not create the x level directory"); - - } else if (!VSI_ISDIR(stat.st_mode)) { - throw CTBException("X level file path is not a directory"); - } - - // Create the filename itself, adding the extension if required - filename += concat(osDirSep, coord->y); - if (extension != NULL) { - filename += "."; - filename += extension; - } - - return filename; -} - /** * Increment a TilerIterator whilst cooperating between threads * @@ -363,6 +318,10 @@ showProgress(int currentIndex, string filename) { return progressFunc(currentIndex / (double) iteratorSize, message.c_str(), NULL); } +int +showProgress(int currentIndex) { + return progressFunc(currentIndex / (double) iteratorSize, NULL, NULL); +} static bool fileExists(const std::string& filename) { @@ -541,7 +500,7 @@ class TerrainMetadata { /// Output GDAL tiles represented by a tiler to a directory static void -buildGDAL(const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { +buildGDAL(GDALSerializer &serializer, const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(command->outputFormat); if (poDriver == NULL) { @@ -553,7 +512,6 @@ buildGDAL(const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *meta } const char *extension = poDriver->GetMetadataItem(GDAL_DMD_EXTENSION); - const string dirname = string(command->outputDir) + osDirSep; i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom, endZoom = (command->endZoom < 0) ? 0 : command->endZoom; @@ -563,38 +521,22 @@ buildGDAL(const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *meta while (!iter.exhausted()) { const TileCoordinate *coordinate = iter.GridIterator::operator*(); - GDALDataset *poDstDS; - const string filename = getTileFilename(coordinate, dirname, extension); if (metadata) metadata->add(tiler.grid(), coordinate); - if( !command->resume || !fileExists(filename) ) { + if (serializer.mustSerializeCoordinate(coordinate)) { GDALTile *tile = *iter; - const string temp_filename = concat(filename, ".tmp"); - poDstDS = poDriver->CreateCopy(temp_filename.c_str(), tile->dataset, FALSE, - command->creationOptions.List(), NULL, NULL ); + serializer.serializeTile(tile, poDriver, extension, command->creationOptions); delete tile; - - // Close the datasets, flushing data to destination - if (poDstDS == NULL) { - throw CTBException("Could not create GDAL tile"); - } - - GDALClose(poDstDS); - - if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) { - throw new CTBException("Could not rename temporary file"); - } } currentIndex = incrementIterator(iter, currentIndex); - showProgress(currentIndex, filename); + showProgress(currentIndex); } } /// Output terrain tiles represented by a tiler to a directory static void -buildTerrain(const TerrainTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { - const string dirname = string(command->outputDir) + osDirSep; +buildTerrain(TerrainSerializer &serializer, const TerrainTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom, endZoom = (command->endZoom < 0) ? 0 : command->endZoom; @@ -605,35 +547,28 @@ buildTerrain(const TerrainTiler &tiler, TerrainBuild *command, TerrainMetadata * while (!iter.exhausted()) { const TileCoordinate *coordinate = iter.GridIterator::operator*(); - const string filename = getTileFilename(coordinate, dirname, "terrain"); if (metadata) metadata->add(tiler.grid(), coordinate); - if( !command->resume || !fileExists(filename) ) { + if (serializer.mustSerializeCoordinate(coordinate)) { TerrainTile *tile = iter.operator*(&reader); - const string temp_filename = concat(filename, ".tmp"); - - tile->writeFile(temp_filename.c_str()); + serializer.serializeTile(tile); delete tile; - - if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) { - throw new CTBException("Could not rename temporary file"); - } } currentIndex = incrementIterator(iter, currentIndex); - showProgress(currentIndex, filename); + showProgress(currentIndex); } } /// Output mesh tiles represented by a tiler to a directory static void -buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata, bool writeVertexNormals = false) { - const string dirname = string(command->outputDir) + osDirSep; +buildMesh(MeshSerializer &serializer, const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata, bool writeVertexNormals = false) { i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom, endZoom = (command->endZoom < 0) ? 0 : command->endZoom; // DEBUG Chunker: #if 0 + const string dirname = string(command->outputDir) + osDirSep; TileCoordinate coordinate(13, 8102, 6047); MeshTile *tile = tiler.createMesh(tiler.dataset(), coordinate); // @@ -660,23 +595,16 @@ buildMesh(const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metada while (!iter.exhausted()) { const TileCoordinate *coordinate = iter.GridIterator::operator*(); - const string filename = getTileFilename(coordinate, dirname, "terrain"); if (metadata) metadata->add(tiler.grid(), coordinate); - if( !command->resume || !fileExists(filename) ) { + if (serializer.mustSerializeCoordinate(coordinate)) { MeshTile *tile = iter.operator*(&reader); - const string temp_filename = concat(filename, ".tmp"); - - tile->writeFile(temp_filename.c_str(), writeVertexNormals); + serializer.serializeTile(tile, writeVertexNormals); delete tile; - - if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) { - throw new CTBException("Could not rename temporary file"); - } } currentIndex = incrementIterator(iter, currentIndex); - showProgress(currentIndex, filename); + showProgress(currentIndex); } } @@ -717,24 +645,30 @@ runTiler(TerrainBuild *command, Grid *grid, TerrainMetadata *metadata) { // Metadata of only this thread, it will be joined to global later TerrainMetadata *threadMetadata = metadata ? new TerrainMetadata() : NULL; + // Choose serializer of tiles (Directory of files, MBTiles store...) + CTBFileTileSerializer serializer(string(command->outputDir) + osDirSep, command->resume); + try { + serializer.startSerialization(); + if (command->metadata) { const RasterTiler tiler(poDataset, *grid, command->tilerOptions); buildMetadata(tiler, command, threadMetadata); } else if (strcmp(command->outputFormat, "Terrain") == 0) { const TerrainTiler tiler(poDataset, *grid); - buildTerrain(tiler, command, threadMetadata); + buildTerrain(serializer, tiler, command, threadMetadata); } else if (strcmp(command->outputFormat, "Mesh") == 0) { const MeshTiler tiler(poDataset, *grid, command->tilerOptions, command->meshQualityFactor); - buildMesh(tiler, command, threadMetadata, command->vertexNormals); + buildMesh(serializer, tiler, command, threadMetadata, command->vertexNormals); } else { // it's a GDAL format const RasterTiler tiler(poDataset, *grid, command->tilerOptions); - buildGDAL(tiler, command, threadMetadata); + buildGDAL(serializer, tiler, command, threadMetadata); } } catch (CTBException &e) { cerr << "Error: " << e.what() << endl; } + serializer.endSerialization(); GDALClose(poDataset); From 3471487aa865f95349d97269334863e99c30b710 Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Thu, 17 May 2018 09:47:28 +0200 Subject: [PATCH 06/26] Add 'octvertexnormals' entry to layer.json file --- tools/ctb-tile.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/tools/ctb-tile.cpp b/tools/ctb-tile.cpp index de3bcbd..5fb03ab 100644 --- a/tools/ctb-tile.cpp +++ b/tools/ctb-tile.cpp @@ -435,7 +435,7 @@ class TerrainMetadata { /// http://help.agi.com/TerrainServer/RESTAPIGuide.html /// Example: /// https://assets.agi.com/stk-terrain/v1/tilesets/world/tiles/layer.json - void writeJsonFile(const std::string &filename, const std::string &datasetName, const std::string &outputFormat = "Terrain", const std::string &profile = "geodetic") const { + void writeJsonFile(const std::string &filename, const std::string &datasetName, const std::string &outputFormat = "Terrain", const std::string &profile = "geodetic", bool writeVertexNormals = false) const { FILE *fp = fopen(filename.c_str(), "w"); if (fp == NULL) { @@ -459,6 +459,9 @@ class TerrainMetadata { } fprintf(fp, " \"attribution\": \"\",\n"); fprintf(fp, " \"schema\": \"tms\",\n"); + if (writeVertexNormals) { + fprintf(fp, " \"extensions\": [ \"octvertexnormals\" ],\n"); + } fprintf(fp, " \"tiles\": [ \"{z}/{x}/{y}.terrain?v={version}\" ],\n"); if (strcmp(profile.c_str(), "geodetic") == 0) { @@ -782,7 +785,7 @@ main(int argc, char *argv[]) { const size_t rfindpos = datasetName.rfind('.'); if (std::string::npos != rfindpos) datasetName = datasetName.erase(rfindpos); - metadata->writeJsonFile(filename, datasetName, std::string(command.outputFormat), std::string(command.profile)); + metadata->writeJsonFile(filename, datasetName, std::string(command.outputFormat), std::string(command.profile), command.vertexNormals); delete metadata; } From 2ebd493f6d9f1f6662d440f95ea763b7f6dce8bc Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Thu, 17 May 2018 13:35:09 +0200 Subject: [PATCH 07/26] Fix available root tiles in layer.json with 'CesiumFriendly' param --- tools/ctb-tile.cpp | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/tools/ctb-tile.cpp b/tools/ctb-tile.cpp index 5fb03ab..2a9e8cb 100644 --- a/tools/ctb-tile.cpp +++ b/tools/ctb-tile.cpp @@ -778,17 +778,6 @@ main(int argc, char *argv[]) { } } - // Write Json metadata file? - if (metadata) { - std::string datasetName(command.getInputFilename()); - datasetName = datasetName.substr(datasetName.find_last_of("/\\") + 1); - const size_t rfindpos = datasetName.rfind('.'); - if (std::string::npos != rfindpos) datasetName = datasetName.erase(rfindpos); - - metadata->writeJsonFile(filename, datasetName, std::string(command.outputFormat), std::string(command.profile), command.vertexNormals); - delete metadata; - } - // CesiumJS friendly? if (command.cesiumFriendly && (strcmp(command.profile, "geodetic") == 0) && command.endZoom <= 0) { std::string dirName0 = string(command.outputDir) + osDirSep + "0" + osDirSep + "0"; @@ -805,6 +794,26 @@ main(int argc, char *argv[]) { VSIMkdir(dirName0.c_str(), 0755); fileCopy(tileName1, tileName0); } + + // Fix available indexes. + if (metadata && metadata->levels.size() > 0) { + TerrainMetadata::LevelInfo &level = metadata->levels.at(0); + level.startX = 0; + level.startY = 0; + level.finalX = 1; + level.finalY = 0; + } + } + + // Write Json metadata file? + if (metadata) { + std::string datasetName(command.getInputFilename()); + datasetName = datasetName.substr(datasetName.find_last_of("/\\") + 1); + const size_t rfindpos = datasetName.rfind('.'); + if (std::string::npos != rfindpos) datasetName = datasetName.erase(rfindpos); + + metadata->writeJsonFile(filename, datasetName, std::string(command.outputFormat), std::string(command.profile), command.vertexNormals); + delete metadata; } return 0; From 4bb78ba2e2c13cb0b4631b355cafb1aa94676f65 Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Thu, 24 May 2018 09:42:42 +0200 Subject: [PATCH 08/26] Create missing root tiles from scratch --- tools/ctb-tile.cpp | 121 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 102 insertions(+), 19 deletions(-) diff --git a/tools/ctb-tile.cpp b/tools/ctb-tile.cpp index 2a9e8cb..af81caf 100644 --- a/tools/ctb-tile.cpp +++ b/tools/ctb-tile.cpp @@ -255,9 +255,9 @@ class TerrainBuild : public Command { * to ensure all tiles are iterated over consecutively. It assumes individual * tile iterators point to the same source GDAL dataset. */ +static int globalIteratorIndex = 0; // keep track of where we are globally template int incrementIterator(T &iter, int currentIndex) { - static int globalIteratorIndex = 0; // keep track of where we are globally static mutex mutex; // ensure iterations occur serially between threads lock_guard lock(mutex); @@ -501,6 +501,71 @@ class TerrainMetadata { } }; +/// Create an empty root temporary elevation file (GTiff) +static std::string +createEmptyRootElevationFile(std::string &fileName, const Grid &grid, const TileCoordinate& coord) { + GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + + if (poDriver == NULL) { + throw CTBException("Could not retrieve GTiff GDAL driver"); + } + + // Create the geo transform for this temporary elevation tile. + // We apply an 1-degree negative offset to avoid problems in borders. + CRSBounds tileBounds = grid.tileBounds(coord); + tileBounds.setMinX(tileBounds.getMinX() + 1); + tileBounds.setMinY(tileBounds.getMinY() + 1); + tileBounds.setMaxX(tileBounds.getMaxX() - 1); + tileBounds.setMaxY(tileBounds.getMaxY() - 1); + const i_tile tileSize = grid.tileSize() - 2; + const double resolution = tileBounds.getWidth() / tileSize; + double adfGeoTransform[6] = { tileBounds.getMinX(), resolution, 0, tileBounds.getMaxY(), 0, -resolution }; + + // Create the spatial reference system for the file + OGRSpatialReference oSRS; + if (oSRS.importFromEPSG(4326) != OGRERR_NONE) { + throw CTBException("Could not create EPSG:4326 spatial reference"); + } + char *pszDstWKT = NULL; + if (oSRS.exportToWkt(&pszDstWKT) != OGRERR_NONE) { + CPLFree(pszDstWKT); + throw CTBException("Could not create EPSG:4326 WKT string"); + } + + // Create the GTiff file + fileName += ".tif"; + GDALDataset *poDataset = poDriver->Create(fileName.c_str(), tileSize, tileSize, 1, GDT_Float32, NULL); + + // Set the projection + if (poDataset->SetProjection(pszDstWKT) != CE_None) { + poDataset->Release(); + CPLFree(pszDstWKT); + throw CTBException("Could not set projection on temporary elevation file"); + } + CPLFree(pszDstWKT); + + // Apply the geo transform + if (poDataset->SetGeoTransform(adfGeoTransform) != CE_None) { + poDataset->Release(); + throw CTBException("Could not set projection on temporary elevation file"); + } + + // Finally write the height data + float *rasterHeights = (float *)CPLCalloc(tileSize * tileSize, sizeof(float)); + GDALRasterBand *heightsBand = poDataset->GetRasterBand(1); + if (heightsBand->RasterIO(GF_Write, 0, 0, tileSize, tileSize, + (void *)rasterHeights, tileSize, tileSize, GDT_Float32, + 0, 0) != CE_None) { + CPLFree(rasterHeights); + throw CTBException("Could not write heights on temporary elevation file"); + } + CPLFree(rasterHeights); + + poDataset->FlushCache(); + poDataset->Release(); + return fileName; +} + /// Output GDAL tiles represented by a tiler to a directory static void buildGDAL(GDALSerializer &serializer, const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { @@ -638,8 +703,8 @@ buildMetadata(const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata * * This function is designed to be run in a separate thread. */ static int -runTiler(TerrainBuild *command, Grid *grid, TerrainMetadata *metadata) { - GDALDataset *poDataset = (GDALDataset *) GDALOpen(command->getInputFilename(), GA_ReadOnly); +runTiler(const char *inputFilename, TerrainBuild *command, Grid *grid, TerrainMetadata *metadata) { + GDALDataset *poDataset = (GDALDataset *) GDALOpen(inputFilename, GA_ReadOnly); if (poDataset == NULL) { cerr << "Error: could not open GDAL dataset" << endl; return 1; @@ -757,9 +822,9 @@ main(int argc, char *argv[]) { // Instantiate the threads using futures from a packaged_task for (int i = 0; i < threadCount ; ++i) { - packaged_task task(runTiler); // wrap the function - tasks.push_back(task.get_future()); // get a future - thread(move(task), &command, &grid, metadata).detach(); // launch on a thread + packaged_task task(runTiler); // wrap the function + tasks.push_back(task.get_future()); // get a future + thread(move(task), command.getInputFilename(), &command, &grid, metadata).detach(); // launch on a thread } // Synchronise the completion of the threads @@ -780,19 +845,37 @@ main(int argc, char *argv[]) { // CesiumJS friendly? if (command.cesiumFriendly && (strcmp(command.profile, "geodetic") == 0) && command.endZoom <= 0) { - std::string dirName0 = string(command.outputDir) + osDirSep + "0" + osDirSep + "0"; - std::string dirName1 = string(command.outputDir) + osDirSep + "0" + osDirSep + "1"; - std::string tileName0 = dirName0 + osDirSep + "0.terrain"; - std::string tileName1 = dirName1 + osDirSep + "0.terrain"; - - if (fileExists(tileName0) && !fileExists(tileName1)) { - VSIMkdir(dirName1.c_str(), 0755); - fileCopy(tileName0, tileName1); - } - else - if (!fileExists(tileName0) && fileExists(tileName1)) { - VSIMkdir(dirName0.c_str(), 0755); - fileCopy(tileName1, tileName0); + + // Create missing root tiles if it is necessary + if (!command.metadata) { + std::string dirName0 = string(command.outputDir) + osDirSep + "0" + osDirSep + "0"; + std::string dirName1 = string(command.outputDir) + osDirSep + "0" + osDirSep + "1"; + std::string tileName0 = dirName0 + osDirSep + "0.terrain"; + std::string tileName1 = dirName1 + osDirSep + "0.terrain"; + + i_zoom missingZoom = 65535; + ctb::TileCoordinate missingTileCoord(missingZoom, 0, 0); + std::string missingTileName; + + if (fileExists(tileName0) && !fileExists(tileName1)) { + VSIMkdir(dirName1.c_str(), 0755); + missingTileCoord = ctb::TileCoordinate(0, 1, 0); + missingTileName = tileName1; + } + else + if (!fileExists(tileName0) && fileExists(tileName1)) { + VSIMkdir(dirName0.c_str(), 0755); + missingTileCoord = ctb::TileCoordinate(0, 0, 0); + missingTileName = tileName0; + } + if (missingTileCoord.zoom != missingZoom) { + globalIteratorIndex = 0; // reset global iterator index + command.startZoom = 0; + command.endZoom = 0; + missingTileName = createEmptyRootElevationFile(missingTileName, grid, missingTileCoord); + runTiler (missingTileName.c_str(), &command, &grid, NULL); + VSIUnlink(missingTileName.c_str()); + } } // Fix available indexes. From ce2da57d00643663b914e1ef2385fb7f200d290b Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Sat, 20 Oct 2018 09:03:39 +0200 Subject: [PATCH 09/26] README adjustments --- README.md | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/README.md b/README.md index 04798cc..ec1d67e 100644 --- a/README.md +++ b/README.md @@ -61,22 +61,26 @@ Usage: ctb-tile [options] GDAL_DATASOURCE Options: - -V, --version output program version - -h, --help output help information - -o, --output-dir specify the output directory for the tiles (defaults to working directory) - -f, --output-format specify the output format for the tiles. This is either `Terrain` (the default) or any format listed by `gdalinfo --formats` - -p, --profile specify the TMS profile for the tiles. This is either `geodetic` (the default) or `mercator` - -c, --thread-count specify the number of threads to use for tile generation. On multicore machines this defaults to the number of CPUs - -t, --tile-size specify the size of the tiles in pixels. This defaults to 65 for terrain tiles and 256 for other GDAL formats - -s, --start-zoom specify the zoom level to start at. This should be greater than the end zoom level - -e, --end-zoom specify the zoom level to end at. This should be less than the start zoom level and >= 0 - -r, --resampling-method specify the raster resampling algorithm. One of: nearest; bilinear; cubic; cubicspline; lanczos; average; mode; max; min; med; q1; q3. Defaults to average. - -n, --creation-option