load correct covering tiles in perspective view

port pointCoordinate from -js to convert screen points to tile
coordinates correctly in perspective view.

8 files changed, 140 insertions, 54 deletions

diff --git a/include/mbgl/util/mat4.hpp b/include/mbgl/util/mat4.hpp
index ea4854306c..9c34332400 100644
--- a/include/mbgl/util/mat4.hpp
+++ b/include/mbgl/util/mat4.hpp
@@ -32,6 +32,7 @@ typedef std::array<double, 16> mat4;
 namespace matrix {
 
 void identity(mat4& out);
+bool invert(mat4& out, mat4& a);
 void ortho(mat4& out, double left, double right, double bottom, double top, double near, double far);
 void perspective(mat4& out, double fovy, double aspect, double near, double far);
 void copy(mat4& out, const mat4& a);
diff --git a/src/mbgl/map/source.cpp b/src/mbgl/map/source.cpp
index a3caab291e..f60ddce573 100644
--- a/src/mbgl/map/source.cpp
+++ b/src/mbgl/map/source.cpp
@@ -9,6 +9,7 @@
 #include <mbgl/storage/response.hpp>
 #include <mbgl/util/math.hpp>
 #include <mbgl/util/box.hpp>
+#include <mbgl/util/tile_coordinate.hpp>
 #include <mbgl/util/mapbox.hpp>
 #include <mbgl/storage/file_source.hpp>
 #include <mbgl/style/style_layer.hpp>
@@ -332,14 +333,14 @@ std::forward_list<TileID> Source::coveringTiles(const TransformState& state) con
 
     // Map four viewport corners to pixel coordinates
     box points = state.cornersToBox(z);
-    const vec2<double>& center = points.center;
+    const TileCoordinate center = state.pointCoordinate({ state.getWidth() / 2.0f, state.getHeight()/ 2.0f }).zoomTo(z);
 
     std::forward_list<TileID> covering_tiles = tileCover(z, points, reparseOverscaled ? actualZ : z);
 
     covering_tiles.sort([&center](const TileID& a, const TileID& b) {
         // Sorts by distance from the box center
-        return std::fabs(a.x - center.x) + std::fabs(a.y - center.y) <
-               std::fabs(b.x - center.x) + std::fabs(b.y - center.y);
+        return std::fabs(a.x - center.column) + std::fabs(a.y - center.row) <
+               std::fabs(b.x - center.column) + std::fabs(b.y - center.row);
     });
 
     return covering_tiles;
diff --git a/src/mbgl/map/transform_state.cpp b/src/mbgl/map/transform_state.cpp
index 1cb1d9f827..c06d944d97 100644
--- a/src/mbgl/map/transform_state.cpp
+++ b/src/mbgl/map/transform_state.cpp
@@ -3,6 +3,8 @@
 #include <mbgl/util/projection.hpp>
 #include <mbgl/util/constants.hpp>
 #include <mbgl/util/box.hpp>
+#include <mbgl/util/tile_coordinate.hpp>
+#include <mbgl/util/interpolate.hpp>
 
 using namespace mbgl;
 
@@ -17,37 +19,37 @@ void TransformState::matrixFor(mat4& matrix, const TileID& id, const int8_t z) c
     matrix::scale(matrix, matrix, s / 4096.0f, s / 4096.0f, 1);
 }
 
