transverse_mercator.cuh

转到此文件的文档。

/*
 * 版权所有 (c) 2023, NVIDIA CORPORATION.
 *
 * 根据 Apache 许可证 2.0 版（“许可证”）获得许可；
 * 除非符合许可证的规定，否则您不得使用此文件。
 * 您可以在以下网址获得许可证副本：
 *
 * https://apache.ac.cn/licenses/LICENSE-2.0
 *
 * 除非适用法律要求或书面同意，否则软件
 * 根据许可证分发，按“原样”提供，
 * 不提供任何明示或默示的担保或条件。
 * 请参阅许可证，了解管理权限和
 * 限制的具体语言。
 */
 
// 此文件中的代码源自 [OSGeo/PROJ 项目](https://github.com/OSGeo/PROJ)
// 并已修改为使用 CUDA 在 GPU 上运行。
//
// PROJ 许可证来自 https://github.com/OSGeo/PROJ/blob/9.2/COPYING:
// 但请注意，此文件源自的文件没有版权声明。
/*
 版权信息可在源代码文件中找到。
 
 --------------
 
 特此授予任何获得此软件及其相关文档文件（“软件”）副本的人，
 免费处理此软件的许可，不受任何限制，包括但不限于
 使用、复制、修改、合并、发布、分发、再许可、
 和/或销售软件副本的权利，并允许接受者按以下条件
 使用本软件：
 上述版权声明和本许可声明应包含在
 
 软件的所有副本或实质部分中。
 本软件按“原样”提供，不提供任何明示
 
 或默示的任何形式的保证，包括但不限于适销性、
 特定用途适用性和非侵权性的保证。在任何情况下，
 作者或版权持有人均不对因软件或
 使用或其他交易引起的任何索赔、损害或其他
 责任，无论是合同行为、侵权行为或其他原因引起的，承担任何责任。
 
 DEALINGS IN THE SOFTWARE.
*/
 
// 以下文本取自 PROJ 文件
// [tmerc.cpp](https://github.com/OSGeo/PROJ/blob/9.2/src/projections/tmerc.cpp)
// 有关原始 PROJ 许可证文本，请参阅 cuspatial 仓库根目录中的文件 LICENSE_PROJ
// license text.
/*****************************************************************************/
//
// 精确横轴墨卡托函数
//
//
// 此文件中的代码主要基于以下程序：
//
// 作者：Knud Poder 和 Karsten Engsager
//
// 基于以下数学理论：R.Koenig 和 K.H. Weise，“Mathematische
// Grundlagen der hoeheren Geodaesie und Kartographie,
// Springer-Verlag, Berlin/Goettingen” 海德堡，1951年。
//
// 经丹麦哥本哈根 Kort og Matrikelstyrelsen (KMS) 参考网络部门许可在此修改和使用
// 
//
/*****************************************************************************/
 
#pragma once
 
#include <cuproj/detail/utility/cuda.hpp>
#include <cuproj/detail/wrap_to_pi.hpp>
#include <cuproj/ellipsoid.hpp>
#include <cuproj/operation/operation.cuh>
#include <cuproj/projection_parameters.hpp>
 
#include <thrust/iterator/transform_iterator.h>
 
#include <assert.h>
 
