From 908570c623d074078a75bd63093fe1002bf4ee1d Mon Sep 17 00:00:00 2001 From: dusan-prascevic Date: Fri, 21 Feb 2025 16:36:16 +0100 Subject: [PATCH] working implementation for forward/backwards transform + shaders, however some issues with display outside of the projection ellipse --- src/core/StelCore.cpp | 3 + src/core/StelProjectorClasses.cpp | 154 +++++++++++++++++------------- src/core/StelProjectorClasses.hpp | 19 +++- 3 files changed, 104 insertions(+), 72 deletions(-) diff --git a/src/core/StelCore.cpp b/src/core/StelCore.cpp index 258bba968d51b..b894b0bd29b74 100644 --- a/src/core/StelCore.cpp +++ b/src/core/StelCore.cpp @@ -432,6 +432,9 @@ StelProjectorP StelCore::getProjection(StelProjector::ModelViewTranformP modelVi case ProjectionHammer: prj = StelProjectorP(new StelProjectorHammer(modelViewTransform)); break; + case ProjectionMollweide: + prj = StelProjectorP(new StelProjectorMollweide(modelViewTransform)); + break; case ProjectionCylinder: prj = StelProjectorP(new StelProjectorCylinder(modelViewTransform)); break; diff --git a/src/core/StelProjectorClasses.cpp b/src/core/StelProjectorClasses.cpp index 7558a21d76624..690ca8c0aa16f 100644 --- a/src/core/StelProjectorClasses.cpp +++ b/src/core/StelProjectorClasses.cpp @@ -552,76 +552,86 @@ QString StelProjectorMollweide::getDescriptionI18() const { return q_("The Mollweide projection is an equal-area map projection introduced by Karl Mollweide in 1805."); } -bool StelProjectorMollweide::forward(Vec3f &v) const { - // Convert from 3D sphere to (x,y) in projection space. +bool StelProjectorMollweide::forward(Vec3f &v) const +{ const float r = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); - // Compute latitude (φ) and longitude (λ) - const float phi = std::asin(v[1] / r); - const float lambda = std::atan2(v[0], -v[2]); - // Solve for theta in: 2θ + sin(2θ) = π sinφ - float theta = phi; // initial guess - const float piSinPhi = static_cast(M_PI) * std::sin(phi); - const int maxIter = 10; - const float tol = 1e-7f; - for (int i = 0; i < maxIter; ++i) { - float f = 2.f*theta + std::sin(2.f*theta) - piSinPhi; - float fprime = 2.f + 2.f*std::cos(2.f*theta); - float dtheta = f / fprime; - theta -= dtheta; - if (std::fabs(dtheta) < tol) - break; + const float lambda = std::atan2(v[0], -v[2]); // Longitude + const float sinDelta = v[1] / r; // Latitude (delta) sine + const float delta = std::asin(sinDelta); + + // Solve 2θ + sin(2θ) = π sinδ using Newton-Raphson + float theta = delta; + for (int i=0; i<5; ++i) { + float f = 2.f*theta + sin(2.f*theta) - M_PIf * sinDelta; + float df = 2.f + 2.f*cos(2.f*theta); + theta -= f / df; } - // Compute projection coordinates: - // widthStretch is a member (as in the Hammer code) - float x = (2.f * static_cast(M_SQRT2) / static_cast(M_PI)) * lambda * std::cos(theta) * widthStretch; - float y = static_cast(M_SQRT2) * std::sin(theta); - v[0] = x; - v[1] = y; - v[2] = r; // store original radius if needed + + const float cosTheta = cos(theta); + v[0] = (2.0f * M_SQRT2 / M_PIf) * lambda * cosTheta * widthStretch; + v[1] = M_SQRT2 * sin(theta); + v[2] = r; return true; } + bool StelProjectorMollweide::backward(Vec3d &v) const { - // Inverse: from projection (x,y) back to sphere. - v[0] /= widthStretch; - // Recover theta from y = √2 sinθ - double theta = std::asin(v[1] / static_cast(M_SQRT2)); - double cosTheta = std::cos(theta); - // Avoid division by zero near the poles - if (std::fabs(cosTheta) < 1e-10) + v[0] /= widthStretch; + const double x = v[0]; + const double y = v[1]; + + // Check if point lies within Mollweide ellipse bounds + const bool ok = (x*x + 4.0*y*y <= 8.0); + + // Compute theta (exit early if invalid) + if (std::abs(y) > M_SQRT2) return false; - // Recover λ from x = (2√2/π) λ cosθ - double lambda = (v[0] * static_cast(M_PI)) / (2.0 * static_cast(M_SQRT2) * cosTheta); - // Recover φ from: 2θ + sin(2θ) = π sinφ => φ = arcsin((2θ + sin2θ)/π) - double phi = std::asin((2.0*theta + std::sin(2.0*theta)) / static_cast(M_PI)); - double cosPhi = std::cos(phi); - // Convert spherical (φ,λ) to 3D vector (using same convention as in the Hammer code) - v[0] = cosPhi * std::sin(lambda); - v[1] = std::sin(phi); - v[2] = - cosPhi * std::cos(lambda); - return true; + const double theta = asin(y / M_SQRT2); + + // Reconstruct spherical coordinates (lambda, delta) + const double cosTheta = std::cos(theta); + const double lambda = (x * M_PI) / (2.0 * M_SQRT2 * cosTheta); + const double sinDelta = (2.0 * theta + sin(2.0*theta)) / M_PI; + const double delta = std::asin(sinDelta); + + // Convert to 3D Cartesian + const double cosDelta = std::cos(delta); + v[0] = cosDelta * std::sin(lambda); + v[1] = sinDelta; + v[2] = -cosDelta * std::cos(lambda); + return ok; } QByteArray StelProjectorMollweide::getForwardTransformShader() const { - return modelViewTransform->getForwardTransformShader() + R"( + return modelViewTransform->getForwardTransformShader() + R"( #line 1 102 // Forward transform shader snippet: uniform float PROJECTOR_FWD_widthStretch; -vec3 projectorForwardTransform(vec3 v) { +vec3 projectorForwardTransform(vec3 v) +{ + const float M_SQRT2 = 1.41421356; + const float PI = 3.14159265; + float widthStretch = PROJECTOR_FWD_widthStretch; + float r = length(v); - float phi = asin(v.y / r); - float lambda = atan(v.x, -v.z); - float theta = phi; - // Iterative solution for theta (10 iterations) - for (int i = 0; i < 10; i++) { - float f = 2.0 * theta + sin(2.0 * theta) - 3.14159265 * sin(phi); - float fprime = 2.0 + 2.0 * cos(2.0 * theta); - theta -= f / fprime; + float lambda = atan(v[0], -v[2]); + float sinDelta = v[1] / r; + float delta = asin(sinDelta); + + // Newton-Raphson iterations for theta + float theta = delta; + for (int i=0; i<5; i++) { + float f = 2.0*theta + sin(2.0*theta) - PI * sinDelta; + float df = 2.0 + 2.0*cos(2.0*theta); + theta = theta - f/df; } - float x = (2.82842712 / 3.14159265) * lambda * cos(theta) * PROJECTOR_FWD_widthStretch; // 2√2 = 2.82842712 - float y = 1.41421356 * sin(theta); // √2 = 1.41421356 - return vec3(x, y, r); + + float cosTheta = cos(theta); + v[0] = (2.0 * M_SQRT2 / PI) * lambda * cosTheta * widthStretch; + v[1] = M_SQRT2 * sin(theta); + v[2] = r; + return v; } #line 1 0 )"; @@ -633,22 +643,28 @@ QByteArray StelProjectorMollweide::getBackwardTransformShader() const #line 1 103 // Backward transform shader snippet: uniform float PROJECTOR_FWD_widthStretch; -vec3 projectorBackwardTransform(vec3 v, out bool ok) { - v.x /= PROJECTOR_FWD_widthStretch; - float theta = asin(v.y / 1.41421356); +vec3 projectorBackwardTransform(vec3 v, out bool ok) +{ + float widthStretch = PROJECTOR_FWD_widthStretch; + v[0] /= widthStretch; + float x = v[0]; + float y = v[1]; + + // Validity check: lies within Mollweide ellipse (x^2 +4y^2 <=8) + ok = (x*x + 4.0*y*y <= 8.0); + + float theta = asin(y / 1.41421356); // y ~ M_SQRT2*sin(theta) float cosTheta = cos(theta); - if (abs(cosTheta) < 0.00001) { - ok = false; - return v; - } - float lambda = (v.x * 3.14159265) / (2.82842712 * cosTheta); - float phi = asin((2.0 * theta + sin(2.0 * theta)) / 3.14159265); - float cosPhi = cos(phi); - float x = cosPhi * sin(lambda); - float y = sin(phi); - float z = - cosPhi * cos(lambda); - ok = true; - return vec3(x, y, z); + float lambda = (x * 3.14159265) / (2.0 * 1.41421356 * cosTheta); + + float sinDelta = (2.0*theta + sin(2.0*theta)) / 3.14159265; + float delta = asin(sinDelta); + + float cosDelta = cos(delta); + v[0] = cosDelta * sin(lambda); + v[1] = sinDelta; + v[2] = -cosDelta * cos(lambda); + return v; } #line 1 0 )"; diff --git a/src/core/StelProjectorClasses.hpp b/src/core/StelProjectorClasses.hpp index aab0244ad8601..057dfc665d3ae 100644 --- a/src/core/StelProjectorClasses.hpp +++ b/src/core/StelProjectorClasses.hpp @@ -113,19 +113,32 @@ class StelProjectorHammer : public StelProjector } }; -class StelProjectorMollweide : public StelProjector -{ +class StelProjectorMollweide : public StelProjector { public: StelProjectorMollweide(ModelViewTranformP func) : StelProjector(func) {} QString getNameI18() const override; QString getDescriptionI18() const override; - float getMaxFov() const override { return 180.f; } + float getMaxFov() const override { return 185.f; } bool forward(Vec3f &v) const override; bool backward(Vec3d &v) const override; + // float fovToViewScalingFactor(float fov) const override; + //float viewScalingFactorToFov(float vsf) const override; + //float deltaZoom(float fov) const override; QByteArray getForwardTransformShader() const override; QByteArray getBackwardTransformShader() const override; protected: +// TODO: take a look at the discontinuities => might be the reason for weird display glitches at high FOV bool hasDiscontinuity() const override { return true; } + bool intersectViewportDiscontinuityInternal(const Vec3d& p1, const Vec3d& p2) const override { + return p1[0] * p2[0] < 0 && !(p1[2] < 0 && p2[2] < 0); + } + bool intersectViewportDiscontinuityInternal(const Vec3d& capN, double capD) const override { + static const SphericalCap cap1(1, 0, 0); + static const SphericalCap cap2(-1, 0, 0); + static const SphericalCap cap3(0, 0, -1); + SphericalCap cap(capN, capD); + return cap.intersects(cap1) && cap.intersects(cap2) && cap.intersects(cap3); + } }; class StelProjectorCylinder : public StelProjector