[ceph.git] / ceph / src / boost / boost / geometry / srs / projections / proj / qsc.hpp

// Boost.Geometry - gis-projections (based on PROJ4)

// Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.

// This file was modified by Oracle on 2017, 2018, 2019.
// Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

// This file is converted from PROJ4, http://trac.osgeo.org/proj
// PROJ4 is originally written by Gerald Evenden (then of the USGS)
// PROJ4 is maintained by Frank Warmerdam
// PROJ4 is converted to Boost.Geometry by Barend Gehrels

// Last updated version of proj: 5.0.0

// Original copyright notice:

// This implements the Quadrilateralized Spherical Cube (QSC) projection.
// Copyright (c) 2011, 2012  Martin Lambers <marlam@marlam.de>

// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.

// The QSC projection was introduced in:
// [OL76]
// E.M. O'Neill and R.E. Laubscher, "Extended Studies of a Quadrilateralized
// Spherical Cube Earth Data Base", Naval Environmental Prediction Research
// Facility Tech. Report NEPRF 3-76 (CSC), May 1976.

// The preceding shift from an ellipsoid to a sphere, which allows to apply
// this projection to ellipsoids as used in the Ellipsoidal Cube Map model,
// is described in
// [LK12]
// M. Lambers and A. Kolb, "Ellipsoidal Cube Maps for Accurate Rendering of
// Planetary-Scale Terrain Data", Proc. Pacfic Graphics (Short Papers), Sep.
// 2012

// You have to choose one of the following projection centers,
// corresponding to the centers of the six cube faces:
// phi0 = 0.0, lam0 = 0.0       ("front" face)
// phi0 = 0.0, lam0 = 90.0      ("right" face)
// phi0 = 0.0, lam0 = 180.0     ("back" face)
// phi0 = 0.0, lam0 = -90.0     ("left" face)
// phi0 = 90.0                  ("top" face)
// phi0 = -90.0                 ("bottom" face)
// Other projection centers will not work!

// In the projection code below, each cube face is handled differently.
// See the computation of the face parameter in the ENTRY0(qsc) function
// and the handling of different face values (FACE_*) in the forward and
// inverse projections.

// Furthermore, the projection is originally only defined for theta angles
// between (-1/4 * PI) and (+1/4 * PI) on the current cube face. This area
// of definition is named AREA_0 in the projection code below. The other
// three areas of a cube face are handled by rotation of AREA_0.

#ifndef BOOST_GEOMETRY_PROJECTIONS_QSC_HPP
#define BOOST_GEOMETRY_PROJECTIONS_QSC_HPP

#include <boost/core/ignore_unused.hpp>
#include <boost/geometry/util/math.hpp>

#include <boost/geometry/srs/projections/impl/base_static.hpp>
#include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
#include <boost/geometry/srs/projections/impl/projects.hpp>
#include <boost/geometry/srs/projections/impl/factory_entry.hpp>

namespace boost { namespace geometry
{

namespace projections
{
    #ifndef DOXYGEN_NO_DETAIL
    namespace detail { namespace qsc
    {

            /* The six cube faces. */
            enum face_type {
                face_front  = 0,
                face_right  = 1,
                face_back   = 2,
                face_left   = 3,
                face_top    = 4,
                face_bottom = 5
            };

            template <typename T>
            struct par_qsc
            {
                T a_squared;
                T b;
                T one_minus_f;
                T one_minus_f_squared;
                face_type face;
            };

            static const double epsilon10 = 1.e-10;

            /* The four areas on a cube face. AREA_0 is the area of definition,
             * the other three areas are counted counterclockwise. */
            enum area_type {
                area_0 = 0,
                area_1 = 1,
                area_2 = 2,
                area_3 = 3
            };

            /* Helper function for forward projection: compute the theta angle
             * and determine the area number. */
            template <typename T>
            inline T qsc_fwd_equat_face_theta(T const& phi, T const& y, T const& x, area_type *area)
            {
                static const T fourth_pi = detail::fourth_pi<T>();
                static const T half_pi = detail::half_pi<T>();
                static const T pi = detail::pi<T>();

                T theta;
                if (phi < epsilon10) {
                    *area = area_0;
                    theta = 0.0;
                } else {
                    theta = atan2(y, x);
                    if (fabs(theta) <= fourth_pi) {
                        *area = area_0;
                    } else if (theta > fourth_pi && theta <= half_pi + fourth_pi) {
                        *area = area_1;
                        theta -= half_pi;
                    } else if (theta > half_pi + fourth_pi || theta <= -(half_pi + fourth_pi)) {
                        *area = area_2;
                        theta = (theta >= 0.0 ? theta - pi : theta + pi);
                    } else {
                        *area = area_3;
                        theta += half_pi;
                    }
                }
                return theta;
            }

            /* Helper function: shift the longitude. */
            template <typename T>
            inline T qsc_shift_lon_origin(T const& lon, T const& offset)
            {
                static const T pi = detail::pi<T>();
                static const T two_pi = detail::two_pi<T>();

                T slon = lon + offset;
                if (slon < -pi) {
                    slon += two_pi;
                } else if (slon > +pi) {
                    slon -= two_pi;
                }
                return slon;
            }

            /* Forward projection, ellipsoid */

            template <typename T, typename Parameters>
            struct base_qsc_ellipsoid
            {
                par_qsc<T> m_proj_parm;

