Apply bilinear interpolation for horizontal shift grids
| 2894 | |
| 2895 | // Apply bilinear interpolation for horizontal shift grids |
| 2896 | static PJ_LP pj_hgrid_interpolate(PJ_LP t, const HorizontalShiftGrid *grid, |
| 2897 | bool compensateNTConvention) { |
| 2898 | PJ_LP val, frct; |
| 2899 | ILP indx; |
| 2900 | int in; |
| 2901 | |
| 2902 | const auto &extent = grid->extentAndRes(); |
| 2903 | t.lam /= extent.resX; |
| 2904 | indx.lam = std::isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam)); |
| 2905 | t.phi /= extent.resY; |
| 2906 | indx.phi = std::isnan(t.phi) ? 0 : (pj_int32)lround(floor(t.phi)); |
| 2907 | |
| 2908 | frct.lam = t.lam - indx.lam; |
| 2909 | frct.phi = t.phi - indx.phi; |
| 2910 | val.lam = val.phi = HUGE_VAL; |
| 2911 | if (indx.lam < 0) { |
| 2912 | if (indx.lam == -1 && frct.lam > 1 - 10 * REL_TOLERANCE_HGRIDSHIFT) { |
| 2913 | ++indx.lam; |
| 2914 | frct.lam = 0.; |
| 2915 | } else |
| 2916 | return val; |
| 2917 | } else if ((in = indx.lam + 1) >= grid->width()) { |
| 2918 | if (in == grid->width() && frct.lam < 10 * REL_TOLERANCE_HGRIDSHIFT) { |
| 2919 | --indx.lam; |
| 2920 | frct.lam = 1.; |
| 2921 | } else |
| 2922 | return val; |
| 2923 | } |
| 2924 | if (indx.phi < 0) { |
| 2925 | if (indx.phi == -1 && frct.phi > 1 - 10 * REL_TOLERANCE_HGRIDSHIFT) { |
| 2926 | ++indx.phi; |
| 2927 | frct.phi = 0.; |
| 2928 | } else |
| 2929 | return val; |
| 2930 | } else if ((in = indx.phi + 1) >= grid->height()) { |
| 2931 | if (in == grid->height() && frct.phi < 10 * REL_TOLERANCE_HGRIDSHIFT) { |
| 2932 | --indx.phi; |
| 2933 | frct.phi = 1.; |
| 2934 | } else |
| 2935 | return val; |
| 2936 | } |
| 2937 | |
| 2938 | float f00Lon = 0, f00Lat = 0; |
| 2939 | float f10Lon = 0, f10Lat = 0; |
| 2940 | float f01Lon = 0, f01Lat = 0; |
| 2941 | float f11Lon = 0, f11Lat = 0; |
| 2942 | if (!grid->valueAt(indx.lam, indx.phi, compensateNTConvention, f00Lon, |
| 2943 | f00Lat) || |
| 2944 | !grid->valueAt(indx.lam + 1, indx.phi, compensateNTConvention, f10Lon, |
| 2945 | f10Lat) || |
| 2946 | !grid->valueAt(indx.lam, indx.phi + 1, compensateNTConvention, f01Lon, |
| 2947 | f01Lat) || |
| 2948 | !grid->valueAt(indx.lam + 1, indx.phi + 1, compensateNTConvention, |
| 2949 | f11Lon, f11Lat)) { |
| 2950 | return val; |
| 2951 | } |
| 2952 | |
| 2953 | double m10 = frct.lam; |