]>
Commit | Line | Data |
---|---|---|
92f5a8d4 | 1 | // Boost.Geometry - gis-projections (based on PROJ4) |
11fdf7f2 TL |
2 | |
3 | // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. | |
4 | ||
92f5a8d4 TL |
5 | // This file was modified by Oracle on 2017, 2018, 2019. |
6 | // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates. | |
11fdf7f2 TL |
7 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle. |
8 | ||
9 | // Use, modification and distribution is subject to the Boost Software License, | |
10 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
11 | // http://www.boost.org/LICENSE_1_0.txt) | |
12 | ||
13 | // This file is converted from PROJ4, http://trac.osgeo.org/proj | |
14 | // PROJ4 is originally written by Gerald Evenden (then of the USGS) | |
15 | // PROJ4 is maintained by Frank Warmerdam | |
16 | // PROJ4 is converted to Boost.Geometry by Barend Gehrels | |
17 | ||
92f5a8d4 | 18 | // Last updated version of proj: 5.0.0 |
11fdf7f2 TL |
19 | |
20 | // Original copyright notice: | |
21 | ||
22 | // Purpose: Implementation of the HEALPix and rHEALPix projections. | |
92f5a8d4 | 23 | // For background see <http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/rhealpix_dggs.pdf>. |
11fdf7f2 | 24 | // Authors: Alex Raichev (raichev@cs.auckland.ac.nz) |
92f5a8d4 | 25 | // Michael Speth (spethm@landcareresearch.co.nz) |
11fdf7f2 | 26 | // Notes: Raichev implemented these projections in Python and |
92f5a8d4 TL |
27 | // Speth translated them into C here. |
28 | ||
11fdf7f2 TL |
29 | // Copyright (c) 2001, Thomas Flemming, tf@ttqv.com |
30 | ||
31 | // Permission is hereby granted, free of charge, to any person obtaining a | |
32 | // copy of this software and associated documentation files (the "Software"), | |
33 | // to deal in the Software without restriction, including without limitation | |
34 | // the rights to use, copy, modify, merge, publish, distribute, sublicense, | |
35 | // and/or sell copies of the Software, and to permit persons to whom the | |
36 | // Software is furnished to do so, subject to the following conditions: | |
37 | ||
38 | // The above copyright notice and this permission notice shall be included | |
39 | // in all copies or substantial portions of the Software. | |
40 | ||
41 | // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS | |
42 | // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | |
43 | // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL | |
44 | // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | |
45 | // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING | |
46 | // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER | |
47 | // DEALINGS IN THE SOFTWARE. | |
48 | ||
92f5a8d4 TL |
49 | #ifndef BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP |
50 | #define BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP | |
11fdf7f2 TL |
51 | |
52 | #include <boost/geometry/srs/projections/impl/base_static.hpp> | |
53 | #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> | |
11fdf7f2 TL |
54 | #include <boost/geometry/srs/projections/impl/factory_entry.hpp> |
55 | #include <boost/geometry/srs/projections/impl/pj_auth.hpp> | |
92f5a8d4 | 56 | #include <boost/geometry/srs/projections/impl/pj_param.hpp> |
11fdf7f2 | 57 | #include <boost/geometry/srs/projections/impl/pj_qsfn.hpp> |
92f5a8d4 | 58 | #include <boost/geometry/srs/projections/impl/projects.hpp> |
11fdf7f2 | 59 | |
92f5a8d4 | 60 | #include <boost/geometry/util/math.hpp> |
11fdf7f2 | 61 | |
92f5a8d4 | 62 | namespace boost { namespace geometry |
11fdf7f2 | 63 | { |
11fdf7f2 TL |
64 | |
65 | namespace projections | |
66 | { | |
67 | #ifndef DOXYGEN_NO_DETAIL | |
68 | namespace detail { namespace healpix | |
69 | { | |
70 | ||
92f5a8d4 TL |
71 | /* Fuzz to handle rounding errors: */ |
72 | static const double epsilon = 1e-15; | |
11fdf7f2 TL |
73 | |
74 | template <typename T> | |
75 | struct par_healpix | |
76 | { | |
92f5a8d4 TL |
77 | T qp; |
78 | detail::apa<T> apa; | |
11fdf7f2 TL |
79 | int north_square; |
80 | int south_square; | |
11fdf7f2 TL |
81 | }; |
82 | ||
11fdf7f2 | 83 | template <typename T> |
92f5a8d4 | 84 | struct cap_map |
11fdf7f2 | 85 | { |
11fdf7f2 | 86 | T x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */ |
92f5a8d4 TL |
87 | int cn; /* An integer 0--3 indicating the position of the polar cap. */ |
88 | enum region_type {north, south, equatorial} region; | |
11fdf7f2 TL |
89 | }; |
90 | template <typename T> | |
92f5a8d4 | 91 | struct point_xy |
11fdf7f2 TL |
92 | { |
93 | T x, y; | |
94 | }; | |
92f5a8d4 TL |
95 | |
96 | /* IDENT, R1, R2, R3, R1 inverse, R2 inverse, R3 inverse:*/ | |
97 | static double rot[7][2][2] = { | |
98 | /* Identity matrix */ | |
99 | {{1, 0},{0, 1}}, | |
100 | /* Matrix for counterclockwise rotation by pi/2: */ | |
101 | {{ 0,-1},{ 1, 0}}, | |
102 | /* Matrix for counterclockwise rotation by pi: */ | |
103 | {{-1, 0},{ 0,-1}}, | |
104 | /* Matrix for counterclockwise rotation by 3*pi/2: */ | |
105 | {{ 0, 1},{-1, 0}}, | |
106 | {{ 0, 1},{-1, 0}}, // 3*pi/2 | |
107 | {{-1, 0},{ 0,-1}}, // pi | |
108 | {{ 0,-1},{ 1, 0}} // pi/2 | |
109 | }; | |
11fdf7f2 TL |
110 | |
111 | /** | |
112 | * Returns the sign of the double. | |
113 | * @param v the parameter whose sign is returned. | |
114 | * @return 1 for positive number, -1 for negative, and 0 for zero. | |
115 | **/ | |
116 | template <typename T> | |
117 | inline T pj_sign (T const& v) | |
118 | { | |
119 | return v > 0 ? 1 : (v < 0 ? -1 : 0); | |
120 | } | |
121 | /** | |
122 | * Return the index of the matrix in {{{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}}}. | |
123 | * @param index ranges from -3 to 3. | |
124 | */ | |
125 | inline int get_rotate_index(int index) | |
126 | { | |
127 | switch(index) { | |
128 | case 0: | |
129 | return 0; | |
130 | case 1: | |
131 | return 1; | |
132 | case 2: | |
133 | return 2; | |
134 | case 3: | |
135 | return 3; | |
136 | case -1: | |
137 | return 4; | |
138 | case -2: | |
139 | return 5; | |
140 | case -3: | |
141 | return 6; | |
142 | } | |
143 | return 0; | |
144 | } | |
145 | /** | |
146 | * Return 1 if point (testx, testy) lies in the interior of the polygon | |
147 | * determined by the vertices in vert, and return 0 otherwise. | |
148 | * See http://paulbourke.net/geometry/polygonmesh/ for more details. | |
149 | * @param nvert the number of vertices in the polygon. | |
150 | * @param vert the (x, y)-coordinates of the polygon's vertices | |
151 | **/ | |
152 | template <typename T> | |
153 | inline int pnpoly(int nvert, T vert[][2], T const& testx, T const& testy) | |
154 | { | |
92f5a8d4 | 155 | int i; |
11fdf7f2 TL |
156 | int counter = 0; |
157 | T xinters; | |
92f5a8d4 TL |
158 | point_xy<T> p1, p2; |
159 | ||
11fdf7f2 TL |
160 | /* Check for boundrary cases */ |
161 | for (i = 0; i < nvert; i++) { | |
162 | if (testx == vert[i][0] && testy == vert[i][1]) { | |
163 | return 1; | |
164 | } | |
165 | } | |
92f5a8d4 | 166 | |
11fdf7f2 TL |
167 | p1.x = vert[0][0]; |
168 | p1.y = vert[0][1]; | |
92f5a8d4 | 169 | |
11fdf7f2 TL |
170 | for (i = 1; i < nvert; i++) { |
171 | p2.x = vert[i % nvert][0]; | |
172 | p2.y = vert[i % nvert][1]; | |
92f5a8d4 TL |
173 | if (testy > (std::min)(p1.y, p2.y) && |
174 | testy <= (std::max)(p1.y, p2.y) && | |
175 | testx <= (std::max)(p1.x, p2.x) && | |
176 | p1.y != p2.y) | |
177 | { | |
178 | xinters = (testy-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x; | |
179 | if (p1.x == p2.x || testx <= xinters) | |
180 | counter++; | |
11fdf7f2 TL |
181 | } |
182 | p1 = p2; | |
183 | } | |
92f5a8d4 | 184 | |
11fdf7f2 TL |
185 | if (counter % 2 == 0) { |
186 | return 0; | |
187 | } else { | |
188 | return 1; | |
189 | } | |
11fdf7f2 TL |
190 | } |
191 | /** | |
192 | * Return 1 if (x, y) lies in (the interior or boundary of) the image of the | |
193 | * HEALPix projection (in case proj=0) or in the image the rHEALPix projection | |
194 | * (in case proj=1), and return 0 otherwise. | |
195 | * @param north_square the position of the north polar square (rHEALPix only) | |
196 | * @param south_square the position of the south polar square (rHEALPix only) | |
197 | **/ | |
198 | template <typename T> | |
199 | inline int in_image(T const& x, T const& y, int proj, int north_square, int south_square) | |
200 | { | |
92f5a8d4 TL |
201 | static const T pi = detail::pi<T>(); |
202 | static const T half_pi = detail::half_pi<T>(); | |
203 | static const T fourth_pi = detail::fourth_pi<T>(); | |
11fdf7f2 TL |
204 | |
205 | if (proj == 0) { | |
206 | T healpixVertsJit[][2] = { | |
92f5a8d4 TL |
207 | {-pi - epsilon, fourth_pi}, |
208 | {-3.0*fourth_pi, half_pi + epsilon}, | |
209 | {-half_pi, fourth_pi + epsilon}, | |
210 | {-fourth_pi, half_pi + epsilon}, | |
211 | {0.0, fourth_pi + epsilon}, | |
212 | {fourth_pi, half_pi + epsilon}, | |
213 | {half_pi, fourth_pi + epsilon}, | |
214 | {3.0*fourth_pi, half_pi + epsilon}, | |
215 | {pi + epsilon, fourth_pi}, | |
216 | {pi + epsilon, -fourth_pi}, | |
217 | {3.0*fourth_pi, -half_pi - epsilon}, | |
218 | {half_pi, -fourth_pi - epsilon}, | |
219 | {fourth_pi, -half_pi - epsilon}, | |
220 | {0.0, -fourth_pi - epsilon}, | |
221 | {-fourth_pi, -half_pi - epsilon}, | |
222 | {-half_pi, -fourth_pi - epsilon}, | |
223 | {-3.0*fourth_pi, -half_pi - epsilon}, | |
224 | {-pi - epsilon, -fourth_pi} | |
11fdf7f2 TL |
225 | }; |
226 | return pnpoly((int)sizeof(healpixVertsJit)/ | |
227 | sizeof(healpixVertsJit[0]), healpixVertsJit, x, y); | |
228 | } else { | |
229 | T rhealpixVertsJit[][2] = { | |
92f5a8d4 TL |
230 | {-pi - epsilon, fourth_pi + epsilon}, |
231 | {-pi + north_square*half_pi - epsilon, fourth_pi + epsilon}, | |
232 | {-pi + north_square*half_pi - epsilon, 3.0*fourth_pi + epsilon}, | |
233 | {-pi + (north_square + 1.0)*half_pi + epsilon, 3.0*fourth_pi + epsilon}, | |
234 | {-pi + (north_square + 1.0)*half_pi + epsilon, fourth_pi + epsilon}, | |
235 | {pi + epsilon, fourth_pi + epsilon}, | |
236 | {pi + epsilon, -fourth_pi - epsilon}, | |
237 | {-pi + (south_square + 1.0)*half_pi + epsilon, -fourth_pi - epsilon}, | |
238 | {-pi + (south_square + 1.0)*half_pi + epsilon, -3.0*fourth_pi - epsilon}, | |
239 | {-pi + south_square*half_pi - epsilon, -3.0*fourth_pi - epsilon}, | |
240 | {-pi + south_square*half_pi - epsilon, -fourth_pi - epsilon}, | |
241 | {-pi - epsilon, -fourth_pi - epsilon} | |
242 | }; | |
243 | ||
11fdf7f2 TL |
244 | return pnpoly((int)sizeof(rhealpixVertsJit)/ |
245 | sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y); | |
246 | } | |
247 | } | |
248 | /** | |
249 | * Return the authalic latitude of latitude alpha (if inverse=0) or | |
250 | * return the approximate latitude of authalic latitude alpha (if inverse=1). | |
251 | * P contains the relavent ellipsoid parameters. | |
252 | **/ | |
253 | template <typename Parameters, typename T> | |
254 | inline T auth_lat(const Parameters& par, const par_healpix<T>& proj_parm, T const& alpha, int inverse) | |
255 | { | |
256 | if (inverse == 0) { | |
257 | /* Authalic latitude. */ | |
258 | T q = pj_qsfn(sin(alpha), par.e, 1.0 - par.es); | |
259 | T qp = proj_parm.qp; | |
260 | T ratio = q/qp; | |
92f5a8d4 TL |
261 | |
262 | if (math::abs(ratio) > 1) { | |
11fdf7f2 TL |
263 | /* Rounding error. */ |
264 | ratio = pj_sign(ratio); | |
265 | } | |
92f5a8d4 | 266 | |
11fdf7f2 TL |
267 | return asin(ratio); |
268 | } else { | |
269 | /* Approximation to inverse authalic latitude. */ | |
270 | return pj_authlat(alpha, proj_parm.apa); | |
271 | } | |
272 | } | |
273 | /** | |
274 | * Return the HEALPix projection of the longitude-latitude point lp on | |
275 | * the unit sphere. | |
276 | **/ | |
277 | template <typename T> | |
278 | inline void healpix_sphere(T const& lp_lam, T const& lp_phi, T& xy_x, T& xy_y) | |
279 | { | |
92f5a8d4 TL |
280 | static const T pi = detail::pi<T>(); |
281 | static const T half_pi = detail::half_pi<T>(); | |
282 | static const T fourth_pi = detail::fourth_pi<T>(); | |
11fdf7f2 TL |
283 | |
284 | T lam = lp_lam; | |
285 | T phi = lp_phi; | |
286 | T phi0 = asin(T(2.0)/T(3.0)); | |
287 | ||
288 | /* equatorial region */ | |
289 | if ( fabsl(phi) <= phi0) { | |
290 | xy_x = lam; | |
92f5a8d4 | 291 | xy_y = 3.0*pi/8.0*sin(phi); |
11fdf7f2 TL |
292 | } else { |
293 | T lamc; | |
92f5a8d4 TL |
294 | T sigma = sqrt(3.0*(1 - math::abs(sin(phi)))); |
295 | T cn = floor(2*lam / pi + 2); | |
11fdf7f2 TL |
296 | if (cn >= 4) { |
297 | cn = 3; | |
298 | } | |
92f5a8d4 | 299 | lamc = -3*fourth_pi + half_pi*cn; |
11fdf7f2 | 300 | xy_x = lamc + (lam - lamc)*sigma; |
92f5a8d4 | 301 | xy_y = pj_sign(phi)*fourth_pi*(2 - sigma); |
11fdf7f2 TL |
302 | } |
303 | return; | |
304 | } | |
305 | /** | |
306 | * Return the inverse of healpix_sphere(). | |
307 | **/ | |
308 | template <typename T> | |
309 | inline void healpix_sphere_inverse(T const& xy_x, T const& xy_y, T& lp_lam, T& lp_phi) | |
310 | { | |
92f5a8d4 TL |
311 | static const T pi = detail::pi<T>(); |
312 | static const T half_pi = detail::half_pi<T>(); | |
313 | static const T fourth_pi = detail::fourth_pi<T>(); | |
11fdf7f2 TL |
314 | |
315 | T x = xy_x; | |
316 | T y = xy_y; | |
92f5a8d4 TL |
317 | T y0 = fourth_pi; |
318 | ||
11fdf7f2 | 319 | /* Equatorial region. */ |
92f5a8d4 | 320 | if (math::abs(y) <= y0) { |
11fdf7f2 | 321 | lp_lam = x; |
92f5a8d4 TL |
322 | lp_phi = asin(8.0*y/(3.0*pi)); |
323 | } else if (fabsl(y) < half_pi) { | |
324 | T cn = floor(2.0*x/pi + 2.0); | |
11fdf7f2 TL |
325 | T xc, tau; |
326 | if (cn >= 4) { | |
327 | cn = 3; | |
328 | } | |
92f5a8d4 TL |
329 | xc = -3.0*fourth_pi + (half_pi)*cn; |
330 | tau = 2.0 - 4.0*fabsl(y)/pi; | |
11fdf7f2 | 331 | lp_lam = xc + (x - xc)/tau; |
92f5a8d4 | 332 | lp_phi = pj_sign(y)*asin(1.0 - math::pow(tau, 2)/3.0); |
11fdf7f2 | 333 | } else { |
92f5a8d4 TL |
334 | lp_lam = -1.0*pi; |
335 | lp_phi = pj_sign(y)*half_pi; | |
11fdf7f2 TL |
336 | } |
337 | return; | |
338 | } | |
339 | /** | |
340 | * Return the vector sum a + b, where a and b are 2-dimensional vectors. | |
341 | * @param ret holds a + b. | |
342 | **/ | |
343 | template <typename T> | |
92f5a8d4 | 344 | inline void vector_add(const T a[2], const T b[2], T ret[2]) |
11fdf7f2 TL |
345 | { |
346 | int i; | |
347 | for(i = 0; i < 2; i++) { | |
348 | ret[i] = a[i] + b[i]; | |
349 | } | |
350 | } | |
351 | /** | |
352 | * Return the vector difference a - b, where a and b are 2-dimensional vectors. | |
353 | * @param ret holds a - b. | |
354 | **/ | |
355 | template <typename T> | |
92f5a8d4 | 356 | inline void vector_sub(const T a[2], const T b[2], T ret[2]) |
11fdf7f2 TL |
357 | { |
358 | int i; | |
359 | for(i = 0; i < 2; i++) { | |
360 | ret[i] = a[i] - b[i]; | |
361 | } | |
362 | } | |
363 | /** | |
364 | * Return the 2 x 1 matrix product a*b, where a is a 2 x 2 matrix and | |
365 | * b is a 2 x 1 matrix. | |
366 | * @param ret holds a*b. | |
367 | **/ | |
92f5a8d4 TL |
368 | template <typename T1, typename T2> |
369 | inline void dot_product(const T1 a[2][2], const T2 b[2], T2 ret[2]) | |
11fdf7f2 TL |
370 | { |
371 | int i, j; | |
372 | int length = 2; | |
373 | for(i = 0; i < length; i++) { | |
374 | ret[i] = 0; | |
375 | for(j = 0; j < length; j++) { | |
376 | ret[i] += a[i][j]*b[j]; | |
377 | } | |
378 | } | |
379 | } | |
380 | /** | |
381 | * Return the number of the polar cap, the pole point coordinates, and | |
382 | * the region that (x, y) lies in. | |
383 | * If inverse=0, then assume (x,y) lies in the image of the HEALPix | |
384 | * projection of the unit sphere. | |
385 | * If inverse=1, then assume (x,y) lies in the image of the | |
386 | * (north_square, south_square)-rHEALPix projection of the unit sphere. | |
387 | **/ | |
388 | template <typename T> | |
92f5a8d4 | 389 | inline cap_map<T> get_cap(T x, T const& y, int north_square, int south_square, |
11fdf7f2 TL |
390 | int inverse) |
391 | { | |
92f5a8d4 TL |
392 | static const T pi = detail::pi<T>(); |
393 | static const T half_pi = detail::half_pi<T>(); | |
394 | static const T fourth_pi = detail::fourth_pi<T>(); | |
11fdf7f2 | 395 | |
92f5a8d4 | 396 | cap_map<T> capmap; |
11fdf7f2 TL |
397 | T c; |
398 | capmap.x = x; | |
399 | capmap.y = y; | |
400 | if (inverse == 0) { | |
92f5a8d4 TL |
401 | if (y > fourth_pi) { |
402 | capmap.region = cap_map<T>::north; | |
403 | c = half_pi; | |
404 | } else if (y < -fourth_pi) { | |
405 | capmap.region = cap_map<T>::south; | |
406 | c = -half_pi; | |
11fdf7f2 | 407 | } else { |
92f5a8d4 | 408 | capmap.region = cap_map<T>::equatorial; |
11fdf7f2 TL |
409 | capmap.cn = 0; |
410 | return capmap; | |
411 | } | |
412 | /* polar region */ | |
92f5a8d4 | 413 | if (x < -half_pi) { |
11fdf7f2 | 414 | capmap.cn = 0; |
92f5a8d4 | 415 | capmap.x = (-3.0*fourth_pi); |
11fdf7f2 | 416 | capmap.y = c; |
92f5a8d4 | 417 | } else if (x >= -half_pi && x < 0) { |
11fdf7f2 | 418 | capmap.cn = 1; |
92f5a8d4 | 419 | capmap.x = -fourth_pi; |
11fdf7f2 | 420 | capmap.y = c; |
92f5a8d4 | 421 | } else if (x >= 0 && x < half_pi) { |
11fdf7f2 | 422 | capmap.cn = 2; |
92f5a8d4 | 423 | capmap.x = fourth_pi; |
11fdf7f2 TL |
424 | capmap.y = c; |
425 | } else { | |
426 | capmap.cn = 3; | |
92f5a8d4 | 427 | capmap.x = 3.0*fourth_pi; |
11fdf7f2 TL |
428 | capmap.y = c; |
429 | } | |
11fdf7f2 | 430 | } else { |
92f5a8d4 TL |
431 | if (y > fourth_pi) { |
432 | capmap.region = cap_map<T>::north; | |
433 | capmap.x = (-3.0*fourth_pi + north_square*half_pi); | |
434 | capmap.y = half_pi; | |
435 | x = x - north_square*half_pi; | |
436 | } else if (y < -fourth_pi) { | |
437 | capmap.region = cap_map<T>::south; | |
438 | capmap.x = (-3.0*fourth_pi + south_square*pi/2); | |
439 | capmap.y = -half_pi; | |
440 | x = x - south_square*half_pi; | |
11fdf7f2 | 441 | } else { |
92f5a8d4 | 442 | capmap.region = cap_map<T>::equatorial; |
11fdf7f2 TL |
443 | capmap.cn = 0; |
444 | return capmap; | |
445 | } | |
446 | /* Polar Region, find the HEALPix polar cap number that | |
447 | x, y moves to when rHEALPix polar square is disassembled. */ | |
92f5a8d4 TL |
448 | if (capmap.region == cap_map<T>::north) { |
449 | if (y >= -x - fourth_pi - epsilon && y < x + 5.0*fourth_pi - epsilon) { | |
11fdf7f2 | 450 | capmap.cn = (north_square + 1) % 4; |
92f5a8d4 | 451 | } else if (y > -x -fourth_pi + epsilon && y >= x + 5.0*fourth_pi - epsilon) { |
11fdf7f2 | 452 | capmap.cn = (north_square + 2) % 4; |
92f5a8d4 | 453 | } else if (y <= -x -fourth_pi + epsilon && y > x + 5.0*fourth_pi + epsilon) { |
11fdf7f2 TL |
454 | capmap.cn = (north_square + 3) % 4; |
455 | } else { | |
456 | capmap.cn = north_square; | |
457 | } | |
92f5a8d4 TL |
458 | } else if (capmap.region == cap_map<T>::south) { |
459 | if (y <= x + fourth_pi + epsilon && y > -x - 5.0*fourth_pi + epsilon) { | |
11fdf7f2 | 460 | capmap.cn = (south_square + 1) % 4; |
92f5a8d4 | 461 | } else if (y < x + fourth_pi - epsilon && y <= -x - 5.0*fourth_pi + epsilon) { |
11fdf7f2 | 462 | capmap.cn = (south_square + 2) % 4; |
92f5a8d4 | 463 | } else if (y >= x + fourth_pi - epsilon && y < -x - 5.0*fourth_pi - epsilon) { |
11fdf7f2 TL |
464 | capmap.cn = (south_square + 3) % 4; |
465 | } else { | |
466 | capmap.cn = south_square; | |
467 | } | |
468 | } | |
11fdf7f2 | 469 | } |
92f5a8d4 | 470 | return capmap; |
11fdf7f2 TL |
471 | } |
472 | /** | |
473 | * Rearrange point (x, y) in the HEALPix projection by | |
474 | * combining the polar caps into two polar squares. | |
475 | * Put the north polar square in position north_square and | |
476 | * the south polar square in position south_square. | |
477 | * If inverse=1, then uncombine the polar caps. | |
478 | * @param north_square integer between 0 and 3. | |
479 | * @param south_square integer between 0 and 3. | |
480 | **/ | |
481 | template <typename T> | |
482 | inline void combine_caps(T& xy_x, T& xy_y, int north_square, int south_square, | |
483 | int inverse) | |
484 | { | |
92f5a8d4 TL |
485 | static const T half_pi = detail::half_pi<T>(); |
486 | static const T fourth_pi = detail::fourth_pi<T>(); | |
11fdf7f2 TL |
487 | |
488 | T v[2]; | |
92f5a8d4 | 489 | T c[2]; |
11fdf7f2 TL |
490 | T vector[2]; |
491 | T v_min_c[2]; | |
492 | T ret_dot[2]; | |
92f5a8d4 TL |
493 | const double (*tmpRot)[2]; |
494 | int pole = 0; | |
495 | ||
496 | cap_map<T> capmap = get_cap(xy_x, xy_y, north_square, south_square, inverse); | |
497 | if (capmap.region == cap_map<T>::equatorial) { | |
11fdf7f2 TL |
498 | xy_x = capmap.x; |
499 | xy_y = capmap.y; | |
500 | return; | |
501 | } | |
92f5a8d4 TL |
502 | |
503 | v[0] = xy_x; v[1] = xy_y; | |
504 | c[0] = capmap.x; c[1] = capmap.y; | |
505 | ||
11fdf7f2 TL |
506 | if (inverse == 0) { |
507 | /* Rotate (xy_x, xy_y) about its polar cap tip and then translate it to | |
508 | north_square or south_square. */ | |
92f5a8d4 TL |
509 | |
510 | if (capmap.region == cap_map<T>::north) { | |
11fdf7f2 | 511 | pole = north_square; |
11fdf7f2 | 512 | tmpRot = rot[get_rotate_index(capmap.cn - pole)]; |
11fdf7f2 TL |
513 | } else { |
514 | pole = south_square; | |
11fdf7f2 | 515 | tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))]; |
11fdf7f2 | 516 | } |
11fdf7f2 TL |
517 | } else { |
518 | /* Inverse function. | |
519 | Unrotate (xy_x, xy_y) and then translate it back. */ | |
92f5a8d4 | 520 | |
11fdf7f2 | 521 | /* disassemble */ |
92f5a8d4 | 522 | if (capmap.region == cap_map<T>::north) { |
11fdf7f2 | 523 | pole = north_square; |
11fdf7f2 | 524 | tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))]; |
11fdf7f2 TL |
525 | } else { |
526 | pole = south_square; | |
11fdf7f2 | 527 | tmpRot = rot[get_rotate_index(capmap.cn - pole)]; |
11fdf7f2 | 528 | } |
11fdf7f2 | 529 | } |
11fdf7f2 | 530 | |
92f5a8d4 TL |
531 | vector_sub(v, c, v_min_c); |
532 | dot_product(tmpRot, v_min_c, ret_dot); | |
11fdf7f2 | 533 | |
92f5a8d4 TL |
534 | { |
535 | T a[2]; | |
536 | /* Workaround cppcheck git issue */ | |
537 | T* pa = a; | |
538 | // TODO: in proj4 5.0.0 this line is used instead | |
539 | //pa[0] = -3.0*fourth_pi + ((inverse == 0) ? 0 : capmap.cn) *half_pi; | |
540 | pa[0] = -3.0*fourth_pi + ((inverse == 0) ? pole : capmap.cn) *half_pi; | |
541 | pa[1] = half_pi; | |
542 | vector_add(ret_dot, a, vector); | |
543 | } | |
544 | ||
545 | xy_x = vector[0]; | |
546 | xy_y = vector[1]; | |
547 | } | |
11fdf7f2 | 548 | |
92f5a8d4 TL |
549 | template <typename T, typename Parameters> |
550 | struct base_healpix_ellipsoid | |
551 | { | |
552 | par_healpix<T> m_proj_parm; | |
11fdf7f2 TL |
553 | |
554 | // FORWARD(e_healpix_forward) ellipsoid | |
555 | // Project coordinates from geographic (lon, lat) to cartesian (x, y) | |
92f5a8d4 | 556 | inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const |
11fdf7f2 | 557 | { |
92f5a8d4 | 558 | lp_lat = auth_lat(par, m_proj_parm, lp_lat, 0); |
11fdf7f2 TL |
559 | return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); |
560 | } | |
561 | ||
562 | // INVERSE(e_healpix_inverse) ellipsoid | |
563 | // Project coordinates from cartesian (x, y) to geographic (lon, lat) | |
92f5a8d4 | 564 | inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const |
11fdf7f2 TL |
565 | { |
566 | /* Check whether (x, y) lies in the HEALPix image. */ | |
567 | if (in_image(xy_x, xy_y, 0, 0, 0) == 0) { | |
568 | lp_lon = HUGE_VAL; | |
569 | lp_lat = HUGE_VAL; | |
92f5a8d4 | 570 | BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); |
11fdf7f2 TL |
571 | } |
572 | healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); | |
92f5a8d4 | 573 | lp_lat = auth_lat(par, m_proj_parm, lp_lat, 1); |
11fdf7f2 TL |
574 | } |
575 | ||
576 | static inline std::string get_name() | |
577 | { | |
578 | return "healpix_ellipsoid"; | |
579 | } | |
580 | ||
581 | }; | |
582 | ||
92f5a8d4 TL |
583 | template <typename T, typename Parameters> |
584 | struct base_healpix_spheroid | |
11fdf7f2 | 585 | { |
92f5a8d4 | 586 | par_healpix<T> m_proj_parm; |
11fdf7f2 TL |
587 | |
588 | // FORWARD(s_healpix_forward) sphere | |
589 | // Project coordinates from geographic (lon, lat) to cartesian (x, y) | |
92f5a8d4 | 590 | inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const |
11fdf7f2 TL |
591 | { |
592 | return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); | |
593 | } | |
594 | ||
595 | // INVERSE(s_healpix_inverse) sphere | |
596 | // Project coordinates from cartesian (x, y) to geographic (lon, lat) | |
92f5a8d4 | 597 | inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const |
11fdf7f2 TL |
598 | { |
599 | /* Check whether (x, y) lies in the HEALPix image */ | |
600 | if (in_image(xy_x, xy_y, 0, 0, 0) == 0) { | |
601 | lp_lon = HUGE_VAL; | |
602 | lp_lat = HUGE_VAL; | |
92f5a8d4 | 603 | BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); |
11fdf7f2 TL |
604 | } |
605 | return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); | |
606 | } | |
607 | ||
608 | static inline std::string get_name() | |
609 | { | |
610 | return "healpix_spheroid"; | |
611 | } | |
612 | ||
613 | }; | |
614 | ||
92f5a8d4 TL |
615 | template <typename T, typename Parameters> |
616 | struct base_rhealpix_ellipsoid | |
11fdf7f2 | 617 | { |
92f5a8d4 | 618 | par_healpix<T> m_proj_parm; |
11fdf7f2 TL |
619 | |
620 | // FORWARD(e_rhealpix_forward) ellipsoid | |
621 | // Project coordinates from geographic (lon, lat) to cartesian (x, y) | |
92f5a8d4 | 622 | inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const |
11fdf7f2 | 623 | { |
92f5a8d4 | 624 | lp_lat = auth_lat(par, m_proj_parm, lp_lat, 0); |
11fdf7f2 TL |
625 | healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); |
626 | combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0); | |
627 | } | |
628 | ||
629 | // INVERSE(e_rhealpix_inverse) ellipsoid | |
630 | // Project coordinates from cartesian (x, y) to geographic (lon, lat) | |
92f5a8d4 | 631 | inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const |
11fdf7f2 TL |
632 | { |
633 | /* Check whether (x, y) lies in the rHEALPix image. */ | |
634 | if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) { | |
635 | lp_lon = HUGE_VAL; | |
636 | lp_lat = HUGE_VAL; | |
92f5a8d4 | 637 | BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); |
11fdf7f2 TL |
638 | } |
639 | combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1); | |
640 | healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); | |
92f5a8d4 | 641 | lp_lat = auth_lat(par, m_proj_parm, lp_lat, 1); |
11fdf7f2 TL |
642 | } |
643 | ||
644 | static inline std::string get_name() | |
645 | { | |
646 | return "rhealpix_ellipsoid"; | |
647 | } | |
648 | ||
649 | }; | |
650 | ||
92f5a8d4 TL |
651 | template <typename T, typename Parameters> |
652 | struct base_rhealpix_spheroid | |
11fdf7f2 | 653 | { |
92f5a8d4 | 654 | par_healpix<T> m_proj_parm; |
11fdf7f2 TL |
655 | |
656 | // FORWARD(s_rhealpix_forward) sphere | |
657 | // Project coordinates from geographic (lon, lat) to cartesian (x, y) | |
92f5a8d4 | 658 | inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const |
11fdf7f2 TL |
659 | { |
660 | healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); | |
661 | combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0); | |
662 | } | |
663 | ||
664 | // INVERSE(s_rhealpix_inverse) sphere | |
665 | // Project coordinates from cartesian (x, y) to geographic (lon, lat) | |
92f5a8d4 | 666 | inline void inv(Parameters const& , T xy_x, T xy_y, T& lp_lon, T& lp_lat) const |
11fdf7f2 TL |
667 | { |
668 | /* Check whether (x, y) lies in the rHEALPix image. */ | |
669 | if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) { | |
670 | lp_lon = HUGE_VAL; | |
671 | lp_lat = HUGE_VAL; | |
92f5a8d4 | 672 | BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); |
11fdf7f2 TL |
673 | } |
674 | combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1); | |
675 | return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); | |
676 | } | |
677 | ||
678 | static inline std::string get_name() | |
679 | { | |
680 | return "rhealpix_spheroid"; | |
681 | } | |
682 | ||
683 | }; | |
684 | ||
685 | // HEALPix | |
686 | template <typename Parameters, typename T> | |
687 | inline void setup_healpix(Parameters& par, par_healpix<T>& proj_parm) | |
688 | { | |
92f5a8d4 TL |
689 | if (par.es != 0.0) { |
690 | proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */ | |
11fdf7f2 TL |
691 | proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */ |
692 | par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */ | |
92f5a8d4 | 693 | pj_calc_ellipsoid_params(par, par.a, par.es); /* Ensure we have a consistent parameter set */ |
11fdf7f2 TL |
694 | } else { |
695 | } | |
696 | } | |
697 | ||
698 | // rHEALPix | |
92f5a8d4 TL |
699 | template <typename Params, typename Parameters, typename T> |
700 | inline void setup_rhealpix(Params const& params, Parameters& par, par_healpix<T>& proj_parm) | |
11fdf7f2 | 701 | { |
92f5a8d4 TL |
702 | proj_parm.north_square = pj_get_param_i<srs::spar::north_square>(params, "north_square", srs::dpar::north_square); |
703 | proj_parm.south_square = pj_get_param_i<srs::spar::south_square>(params, "south_square", srs::dpar::south_square); | |
11fdf7f2 | 704 | /* Check for valid north_square and south_square inputs. */ |
92f5a8d4 TL |
705 | if ((proj_parm.north_square < 0) || (proj_parm.north_square > 3)) { |
706 | BOOST_THROW_EXCEPTION( projection_exception(error_axis) ); | |
11fdf7f2 | 707 | } |
92f5a8d4 TL |
708 | if ((proj_parm.south_square < 0) || (proj_parm.south_square > 3)) { |
709 | BOOST_THROW_EXCEPTION( projection_exception(error_axis) ); | |
11fdf7f2 | 710 | } |
92f5a8d4 TL |
711 | if (par.es != 0.0) { |
712 | proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */ | |
11fdf7f2 TL |
713 | proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */ |
714 | par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */ | |
92f5a8d4 TL |
715 | // TODO: why not the same as in healpix? |
716 | //pj_calc_ellipsoid_params(par, par.a, par.es); | |
11fdf7f2 TL |
717 | par.ra = 1.0/par.a; |
718 | } else { | |
719 | } | |
720 | } | |
721 | ||
722 | }} // namespace detail::healpix | |
723 | #endif // doxygen | |
724 | ||
725 | /*! | |
726 | \brief HEALPix projection | |
727 | \ingroup projections | |
728 | \tparam Geographic latlong point type | |
729 | \tparam Cartesian xy point type | |
730 | \tparam Parameters parameter type | |
731 | \par Projection characteristics | |
732 | - Spheroid | |
733 | - Ellipsoid | |
734 | \par Example | |
735 | \image html ex_healpix.gif | |
736 | */ | |
92f5a8d4 TL |
737 | template <typename T, typename Parameters> |
738 | struct healpix_ellipsoid : public detail::healpix::base_healpix_ellipsoid<T, Parameters> | |
11fdf7f2 | 739 | { |
92f5a8d4 TL |
740 | template <typename Params> |
741 | inline healpix_ellipsoid(Params const& , Parameters & par) | |
11fdf7f2 | 742 | { |
92f5a8d4 | 743 | detail::healpix::setup_healpix(par, this->m_proj_parm); |
11fdf7f2 TL |
744 | } |
745 | }; | |
746 | ||
747 | /*! | |
748 | \brief HEALPix projection | |
749 | \ingroup projections | |
750 | \tparam Geographic latlong point type | |
751 | \tparam Cartesian xy point type | |
752 | \tparam Parameters parameter type | |
753 | \par Projection characteristics | |
754 | - Spheroid | |
755 | - Ellipsoid | |
756 | \par Example | |
757 | \image html ex_healpix.gif | |
758 | */ | |
92f5a8d4 TL |
759 | template <typename T, typename Parameters> |
760 | struct healpix_spheroid : public detail::healpix::base_healpix_spheroid<T, Parameters> | |
11fdf7f2 | 761 | { |
92f5a8d4 TL |
762 | template <typename Params> |
763 | inline healpix_spheroid(Params const& , Parameters & par) | |
11fdf7f2 | 764 | { |
92f5a8d4 | 765 | detail::healpix::setup_healpix(par, this->m_proj_parm); |
11fdf7f2 TL |
766 | } |
767 | }; | |
768 | ||
769 | /*! | |
770 | \brief rHEALPix projection | |
771 | \ingroup projections | |
772 | \tparam Geographic latlong point type | |
773 | \tparam Cartesian xy point type | |
774 | \tparam Parameters parameter type | |
775 | \par Projection characteristics | |
776 | - Spheroid | |
777 | - Ellipsoid | |
778 | \par Projection parameters | |
779 | - north_square (integer) | |
780 | - south_square (integer) | |
781 | \par Example | |
782 | \image html ex_rhealpix.gif | |
783 | */ | |
92f5a8d4 TL |
784 | template <typename T, typename Parameters> |
785 | struct rhealpix_ellipsoid : public detail::healpix::base_rhealpix_ellipsoid<T, Parameters> | |
11fdf7f2 | 786 | { |
92f5a8d4 TL |
787 | template <typename Params> |
788 | inline rhealpix_ellipsoid(Params const& params, Parameters & par) | |
11fdf7f2 | 789 | { |
92f5a8d4 | 790 | detail::healpix::setup_rhealpix(params, par, this->m_proj_parm); |
11fdf7f2 TL |
791 | } |
792 | }; | |
793 | ||
794 | /*! | |
795 | \brief rHEALPix projection | |
796 | \ingroup projections | |
797 | \tparam Geographic latlong point type | |
798 | \tparam Cartesian xy point type | |
799 | \tparam Parameters parameter type | |
800 | \par Projection characteristics | |
801 | - Spheroid | |
802 | - Ellipsoid | |
803 | \par Projection parameters | |
804 | - north_square (integer) | |
805 | - south_square (integer) | |
806 | \par Example | |
807 | \image html ex_rhealpix.gif | |
808 | */ | |
92f5a8d4 TL |
809 | template <typename T, typename Parameters> |
810 | struct rhealpix_spheroid : public detail::healpix::base_rhealpix_spheroid<T, Parameters> | |
11fdf7f2 | 811 | { |
92f5a8d4 TL |
812 | template <typename Params> |
813 | inline rhealpix_spheroid(Params const& params, Parameters & par) | |
11fdf7f2 | 814 | { |
92f5a8d4 | 815 | detail::healpix::setup_rhealpix(params, par, this->m_proj_parm); |
11fdf7f2 TL |
816 | } |
817 | }; | |
818 | ||
819 | #ifndef DOXYGEN_NO_DETAIL | |
820 | namespace detail | |
821 | { | |
822 | ||
823 | // Static projection | |
92f5a8d4 TL |
824 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_healpix, healpix_spheroid, healpix_ellipsoid) |
825 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_rhealpix, rhealpix_spheroid, rhealpix_ellipsoid) | |
11fdf7f2 TL |
826 | |
827 | // Factory entry(s) | |
92f5a8d4 TL |
828 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(healpix_entry, healpix_spheroid, healpix_ellipsoid) |
829 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(rhealpix_entry, rhealpix_spheroid, rhealpix_ellipsoid) | |
830 | ||
831 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(healpix_init) | |
11fdf7f2 | 832 | { |
92f5a8d4 TL |
833 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(healpix, healpix_entry) |
834 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(rhealpix, rhealpix_entry) | |
11fdf7f2 TL |
835 | } |
836 | ||
837 | } // namespace detail | |
838 | #endif // doxygen | |
839 | ||
840 | } // namespace projections | |
841 | ||
842 | }} // namespace boost::geometry | |
843 | ||
844 | #endif // BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP | |
845 |