                // FORWARD(e_forward)
                // Project coordinates from geographic (lon, lat) to cartesian (x, y)
                inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
                {
                    static const T fourth_pi = detail::fourth_pi<T>();
                    static const T half_pi = detail::half_pi<T>();
                    static const T pi = detail::pi<T>();

                        T lat, lon;
                        T theta, phi;
                        T t, mu; /* nu; */ 
                        area_type area;

                        /* Convert the geodetic latitude to a geocentric latitude.
                         * This corresponds to the shift from the ellipsoid to the sphere
                         * described in [LK12]. */
                        if (par.es != 0.0) {
                            lat = atan(this->m_proj_parm.one_minus_f_squared * tan(lp_lat));
                        } else {
                            lat = lp_lat;
                        }

                        /* Convert the input lat, lon into theta, phi as used by QSC.
                         * This depends on the cube face and the area on it.
                         * For the top and bottom face, we can compute theta and phi
                         * directly from phi, lam. For the other faces, we must use
                         * unit sphere cartesian coordinates as an intermediate step. */
                        lon = lp_lon;
                        if (this->m_proj_parm.face == face_top) {
                            phi = half_pi - lat;
                            if (lon >= fourth_pi && lon <= half_pi + fourth_pi) {
                                area = area_0;
                                theta = lon - half_pi;
                            } else if (lon > half_pi + fourth_pi || lon <= -(half_pi + fourth_pi)) {
                                area = area_1;
                                theta = (lon > 0.0 ? lon - pi : lon + pi);
                            } else if (lon > -(half_pi + fourth_pi) && lon <= -fourth_pi) {
                                area = area_2;
                                theta = lon + half_pi;
                            } else {
                                area = area_3;
                                theta = lon;
                            }
                        } else if (this->m_proj_parm.face == face_bottom) {
                            phi = half_pi + lat;
                            if (lon >= fourth_pi && lon <= half_pi + fourth_pi) {
                                area = area_0;
                                theta = -lon + half_pi;
                            } else if (lon < fourth_pi && lon >= -fourth_pi) {
                                area = area_1;
                                theta = -lon;
                            } else if (lon < -fourth_pi && lon >= -(half_pi + fourth_pi)) {
                                area = area_2;
                                theta = -lon - half_pi;
                            } else {
                                area = area_3;
                                theta = (lon > 0.0 ? -lon + pi : -lon - pi);
                            }
                        } else {
                            T q, r, s;
                            T sinlat, coslat;
                            T sinlon, coslon;

                            if (this->m_proj_parm.face == face_right) {
                                lon = qsc_shift_lon_origin(lon, +half_pi);
                            } else if (this->m_proj_parm.face == face_back) {
                                lon = qsc_shift_lon_origin(lon, +pi);
                            } else if (this->m_proj_parm.face == face_left) {
                                lon = qsc_shift_lon_origin(lon, -half_pi);
                            }
                            sinlat = sin(lat);
                            coslat = cos(lat);
                            sinlon = sin(lon);
                            coslon = cos(lon);
                            q = coslat * coslon;
                            r = coslat * sinlon;
                            s = sinlat;

                            if (this->m_proj_parm.face == face_front) {
                                phi = acos(q);
                                theta = qsc_fwd_equat_face_theta(phi, s, r, &area);
                            } else if (this->m_proj_parm.face == face_right) {
                                phi = acos(r);
                                theta = qsc_fwd_equat_face_theta(phi, s, -q, &area);
                            } else if (this->m_proj_parm.face == face_back) {
                                phi = acos(-q);
                                theta = qsc_fwd_equat_face_theta(phi, s, -r, &area);
                            } else if (this->m_proj_parm.face == face_left) {
                                phi = acos(-r);
                                theta = qsc_fwd_equat_face_theta(phi, s, q, &area);
                            } else {
                                /* Impossible */
                                phi = theta = 0.0;
                                area = area_0;
                            }
                        }

                        /* Compute mu and nu for the area of definition.
                         * For mu, see Eq. (3-21) in [OL76], but note the typos:
                         * compare with Eq. (3-14). For nu, see Eq. (3-38). */
                        mu = atan((12.0 / pi) * (theta + acos(sin(theta) * cos(fourth_pi)) - half_pi));
                        // TODO: (cos(mu) * cos(mu)) could be replaced with sqr(cos(mu))
                        t = sqrt((1.0 - cos(phi)) / (cos(mu) * cos(mu)) / (1.0 - cos(atan(1.0 / cos(theta)))));
                        /* nu = atan(t);        We don't really need nu, just t, see below. */

                        /* Apply the result to the real area. */
                        if (area == area_1) {
                            mu += half_pi;
                        } else if (area == area_2) {
                            mu += pi;
                        } else if (area == area_3) {
                            mu += half_pi + pi;
                        }

                        /* Now compute x, y from mu and nu */
                        /* t = tan(nu); */
                        xy_x = t * cos(mu);
                        xy_y = t * sin(mu);
                }
                /* Inverse projection, ellipsoid */

