Add support for proper lunar mesh instead of the ellipse

.appveyor.yml

@@ -42,8 +42,10 @@ environment:
     VSPATH: C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\SDK\ScopeCppSDK\VC\bin
     QT_VERSION_MAJOR: 5
     QT_BASEDIR: C:\Qt\5.12\msvc2017_64
-    CMAKE_GENERATOR: Visual Studio 15 2017 Win64
-    CMAKE_ARGS: --
+    CMAKE_GENERATOR: Visual Studio 15 2017
+    CMAKE_ARGS: -A x64
+    cmakeURL: https://github.com/Kitware/CMake/releases/download/v3.18.6/cmake-3.18.6-win64-x64.zip
+    cmakeBaseName: cmake-3.18.6-win64-x64
     exiv2url: https://github.com/10110111/exiv2/releases/download/ver0.28.0-final/exiv2-0.28.0-2017msvc64.zip
     exiv2baseName: exiv2-0.28.0-2017msvc64
     scConverterEnabled: false
@@ -113,6 +115,10 @@ before_build:
 
   - ps: if($env:exiv2url -ne $null) { appveyor DownloadFile $env:exiv2url -FileName c:\$env:exiv2baseName.zip }
   - ps: if($env:exiv2url -ne $null) { 7z e c:\$env:exiv2baseName.zip -spf -oc:\ }
+  - ps: if($env:cmakeURL -ne $null) { appveyor DownloadFile $env:cmakeURL -FileName c:\$env:cmakeBaseName.zip }
+  - ps: if($env:cmakeURL -ne $null) { 7z e c:\$env:cmakeBaseName.zip -spf -oc:\ }
+  - ps: if($env:cmakeURL -ne $null) { $env:PATH="c:\$env:cmakeBaseName\bin;$env:PATH" }
+
   - if [%PUBLISH_BINARY%]==[true] appveyor DownloadFile https://github.com/Stellarium/stellarium-data/releases/download/guide/guide.pdf -FileName c:\stellarium\guide\guide.pdf
 
   - ps: if($env:INSTALL_CONVERTER_DEPS -eq "true") { appveyor DownloadFile $env:gettextURL -FileName c:\$env:gettextBaseName.zip }

CMakeLists.txt

@@ -1,4 +1,4 @@
-CMAKE_MINIMUM_REQUIRED(VERSION 3.16)
+CMAKE_MINIMUM_REQUIRED(VERSION 3.18)
 # CMake 3.20+ is required if x-compiling for Windows on ARM
 MESSAGE(STATUS "Found CMake ${CMAKE_VERSION}")
 

data/shaders/planet.frag

@@ -35,6 +35,7 @@ uniform mediump float skyBrightness;
 uniform bool hasAtmosphere;
 uniform bool hasNormalMap;
 uniform bool hasHorizonMap;
+uniform bool texCoordsFromFragment;
 
 uniform int shadowCount;
 uniform highp mat4 shadowData;
@@ -171,6 +172,15 @@ lowp float outgasFactor(in mediump vec3 normal, in highp vec3 lightDir, in mediu
     return opac;
 }
 
+vec4 sampleTexDxDy(sampler2D sampler, vec2 texCoord, vec2 texDx, vec2 texDy)
+{
+#ifdef textureGrad_SUPPORTED
+    if(texCoordsFromFragment)
+        return textureGrad(sampler, texCoord, texDx, texDy);
+#endif
+    return texture2D(sampler, texCoord);
+}
+
 void main()
 {
 #ifndef IS_MOON
@@ -288,9 +298,36 @@ void main()
     }
 
 #ifdef IS_MOON
