/src/PROJ/src/projections/s2.cpp

Source
/******************************************************************************
 * Project:  PROJ
 * Purpose:  Implementing the S2 family of projections in PROJ
 * Author:   Marcus Elia, <marcus at geopi.pe>
 *
 ******************************************************************************
 * Copyright (c) 2021, Marcus Elia, <marcus at geopi.pe>
 *
 * 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.
 *****************************************************************************
 *
 * This implements the S2 projection.  This code is similar
 * to the QSC projection.
 *
 *
 * 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!
 *
 * You must also choose which conversion from UV to ST coordinates
 * is used (linear, quadratic, tangent). Read about them in
 * https://github.com/google/s2geometry/blob/0c4c460bdfe696da303641771f9def900b3e440f/src/s2/s2coords.h
 * The S2 projection functions are adapted from the above link and the S2
 * Math util functions are adapted from
 * https://github.com/google/s2geometry/blob/0c4c460bdfe696da303641771f9def900b3e440f/src/s2/util/math/vector.h
 ****************************************************************************/

/* enable predefined math constants M_* for MS Visual Studio */
#if defined(_MSC_VER) || defined(_WIN32)
#ifndef _USE_MATH_DEFINES
#define _USE_MATH_DEFINES
#endif
#endif

#include <cmath>
#include <cstdint>
#include <errno.h>

#include "proj.h"
#include "proj_internal.h"

/* The six cube faces. */
namespace { // anonymous namespace
enum Face {
    FACE_FRONT = 0,
    FACE_RIGHT = 1,
    FACE_TOP = 2,
    FACE_BACK = 3,
    FACE_LEFT = 4,
    FACE_BOTTOM = 5
};
} // anonymous namespace

enum S2ProjectionType { Linear, Quadratic, Tangent, NoUVtoST };
static std::map<std::string, S2ProjectionType> stringToS2ProjectionType{
    {"linear", Linear},
    {"quadratic", Quadratic},
    {"tangent", Tangent},
    {"none", NoUVtoST}};

namespace { // anonymous namespace
struct pj_s2 {
    enum Face face;
    double a_squared;
    double one_minus_f;
    double one_minus_f_squared;
    S2ProjectionType UVtoST;
};
} // anonymous namespace
PROJ_HEAD(s2, "S2") "\n\tMisc, Sph&Ell";

/* The four areas on a cube face. AREA_0 is the area of definition,
 * the other three areas are counted counterclockwise. */
namespace { // anonymous namespace
enum Area { AREA_0 = 0, AREA_1 = 1, AREA_2 = 2, AREA_3 = 3 };
} // anonymous namespace

// =================================================
//
//                  S2 Math Util
//
// =================================================
static PJ_XYZ Abs(const PJ_XYZ &p) {
    return {std::fabs(p.x), std::fabs(p.y), std::fabs(p.z)};
}
// return the index of the largest component (fabs)
static int LargestAbsComponent(const PJ_XYZ &p) {
    PJ_XYZ temp = Abs(p);
    return temp.x > temp.y ? temp.x > temp.z ? 0 : 2 : temp.y > temp.z ? 1 : 2;
}

// =================================================
//
//              S2 Projection Functions
//
// =================================================

// Unfortunately, tan(M_PI_4) is slightly less than 1.0.  This isn't due to
// a flaw in the implementation of tan(), it's because the derivative of
// tan(x) at x=pi/4 is 2, and it happens that the two adjacent floating
// point numbers on either side of the infinite-precision value of pi/4 have
// tangents that are slightly below and slightly above 1.0 when rounded to
// the nearest double-precision result.
static double STtoUV(double s, S2ProjectionType s2_projection) {
    switch (s2_projection) {
    case Linear:
        return 2 * s - 1;
        break;
    case Quadratic:
        if (s >= 0.5)
            return (1 / 3.) * (4 * s * s - 1);
        else
            return (1 / 3.) * (1 - 4 * (1 - s) * (1 - s));
        break;
    case Tangent:
        s = std::tan(M_PI_2 * s - M_PI_4);
        return s +
               (1.0 / static_cast<double>(static_cast<std::int64_t>(1) << 53)) *
                   s;
        break;
    default:
        return s;
    }
}

