/src/PROJ/src/projections/qsc.cpp

Source (jump to first uncovered line)
/*
 * This implements the Quadrilateralized Spherical Cube (QSC) projection.
 *
 * Copyright (c) 2011, 2012  Martin Lambers <marlam@marlam.de>
 *
 * 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. Pacific 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 PROJECTION(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 pj_qsc_ns::AREA_0 in the projection code below. The
 * other three areas of a cube face are handled by rotation of
 * pj_qsc_ns::AREA_0.
 */

#include <errno.h>
#include <math.h>

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

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

namespace { // anonymous namespace
struct pj_qsc_data {
    enum pj_qsc_ns::Face face;
    double a_squared;
    double b;
    double one_minus_f;
    double one_minus_f_squared;
};
} // anonymous namespace
PROJ_HEAD(qsc, "Quadrilateralized Spherical Cube") "\n\tAzi, Sph";

#define EPS10 1.e-10

/* The four areas on a cube face. AREA_0 is the area of definition,
 * the other three areas are counted counterclockwise. */
namespace pj_qsc_ns {
enum Area { 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. */
static double qsc_fwd_equat_face_theta(double phi, double y, double x,
                                       enum pj_qsc_ns::Area *area) {
    double theta;
    if (phi < EPS10) {
        *area = pj_qsc_ns::AREA_0;
        theta = 0.0;
    } else {
        theta = atan2(y, x);
        if (fabs(theta) <= M_FORTPI) {
            *area = pj_qsc_ns::AREA_0;
        } else if (theta > M_FORTPI && theta <= M_HALFPI + M_FORTPI) {
            *area = pj_qsc_ns::AREA_1;
            theta -= M_HALFPI;
        } else if (theta > M_HALFPI + M_FORTPI ||
                   theta <= -(M_HALFPI + M_FORTPI)) {
            *area = pj_qsc_ns::AREA_2;
            theta = (theta >= 0.0 ? theta - M_PI : theta + M_PI);
        } else {
            *area = pj_qsc_ns::AREA_3;
            theta += M_HALFPI;
        }
    }
    return theta;
}

/* Helper function: shift the longitude. */
static double qsc_shift_longitude_origin(double longitude, double offset) {
    double slon = longitude + offset;
    if (slon < -M_PI) {
        slon += M_TWOPI;
    } else if (slon > +M_PI) {
        slon -= M_TWOPI;
    }
    return slon;
}

static PJ_XY qsc_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_qsc_data *Q = static_cast<struct pj_qsc_data *>(P->opaque);
    double lat, longitude;
    double theta, phi;
    double t, mu; /* nu; */
    enum pj_qsc_ns::Area area;

