Cache first row and column of pixels from neighbors. (#32)

* Cache first row and column of pixels from neighbors.

In a large prominence run, this can remove the need for 3/4
of the disk reads.

Also iterate east to west to get more cache hits, as we
read our eastern neighbor for pixels.

* Option to filter points *outside* a polygon.

code/filter_points.cpp

@@ -47,12 +47,14 @@ static void usage() {
   printf("  Options:\n");
   printf("  -a longitude    Add 360 to any longitudes < this value in polygon.\n");
   printf("                  Allows polygon to cross antimeridian (negative longitude gets +360)\n");
+  printf("  -x              Print only lines *outside* the polygon instead\n");
   printf("\n");
   exit(1);
 }
 
 int main(int argc, char **argv) {
   double wrapLongitude = -180;
+  bool includeInside = true;
 
   // Parse options
   START_EASYLOGGINGPP(argc, argv);
@@ -63,11 +65,14 @@ int main(int argc, char **argv) {
     {"v", required_argument, nullptr, 0},
     {nullptr, 0, 0, 0},
   };
-  while ((ch = getopt_long(argc, argv, "a:", long_options, nullptr)) != -1) {
+  while ((ch = getopt_long(argc, argv, "a:x", long_options, nullptr)) != -1) {
     switch (ch) {
     case 'a':
       wrapLongitude = atof(optarg);
       break;
+    case 'x':
+      includeInside = false;
+      break;
     }
   }
   argc -= optind;
@@ -115,15 +120,18 @@ int main(int argc, char **argv) {
     double lng = stod(elements[1]);
     LatLng point(lat, lng);
 
-    if (filter.isPointInside(point)) {
+    bool isInside = filter.isPointInside(point);
+    if (isInside) {
       VLOG(1) << "Point is in polygon: " <<  lat << ", " <<  lng;
       numPointsInPolygon += 1;
-
-      outputFile << line << std::endl;
     } else {
       VLOG(1) << "Point is not in polygon: " <<  lat << ", " <<  lng;
       numPointsNotInPolygon += 1;
     }
+
+    if ((isInside && includeInside) || (!isInside && !includeInside)) {
+      outputFile << line << std::endl;
+    }
   }
 
   outputFile.close();

code/prominence.cpp

@@ -204,8 +204,10 @@ int main(int argc, char **argv) {
   // Use double precision to avoid accumulating floating-point error during loop
   double lat = bounds[0];
   while (lat < bounds[1]) {
-    double lng = bounds[2];
-    while (lng < bounds[3]) {
+    // Go east-to-west to take advantage of some tile caching, since we need pixels
+    // from our eastern neighbor.
+    double lng = bounds[3];
+    while (lng >= bounds[2]) {
       // Allow specifying longitude ranges that span the antimeridian (lng > 180)
       auto wrappedLng = lng;
       if (!fileFormat.isUtm() && wrappedLng >= 180) {
@@ -227,7 +229,7 @@ int main(int argc, char **argv) {
             }));
       }
 
-      lng += fileFormat.degreesAcross();
+      lng -= fileFormat.degreesAcross();
     }
     lat += fileFormat.degreesAcross();
   }

code/tile_cache.cpp

@@ -34,9 +34,20 @@ using std::string;
 TileCache::TileCache(TileLoadingPolicy *policy, int maxEntries)
     : mCache(maxEntries),
       mLoadingPolicy(policy) {
+  // Allow loading policy to use the cache.  This is a circular dependency,
+  // but allows for a nice optimization where the cache can store pixels
+  // for neighboring tiles that the loading policy needs (for prominence).
+  policy->setTileCache(this);
 }
 
 TileCache::~TileCache() {
+  // Clear cached first rows and columns
+  for (auto row : mFirstRows) {
+    delete [] row.second;
+  }
+  for (auto col : mFirstCols) {
+    delete [] col.second;
+  }
 }
 
 Tile *TileCache::getOrLoad(double minLat, double minLng,
@@ -125,6 +136,19 @@ Tile *TileCache::loadWithoutCaching(double minLat, double minLng,
     }
   }
 
+  // As an optimization, remember the first row and column
+  auto key = makeCacheKey(minLat, minLng);
+  Elevation *firstRow = new Elevation[tile->width()];
+  Elevation *firstCol = new Elevation[tile->height()];
+  for (auto i = 0; i < tile->width(); ++i) {
+    firstRow[i] = tile->get(i, 0);
+  }
+  for (auto i = 0; i < tile->height(); ++i) {
+    firstCol[i] = tile->get(0, i);
+  }
+  mFirstRows[key] = firstRow;
+  mFirstCols[key] = firstCol;
+
   VLOG(1) << "Loaded tile at " << minLat << " " << minLng << " with max elevation " << tile->maxElevation();
 
   return tile;
@@ -148,6 +172,32 @@ bool TileCache::getMaxElevation(double lat, double lng, Elevation *elev) {
   return retval;
 }
 
+bool TileCache::getFirstRow(double lat, double lng, Elevation **elevs) {
+  return getFirstRowOrCol(lat, lng, mFirstRows, elevs);
+}
+
+bool TileCache::getFirstColumn(double lat, double lng, Elevation **elevs) {
+  return getFirstRowOrCol(lat, lng, mFirstCols, elevs);
+}
+
+bool TileCache::getFirstRowOrCol(double lat, double lng, const CachedElevations &cache, Elevation **elevs) {
+  assert(elevs != nullptr);
+
+  bool retval = true;
+  int key = makeCacheKey(lat, lng);
+
+  mLock.lock();
+  auto it = cache.find(key);
+  if (it == cache.end()) {
+    retval = false;
+  } else {
+    *elevs = it->second;
+  }
+  mLock.unlock();
+
+  return retval;
+}
+
 int TileCache::makeCacheKey(double minLat, double minLng) const {
   // Round to some reasonable precision
   int latKey = static_cast<int>(minLat * 100);

code/tile_cache.h

@@ -52,7 +52,15 @@ class TileCache {
   // If we've ever loaded the tile with the given minimum lat/lng, set elev to its maximum
   // elevation and return true, otherwise return false.
   bool getMaxElevation(double lat, double lng, Elevation *elev);
-
+
+  // If we've ever loaded the tile with the given minimum lat/lng, set elevs to its
+  // first row and return true, otherwise return false.
+  bool getFirstRow(double lat, double lng, Elevation **elevs);
+
+  // If we've ever loaded the tile with the given minimum lat/lng, set elevs to its
+  // first column and return true, otherwise return false.
+  bool getFirstColumn(double lat, double lng, Elevation **elevs);
+
 private:
 
   Lock mLock;
@@ -61,9 +69,17 @@ class TileCache {
   // Map of encoded lat/lng to max elevation in that tile
   std::unordered_map<int, Elevation> mMaxElevations;
 
+  // Map of encoded lat/lng to first row/col in that tile.
+  // This is an optimization for the prominence algorithm.
+  typedef std::unordered_map<int, Elevation *> CachedElevations;
+  CachedElevations mFirstRows;
+  CachedElevations mFirstCols;
+
   Tile *loadInternal(double minLat, double minLng) const;
 
   int makeCacheKey(double minLat, double minLng) const;
+
+  bool getFirstRowOrCol(double lat, double lng, const CachedElevations &cache, Elevation **elevs);
 };
 
 #endif  // _TILE_CACHE_H_

code/tile_loading_policy.cpp

@@ -28,7 +28,7 @@
 #include "flt_loader.h"
 #include "glo_loader.h"
 #include "hgt_loader.h"
-#include "tile.h"
+#include "tile_cache.h"
 #include "easylogging++.h"
 
 using std::string;
@@ -38,7 +38,8 @@ BasicTileLoadingPolicy::BasicTileLoadingPolicy(const string &directory,
     : mDirectory(directory),
       mFileFormat(format),
       mNeighborEdgeLoadingEnabled(false),
-      mUtmZone(-1) {
+      mUtmZone(-1),
+      mTileCache(nullptr) {
 }
 
 void BasicTileLoadingPolicy::enableNeighborEdgeLoading(bool enabled) {
@@ -49,6 +50,10 @@ void BasicTileLoadingPolicy::setUtmZone(int utmZone) {
   mUtmZone = utmZone;
 }
 
+void BasicTileLoadingPolicy::setTileCache(TileCache *cache) {
+  mTileCache = cache;
+}
+
 Tile *BasicTileLoadingPolicy::loadTile(double minLat, double minLng) const {
   Tile *tile = loadInternal(minLat, minLng);
   if (tile == nullptr) {
@@ -191,33 +196,60 @@ Tile *BasicTileLoadingPolicy::appendPixelsFromNeighbors(Tile *tile, double minLa
     }
   }
 
+  Elevation *firstRow, *firstColumn;
+
   auto tileSpan = mFileFormat.degreesAcross();
   auto bottomLat = minLat - tileSpan;
-  Tile *neighbor = loadInternal(bottomLat, minLng);  // bottom neighbor
-  if (neighbor != nullptr) {
+
+  // The tile cache stores the first row and column of each tile it has loaded
+  // so that we can get them here without having to load the entire neighboring
+  // tile from disk again.
+  if (mTileCache != nullptr &&
+      mTileCache->getFirstRow(bottomLat, minLng, &firstRow)) {
     for (int i = 0; i < oldWidth; ++i) {
-      samples[oldHeight * newWidth + i] = neighbor->get(i, 0);
+      samples[oldHeight * newWidth + i] = firstRow[i];
+    }
+  } else {
+    Tile *neighbor = loadInternal(bottomLat, minLng);  // bottom neighbor
+    if (neighbor != nullptr) {
+      for (int i = 0; i < oldWidth; ++i) {
+        samples[oldHeight * newWidth + i] = neighbor->get(i, 0);
+      }
+      delete neighbor;
     }
-    delete neighbor;
   }
 
   auto rightLng = minLng + tileSpan;
   if (rightLng == 180) {
     rightLng = -180;  // antimeridian
   }
-  neighbor = loadInternal(minLat, rightLng);  // right neighbor
-  if (neighbor != nullptr) {
+
+  if (mTileCache != nullptr &&
+      mTileCache->getFirstColumn(minLat, rightLng, &firstColumn)) {
     for (int i = 0; i < oldHeight; ++i) {
-      samples[i * newWidth + oldWidth] = neighbor->get(0, i);
+      samples[i * newWidth + oldWidth] = firstColumn[i];
+    }
+  } else {
+    Tile *neighbor = loadInternal(minLat, rightLng);  // right neighbor
+    if (neighbor != nullptr) {
+      for (int i = 0; i < oldHeight; ++i) {
+        samples[i * newWidth + oldWidth] = neighbor->get(0, i);
+      }
+      delete neighbor;
     }
-    delete neighbor;
   }
 
   // Have to get a single corner pixel from the SE neighbor--annoying
-  neighbor = loadInternal(bottomLat, rightLng);
-  if (neighbor != nullptr) {
-    samples[newHeight * newWidth - 1] = neighbor->get(0, 0);
-    delete neighbor;
+  // Optimization: maybe we have this pixel cached.
+  if (mTileCache != nullptr &&
+      mTileCache->getFirstRow(bottomLat, rightLng, &firstRow)) {
+    samples[newHeight * newWidth - 1] = firstRow[0];
+  } else {
+    Tile *neighbor = loadInternal(bottomLat, rightLng);
+    if (neighbor != nullptr) {
+      samples[newHeight * newWidth - 1] = neighbor->get(0, 0);
+      delete neighbor;
+    }
   }
 
   Tile *newTile = new Tile(newWidth, newHeight, samples);

code/tile_loading_policy.h

@@ -31,11 +31,16 @@
 
 #include <string>
 
+class TileCache;
+
 // Responsible for loading a tile given location.
 
 class TileLoadingPolicy {
 public:
   virtual Tile *loadTile(double minLat, double minLng) const = 0;
+
+  // Optionally get access to the tile cache
+  virtual void setTileCache(TileCache *cache) = 0;
 };
 
 
@@ -60,11 +65,14 @@ class BasicTileLoadingPolicy : public TileLoadingPolicy {
 
   virtual Tile *loadTile(double minLat, double minLng) const;
 
+  virtual void setTileCache(TileCache *cache);
+
 private:
   std::string mDirectory;  // Directory for loading tiles
   FileFormat mFileFormat;  
   bool mNeighborEdgeLoadingEnabled;
   int mUtmZone;
+  TileCache *mTileCache;
 
   // Load tile without modifications
   Tile *loadInternal(double minLat, double minLng) const;