static double UVtoST(double u, S2ProjectionType s2_projection) {
    double ret = u;
    switch (s2_projection) {
    case Linear:
        ret = 0.5 * (u + 1);
        break;
    case Quadratic:
        if (u >= 0)
            ret = 0.5 * std::sqrt(1 + 3 * u);
        else
            ret = 1 - 0.5 * std::sqrt(1 - 3 * u);
        break;
    case Tangent: {
        volatile double a = std::atan(u);
        ret = (2 * M_1_PI) * (a + M_PI_4);
        break;
    }
    case NoUVtoST:
        break;
    }
    return ret;
}

inline PJ_XYZ FaceUVtoXYZ(int face, double u, double v) {
    switch (face) {
    case 0:
        return {1, u, v};
    case 1:
        return {-u, 1, v};
    case 2:
        return {-u, -v, 1};
    case 3:
        return {-1, -v, -u};
    case 4:
        return {v, -1, -u};
    default:
        return {v, u, -1};
    }
}

inline PJ_XYZ FaceUVtoXYZ(int face, const PJ_XY &uv) {
    return FaceUVtoXYZ(face, uv.x, uv.y);
}

inline void ValidFaceXYZtoUV(int face, const PJ_XYZ &p, double *pu,
                             double *pv) {
    switch (face) {
    case 0:
        *pu = p.y / p.x;
        *pv = p.z / p.x;
        break;
    case 1:
        *pu = -p.x / p.y;
        *pv = p.z / p.y;
        break;
    case 2:
        *pu = -p.x / p.z;
        *pv = -p.y / p.z;
        break;
    case 3:
        *pu = p.z / p.x;
        *pv = p.y / p.x;
        break;
    case 4:
        *pu = p.z / p.y;
        *pv = -p.x / p.y;
        break;
    default:
        *pu = -p.y / p.z;
        *pv = -p.x / p.z;
        break;
    }
}

inline void ValidFaceXYZtoUV(int face, const PJ_XYZ &p, PJ_XY *puv) {
    ValidFaceXYZtoUV(face, p, &(*puv).x, &(*puv).y);
}

inline int GetFace(const PJ_XYZ &p) {
    int face = LargestAbsComponent(p);
    double pFace;
    switch (face) {
    case 0:
        pFace = p.x;
        break;
    case 1:
        pFace = p.y;
        break;
    default:
        pFace = p.z;
        break;
    }
    if (pFace < 0)
        face += 3;
    return face;
}

inline int XYZtoFaceUV(const PJ_XYZ &p, double *pu, double *pv) {
    int face = GetFace(p);
    ValidFaceXYZtoUV(face, p, pu, pv);
    return face;
}

inline int XYZtoFaceUV(const PJ_XYZ &p, PJ_XY *puv) {
    return XYZtoFaceUV(p, &(*puv).x, &(*puv).y);
}

inline bool FaceXYZtoUV(int face, const PJ_XYZ &p, double *pu, double *pv) {
    double pFace;
    switch (face) {
    case 0:
        pFace = p.x;
        break;
    case 1:
        pFace = p.y;
        break;
    case 2:
        pFace = p.z;
        break;
    case 3:
        pFace = p.x;
        break;
    case 4:
        pFace = p.y;
        break;
    default:
        pFace = p.z;
        break;
    }
    if (face < 3) {
        if (pFace <= 0)
            return false;
    } else {
        if (pFace >= 0)
            return false;
    }
    ValidFaceXYZtoUV(face, p, pu, pv);
    return true;
}

inline bool FaceXYZtoUV(int face, const PJ_XYZ &p, PJ_XY *puv) {
    return FaceXYZtoUV(face, p, &(*puv).x, &(*puv).y);
}