                // INVERSE(e_inverse)
                // Project coordinates from cartesian (x, y) to geographic (lon, lat)
                inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
                {
                    static const T half_pi = detail::half_pi<T>();
                    static const T pi = detail::pi<T>();

                        T mu, nu, cosmu, tannu;
                        T tantheta, theta, cosphi, phi;
                        T t;
                        int area;

                        /* Convert the input x, y to the mu and nu angles as used by QSC.
                         * This depends on the area of the cube face. */
                        nu = atan(sqrt(xy_x * xy_x + xy_y * xy_y));
                        mu = atan2(xy_y, xy_x);
                        if (xy_x >= 0.0 && xy_x >= fabs(xy_y)) {
                            area = area_0;
                        } else if (xy_y >= 0.0 && xy_y >= fabs(xy_x)) {
                            area = area_1;
                            mu -= half_pi;
                        } else if (xy_x < 0.0 && -xy_x >= fabs(xy_y)) {
                            area = area_2;
                            mu = (mu < 0.0 ? mu + pi : mu - pi);
                        } else {
                            area = area_3;
                            mu += half_pi;
                        }

                        /* Compute phi and theta for the area of definition.
                         * The inverse projection is not described in the original paper, but some
                         * good hints can be found here (as of 2011-12-14):
                         * http://fits.gsfc.nasa.gov/fitsbits/saf.93/saf.9302
                         * (search for "Message-Id: <9302181759.AA25477 at fits.cv.nrao.edu>") */
                        t = (pi / 12.0) * tan(mu);
                        tantheta = sin(t) / (cos(t) - (1.0 / sqrt(2.0)));
                        theta = atan(tantheta);
                        cosmu = cos(mu);
                        tannu = tan(nu);
                        cosphi = 1.0 - cosmu * cosmu * tannu * tannu * (1.0 - cos(atan(1.0 / cos(theta))));
                        if (cosphi < -1.0) {
                            cosphi = -1.0;
                        } else if (cosphi > +1.0) {
                            cosphi = +1.0;
                        }

                        /* Apply the result to the real area on the cube face.
                         * For the top and bottom face, we can compute phi and lam directly.
                         * For the other faces, we must use unit sphere cartesian coordinates
                         * as an intermediate step. */
                        if (this->m_proj_parm.face == face_top) {
                            phi = acos(cosphi);
                            lp_lat = half_pi - phi;
                            if (area == area_0) {
                                lp_lon = theta + half_pi;
                            } else if (area == area_1) {
                                lp_lon = (theta < 0.0 ? theta + pi : theta - pi);
                            } else if (area == area_2) {
                                lp_lon = theta - half_pi;
                            } else /* area == AREA_3 */ {
                                lp_lon = theta;
                            }
                        } else if (this->m_proj_parm.face == face_bottom) {
                            phi = acos(cosphi);
                            lp_lat = phi - half_pi;
                            if (area == area_0) {
                                lp_lon = -theta + half_pi;
                            } else if (area == area_1) {
                                lp_lon = -theta;
                            } else if (area == area_2) {
                                lp_lon = -theta - half_pi;
                            } else /* area == area_3 */ {
                                lp_lon = (theta < 0.0 ? -theta - pi : -theta + pi);
                            }
                        } else {
                            /* Compute phi and lam via cartesian unit sphere coordinates. */
                            T q, r, s;
                            q = cosphi;
                            t = q * q;
                            if (t >= 1.0) {
                                s = 0.0;
                            } else {
                                s = sqrt(1.0 - t) * sin(theta);
                            }
                            t += s * s;
                            if (t >= 1.0) {
                                r = 0.0;
                            } else {
                                r = sqrt(1.0 - t);
                            }
                            /* Rotate q,r,s into the correct area. */
                            if (area == area_1) {
                                t = r;
                                r = -s;
                                s = t;
                            } else if (area == area_2) {
                                r = -r;
                                s = -s;
                            } else if (area == area_3) {
                                t = r;
                                r = s;
                                s = -t;
                            }
                            /* Rotate q,r,s into the correct cube face. */
                            if (this->m_proj_parm.face == face_right) {
                                t = q;
                                q = -r;
                                r = t;
                            } else if (this->m_proj_parm.face == face_back) {
                                q = -q;
                                r = -r;
                            } else if (this->m_proj_parm.face == face_left) {
                                t = q;
                                q = r;
                                r = -t;
                            }
                            /* Now compute phi and lam from the unit sphere coordinates. */
                            lp_lat = acos(-s) - half_pi;
                            lp_lon = atan2(r, q);
                            if (this->m_proj_parm.face == face_right) {
                                lp_lon = qsc_shift_lon_origin(lp_lon, -half_pi);
                            } else if (this->m_proj_parm.face == face_back) {
                                lp_lon = qsc_shift_lon_origin(lp_lon, -pi);
                            } else if (this->m_proj_parm.face == face_left) {
                                lp_lon = qsc_shift_lon_origin(lp_lon, +half_pi);
                            }
                        }

                        /* Apply the shift from the sphere to the ellipsoid as described
                         * in [LK12]. */
                        if (par.es != 0.0) {
                            int invert_sign;
                            T tanphi, xa;
                            invert_sign = (lp_lat < 0.0 ? 1 : 0);
                            tanphi = tan(lp_lat);
                            xa = this->m_proj_parm.b / sqrt(tanphi * tanphi + this->m_proj_parm.one_minus_f_squared);
                            lp_lat = atan(sqrt(par.a * par.a - xa * xa) / (this->m_proj_parm.one_minus_f * xa));
                            if (invert_sign) {
                                lp_lat = -lp_lat;
                            }
                        }
                }

