Add MBTiles output and support for skipping empty tiles

src/BoundingSphere.hpp

@@ -0,0 +1,184 @@
+#ifndef BBSPHERE_HPP
+#define BBSPHERE_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <[email protected]>
+ *
+ * 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 <[email protected]>
+ */
+
+#include <vector>
+#include <limits>
+
+#include "Coordinate3D.hpp"
+#include "types.hpp"
+
+namespace ctb {
+  template <class T> class BoundingSphere;
+  template <class T> class BoundingBox;
+}
+
+/// A spherical bounding region which is defined by a center point and a radius
+template <class T>
+class ctb::BoundingSphere {
+public:
+  Coordinate3D<T> 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<Coordinate3D<T>> &points) {
+    fromPoints(points);
+  }
+
+  /// Calculate the center and radius from the specified point stream
+  /// Based on Ritter's algorithm
+  void fromPoints(const std::vector<Coordinate3D<T>> &points) {
+    const T MAX =  std::numeric_limits<T>::infinity();
+    const T MIN = -std::numeric_limits<T>::infinity();
+
+    Coordinate3D<T> minPointX(MAX, MAX, MAX);
+    Coordinate3D<T> minPointY(MAX, MAX, MAX);
+    Coordinate3D<T> minPointZ(MAX, MAX, MAX);
+    Coordinate3D<T> maxPointX(MIN, MIN, MIN);
+    Coordinate3D<T> maxPointY(MIN, MIN, MIN);
+    Coordinate3D<T> 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<T> &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<T> diameter1 = minPointX;
+    Coordinate3D<T> diameter2 = maxPointX;
+    T maxSpan = xSpan;
+    if (ySpan > maxSpan) {
+      diameter1 = minPointY;
+      diameter2 = maxPointY;
+      maxSpan = ySpan;
+    }
+    if (zSpan > maxSpan) {
+      diameter1 = minPointZ;
+      diameter2 = maxPointZ;
+      maxSpan = zSpan;
+    }
+
+    Coordinate3D<T> ritterCenter = Coordinate3D<T>(
+      (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<T> minBoxPt(minPointX.x, minPointY.y, minPointZ.z);
+    Coordinate3D<T> maxBoxPt(maxPointX.x, maxPointY.y, maxPointZ.z);
+    Coordinate3D<T> naiveCenter = (minBoxPt + maxBoxPt) * 0.5;
+    T naiveRadius = 0;
+
+    for (int i = 0, icount = points.size(); i < icount; i++) {
+      const Coordinate3D<T> &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 T>
+class ctb::BoundingBox {
+public:
+  Coordinate3D<T> min; ///< The min coordinate of the BoundingBox
+  Coordinate3D<T> max; ///< The max coordinate of the BoundingBox
+
+  /// Create an empty BoundingBox
+  BoundingBox() {
+  }
+  /// Create a BoundingBox from the specified point stream
+  BoundingBox(const std::vector<Coordinate3D<T>> &points) {
+    fromPoints(points);
+  }
+
+  /// Calculate the BBOX from the specified point stream
+  void fromPoints(const std::vector<Coordinate3D<T>> &points) {
+    const T MAX =  std::numeric_limits<T>::infinity();
+    const T MIN = -std::numeric_limits<T>::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<T> &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 */

src/CMakeLists.txt

@@ -22,20 +22,29 @@ 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})
 
 # 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

src/Coordinate3D.hpp

@@ -0,0 +1,154 @@
+#ifndef COORDINATE3D_HPP
+#define COORDINATE3D_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <[email protected]>
+ *
+ * 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 <vector>
+#include <cmath>
+
+/**
+ * @file Coordinate3D.hpp
+ * @brief This declares and defines the `Coordinate3D` class
+ * @author Alvaro Huarte <[email protected]>
+ */
+
+namespace ctb {
+  template <class T> class Coordinate3D;
+}
+
+/// A representation of a 3-dimensional point coordinate
+template <class T>
+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<T> cross(const Coordinate3D<T> &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<T> &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<T> add(const Coordinate3D<T> &p1, const Coordinate3D<T> &p2) {
+    return p1 + p2;
+  }
+  inline static Coordinate3D<T> subtract(const Coordinate3D<T> &p1, const Coordinate3D<T> &p2) {
+    return p1 - p2;
+  }
+  inline static T distanceSquared(const Coordinate3D<T> &p1, const Coordinate3D<T> &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<T> &p1, const Coordinate3D<T> &p2) {
+    return std::sqrt(distanceSquared(p1, p2));
+  }
+  inline Coordinate3D<T> normalize(void) const {
+    T mgn = magnitude();
+    return Coordinate3D(x / mgn, y / mgn, z / mgn);
+  }
+};
+
+#endif /* COORDINATE3D_HPP */