namespace cuproj {

 
namespace detail {

template <typename T>
inline static CUPROJ_HOST_DEVICE T gatg(T const* p1, int len_p1, T B, T cos_2B, T sin_2B)
{
 T h = 0, h1, h2 = 0;
 
 T const two_cos_2B = 2 * cos_2B;
 T const* p = p1 + len_p1;
 h1 = *--p;
 while (p - p1) {
 h = -h2 + two_cos_2B * h1 + *--p;
 h2 = h1;
 h1 = h;
  }
 return (B + h * sin_2B);
}

template <typename T>
inline static CUPROJ_HOST_DEVICE T clenshaw_complex(
 T const* a, int size, T sin_arg_r, T cos_arg_r, T sinh_arg_i, T cosh_arg_i, T* R, T* I)
{
 T r, i, hr, hr1, hr2, hi, hi1, hi2;
 
 // 参数
 T const* p = a + size;
 r = 2 * cos_arg_r * cosh_arg_i;
 i = -2 * sin_arg_r * sinh_arg_i;
 
 // 求和循环
 hi1 = hr1 = hi = 0;
 hr = *--p;
 for (; a - p;) {
 hr2 = hr1;
 hi2 = hi1;
 hr1 = hr;
 hi1 = hi;
 hr = -hr2 + r * hr1 - i * hi1 + *--p;
 hi = -hi2 + i * hr1 + r * hi1;
  }
 
 r = sin_arg_r * cosh_arg_i;
 i = cos_arg_r * sinh_arg_i;
 *R = r * hr - i * hi;
 *I = r * hi + i * hr;
 return *R;
}

template <typename T>
static CUPROJ_HOST_DEVICE T clenshaw_real(T const* a, int size, T arg_r)
{
 T r, hr, hr1, hr2, cos_arg_r;
 
 T const* p = a + size;
 cos_arg_r = cos(arg_r);
 r = 2 * cos_arg_r;
 
 // 求和循环
 hr1 = 0;
 hr = *--p;
 for (; a - p;) {
 hr2 = hr1;
 hr1 = hr;
 hr = -hr2 + r * hr1 + *--p;
  }
 return sin(arg_r) * hr;
}
} // namespace detail
 
template <typename Coordinate, typename T = typename Coordinate::value_type>
class transverse_mercator : operation<Coordinate> {
 public
 static constexpr int ETMERC_ORDER = 6;

 CUPROJ_HOST_DEVICE transverse_mercator(projection_parameters<T> const& params) : params_(params)
  {
  }

 CUPROJ_HOST_DEVICE Coordinate operator()(Coordinate const& coord, direction dir) const
  {
 if (dir == direction::FORWARD)
 return forward(coord);
 else
 return inverse(coord);
  }
 
 static constexpr T utm_central_meridian = T{500000}; // UTM 区域的假东值中心
 // UTM 区域的假北值中心
 static constexpr T utm_central_parallel(hemisphere h)
  {
 return (h == hemisphere::NORTH) ? T{0} : T{10000000};
  }

