projections/proj/qsc.hpp

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2017, 2018, 2019.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017
0018 // Last updated version of proj: 5.0.0
0019
0020 // Original copyright notice:
0021
0022 // This implements the Quadrilateralized Spherical Cube (QSC) projection.
0023 // Copyright (c) 2011, 2012  Martin Lambers <marlam@marlam.de>
0024
0025 // Permission is hereby granted, free of charge, to any person obtaining a
0026 // copy of this software and associated documentation files (the "Software"),
0027 // to deal in the Software without restriction, including without limitation
0028 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0029 // and/or sell copies of the Software, and to permit persons to whom the
0030 // Software is furnished to do so, subject to the following conditions:
0031
0032 // The above copyright notice and this permission notice shall be included
0033 // in all copies or substantial portions of the Software.
0034
0035 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0036 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0037 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0038 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0039 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0040 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0041 // DEALINGS IN THE SOFTWARE.
0042
0043 // The QSC projection was introduced in:
0044 // [OL76]
0045 // E.M. O'Neill and R.E. Laubscher, "Extended Studies of a Quadrilateralized
0046 // Spherical Cube Earth Data Base", Naval Environmental Prediction Research
0047 // Facility Tech. Report NEPRF 3-76 (CSC), May 1976.
0048
0049 // The preceding shift from an ellipsoid to a sphere, which allows to apply
0050 // this projection to ellipsoids as used in the Ellipsoidal Cube Map model,
0051 // is described in
0052 // [LK12]
0053 // M. Lambers and A. Kolb, "Ellipsoidal Cube Maps for Accurate Rendering of
0054 // Planetary-Scale Terrain Data", Proc. Pacfic Graphics (Short Papers), Sep.
0055 // 2012
0056
0057 // You have to choose one of the following projection centers,
0058 // corresponding to the centers of the six cube faces:
0059 // phi0 = 0.0, lam0 = 0.0       ("front" face)
0060 // phi0 = 0.0, lam0 = 90.0      ("right" face)
0061 // phi0 = 0.0, lam0 = 180.0     ("back" face)
0062 // phi0 = 0.0, lam0 = -90.0     ("left" face)
0063 // phi0 = 90.0                  ("top" face)
0064 // phi0 = -90.0                 ("bottom" face)
0065 // Other projection centers will not work!
0066
0067 // In the projection code below, each cube face is handled differently.
0068 // See the computation of the face parameter in the ENTRY0(qsc) function
0069 // and the handling of different face values (FACE_*) in the forward and
0070 // inverse projections.
0071
0072 // Furthermore, the projection is originally only defined for theta angles
0073 // between (-1/4 * PI) and (+1/4 * PI) on the current cube face. This area
0074 // of definition is named AREA_0 in the projection code below. The other
0075 // three areas of a cube face are handled by rotation of AREA_0.
0076
0077 #ifndef BOOST_GEOMETRY_PROJECTIONS_QSC_HPP
0078 #define BOOST_GEOMETRY_PROJECTIONS_QSC_HPP
0079
0080 #include <boost/core/ignore_unused.hpp>
0081 #include <boost/geometry/util/math.hpp>
0082
0083 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0084 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0085 #include <boost/geometry/srs/projections/impl/projects.hpp>
0086 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0087
0088 namespace boost { namespace geometry
0089 {
0090
0091 namespace projections
0092 {
0093     #ifndef DOXYGEN_NO_DETAIL
0094     namespace detail { namespace qsc
0095     {
0096
0097             /* The six cube faces. */
0098             enum face_type {
0099                 face_front  = 0,
0100                 face_right  = 1,
0101                 face_back   = 2,
0102                 face_left   = 3,
0103                 face_top    = 4,
0104                 face_bottom = 5
0105             };
0106
0107             template <typename T>
0108             struct par_qsc
0109             {
0110                 T a_squared;
0111                 T b;
0112                 T one_minus_f;
0113                 T one_minus_f_squared;
0114                 face_type face;
0115             };
0116
0117             static const double epsilon10 = 1.e-10;
0118
0119             /* The four areas on a cube face. AREA_0 is the area of definition,
0120              * the other three areas are counted counterclockwise. */
0121             enum area_type {
0122                 area_0 = 0,
0123                 area_1 = 1,
0124                 area_2 = 2,
0125                 area_3 = 3
0126             };
0127
0128             /* Helper function for forward projection: compute the theta angle
0129              * and determine the area number. */
0130             template <typename T>
0131             inline T qsc_fwd_equat_face_theta(T const& phi, T const& y, T const& x, area_type *area)
0132             {
0133                 static const T fourth_pi = detail::fourth_pi<T>();
0134                 static const T half_pi = detail::half_pi<T>();
0135                 static const T pi = detail::pi<T>();
0136
0137                 T theta;
0138                 if (phi < epsilon10) {
0139                     *area = area_0;
0140                     theta = 0.0;
0141                 } else {
0142                     theta = atan2(y, x);
0143                     if (fabs(theta) <= fourth_pi) {
0144                         *area = area_0;
0145                     } else if (theta > fourth_pi && theta <= half_pi + fourth_pi) {
0146                         *area = area_1;
0147                         theta -= half_pi;
0148                     } else if (theta > half_pi + fourth_pi || theta <= -(half_pi + fourth_pi)) {
0149                         *area = area_2;
0150                         theta = (theta >= 0.0 ? theta - pi : theta + pi);
0151                     } else {
0152                         *area = area_3;
0153                         theta += half_pi;
0154                     }
0155                 }
0156                 return theta;
0157             }
0158
0159             /* Helper function: shift the longitude. */
0160             template <typename T>
0161             inline T qsc_shift_lon_origin(T const& lon, T const& offset)
0162             {
0163                 static const T pi = detail::pi<T>();
0164                 static const T two_pi = detail::two_pi<T>();
0165
0166                 T slon = lon + offset;
0167                 if (slon < -pi) {
0168                     slon += two_pi;
0169                 } else if (slon > +pi) {
0170                     slon -= two_pi;
0171                 }
0172                 return slon;
0173             }
0174
0175             /* Forward projection, ellipsoid */
0176
0177             template <typename T, typename Parameters>
0178             struct base_qsc_ellipsoid
0179             {
0180                 par_qsc<T> m_proj_parm;
0181
0182                 // FORWARD(e_forward)
0183                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0184                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0185                 {
0186                     static const T fourth_pi = detail::fourth_pi<T>();
0187                     static const T half_pi = detail::half_pi<T>();
0188                     static const T pi = detail::pi<T>();
0189
0190                         T lat, lon;
0191                         T theta, phi;
0192                         T t, mu; /* nu; */
0193                         area_type area;
0194
0195                         /* Convert the geodetic latitude to a geocentric latitude.
0196                          * This corresponds to the shift from the ellipsoid to the sphere
0197                          * described in [LK12]. */
0198                         if (par.es != 0.0) {
0199                             lat = atan(this->m_proj_parm.one_minus_f_squared * tan(lp_lat));
0200                         } else {
0201                             lat = lp_lat;
0202                         }
0203
0204                         /* Convert the input lat, lon into theta, phi as used by QSC.
0205                          * This depends on the cube face and the area on it.
0206                          * For the top and bottom face, we can compute theta and phi
0207                          * directly from phi, lam. For the other faces, we must use
0208                          * unit sphere cartesian coordinates as an intermediate step. */
0209                         lon = lp_lon;
0210                         if (this->m_proj_parm.face == face_top) {
0211                             phi = half_pi - lat;
0212                             if (lon >= fourth_pi && lon <= half_pi + fourth_pi) {
0213                                 area = area_0;
0214                                 theta = lon - half_pi;
0215                             } else if (lon > half_pi + fourth_pi || lon <= -(half_pi + fourth_pi)) {
0216                                 area = area_1;
0217                                 theta = (lon > 0.0 ? lon - pi : lon + pi);
0218                             } else if (lon > -(half_pi + fourth_pi) && lon <= -fourth_pi) {
0219                                 area = area_2;
0220                                 theta = lon + half_pi;
0221                             } else {
0222                                 area = area_3;
0223                                 theta = lon;
0224                             }
0225                         } else if (this->m_proj_parm.face == face_bottom) {
0226                             phi = half_pi + lat;
0227                             if (lon >= fourth_pi && lon <= half_pi + fourth_pi) {
0228                                 area = area_0;
0229                                 theta = -lon + half_pi;
0230                             } else if (lon < fourth_pi && lon >= -fourth_pi) {
0231                                 area = area_1;
0232                                 theta = -lon;
0233                             } else if (lon < -fourth_pi && lon >= -(half_pi + fourth_pi)) {
0234                                 area = area_2;
0235                                 theta = -lon - half_pi;
0236                             } else {
0237                                 area = area_3;
0238                                 theta = (lon > 0.0 ? -lon + pi : -lon - pi);
0239                             }
0240                         } else {
0241                             T q, r, s;
0242                             T sinlat, coslat;
0243                             T sinlon, coslon;
0244
0245                             if (this->m_proj_parm.face == face_right) {
0246                                 lon = qsc_shift_lon_origin(lon, +half_pi);
0247                             } else if (this->m_proj_parm.face == face_back) {
0248                                 lon = qsc_shift_lon_origin(lon, +pi);
0249                             } else if (this->m_proj_parm.face == face_left) {
0250                                 lon = qsc_shift_lon_origin(lon, -half_pi);
0251                             }
0252                             sinlat = sin(lat);
0253                             coslat = cos(lat);
0254                             sinlon = sin(lon);
0255                             coslon = cos(lon);
0256                             q = coslat * coslon;
0257                             r = coslat * sinlon;
0258                             s = sinlat;
0259
0260                             if (this->m_proj_parm.face == face_front) {
0261                                 phi = acos(q);
0262                                 theta = qsc_fwd_equat_face_theta(phi, s, r, &area);
0263                             } else if (this->m_proj_parm.face == face_right) {
0264                                 phi = acos(r);
0265                                 theta = qsc_fwd_equat_face_theta(phi, s, -q, &area);
0266                             } else if (this->m_proj_parm.face == face_back) {
0267                                 phi = acos(-q);
0268                                 theta = qsc_fwd_equat_face_theta(phi, s, -r, &area);
0269                             } else if (this->m_proj_parm.face == face_left) {
0270                                 phi = acos(-r);
0271                                 theta = qsc_fwd_equat_face_theta(phi, s, q, &area);
0272                             } else {
0273                                 /* Impossible */
0274                                 phi = theta = 0.0;
0275                                 area = area_0;
0276                             }
0277                         }
0278
0279                         /* Compute mu and nu for the area of definition.
0280                          * For mu, see Eq. (3-21) in [OL76], but note the typos:
0281                          * compare with Eq. (3-14). For nu, see Eq. (3-38). */
0282                         mu = atan((12.0 / pi) * (theta + acos(sin(theta) * cos(fourth_pi)) - half_pi));
0283                         // TODO: (cos(mu) * cos(mu)) could be replaced with sqr(cos(mu))
0284                         t = sqrt((1.0 - cos(phi)) / (cos(mu) * cos(mu)) / (1.0 - cos(atan(1.0 / cos(theta)))));
0285                         /* nu = atan(t);        We don't really need nu, just t, see below. */
0286
0287                         /* Apply the result to the real area. */
0288                         if (area == area_1) {
0289                             mu += half_pi;
0290                         } else if (area == area_2) {
0291                             mu += pi;
0292                         } else if (area == area_3) {
0293                             mu += half_pi + pi;
0294                         }
0295
0296                         /* Now compute x, y from mu and nu */
0297                         /* t = tan(nu); */
0298                         xy_x = t * cos(mu);
0299                         xy_y = t * sin(mu);
0300                 }
0301                 /* Inverse projection, ellipsoid */
0302
0303                 // INVERSE(e_inverse)
0304                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0305                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0306                 {
0307                     static const T half_pi = detail::half_pi<T>();
0308                     static const T pi = detail::pi<T>();
0309
0310                         T mu, nu, cosmu, tannu;
0311                         T tantheta, theta, cosphi, phi;
0312                         T t;
0313                         int area;
0314
0315                         /* Convert the input x, y to the mu and nu angles as used by QSC.
0316                          * This depends on the area of the cube face. */
0317                         nu = atan(sqrt(xy_x * xy_x + xy_y * xy_y));
0318                         mu = atan2(xy_y, xy_x);
0319                         if (xy_x >= 0.0 && xy_x >= fabs(xy_y)) {
0320                             area = area_0;
0321                         } else if (xy_y >= 0.0 && xy_y >= fabs(xy_x)) {
0322                             area = area_1;
0323                             mu -= half_pi;
0324                         } else if (xy_x < 0.0 && -xy_x >= fabs(xy_y)) {
0325                             area = area_2;
0326                             mu = (mu < 0.0 ? mu + pi : mu - pi);
0327                         } else {
0328                             area = area_3;
0329                             mu += half_pi;
0330                         }
0331
0332                         /* Compute phi and theta for the area of definition.
0333                          * The inverse projection is not described in the original paper, but some
0334                          * good hints can be found here (as of 2011-12-14):
0335                          * http://fits.gsfc.nasa.gov/fitsbits/saf.93/saf.9302
0336                          * (search for "Message-Id: <9302181759.AA25477 at fits.cv.nrao.edu>") */
0337                         t = (pi / 12.0) * tan(mu);
0338                         tantheta = sin(t) / (cos(t) - (1.0 / sqrt(2.0)));
0339                         theta = atan(tantheta);
0340                         cosmu = cos(mu);
0341                         tannu = tan(nu);
0342                         cosphi = 1.0 - cosmu * cosmu * tannu * tannu * (1.0 - cos(atan(1.0 / cos(theta))));
0343                         if (cosphi < -1.0) {
0344                             cosphi = -1.0;
0345                         } else if (cosphi > +1.0) {
0346                             cosphi = +1.0;
0347                         }
0348
0349                         /* Apply the result to the real area on the cube face.
0350                          * For the top and bottom face, we can compute phi and lam directly.
0351                          * For the other faces, we must use unit sphere cartesian coordinates
0352                          * as an intermediate step. */
0353                         if (this->m_proj_parm.face == face_top) {
0354                             phi = acos(cosphi);
0355                             lp_lat = half_pi - phi;
0356                             if (area == area_0) {
0357                                 lp_lon = theta + half_pi;
0358                             } else if (area == area_1) {
0359                                 lp_lon = (theta < 0.0 ? theta + pi : theta - pi);
0360                             } else if (area == area_2) {
0361                                 lp_lon = theta - half_pi;
0362                             } else /* area == AREA_3 */ {
0363                                 lp_lon = theta;
0364                             }
0365                         } else if (this->m_proj_parm.face == face_bottom) {
0366                             phi = acos(cosphi);
0367                             lp_lat = phi - half_pi;
0368                             if (area == area_0) {
0369                                 lp_lon = -theta + half_pi;
0370                             } else if (area == area_1) {
0371                                 lp_lon = -theta;
0372                             } else if (area == area_2) {
0373                                 lp_lon = -theta - half_pi;
0374                             } else /* area == area_3 */ {
0375                                 lp_lon = (theta < 0.0 ? -theta - pi : -theta + pi);
0376                             }
0377                         } else {
0378                             /* Compute phi and lam via cartesian unit sphere coordinates. */
0379                             T q, r, s;
0380                             q = cosphi;
0381                             t = q * q;
0382                             if (t >= 1.0) {
0383                                 s = 0.0;
0384                             } else {
0385                                 s = sqrt(1.0 - t) * sin(theta);
0386                             }
0387                             t += s * s;
0388                             if (t >= 1.0) {
0389                                 r = 0.0;
0390                             } else {
0391                                 r = sqrt(1.0 - t);
0392                             }
0393                             /* Rotate q,r,s into the correct area. */
0394                             if (area == area_1) {
0395                                 t = r;
0396                                 r = -s;
0397                                 s = t;
0398                             } else if (area == area_2) {
0399                                 r = -r;
0400                                 s = -s;
0401                             } else if (area == area_3) {
0402                                 t = r;
0403                                 r = s;
0404                                 s = -t;
0405                             }
0406                             /* Rotate q,r,s into the correct cube face. */
0407                             if (this->m_proj_parm.face == face_right) {
0408                                 t = q;
0409                                 q = -r;
0410                                 r = t;
0411                             } else if (this->m_proj_parm.face == face_back) {
0412                                 q = -q;
0413                                 r = -r;
0414                             } else if (this->m_proj_parm.face == face_left) {
0415                                 t = q;
0416                                 q = r;
0417                                 r = -t;
0418                             }
0419                             /* Now compute phi and lam from the unit sphere coordinates. */
0420                             lp_lat = acos(-s) - half_pi;
0421                             lp_lon = atan2(r, q);
0422                             if (this->m_proj_parm.face == face_right) {
0423                                 lp_lon = qsc_shift_lon_origin(lp_lon, -half_pi);
0424                             } else if (this->m_proj_parm.face == face_back) {
0425                                 lp_lon = qsc_shift_lon_origin(lp_lon, -pi);
0426                             } else if (this->m_proj_parm.face == face_left) {
0427                                 lp_lon = qsc_shift_lon_origin(lp_lon, +half_pi);
0428                             }
0429                         }
0430
0431                         /* Apply the shift from the sphere to the ellipsoid as described
0432                          * in [LK12]. */
0433                         if (par.es != 0.0) {
0434                             int invert_sign;
0435                             T tanphi, xa;
0436                             invert_sign = (lp_lat < 0.0 ? 1 : 0);
0437                             tanphi = tan(lp_lat);
0438                             xa = this->m_proj_parm.b / sqrt(tanphi * tanphi + this->m_proj_parm.one_minus_f_squared);
0439                             lp_lat = atan(sqrt(par.a * par.a - xa * xa) / (this->m_proj_parm.one_minus_f * xa));
0440                             if (invert_sign) {
0441                                 lp_lat = -lp_lat;
0442                             }
0443                         }
0444                 }
0445
0446                 static inline std::string get_name()
0447                 {
0448                     return "qsc_ellipsoid";
0449                 }
0450
0451             };
0452
0453             // Quadrilateralized Spherical Cube
0454             template <typename Parameters, typename T>
0455             inline void setup_qsc(Parameters const& par, par_qsc<T>& proj_parm)
0456             {
0457                 static const T fourth_pi = detail::fourth_pi<T>();
0458                 static const T half_pi = detail::half_pi<T>();
0459
0460                 /* Determine the cube face from the center of projection. */
0461                 if (par.phi0 >= half_pi - fourth_pi / 2.0) {
0462                     proj_parm.face = face_top;
0463                 } else if (par.phi0 <= -(half_pi - fourth_pi / 2.0)) {
0464                     proj_parm.face = face_bottom;
0465                 } else if (fabs(par.lam0) <= fourth_pi) {
0466                     proj_parm.face = face_front;
0467                 } else if (fabs(par.lam0) <= half_pi + fourth_pi) {
0468                     proj_parm.face = (par.lam0 > 0.0 ? face_right : face_left);
0469                 } else {
0470                     proj_parm.face = face_back;
0471                 }
0472                 /* Fill in useful values for the ellipsoid <-> sphere shift
0473                  * described in [LK12]. */
0474                 if (par.es != 0.0) {
0475                     proj_parm.a_squared = par.a * par.a;
0476                     proj_parm.b = par.a * sqrt(1.0 - par.es);
0477                     proj_parm.one_minus_f = 1.0 - (par.a - proj_parm.b) / par.a;
0478                     proj_parm.one_minus_f_squared = proj_parm.one_minus_f * proj_parm.one_minus_f;
0479                 }
0480             }
0481
0482     }} // namespace detail::qsc
0483     #endif // doxygen
0484
0485     /*!
0486         \brief Quadrilateralized Spherical Cube projection
0487         \ingroup projections
0488         \tparam Geographic latlong point type
0489         \tparam Cartesian xy point type
0490         \tparam Parameters parameter type
0491         \par Projection characteristics
0492          - Azimuthal
0493          - Spheroid
0494         \par Example
0495         \image html ex_qsc.gif
0496     */
0497     template <typename T, typename Parameters>
0498     struct qsc_ellipsoid : public detail::qsc::base_qsc_ellipsoid<T, Parameters>
0499     {
0500         template <typename Params>
0501         inline qsc_ellipsoid(Params const& , Parameters const& par)
0502         {
0503             detail::qsc::setup_qsc(par, this->m_proj_parm);
0504         }
0505     };
0506
0507     #ifndef DOXYGEN_NO_DETAIL
0508     namespace detail
0509     {
0510
0511         // Static projection
0512         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_qsc, qsc_ellipsoid)
0513
0514         // Factory entry(s)
0515         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(qsc_entry, qsc_ellipsoid)
0516
0517         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(qsc_init)
0518         {
0519             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(qsc, qsc_entry)
0520         }
0521
0522     } // namespace detail
0523     #endif // doxygen
0524
0525 } // namespace projections
0526
0527 }} // namespace boost::geometry
0528
0529 #endif // BOOST_GEOMETRY_PROJECTIONS_QSC_HPP
0530