    /* 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 input lat, longitude 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. */
    longitude = lp.lam;
    if (Q->face == pj_qsc_ns::FACE_TOP) {
        phi = M_HALFPI - lat;
        if (longitude >= M_FORTPI && longitude <= M_HALFPI + M_FORTPI) {
            area = pj_qsc_ns::AREA_0;
            theta = longitude - M_HALFPI;
        } else if (longitude > M_HALFPI + M_FORTPI ||
                   longitude <= -(M_HALFPI + M_FORTPI)) {
            area = pj_qsc_ns::AREA_1;
            theta = (longitude > 0.0 ? longitude - M_PI : longitude + M_PI);
        } else if (longitude > -(M_HALFPI + M_FORTPI) &&
                   longitude <= -M_FORTPI) {
            area = pj_qsc_ns::AREA_2;
            theta = longitude + M_HALFPI;
        } else {
            area = pj_qsc_ns::AREA_3;
            theta = longitude;
        }
    } else if (Q->face == pj_qsc_ns::FACE_BOTTOM) {
        phi = M_HALFPI + lat;
        if (longitude >= M_FORTPI && longitude <= M_HALFPI + M_FORTPI) {
            area = pj_qsc_ns::AREA_0;
            theta = -longitude + M_HALFPI;
        } else if (longitude < M_FORTPI && longitude >= -M_FORTPI) {
            area = pj_qsc_ns::AREA_1;
            theta = -longitude;
        } else if (longitude < -M_FORTPI &&
                   longitude >= -(M_HALFPI + M_FORTPI)) {
            area = pj_qsc_ns::AREA_2;
            theta = -longitude - M_HALFPI;
        } else {
            area = pj_qsc_ns::AREA_3;
            theta = (longitude > 0.0 ? -longitude + M_PI : -longitude - M_PI);
        }
    } else {
        double q, r, s;
        double sinlat, coslat;
        double sinlon, coslon;

        if (Q->face == pj_qsc_ns::FACE_RIGHT) {
            longitude = qsc_shift_longitude_origin(longitude, +M_HALFPI);
        } else if (Q->face == pj_qsc_ns::FACE_BACK) {
            longitude = qsc_shift_longitude_origin(longitude, +M_PI);
        } else if (Q->face == pj_qsc_ns::FACE_LEFT) {
            longitude = qsc_shift_longitude_origin(longitude, -M_HALFPI);
        }
        sinlat = sin(lat);
        coslat = cos(lat);
        sinlon = sin(longitude);
        coslon = cos(longitude);
        q = coslat * coslon;
        r = coslat * sinlon;
        s = sinlat;

        if (Q->face == pj_qsc_ns::FACE_FRONT) {
            phi = acos(q);
            theta = qsc_fwd_equat_face_theta(phi, s, r, &area);
        } else if (Q->face == pj_qsc_ns::FACE_RIGHT) {
            phi = acos(r);
            theta = qsc_fwd_equat_face_theta(phi, s, -q, &area);
        } else if (Q->face == pj_qsc_ns::FACE_BACK) {
            phi = acos(-q);
            theta = qsc_fwd_equat_face_theta(phi, s, -r, &area);
        } else if (Q->face == pj_qsc_ns::FACE_LEFT) {
            phi = acos(-r);
            theta = qsc_fwd_equat_face_theta(phi, s, q, &area);
        } else {
            /* Impossible */
            phi = theta = 0.0;
            area = pj_qsc_ns::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 / M_PI) *
              (theta + acos(sin(theta) * cos(M_FORTPI)) - M_HALFPI));
    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 == pj_qsc_ns::AREA_1) {
        mu += M_HALFPI;
    } else if (area == pj_qsc_ns::AREA_2) {
        mu += M_PI;
    } else if (area == pj_qsc_ns::AREA_3) {
        mu += M_PI_HALFPI;
    }

    /* Now compute x, y from mu and nu */
    /* t = tan(nu); */
    xy.x = t * cos(mu);
    xy.y = t * sin(mu);
    return xy;
}

