Apply bilinear interpolation for horizontal shift grids
| 3323 | |
| 3324 | // Apply bilinear interpolation for horizontal shift grids |
| 3325 | static PJ_LP pj_hgrid_interpolate(PJ_LP t, const HorizontalShiftGrid *grid, |
| 3326 | bool compensateNTConvention) { |
| 3327 | PJ_LP val, frct; |
| 3328 | ILP indx; |
| 3329 | int in; |
| 3330 | |
| 3331 | const auto &extent = grid->extentAndRes(); |
| 3332 | t.lam /= extent.resX; |
| 3333 | indx.lam = std::isnan(t.lam) ? 0 : (int32_t)lround(floor(t.lam)); |
| 3334 | t.phi /= extent.resY; |
| 3335 | indx.phi = std::isnan(t.phi) ? 0 : (int32_t)lround(floor(t.phi)); |
| 3336 | |
| 3337 | frct.lam = t.lam - indx.lam; |
| 3338 | frct.phi = t.phi - indx.phi; |
| 3339 | val.lam = val.phi = HUGE_VAL; |
| 3340 | if (indx.lam < 0) { |
| 3341 | if (indx.lam == -1 && frct.lam > 1 - 10 * REL_TOLERANCE_HGRIDSHIFT) { |
| 3342 | ++indx.lam; |
| 3343 | frct.lam = 0.; |
| 3344 | } else |
| 3345 | return val; |
| 3346 | } else if ((in = indx.lam + 1) >= grid->width()) { |
| 3347 | if (in == grid->width() && frct.lam < 10 * REL_TOLERANCE_HGRIDSHIFT) { |
| 3348 | --indx.lam; |
| 3349 | frct.lam = 1.; |
| 3350 | } else |
| 3351 | return val; |
| 3352 | } |
| 3353 | if (indx.phi < 0) { |
| 3354 | if (indx.phi == -1 && frct.phi > 1 - 10 * REL_TOLERANCE_HGRIDSHIFT) { |
| 3355 | ++indx.phi; |
| 3356 | frct.phi = 0.; |
| 3357 | } else |
| 3358 | return val; |
| 3359 | } else if ((in = indx.phi + 1) >= grid->height()) { |
| 3360 | if (in == grid->height() && frct.phi < 10 * REL_TOLERANCE_HGRIDSHIFT) { |
| 3361 | --indx.phi; |
| 3362 | frct.phi = 1.; |
| 3363 | } else |
| 3364 | return val; |
| 3365 | } |
| 3366 | |
| 3367 | float f00Long = 0, f00Lat = 0; |
| 3368 | float f10Long = 0, f10Lat = 0; |
| 3369 | float f01Long = 0, f01Lat = 0; |
| 3370 | float f11Long = 0, f11Lat = 0; |
| 3371 | if (!grid->valueAt(indx.lam, indx.phi, compensateNTConvention, f00Long, |
| 3372 | f00Lat) || |
| 3373 | !grid->valueAt(indx.lam + 1, indx.phi, compensateNTConvention, f10Long, |
| 3374 | f10Lat) || |
| 3375 | !grid->valueAt(indx.lam, indx.phi + 1, compensateNTConvention, f01Long, |
| 3376 | f01Lat) || |
| 3377 | !grid->valueAt(indx.lam + 1, indx.phi + 1, compensateNTConvention, |
| 3378 | f11Long, f11Lat)) { |
| 3379 | return val; |
| 3380 | } |
| 3381 | |
| 3382 | double m10 = frct.lam; |
no test coverage detected