Fix normals map of the Moon, increase its resolution and implement sh…

…adows

With the old normal map, near the poles we had dark-shaded parts of
craters that should actually be brighter — as if the Sun is not where it
is, and this behavior inconsistently depended on longitude... The new
texture added in this commit fixes behavior of shadows.

Now the resolution of the normal map is the half that of the reflectance
texture, 2048×1024. Further increase leads to excessive noise due to
lack of proper mip mapping (and its non-triviality for normal maps), so
this looks like a good compromise between the shapes of the craters and
background noise.

Additionally, a new kind of texture is added: horizon map. This texture
contains elevation of the horizon in four directions: North, East, South, West.
These elevations take terrain (mountains, craters etc.) surrounding each
point of the map into account.
These data are used in the fragment shader to determine whether the Sun is
visible or is hidden behind the terrain.

The data for the normals and horizon maps were generated using the
elevation data in the CGI Moon Kit:
https://svs.gsfc.nasa.gov/cgi-bin/details.cgi?aid=4720 .

data/shaders/planet.frag

@@ -24,6 +24,8 @@
 VARYING mediump vec2 texc; //texture coord
 VARYING highp vec3 P; //original vertex pos in model space
 
+const highp float PI = 3.14159265;
+
 uniform sampler2D tex;
 uniform mediump vec2 poleLat; //latitudes of pole caps, in terms of texture coordinate. x>0...north, y<1...south. 
 uniform mediump vec3 ambientLight;
@@ -63,6 +65,7 @@ VARYING highp vec4 shadowCoord;
     uniform sampler2D earthShadow;
     uniform mediump float eclipsePush;
     uniform sampler2D normalMap;
+	uniform sampler2D horizonMap;
 
     VARYING highp vec3 normalX;
     VARYING highp vec3 normalY;
@@ -263,6 +266,53 @@ void main()
     mediump vec3 normal = texture2D(normalMap, texc).rgb-vec3(0.5, 0.5, 0);
     normal = normalize(normalX*normal.x+normalY*normal.y+normalZ*normal.z);
     // normal now contains the real surface normal taking normal map into account
+
+	mediump float horizonShadowCoefficient = 1;
+	{
+		// Check whether the fragment is in the shadow of surrounding mountains or the horizon
+		mediump vec3 lonDir = normalX;
+		mediump vec3 northDir = normalY;
+		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 sinHorizElevs = sign(horizonElevSample) * horizonElevSample * horizonElevSample;
+		mediump float sinHorizElevLeft, sinHorizElevRight;
+		mediump float alpha;
+		if(sunAzimuth >= PI/2)
+		{
+			// Sun is between East and South
+			sinHorizElevLeft = sinHorizElevs[1];
+			sinHorizElevRight = sinHorizElevs[2];
+			alpha = (sunAzimuth - PI/2) / (PI/2);
+		}
+		else if(sunAzimuth >= 0)
+		{
+			// Sun is between North and East
+			sinHorizElevLeft = sinHorizElevs[0];
+			sinHorizElevRight = sinHorizElevs[1];
+			alpha = sunAzimuth / (PI/2);
+		}
+		else if(sunAzimuth <= -PI/2)
+		{
+			// Sun is between South and West
+			sinHorizElevLeft = sinHorizElevs[2];
+			sinHorizElevRight = sinHorizElevs[3];
+			alpha = (sunAzimuth + PI) / (PI/2);
+		}
+		else
+		{
+			// Sun is between West and North
+			sinHorizElevLeft = sinHorizElevs[3];
+			sinHorizElevRight = sinHorizElevs[0];
+			alpha = (sunAzimuth + PI/2) / (PI/2);
+		}
+		mediump float horizElevLeft = asin(sinHorizElevLeft);
+		mediump float horizElevRight = asin(sinHorizElevRight);
+		mediump float horizElev = horizElevLeft + (horizElevRight-horizElevLeft)*alpha;
+		if(sinSunElevation < sin(horizElev))
+			horizonShadowCoefficient = 0;
+	}
 #else
     // important to normalize here again
     mediump vec3 normal = normalize(normalVS);
@@ -271,6 +321,9 @@ void main()
     // Use an Oren-Nayar model for rough surfaces
     // Ref: http://content.gpwiki.org/index.php/D3DBook:(Lighting)_Oren-Nayar
     lum = orenNayar(normal, lightDirection, eyeDirection, orenNayarParameters.x, orenNayarParameters.y, orenNayarParameters.z, orenNayarParameters.w);
+#endif
+#ifdef IS_MOON
+	lum *= horizonShadowCoefficient;
 #endif
     //calculate pseudo-outgassing/rim-lighting effect
     lowp float outgas = 0.0;

data/ssystem_major.ini

@@ -1089,6 +1089,7 @@ coord_func=lunar_special
 landscape=moon
 name=Moon
 normals_map=moon_normals.png