+    mediump vec2 moonTexCoord = texc;
+    mediump vec2 texDx = vec2(0), texDy = vec2(0);
+    if(texCoordsFromFragment)
+    {
+        moonTexCoord = vec2(atan(normalZ.x, -normalZ.y)/(2.*PI)+0.5, asin(normalize(normalZ).z)/PI+0.5);
+
+#ifdef textureGrad_SUPPORTED
+        // The usual automatic computation of derivatives of texture coordinates
+        // breaks down at the discontinuity of atan, resulting in choosing the most
+        // minified mip level instead of the correct one, which looks as a seam on
+        // the screen. Thus, we need to compute them in a custom way, treating atan
+        // as a (continuous) multivalued function. We differentiate
+        // atan(normalZ.x(x,y), -normalZ.y(x,y)) with respect to x and y and yield
+        // gradLongitude vector.
+        vec2 gradNormalZX = vec2(dFdx(normalZ.x), dFdy(normalZ.x));
+        vec2 gradNormalZY = vec2(dFdx(-normalZ.y), dFdy(-normalZ.y));
+        vec2 gradLongitude = vec2(-normalZ.y*gradNormalZX.s-normalZ.x*gradNormalZY.s,
+                                  -normalZ.y*gradNormalZX.t-normalZ.x*gradNormalZY.t)
+                                                       /
+                                             dot(normalZ, normalZ);
+        float texTdx = dFdx(moonTexCoord.t);
+        float texTdy = dFdy(moonTexCoord.t);
+        texDx = vec2(gradLongitude.s/(2.*PI), texTdx);
+        texDy = vec2(gradLongitude.t/(2.*PI), texTdy);
+#endif
+    }
+
     mediump vec3 normal;
     if(hasNormalMap)
-        normal = texture2D(normalMap, texc).rgb-vec3(0.5);
+        normal = sampleTexDxDy(normalMap, moonTexCoord, texDx, texDy).rgb-vec3(0.5);
     else
         normal = vec3(0,0,1);
 
@@ -306,7 +343,7 @@ void main()
         mediump vec3 zenith = normalZ;
         mediump float sunAzimuth = atan(dot(lightDirection,lonDir), dot(lightDirection,northDir));
         mediump float sinSunElevation = dot(zenith, lightDirection);
-        mediump vec4 horizonElevSample = (texture2D(horizonMap, texc) - 0.5) * 2.;
+        mediump vec4 horizonElevSample = (sampleTexDxDy(horizonMap, moonTexCoord, texDx, texDy) - 0.5) * 2.;
         mediump vec4 sinHorizElevs = sign(horizonElevSample) * horizonElevSample * horizonElevSample;
         mediump float sinHorizElevLeft, sinHorizElevRight;
         mediump float alpha;
@@ -390,7 +427,10 @@ void main()
     //litColor.xyz = clamp( litColor.xyz + vec3(outgas), 0.0, 1.0);
 
     lowp vec4 texColor;
-#ifdef RINGS_SUPPORT
+#ifdef IS_MOON
+    texColor = sampleTexDxDy(tex, moonTexCoord, texDx, texDy);
+#else
+# ifdef RINGS_SUPPORT
     if(isRing)
     {
         float radius = length(texc);
@@ -402,10 +442,11 @@ void main()
             texColor = vec4(0);
     }
     else
-#endif
+# endif // RINGS_SUPPORT
     {
         texColor = texture2D(tex, texc);
     }
+#endif // IS_MOON
 
     texColor.rgb = srgbToLinear(texColor.rgb);
 

data/shaders/planet.vert

@@ -23,7 +23,9 @@
 
 ATTRIBUTE highp vec4 vertex; // vertex projected by CPU
 ATTRIBUTE mediump vec2 texCoord;
+#ifndef PROJECTOR_PRESENT
 ATTRIBUTE highp vec4 unprojectedVertex; //original vertex coordinate (in km for OBJ models, in AU otherwise)
+#endif
 #ifdef IS_OBJ
     ATTRIBUTE mediump vec3 normalIn; // OBJs have pre-calculated normals
 #endif
@@ -49,9 +51,18 @@ VARYING highp vec3 P; //original unprojected position (in AU)
     uniform highp vec3 lightDirection;
 #endif
 
