Merge pull request #2787 from rouault/vgrid_perf_improvements

GeoTIFF grid reading: perf improvements (fixes #2785)

8 files changed, 252 insertions, 187 deletions

diff --git a/include/proj/internal/lru_cache.hpp b/include/proj/internal/lru_cache.hpp
index b7aff6b9..b2e997b3 100644
--- a/include/proj/internal/lru_cache.hpp
+++ b/include/proj/internal/lru_cache.hpp
@@ -160,6 +160,21 @@ class Cache {
     keys_.splice(keys_.begin(), keys_, iter->second);
     return iter->second->value;
   }
+
+  /**
+   *	The const reference returned here is only
+   *    guaranteed to be valid till the next insert/delete
+   */
+  const Value* getPtr(const Key& k) {
+    Guard g(lock_);
+    const auto iter = cache_.find(k);
+    if (iter == cache_.end()) {
+      return nullptr;
+    }
+    keys_.splice(keys_.begin(), keys_, iter->second);
+    return &(iter->second->value);
+  }
+
   /**
    * returns a copy of the stored object (if found)
    */
diff --git a/src/fwd.cpp b/src/fwd.cpp
index 97ba5999..8583a6e7 100644
--- a/src/fwd.cpp
+++ b/src/fwd.cpp
@@ -38,9 +38,12 @@
 #define OUTPUT_UNITS P->right
 
 
-static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
+static void fwd_prepare (PJ *P, PJ_COORD& coo) {
     if (HUGE_VAL==coo.v[0] || HUGE_VAL==coo.v[1] || HUGE_VAL==coo.v[2])
-        return proj_coord_error ();
+    {
+        coo = proj_coord_error ();
+        return;
+    }
 
     /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
     if (HUGE_VAL==coo.v[2] && P->helmert) coo.v[2] = 0.0;
@@ -56,13 +59,15 @@ static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
         {
             proj_log_error(P, _("Invalid latitude"));
             proj_errno_set (P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
-            return proj_coord_error ();
+            coo = proj_coord_error ();
+            return;
         }
         if (coo.lp.lam > 10  ||  coo.lp.lam < -10)
         {
             proj_log_error(P, _("Invalid longitude"));
             proj_errno_set (P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
-            return proj_coord_error ();
+            coo = proj_coord_error ();
+            return;
         }
 
 
@@ -89,7 +94,7 @@ static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
             coo = proj_trans (P->cart,       PJ_INV, coo); /* Go back to angular using local ellps */
         }
         if (coo.lp.lam==HUGE_VAL)
-            return coo;
+            return;
         if (P->vgridshift)
             coo = proj_trans (P->vgridshift, PJ_FWD, coo); /* Go orthometric from geometric */
 
@@ -100,18 +105,18 @@ static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
         if (0==P->over)
             coo.lp.lam = adjlon(coo.lp.lam);
 
-        return coo;
+        return;
     }
 
 
     /* We do not support gridshifts on cartesian input */
     if (INPUT_UNITS==PJ_IO_UNITS_CARTESIAN && P->helmert)
-            return proj_trans (P->helmert, PJ_INV, coo);
-    return coo;
+            coo = proj_trans (P->helmert, PJ_INV, coo);
+    return;
 }
 
 
-static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
+static void fwd_finalize (PJ *P, PJ_COORD& coo) {
 
     switch (OUTPUT_UNITS) {
 
@@ -161,29 +166,28 @@ static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
 
     if (P->axisswap)
         coo = proj_trans (P->axisswap, PJ_FWD, coo);
-
-    return coo;
 }
 
 
-static PJ_COORD error_or_coord(PJ *P, PJ_COORD coord, int last_errno) {
-    if (proj_errno(P))
+static inline PJ_COORD error_or_coord(PJ *P, PJ_COORD coord, int last_errno) {
+    if (P->ctx->last_errno)
         return proj_coord_error();
 
-    proj_errno_restore(P, last_errno);
+    P->ctx->last_errno = last_errno;
+
     return coord;
 }
 
 
 PJ_XY pj_fwd(PJ_LP lp, PJ *P) {
-    int last_errno;
     PJ_COORD coo = {{0,0,0,0}};
     coo.lp = lp;
 
-    last_errno = proj_errno_reset(P);
+    const int last_errno = P->ctx->last_errno;
+    P->ctx->last_errno = 0;
 
     if (!P->skip_fwd_prepare)
-        coo = fwd_prepare (P, coo);
+        fwd_prepare (P, coo);
     if (HUGE_VAL==coo.v[0] || HUGE_VAL==coo.v[1])
         return proj_coord_error ().xy;
 
@@ -202,7 +206,7 @@ PJ_XY pj_fwd(PJ_LP lp, PJ *P) {
         return proj_coord_error ().xy;
 
     if (!P->skip_fwd_finalize)
-        coo = fwd_finalize (P, coo);
+        fwd_finalize (P, coo);
 
     return error_or_coord(P, coo, last_errno).xy;
 }
@@ -210,14 +214,14 @@ PJ_XY pj_fwd(PJ_LP lp, PJ *P) {
 
 
 PJ_XYZ pj_fwd3d(PJ_LPZ lpz, PJ *P) {
-    int last_errno;
     PJ_COORD coo = {{0,0,0,0}};
     coo.lpz = lpz;
 
-    last_errno = proj_errno_reset(P);
+    const int last_errno = P->ctx->last_errno;
+    P->ctx->last_errno = 0;
 
     if (!P->skip_fwd_prepare)
-        coo = fwd_prepare (P, coo);
+        fwd_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
         return proj_coord_error ().xyz;
 
@@ -236,7 +240,7 @@ PJ_XYZ pj_fwd3d(PJ_LPZ lpz, PJ *P) {
         return proj_coord_error ().xyz;
 
     if (!P->skip_fwd_finalize)
-        coo = fwd_finalize (P, coo);
+        fwd_finalize (P, coo);
 
     return error_or_coord(P, coo, last_errno).xyz;
 }
@@ -244,10 +248,12 @@ PJ_XYZ pj_fwd3d(PJ_LPZ lpz, PJ *P) {
 
 
 PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P) {
-    int last_errno = proj_errno_reset(P);
+
+    const int last_errno = P->ctx->last_errno;
+    P->ctx->last_errno = 0;
 
     if (!P->skip_fwd_prepare)
-        coo = fwd_prepare (P, coo);
+        fwd_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
         return proj_coord_error ();
 
@@ -266,7 +272,7 @@ PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P) {
         return proj_coord_error ();
 
     if (!P->skip_fwd_finalize)
-        coo = fwd_finalize (P, coo);
+        fwd_finalize (P, coo);
 
     return error_or_coord(P, coo, last_errno);
 }
diff --git a/src/grids.cpp b/src/grids.cpp
index 3bde3cad..c49cbb8f 100644
--- a/src/grids.cpp
+++ b/src/grids.cpp
@@ -76,6 +76,13 @@ static void swap_words(void *dataIn, size_t word_size, size_t word_count)
 
 // ---------------------------------------------------------------------------
 
+void ExtentAndRes::computeInvRes() {
+    invResX = 1.0 / resX;
+    invResY = 1.0 / resY;
+}
+
+// ---------------------------------------------------------------------------
+
 bool ExtentAndRes::fullWorldLongitude() const {
     return isGeographic && east - west + resX >= 2 * M_PI - 1e-10;
 }
@@ -126,6 +133,7 @@ static ExtentAndRes globalExtent() {
     extent.north = M_PI / 2;
     extent.resX = M_PI;
     extent.resY = M_PI / 2;
+    extent.computeInvRes();
     return extent;
 }
 
@@ -249,6 +257,7 @@ GTXVerticalShiftGrid *GTXVerticalShiftGrid::open(PJ_CONTEXT *ctx,
     extent.resY = ystep * DEG_TO_RAD;
     extent.east = (xorigin + xstep * (columns - 1)) * DEG_TO_RAD;
     extent.north = (yorigin + ystep * (rows - 1)) * DEG_TO_RAD;
+    extent.computeInvRes();
 
     return new GTXVerticalShiftGrid(ctx, std::move(fp), name, columns, rows,
                                     extent);
@@ -328,54 +337,30 @@ class BlockCache {
   public:
     void insert(uint32_t ifdIdx, uint32_t blockNumber,
                 const std::vector<unsigned char> &data);
-    std::shared_ptr<std::vector<unsigned char>> get(uint32_t ifdIdx,
-                                                    uint32_t blockNumber);
+    const std::vector<unsigned char> *get(uint32_t ifdIdx,
+                                          uint32_t blockNumber);
 
   private:
-    struct Key {
-        uint32_t ifdIdx;
-        uint32_t blockNumber;
-
-        Key(uint32_t ifdIdxIn, uint32_t blockNumberIn)
-            : ifdIdx(ifdIdxIn), blockNumber(blockNumberIn) {}
-        bool operator==(const Key &other) const {
-            return ifdIdx == other.ifdIdx && blockNumber == other.blockNumber;
-        }
-    };
-
-    struct KeyHasher {
-        std::size_t operator()(const Key &k) const {
-            return k.ifdIdx ^ (k.blockNumber << 16) ^ (k.blockNumber >> 16);
-        }
-    };
+    typedef uint64_t Key;
 
     static constexpr int NUM_BLOCKS_AT_CROSSING_TILES = 4;
     static constexpr int MAX_SAMPLE_COUNT = 3;
-    lru11::Cache<
-        Key, std::shared_ptr<std::vector<unsigned char>>, lru11::NullLock,
-        std::unordered_map<
-            Key,
-            typename std::list<lru11::KeyValuePair<
-                Key, std::shared_ptr<std::vector<unsigned char>>>>::iterator,
-            KeyHasher>>
-        cache_{NUM_BLOCKS_AT_CROSSING_TILES * MAX_SAMPLE_COUNT};
+    lru11::Cache<Key, std::vector<unsigned char>, lru11::NullLock> cache_{
+        NUM_BLOCKS_AT_CROSSING_TILES * MAX_SAMPLE_COUNT};
 };
 
 // ---------------------------------------------------------------------------
 
 void BlockCache::insert(uint32_t ifdIdx, uint32_t blockNumber,
                         const std::vector<unsigned char> &data) {
-    cache_.insert(Key(ifdIdx, blockNumber),
-                  std::make_shared<std::vector<unsigned char>>(data));
+    cache_.insert((static_cast<uint64_t>(ifdIdx) << 32) | blockNumber, data);
 }
 
 // ---------------------------------------------------------------------------
 
-std::shared_ptr<std::vector<unsigned char>>
-BlockCache::get(uint32_t ifdIdx, uint32_t blockNumber) {
-    std::shared_ptr<std::vector<unsigned char>> ret;
-    cache_.tryGet(Key(ifdIdx, blockNumber), ret);
-    return ret;
+const std::vector<unsigned char> *BlockCache::get(uint32_t ifdIdx,
+                                                  uint32_t blockNumber) {
+    return cache_.getPtr((static_cast<uint64_t>(ifdIdx) << 32) | blockNumber);
 }
 
 // ---------------------------------------------------------------------------
@@ -388,7 +373,7 @@ class GTiffGrid : public Grid {
     uint32_t m_ifdIdx;
     TIFFDataType m_dt;
     uint16_t m_samplesPerPixel;
-    uint16_t m_planarConfig;
+    uint16_t m_planarConfig; // set to -1 if m_samplesPerPixel == 1
     bool m_bottomUp;
     toff_t m_dirOffset;
     bool m_tiled;
@@ -397,18 +382,19 @@ class GTiffGrid : public Grid {
     mutable std::vector<unsigned char> m_buffer{};
     unsigned m_blocksPerRow = 0;
     unsigned m_blocksPerCol = 0;
-    std::map<int, double> m_mapOffset{};
-    std::map<int, double> m_mapScale{};
+    unsigned m_blocks = 0;
+    std::vector<double> m_adfOffset{};
+    std::vector<double> m_adfScale{};
     std::map<std::pair<int, std::string>, std::string> m_metadata{};
     bool m_hasNodata = false;
+    bool m_blockIs256Pixel = false;
+    bool m_isSingleBlock = false;
     float m_noData = 0.0f;
     uint32_t m_subfileType = 0;
 
     GTiffGrid(const GTiffGrid &) = delete;
     GTiffGrid &operator=(const GTiffGrid &) = delete;
 
-    void getScaleOffset(double &scale, double &offset, uint16_t sample) const;
-
     template <class T>
     float readValue(const std::vector<unsigned char> &buffer,
                     uint32_t offsetInBlock, uint16_t sample) const;
@@ -446,7 +432,9 @@ GTiffGrid::GTiffGrid(PJ_CONTEXT *ctx, TIFF *hTIFF, BlockCache &cache, File *fp,
                      uint16_t planarConfig, bool bottomUpIn)
     : Grid(nameIn, widthIn, heightIn, extentIn), m_ctx(ctx), m_hTIFF(hTIFF),
       m_cache(cache), m_fp(fp), m_ifdIdx(ifdIdx), m_dt(dtIn),
-      m_samplesPerPixel(samplesPerPixelIn), m_planarConfig(planarConfig),
+      m_samplesPerPixel(samplesPerPixelIn),
+      m_planarConfig(samplesPerPixelIn == 1 ? static_cast<uint16_t>(-1)
+                                            : planarConfig),
       m_bottomUp(bottomUpIn), m_dirOffset(TIFFCurrentDirOffset(hTIFF)),
       m_tiled(TIFFIsTiled(hTIFF) != 0) {
 
@@ -460,10 +448,15 @@ GTiffGrid::GTiffGrid(PJ_CONTEXT *ctx, TIFF *hTIFF, BlockCache &cache, File *fp,
             m_blockHeight = m_height;
     }
 
+    m_blockIs256Pixel = (m_blockWidth == 256) && (m_blockHeight == 256);
+    m_isSingleBlock = (m_blockWidth == static_cast<uint32_t>(m_width)) &&
+                      (m_blockHeight == static_cast<uint32_t>(m_height));
+
     TIFFGetField(m_hTIFF, TIFFTAG_SUBFILETYPE, &m_subfileType);
 
     m_blocksPerRow = (m_width + m_blockWidth - 1) / m_blockWidth;
     m_blocksPerCol = (m_height + m_blockHeight - 1) / m_blockHeight;
+    m_blocks = m_blocksPerRow * m_blocksPerCol;
 
     const char *text = nullptr;
     // Poor-man XML parsing of TIFFTAG_GDAL_METADATA tag. Hopefully good
@@ -515,16 +508,26 @@ GTiffGrid::GTiffGrid(PJ_CONTEXT *ctx, TIFF *hTIFF, BlockCache &cache, File *fp,
                     break;
                 const auto role = tag.substr(rolePos, endQuote - rolePos);
                 if (role == "offset") {
-                    if (sample >= 0) {
+                    if (sample >= 0 &&
+                        static_cast<unsigned>(sample) <= m_samplesPerPixel) {
                         try {
-                            m_mapOffset[sample] = c_locale_stod(value);
+                            if (m_adfOffset.empty()) {
+                                m_adfOffset.resize(m_samplesPerPixel);
+                                m_adfScale.resize(m_samplesPerPixel, 1);
+                            }
+                            m_adfOffset[sample] = c_locale_stod(value);
                         } catch (const std::exception &) {
                         }
                     }
                 } else if (role == "scale") {
-                    if (sample >= 0) {
+                    if (sample >= 0 &&
+                        static_cast<unsigned>(sample) <= m_samplesPerPixel) {
                         try {
-                            m_mapScale[sample] = c_locale_stod(value);
+                            if (m_adfOffset.empty()) {
+                                m_adfOffset.resize(m_samplesPerPixel);
+                                m_adfScale.resize(m_samplesPerPixel, 1);
+                            }
+                            m_adfScale[sample] = c_locale_stod(value);
                         } catch (const std::exception &) {
                         }
                     }
@@ -555,33 +558,16 @@ GTiffGrid::~GTiffGrid() = default;
 
 // ---------------------------------------------------------------------------
 
-void GTiffGrid::getScaleOffset(double &scale, double &offset,
-                               uint16_t sample) const {
-    {
-        auto iter = m_mapScale.find(sample);
-        if (iter != m_mapScale.end())
-            scale = iter->second;
-    }
-
-    {
-        auto iter = m_mapOffset.find(sample);
-        if (iter != m_mapOffset.end())
-            offset = iter->second;
-    }
-}
-
-// ---------------------------------------------------------------------------
-
 template <class T>
 float GTiffGrid::readValue(const std::vector<unsigned char> &buffer,
                            uint32_t offsetInBlock, uint16_t sample) const {
     const auto ptr = reinterpret_cast<const T *>(buffer.data());
     assert(offsetInBlock < buffer.size() / sizeof(T));
     const auto val = ptr[offsetInBlock];
-    if (!m_hasNodata || static_cast<float>(val) != m_noData) {
-        double scale = 1;
-        double offset = 0;
-        getScaleOffset(scale, offset, sample);
+    if ((!m_hasNodata || static_cast<float>(val) != m_noData) &&
+        sample < m_adfScale.size()) {
+        double scale = m_adfScale[sample];
+        double offset = m_adfOffset[sample];
         return static_cast<float>(val * scale + offset);
     } else {
         return static_cast<float>(val);
@@ -595,27 +581,41 @@ bool GTiffGrid::valueAt(uint16_t sample, int x, int yFromBottom,
     assert(x >= 0 && yFromBottom >= 0 && x < m_width && yFromBottom < m_height);
     assert(sample < m_samplesPerPixel);
 
-    const int blockX = x / m_blockWidth;
-
     // All non-TIFF grids have the first rows in the file being the one
     // corresponding to the southern-most row. In GeoTIFF, the convention is
     // *generally* different (when m_bottomUp == false), TIFF being an
     // image-oriented image. If m_bottomUp == true, then we had GeoTIFF hints
     // that the first row of the image is the southern-most.
     const int yTIFF = m_bottomUp ? yFromBottom : m_height - 1 - yFromBottom;
-    const int blockY = yTIFF / m_blockHeight;
 
-    uint32_t blockId = blockY * m_blocksPerRow + blockX;
+    int blockXOff;
+    int blockYOff;
+    uint32_t blockId;
+
+    if (m_blockIs256Pixel) {
+        const int blockX = x / 256;
+        blockXOff = x % 256;
+        const int blockY = yTIFF / 256;
+        blockYOff = yTIFF % 256;
+        blockId = blockY * m_blocksPerRow + blockX;
+    } else if (m_isSingleBlock) {
+        blockXOff = x;
+        blockYOff = yTIFF;
+        blockId = 0;
+    } else {
+        const int blockX = x / m_blockWidth;
+        blockXOff = x % m_blockWidth;
+        const int blockY = yTIFF / m_blockHeight;
+        blockYOff = yTIFF % m_blockHeight;
+        blockId = blockY * m_blocksPerRow + blockX;
+    }
+
     if (m_planarConfig == PLANARCONFIG_SEPARATE) {
-        blockId += sample * m_blocksPerCol * m_blocksPerRow;
+        blockId += sample * m_blocks;
     }
 
-    auto cachedBuffer = m_cache.get(m_ifdIdx, blockId);
-    std::vector<unsigned char> *pBuffer = &m_buffer;
-    if (cachedBuffer != nullptr) {
-        // Safe as we don't access the cache before pBuffer is used
-        pBuffer = cachedBuffer.get();
-    } else {
+    const std::vector<unsigned char> *pBuffer = m_cache.get(m_ifdIdx, blockId);
+    if (pBuffer == nullptr) {
         if (TIFFCurrentDirOffset(m_hTIFF) != m_dirOffset &&
             !TIFFSetSubDirectory(m_hTIFF, m_dirOffset)) {
             return false;
@@ -643,6 +643,7 @@ bool GTiffGrid::valueAt(uint16_t sample, int x, int yFromBottom,
             }
         }
 
+        pBuffer = &m_buffer;
         try {
             m_cache.insert(m_ifdIdx, blockId, m_buffer);
         } catch (const std::exception &e) {
@@ -651,8 +652,11 @@ bool GTiffGrid::valueAt(uint16_t sample, int x, int yFromBottom,
         }
     }
 
-    uint32_t offsetInBlock =
-        (x % m_blockWidth) + (yTIFF % m_blockHeight) * m_blockWidth;
+    uint32_t offsetInBlock;
+    if (m_blockIs256Pixel)
+        offsetInBlock = blockXOff + blockYOff * 256U;
+    else
+        offsetInBlock = blockXOff + blockYOff * m_blockWidth;
     if (m_planarConfig == PLANARCONFIG_CONTIG)
         offsetInBlock = offsetInBlock * m_samplesPerPixel + sample;
 
@@ -1053,6 +1057,7 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() {
     extent.resY = fabs(vRes) * mulFactor;
     extent.east = (west + hRes * (width - 1)) * mulFactor;
     extent.south = (north - vRes * (height - 1)) * mulFactor;
+    extent.computeInvRes();
 
     if (vRes < 0) {
         std::swap(extent.north, extent.south);
@@ -1451,7 +1456,7 @@ const VerticalShiftGrid *VerticalShiftGrid::gridAt(double lon,
 const VerticalShiftGrid *VerticalShiftGridSet::gridAt(double lon,
                                                       double lat) const {
     for (const auto &grid : m_grids) {
-        if (dynamic_cast<NullVerticalShiftGrid *>(grid.get())) {
+        if (grid->isNullGrid()) {
             return grid.get();
         }
         const auto &extent = grid->extentAndRes();
@@ -1604,6 +1609,8 @@ NTv1Grid *NTv1Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
     extent.north = to_double(header + 40) * DEG_TO_RAD;
     extent.resX = to_double(header + 104) * DEG_TO_RAD;
     extent.resY = to_double(header + 88) * DEG_TO_RAD;
+    extent.computeInvRes();
+
     if (!(fabs(extent.west) <= 4 * M_PI && fabs(extent.east) <= 4 * M_PI &&
           fabs(extent.north) <= M_PI + 1e-5 &&
           fabs(extent.south) <= M_PI + 1e-5 && extent.west < extent.east &&
@@ -1616,9 +1623,9 @@ NTv1Grid *NTv1Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
         return nullptr;
     }
     const int columns = static_cast<int>(
-        fabs((extent.east - extent.west) / extent.resX + 0.5) + 1);
+        fabs((extent.east - extent.west) * extent.invResX + 0.5) + 1);
     const int rows = static_cast<int>(
-        fabs((extent.north - extent.south) / extent.resY + 0.5) + 1);
+        fabs((extent.north - extent.south) * extent.invResY + 0.5) + 1);
 
     return new NTv1Grid(ctx, std::move(fp), filename, columns, rows, extent);
 }
@@ -1948,6 +1955,7 @@ std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx,
             to_double(header + OFFSET_SOUTH_LAT + 16 * 4) * DEG_TO_RAD / 3600.0;
         extent.resX =
             to_double(header + OFFSET_SOUTH_LAT + 16 * 5) * DEG_TO_RAD / 3600.0;
+        extent.computeInvRes();
 
         if (!(fabs(extent.west) <= 4 * M_PI && fabs(extent.east) <= 4 * M_PI &&
               fabs(extent.north) <= M_PI + 1e-5 &&
@@ -1961,9 +1969,9 @@ std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx,
             return nullptr;
         }
         const int columns = static_cast<int>(
-            fabs((extent.east - extent.west) / extent.resX + 0.5) + 1);
+            fabs((extent.east - extent.west) * extent.invResX + 0.5) + 1);
         const int rows = static_cast<int>(
-            fabs((extent.north - extent.south) / extent.resY + 0.5) + 1);
+            fabs((extent.north - extent.south) * extent.invResY + 0.5) + 1);
 
         pj_log(ctx, PJ_LOG_TRACE,
                "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", gridName.c_str(),
@@ -2429,7 +2437,7 @@ const HorizontalShiftGrid *HorizontalShiftGrid::gridAt(double lon,
 const HorizontalShiftGrid *HorizontalShiftGridSet::gridAt(double lon,
                                                           double lat) const {
     for (const auto &grid : m_grids) {
-        if (dynamic_cast<NullHorizontalShiftGrid *>(grid.get())) {
+        if (grid->isNullGrid()) {
             return grid.get();
         }
         const auto &extent = grid->extentAndRes();
@@ -2760,7 +2768,7 @@ const GenericShiftGrid *GenericShiftGrid::gridAt(double x, double y) const {
 
 const GenericShiftGrid *GenericShiftGridSet::gridAt(double x, double y) const {
     for (const auto &grid : m_grids) {
-        if (dynamic_cast<NullGenericShiftGrid *>(grid.get())) {
+        if (grid->isNullGrid()) {
             return grid.get();
         }
         const auto &extent = grid->extentAndRes();
@@ -3196,7 +3204,7 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids,
     }
 
     /* Interpolation of a location within the grid */
-    double grid_x = (input.lam - extent.west) / extent.resX;
+    double grid_x = (input.lam - extent.west) * extent.invResX;
     if (input.lam < extent.west) {
         if (extent.fullWorldLongitude()) {
             // The first fmod goes to ]-lim, lim[ range
@@ -3205,7 +3213,7 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids,
                               grid->width(),
                           grid->width());
         } else {
-            grid_x = (input.lam + 2 * M_PI - extent.west) / extent.resX;
+            grid_x = (input.lam + 2 * M_PI - extent.west) * extent.invResX;
         }
     } else if (input.lam > extent.east) {
         if (extent.fullWorldLongitude()) {
@@ -3215,10 +3223,10 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids,
                               grid->width(),
                           grid->width());
         } else {
-            grid_x = (input.lam - 2 * M_PI - extent.west) / extent.resX;
+            grid_x = (input.lam - 2 * M_PI - extent.west) * extent.invResX;
         }
     }
-    double grid_y = (input.phi - extent.south) / extent.resY;
+    double grid_y = (input.phi - extent.south) * extent.invResY;
     int grid_ix = static_cast<int>(lround(floor(grid_x)));
     if (!(grid_ix >= 0 && grid_ix < grid->width())) {
         // in the unlikely case we end up here...
@@ -3262,39 +3270,60 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids,
         return HUGE_VAL;
     }
 
-    double total_weight = 0.0;
-    int n_weights = 0;
-    double value = 0.0f;
+    double value = 0.0;
 
-    if (!grid->isNodata(value_a, vmultiplier)) {
-        double weight = (1.0 - grid_x) * (1.0 - grid_y);
-        value += value_a * weight;
-        total_weight += weight;
-        n_weights++;
-    }
-    if (!grid->isNodata(value_b, vmultiplier)) {
-        double weight = (grid_x) * (1.0 - grid_y);
-        value += value_b * weight;
-        total_weight += weight;
-        n_weights++;
-    }
-    if (!grid->isNodata(value_c, vmultiplier)) {
-        double weight = (1.0 - grid_x) * (grid_y);
-        value += value_c * weight;
-        total_weight += weight;
-        n_weights++;
-    }
-    if (!grid->isNodata(value_d, vmultiplier)) {
-        double weight = (grid_x) * (grid_y);
-        value += value_d * weight;
-        total_weight += weight;
-        n_weights++;
-    }
-    if (n_weights == 0) {
+    const double grid_x_y = grid_x * grid_y;
+    const bool a_valid = !grid->isNodata(value_a, vmultiplier);
+    const bool b_valid = !grid->isNodata(value_b, vmultiplier);
+    const bool c_valid = !grid->isNodata(value_c, vmultiplier);
+    const bool d_valid = !grid->isNodata(value_d, vmultiplier);
+    const int countValid =
+        static_cast<int>(a_valid) + static_cast<int>(b_valid) +
+        static_cast<int>(c_valid) + static_cast<int>(d_valid);
+    if (countValid == 4) {
+        {
+            double weight = 1.0 - grid_x - grid_y + grid_x_y;
+            value = value_a * weight;
+        }
+        {
+            double weight = grid_x - grid_x_y;
+            value += value_b * weight;
+        }
+        {
+            double weight = grid_y - grid_x_y;
+            value += value_c * weight;
+        }
+        {
+            double weight = grid_x_y;
+            value += value_d * weight;
+        }
+    } else if (countValid == 0) {
         proj_context_errno_set(ctx, PROJ_ERR_COORD_TRANSFM_GRID_AT_NODATA);
         value = HUGE_VAL;
-    } else if (n_weights != 4)
+    } else {
+        double total_weight = 0.0;
+        if (a_valid) {
+            double weight = 1.0 - grid_x - grid_y + grid_x_y;
+            value = value_a * weight;
+            total_weight = weight;
+        }
+        if (b_valid) {
+            double weight = grid_x - grid_x_y;
+            value += value_b * weight;
+            total_weight += weight;
+        }
+        if (c_valid) {
+            double weight = grid_y - grid_x_y;
+            value += value_c * weight;
+            total_weight += weight;
+        }
+        if (d_valid) {
+            double weight = grid_x_y;
+            value += value_d * weight;
+            total_weight += weight;
+        }
         value /= total_weight;
+    }
 
     return value * vmultiplier;
 }
@@ -3364,8 +3393,10 @@ double pj_vgrid_value(PJ *P, const ListOfVGrids &grids, PJ_LP lp,
     double value;
 
     value = read_vgrid_value(P->ctx, grids, lp, vmultiplier);
-    proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam * RAD_TO_DEG,
-                   lp.phi * RAD_TO_DEG, value);
+    if (pj_log_active(P->ctx, PJ_LOG_TRACE)) {
+        proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f",
+                       lp.lam * RAD_TO_DEG, lp.phi * RAD_TO_DEG, value);
+    }
 
     return value;
 }
@@ -3413,14 +3444,14 @@ bool pj_bilinear_interpolation_three_samples(
     // by identifying the lower-left x,y of it (ix, iy), and the upper-right
     // (ix2, iy2)
 
-    double grid_x = (lp.lam - extent.west) / extent.resX;
+    double grid_x = (lp.lam - extent.west) * extent.invResX;
     // Special case for grids with world extent, and dealing with wrap-around
     if (lp.lam < extent.west) {
-        grid_x = (lp.lam + 2 * M_PI - extent.west) / extent.resX;
+        grid_x = (lp.lam + 2 * M_PI - extent.west) * extent.invResX;
     } else if (lp.lam > extent.east) {
-        grid_x = (lp.lam - 2 * M_PI - extent.west) / extent.resX;
+        grid_x = (lp.lam - 2 * M_PI - extent.west) * extent.invResX;
     }
-    double grid_y = (lp.phi - extent.south) / extent.resY;
+    double grid_y = (lp.phi - extent.south) * extent.invResY;
     int ix = static_cast<int>(grid_x);
     int iy = static_cast<int>(grid_y);
     int ix2 = std::min(ix + 1, grid->width() - 1);
diff --git a/src/grids.hpp b/src/grids.hpp
index d060fc95..459bde07 100644
--- a/src/grids.hpp
+++ b/src/grids.hpp
@@ -45,6 +45,10 @@ struct ExtentAndRes {
     double north;      // in radian for geographic, in CRS units otherwise
     double resX;       // in radian for geographic, in CRS units otherwise
     double resY;       // in radian for geographic, in CRS units otherwise
+    double invResX;    // = 1 / resX;
+    double invResY;    // = 1 / resY;
+
+    void computeInvRes();
 
     bool fullWorldLongitude() const;
     bool contains(const ExtentAndRes &other) const;
diff --git a/src/inv.cpp b/src/inv.cpp
index b4a42189..8925f0e9 100644
--- a/src/inv.cpp
+++ b/src/inv.cpp
@@ -36,10 +36,11 @@
 #define INPUT_UNITS  P->right
 #define OUTPUT_UNITS P->left
 
-static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
+static void inv_prepare (PJ *P, PJ_COORD& coo) {
     if (coo.v[0] == HUGE_VAL || coo.v[1] == HUGE_VAL || coo.v[2] == HUGE_VAL) {
         proj_errno_set (P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
-        return proj_coord_error ();
+        coo = proj_coord_error ();
+        return;
     }
 
     /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
@@ -52,10 +53,10 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
     /* Handle remaining possible input types */
     switch (INPUT_UNITS) {
     case PJ_IO_UNITS_WHATEVER:
-        return coo;
+        break;
 
     case PJ_IO_UNITS_DEGREES:
-        return coo;
+        break;
 
     /* de-scale and de-offset */
     case PJ_IO_UNITS_CARTESIAN:
@@ -65,8 +66,7 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
         if (P->is_geocent) {
             coo = proj_trans (P->cart, PJ_INV, coo);
         }
-
-        return coo;
+        break;
 
     case PJ_IO_UNITS_PROJECTED:
     case PJ_IO_UNITS_CLASSIC:
@@ -74,7 +74,7 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
         coo.xyz.y = P->to_meter  * coo.xyz.y - P->y0;
         coo.xyz.z = P->vto_meter * coo.xyz.z - P->z0;
         if (INPUT_UNITS==PJ_IO_UNITS_PROJECTED)
-            return coo;
+            return;
 
         /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis  */
         /* Multiplying by ra, rather than dividing by a because the CalCOFI projection       */
@@ -82,23 +82,20 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
         /* (CalCOFI avoids further scaling by stomping - but a better solution is possible)  */
         coo.xyz.x *= P->ra;
         coo.xyz.y *= P->ra;
-        return coo;
+        break;
 
     case PJ_IO_UNITS_RADIANS:
         coo.lpz.z = P->vto_meter * coo.lpz.z - P->z0;
         break;
     }
-
-    /* Should not happen, so we could return pj_coord_err here */
-    return coo;
 }
 
 
 
-static PJ_COORD inv_finalize (PJ *P, PJ_COORD coo) {
+static void inv_finalize (PJ *P, PJ_COORD& coo) {
     if (coo.xyz.x == HUGE_VAL) {
         proj_errno_set (P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
-        return proj_coord_error ();
+        coo = proj_coord_error ();
     }
 
     if (OUTPUT_UNITS==PJ_IO_UNITS_RADIANS) {
@@ -113,7 +110,7 @@ static PJ_COORD inv_finalize (PJ *P, PJ_COORD coo) {
         if (P->vgridshift)
             coo = proj_trans (P->vgridshift, PJ_INV, coo); /* Go geometric from orthometric */
         if (coo.lp.lam==HUGE_VAL)
-            return coo;
+            return;
         if (P->hgridshift)
             coo = proj_trans (P->hgridshift, PJ_FWD, coo);
         else if (P->helmert || (P->cart_wgs84 != nullptr && P->cart != nullptr)) {
@@ -123,35 +120,32 @@ static PJ_COORD inv_finalize (PJ *P, PJ_COORD coo) {
             coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
         }
         if (coo.lp.lam==HUGE_VAL)
-            return coo;
+            return;
 
         /* If input latitude was geocentrical, convert back to geocentrical */
         if (P->geoc)
             coo = pj_geocentric_latitude (P, PJ_FWD, coo);
     }
-
-    return coo;
 }
 
-
-static PJ_COORD error_or_coord(PJ *P, PJ_COORD coord, int last_errno) {
-    if (proj_errno(P))
+static inline PJ_COORD error_or_coord(PJ *P, PJ_COORD coord, int last_errno) {
+    if (P->ctx->last_errno)
         return proj_coord_error();
 
-    proj_errno_restore(P, last_errno);
+    P->ctx->last_errno = last_errno;
+
     return coord;
 }
 
-
 PJ_LP pj_inv(PJ_XY xy, PJ *P) {
-    int last_errno;
     PJ_COORD coo = {{0,0,0,0}};
     coo.xy = xy;
 
-    last_errno = proj_errno_reset(P);
+    const int last_errno = P->ctx->last_errno;
+    P->ctx->last_errno = 0;
 
     if (!P->skip_inv_prepare)
-        coo = inv_prepare (P, coo);
+        inv_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
         return proj_coord_error ().lp;
 
@@ -170,7 +164,7 @@ PJ_LP pj_inv(PJ_XY xy, PJ *P) {
         return proj_coord_error ().lp;
 
     if (!P->skip_inv_finalize)
-        coo = inv_finalize (P, coo);
+        inv_finalize (P, coo);
 
     return error_or_coord(P, coo, last_errno).lp;
 }
@@ -178,14 +172,14 @@ PJ_LP pj_inv(PJ_XY xy, PJ *P) {
 
 
 PJ_LPZ pj_inv3d (PJ_XYZ xyz, PJ *P) {
-    int last_errno;
     PJ_COORD coo = {{0,0,0,0}};
     coo.xyz = xyz;
 
-    last_errno = proj_errno_reset(P);
+    const int last_errno = P->ctx->last_errno;
+    P->ctx->last_errno = 0;
 
     if (!P->skip_inv_prepare)
-        coo = inv_prepare (P, coo);
+         inv_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
         return proj_coord_error ().lpz;
 
@@ -204,7 +198,7 @@ PJ_LPZ pj_inv3d (PJ_XYZ xyz, PJ *P) {
         return proj_coord_error ().lpz;
 
     if (!P->skip_inv_finalize)
-        coo = inv_finalize (P, coo);
+        inv_finalize (P, coo);
 
     return error_or_coord(P, coo, last_errno).lpz;
 }
@@ -212,10 +206,12 @@ PJ_LPZ pj_inv3d (PJ_XYZ xyz, PJ *P) {
 
 
 PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
-    int last_errno = proj_errno_reset(P);
+
+    const int last_errno = P->ctx->last_errno;
+    P->ctx->last_errno = 0;
 
     if (!P->skip_inv_prepare)
-        coo = inv_prepare (P, coo);
+        inv_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
         return proj_coord_error ();
 
@@ -234,7 +230,7 @@ PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
         return proj_coord_error ();
 
     if (!P->skip_inv_finalize)
-        coo = inv_finalize (P, coo);
+        inv_finalize (P, coo);
 
     return error_or_coord(P, coo, last_errno);
 }
diff --git a/src/log.cpp b/src/log.cpp
index 4a2cc73d..9c96f53a 100644
--- a/src/log.cpp
+++ b/src/log.cpp
@@ -46,25 +46,36 @@ void pj_stderr_logger( void *app_data, int level, const char *msg )
 }
 
 /************************************************************************/
-/*                               pj_vlog()                              */
+/*                           pj_log_active()                            */
 /************************************************************************/
 
-static
-void pj_vlog( PJ_CONTEXT *ctx, int level, const PJ* P, const char *fmt, va_list args )
-
+bool pj_log_active( PJ_CONTEXT *ctx, int level )
 {
-    char *msg_buf;
     int debug_level = ctx->debug_level;
     int shutup_unless_errno_set = debug_level < 0;
 
     /* For negative debug levels, we first start logging when errno is set */
     if (ctx->last_errno==0 && shutup_unless_errno_set)
-        return;
+        return false;
 
     if (debug_level < 0)
         debug_level = -debug_level;
 
     if( level > debug_level )
+        return false;
+    return true;
+}
+
+/************************************************************************/
+/*                               pj_vlog()                              */
+/************************************************************************/
+
+static
+void pj_vlog( PJ_CONTEXT *ctx, int level, const PJ* P, const char *fmt, va_list args )
+
+{
+    char *msg_buf;
+    if( !pj_log_active(ctx, level) )
         return;
 
     constexpr size_t BUF_SIZE = 100000;
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 582bb3c5..6edb6aec 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -918,6 +918,7 @@ void pj_acquire_lock(void);
 void pj_release_lock(void);
 void pj_cleanup_lock(void);
 
+bool pj_log_active( PJ_CONTEXT *ctx, int level );
 void pj_log( PJ_CONTEXT * ctx, int level, const char *fmt, ... );
 void pj_stderr_logger( void *, int, const char * );
 
diff --git a/src/proj_symbol_rename.h b/src/proj_symbol_rename.h
index 7fbe4242..6f9e8ab6 100644
--- a/src/proj_symbol_rename.h
+++ b/src/proj_symbol_rename.h
@@ -71,6 +71,7 @@
 #define pj_is_latlong internal_pj_is_latlong
 #define pj_latlong_from_proj internal_pj_latlong_from_proj
 #define pj_log internal_pj_log
+#define pj_log_active internal_pj_log_active
 #define pj_malloc internal_pj_malloc
 #define pj_open_lib internal_pj_open_lib
 #define pj_pr_list internal_pj_pr_list