+horizon_map=moon_horizon.png
 orbit_visualization_period=27.321582
 parent=Earth
 radius=1737.4

src/core/StelObserver.cpp

@@ -52,7 +52,7 @@ class ArtificialPlanet : public Planet
 };
 
 ArtificialPlanet::ArtificialPlanet(const PlanetP& orig) :
-		Planet("art", 0, 0, Vec3f(0,0,0), 0, 0, "", "", "", Q_NULLPTR, Q_NULLPTR, Q_NULLPTR, false, true, false, false, "artificial"),
+		Planet("art", 0, 0, Vec3f(0,0,0), 0, 0, "", "", "", "", Q_NULLPTR, Q_NULLPTR, Q_NULLPTR, false, true, false, false, "artificial"),
 		dest(Q_NULLPTR), orig_name(orig->getEnglishName()), orig_name_i18n(orig->getNameI18n())
 {
 	// set parent = sun:

src/core/modules/Comet.cpp

@@ -75,6 +75,7 @@ Comet::Comet(const QString& englishName,
 		  roughness,
 		  atexMapName,
 		  "", // no normalmap.
+		  "", // no horizon map
 		  aobjModelName,
 		  coordFunc,
 		  orbitPtr,

src/core/modules/MinorPlanet.cpp

@@ -36,6 +36,7 @@ MinorPlanet::MinorPlanet(const QString& englishName,
 			 float roughness,
 			 const QString& atexMapName,
 			 const QString& anormalMapName,
+			 const QString& ahorizonMapName,
 			 const QString& aobjModelName,
 			 posFuncType coordFunc,
 			 KeplerOrbit* orbitPtr,
@@ -51,6 +52,7 @@ MinorPlanet::MinorPlanet(const QString& englishName,
 		  roughness,
 		  atexMapName,
 		  anormalMapName,
+		  ahorizonMapName,
 		  aobjModelName,
 		  coordFunc,
 		  orbitPtr,

src/core/modules/MinorPlanet.hpp

@@ -46,6 +46,7 @@ class MinorPlanet : public Planet
 		    float roughness,
 		    const QString& texMapName,
 		    const QString& normalMapName,
+			const QString& horizonMapName,
 		    const QString& objModelName,
 		    posFuncType _coordFunc,
 		    KeplerOrbit *orbitPtr,

src/core/modules/Planet.cpp

@@ -208,6 +208,7 @@ Planet::Planet(const QString& englishName,
 	       float roughness,
 	       const QString& atexMapName,
 	       const QString& anormalMapName,
+		   const QString& ahorizonMapName,
 	       const QString& aobjModelName,
 	       posFuncType coordFunc,
 	       Orbit* anOrbitPtr,
@@ -226,6 +227,7 @@ Planet::Planet(const QString& englishName,
 	  nativeName(""),
 	  texMapName(atexMapName),
 	  normalMapName(anormalMapName),
+	  horizonMapName(ahorizonMapName),
 	  siderealPeriod(0.),
 	  equatorialRadius(radius),
 	  oneMinusOblateness(1.0-oblateness),
@@ -297,6 +299,16 @@ Planet::Planet(const QString& englishName,
 			normalMapFileOrig = normalMapFile;
 		}
 	}
+	if(!horizonMapName.isEmpty())
+	{
+		// TODO: use StelFileMgr::findFileInAllPaths() after introducing an Add-On Manager
+		QString horizonMapFile = StelFileMgr::findFile("textures/"+horizonMapName, StelFileMgr::File);
+		if (!horizonMapFile.isEmpty())
+		{
+			horizonMap = StelApp::getInstance().getTextureManager().createTextureThread(horizonMapFile, StelTexture::StelTextureParams(true, GL_LINEAR, GL_REPEAT));
+			horizonMapFileOrig = horizonMapFile;
+		}
+	}
 	//the OBJ is lazily loaded when first required
 	if(!aobjModelName.isEmpty())
 	{
@@ -354,6 +366,9 @@ void Planet::resetTextures()
 	// restore normal map
 	if (!normalMapFileOrig.isEmpty())
 		normalMap = StelApp::getInstance().getTextureManager().createTextureThread(normalMapFileOrig, StelTexture::StelTextureParams(true, GL_LINEAR, GL_REPEAT));
+
+	if (!horizonMapFileOrig.isEmpty())
+		horizonMap = StelApp::getInstance().getTextureManager().createTextureThread(horizonMapFileOrig, StelTexture::StelTextureParams(true, GL_LINEAR, GL_REPEAT));
 }
 
 void Planet::replaceTexture(const QString &texName)
@@ -2968,6 +2983,7 @@ void Planet::PlanetShaderVars::initLocations(QOpenGLShaderProgram* p)
 	GL(earthShadow = p->uniformLocation("earthShadow"));
 	GL(eclipsePush = p->uniformLocation("eclipsePush"));
 	GL(normalMap = p->uniformLocation("normalMap"));
+	GL(horizonMap = p->uniformLocation("horizonMap"));
 
 	// Rings-specific variables
 	GL(isRing = p->uniformLocation("isRing"));
@@ -3979,6 +3995,8 @@ void Planet::drawSphere(StelPainter* painter, float screenRd, bool drawOnlyRing)
 
 			GL(moonShaderProgram->setUniformValue(moonShaderVars.eclipsePush, push)); // constant for now...
 		}
+		GL(horizonMap->bind(4));
+		GL(moonShaderProgram->setUniformValue(moonShaderVars.horizonMap, 4));
 	}
 
 	if (englishName=="Mars")

src/core/modules/Planet.hpp

@@ -165,6 +165,7 @@ class Planet : public StelObject
 	       float roughness,
 	       const QString& texMapName,
 	       const QString& normalMapName,
+	       const QString& ahorizonMapName,
 	       const QString& objModelName,
 	       posFuncType _coordFunc,
 	       Orbit *anOrbitPtr,
@@ -698,6 +699,7 @@ class Planet : public StelObject
 	QString nativeNameMeaningI18n;   // Can be used in a skyculture
 	QString texMapName;              // Texture file path
 	QString normalMapName;           // Texture file path
+	QString horizonMapName;          // Texture file path
 	RotationElements re;             // Rotation and axis orientation parameters
 	double siderealPeriod;           // sidereal period (Planet year or a moon's sidereal month) [earth days]
 	double equatorialRadius;         // Planet's equatorial radius in AU
@@ -732,6 +734,7 @@ class Planet : public StelObject
 					 // i.e. angle between ascending node of body equator w.r.t. ICRF equator and its prime meridian.
 	StelTextureSP texMap;            // Planet map texture
 	StelTextureSP normalMap;         // Planet normal map texture
+	StelTextureSP horizonMap;        // Planet horizon map texture
 
 	PlanetOBJModel* objModel;               // Planet model (when it has been loaded)
 	QFuture<PlanetOBJModel*>* objModelLoader;// For async loading of the OBJ file
@@ -782,6 +785,7 @@ class Planet : public StelObject
 	// File path for texture and normal map; both variables used for saving original names of files
 	QString texMapFileOrig;
 	QString normalMapFileOrig;
+	QString horizonMapFileOrig;
 
 	std::unique_ptr<QOpenGLVertexArrayObject> sphereVAO;
 	std::unique_ptr<QOpenGLVertexArrayObject> ringsVAO;
@@ -823,6 +827,7 @@ class Planet : public StelObject
 		int earthShadow;
 		int eclipsePush; // apparent brightness push for partial Lunar Eclipse (make bright rim overbright)
 		int normalMap;
+		int horizonMap;
 
 		// Rings-specific variables
 		int isRing;

src/core/modules/SolarSystem.cpp

@@ -984,6 +984,7 @@ bool SolarSystem::loadPlanets(const QString& filePath)
 
 			const bool hidden = pd.value(secname+"/hidden", false).toBool();
 			const QString normalMapName = ( hidden ? "" : englishName.toLower().append("_normals.png")); // no normal maps for invisible objects!
+			const QString horizonMapName = ( hidden ? "" : englishName.toLower().append("_normals.png")); // no normal maps for invisible objects!
 
 			newP = PlanetP(new MinorPlanet(englishName,
 						    pd.value(secname+"/radius", 1.0).toDouble()/AU,
@@ -993,6 +994,7 @@ bool SolarSystem::loadPlanets(const QString& filePath)
 						    pd.value(secname+"/roughness",0.9f).toFloat(),
 						    pd.value(secname+"/tex_map", "nomap.png").toString(),
 						    pd.value(secname+"/normals_map", normalMapName).toString(),
+						    pd.value(secname+"/horizon_map", horizonMapName).toString(),
 						    pd.value(secname+"/model").toString(),
 						    posfunc,
 						    static_cast<KeplerOrbit*>(orbitPtr), // the KeplerOrbit object created previously
@@ -1072,6 +1074,7 @@ bool SolarSystem::loadPlanets(const QString& filePath)
 					       pd.value(secname+"/roughness",0.9f).toFloat(),
 					       pd.value(secname+"/tex_map", "nomap.png").toString(),
 					       pd.value(secname+"/normals_map", englishName.toLower().append("_normals.png")).toString(),
+					       pd.value(secname+"/horizon_map", englishName.toLower().append("_horizon.png")).toString(),
 					       pd.value(secname+"/model").toString(),
 					       posfunc,
 					       static_cast<KeplerOrbit*>(orbitPtr), // This remains Q_NULLPTR for the major planets, or has a KeplerOrbit for planet moons.