 projection_parameters<T> setup(projection_parameters<T> const& input_params)
  {
 params_ = input_params;
 
 // 这样我们就不必在每个地方限定类名。
 auto& params = params_;
 auto& tmerc_params = params_.tmerc_params_;
 auto& ellipsoid = params_.ellipsoid_;
 
 assert(ellipsoid.es > 0);
 
 params.x0 = utm_central_meridian;
 params.y0 = utm_central_parallel(params.utm_hemisphere_);
 
 if (params.utm_zone_ > 0 && params.utm_zone_ <= 60) {
 --params.utm_zone_;
 } else {
 params.utm_zone_ = lround((floor((detail::wrap_to_pi(params.lam0_) + M_PI) * 30. / M_PI)));
 params.utm_zone_ = std::clamp(params.utm_zone_, 0, 59);
    }
 
 params.lam0_ = (params.utm_zone_ + .5) * M_PI / 30. - M_PI;
 params.k0 = T{0.9996};
 params.phi0 = T{0};
 
 // 第三扁率
 T const n = ellipsoid.n;
 T np = n;
 
 // 三角级数系数 GEO <-> GAUSS
 // cgb := 高斯坐标 -> 大地坐标, KW p190 - 191 (61) - (62)
 // cbg := 大地坐标 -> 高斯坐标, KW p186 - 187 (51) - (52)
 // ETMERC_ORDER = 6阶：Engsager 和 Poder：ICC2007
 
 tmerc_params.cgb[0] =
 n *
 (2 + n * (-2 / 3.0 + n * (-2 + n * (116 / 45.0 + n * (26 / 45.0 + n * (-2854 / 675.0))))));
 tmerc_params.cbg[0] =
 n * (-2 + n * (2 / 3.0 +
 n * (4 / 3.0 + n * (-82 / 45.0 + n * (32 / 45.0 + n * (4642 / 4725.0))))));
 np *= n;
 tmerc_params.cgb[1] =
 np * (7 / 3.0 + n * (-8 / 5.0 + n * (-227 / 45.0 + n * (2704 / 315.0 + n * (2323 / 945.0)))));
 tmerc_params.cbg[1] =
 np * (5 / 3.0 + n * (-16 / 15.0 + n * (-13 / 9.0 + n * (904 / 315.0 + n * (-1522 / 945.0)))));
 np *= n;
 // n^5 系数从 1262/105 修正为 -1262/105
 tmerc_params.cgb[2] =
 np * (56 / 15.0 + n * (-136 / 35.0 + n * (-1262 / 105.0 + n * (73814 / 2835.0))));
 tmerc_params.cbg[2] =
 np * (-26 / 15.0 + n * (34 / 21.0 + n * (8 / 5.0 + n * (-12686 / 2835.0))));
 np *= n;
 // n^5 系数从 322/35 修正为 332/35
 tmerc_params.cgb[3] = np * (4279 / 630.0 + n * (-332 / 35.0 + n * (-399572 / 14175.0)));
 tmerc_params.cbg[3] = np * (1237 / 630.0 + n * (-12 / 5.0 + n * (-24832 / 14175.0)));
 np *= n;
 tmerc_params.cgb[4] = np * (4174 / 315.0 + n * (-144838 / 6237.0));
 tmerc_params.cbg[4] = np * (-734 / 315.0 + n * (109598 / 31185.0));
 np *= n;
 tmerc_params.cgb[5] = np * (601676 / 22275.0);
 tmerc_params.cbg[5] = np * (444337 / 155925.0);
 
 // 投影常数
 // 横轴墨卡托 (UTM, ITM 等)
 np = n * n;
 // 标准子午线象限, K&W p.50 (96), p.19 (38b), p.5 (2)
 tmerc_params.Qn = params.k0 / (1 + n) * (1 + np * (1 / 4.0 + np * (1 / 64.0 + np / 256.0)));
 // 三角级数系数
 // utg := 椭球北、东 -> 球北、东, KW p194 (65)
 // gtu := 球北、东 -> 椭球北、东, KW p196 (69)
 tmerc_params.utg[0] =
 n * (-0.5 +
 n * (2 / 3.0 +
 n * (-37 / 96.0 + n * (1 / 360.0 + n * (81 / 512.0 + n * (-96199 / 604800.0))))));
 tmerc_params.gtu[0] =
 n * (0.5 + n * (-2 / 3.0 + n * (5 / 16.0 + n * (41 / 180.0 +
 n * (-127 / 288.0 + n * (7891 / 37800.0))))));
 tmerc_params.utg[1] =
 np * (-1 / 48.0 +
 n * (-1 / 15.0 + n * (437 / 1440.0 + n * (-46 / 105.0 + n * (1118711 / 3870720.0)))));
 tmerc_params.gtu[1] =
 np * (13 / 48.0 +
 n * (-3 / 5.0 + n * (557 / 1440.0 + n * (281 / 630.0 + n * (-1983433 / 1935360.0)))));
 np *= n;
 tmerc_params.utg[2] =
 np * (-17 / 480.0 + n * (37 / 840.0 + n * (209 / 4480.0 + n * (-5569 / 90720.0))));
 tmerc_params.gtu[2] =
 np * (61 / 240.0 + n * (-103 / 140.0 + n * (15061 / 26880.0 + n * (167603 / 181440.0))));
 np *= n;
 tmerc_params.utg[3] = np * (-4397 / 161280.0 + n * (11 / 504.0 + n * (830251 / 7257600.0)));
 tmerc_params.gtu[3] = np * (49561 / 161280.0 + n * (-179 / 168.0 + n * (6601661 / 7257600.0)));
 np *= n;
 tmerc_params.utg[4] = np * (-4583 / 161280.0 + n * (108847 / 3991680.0));
 tmerc_params.gtu[4] = np * (34729 / 80640.0 + n * (-3418889 / 1995840.0));
 np *= n;
 tmerc_params.utg[5] = np * (-20648693 / 638668800.0);
 tmerc_params.gtu[5] = np * (212378941 / 319334400.0);
 
 // 原点纬度的高斯纬度值
 T const Z = detail::gatg(
 tmerc_params.cbg, ETMERC_ORDER, params.phi0, cos(2 * params.phi0), sin(2 * params.phi0));
 
 // 原点北值减去原点纬度的真实北值
 // 即 真实北值 = N - P->Zb
 tmerc_params.Zb =
 -tmerc_params.Qn * (Z + detail::clenshaw_real(tmerc_params.gtu, ETMERC_ORDER, 2 * Z));
 
 return params;
  }
 