static PJ_LP qsc_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_qsc_data *Q = static_cast<struct pj_qsc_data *>(P->opaque);
    double mu, nu, cosmu, tannu;
    double tantheta, theta, cosphi, phi;
    double 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 = pj_qsc_ns::AREA_0;
    } else if (xy.y >= 0.0 && xy.y >= fabs(xy.x)) {
        area = pj_qsc_ns::AREA_1;
        mu -= M_HALFPI;
    } else if (xy.x < 0.0 && -xy.x >= fabs(xy.y)) {
        area = pj_qsc_ns::AREA_2;
        mu = (mu < 0.0 ? mu + M_PI : mu - M_PI);
    } else {
        area = pj_qsc_ns::AREA_3;
        mu += M_HALFPI;
    }

    /* 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 = (M_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 (Q->face == pj_qsc_ns::FACE_TOP) {
        phi = acos(cosphi);
        lp.phi = M_HALFPI - phi;
        if (area == pj_qsc_ns::AREA_0) {
            lp.lam = theta + M_HALFPI;
        } else if (area == pj_qsc_ns::AREA_1) {
            lp.lam = (theta < 0.0 ? theta + M_PI : theta - M_PI);
        } else if (area == pj_qsc_ns::AREA_2) {
            lp.lam = theta - M_HALFPI;
        } else /* area == pj_qsc_ns::AREA_3 */ {
            lp.lam = theta;
        }
    } else if (Q->face == pj_qsc_ns::FACE_BOTTOM) {
        phi = acos(cosphi);
        lp.phi = phi - M_HALFPI;
        if (area == pj_qsc_ns::AREA_0) {
            lp.lam = -theta + M_HALFPI;
        } else if (area == pj_qsc_ns::AREA_1) {
            lp.lam = -theta;
        } else if (area == pj_qsc_ns::AREA_2) {
            lp.lam = -theta - M_HALFPI;
        } else /* area == pj_qsc_ns::AREA_3 */ {
            lp.lam = (theta < 0.0 ? -theta - M_PI : -theta + M_PI);
        }
    } else {
        /* Compute phi and lam via cartesian unit sphere coordinates. */
        double 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 == pj_qsc_ns::AREA_1) {
            t = r;
            r = -s;
            s = t;
        } else if (area == pj_qsc_ns::AREA_2) {
            r = -r;
            s = -s;
        } else if (area == pj_qsc_ns::AREA_3) {
            t = r;
            r = s;
            s = -t;
        }
        /* Rotate q,r,s into the correct cube face. */
        if (Q->face == pj_qsc_ns::FACE_RIGHT) {
            t = q;
            q = -r;
            r = t;
        } else if (Q->face == pj_qsc_ns::FACE_BACK) {
            q = -q;
            r = -r;
        } else if (Q->face == pj_qsc_ns::FACE_LEFT) {
            t = q;
            q = r;
            r = -t;
        }
        /* Now compute phi and lam from the unit sphere coordinates. */
        lp.phi = acos(-s) - M_HALFPI;
        lp.lam = atan2(r, q);
        if (Q->face == pj_qsc_ns::FACE_RIGHT) {
            lp.lam = qsc_shift_longitude_origin(lp.lam, -M_HALFPI);
        } else if (Q->face == pj_qsc_ns::FACE_BACK) {
            lp.lam = qsc_shift_longitude_origin(lp.lam, -M_PI);
        } else if (Q->face == pj_qsc_ns::FACE_LEFT) {
            lp.lam = qsc_shift_longitude_origin(lp.lam, +M_HALFPI);
        }
    }

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

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

    P->inv = qsc_e_inverse;
    P->fwd = qsc_e_forward;
    /* Determine the cube face from the center of projection. */
    if (P->phi0 >= M_HALFPI - M_FORTPI / 2.0) {
        Q->face = pj_qsc_ns::FACE_TOP;
    } else if (P->phi0 <= -(M_HALFPI - M_FORTPI / 2.0)) {
        Q->face = pj_qsc_ns::FACE_BOTTOM;
    } else if (fabs(P->lam0) <= M_FORTPI) {
        Q->face = pj_qsc_ns::FACE_FRONT;
    } else if (fabs(P->lam0) <= M_HALFPI + M_FORTPI) {
        Q->face =
            (P->lam0 > 0.0 ? pj_qsc_ns::FACE_RIGHT : pj_qsc_ns::FACE_LEFT);
    } else {
        Q->face = pj_qsc_ns::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->b = P->a * sqrt(1.0 - P->es);
        Q->one_minus_f = 1.0 - (P->a - Q->b) / P->a;
        Q->one_minus_f_squared = Q->one_minus_f * Q->one_minus_f;
    }

    return P;
}

#undef EPS10

Coverage Report

Created: 2025-08-26 07:08