                static inline std::string get_name()
                {
                    return "qsc_ellipsoid";
                }

            };

            // Quadrilateralized Spherical Cube
            template <typename Parameters, typename T>
            inline void setup_qsc(Parameters const& par, par_qsc<T>& proj_parm)
            {
                static const T fourth_pi = detail::fourth_pi<T>();
                static const T half_pi = detail::half_pi<T>();

                /* Determine the cube face from the center of projection. */
                if (par.phi0 >= half_pi - fourth_pi / 2.0) {
                    proj_parm.face = face_top;
                } else if (par.phi0 <= -(half_pi - fourth_pi / 2.0)) {
                    proj_parm.face = face_bottom;
                } else if (fabs(par.lam0) <= fourth_pi) {
                    proj_parm.face = face_front;
                } else if (fabs(par.lam0) <= half_pi + fourth_pi) {
                    proj_parm.face = (par.lam0 > 0.0 ? face_right : face_left);
                } else {
                    proj_parm.face = face_back;
                }
                /* Fill in useful values for the ellipsoid <-> sphere shift
                 * described in [LK12]. */
                if (par.es != 0.0) {
                    proj_parm.a_squared = par.a * par.a;
                    proj_parm.b = par.a * sqrt(1.0 - par.es);
                    proj_parm.one_minus_f = 1.0 - (par.a - proj_parm.b) / par.a;
                    proj_parm.one_minus_f_squared = proj_parm.one_minus_f * proj_parm.one_minus_f;
                }
            }

    }} // namespace detail::qsc
    #endif // doxygen

    /*!
        \brief Quadrilateralized Spherical Cube projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Azimuthal
         - Spheroid
        \par Example
        \image html ex_qsc.gif
    */
    template <typename T, typename Parameters>
    struct qsc_ellipsoid : public detail::qsc::base_qsc_ellipsoid<T, Parameters>
    {
        template <typename Params>
        inline qsc_ellipsoid(Params const& , Parameters const& par)
        {
            detail::qsc::setup_qsc(par, this->m_proj_parm);
        }
    };

    #ifndef DOXYGEN_NO_DETAIL
    namespace detail
    {

        // Static projection
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_qsc, qsc_ellipsoid)

        // Factory entry(s)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(qsc_entry, qsc_ellipsoid)

        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(qsc_init)
        {
            BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(qsc, qsc_entry)
        }

    } // namespace detail
    #endif // doxygen

} // namespace projections

}} // namespace boost::geometry

#endif // BOOST_GEOMETRY_PROJECTIONS_QSC_HPP
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.
92f5a8d4 TL	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	// This implements the Quadrilateralized Spherical Cube (QSC) projection.
	23	// Copyright (c) 2011, 2012 Martin Lambers <marlam@marlam.de>
92f5a8d4 TL	24
	25	// Permission is hereby granted, free of charge, to any person obtaining a
	26	// copy of this software and associated documentation files (the "Software"),
	27	// to deal in the Software without restriction, including without limitation
	28	// the rights to use, copy, modify, merge, publish, distribute, sublicense,
	29	// and/or sell copies of the Software, and to permit persons to whom the
	30	// Software is furnished to do so, subject to the following conditions:
	31
	32	// The above copyright notice and this permission notice shall be included
	33	// in all copies or substantial portions of the Software.
	34
	35	// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	36	// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	37	// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
	38	// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	39	// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
	40	// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
	41	// DEALINGS IN THE SOFTWARE.
	42
11fdf7f2 TL	43	// The QSC projection was introduced in:
	44	// [OL76]
	45	// E.M. O'Neill and R.E. Laubscher, "Extended Studies of a Quadrilateralized
	46	// Spherical Cube Earth Data Base", Naval Environmental Prediction Research
	47	// Facility Tech. Report NEPRF 3-76 (CSC), May 1976.
92f5a8d4	48
11fdf7f2 TL	49	// The preceding shift from an ellipsoid to a sphere, which allows to apply
	50	// this projection to ellipsoids as used in the Ellipsoidal Cube Map model,
	51	// is described in
	52	// [LK12]
	53	// M. Lambers and A. Kolb, "Ellipsoidal Cube Maps for Accurate Rendering of
	54	// Planetary-Scale Terrain Data", Proc. Pacfic Graphics (Short Papers), Sep.
	55	// 2012
92f5a8d4	56
11fdf7f2 TL	57	// You have to choose one of the following projection centers,
	58	// corresponding to the centers of the six cube faces:
	59	// phi0 = 0.0, lam0 = 0.0 ("front" face)
	60	// phi0 = 0.0, lam0 = 90.0 ("right" face)
	61	// phi0 = 0.0, lam0 = 180.0 ("back" face)
	62	// phi0 = 0.0, lam0 = -90.0 ("left" face)
	63	// phi0 = 90.0 ("top" face)
	64	// phi0 = -90.0 ("bottom" face)
	65	// Other projection centers will not work!
92f5a8d4	66
11fdf7f2 TL	67	// In the projection code below, each cube face is handled differently.
	68	// See the computation of the face parameter in the ENTRY0(qsc) function
	69	// and the handling of different face values (FACE_*) in the forward and
	70	// inverse projections.