 private
 CUPROJ_HOST_DEVICE Coordinate forward(Coordinate const& coord) const
  {
 // 这样我们就不必在每个地方限定类名。
 auto& tmerc_params = this->params_.tmerc_params_;
 auto& ellipsoid = this->params_.ellipsoid_;
 
 // 椭球纬度、经度 -> 高斯纬度、经度
 T Cn =
 detail::gatg(tmerc_params.cbg, ETMERC_ORDER, coord.y, cos(2 * coord.y), sin(2 * coord.y));
 
 // 高斯纬度、经度 -> 复球纬度
 T const sin_Cn = sin(Cn);
 T const cos_Cn = cos(Cn);
 T const sin_Ce = sin(coord.x);
 T const cos_Ce = cos(coord.x);
 
 T const cos_Cn_cos_Ce = cos_Cn * cos_Ce;
 Cn = atan2(sin_Cn, cos_Cn_cos_Ce);
 
 T const inv_denom_tan_Ce = 1. / hypot(sin_Cn, cos_Cn_cos_Ce);
 T const tan_Ce = sin_Ce * cos_Cn * inv_denom_tan_Ce;
 
 // 上述变体：未发现有可测量的加速
 // T const sin_Ce_cos_Cn = sin_Ce*cos_Cn;
 // T const denom = sqrt(1 - sin_Ce_cos_Cn * sin_Ce_cos_Cn);
 // T const tan_Ce = sin_Ce_cos_Cn / denom;
 
 // 复球北、东 -> 椭球标准北、东
 T Ce = asinh(tan_Ce); /* 替换为：Ce = log(tan(FORTPI + Ce*0.5)); */
 
 // 未优化版本：
 // T const sin_arg_r = sin(2*Cn);
 // T const cos_arg_r = cos(2*Cn);
    //
 // 已知：
 // sin(2 * Cn) = 2 sin(Cn) cos(Cn)
 // sin(atan(y)) = y / sqrt(1 + y^2)
 // cos(atan(y)) = 1 / sqrt(1 + y^2)
 // ==> sin(2 * Cn) = 2 tan_Cn / (1 + tan_Cn^2)
    //
 // cos(2 * Cn) = 2cos^2(Cn) - 1
 // = 2 / (1 + tan_Cn^2) - 1
    //
 T const two_inv_denom_tan_Ce = 2 * inv_denom_tan_Ce;
 T const two_inv_denom_tan_Ce_square = two_inv_denom_tan_Ce * inv_denom_tan_Ce;
 T const tmp_r = cos_Cn_cos_Ce * two_inv_denom_tan_Ce_square;
 T const sin_arg_r = sin_Cn * tmp_r;
 T const cos_arg_r = cos_Cn_cos_Ce * tmp_r - 1;
 
 // 未优化版本：
 // T const sinh_arg_i = sinh(2*Ce);
 // T const cosh_arg_i = cosh(2*Ce);
    //
 // 已知
 // sinh(2 * Ce) = 2 sinh(Ce) cosh(Ce)
 // sinh(asinh(y)) = y
 // cosh(asinh(y)) = sqrt(1 + y^2)
 // ==> sinh(2 * Ce) = 2 tan_Ce sqrt(1 + tan_Ce^2)
    //
 // cosh(2 * Ce) = 2cosh^2(Ce) - 1
 // = 2 * (1 + tan_Ce^2) - 1
    //
 // 且 1+tan_Ce^2 = 1 + sin_Ce^2 * cos_Cn^2 / (sin_Cn^2 + cos_Cn^2 *
 // cos_Ce^2) = (sin_Cn^2 + cos_Cn^2 * cos_Ce^2 + sin_Ce^2 * cos_Cn^2) /
 // (sin_Cn^2 + cos_Cn^2 * cos_Ce^2) = 1. / (sin_Cn^2 + cos_Cn^2 * cos_Ce^2)
 // = inv_denom_tan_Ce^2
 
 T const sinh_arg_i = tan_Ce * two_inv_denom_tan_Ce;
 T const cosh_arg_i = two_inv_denom_tan_Ce_square - 1;
 
 T dCn, dCe;
 Cn += detail::clenshaw_complex(
 tmerc_params.gtu, ETMERC_ORDER, sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i, &dCn, &dCe);
 
 Ce += dCe;
 
 CUPROJ_HOST_DEVICE_EXPECTS(fabs(Ce) <= 2.623395162778, // 值来自 PROJ
 "坐标变换超出投影域");
 Coordinate xy{0.0, 0.0};
 xy.y = tmerc_params.Qn * Cn + tmerc_params.Zb; // 北值
 xy.x = tmerc_params.Qn * Ce; // 东值
 
 return xy;
  }

