diff options
Diffstat (limited to 'boost/geometry/srs/projections/proj/healpix.hpp')
-rw-r--r-- | boost/geometry/srs/projections/proj/healpix.hpp | 508 |
1 files changed, 253 insertions, 255 deletions
diff --git a/boost/geometry/srs/projections/proj/healpix.hpp b/boost/geometry/srs/projections/proj/healpix.hpp index d2f5a8c081..4615448192 100644 --- a/boost/geometry/srs/projections/proj/healpix.hpp +++ b/boost/geometry/srs/projections/proj/healpix.hpp @@ -1,8 +1,4 @@ -#ifndef BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP -#define BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP - -// Boost.Geometry - extensions-gis-projections (based on PROJ4) -// This file is automatically generated. DO NOT EDIT. +// Boost.Geometry - gis-projections (based on PROJ4) // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. @@ -19,16 +15,17 @@ // PROJ4 is maintained by Frank Warmerdam // PROJ4 is converted to Boost.Geometry by Barend Gehrels -// Last updated version of proj: 4.9.1 +// Last updated version of proj: 5.0.0 // Original copyright notice: // Purpose: Implementation of the HEALPix and rHEALPix projections. -// For background see <http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/rhealpix_dggs.pdf>. +// For background see <http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/rhealpix_dggs.pdf>. // Authors: Alex Raichev (raichev@cs.auckland.ac.nz) -// Michael Speth (spethm@landcareresearch.co.nz) +// Michael Speth (spethm@landcareresearch.co.nz) // Notes: Raichev implemented these projections in Python and -// Speth translated them into C here. +// Speth translated them into C here. + // Copyright (c) 2001, Thomas Flemming, tf@ttqv.com // Permission is hereby granted, free of charge, to any person obtaining a @@ -49,6 +46,9 @@ // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER // DEALINGS IN THE SOFTWARE. +#ifndef BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP +#define BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP + #include <boost/geometry/util/math.hpp> #include <boost/geometry/srs/projections/impl/base_static.hpp> @@ -63,8 +63,8 @@ namespace boost { namespace geometry namespace srs { namespace par4 { - struct healpix {}; - struct rhealpix {}; + struct healpix {}; // HEALPix + struct rhealpix {}; // rHEALPix }} //namespace srs::par4 @@ -74,7 +74,8 @@ namespace projections namespace detail { namespace healpix { - static const double EPS = 1e-15; + /* Fuzz to handle rounding errors: */ + static const double epsilon = 1e-15; template <typename T> struct par_healpix @@ -82,28 +83,36 @@ namespace projections int north_square; int south_square; T qp; - T apa[APA_SIZE]; + detail::apa<T> apa; }; - /* Matrix for counterclockwise rotation by pi/2: */ - /* Matrix for counterclockwise rotation by pi: */ - /* Matrix for counterclockwise rotation by 3*pi/2: */ - /* Identity matrix */ - /* IDENT, R1, R2, R3, R1 inverse, R2 inverse, R3 inverse:*/ - /* Fuzz to handle rounding errors: */ template <typename T> - struct CapMap + struct cap_map { int cn; /* An integer 0--3 indicating the position of the polar cap. */ T x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */ - enum Region {north, south, equatorial} region; + enum region_type {north, south, equatorial} region; }; template <typename T> - struct Point + struct point_xy { T x, y; }; - static double rot[7][2][2] = {{{1, 0},{0, 1}}, {{ 0,-1},{ 1, 0}}, {{-1, 0},{ 0,-1}}, {{ 0, 1},{-1, 0}}, {{ 0, 1},{-1, 0}}, {{-1, 0},{ 0,-1}}, {{ 0,-1},{ 1, 0}}}; + + /* IDENT, R1, R2, R3, R1 inverse, R2 inverse, R3 inverse:*/ + static double rot[7][2][2] = { + /* Identity matrix */ + {{1, 0},{0, 1}}, + /* Matrix for counterclockwise rotation by pi/2: */ + {{ 0,-1},{ 1, 0}}, + /* Matrix for counterclockwise rotation by pi: */ + {{-1, 0},{ 0,-1}}, + /* Matrix for counterclockwise rotation by 3*pi/2: */ + {{ 0, 1},{-1, 0}}, + {{ 0, 1},{-1, 0}}, // 3*pi/2 + {{-1, 0},{ 0,-1}}, // pi + {{ 0,-1},{ 1, 0}} // pi/2 + }; /** * Returns the sign of the double. @@ -149,41 +158,41 @@ namespace projections template <typename T> inline int pnpoly(int nvert, T vert[][2], T const& testx, T const& testy) { - int i, c = 0; + int i; int counter = 0; T xinters; - Point<T> p1, p2; + point_xy<T> p1, p2; + /* Check for boundrary cases */ for (i = 0; i < nvert; i++) { if (testx == vert[i][0] && testy == vert[i][1]) { return 1; } } + p1.x = vert[0][0]; p1.y = vert[0][1]; + for (i = 1; i < nvert; i++) { p2.x = vert[i % nvert][0]; p2.y = vert[i % nvert][1]; - if (testy > (std::min)(p1.y, p2.y)) { - if (testy <= (std::max)(p1.y, p2.y)) { - if (testx <= (std::max)(p1.x, p2.x)) { - if (p1.y != p2.y) { - xinters = (testy-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x; - if (p1.x == p2.x || testx <= xinters) { - counter++; - } - } - } - } + if (testy > (std::min)(p1.y, p2.y) && + testy <= (std::max)(p1.y, p2.y) && + testx <= (std::max)(p1.x, p2.x) && + p1.y != p2.y) + { + xinters = (testy-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x; + if (p1.x == p2.x || testx <= xinters) + counter++; } p1 = p2; } + if (counter % 2 == 0) { return 0; } else { return 1; } - return c; } /** * Return 1 if (x, y) lies in (the interior or boundary of) the image of the @@ -195,45 +204,49 @@ namespace projections template <typename T> inline int in_image(T const& x, T const& y, int proj, int north_square, int south_square) { - static const T ONEPI = detail::ONEPI<T>(); + static const T pi = detail::pi<T>(); + static const T half_pi = detail::half_pi<T>(); + static const T fourth_pi = detail::fourth_pi<T>(); if (proj == 0) { T healpixVertsJit[][2] = { - {-1.0*ONEPI- EPS, ONEPI/4.0}, - {-3.0*ONEPI/4.0, ONEPI/2.0 + EPS}, - {-1.0*ONEPI/2.0, ONEPI/4.0 + EPS}, - {-1.0*ONEPI/4.0, ONEPI/2.0 + EPS}, - {0.0, ONEPI/4.0 + EPS}, - {ONEPI/4.0, ONEPI/2.0 + EPS}, - {ONEPI/2.0, ONEPI/4.0 + EPS}, - {3.0*ONEPI/4.0, ONEPI/2.0 + EPS}, - {ONEPI+ EPS, ONEPI/4.0}, - {ONEPI+ EPS, -1.0*ONEPI/4.0}, - {3.0*ONEPI/4.0, -1.0*ONEPI/2.0 - EPS}, - {ONEPI/2.0, -1.0*ONEPI/4.0 - EPS}, - {ONEPI/4.0, -1.0*ONEPI/2.0 - EPS}, - {0.0, -1.0*ONEPI/4.0 - EPS}, - {-1.0*ONEPI/4.0, -1.0*ONEPI/2.0 - EPS}, - {-1.0*ONEPI/2.0, -1.0*ONEPI/4.0 - EPS}, - {-3.0*ONEPI/4.0, -1.0*ONEPI/2.0 - EPS}, - {-1.0*ONEPI - EPS, -1.0*ONEPI/4.0} + {-pi - epsilon, fourth_pi}, + {-3.0*fourth_pi, half_pi + epsilon}, + {-half_pi, fourth_pi + epsilon}, + {-fourth_pi, half_pi + epsilon}, + {0.0, fourth_pi + epsilon}, + {fourth_pi, half_pi + epsilon}, + {half_pi, fourth_pi + epsilon}, + {3.0*fourth_pi, half_pi + epsilon}, + {pi + epsilon, fourth_pi}, + {pi + epsilon, -fourth_pi}, + {3.0*fourth_pi, -half_pi - epsilon}, + {half_pi, -fourth_pi - epsilon}, + {fourth_pi, -half_pi - epsilon}, + {0.0, -fourth_pi - epsilon}, + {-fourth_pi, -half_pi - epsilon}, + {-half_pi, -fourth_pi - epsilon}, + {-3.0*fourth_pi, -half_pi - epsilon}, + {-pi - epsilon, -fourth_pi} }; return pnpoly((int)sizeof(healpixVertsJit)/ sizeof(healpixVertsJit[0]), healpixVertsJit, x, y); } else { T rhealpixVertsJit[][2] = { - {-1.0*ONEPI - EPS, ONEPI/4.0 + EPS}, - {-1.0*ONEPI + north_square*ONEPI/2.0- EPS, ONEPI/4.0 + EPS}, - {-1.0*ONEPI + north_square*ONEPI/2.0- EPS, 3*ONEPI/4.0 + EPS}, - {-1.0*ONEPI + (north_square + 1.0)*ONEPI/2.0 + EPS, 3*ONEPI/4.0 + EPS}, - {-1.0*ONEPI + (north_square + 1.0)*ONEPI/2.0 + EPS, ONEPI/4.0 + EPS}, - {ONEPI + EPS, ONEPI/4.0 + EPS}, - {ONEPI + EPS, -1.0*ONEPI/4.0 - EPS}, - {-1.0*ONEPI + (south_square + 1.0)*ONEPI/2.0 + EPS, -1.0*ONEPI/4.0 - EPS}, - {-1.0*ONEPI + (south_square + 1.0)*ONEPI/2.0 + EPS, -3.0*ONEPI/4.0 - EPS}, - {-1.0*ONEPI + south_square*ONEPI/2.0 - EPS, -3.0*ONEPI/4.0 - EPS}, - {-1.0*ONEPI + south_square*ONEPI/2.0 - EPS, -1.0*ONEPI/4.0 - EPS}, - {-1.0*ONEPI - EPS, -1.0*ONEPI/4.0 - EPS}}; + {-pi - epsilon, fourth_pi + epsilon}, + {-pi + north_square*half_pi - epsilon, fourth_pi + epsilon}, + {-pi + north_square*half_pi - epsilon, 3.0*fourth_pi + epsilon}, + {-pi + (north_square + 1.0)*half_pi + epsilon, 3.0*fourth_pi + epsilon}, + {-pi + (north_square + 1.0)*half_pi + epsilon, fourth_pi + epsilon}, + {pi + epsilon, fourth_pi + epsilon}, + {pi + epsilon, -fourth_pi - epsilon}, + {-pi + (south_square + 1.0)*half_pi + epsilon, -fourth_pi - epsilon}, + {-pi + (south_square + 1.0)*half_pi + epsilon, -3.0*fourth_pi - epsilon}, + {-pi + south_square*half_pi - epsilon, -3.0*fourth_pi - epsilon}, + {-pi + south_square*half_pi - epsilon, -fourth_pi - epsilon}, + {-pi - epsilon, -fourth_pi - epsilon} + }; + return pnpoly((int)sizeof(rhealpixVertsJit)/ sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y); } @@ -251,10 +264,12 @@ namespace projections T q = pj_qsfn(sin(alpha), par.e, 1.0 - par.es); T qp = proj_parm.qp; T ratio = q/qp; - if (fabsl(ratio) > 1) { + + if (math::abs(ratio) > 1) { /* Rounding error. */ ratio = pj_sign(ratio); } + return asin(ratio); } else { /* Approximation to inverse authalic latitude. */ @@ -268,7 +283,9 @@ namespace projections template <typename T> inline void healpix_sphere(T const& lp_lam, T const& lp_phi, T& xy_x, T& xy_y) { - static const T ONEPI = detail::ONEPI<T>(); + static const T pi = detail::pi<T>(); + static const T half_pi = detail::half_pi<T>(); + static const T fourth_pi = detail::fourth_pi<T>(); T lam = lp_lam; T phi = lp_phi; @@ -277,17 +294,17 @@ namespace projections /* equatorial region */ if ( fabsl(phi) <= phi0) { xy_x = lam; - xy_y = 3.0*ONEPI/8.0*sin(phi); + xy_y = 3.0*pi/8.0*sin(phi); } else { T lamc; - T sigma = sqrt(3.0*(1 - fabsl(sin(phi)))); - T cn = floor(2*lam / ONEPI + 2); + T sigma = sqrt(3.0*(1 - math::abs(sin(phi)))); + T cn = floor(2*lam / pi + 2); if (cn >= 4) { cn = 3; } - lamc = -3*ONEPI/4 + (ONEPI/2)*cn; + lamc = -3*fourth_pi + half_pi*cn; xy_x = lamc + (lam - lamc)*sigma; - xy_y = pj_sign(phi)*ONEPI/4*(2 - sigma); + xy_y = pj_sign(phi)*fourth_pi*(2 - sigma); } return; } @@ -297,28 +314,31 @@ namespace projections template <typename T> inline void healpix_sphere_inverse(T const& xy_x, T const& xy_y, T& lp_lam, T& lp_phi) { - static const T ONEPI = detail::ONEPI<T>(); + static const T pi = detail::pi<T>(); + static const T half_pi = detail::half_pi<T>(); + static const T fourth_pi = detail::fourth_pi<T>(); T x = xy_x; T y = xy_y; - T y0 = ONEPI/4.0; + T y0 = fourth_pi; + /* Equatorial region. */ - if (fabsl(y) <= y0) { + if (math::abs(y) <= y0) { lp_lam = x; - lp_phi = asin(8.0*y/(3.0*ONEPI)); - } else if (fabsl(y) < ONEPI/2.0) { - T cn = floor(2.0*x/ONEPI + 2.0); + lp_phi = asin(8.0*y/(3.0*pi)); + } else if (fabsl(y) < half_pi) { + T cn = floor(2.0*x/pi + 2.0); T xc, tau; if (cn >= 4) { cn = 3; } - xc = -3.0*ONEPI/4.0 + (ONEPI/2.0)*cn; - tau = 2.0 - 4.0*fabsl(y)/ONEPI; + xc = -3.0*fourth_pi + (half_pi)*cn; + tau = 2.0 - 4.0*fabsl(y)/pi; lp_lam = xc + (x - xc)/tau; - lp_phi = pj_sign(y)*asin(1.0 - pow(tau , 2.0)/3.0); + lp_phi = pj_sign(y)*asin(1.0 - math::pow(tau, 2)/3.0); } else { - lp_lam = -1.0*ONEPI; - lp_phi = pj_sign(y)*ONEPI/2.0; + lp_lam = -1.0*pi; + lp_phi = pj_sign(y)*half_pi; } return; } @@ -351,8 +371,8 @@ namespace projections * b is a 2 x 1 matrix. * @param ret holds a*b. **/ - template <typename T> - inline void dot_product(T a[2][2], T b[2], T *ret) + template <typename T1, typename T2> + inline void dot_product(T1 a[2][2], T2 b[2], T2 *ret) { int i, j; int length = 2; @@ -372,89 +392,88 @@ namespace projections * (north_square, south_square)-rHEALPix projection of the unit sphere. **/ template <typename T> - inline CapMap<T> get_cap(T x, T const& y, int north_square, int south_square, + inline cap_map<T> get_cap(T x, T const& y, int north_square, int south_square, int inverse) { - static const T ONEPI = detail::ONEPI<T>(); + static const T pi = detail::pi<T>(); + static const T half_pi = detail::half_pi<T>(); + static const T fourth_pi = detail::fourth_pi<T>(); - CapMap<T> capmap; + cap_map<T> capmap; T c; capmap.x = x; capmap.y = y; if (inverse == 0) { - if (y > ONEPI/4.0) { - capmap.region = CapMap<T>::north; - c = ONEPI/2.0; - } else if (y < -1*ONEPI/4.0) { - capmap.region = CapMap<T>::south; - c = -1*ONEPI/2.0; + if (y > fourth_pi) { + capmap.region = cap_map<T>::north; + c = half_pi; + } else if (y < -fourth_pi) { + capmap.region = cap_map<T>::south; + c = -half_pi; } else { - capmap.region = CapMap<T>::equatorial; + capmap.region = cap_map<T>::equatorial; capmap.cn = 0; return capmap; } /* polar region */ - if (x < -1*ONEPI/2.0) { + if (x < -half_pi) { capmap.cn = 0; - capmap.x = (-1*3.0*ONEPI/4.0); + capmap.x = (-3.0*fourth_pi); capmap.y = c; - } else if (x >= -1*ONEPI/2.0 && x < 0) { + } else if (x >= -half_pi && x < 0) { capmap.cn = 1; - capmap.x = -1*ONEPI/4.0; + capmap.x = -fourth_pi; capmap.y = c; - } else if (x >= 0 && x < ONEPI/2.0) { + } else if (x >= 0 && x < half_pi) { capmap.cn = 2; - capmap.x = ONEPI/4.0; + capmap.x = fourth_pi; capmap.y = c; } else { capmap.cn = 3; - capmap.x = 3.0*ONEPI/4.0; + capmap.x = 3.0*fourth_pi; capmap.y = c; } - return capmap; } else { - T eps; - if (y > ONEPI/4.0) { - capmap.region = CapMap<T>::north; - capmap.x = (-3.0*ONEPI/4.0 + north_square*ONEPI/2.0); - capmap.y = ONEPI/2.0; - x = x - north_square*ONEPI/2.0; - } else if (y < -1*ONEPI/4.0) { - capmap.region = CapMap<T>::south; - capmap.x = (-3.0*ONEPI/4.0 + south_square*ONEPI/2); - capmap.y = -1*ONEPI/2.0; - x = x - south_square*ONEPI/2.0; + if (y > fourth_pi) { + capmap.region = cap_map<T>::north; + capmap.x = (-3.0*fourth_pi + north_square*half_pi); + capmap.y = half_pi; + x = x - north_square*half_pi; + } else if (y < -fourth_pi) { + capmap.region = cap_map<T>::south; + capmap.x = (-3.0*fourth_pi + south_square*pi/2); + capmap.y = -half_pi; + x = x - south_square*half_pi; } else { - capmap.region = CapMap<T>::equatorial; + capmap.region = cap_map<T>::equatorial; capmap.cn = 0; return capmap; } /* Polar Region, find the HEALPix polar cap number that x, y moves to when rHEALPix polar square is disassembled. */ - eps = 1e-15; /* Kludge. Fuzz to avoid some rounding errors. */ - if (capmap.region == CapMap<T>::north) { - if (y >= -1*x - ONEPI/4.0 - eps && y < x + 5.0*ONEPI/4.0 - eps) { + if (capmap.region == cap_map<T>::north) { + if (y >= -x - fourth_pi - epsilon && y < x + 5.0*fourth_pi - epsilon) { capmap.cn = (north_square + 1) % 4; - } else if (y > -1*x -1*ONEPI/4.0 + eps && y >= x + 5.0*ONEPI/4.0 - eps) { + } else if (y > -x -fourth_pi + epsilon && y >= x + 5.0*fourth_pi - epsilon) { capmap.cn = (north_square + 2) % 4; - } else if (y <= -1*x -1*ONEPI/4.0 + eps && y > x + 5.0*ONEPI/4.0 + eps) { + } else if (y <= -x -fourth_pi + epsilon && y > x + 5.0*fourth_pi + epsilon) { capmap.cn = (north_square + 3) % 4; } else { capmap.cn = north_square; } - } else if (capmap.region == CapMap<T>::south) { - if (y <= x + ONEPI/4.0 + eps && y > -1*x - 5.0*ONEPI/4 + eps) { + } else if (capmap.region == cap_map<T>::south) { + if (y <= x + fourth_pi + epsilon && y > -x - 5.0*fourth_pi + epsilon) { capmap.cn = (south_square + 1) % 4; - } else if (y < x + ONEPI/4.0 - eps && y <= -1*x - 5.0*ONEPI/4.0 + eps) { + } else if (y < x + fourth_pi - epsilon && y <= -x - 5.0*fourth_pi + epsilon) { capmap.cn = (south_square + 2) % 4; - } else if (y >= x + ONEPI/4.0 - eps && y < -1*x - 5.0*ONEPI/4.0 - eps) { + } else if (y >= x + fourth_pi - epsilon && y < -x - 5.0*fourth_pi - epsilon) { capmap.cn = (south_square + 3) % 4; } else { capmap.cn = south_square; } } - return capmap; } + return capmap; } /** * Rearrange point (x, y) in the HEALPix projection by @@ -469,95 +488,84 @@ namespace projections inline void combine_caps(T& xy_x, T& xy_y, int north_square, int south_square, int inverse) { - static const T ONEPI = detail::ONEPI<T>(); + static const T half_pi = detail::half_pi<T>(); + static const T fourth_pi = detail::fourth_pi<T>(); T v[2]; - T a[2]; + T c[2]; T vector[2]; T v_min_c[2]; T ret_dot[2]; - CapMap<T> capmap = get_cap(xy_x, xy_y, north_square, south_square, inverse); - if (capmap.region == CapMap<T>::equatorial) { + const double (*tmpRot)[2]; + int pole = 0; + + cap_map<T> capmap = get_cap(xy_x, xy_y, north_square, south_square, inverse); + if (capmap.region == cap_map<T>::equatorial) { xy_x = capmap.x; xy_y = capmap.y; return; } - v[0] = xy_x; - v[1] = xy_y; + + v[0] = xy_x; v[1] = xy_y; + c[0] = capmap.x; c[1] = capmap.y; + if (inverse == 0) { /* Rotate (xy_x, xy_y) about its polar cap tip and then translate it to north_square or south_square. */ - int pole = 0; - T (*tmpRot)[2]; - T c[2] = {capmap.x, capmap.y}; - if (capmap.region == CapMap<T>::north) { + + if (capmap.region == cap_map<T>::north) { pole = north_square; - a[0] = (-3.0*ONEPI/4.0 + pole*ONEPI/2); - a[1] = (ONEPI/2.0 + pole*0); tmpRot = rot[get_rotate_index(capmap.cn - pole)]; - vector_sub(v, c, v_min_c); - dot_product(tmpRot, v_min_c, ret_dot); - vector_add(ret_dot, a, vector); } else { pole = south_square; - a[0] = (-3.0*ONEPI/4.0 + pole*ONEPI/2); - a[1] = (ONEPI/-2.0 + pole*0); tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))]; - vector_sub(v, c, v_min_c); - dot_product(tmpRot, v_min_c, ret_dot); - vector_add(ret_dot, a, vector); } - xy_x = vector[0]; - xy_y = vector[1]; - return; } else { /* Inverse function. Unrotate (xy_x, xy_y) and then translate it back. */ - int pole = 0; - T (*tmpRot)[2]; - T c[2] = {capmap.x, capmap.y}; + /* disassemble */ - if (capmap.region == CapMap<T>::north) { + if (capmap.region == cap_map<T>::north) { pole = north_square; - a[0] = (-3.0*ONEPI/4.0 + capmap.cn*ONEPI/2); - a[1] = (ONEPI/2.0 + capmap.cn*0); tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))]; - vector_sub(v, c, v_min_c); - dot_product(tmpRot, v_min_c, ret_dot); - vector_add(ret_dot, a, vector); } else { pole = south_square; - a[0] = (-3.0*ONEPI/4.0 + capmap.cn*ONEPI/2); - a[1] = (ONEPI/-2.0 + capmap.cn*0); tmpRot = rot[get_rotate_index(capmap.cn - pole)]; - vector_sub(v, c, v_min_c); - dot_product(tmpRot, v_min_c, ret_dot); - vector_add(ret_dot, a, vector); } - xy_x = vector[0]; - xy_y = vector[1]; - return; } + + vector_sub(v, c, v_min_c); + dot_product(tmpRot, v_min_c, ret_dot); + + { + T a[2]; + /* Workaround cppcheck git issue */ + T* pa = a; + // TODO: in proj4 5.0.0 this line is used instead + //pa[0] = -3.0*fourth_pi + ((inverse == 0) ? 0 : capmap.cn) *half_pi; + pa[0] = -3.0*fourth_pi + ((inverse == 0) ? pole : capmap.cn) *half_pi; + pa[1] = half_pi; + vector_add(ret_dot, a, vector); + } + + xy_x = vector[0]; + xy_y = vector[1]; } // template class, using CRTP to implement forward/inverse - template <typename CalculationType, typename Parameters> - struct base_healpix_ellipsoid : public base_t_fi<base_healpix_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters> + template <typename T, typename Parameters> + struct base_healpix_ellipsoid + : public base_t_fi<base_healpix_ellipsoid<T, Parameters>, T, Parameters> { - - typedef CalculationType geographic_type; - typedef CalculationType cartesian_type; - - par_healpix<CalculationType> m_proj_parm; + par_healpix<T> m_proj_parm; inline base_healpix_ellipsoid(const Parameters& par) - : base_t_fi<base_healpix_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters>(*this, par) {} + : base_t_fi<base_healpix_ellipsoid<T, Parameters>, T, Parameters>(*this, par) + {} // FORWARD(e_healpix_forward) ellipsoid // Project coordinates from geographic (lon, lat) to cartesian (x, y) - inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const { lp_lat = auth_lat(this->params(), m_proj_parm, lp_lat, 0); return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); @@ -565,13 +573,13 @@ namespace projections // INVERSE(e_healpix_inverse) ellipsoid // Project coordinates from cartesian (x, y) to geographic (lon, lat) - inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const { /* Check whether (x, y) lies in the HEALPix image. */ if (in_image(xy_x, xy_y, 0, 0, 0) == 0) { lp_lon = HUGE_VAL; lp_lat = HUGE_VAL; - BOOST_THROW_EXCEPTION( projection_exception(-15) ); + BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); } healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); lp_lat = auth_lat(this->params(), m_proj_parm, lp_lat, 1); @@ -585,36 +593,32 @@ namespace projections }; // template class, using CRTP to implement forward/inverse - template <typename CalculationType, typename Parameters> - struct base_healpix_spheroid : public base_t_fi<base_healpix_spheroid<CalculationType, Parameters>, - CalculationType, Parameters> + template <typename T, typename Parameters> + struct base_healpix_spheroid + : public base_t_fi<base_healpix_spheroid<T, Parameters>, T, Parameters> { - - typedef CalculationType geographic_type; - typedef CalculationType cartesian_type; - - par_healpix<CalculationType> m_proj_parm; + par_healpix<T> m_proj_parm; inline base_healpix_spheroid(const Parameters& par) - : base_t_fi<base_healpix_spheroid<CalculationType, Parameters>, - CalculationType, Parameters>(*this, par) {} + : base_t_fi<base_healpix_spheroid<T, Parameters>, T, Parameters>(*this, par) + {} // FORWARD(s_healpix_forward) sphere // Project coordinates from geographic (lon, lat) to cartesian (x, y) - inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const { return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); } // INVERSE(s_healpix_inverse) sphere // Project coordinates from cartesian (x, y) to geographic (lon, lat) - inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const { /* Check whether (x, y) lies in the HEALPix image */ if (in_image(xy_x, xy_y, 0, 0, 0) == 0) { lp_lon = HUGE_VAL; lp_lat = HUGE_VAL; - BOOST_THROW_EXCEPTION( projection_exception(-15) ); + BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); } return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); } @@ -627,23 +631,19 @@ namespace projections }; // template class, using CRTP to implement forward/inverse - template <typename CalculationType, typename Parameters> - struct base_rhealpix_ellipsoid : public base_t_fi<base_rhealpix_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters> + template <typename T, typename Parameters> + struct base_rhealpix_ellipsoid + : public base_t_fi<base_rhealpix_ellipsoid<T, Parameters>, T, Parameters> { - - typedef CalculationType geographic_type; - typedef CalculationType cartesian_type; - - par_healpix<CalculationType> m_proj_parm; + par_healpix<T> m_proj_parm; inline base_rhealpix_ellipsoid(const Parameters& par) - : base_t_fi<base_rhealpix_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters>(*this, par) {} + : base_t_fi<base_rhealpix_ellipsoid<T, Parameters>, T, Parameters>(*this, par) + {} // FORWARD(e_rhealpix_forward) ellipsoid // Project coordinates from geographic (lon, lat) to cartesian (x, y) - inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const { lp_lat = auth_lat(this->params(), m_proj_parm, lp_lat, 0); healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); @@ -652,13 +652,13 @@ namespace projections // INVERSE(e_rhealpix_inverse) ellipsoid // Project coordinates from cartesian (x, y) to geographic (lon, lat) - inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const { /* Check whether (x, y) lies in the rHEALPix image. */ if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) { lp_lon = HUGE_VAL; lp_lat = HUGE_VAL; - BOOST_THROW_EXCEPTION( projection_exception(-15) ); + BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); } combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1); healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); @@ -673,23 +673,19 @@ namespace projections }; // template class, using CRTP to implement forward/inverse - template <typename CalculationType, typename Parameters> - struct base_rhealpix_spheroid : public base_t_fi<base_rhealpix_spheroid<CalculationType, Parameters>, - CalculationType, Parameters> + template <typename T, typename Parameters> + struct base_rhealpix_spheroid + : public base_t_fi<base_rhealpix_spheroid<T, Parameters>, T, Parameters> { - - typedef CalculationType geographic_type; - typedef CalculationType cartesian_type; - - par_healpix<CalculationType> m_proj_parm; + par_healpix<T> m_proj_parm; inline base_rhealpix_spheroid(const Parameters& par) - : base_t_fi<base_rhealpix_spheroid<CalculationType, Parameters>, - CalculationType, Parameters>(*this, par) {} + : base_t_fi<base_rhealpix_spheroid<T, Parameters>, T, Parameters>(*this, par) + {} // FORWARD(s_rhealpix_forward) sphere // Project coordinates from geographic (lon, lat) to cartesian (x, y) - inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const { healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0); @@ -697,13 +693,13 @@ namespace projections // INVERSE(s_rhealpix_inverse) sphere // Project coordinates from cartesian (x, y) to geographic (lon, lat) - inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const { /* Check whether (x, y) lies in the rHEALPix image. */ if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) { lp_lon = HUGE_VAL; lp_lat = HUGE_VAL; - BOOST_THROW_EXCEPTION( projection_exception(-15) ); + BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); } combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1); return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); @@ -720,11 +716,11 @@ namespace projections template <typename Parameters, typename T> inline void setup_healpix(Parameters& par, par_healpix<T>& proj_parm) { - if (par.es) { - pj_authset(par.es, proj_parm.apa); /* For auth_lat(). */ + if (par.es != 0.0) { + proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */ proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */ par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */ - par.ra = 1.0/par.a; + pj_calc_ellipsoid_params(par, par.a, par.es); /* Ensure we have a consistent parameter set */ } else { } } @@ -733,19 +729,21 @@ namespace projections template <typename Parameters, typename T> inline void setup_rhealpix(Parameters& par, par_healpix<T>& proj_parm) { - proj_parm.north_square = pj_param(par.params,"inorth_square").i; - proj_parm.south_square = pj_param(par.params,"isouth_square").i; + proj_parm.north_square = pj_get_param_i(par.params, "north_square"); + proj_parm.south_square = pj_get_param_i(par.params, "south_square"); /* Check for valid north_square and south_square inputs. */ if (proj_parm.north_square < 0 || proj_parm.north_square > 3) { - BOOST_THROW_EXCEPTION( projection_exception(-47) ); + BOOST_THROW_EXCEPTION( projection_exception(error_axis) ); } if (proj_parm.south_square < 0 || proj_parm.south_square > 3) { - BOOST_THROW_EXCEPTION( projection_exception(-47) ); + BOOST_THROW_EXCEPTION( projection_exception(error_axis) ); } - if (par.es) { - pj_authset(par.es, proj_parm.apa); /* For auth_lat(). */ + if (par.es != 0.0) { + proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */ proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */ par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */ + // TODO: why not the same as in healpix? + //pj_calc_ellipsoid_params(par, par.a, par.es); par.ra = 1.0/par.a; } else { } @@ -766,10 +764,10 @@ namespace projections \par Example \image html ex_healpix.gif */ - template <typename CalculationType, typename Parameters> - struct healpix_ellipsoid : public detail::healpix::base_healpix_ellipsoid<CalculationType, Parameters> + template <typename T, typename Parameters> + struct healpix_ellipsoid : public detail::healpix::base_healpix_ellipsoid<T, Parameters> { - inline healpix_ellipsoid(const Parameters& par) : detail::healpix::base_healpix_ellipsoid<CalculationType, Parameters>(par) + inline healpix_ellipsoid(const Parameters& par) : detail::healpix::base_healpix_ellipsoid<T, Parameters>(par) { detail::healpix::setup_healpix(this->m_par, this->m_proj_parm); } @@ -787,10 +785,10 @@ namespace projections \par Example \image html ex_healpix.gif */ - template <typename CalculationType, typename Parameters> - struct healpix_spheroid : public detail::healpix::base_healpix_spheroid<CalculationType, Parameters> + template <typename T, typename Parameters> + struct healpix_spheroid : public detail::healpix::base_healpix_spheroid<T, Parameters> { - inline healpix_spheroid(const Parameters& par) : detail::healpix::base_healpix_spheroid<CalculationType, Parameters>(par) + inline healpix_spheroid(const Parameters& par) : detail::healpix::base_healpix_spheroid<T, Parameters>(par) { detail::healpix::setup_healpix(this->m_par, this->m_proj_parm); } @@ -811,10 +809,10 @@ namespace projections \par Example \image html ex_rhealpix.gif */ - template <typename CalculationType, typename Parameters> - struct rhealpix_ellipsoid : public detail::healpix::base_rhealpix_ellipsoid<CalculationType, Parameters> + template <typename T, typename Parameters> + struct rhealpix_ellipsoid : public detail::healpix::base_rhealpix_ellipsoid<T, Parameters> { - inline rhealpix_ellipsoid(const Parameters& par) : detail::healpix::base_rhealpix_ellipsoid<CalculationType, Parameters>(par) + inline rhealpix_ellipsoid(const Parameters& par) : detail::healpix::base_rhealpix_ellipsoid<T, Parameters>(par) { detail::healpix::setup_rhealpix(this->m_par, this->m_proj_parm); } @@ -835,10 +833,10 @@ namespace projections \par Example \image html ex_rhealpix.gif */ - template <typename CalculationType, typename Parameters> - struct rhealpix_spheroid : public detail::healpix::base_rhealpix_spheroid<CalculationType, Parameters> + template <typename T, typename Parameters> + struct rhealpix_spheroid : public detail::healpix::base_rhealpix_spheroid<T, Parameters> { - inline rhealpix_spheroid(const Parameters& par) : detail::healpix::base_rhealpix_spheroid<CalculationType, Parameters>(par) + inline rhealpix_spheroid(const Parameters& par) : detail::healpix::base_rhealpix_spheroid<T, Parameters>(par) { detail::healpix::setup_rhealpix(this->m_par, this->m_proj_parm); } @@ -853,37 +851,37 @@ namespace projections BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::rhealpix, rhealpix_spheroid, rhealpix_ellipsoid) // Factory entry(s) - template <typename CalculationType, typename Parameters> - class healpix_entry : public detail::factory_entry<CalculationType, Parameters> + template <typename T, typename Parameters> + class healpix_entry : public detail::factory_entry<T, Parameters> { public : - virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const + virtual base_v<T, Parameters>* create_new(const Parameters& par) const { if (par.es) - return new base_v_fi<healpix_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par); + return new base_v_fi<healpix_ellipsoid<T, Parameters>, T, Parameters>(par); else - return new base_v_fi<healpix_spheroid<CalculationType, Parameters>, CalculationType, Parameters>(par); + return new base_v_fi<healpix_spheroid<T, Parameters>, T, Parameters>(par); } }; - template <typename CalculationType, typename Parameters> - class rhealpix_entry : public detail::factory_entry<CalculationType, Parameters> + template <typename T, typename Parameters> + class rhealpix_entry : public detail::factory_entry<T, Parameters> { public : - virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const + virtual base_v<T, Parameters>* create_new(const Parameters& par) const { if (par.es) - return new base_v_fi<rhealpix_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par); + return new base_v_fi<rhealpix_ellipsoid<T, Parameters>, T, Parameters>(par); else - return new base_v_fi<rhealpix_spheroid<CalculationType, Parameters>, CalculationType, Parameters>(par); + return new base_v_fi<rhealpix_spheroid<T, Parameters>, T, Parameters>(par); } }; - template <typename CalculationType, typename Parameters> - inline void healpix_init(detail::base_factory<CalculationType, Parameters>& factory) + template <typename T, typename Parameters> + inline void healpix_init(detail::base_factory<T, Parameters>& factory) { - factory.add_to_factory("healpix", new healpix_entry<CalculationType, Parameters>); - factory.add_to_factory("rhealpix", new rhealpix_entry<CalculationType, Parameters>); + factory.add_to_factory("healpix", new healpix_entry<T, Parameters>); + factory.add_to_factory("rhealpix", new rhealpix_entry<T, Parameters>); } } // namespace detail |