summaryrefslogtreecommitdiff
path: root/boost/geometry/srs/projections/proj/healpix.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/geometry/srs/projections/proj/healpix.hpp')
-rw-r--r--boost/geometry/srs/projections/proj/healpix.hpp508
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
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-18 01:52:32 +0000