Line	Count	Source (jump to first uncovered line)
1		/*
2		* This implements the Quadrilateralized Spherical Cube (QSC) projection.
3		*
4		* Copyright (c) 2011, 2012 Martin Lambers <marlam@marlam.de>
5		*
6		* The QSC projection was introduced in:
7		* [OL76]
8		* E.M. O'Neill and R.E. Laubscher, "Extended Studies of a Quadrilateralized
9		* Spherical Cube Earth Data Base", Naval Environmental Prediction Research
10		* Facility Tech. Report NEPRF 3-76 (CSC), May 1976.
11		*
12		* The preceding shift from an ellipsoid to a sphere, which allows to apply
13		* this projection to ellipsoids as used in the Ellipsoidal Cube Map model,
14		* is described in
15		* [LK12]
16		* M. Lambers and A. Kolb, "Ellipsoidal Cube Maps for Accurate Rendering of
17		* Planetary-Scale Terrain Data", Proc. Pacific Graphics (Short Papers), Sep.
18		* 2012
19		*
20		* You have to choose one of the following projection centers,
21		* corresponding to the centers of the six cube faces:
22		* phi0 = 0.0, lam0 = 0.0 ("front" face)
23		* phi0 = 0.0, lam0 = 90.0 ("right" face)
24		* phi0 = 0.0, lam0 = 180.0 ("back" face)
25		* phi0 = 0.0, lam0 = -90.0 ("left" face)
26		* phi0 = 90.0 ("top" face)
27		* phi0 = -90.0 ("bottom" face)
28		* Other projection centers will not work!
29		*
30		* In the projection code below, each cube face is handled differently.
31		* See the computation of the face parameter in the PROJECTION(qsc) function
32		* and the handling of different face values (FACE_*) in the forward and
33		* inverse projections.
34		*
35		* Furthermore, the projection is originally only defined for theta angles
36		* between (-1/4 * PI) and (+1/4 * PI) on the current cube face. This area
37		* of definition is named pj_qsc_ns::AREA_0 in the projection code below. The
38		* other three areas of a cube face are handled by rotation of
39		* pj_qsc_ns::AREA_0.
40		*/
41
42		#include <errno.h>
43		#include <math.h>
44
45		#include "proj.h"
46		#include "proj_internal.h"
47
48		/* The six cube faces. */
49		namespace pj_qsc_ns {
50		enum Face {
51		FACE_FRONT = 0,
52		FACE_RIGHT = 1,
53		FACE_BACK = 2,
54		FACE_LEFT = 3,
55		FACE_TOP = 4,
56		FACE_BOTTOM = 5
57		};
58		}
59
60		namespace { // anonymous namespace
61		struct pj_qsc_data {
62		enum pj_qsc_ns::Face face;
63		double a_squared;
64		double b;
65		double one_minus_f;
66		double one_minus_f_squared;
67		};
68		} // anonymous namespace
69		PROJ_HEAD(qsc, "Quadrilateralized Spherical Cube") "\n\tAzi, Sph";
70
71	201k	#define EPS10 1.e-10
72
73		/* The four areas on a cube face. AREA_0 is the area of definition,
74		* the other three areas are counted counterclockwise. */
75		namespace pj_qsc_ns {
76		enum Area { AREA_0 = 0, AREA_1 = 1, AREA_2 = 2, AREA_3 = 3 };
77		}
78
79		/* Helper function for forward projection: compute the theta angle
80		* and determine the area number. */
81		static double qsc_fwd_equat_face_theta(double phi, double y, double x,
82	201k	enum pj_qsc_ns::Area *area) {
83	201k	double theta;
84	201k	if (phi < EPS10) {
85	0	*area = pj_qsc_ns::AREA_0;
86	0	theta = 0.0;
87	201k	} else {
88	201k	theta = atan2(y, x);
89	201k	if (fabs(theta) <= M_FORTPI) {
90	0	*area = pj_qsc_ns::AREA_0;
91	201k	} else if (theta > M_FORTPI && theta <= M_HALFPI + M_FORTPI) {
92	89.7k	*area = pj_qsc_ns::AREA_1;
93	89.7k	theta -= M_HALFPI;
94	111k	} else if (theta > M_HALFPI + M_FORTPI \|\|
95	111k	theta <= -(M_HALFPI + M_FORTPI)) {
96	111k	*area = pj_qsc_ns::AREA_2;
97	111k	theta = (theta >= 0.0 ? theta - M_PI : theta + M_PI);
98	111k	} else {
99	0	*area = pj_qsc_ns::AREA_3;
100	0	theta += M_HALFPI;
101	0	}
102	201k	}
103	201k	return theta;
104	201k	}
105
106		/* Helper function: shift the longitude. */
107	30.2k	static double qsc_shift_longitude_origin(double longitude, double offset) {
108	30.2k	double slon = longitude + offset;
109	30.2k	if (slon < -M_PI) {
110	0	slon += M_TWOPI;
111	30.2k	} else if (slon > +M_PI) {
112	2.74k	slon -= M_TWOPI;
113	2.74k	}
114	30.2k	return slon;
115	30.2k	}
116
117	202k	static PJ_XY qsc_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
118	202k	PJ_XY xy = {0.0, 0.0};
119	202k	struct pj_qsc_data Q = static_cast<struct pj_qsc_data >(P->opaque);
120	202k	double lat, longitude;
121	202k	double theta, phi;
122	202k	double t, mu; /* nu; */
123	202k	enum pj_qsc_ns::Area area;
124
125		/* Convert the geodetic latitude to a geocentric latitude.
126		* This corresponds to the shift from the ellipsoid to the sphere
127		* described in [LK12]. */
128	202k	if (P->es != 0.0) {
129	202k	lat = atan(Q->one_minus_f_squared * tan(lp.phi));
130	202k	} else {
131	0	lat = lp.phi;
132	0	}
133
134		/* Convert the input lat, longitude into theta, phi as used by QSC.
135		* This depends on the cube face and the area on it.
136		* For the top and bottom face, we can compute theta and phi
137		* directly from phi, lam. For the other faces, we must use
138		* unit sphere cartesian coordinates as an intermediate step. */
139	202k	longitude = lp.lam;
140	202k	if (Q->face == pj_qsc_ns::FACE_TOP) {
141	588	phi = M_HALFPI - lat;
142	588	if (longitude >= M_FORTPI && longitude <= M_HALFPI + M_FORTPI) {
143	0	area = pj_qsc_ns::AREA_0;
144	0	theta = longitude - M_HALFPI;
145	588	} else if (longitude > M_HALFPI + M_FORTPI \|\|
146	588	longitude <= -(M_HALFPI + M_FORTPI)) {
147	0	area = pj_qsc_ns::AREA_1;
148	0	theta = (longitude > 0.0 ? longitude - M_PI : longitude + M_PI);
149	588	} else if (longitude > -(M_HALFPI + M_FORTPI) &&
150	588	longitude <= -M_FORTPI) {
151	0	area = pj_qsc_ns::AREA_2;
152	0	theta = longitude + M_HALFPI;
153	588	} else {
154	588	area = pj_qsc_ns::AREA_3;
155	588	theta = longitude;
156	588	}
157	201k	} else if (Q->face == pj_qsc_ns::FACE_BOTTOM) {
158	0	phi = M_HALFPI + lat;
159	0	if (longitude >= M_FORTPI && longitude <= M_HALFPI + M_FORTPI) {
160	0	area = pj_qsc_ns::AREA_0;
161	0	theta = -longitude + M_HALFPI;
162	0	} else if (longitude < M_FORTPI && longitude >= -M_FORTPI) {
163	0	area = pj_qsc_ns::AREA_1;
164	0	theta = -longitude;
165	0	} else if (longitude < -M_FORTPI &&
166	0	longitude >= -(M_HALFPI + M_FORTPI)) {
167	0	area = pj_qsc_ns::AREA_2;
168	0	theta = -longitude - M_HALFPI;
169	0	} else {
170	0	area = pj_qsc_ns::AREA_3;
171	0	theta = (longitude > 0.0 ? -longitude + M_PI : -longitude - M_PI);
172	0	}
173	201k	} else {
174	201k	double q, r, s;
175	201k	double sinlat, coslat;
176	201k	double sinlon, coslon;
177
178	201k	if (Q->face == pj_qsc_ns::FACE_RIGHT) {
179	5.12k	longitude = qsc_shift_longitude_origin(longitude, +M_HALFPI);
180	196k	} else if (Q->face == pj_qsc_ns::FACE_BACK) {
181	25.1k	longitude = qsc_shift_longitude_origin(longitude, +M_PI);
182	171k	} else if (Q->face == pj_qsc_ns::FACE_LEFT) {
183	0	longitude = qsc_shift_longitude_origin(longitude, -M_HALFPI);
184	0	}
185	201k	sinlat = sin(lat);
186	201k	coslat = cos(lat);
187	201k	sinlon = sin(longitude);
188	201k	coslon = cos(longitude);
189	201k	q = coslat * coslon;
190	201k	r = coslat * sinlon;
191	201k	s = sinlat;
192
193	201k	if (Q->face == pj_qsc_ns::FACE_FRONT) {
194	171k	phi = acos(q);
195	171k	theta = qsc_fwd_equat_face_theta(phi, s, r, &area);
196	171k	} else if (Q->face == pj_qsc_ns::FACE_RIGHT) {
197	5.12k	phi = acos(r);
198	5.12k	theta = qsc_fwd_equat_face_theta(phi, s, -q, &area);
199	25.1k	} else if (Q->face == pj_qsc_ns::FACE_BACK) {
200	25.1k	phi = acos(-q);
201	25.1k	theta = qsc_fwd_equat_face_theta(phi, s, -r, &area);
202	25.1k	} else if (Q->face == pj_qsc_ns::FACE_LEFT) {
203	0	phi = acos(-r);
204	0	theta = qsc_fwd_equat_face_theta(phi, s, q, &area);
205	0	} else {
206		/* Impossible */
207	0	phi = theta = 0.0;
208	0	area = pj_qsc_ns::AREA_0;
209	0	}
210	201k	}
211
212		/* Compute mu and nu for the area of definition.
213		* For mu, see Eq. (3-21) in [OL76], but note the typos:
214		* compare with Eq. (3-14). For nu, see Eq. (3-38). */
215	202k	mu = atan((12.0 / M_PI) *
216	202k	(theta + acos(sin(theta) * cos(M_FORTPI)) - M_HALFPI));
217	202k	t = sqrt((1.0 - cos(phi)) / (cos(mu) * cos(mu)) /
218	202k	(1.0 - cos(atan(1.0 / cos(theta)))));
219		/* nu = atan(t); We don't really need nu, just t, see below. */
220
221		/* Apply the result to the real area. */
222	202k	if (area == pj_qsc_ns::AREA_1) {
223	89.7k	mu += M_HALFPI;
224	112k	} else if (area == pj_qsc_ns::AREA_2) {
225	111k	mu += M_PI;
226	111k	} else if (area == pj_qsc_ns::AREA_3) {
227	588	mu += M_PI_HALFPI;
228	588	}
229
230		/* Now compute x, y from mu and nu */
231		/* t = tan(nu); */
232	202k	xy.x = t * cos(mu);
233	202k	xy.y = t * sin(mu);
234	202k	return xy;
235	202k	}
236
237	0	static PJ_LP qsc_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
238	0	PJ_LP lp = {0.0, 0.0};
239	0	struct pj_qsc_data Q = static_cast<struct pj_qsc_data >(P->opaque);
240	0	double mu, nu, cosmu, tannu;
241	0	double tantheta, theta, cosphi, phi;
242	0	double t;
243	0	int area;
244
245		/* Convert the input x, y to the mu and nu angles as used by QSC.
246		* This depends on the area of the cube face. */
247	0	nu = atan(sqrt(xy.x * xy.x + xy.y * xy.y));
248	0	mu = atan2(xy.y, xy.x);
249	0	if (xy.x >= 0.0 && xy.x >= fabs(xy.y)) {
250	0	area = pj_qsc_ns::AREA_0;
251	0	} else if (xy.y >= 0.0 && xy.y >= fabs(xy.x)) {
252	0	area = pj_qsc_ns::AREA_1;
253	0	mu -= M_HALFPI;
254	0	} else if (xy.x < 0.0 && -xy.x >= fabs(xy.y)) {
255	0	area = pj_qsc_ns::AREA_2;
256	0	mu = (mu < 0.0 ? mu + M_PI : mu - M_PI);
257	0	} else {
258	0	area = pj_qsc_ns::AREA_3;
259	0	mu += M_HALFPI;
260	0	}
261
262		/* Compute phi and theta for the area of definition.
263		* The inverse projection is not described in the original paper, but some
264		* good hints can be found here (as of 2011-12-14):
265		* http://fits.gsfc.nasa.gov/fitsbits/saf.93/saf.9302
266		* (search for "Message-Id: <9302181759.AA25477 at fits.cv.nrao.edu>") */
267	0	t = (M_PI / 12.0) * tan(mu);
268	0	tantheta = sin(t) / (cos(t) - (1.0 / sqrt(2.0)));
269	0	theta = atan(tantheta);
270	0	cosmu = cos(mu);
271	0	tannu = tan(nu);
272	0	cosphi = 1.0 - cosmu * cosmu * tannu * tannu *
273	0	(1.0 - cos(atan(1.0 / cos(theta))));
274	0	if (cosphi < -1.0) {
275	0	cosphi = -1.0;
276	0	} else if (cosphi > +1.0) {
277	0	cosphi = +1.0;
278	0	}
279
280		/* Apply the result to the real area on the cube face.
281		* For the top and bottom face, we can compute phi and lam directly.
282		* For the other faces, we must use unit sphere cartesian coordinates
283		* as an intermediate step. */
284	0	if (Q->face == pj_qsc_ns::FACE_TOP) {
285	0	phi = acos(cosphi);
286	0	lp.phi = M_HALFPI - phi;
287	0	if (area == pj_qsc_ns::AREA_0) {
288	0	lp.lam = theta + M_HALFPI;
289	0	} else if (area == pj_qsc_ns::AREA_1) {
290	0	lp.lam = (theta < 0.0 ? theta + M_PI : theta - M_PI);
291	0	} else if (area == pj_qsc_ns::AREA_2) {
292	0	lp.lam = theta - M_HALFPI;
293	0	} else /* area == pj_qsc_ns::AREA_3 */ {
294	0	lp.lam = theta;
295	0	}
296	0	} else if (Q->face == pj_qsc_ns::FACE_BOTTOM) {
297	0	phi = acos(cosphi);
298	0	lp.phi = phi - M_HALFPI;
299	0	if (area == pj_qsc_ns::AREA_0) {
300	0	lp.lam = -theta + M_HALFPI;
301	0	} else if (area == pj_qsc_ns::AREA_1) {
302	0	lp.lam = -theta;
303	0	} else if (area == pj_qsc_ns::AREA_2) {
304	0	lp.lam = -theta - M_HALFPI;
305	0	} else /* area == pj_qsc_ns::AREA_3 */ {
306	0	lp.lam = (theta < 0.0 ? -theta - M_PI : -theta + M_PI);
307	0	}
308	0	} else {
309		/* Compute phi and lam via cartesian unit sphere coordinates. */
310	0	double q, r, s;
311	0	q = cosphi;
312	0	t = q * q;
313	0	if (t >= 1.0) {
314	0	s = 0.0;
315	0	} else {
316	0	s = sqrt(1.0 - t) * sin(theta);
317	0	}
318	0	t += s * s;
319	0	if (t >= 1.0) {
320	0	r = 0.0;
321	0	} else {
322	0	r = sqrt(1.0 - t);
323	0	}
324		/* Rotate q,r,s into the correct area. */
325	0	if (area == pj_qsc_ns::AREA_1) {
326	0	t = r;
327	0	r = -s;
328	0	s = t;
329	0	} else if (area == pj_qsc_ns::AREA_2) {
330	0	r = -r;
331	0	s = -s;
332	0	} else if (area == pj_qsc_ns::AREA_3) {
333	0	t = r;
334	0	r = s;
335	0	s = -t;
336	0	}
337		/* Rotate q,r,s into the correct cube face. */
338	0	if (Q->face == pj_qsc_ns::FACE_RIGHT) {
339	0	t = q;
340	0	q = -r;
341	0	r = t;
342	0	} else if (Q->face == pj_qsc_ns::FACE_BACK) {
343	0	q = -q;
344	0	r = -r;
345	0	} else if (Q->face == pj_qsc_ns::FACE_LEFT) {
346	0	t = q;
347	0	q = r;
348	0	r = -t;
349	0	}
350		/* Now compute phi and lam from the unit sphere coordinates. */
351	0	lp.phi = acos(-s) - M_HALFPI;
352	0	lp.lam = atan2(r, q);
353	0	if (Q->face == pj_qsc_ns::FACE_RIGHT) {
354	0	lp.lam = qsc_shift_longitude_origin(lp.lam, -M_HALFPI);
355	0	} else if (Q->face == pj_qsc_ns::FACE_BACK) {
356	0	lp.lam = qsc_shift_longitude_origin(lp.lam, -M_PI);
357	0	} else if (Q->face == pj_qsc_ns::FACE_LEFT) {
358	0	lp.lam = qsc_shift_longitude_origin(lp.lam, +M_HALFPI);
359	0	}
360	0	}
361
362		/* Apply the shift from the sphere to the ellipsoid as described
363		* in [LK12]. */
364	0	if (P->es != 0.0) {
365	0	int invert_sign;
366	0	double tanphi, xa;
367	0	invert_sign = (lp.phi < 0.0 ? 1 : 0);
368	0	tanphi = tan(lp.phi);
369	0	xa = Q->b / sqrt(tanphi * tanphi + Q->one_minus_f_squared);
370	0	lp.phi = atan(sqrt(P->a * P->a - xa * xa) / (Q->one_minus_f * xa));
371	0	if (invert_sign) {
372	0	lp.phi = -lp.phi;
373	0	}
374	0	}
375	0	return lp;
376	0	}
377
378	1.94k	PJ *PJ_PROJECTION(qsc) {
379	1.94k	struct pj_qsc_data Q = static_cast<struct pj_qsc_data >(
380	1.94k	calloc(1, sizeof(struct pj_qsc_data)));
381	1.94k	if (nullptr == Q)
382	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
383	1.94k	P->opaque = Q;
384
385	1.94k	P->inv = qsc_e_inverse;
386	1.94k	P->fwd = qsc_e_forward;
387		/* Determine the cube face from the center of projection. */
388	1.94k	if (P->phi0 >= M_HALFPI - M_FORTPI / 2.0) {
389	25	Q->face = pj_qsc_ns::FACE_TOP;
390	1.91k	} else if (P->phi0 <= -(M_HALFPI - M_FORTPI / 2.0)) {
391	36	Q->face = pj_qsc_ns::FACE_BOTTOM;
392	1.87k	} else if (fabs(P->lam0) <= M_FORTPI) {
393	1.79k	Q->face = pj_qsc_ns::FACE_FRONT;
394	1.79k	} else if (fabs(P->lam0) <= M_HALFPI + M_FORTPI) {
395	38	Q->face =
396	38	(P->lam0 > 0.0 ? pj_qsc_ns::FACE_RIGHT : pj_qsc_ns::FACE_LEFT);
397	46	} else {
398	46	Q->face = pj_qsc_ns::FACE_BACK;
399	46	}
400		/* Fill in useful values for the ellipsoid <-> sphere shift
401		* described in [LK12]. */
402	1.94k	if (P->es != 0.0) {
403	1.92k	Q->a_squared = P->a * P->a;
404	1.92k	Q->b = P->a * sqrt(1.0 - P->es);
405	1.92k	Q->one_minus_f = 1.0 - (P->a - Q->b) / P->a;
406	1.92k	Q->one_minus_f_squared = Q->one_minus_f * Q->one_minus_f;
407	1.92k	}
408
409	1.94k	return P;
410	1.94k	}
411
412		#undef EPS10