// This function inverts ValidFaceXYZtoUV()
inline bool UVtoSphereXYZ(int face, double u, double v, PJ_XYZ *xyz) {
    double major_coord = 1 / sqrt(1 + u * u + v * v);
    double minor_coord_1 = u * major_coord;
    double minor_coord_2 = v * major_coord;

    switch (face) {
    case 0:
        xyz->x = major_coord;
        xyz->y = minor_coord_1;
        xyz->z = minor_coord_2;
        break;
    case 1:
        xyz->x = -minor_coord_1;
        xyz->y = major_coord;
        xyz->z = minor_coord_2;
        break;
    case 2:
        xyz->x = -minor_coord_1;
        xyz->y = -minor_coord_2;
        xyz->z = major_coord;
        break;
    case 3:
        xyz->x = -major_coord;
        xyz->y = -minor_coord_2;
        xyz->z = -minor_coord_1;
        break;
    case 4:
        xyz->x = minor_coord_2;
        xyz->y = -major_coord;
        xyz->z = -minor_coord_1;
        break;
    default:
        xyz->x = minor_coord_2;
        xyz->y = minor_coord_1;
        xyz->z = -major_coord;
        break;
    }

    return true;
}

// ============================================
//
//      The Forward and Inverse Functions
//
// ============================================
static PJ_XY s2_forward(PJ_LP lp, PJ *P) {
    struct pj_s2 *Q = static_cast<struct pj_s2 *>(P->opaque);
    double lat;

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

    // Convert the lat/long to x,y,z on the unit sphere
    double x, y, z;
    double sinlat, coslat;
    double sinlon, coslon;

    sinlat = sin(lat);
    coslat = cos(lat);
    sinlon = sin(lp.lam);
    coslon = cos(lp.lam);
    x = coslat * coslon;
    y = coslat * sinlon;
    z = sinlat;

    PJ_XYZ spherePoint{x, y, z};
    PJ_XY uvCoords;

    ValidFaceXYZtoUV(Q->face, spherePoint, &uvCoords.x, &uvCoords.y);
    double s = UVtoST(uvCoords.x, Q->UVtoST);
    double t = UVtoST(uvCoords.y, Q->UVtoST);
    return {s, t};
}

static PJ_LP s2_inverse(PJ_XY xy, PJ *P) {
    PJ_LP lp = {0.0, 0.0};
    struct pj_s2 *Q = static_cast<struct pj_s2 *>(P->opaque);

    // Do the S2 projections to get from s,t to u,v to x,y,z
    double u = STtoUV(xy.x, Q->UVtoST);
    double v = STtoUV(xy.y, Q->UVtoST);

    PJ_XYZ sphereCoords;
    UVtoSphereXYZ(Q->face, u, v, &sphereCoords);
    double q = sphereCoords.x;
    double r = sphereCoords.y;
    double s = sphereCoords.z;

    // Get the spherical angles from the x y z
    lp.phi = acos(-s) - M_HALFPI;
    lp.lam = atan2(r, q);

    /* Apply the shift from the sphere to the ellipsoid as described
     * in [LK12]. */
    if (P->es != 0.0) {
        int invert_sign;
        volatile double tanphi, xa;
        invert_sign = (lp.phi < 0.0 ? 1 : 0);
        tanphi = tan(lp.phi);
        xa = P->b / sqrt(tanphi * tanphi + Q->one_minus_f_squared);
        lp.phi = atan(sqrt(Q->a_squared - xa * xa) / (Q->one_minus_f * xa));
        if (invert_sign) {
            lp.phi = -lp.phi;
        }
    }

    return lp;
}