-box TransformState::cornersToBox(uint32_t z) const {
-    const double ref_scale = std::pow(2, z);
-
-    const double angle_sin = std::sin(-angle);
-    const double angle_cos = std::cos(-angle);
-
-    const double w_2 = static_cast<double>(width) / 2.0;
-    const double h_2 = static_cast<double>(height) / 2.0;
-    const double ss_0 = scale * util::tileSize;
-    const double ss_1 = ref_scale / ss_0;
-    const double ss_2 = ss_0 / 2.0;
-
-    // Calculate the corners of the map view. The resulting coordinates will be
-    // in fractional tile coordinates.
-    box b;
+void TransformState::getProjMatrix(mat4& projMatrix) const {
+    double halfFov = std::atan(0.5 / getAltitude());
+    double topHalfSurfaceDistance = std::sin(halfFov) * getAltitude() /
+        std::sin(M_PI / 2.0f - getPitch() - halfFov);
+    // Calculate z value of the farthest fragment that should be rendered.
+    double farZ = std::cos(M_PI / 2.0f - getPitch()) * topHalfSurfaceDistance + getAltitude();
 
-    b.tl.x = ((-w_2) * angle_cos - (-h_2) * angle_sin + ss_2 - x) * ss_1;
-    b.tl.y = ((-w_2) * angle_sin + (-h_2) * angle_cos + ss_2 - y) * ss_1;
+    matrix::perspective(projMatrix, 2.0f * std::atan((getHeight() / 2.0f) / getAltitude()),
+            double(getWidth()) / getHeight(), 0.1, farZ);
 
-    b.tr.x = ((+w_2) * angle_cos - (-h_2) * angle_sin + ss_2 - x) * ss_1;
-    b.tr.y = ((+w_2) * angle_sin + (-h_2) * angle_cos + ss_2 - y) * ss_1;
+    matrix::translate(projMatrix, projMatrix, 0, 0, -getAltitude());
 
-    b.bl.x = ((-w_2) * angle_cos - (+h_2) * angle_sin + ss_2 - x) * ss_1;
-    b.bl.y = ((-w_2) * angle_sin + (+h_2) * angle_cos + ss_2 - y) * ss_1;
+    // After the rotateX, z values are in pixel units. Convert them to
+    // altitude unites. 1 altitude unit = the screen height.
+    matrix::scale(projMatrix, projMatrix, 1, -1, 1.0f / getHeight());
 
-    b.br.x = ((+w_2) * angle_cos - (+h_2) * angle_sin + ss_2 - x) * ss_1;
-    b.br.y = ((+w_2) * angle_sin + (+h_2) * angle_cos + ss_2 - y) * ss_1;
+    matrix::rotate_x(projMatrix, projMatrix, getPitch());
+    matrix::rotate_z(projMatrix, projMatrix, getAngle());
 
-    b.center.x = (ss_2 - x) * ss_1;
-    b.center.y = (ss_2 - y) * ss_1;
+    matrix::translate(projMatrix, projMatrix, pixel_x() - getWidth() / 2.0f,
+            pixel_y() - getWidth() / 2.0f, 0);
+}
 
+box TransformState::cornersToBox(uint32_t z) const {
+    double w = width;
+    double h = height;
+    box b(
+    pointCoordinate({ 0, 0 }).zoomTo(z),
+    pointCoordinate({ w, 0 }).zoomTo(z),
+    pointCoordinate({ w, h }).zoomTo(z),
+    pointCoordinate({ 0, h }).zoomTo(z));
     return b;
 }
 
@@ -240,6 +242,61 @@ const LatLng TransformState::latLngForPixel(const vec2<double> pixel) const {
     return Projection::latLngForProjectedMeters(givenMeters);
 }
 
+TileCoordinate TransformState::pointCoordinate(const vec2<double> point) const {
+
+    float targetZ = 0;
+    float tileZoom = std::floor(getZoom());
+
+    mat4 mat = coordinatePointMatrix(tileZoom);
+
+    mat4 inverted;
+    bool err = matrix::invert(inverted, mat);
+
+    if (err) throw std::runtime_error("failed to invert coordinatePointMatrix");
+
+    // since we don't know the correct projected z value for the point,
+    // unproject two points to get a line and then find the point on that
+    // line with z=0
+
+    matrix::vec4 coord0;
+    matrix::vec4 coord1;
+    matrix::vec4 point0 = {{ point.x, point.y, 0, 1 }};
+    matrix::vec4 point1 = {{ point.x, point.y, 1, 1 }};
+    matrix::transformMat4(coord0, point0, inverted);
+    matrix::transformMat4(coord1, point1, inverted);
+
+    float w0 = coord0[3];
+    float w1 = coord1[3];
+    float x0 = coord0[0] / w0;
+    float x1 = coord1[0] / w1;
+    float y0 = coord0[1] / w0;
+    float y1 = coord1[1] / w1;
+    float z0 = coord0[2] / w0;
+    float z1 = coord1[2] / w1;
+
+
+    float t = z0 == z1 ? 0 : (targetZ - z0) / (z1 - z0);
+  
+    return { util::interpolate(x0, x1, t), util::interpolate(y0, y1, t), tileZoom };
+}
+
+mat4 TransformState::coordinatePointMatrix(float z) const {
+    mat4 proj;
+    getProjMatrix(proj);
+    float s = util::tileSize * scale / std::pow(2, z);
+    matrix::scale(proj, proj, s , s, 1);
+    matrix::multiply(proj, getPixelMatrix(), proj);
+    return proj;
+}
+
+mat4 TransformState::getPixelMatrix() const {
+    mat4 m;
+    matrix::identity(m);
+    matrix::scale(m, m, width / 2.0f, -height / 2.0f, 1);
+    matrix::translate(m, m, 1, -1, 0);
+    return m;
+}
+
 
 #pragma mark - Changing
 
diff --git a/src/mbgl/map/transform_state.hpp b/src/mbgl/map/transform_state.hpp
index 92e2fb8097..5b648ad07d 100644
--- a/src/mbgl/map/transform_state.hpp
+++ b/src/mbgl/map/transform_state.hpp
@@ -1,10 +1,11 @@
 #ifndef MBGL_MAP_TRANSFORM_STATE
 #define MBGL_MAP_TRANSFORM_STATE
 
-#include <mbgl/util/mat4.hpp>
 #include <mbgl/util/geo.hpp>
 #include <mbgl/util/constants.hpp>
 #include <mbgl/util/vec.hpp>
+#include <mbgl/util/mat4.hpp>
+#include <mbgl/util/vec4.hpp>
 
 #include <cstdint>
 #include <array>
@@ -14,6 +15,7 @@ namespace mbgl {
 
 class TileID;
 struct box;
+struct TileCoordinate;
 
 class TransformState {
     friend class Transform;
@@ -21,6 +23,7 @@ class TransformState {
 public:
     // Matrix
     void matrixFor(mat4& matrix, const TileID& id, const int8_t z) const;
+    void getProjMatrix(mat4& matrix) const;
     box cornersToBox(uint32_t z) const;
 
     // Dimensions
@@ -62,6 +65,9 @@ public:
     double pixel_x() const;
     double pixel_y() const;
 
+    // Conversions
+    TileCoordinate pointCoordinate(const vec2<double> point) const;
+
 private:
     void constrain(double& scale, double& y) const;
 
@@ -72,6 +78,9 @@ private:
     // logical dimensions
     uint16_t width = 0, height = 0;
 
+    mat4 coordinatePointMatrix(float z) const;
+    mat4 getPixelMatrix() const;
+
 private:
     // animation state
     bool rotating = false;
diff --git a/src/mbgl/renderer/painter.cpp b/src/mbgl/renderer/painter.cpp
index 7f06b2345c..74268a4017 100644
--- a/src/mbgl/renderer/painter.cpp
+++ b/src/mbgl/renderer/painter.cpp
@@ -133,26 +133,7 @@ void Painter::lineWidth(float line_width) {
 
 void Painter::changeMatrix() {
 
-    double halfFov = std::atan(0.5 / state.getAltitude());
-    double topHalfSurfaceDistance = std::sin(halfFov) * state.getAltitude() /
-        std::sin(M_PI / 2.0f - state.getPitch() - halfFov);
-    // Calculate z value of the farthest fragment that should be rendered.
-    double farZ = std::cos(M_PI / 2.0f - state.getPitch()) * topHalfSurfaceDistance + state.getAltitude();
-
-    matrix::perspective(projMatrix, 2.0f * std::atan((state.getHeight() / 2.0f) / state.getAltitude()),
-            double(state.getWidth()) / state.getHeight(), 0.1, farZ);
-
-    matrix::translate(projMatrix, projMatrix, 0, 0, -state.getAltitude());
-
-    // After the rotateX, z values are in pixel units. Convert them to
-    // altitude unites. 1 altitude unit = the screen height.
-    matrix::scale(projMatrix, projMatrix, 1, -1, 1.0f / state.getHeight());
-
-    matrix::rotate_x(projMatrix, projMatrix, state.getPitch());
-    matrix::rotate_z(projMatrix, projMatrix, state.getAngle());
-
-    matrix::translate(projMatrix, projMatrix, state.pixel_x() - state.getWidth() / 2.0f,
-            state.pixel_y() - state.getWidth() / 2.0f, 0);
+    state.getProjMatrix(projMatrix);
 
     // The extrusion matrix.
     matrix::ortho(extrudeMatrix, 0, state.getWidth(), state.getHeight(), 0, 0, -1);
diff --git a/src/mbgl/util/box.hpp b/src/mbgl/util/box.hpp
index 55a5d46fbc..f03adc2e50 100644
--- a/src/mbgl/util/box.hpp
+++ b/src/mbgl/util/box.hpp
@@ -1,13 +1,14 @@
 #ifndef MBGL_UTIL_BOX
 #define MBGL_UTIL_BOX
 
-#include <mbgl/util/vec.hpp>
+#include <mbgl/util/tile_coordinate.hpp>
 
 namespace mbgl {
 
 struct box {
-    vec2<double> tl, tr, bl, br;
-    vec2<double> center;
+    box(TileCoordinate tl_, TileCoordinate tr_, TileCoordinate br_, TileCoordinate bl_) :
+        tl(tl_), tr(tr_), br(br_), bl(bl_) {}
+    TileCoordinate tl, tr, br, bl;
 };
 
 }
diff --git a/src/mbgl/util/tile_coordinate.hpp b/src/mbgl/util/tile_coordinate.hpp
new file mode 100644
index 0000000000..0f5f0dcb8b
--- /dev/null
+++ b/src/mbgl/util/tile_coordinate.hpp
@@ -0,0 +1,29 @@
+#ifndef MBGL_UTIL_TILE_COORDINATE
+#define MBGL_UTIL_TILE_COORDINATE
+
+#include <mbgl/util/vec.hpp>
+
+namespace mbgl {
+
+struct TileCoordinate {
+    float column;
+    float row;
+    float zoom;
+
+    TileCoordinate(float column_, float row_, float zoom_) :
+        column(column_), row(row_), zoom(zoom_) {}
+
+    TileCoordinate zoomTo(float targetZoom) {
+        float scale = std::pow(2, targetZoom - zoom);
+        return { column * scale, row * scale, targetZoom };
+    }
+
+    TileCoordinate operator-(TileCoordinate c) {
+        c = c.zoomTo(zoom);
+        return { column - c.column, row - c.row, zoom };
+    };
+};
+
+}
+
+#endif
diff --git a/src/mbgl/util/tile_cover.cpp b/src/mbgl/util/tile_cover.cpp
index 9e444fec12..b7937bf849 100644
--- a/src/mbgl/util/tile_cover.cpp
+++ b/src/mbgl/util/tile_cover.cpp
@@ -1,6 +1,7 @@
 #include <mbgl/util/tile_cover.hpp>
 #include <mbgl/util/vec.hpp>
 #include <mbgl/util/box.hpp>
+#include <mbgl/util/tile_coordinate.hpp>
 
 namespace mbgl {
 
@@ -78,13 +79,19 @@ std::forward_list<TileID> tileCover(int8_t z, const mbgl::box &bounds, int8_t ac
         }
     };
 
+    mbgl::vec2<double> tl = { bounds.tl.column, bounds.tl.row };
+    mbgl::vec2<double> tr = { bounds.tr.column, bounds.tr.row };
+    mbgl::vec2<double> br = { bounds.br.column, bounds.br.row };
+    mbgl::vec2<double> bl = { bounds.bl.column, bounds.bl.row };
+
     // Divide the screen up in two triangles and scan each of them:
     // \---+
     // | \ |
     // +---\.
-    scanTriangle(bounds.tl, bounds.tr, bounds.br, 0, tiles, scanLine);
-    scanTriangle(bounds.br, bounds.bl, bounds.tl, 0, tiles, scanLine);
+    scanTriangle(tl, tr, br, 0, tiles, scanLine);
+    scanTriangle(br, bl, tl, 0, tiles, scanLine);
 
+    t.sort();
     t.unique();
 
     return t;