92f5a8d4	71
11fdf7f2 TL	72	// Furthermore, the projection is originally only defined for theta angles
	73	// between (-1/4 * PI) and (+1/4 * PI) on the current cube face. This area
	74	// of definition is named AREA_0 in the projection code below. The other
	75	// three areas of a cube face are handled by rotation of AREA_0.
	76
92f5a8d4 TL	77	#ifndef BOOST_GEOMETRY_PROJECTIONS_QSC_HPP
92f5a8d4 TL	78	#define BOOST_GEOMETRY_PROJECTIONS_QSC_HPP
11fdf7f2 TL	79
	80	#include <boost/core/ignore_unused.hpp>
	81	#include <boost/geometry/util/math.hpp>
	82
	83	#include <boost/geometry/srs/projections/impl/base_static.hpp>
	84	#include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
	85	#include <boost/geometry/srs/projections/impl/projects.hpp>
	86	#include <boost/geometry/srs/projections/impl/factory_entry.hpp>
	87
	88	namespace boost { namespace geometry
	89	{
	90
11fdf7f2 TL	91	namespace projections
	92	{
	93	#ifndef DOXYGEN_NO_DETAIL
	94	namespace detail { namespace qsc
	95	{
92f5a8d4 TL	96
	97	/* The six cube faces. */
	98	enum face_type {
	99	face_front = 0,
	100	face_right = 1,
	101	face_back = 2,
	102	face_left = 3,
	103	face_top = 4,
	104	face_bottom = 5
	105	};
11fdf7f2 TL	106
	107	template <typename T>
	108	struct par_qsc
	109	{
92f5a8d4 TL	110	T a_squared;
	111	T b;
	112	T one_minus_f;
	113	T one_minus_f_squared;
	114	face_type face;
11fdf7f2 TL	115	};
11fdf7f2 TL	116
92f5a8d4	117	static const double epsilon10 = 1.e-10;
11fdf7f2 TL	118
	119	/* The four areas on a cube face. AREA_0 is the area of definition,
	120	* the other three areas are counted counterclockwise. */
92f5a8d4 TL	121	enum area_type {
	122	area_0 = 0,
	123	area_1 = 1,
	124	area_2 = 2,
	125	area_3 = 3
	126	};
11fdf7f2 TL	127
	128	/* Helper function for forward projection: compute the theta angle
	129	* and determine the area number. */
	130	template <typename T>
92f5a8d4	131	inline T qsc_fwd_equat_face_theta(T const& phi, T const& y, T const& x, area_type *area)
11fdf7f2	132	{
92f5a8d4 TL	133	static const T fourth_pi = detail::fourth_pi<T>();
	134	static const T half_pi = detail::half_pi<T>();
	135	static const T pi = detail::pi<T>();
	136
	137	T theta;
	138	if (phi < epsilon10) {
	139	*area = area_0;
	140	theta = 0.0;
	141	} else {
	142	theta = atan2(y, x);
	143	if (fabs(theta) <= fourth_pi) {
	144	*area = area_0;
	145	} else if (theta > fourth_pi && theta <= half_pi + fourth_pi) {
	146	*area = area_1;
	147	theta -= half_pi;
	148	} else if (theta > half_pi + fourth_pi \|\| theta <= -(half_pi + fourth_pi)) {
	149	*area = area_2;
	150	theta = (theta >= 0.0 ? theta - pi : theta + pi);
11fdf7f2	151	} else {
92f5a8d4 TL	152	*area = area_3;
92f5a8d4 TL	153	theta += half_pi;
11fdf7f2	154	}
92f5a8d4 TL	155	}
92f5a8d4 TL	156	return theta;
11fdf7f2 TL	157	}
	158
	159	/* Helper function: shift the longitude. */
	160	template <typename T>
	161	inline T qsc_shift_lon_origin(T const& lon, T const& offset)
	162	{
92f5a8d4 TL	163	static const T pi = detail::pi<T>();
92f5a8d4 TL	164	static const T two_pi = detail::two_pi<T>();
11fdf7f2 TL	165
11fdf7f2 TL	166	T slon = lon + offset;
92f5a8d4 TL	167	if (slon < -pi) {
	168	slon += two_pi;
	169	} else if (slon > +pi) {
	170	slon -= two_pi;
11fdf7f2 TL	171	}
	172	return slon;
	173	}
	174
	175	/* Forward projection, ellipsoid */
	176
92f5a8d4 TL	177	template <typename T, typename Parameters>
92f5a8d4 TL	178	struct base_qsc_ellipsoid
11fdf7f2	179	{
92f5a8d4	180	par_qsc<T> m_proj_parm;
11fdf7f2 TL	181
	182	// FORWARD(e_forward)
	183	// Project coordinates from geographic (lon, lat) to cartesian (x, y)
92f5a8d4	184	inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
11fdf7f2	185	{
92f5a8d4 TL	186	static const T fourth_pi = detail::fourth_pi<T>();
	187	static const T half_pi = detail::half_pi<T>();
	188	static const T pi = detail::pi<T>();
	189
	190	T lat, lon;
	191	T theta, phi;
	192	T t, mu; /* nu; */
	193	area_type area;
11fdf7f2 TL	194
	195	/* Convert the geodetic latitude to a geocentric latitude.
	196	* This corresponds to the shift from the ellipsoid to the sphere
	197	* described in [LK12]. */
92f5a8d4	198	if (par.es != 0.0) {
11fdf7f2 TL	199	lat = atan(this->m_proj_parm.one_minus_f_squared * tan(lp_lat));
	200	} else {
	201	lat = lp_lat;
	202	}
	203
	204	/* Convert the input lat, lon into theta, phi as used by QSC.
	205	* This depends on the cube face and the area on it.
	206	* For the top and bottom face, we can compute theta and phi
	207	* directly from phi, lam. For the other faces, we must use
	208	* unit sphere cartesian coordinates as an intermediate step. */
	209	lon = lp_lon;
92f5a8d4 TL	210	if (this->m_proj_parm.face == face_top) {
	211	phi = half_pi - lat;
	212	if (lon >= fourth_pi && lon <= half_pi + fourth_pi) {
	213	area = area_0;
	214	theta = lon - half_pi;
	215	} else if (lon > half_pi + fourth_pi \|\| lon <= -(half_pi + fourth_pi)) {
	216	area = area_1;
	217	theta = (lon > 0.0 ? lon - pi : lon + pi);
	218	} else if (lon > -(half_pi + fourth_pi) && lon <= -fourth_pi) {
	219	area = area_2;
	220	theta = lon + half_pi;
	221	} else {
	222	area = area_3;
	223	theta = lon;
	224	}
	225	} else if (this->m_proj_parm.face == face_bottom) {
	226	phi = half_pi + lat;
	227	if (lon >= fourth_pi && lon <= half_pi + fourth_pi) {
	228	area = area_0;
	229	theta = -lon + half_pi;
	230	} else if (lon < fourth_pi && lon >= -fourth_pi) {
	231	area = area_1;
	232	theta = -lon;
	233	} else if (lon < -fourth_pi && lon >= -(half_pi + fourth_pi)) {
	234	area = area_2;
	235	theta = -lon - half_pi;
	236	} else {
	237	area = area_3;
	238	theta = (lon > 0.0 ? -lon + pi : -lon - pi);
	239	}
	240	} else {
	241	T q, r, s;
	242	T sinlat, coslat;
	243	T sinlon, coslon;
	244
	245	if (this->m_proj_parm.face == face_right) {
	246	lon = qsc_shift_lon_origin(lon, +half_pi);
	247	} else if (this->m_proj_parm.face == face_back) {
	248	lon = qsc_shift_lon_origin(lon, +pi);
	249	} else if (this->m_proj_parm.face == face_left) {
	250	lon = qsc_shift_lon_origin(lon, -half_pi);
11fdf7f2 TL	251	}
	252	sinlat = sin(lat);
	253	coslat = cos(lat);
	254	sinlon = sin(lon);
	255	coslon = cos(lon);
	256	q = coslat * coslon;
	257	r = coslat * sinlon;
	258	s = sinlat;
92f5a8d4 TL	259
	260	if (this->m_proj_parm.face == face_front) {
	261	phi = acos(q);
	262	theta = qsc_fwd_equat_face_theta(phi, s, r, &area);
	263	} else if (this->m_proj_parm.face == face_right) {
	264	phi = acos(r);
	265	theta = qsc_fwd_equat_face_theta(phi, s, -q, &area);
	266	} else if (this->m_proj_parm.face == face_back) {
	267	phi = acos(-q);
	268	theta = qsc_fwd_equat_face_theta(phi, s, -r, &area);
	269	} else if (this->m_proj_parm.face == face_left) {
	270	phi = acos(-r);
	271	theta = qsc_fwd_equat_face_theta(phi, s, q, &area);
11fdf7f2	272	} else {
92f5a8d4 TL	273	/* Impossible */
	274	phi = theta = 0.0;
	275	area = area_0;
11fdf7f2 TL	276	}
	277	}
	278
	279	/* Compute mu and nu for the area of definition.
	280	* For mu, see Eq. (3-21) in [OL76], but note the typos:
	281	* compare with Eq. (3-14). For nu, see Eq. (3-38). */
92f5a8d4 TL	282	mu = atan((12.0 / pi) * (theta + acos(sin(theta) * cos(fourth_pi)) - half_pi));
92f5a8d4 TL	283	// TODO: (cos(mu) * cos(mu)) could be replaced with sqr(cos(mu))
11fdf7f2 TL	284	t = sqrt((1.0 - cos(phi)) / (cos(mu) * cos(mu)) / (1.0 - cos(atan(1.0 / cos(theta)))));
	285	/* nu = atan(t); We don't really need nu, just t, see below. */
	286
	287	/* Apply the result to the real area. */
92f5a8d4 TL	288	if (area == area_1) {
	289	mu += half_pi;
	290	} else if (area == area_2) {
	291	mu += pi;
	292	} else if (area == area_3) {
	293	mu += half_pi + pi;
11fdf7f2 TL	294	}
	295
	296	/* Now compute x, y from mu and nu */
	297	/* t = tan(nu); */
	298	xy_x = t * cos(mu);
	299	xy_y = t * sin(mu);
11fdf7f2 TL	300	}
	301	/* Inverse projection, ellipsoid */
	302
	303	// INVERSE(e_inverse)
	304	// Project coordinates from cartesian (x, y) to geographic (lon, lat)
92f5a8d4	305	inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
11fdf7f2	306	{
92f5a8d4 TL	307	static const T half_pi = detail::half_pi<T>();
92f5a8d4 TL	308	static const T pi = detail::pi<T>();
11fdf7f2	309
92f5a8d4 TL	310	T mu, nu, cosmu, tannu;
	311	T tantheta, theta, cosphi, phi;
	312	T t;
11fdf7f2 TL	313	int area;
	314
	315	/* Convert the input x, y to the mu and nu angles as used by QSC.
	316	* This depends on the area of the cube face. */
	317	nu = atan(sqrt(xy_x * xy_x + xy_y * xy_y));
	318	mu = atan2(xy_y, xy_x);
	319	if (xy_x >= 0.0 && xy_x >= fabs(xy_y)) {
92f5a8d4	320	area = area_0;
11fdf7f2	321	} else if (xy_y >= 0.0 && xy_y >= fabs(xy_x)) {
92f5a8d4 TL	322	area = area_1;
92f5a8d4 TL	323	mu -= half_pi;
11fdf7f2	324	} else if (xy_x < 0.0 && -xy_x >= fabs(xy_y)) {
92f5a8d4 TL	325	area = area_2;
92f5a8d4 TL	326	mu = (mu < 0.0 ? mu + pi : mu - pi);
11fdf7f2	327	} else {
92f5a8d4 TL	328	area = area_3;
92f5a8d4 TL	329	mu += half_pi;
11fdf7f2 TL	330	}
	331
	332	/* Compute phi and theta for the area of definition.
	333	* The inverse projection is not described in the original paper, but some
	334	* good hints can be found here (as of 2011-12-14):
	335	* http://fits.gsfc.nasa.gov/fitsbits/saf.93/saf.9302
	336	* (search for "Message-Id: <9302181759.AA25477 at fits.cv.nrao.edu>") */
92f5a8d4	337	t = (pi / 12.0) * tan(mu);
11fdf7f2 TL	338	tantheta = sin(t) / (cos(t) - (1.0 / sqrt(2.0)));
	339	theta = atan(tantheta);
	340	cosmu = cos(mu);
	341	tannu = tan(nu);
	342	cosphi = 1.0 - cosmu * cosmu * tannu * tannu * (1.0 - cos(atan(1.0 / cos(theta))));
	343	if (cosphi < -1.0) {
	344	cosphi = -1.0;
	345	} else if (cosphi > +1.0) {
	346	cosphi = +1.0;
	347	}
	348
	349	/* Apply the result to the real area on the cube face.
	350	* For the top and bottom face, we can compute phi and lam directly.
	351	* For the other faces, we must use unit sphere cartesian coordinates
	352	* as an intermediate step. */
92f5a8d4	353	if (this->m_proj_parm.face == face_top) {
11fdf7f2	354	phi = acos(cosphi);
92f5a8d4 TL	355	lp_lat = half_pi - phi;
	356	if (area == area_0) {
	357	lp_lon = theta + half_pi;
	358	} else if (area == area_1) {
	359	lp_lon = (theta < 0.0 ? theta + pi : theta - pi);
	360	} else if (area == area_2) {
	361	lp_lon = theta - half_pi;
11fdf7f2 TL	362	} else /* area == AREA_3 */ {
	363	lp_lon = theta;
	364	}
92f5a8d4	365	} else if (this->m_proj_parm.face == face_bottom) {
11fdf7f2	366	phi = acos(cosphi);
92f5a8d4 TL	367	lp_lat = phi - half_pi;
	368	if (area == area_0) {
	369	lp_lon = -theta + half_pi;
	370	} else if (area == area_1) {
11fdf7f2	371	lp_lon = -theta;
92f5a8d4 TL	372	} else if (area == area_2) {
	373	lp_lon = -theta - half_pi;
	374	} else /* area == area_3 */ {
	375	lp_lon = (theta < 0.0 ? -theta - pi : -theta + pi);
11fdf7f2 TL	376	}
	377	} else {
	378	/* Compute phi and lam via cartesian unit sphere coordinates. */
92f5a8d4	379	T q, r, s;
11fdf7f2 TL	380	q = cosphi;
	381	t = q * q;
	382	if (t >= 1.0) {
	383	s = 0.0;
	384	} else {
	385	s = sqrt(1.0 - t) * sin(theta);
	386	}
	387	t += s * s;
	388	if (t >= 1.0) {
	389	r = 0.0;
	390	} else {
	391	r = sqrt(1.0 - t);
	392	}
	393	/* Rotate q,r,s into the correct area. */
92f5a8d4	394	if (area == area_1) {
11fdf7f2 TL	395	t = r;
	396	r = -s;
	397	s = t;
92f5a8d4	398	} else if (area == area_2) {
11fdf7f2 TL	399	r = -r;
11fdf7f2 TL	400	s = -s;
92f5a8d4	401	} else if (area == area_3) {
11fdf7f2 TL	402	t = r;
	403	r = s;
	404	s = -t;
	405	}
	406	/* Rotate q,r,s into the correct cube face. */
92f5a8d4	407	if (this->m_proj_parm.face == face_right) {
11fdf7f2 TL	408	t = q;
	409	q = -r;
	410	r = t;
92f5a8d4	411	} else if (this->m_proj_parm.face == face_back) {
11fdf7f2 TL	412	q = -q;
11fdf7f2 TL	413	r = -r;
92f5a8d4	414	} else if (this->m_proj_parm.face == face_left) {
11fdf7f2 TL	415	t = q;
	416	q = r;
	417	r = -t;
	418	}
	419	/* Now compute phi and lam from the unit sphere coordinates. */
92f5a8d4	420	lp_lat = acos(-s) - half_pi;
11fdf7f2	421	lp_lon = atan2(r, q);
92f5a8d4 TL	422	if (this->m_proj_parm.face == face_right) {
	423	lp_lon = qsc_shift_lon_origin(lp_lon, -half_pi);
	424	} else if (this->m_proj_parm.face == face_back) {
	425	lp_lon = qsc_shift_lon_origin(lp_lon, -pi);
	426	} else if (this->m_proj_parm.face == face_left) {
	427	lp_lon = qsc_shift_lon_origin(lp_lon, +half_pi);
11fdf7f2 TL	428	}
	429	}
	430
	431	/* Apply the shift from the sphere to the ellipsoid as described
	432	* in [LK12]. */
92f5a8d4	433	if (par.es != 0.0) {
11fdf7f2	434	int invert_sign;
92f5a8d4	435	T tanphi, xa;
11fdf7f2 TL	436	invert_sign = (lp_lat < 0.0 ? 1 : 0);
	437	tanphi = tan(lp_lat);
	438	xa = this->m_proj_parm.b / sqrt(tanphi * tanphi + this->m_proj_parm.one_minus_f_squared);
92f5a8d4	439	lp_lat = atan(sqrt(par.a * par.a - xa * xa) / (this->m_proj_parm.one_minus_f * xa));
11fdf7f2 TL	440	if (invert_sign) {
	441	lp_lat = -lp_lat;
	442	}
	443	}
	444	}
	445
	446	static inline std::string get_name()
	447	{
	448	return "qsc_ellipsoid";
	449	}
	450
	451	};
	452
	453	// Quadrilateralized Spherical Cube
	454	template <typename Parameters, typename T>
92f5a8d4	455	inline void setup_qsc(Parameters const& par, par_qsc<T>& proj_parm)
11fdf7f2	456	{
92f5a8d4 TL	457	static const T fourth_pi = detail::fourth_pi<T>();
	458	static const T half_pi = detail::half_pi<T>();
	459
	460	/* Determine the cube face from the center of projection. */
	461	if (par.phi0 >= half_pi - fourth_pi / 2.0) {
	462	proj_parm.face = face_top;
	463	} else if (par.phi0 <= -(half_pi - fourth_pi / 2.0)) {
	464	proj_parm.face = face_bottom;
	465	} else if (fabs(par.lam0) <= fourth_pi) {
	466	proj_parm.face = face_front;
	467	} else if (fabs(par.lam0) <= half_pi + fourth_pi) {
	468	proj_parm.face = (par.lam0 > 0.0 ? face_right : face_left);
	469	} else {
	470	proj_parm.face = face_back;
	471	}
	472	/* Fill in useful values for the ellipsoid <-> sphere shift
	473	* described in [LK12]. */
	474	if (par.es != 0.0) {
	475	proj_parm.a_squared = par.a * par.a;
	476	proj_parm.b = par.a * sqrt(1.0 - par.es);
	477	proj_parm.one_minus_f = 1.0 - (par.a - proj_parm.b) / par.a;
	478	proj_parm.one_minus_f_squared = proj_parm.one_minus_f * proj_parm.one_minus_f;
	479	}
11fdf7f2 TL	480	}
	481
	482	}} // namespace detail::qsc
	483	#endif // doxygen
	484
	485	/*!
	486	\brief Quadrilateralized Spherical Cube projection
	487	\ingroup projections
	488	\tparam Geographic latlong point type
	489	\tparam Cartesian xy point type
	490	\tparam Parameters parameter type
	491	\par Projection characteristics
	492	- Azimuthal
	493	- Spheroid
	494	\par Example
	495	\image html ex_qsc.gif
	496	*/
92f5a8d4 TL	497	template <typename T, typename Parameters>
92f5a8d4 TL	498	struct qsc_ellipsoid : public detail::qsc::base_qsc_ellipsoid<T, Parameters>
11fdf7f2	499	{
92f5a8d4 TL	500	template <typename Params>
92f5a8d4 TL	501	inline qsc_ellipsoid(Params const& , Parameters const& par)
11fdf7f2	502	{
92f5a8d4	503	detail::qsc::setup_qsc(par, this->m_proj_parm);
11fdf7f2 TL	504	}
	505	};
	506
	507	#ifndef DOXYGEN_NO_DETAIL
	508	namespace detail
	509	{
	510
	511	// Static projection
92f5a8d4	512	BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_qsc, qsc_ellipsoid)
11fdf7f2 TL	513
11fdf7f2 TL	514	// Factory entry(s)
92f5a8d4	515	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(qsc_entry, qsc_ellipsoid)
11fdf7f2	516
92f5a8d4	517	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(qsc_init)
11fdf7f2	518	{
92f5a8d4	519	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(qsc, qsc_entry)
11fdf7f2 TL	520	}
	521
	522	} // namespace detail
	523	#endif // doxygen
	524
	525	} // namespace projections
	526
	527	}} // namespace boost::geometry
	528
	529	#endif // BOOST_GEOMETRY_PROJECTIONS_QSC_HPP
	530