PJ *PJ_PROJECTION(s2) {
    struct pj_s2 *Q =
        static_cast<struct pj_s2 *>(calloc(1, sizeof(struct pj_s2)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;

    /* Determine which UVtoST function is to be used */
    PROJVALUE maybeUVtoST = pj_param(P->ctx, P->params, "sUVtoST");
    if (nullptr != maybeUVtoST.s) {
        try {
            Q->UVtoST = stringToS2ProjectionType.at(maybeUVtoST.s);
        } catch (const std::out_of_range &) {
            proj_log_error(
                P, _("Invalid value for s2 parameter: should be linear, "
                     "quadratic, tangent, or none."));
            return pj_default_destructor(P,
                                         PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
    } else {
        Q->UVtoST = Quadratic;
    }

    P->left = PJ_IO_UNITS_RADIANS;
    P->right = PJ_IO_UNITS_PROJECTED;
    P->from_greenwich = -P->lam0;

    P->inv = s2_inverse;
    P->fwd = s2_forward;

    /* Determine the cube face from the center of projection. */
    if (P->phi0 >= M_HALFPI - M_FORTPI / 2.0) {
        Q->face = FACE_TOP;
    } else if (P->phi0 <= -(M_HALFPI - M_FORTPI / 2.0)) {
        Q->face = FACE_BOTTOM;
    } else if (fabs(P->lam0) <= M_FORTPI) {
        Q->face = FACE_FRONT;
    } else if (fabs(P->lam0) <= M_HALFPI + M_FORTPI) {
        Q->face = (P->lam0 > 0.0 ? FACE_RIGHT : FACE_LEFT);
    } else {
        Q->face = FACE_BACK;
    }
    /* Fill in useful values for the ellipsoid <-> sphere shift
     * described in [LK12]. */
    if (P->es != 0.0) {
        Q->a_squared = P->a * P->a;
        Q->one_minus_f = 1.0 - (P->a - P->b) / P->a;
        Q->one_minus_f_squared = Q->one_minus_f * Q->one_minus_f;
    }

    return P;
}

Coverage Report

Created: 2025-10-13 06:59

Line	Count	Source
1		/******************************************************************************
2		* Project: PROJ
3		* Purpose: Implementing the S2 family of projections in PROJ
4		* Author: Marcus Elia, <marcus at geopi.pe>
5		*
6		******************************************************************************
7		* Copyright (c) 2021, Marcus Elia, <marcus at geopi.pe>
8		*
9		* Permission is hereby granted, free of charge, to any person obtaining a
10		* copy of this software and associated documentation files (the "Software"),
11		* to deal in the Software without restriction, including without limitation
12		* the rights to use, copy, modify, merge, publish, distribute, sublicense,
13		* and/or sell copies of the Software, and to permit persons to whom the
14		* Software is furnished to do so, subject to the following conditions:
15		*
16		* The above copyright notice and this permission notice shall be included
17		* in all copies or substantial portions of the Software.
18		*
19		* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
20		* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21		* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
22		* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23		* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
24		* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
25		* DEALINGS IN THE SOFTWARE.
26		*****************************************************************************
27		*
28		* This implements the S2 projection. This code is similar
29		* to the QSC projection.
30		*
31		*
32		* You have to choose one of the following projection centers,
33		* corresponding to the centers of the six cube faces:
34		* phi0 = 0.0, lam0 = 0.0 ("front" face)
35		* phi0 = 0.0, lam0 = 90.0 ("right" face)
36		* phi0 = 0.0, lam0 = 180.0 ("back" face)
37		* phi0 = 0.0, lam0 = -90.0 ("left" face)
38		* phi0 = 90.0 ("top" face)
39		* phi0 = -90.0 ("bottom" face)
40		* Other projection centers will not work!
41		*
42		* You must also choose which conversion from UV to ST coordinates
43		* is used (linear, quadratic, tangent). Read about them in
44		* https://github.com/google/s2geometry/blob/0c4c460bdfe696da303641771f9def900b3e440f/src/s2/s2coords.h
45		* The S2 projection functions are adapted from the above link and the S2
46		* Math util functions are adapted from
47		* https://github.com/google/s2geometry/blob/0c4c460bdfe696da303641771f9def900b3e440f/src/s2/util/math/vector.h
48		****************************************************************************/
49
50		/* enable predefined math constants M_* for MS Visual Studio */
51		#if defined(_MSC_VER) \|\| defined(_WIN32)
52		#ifndef _USE_MATH_DEFINES
53		#define _USE_MATH_DEFINES
54		#endif
55		#endif
56
57		#include <cmath>
58		#include <cstdint>
59		#include <errno.h>
60
61		#include "proj.h"
62		#include "proj_internal.h"
63
64		/* The six cube faces. */
65		namespace { // anonymous namespace
66		enum Face {
67		FACE_FRONT = 0,
68		FACE_RIGHT = 1,
69		FACE_TOP = 2,
70		FACE_BACK = 3,
71		FACE_LEFT = 4,
72		FACE_BOTTOM = 5
73		};
74		} // anonymous namespace
75
76		enum S2ProjectionType { Linear, Quadratic, Tangent, NoUVtoST };
77		static std::map<std::string, S2ProjectionType> stringToS2ProjectionType{
78		{"linear", Linear},
79		{"quadratic", Quadratic},
80		{"tangent", Tangent},
81		{"none", NoUVtoST}};
82
83		namespace { // anonymous namespace
84		struct pj_s2 {
85		enum Face face;
86		double a_squared;
87		double one_minus_f;
88		double one_minus_f_squared;
89		S2ProjectionType UVtoST;
90		};
91		} // anonymous namespace
92		PROJ_HEAD(s2, "S2") "\n\tMisc, Sph&Ell";
93
94		/* The four areas on a cube face. AREA_0 is the area of definition,
95		* the other three areas are counted counterclockwise. */
96		namespace { // anonymous namespace
97		enum Area { AREA_0 = 0, AREA_1 = 1, AREA_2 = 2, AREA_3 = 3 };
98		} // anonymous namespace
99
100		// =================================================
101		//
102		// S2 Math Util
103		//
104		// =================================================
105	0	static PJ_XYZ Abs(const PJ_XYZ &p) {
106	0	return {std::fabs(p.x), std::fabs(p.y), std::fabs(p.z)};
107	0	}
108		// return the index of the largest component (fabs)
109	0	static int LargestAbsComponent(const PJ_XYZ &p) {
110	0	PJ_XYZ temp = Abs(p);
111	0	return temp.x > temp.y ? temp.x > temp.z ? 0 : 2 : temp.y > temp.z ? 1 : 2;
112	0	}
113
114		// =================================================
115		//
116		// S2 Projection Functions
117		//
118		// =================================================
119
120		// Unfortunately, tan(M_PI_4) is slightly less than 1.0. This isn't due to
121		// a flaw in the implementation of tan(), it's because the derivative of
122		// tan(x) at x=pi/4 is 2, and it happens that the two adjacent floating
123		// point numbers on either side of the infinite-precision value of pi/4 have
124		// tangents that are slightly below and slightly above 1.0 when rounded to
125		// the nearest double-precision result.
126	0	static double STtoUV(double s, S2ProjectionType s2_projection) {
127	0	switch (s2_projection) {
128	0	case Linear:
129	0	return 2 * s - 1;
130	0	break;
131	0	case Quadratic:
132	0	if (s >= 0.5)
133	0	return (1 / 3.) * (4 * s * s - 1);
134	0	else
135	0	return (1 / 3.) * (1 - 4 * (1 - s) * (1 - s));
136	0	break;
137	0	case Tangent:
138	0	s = std::tan(M_PI_2 * s - M_PI_4);
139	0	return s +
140	0	(1.0 / static_cast<double>(static_cast<std::int64_t>(1) << 53)) *
141	0	s;
142	0	break;
143	0	default:
144	0	return s;
145	0	}
146	0	}
147
148	28.7k	static double UVtoST(double u, S2ProjectionType s2_projection) {
149	28.7k	double ret = u;
150	28.7k	switch (s2_projection) {
151	0	case Linear:
152	0	ret = 0.5 * (u + 1);
153	0	break;
154	28.7k	case Quadratic:
155	28.7k	if (u >= 0)
156	14.3k	ret = 0.5 * std::sqrt(1 + 3 * u);
157	14.4k	else
158	14.4k	ret = 1 - 0.5 * std::sqrt(1 - 3 * u);
159	28.7k	break;
160	0	case Tangent: {
161	0	volatile double a = std::atan(u);
162	0	ret = (2 * M_1_PI) * (a + M_PI_4);
163	0	break;
164	0	}
165	0	case NoUVtoST:
166	0	break;
167	28.7k	}
168	28.7k	return ret;
169	28.7k	}
170
171	0	inline PJ_XYZ FaceUVtoXYZ(int face, double u, double v) {
172	0	switch (face) {
173	0	case 0:
174	0	return {1, u, v};
175	0	case 1:
176	0	return {-u, 1, v};
177	0	case 2:
178	0	return {-u, -v, 1};
179	0	case 3:
180	0	return {-1, -v, -u};
181	0	case 4:
182	0	return {v, -1, -u};
183	0	default:
184	0	return {v, u, -1};
185	0	}
186	0	}
187
188	0	inline PJ_XYZ FaceUVtoXYZ(int face, const PJ_XY &uv) {
189	0	return FaceUVtoXYZ(face, uv.x, uv.y);
190	0	}
191
192		inline void ValidFaceXYZtoUV(int face, const PJ_XYZ &p, double *pu,
193	14.3k	double *pv) {
194	14.3k	switch (face) {
195	14.3k	case 0:
196	14.3k	*pu = p.y / p.x;
197	14.3k	*pv = p.z / p.x;
198	14.3k	break;
199	0	case 1:
200	0	*pu = -p.x / p.y;
201	0	*pv = p.z / p.y;
202	0	break;
203	0	case 2:
204	0	*pu = -p.x / p.z;
205	0	*pv = -p.y / p.z;
206	0	break;
207	0	case 3:
208	0	*pu = p.z / p.x;
209	0	*pv = p.y / p.x;
210	0	break;
211	0	case 4:
212	0	*pu = p.z / p.y;
213	0	*pv = -p.x / p.y;
214	0	break;
215	0	default:
216	0	*pu = -p.y / p.z;
217	0	*pv = -p.x / p.z;
218	0	break;
219	14.3k	}
220	14.3k	}
221
222	0	inline void ValidFaceXYZtoUV(int face, const PJ_XYZ &p, PJ_XY *puv) {
223	0	ValidFaceXYZtoUV(face, p, &(puv).x, &(puv).y);
224	0	}
225
226	0	inline int GetFace(const PJ_XYZ &p) {
227	0	int face = LargestAbsComponent(p);
228	0	double pFace;
229	0	switch (face) {
230	0	case 0:
231	0	pFace = p.x;
232	0	break;
233	0	case 1:
234	0	pFace = p.y;
235	0	break;
236	0	default:
237	0	pFace = p.z;
238	0	break;
239	0	}
240	0	if (pFace < 0)
241	0	face += 3;
242	0	return face;
243	0	}
244
245	0	inline int XYZtoFaceUV(const PJ_XYZ &p, double pu, double pv) {
246	0	int face = GetFace(p);
247	0	ValidFaceXYZtoUV(face, p, pu, pv);
248	0	return face;
249	0	}
250
251	0	inline int XYZtoFaceUV(const PJ_XYZ &p, PJ_XY *puv) {
252	0	return XYZtoFaceUV(p, &(puv).x, &(puv).y);
253	0	}
254
255	0	inline bool FaceXYZtoUV(int face, const PJ_XYZ &p, double pu, double pv) {
256	0	double pFace;
257	0	switch (face) {
258	0	case 0:
259	0	pFace = p.x;
260	0	break;
261	0	case 1:
262	0	pFace = p.y;
263	0	break;
264	0	case 2:
265	0	pFace = p.z;
266	0	break;
267	0	case 3:
268	0	pFace = p.x;
269	0	break;
270	0	case 4:
271	0	pFace = p.y;
272	0	break;
273	0	default:
274	0	pFace = p.z;
275	0	break;
276	0	}
277	0	if (face < 3) {
278	0	if (pFace <= 0)
279	0	return false;
280	0	} else {
281	0	if (pFace >= 0)
282	0	return false;
283	0	}
284	0	ValidFaceXYZtoUV(face, p, pu, pv);
285	0	return true;
286	0	}
287
288	0	inline bool FaceXYZtoUV(int face, const PJ_XYZ &p, PJ_XY *puv) {
289	0	return FaceXYZtoUV(face, p, &(puv).x, &(puv).y);
290	0	}
291
292		// This function inverts ValidFaceXYZtoUV()
293	0	inline bool UVtoSphereXYZ(int face, double u, double v, PJ_XYZ *xyz) {
294	0	double major_coord = 1 / sqrt(1 + u * u + v * v);
295	0	double minor_coord_1 = u * major_coord;
296	0	double minor_coord_2 = v * major_coord;
297
298	0	switch (face) {
299	0	case 0:
300	0	xyz->x = major_coord;
301	0	xyz->y = minor_coord_1;
302	0	xyz->z = minor_coord_2;
303	0	break;
304	0	case 1:
305	0	xyz->x = -minor_coord_1;
306	0	xyz->y = major_coord;
307	0	xyz->z = minor_coord_2;
308	0	break;
309	0	case 2:
310	0	xyz->x = -minor_coord_1;
311	0	xyz->y = -minor_coord_2;
312	0	xyz->z = major_coord;
313	0	break;
314	0	case 3:
315	0	xyz->x = -major_coord;
316	0	xyz->y = -minor_coord_2;
317	0	xyz->z = -minor_coord_1;
318	0	break;
319	0	case 4:
320	0	xyz->x = minor_coord_2;
321	0	xyz->y = -major_coord;
322	0	xyz->z = -minor_coord_1;
323	0	break;
324	0	default:
325	0	xyz->x = minor_coord_2;
326	0	xyz->y = minor_coord_1;
327	0	xyz->z = -major_coord;
328	0	break;
329	0	}
330
331	0	return true;
332	0	}
333
334		// ============================================
335		//
336		// The Forward and Inverse Functions
337		//
338		// ============================================
339	14.3k	static PJ_XY s2_forward(PJ_LP lp, PJ *P) {
340	14.3k	struct pj_s2 Q = static_cast<struct pj_s2 >(P->opaque);
341	14.3k	double lat;
342
343		/* Convert the geodetic latitude to a geocentric latitude.
344		* This corresponds to the shift from the ellipsoid to the sphere
345		* described in [LK12]. */
346	14.3k	if (P->es != 0.0) {
347	14.3k	lat = atan(Q->one_minus_f_squared * tan(lp.phi));
348	14.3k	} else {
349	0	lat = lp.phi;
350	0	}
351
352		// Convert the lat/long to x,y,z on the unit sphere
353	14.3k	double x, y, z;
354	14.3k	double sinlat, coslat;
355	14.3k	double sinlon, coslon;
356
357	14.3k	sinlat = sin(lat);
358	14.3k	coslat = cos(lat);
359	14.3k	sinlon = sin(lp.lam);
360	14.3k	coslon = cos(lp.lam);
361	14.3k	x = coslat * coslon;
362	14.3k	y = coslat * sinlon;
363	14.3k	z = sinlat;
364
365	14.3k	PJ_XYZ spherePoint{x, y, z};
366	14.3k	PJ_XY uvCoords;
367
368	14.3k	ValidFaceXYZtoUV(Q->face, spherePoint, &uvCoords.x, &uvCoords.y);
369	14.3k	double s = UVtoST(uvCoords.x, Q->UVtoST);
370	14.3k	double t = UVtoST(uvCoords.y, Q->UVtoST);
371	14.3k	return {s, t};
372	14.3k	}
373
374	0	static PJ_LP s2_inverse(PJ_XY xy, PJ *P) {
375	0	PJ_LP lp = {0.0, 0.0};
376	0	struct pj_s2 Q = static_cast<struct pj_s2 >(P->opaque);
377
378		// Do the S2 projections to get from s,t to u,v to x,y,z
379	0	double u = STtoUV(xy.x, Q->UVtoST);
380	0	double v = STtoUV(xy.y, Q->UVtoST);
381
382	0	PJ_XYZ sphereCoords;
383	0	UVtoSphereXYZ(Q->face, u, v, &sphereCoords);
384	0	double q = sphereCoords.x;
385	0	double r = sphereCoords.y;
386	0	double s = sphereCoords.z;
387
388		// Get the spherical angles from the x y z
389	0	lp.phi = acos(-s) - M_HALFPI;
390	0	lp.lam = atan2(r, q);
391
392		/* Apply the shift from the sphere to the ellipsoid as described
393		* in [LK12]. */
394	0	if (P->es != 0.0) {
395	0	int invert_sign;
396	0	volatile double tanphi, xa;
397	0	invert_sign = (lp.phi < 0.0 ? 1 : 0);
398	0	tanphi = tan(lp.phi);
399	0	xa = P->b / sqrt(tanphi * tanphi + Q->one_minus_f_squared);
400	0	lp.phi = atan(sqrt(Q->a_squared - xa * xa) / (Q->one_minus_f * xa));
401	0	if (invert_sign) {
402	0	lp.phi = -lp.phi;
403	0	}
404	0	}
405
406	0	return lp;
407	0	}
408
409	445	PJ *PJ_PROJECTION(s2) {
410	445	struct pj_s2 *Q =
411	445	static_cast<struct pj_s2 *>(calloc(1, sizeof(struct pj_s2)));
412	445	if (nullptr == Q)
413	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
414	445	P->opaque = Q;
415
416		/* Determine which UVtoST function is to be used */
417	445	PROJVALUE maybeUVtoST = pj_param(P->ctx, P->params, "sUVtoST");
418	445	if (nullptr != maybeUVtoST.s) {
419	32	try {
420	32	Q->UVtoST = stringToS2ProjectionType.at(maybeUVtoST.s);
421	32	} catch (const std::out_of_range &) {
422	29	proj_log_error(
423	29	P, _("Invalid value for s2 parameter: should be linear, "
424	29	"quadratic, tangent, or none."));
425	29	return pj_default_destructor(P,
426	29	PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
427	29	}
428	413	} else {
429	413	Q->UVtoST = Quadratic;
430	413	}
431
432	416	P->left = PJ_IO_UNITS_RADIANS;
433	416	P->right = PJ_IO_UNITS_PROJECTED;
434	416	P->from_greenwich = -P->lam0;
435
436	416	P->inv = s2_inverse;
437	416	P->fwd = s2_forward;
438
439		/* Determine the cube face from the center of projection. */
440	416	if (P->phi0 >= M_HALFPI - M_FORTPI / 2.0) {
441	11	Q->face = FACE_TOP;
442	405	} else if (P->phi0 <= -(M_HALFPI - M_FORTPI / 2.0)) {
443	37	Q->face = FACE_BOTTOM;
444	368	} else if (fabs(P->lam0) <= M_FORTPI) {
445	315	Q->face = FACE_FRONT;
446	315	} else if (fabs(P->lam0) <= M_HALFPI + M_FORTPI) {
447	26	Q->face = (P->lam0 > 0.0 ? FACE_RIGHT : FACE_LEFT);
448	27	} else {
449	27	Q->face = FACE_BACK;
450	27	}
451		/* Fill in useful values for the ellipsoid <-> sphere shift
452		* described in [LK12]. */
453	416	if (P->es != 0.0) {
454	316	Q->a_squared = P->a * P->a;
455	316	Q->one_minus_f = 1.0 - (P->a - P->b) / P->a;
456	316	Q->one_minus_f_squared = Q->one_minus_f * Q->one_minus_f;
457	316	}
458
459	416	return P;
460	445	}