+#ifdef PROJECTOR_PRESENT
+uniform float sphereScale;
+#endif
+
 void main()
 {
+#ifdef PROJECTOR_PRESENT
+	highp vec4 unprojectedVertex = vertex;
+	gl_Position = projectionMatrix*vec4(project(sphereScale*vertex.xyz), 1);
+#else
     gl_Position = projectionMatrix * vertex;
+#endif
     texc = texCoord;
 
 #ifdef SHADOWMAP

models/CMakeLists.txt

@@ -1,3 +1,4 @@
+file(ARCHIVE_EXTRACT INPUT ${CMAKE_CURRENT_SOURCE_DIR}/moon-vertices-indices.bin.7z DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 
 ########### install files ###############
 
@@ -9,6 +10,7 @@ INSTALL(DIRECTORY ./ DESTINATION ${SDATALOC}/models
 	PATTERN "*.obj"
 	PATTERN "*.gz"
 	PATTERN "CMakeFiles" EXCLUDE )
+INSTALL(FILES ${CMAKE_CURRENT_BINARY_DIR}/moon-vertices-indices.bin DESTINATION ${SDATALOC}/models)
 
 INSTALL(DIRECTORY blender/ DESTINATION ${SDATALOC}/models/blender
 	FILES_MATCHING 

src/core/RefractionExtinction.cpp

@@ -282,13 +282,30 @@ QByteArray Refraction::getForwardTransformShader() const
 {
 	return QByteArray(1+R"(
 uniform float REFRACTION_press_temp_corr;
-vec3 innerRefractionForward(vec3 altAzPos)
+// Workaround for Intel's unusable implementation of sin & cos, which leads to broken geometry of the Moon every 0.9° of elevation.
+highp float REFRACTION_sin(highp float x)
 {
 	const float PI = 3.14159265;
-	const float M_180_PI = 180./PI;
-	const float M_PI_180 = PI/180.;
-	const float MIN_GEO_ALTITUDE_DEG=@MIN_GEO_ALTITUDE_DEG@;
-	const float TRANSITION_WIDTH_GEO_DEG=@TRANSITION_WIDTH_GEO_DEG@;
+	x = mod(x+PI, 2.0*PI)-PI;
+	return x*(0.999999599920672 + x*x*(-0.166665526354071 + x*x*(0.00833240298869917 + x*x*(-0.0001980863334175 + x*x*(2.69971463693744e-6 - 2.03622449118901e-8*x*x)))));
+}
+// Workaround for Intel's and AMD's unusable implementation of asin, which leads to time-dependent shifts of the Moon from its refracted positions.
+highp float REFRACTION_asin(highp float x)
+{
+	float sign = x < 0. ? -1. : 1.;
+	if(x < 0.) x = -x;
+	x = 2. * sqrt(1. - x) - 1.;
+	float v = 0.848061881596496 + x*(-0.755929497461161 + x*(-0.0539853235799928 + x*(-0.025701166121295 + x*(-0.00723634508887381 +
+	           x*(-0.00314276748524976 + x*(-0.00101365621070969 + x*(-0.000411380592372285 + x*(-0.000428167561924693 - 0.000213293541786718*x))))))));
+	return v * sign;
+}
+highp vec3 innerRefractionForward(highp vec3 altAzPos)
+{
+	const highp float PI = 3.14159265;
+	const highp float M_180_PI = 180./PI;
+	const highp float M_PI_180 = PI/180.;
+	const highp float MIN_GEO_ALTITUDE_DEG=@MIN_GEO_ALTITUDE_DEG@;
+	const highp float TRANSITION_WIDTH_GEO_DEG=@TRANSITION_WIDTH_GEO_DEG@;
 
 	float press_temp_corr = REFRACTION_press_temp_corr;
 
@@ -299,32 +316,32 @@ vec3 innerRefractionForward(vec3 altAzPos)
 		return altAzPos;
 	}
 
-	float sinGeo = altAzPos[2]/len;
-	float geom_alt_rad = asin(sinGeo);
-	float geom_alt_deg = M_180_PI*geom_alt_rad;
+	highp float sinGeo = altAzPos[2]/len;
+	highp float geom_alt_rad = REFRACTION_asin(sinGeo);
+	highp float geom_alt_deg = M_180_PI*geom_alt_rad;
 	if (geom_alt_deg > MIN_GEO_ALTITUDE_DEG)
 	{
 		// refraction from Saemundsson, S&T1986 p70 / in Meeus, Astr.Alg.
-		float r=press_temp_corr * ( 1.02 / tan((geom_alt_deg+10.3/(geom_alt_deg+5.11))*M_PI_180) + 0.0019279);
+		highp float r=press_temp_corr * ( 1.02 / tan((geom_alt_deg+10.3/(geom_alt_deg+5.11))*M_PI_180) + 0.0019279);
 		geom_alt_deg += r;
 		if (geom_alt_deg > 90.)
 			geom_alt_deg=90.;
 	}
 	else if(geom_alt_deg>MIN_GEO_ALTITUDE_DEG-TRANSITION_WIDTH_GEO_DEG)
 	{
 		// Avoids the jump below -5 by interpolating linearly between MIN_GEO_ALTITUDE_DEG and bottom of transition zone
-		float r_m5=press_temp_corr * ( 1.02 / tan((MIN_GEO_ALTITUDE_DEG+10.3/(MIN_GEO_ALTITUDE_DEG+5.11))*M_PI_180) + 0.0019279);
+		highp float r_m5=press_temp_corr * ( 1.02 / tan((MIN_GEO_ALTITUDE_DEG+10.3/(MIN_GEO_ALTITUDE_DEG+5.11))*M_PI_180) + 0.0019279);
 		geom_alt_deg += r_m5*(geom_alt_deg-(MIN_GEO_ALTITUDE_DEG-TRANSITION_WIDTH_GEO_DEG))/TRANSITION_WIDTH_GEO_DEG;
 	}
 	else return altAzPos;
 	// At this point we have corrected geometric altitude. Note that if we just change altAzPos[2], we would change vector length, so this would change our angles.
 	// We have to shorten X,Y components of the vector as well by the change in cosines of altitude, or (sqrt(1-sin(alt))
 
-	float refr_alt_rad=geom_alt_deg*M_PI_180;
-	float sinRef=sin(refr_alt_rad);
+	highp float refr_alt_rad=geom_alt_deg*M_PI_180;
+	highp float sinRef = REFRACTION_sin(refr_alt_rad);
 
 	// FIXME: do we really need double's mantissa length here as a comment in the C++ code says?
-	float shortenxy = abs(sinGeo)>=1.0 ? 1.0 : sqrt((1.-sinRef*sinRef)/(1.-sinGeo*sinGeo));
+	highp float shortenxy = abs(sinGeo)>=1.0 ? 1.0 : sqrt((1.-sinRef*sinRef)/(1.-sinGeo*sinGeo));
 
 	altAzPos[0]*=shortenxy;
 	altAzPos[1]*=shortenxy;
@@ -333,18 +350,18 @@ vec3 innerRefractionForward(vec3 altAzPos)
 	return altAzPos;
 }
 
-uniform mat4 REFRACTION_preTransfoMat;
-uniform mat4 REFRACTION_postTransfoMat;
+uniform highp mat4 REFRACTION_preTransfoMat;
+uniform highp mat4 REFRACTION_postTransfoMat;
 
-vec3 vertexToAltAzPos(vec3 v)
+highp vec3 vertexToAltAzPos(highp vec3 v)
 {
-	vec3 altAzPosNotRefracted = (REFRACTION_preTransfoMat * vec4(v,1)).xyz;
+	highp vec3 altAzPosNotRefracted = (REFRACTION_preTransfoMat * vec4(v,1)).xyz;
 	return innerRefractionForward(altAzPosNotRefracted);
 }
 
-vec3 modelViewForwardTransform(vec3 v)
+highp vec3 modelViewForwardTransform(highp vec3 v)
 {
-	vec3 altAzPos = vertexToAltAzPos(v);
+	highp vec3 altAzPos = vertexToAltAzPos(v);
 	return (REFRACTION_postTransfoMat * vec4(altAzPos, 1)).xyz;
 }
 

src/core/StelProjector.cpp

@@ -92,11 +92,11 @@ QByteArray StelProjector::Mat4dTransform::getForwardTransformShader() const
 #line 1 104
 uniform mat4 PROJECTOR_vertexToAltAzPosMatrix;
 uniform mat4 PROJECTOR_modelViewMatrix;
-vec3 modelViewForwardTransform(vec3 v)
+vec3 modelViewForwardTransform(highp vec3 v)
 {
 	return (PROJECTOR_modelViewMatrix * vec4(v,1)).xyz;
 }
-vec3 vertexToAltAzPos(vec3 v)
+vec3 vertexToAltAzPos(highp vec3 v)
 {
 	return (PROJECTOR_vertexToAltAzPosMatrix * vec4(v,1)).xyz;
 }
@@ -466,21 +466,21 @@ QByteArray StelProjector::getProjectShader() const
 {
 	static const char*const projectSrc = 1+R"(
 #line 1 106
-uniform vec2 PROJECTOR_viewportCenter;
-uniform vec2 PROJECTOR_flip;
-uniform float PROJECTOR_pixelPerRad;
-uniform float PROJECTOR_zNear;
-uniform float PROJECTOR_oneOverZNearMinusZFar;
-vec3 project(vec3 modelSpacePoint)
-{
-	float flipHorz = PROJECTOR_flip.x, flipVert = PROJECTOR_flip.y;
-	vec2  viewportCenter = PROJECTOR_viewportCenter;
-	float pixelPerRad = PROJECTOR_pixelPerRad;
-	float zNear = PROJECTOR_zNear;
-	float oneOverZNearMinusZFar = PROJECTOR_oneOverZNearMinusZFar;
-
-	vec3 worldSpacePoint = modelViewForwardTransform(modelSpacePoint);
-	vec3 v = projectorForwardTransform(worldSpacePoint);
+uniform highp vec2 PROJECTOR_viewportCenter;
+uniform highp vec2 PROJECTOR_flip;
+uniform highp float PROJECTOR_pixelPerRad;
+uniform highp float PROJECTOR_zNear;
+uniform highp float PROJECTOR_oneOverZNearMinusZFar;
+highp vec3 project(highp vec3 modelSpacePoint)
+{
+	highp float flipHorz = PROJECTOR_flip.x, flipVert = PROJECTOR_flip.y;
+	highp vec2  viewportCenter = PROJECTOR_viewportCenter;
+	highp float pixelPerRad = PROJECTOR_pixelPerRad;
+	highp float zNear = PROJECTOR_zNear;
+	highp float oneOverZNearMinusZFar = PROJECTOR_oneOverZNearMinusZFar;
+
+	highp vec3 worldSpacePoint = modelViewForwardTransform(modelSpacePoint);
+	highp vec3 v = projectorForwardTransform(worldSpacePoint);
 	return vec3(viewportCenter[0] + flipHorz * pixelPerRad * v[0],
 	            viewportCenter[1] + flipVert * pixelPerRad * v[1],
 	            (v[2] - zNear) * oneOverZNearMinusZFar);
@@ -504,25 +504,25 @@ QByteArray StelProjector::getUnProjectShader() const
 {
 	static const char*const projectSrc = 1+R"(
 #line 1 107
-uniform vec2 PROJECTOR_viewportCenter;
-uniform vec2 PROJECTOR_flip;
-uniform float PROJECTOR_pixelPerRad;
-vec3 winPosToWorldPos(float x, float y, out bool ok)
+uniform highp vec2 PROJECTOR_viewportCenter;
+uniform highp vec2 PROJECTOR_flip;
+uniform highp float PROJECTOR_pixelPerRad;
+highp vec3 winPosToWorldPos(highp float x, highp float y, out bool ok)
 {
-	float flipHorz = PROJECTOR_flip.x, flipVert = PROJECTOR_flip.y;
-	vec2  viewportCenter = PROJECTOR_viewportCenter;
-	float pixelPerRad = PROJECTOR_pixelPerRad;
+	highp float flipHorz = PROJECTOR_flip.x, flipVert = PROJECTOR_flip.y;
+	highp vec2  viewportCenter = PROJECTOR_viewportCenter;
+	highp float pixelPerRad = PROJECTOR_pixelPerRad;
 
-	vec3 v;
+	highp vec3 v;
 	v[0] = flipHorz * (x - viewportCenter[0]) / pixelPerRad;
 	v[1] = flipVert * (y - viewportCenter[1]) / pixelPerRad;
 	v[2] = 0.;
 	return projectorBackwardTransform(v, ok);
 }
 
-vec3 unProject(float x, float y, out bool ok)
+highp vec3 unProject(highp float x, highp float y, out bool ok)
 {
-	vec3 worldPos = winPosToWorldPos(x, y, ok);
+	highp vec3 worldPos = winPosToWorldPos(x, y, ok);
 	// Even when the reprojected point comes from a region of the screen,
 	// where nothing is projected to (rval=false), we finish reprojecting.
 	// This looks good for atmosphere rendering, and it helps avoiding