diff --git a/include/proj/internal/lru_cache.hpp b/include/proj/internal/lru_cache.hpp index b7aff6b978..b2e997b36c 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 97ba59996c..8583a6e734 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 3bde3cad5c..c49cbb8f91 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 &data); - std::shared_ptr> get(uint32_t ifdIdx, - uint32_t blockNumber); + const std::vector *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>, lru11::NullLock, - std::unordered_map< - Key, - typename std::list>>>::iterator, - KeyHasher>> - cache_{NUM_BLOCKS_AT_CROSSING_TILES * MAX_SAMPLE_COUNT}; + lru11::Cache, lru11::NullLock> cache_{ + NUM_BLOCKS_AT_CROSSING_TILES * MAX_SAMPLE_COUNT}; }; // --------------------------------------------------------------------------- void BlockCache::insert(uint32_t ifdIdx, uint32_t blockNumber, const std::vector &data) { - cache_.insert(Key(ifdIdx, blockNumber), - std::make_shared>(data)); + cache_.insert((static_cast(ifdIdx) << 32) | blockNumber, data); } // --------------------------------------------------------------------------- -std::shared_ptr> -BlockCache::get(uint32_t ifdIdx, uint32_t blockNumber) { - std::shared_ptr> ret; - cache_.tryGet(Key(ifdIdx, blockNumber), ret); - return ret; +const std::vector *BlockCache::get(uint32_t ifdIdx, + uint32_t blockNumber) { + return cache_.getPtr((static_cast(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 m_buffer{}; unsigned m_blocksPerRow = 0; unsigned m_blocksPerCol = 0; - std::map m_mapOffset{}; - std::map m_mapScale{}; + unsigned m_blocks = 0; + std::vector m_adfOffset{}; + std::vector m_adfScale{}; std::map, 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 float readValue(const std::vector &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(-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(m_width)) && + (m_blockHeight == static_cast(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(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(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 float GTiffGrid::readValue(const std::vector &buffer, uint32_t offsetInBlock, uint16_t sample) const { const auto ptr = reinterpret_cast(buffer.data()); assert(offsetInBlock < buffer.size() / sizeof(T)); const auto val = ptr[offsetInBlock]; - if (!m_hasNodata || static_cast(val) != m_noData) { - double scale = 1; - double offset = 0; - getScaleOffset(scale, offset, sample); + if ((!m_hasNodata || static_cast(val) != m_noData) && + sample < m_adfScale.size()) { + double scale = m_adfScale[sample]; + double offset = m_adfOffset[sample]; return static_cast(val * scale + offset); } else { return static_cast(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 *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 *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 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(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 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 fp, return nullptr; } const int columns = static_cast( - 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( - 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::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::open(PJ_CONTEXT *ctx, return nullptr; } const int columns = static_cast( - 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( - 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(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(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(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(a_valid) + static_cast(b_valid) + + static_cast(c_valid) + static_cast(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(grid_x); int iy = static_cast(grid_y); int ix2 = std::min(ix + 1, grid->width() - 1); diff --git a/src/grids.hpp b/src/grids.hpp index d060fc9534..459bde0718 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 b4a4218994..8925f0e979 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 4a2cc73dd4..9c96f53a06 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 582bb3c5fe..6edb6aec04 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 7fbe4242f4..6f9e8ab694 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