+namespace {
+/// map p onto cube, s gives the surface, u and v the position in [-1,1]
+void cubemap(const glm::dvec3 &p, int &s, double &u, double &v) noexcept {
+ const glm::dvec3 p_abs(abs(p));
+ const glm::bvec3 p_pos(greaterThan(p, glm::dvec3(0.0)));
+ double max_axis = 0.0;
+
+ if (p_pos.x && p_abs.x >= p_abs.y && p_abs.x >= p_abs.z) {
+ max_axis = p_abs.x;
+ u = p.y;
+ v = p.z;
+ s = 1;
+ }
+ if (!p_pos.x && p_abs.x >= p_abs.y && p_abs.x >= p_abs.z) {
+ max_axis = p_abs.x;
+ u = -p.y;
+ v = -p.z;
+ s = 4;
+ }
+ if (p_pos.y && p_abs.y >= p_abs.x && p_abs.y >= p_abs.z) {
+ max_axis = p_abs.y;
+ u = p.z;
+ v = p.x;
+ s = 2;
+ }
+ if (!p_pos.y && p_abs.y >= p_abs.x && p_abs.y >= p_abs.z) {
+ max_axis = p_abs.y;
+ u = -p.z;
+ v = -p.x;
+ s = 5;
+ }
+ if (p_pos.z && p_abs.z >= p_abs.x && p_abs.z >= p_abs.y) {
+ max_axis = p_abs.z;
+ u = p.x;
+ v = p.y;
+ s = 0;
+ }
+ if (!p_pos.z && p_abs.z >= p_abs.x && p_abs.z >= p_abs.y) {
+ max_axis = p_abs.z;
+ u = -p.x;
+ v = -p.y;
+ s = 3;
+ }
+ u /= max_axis;
+ v /= max_axis;
+}
+/// get p from cube, s being surface, u and v the position in [-1,1],
+/// gives a vector from the center to the surface
+glm::dvec3 cubeunmap(int s, double u, double v) {
+ switch (s) {
+ default:
+ case 0: return glm::dvec3(u, v, 1.0); // +Z
+ case 1: return glm::dvec3(1.0, u, v); // +X
+ case 2: return glm::dvec3(v, 1.0, u); // +Y
+ case 3: return glm::dvec3(-u, -v, -1.0); // -Z
+ case 4: return glm::dvec3(-1.0, -u, -v); // -X
+ case 5: return glm::dvec3(-v, -1.0, -u); // -Y
+ };
+}