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..c0463d5 --- /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 unsigned short int 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..4077a2b 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,9 @@ class TerrainBuild : public Command { startZoom(-1), endZoom(-1), verbosity(1), - resume(false) + resume(false), + meshQualityFactor(1.0), + metadata(false) {} void @@ -198,6 +201,16 @@ 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; + } + const char *outputDir, *outputFormat, *profile; @@ -212,6 +225,9 @@ class TerrainBuild : public Command { CPLStringList creationOptions; TilerOptions tilerOptions; + + double meshQualityFactor; + bool metadata; }; /** @@ -423,6 +439,177 @@ buildTerrain(const TerrainTiler &tiler, TerrainBuild *command) { } } +/// Output mesh tiles represented by a tiler to a directory +static void +buildMesh(const MeshTiler &tiler, TerrainBuild *command) { + 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( !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); + } +} + +/// 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 +static void +buildMetadata(const RasterTiler &tiler, TerrainBuild *command) { + const string dirname = string(command->outputDir) + osDirSep; + i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom, + endZoom = (command->endZoom < 0) ? 0 : command->endZoom; + + // Defines the tiles that are available in a level of 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 addTileCoordinate(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); + } + }; + + 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); + + const std::string filename = concat(dirname, "layer.json"); + CRSBounds bounds; + std::vector levels; + for (size_t i = 0; i <= startZoom; i++) { + levels.push_back(LevelInfo()); + } + + RasterIterator iter(tiler, startZoom, endZoom); + int currentIndex = incrementIterator(iter, 0); + setIteratorSize(iter); + + while (!iter.exhausted()) { + const TileCoordinate *coordinate = iter.GridIterator::operator*(); + + CRSBounds tileBounds = tiler.grid().tileBounds(*coordinate); + LevelInfo &level = levels[coordinate->zoom]; + level.addTileCoordinate(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())); + } + currentIndex = incrementIterator(iter, currentIndex); + showProgress(currentIndex, filename); + } + + // Write metadata file + 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(command->outputFormat, "Terrain") == 0) { + fprintf(fp, " \"format\": \"heightmap-1.0\",\n"); + } + else if (strcmp(command->outputFormat, "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"); + fprintf(fp, " \"projection\": \"EPSG:4326\",\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); +} + /** * Perform a tile building operation * @@ -437,9 +624,15 @@ runTiler(TerrainBuild *command, Grid *grid) { } try { - if (strcmp(command->outputFormat, "Terrain") == 0) { + if (command->metadata) { + const RasterTiler tiler(poDataset, *grid, command->tilerOptions); + buildMetadata(tiler, command); + } else if (strcmp(command->outputFormat, "Terrain") == 0) { const TerrainTiler tiler(poDataset, *grid); buildTerrain(tiler, command); + } else if (strcmp(command->outputFormat, "Mesh") == 0) { + const MeshTiler tiler(poDataset, *grid, command->tilerOptions, command->meshQualityFactor); + buildMesh(tiler, command); } else { // it's a GDAL format const RasterTiler tiler(poDataset, *grid, command->tilerOptions); buildGDAL(tiler, command); @@ -460,7 +653,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 +664,8 @@ 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("-q", "--quiet", "only output errors", TerrainBuild::setQuiet); command.option("-v", "--verbose", "be more noisy", TerrainBuild::setVerbose); @@ -513,6 +708,7 @@ main(int argc, char *argv[]) { // Run the tilers in separate threads vector> tasks; int threadCount = (command.threadCount > 0) ? command.threadCount : CPLGetNumCPUs(); + if (command.metadata) threadCount = 1; // Instantiate the threads using futures from a packaged_task for (int i = 0; i < threadCount ; ++i) {