 CUPROJ_HOST_DEVICE Coordinate inverse(Coordinate const& coord) const
  {
 // 这样我们就不必在每个地方限定类名。
 auto& tmerc_params = this->params_.tmerc_params_;
 auto& ellipsoid = this->params_.ellipsoid_;
 
 // 标准化北值、东值
 T Cn = (coord.y - tmerc_params.Zb) / tmerc_params.Qn;
 T Ce = coord.x / tmerc_params.Qn;
 
 CUPROJ_HOST_DEVICE_EXPECTS(fabs(Ce) <= 2.623395162778, // 值来自 PROJ
 "坐标变换超出投影域");
 
 // 标准化北值、东值 -> 复球纬度、经度
 T const sin_arg_r = sin(2 * Cn);
 T const cos_arg_r = cos(2 * Cn);
 
 // T const sinh_arg_i = sinh(2*Ce);
 // T const cosh_arg_i = cosh(2*Ce);
 T const exp_2_Ce = exp(2 * Ce);
 T const half_inv_exp_2_Ce = T{0.5} / exp_2_Ce;
 T const sinh_arg_i = T{0.5} * exp_2_Ce - half_inv_exp_2_Ce;
 T const cosh_arg_i = T{0.5} * exp_2_Ce + half_inv_exp_2_Ce;
 
 T dCn_ignored, dCe;
 Cn += detail::clenshaw_complex(tmerc_params.utg,
 ETMERC_ORDER,
 sin_arg_r,
 cos_arg_r,
 sinh_arg_i,
 cosh_arg_i,
 &dCn_ignored,
 &dCe);
 Ce += dCe;
 
 // 复球纬度 -> 高斯纬度、经度
 T const sin_Cn = sin(Cn);
 T const cos_Cn = cos(Cn);
 
#if 0
 // 未优化版本：
 T sin_Ce, cos_Ce;
 Ce = atan (sinh (Ce)); // 替换为：Ce = 2*(atan(exp(Ce)) - FORTPI);
 sin_Ce = sin (Ce);
 cos_Ce = cos (Ce);
 Ce = atan2 (sin_Ce, cos_Ce*cos_Cn);
 Cn = atan2 (sin_Cn*cos_Ce, hypot (sin_Ce, cos_Ce*cos_Cn));
#else
 // 可以将 Ce = atan2(...) 的两个参数都除以 cos_Ce，
 // 得到：Ce = atan2 (tan_Ce, cos_Cn) = atan2(sinh(Ce), cos_Cn)
    //
 // Cn = atan2(...) 也类似
 // Cn = atan2 (sin_Cn, hypot (sin_Ce, cos_Ce*cos_Cn)/cos_Ce)
 // = atan2 (sin_Cn, hypot (sin_Ce/cos_Ce, cos_Cn))
 // = atan2 (sin_Cn, hypot (tan_Ce, cos_Cn))
 // = atan2 (sin_Cn, hypot (sinhCe, cos_Cn))
 T const sinhCe = sinh(Ce);
 Ce = atan2(sinhCe, cos_Cn);
 T const modulus_Ce = hypot(sinhCe, cos_Cn);
 Cn = atan2(sin_Cn, modulus_Ce);
#endif
 
 // 高斯纬度、经度 -> 椭球纬度、经度
 
 // cos(2*Cn) 和 sin(2*Cn) 计算的优化
 T const tmp = 2 * modulus_Ce / (sinhCe * sinhCe + 1);
 T const sin_2_Cn = sin_Cn * tmp;
 T const cos_2_Cn = tmp * modulus_Ce - 1.;
 // T const cos_2_Cn = cos(2 * Cn);
 // T const sin_2_Cn = sin(2 * Cn);
 
 return Coordinate{Ce, detail::gatg(tmerc_params.cgb, ETMERC_ORDER, Cn, cos_2_Cn, sin_2_Cn)};
  }

 T lam0() const { return this->params_.lam0_; }
 
 private
 projection_parameters<T> params_{};
};

 
} // namespace cuproj