Merge branch 'askaniys-stars' of github.com:CelestiaProject/Celestia …

…into askaniys-stars

shaders/star_frag.glsl

@@ -1,59 +1,59 @@
 const float degree_per_px = 0.05;
 const float br_limit = 1.0 / (255.0 * 12.92);
+// empirical constants
+const float a = 0.123;
+const float k = 0.0016;
 
 varying vec3 v_color;
-varying vec3 v_max_tetha_hk;
+varying float max_theta;
 varying float pointSize;
+varying float br;
 
-float psf_square(float theta, float min_theta, float max_theta, float h, float k, float b)
+// py: def PSF_Bounded(theta: float, max_theta: float, br_center: float):
+// max_theta is common for all pixels so it's set via `varying` in the vertex shader
+float psf_bounded(float theta, float br_center)
 {
-    // Human eye's point source function, optimized to fit a square.
-    // Lower limit on brightness and angular size: 1 Vega and 0.05 degrees per pixel.
-    // No upper limits.
+    // Human eye's point source function from the research by Greg Spencer et al., optimized to fit a square.
+    // Lower limit on brightness and angular size: 1 Vega and 0.05 degrees per pixel. No upper limits.
 
-    if (theta < min_theta)
-        return 1.0; // overexposed
+    // py: if theta == 0:
+    if (theta == 0)
+    // py:  return br_center
+        return br_center; // the center is always overexposed (zero division error)
 
+    // py:  elif theta < max_theta:
     if (theta < max_theta)
     {
-        float brackets = b / (theta - h) - 1.0;
-        return brackets * brackets / k;
+        // py: brackets = max_theta / theta - 1
+        float brackets = max_theta / theta - 1.0;
+        // py: return k * brackets * brackets
+        return k * brackets * brackets;
     }
 
+    // py: return 0. # after max_theta function starts to grow again
     return 0.0; // after max_theta function starts to grow again
 }
 
-/*
-float psf_fullscreen(float theta, float min_theta):
-{
-    // Human eye's point source function, optimized to be a full-screen shader.
-    // The price to pay for simplification is a brightness reduction compared to the original PSF.
-
-    if (theta2 < min_theta)
-        return 1; // overexposed
-
-    return 4.43366571e-6 / theta;
-}
-*/
-
 void main(void)
 {
-    float max_theta = v_max_tetha_hk.x;
     if (max_theta == -1.0)
     {
+        // just a one pixel star
+        // py: arr[center[1], center[0]] += scaled_color
         gl_FragColor = vec4(v_color, 1.0);
     }
     else
     {
-        float h = v_max_tetha_hk.y;
-        float k = v_max_tetha_hk.z;
-
-        float b = max_theta - h;
-        float min_theta = h + b / (sqrt(k) + 1.0);
-
+        // Option 2: glare square render
+        // py: theta = np.sqrt(xx*xx + yy*yy) * degree_per_px # array of distances to the center
+        // in fragment shader all points have virtual dimension 1x1, so gl_PointCoord has a value from [0; 1]
         vec2 offset = (gl_PointCoord.xy - vec2(0.5)) * pointSize;
         float theta = length(offset) * degree_per_px;
-
-        gl_FragColor = vec4(v_color * psf_square(theta, min_theta, max_theta, h, k, b), 1.0);
+        // py: glow_bw = PSF_Bounded(theta, max_theta, br0) # in the [0, 1] range, like in Celestia
+        float glow_bw = psf_bounded(theta, br);
+        // py: glow_colored = color * np.repeat(np.expand_dims(glow_bw, axis=2), 3, axis=2) # scaling
+        vec3 glow_colored = v_color * glow_bw;
+        // py: arr[center[1]+y_min:center[1]+y_max, center[0]+x_min:center[0]+x_max] += glow_colored
+        gl_FragColor = vec4(glow_colored, 1.0);
     }
 }

shaders/star_vert.glsl

@@ -1,41 +1,90 @@
 
 const float degree_per_px = 0.05;
-const float br_limit = 1.0 / (255.0 * 12.92);
+// empirical constants
+const float a = 0.123;
+const float k = 0.0016;
+const float exposure = 1.0;
 
-varying vec3 v_color; // 12
-varying vec3 v_max_tetha_hk; // 24
-varying float pointSize; // 28
+// py: max_square_size = 512 # px
+const float max_square_size = 256.0;
+// py: max_br = (degree_per_px * max_square_size / algorithms.a)**2 / (2*np.pi)
+const float max_br = pow((degree_per_px * max_square_size / a), 2.0) / (2.0 * 3.141592653);
 
-uniform vec2 viewportSize;
+varying vec3 v_color;
+varying float max_theta;
+varying float pointSize;
+varying float br;
 
 attribute vec4 in_Position;
 attribute vec3 in_Color;
-attribute float in_PointSize; // scaled brightness measured in Vegas
+attribute float in_PointSize;
 
+const float color_saturation_limit = 0.1; // The ratio of the minimum color component to the maximum
+
+//! Normalizes the color by its green value and corrects extreme saturation
+// py: def green_normalization(color: np.ndarray):
+vec3 green_normalization(vec3 color)
+{
+    // py: color /= color.max()
+    // color /= max(color.r, max(color.g, color.b)); // we do this in XYZRGBConverter::convertUnnormalized()
+
+    // py: delta = color_saturation_limit - color.min()
+    float delta = color_saturation_limit - min(color.r, min(color.g, color.b));
+
+    // py: if delta > 0:
+    if (delta > 0)
+    {
+        // py: color += delta * (1-color)**2 # desaturating to the saturation limit
+        vec3 diff = vec3(1.0) - color;
+        color += delta * (diff * diff); // desaturating to the saturation limit
+    }
+    // py: return color / color[1]
+    return color / color.g;
+}
 
 void main(void)
 {
-    float linearBr = pow(10.0, 0.4f * in_PointSize) * br_limit;
-    vec3 color0 = in_Color * (linearBr / in_Color.g); // scaling on brightness and normalizing by green channel
+    // py: linear_br = 10**(-0.4 * star_mag) * exposure # scaled brightness measured in Vegas
+    // here i use +0.4 because `in_PointSize` is not the actual magnitude but `faintest` - `actual`
+    float br0 = pow(10.0, 0.4 * in_PointSize) * exposure;
+
+    // py: color = auxiliary.green_normalization(color0)
+    vec3 color = green_normalization(in_Color);
 
-    vec3 check_vec = step(vec3(1.0), color0); // step(edge, x) - For element i of the return value, 0.0 is returned if x[i] < edge[i], and 1.0 is returned otherwise.
-    float check = check_vec.x + check_vec.y + check_vec.z;
-    if (check == 0.0)
+    // py: scaled_color = color * br0
+    vec3 scaled_color = color * br0;
+
+    // py: if np.all(scaled_color < 1):
+    if (all(lessThan(scaled_color, vec3(1.0))))
     {
+        // we set color in the fragment shader (using v_color) so here we just set point size to 1px
         pointSize = 1.0;
-        v_max_tetha_hk = vec3(-1.0);
+        // use max_theta == -1 as as signal that the point size is 1px
+        max_theta = -1.0;
+
+        v_color = scaled_color;
     }
     else
     {
-        float max_theta = 0.33435822702992773 * sqrt(max(color0.r, max(color0.g, color0.b))); // glare radius
-        pointSize =  2.0 * ceil(max_theta / degree_per_px);
+        // py: br = np.arctan(br0 / max_br) * max_br # dimmed brightness
+        br = atan(br0 / max_br) * max_br; // dimmed brightness
+        // py: max_theta = a * np.sqrt(br) # glow radius
+        max_theta = a * sqrt(br); // glow radius
+        // py: half_sq = floor(max_theta / degree_per_px)
+        float half_sq = max_theta / degree_per_px;
+        // py: x_min = -min(half_sq, center[0])
+        // py: x_max = min(half_sq+1, width-center[0])
+        // py: y_min = -min(half_sq, center[1])
+        // py: y_max = min(half_sq+1, hight-center[1])
+        // py: x = np.arange(x_min, x_max)
+        // py: y = np.arange(y_min, y_max)
+        // py: xx, yy = np.meshgrid(x, y)
+        // we just set a point size. all iteration over every px is done in the fragment shader
+        pointSize = 2.0 * half_sq - 1.0;
 
-        float h = 0.0082234880783653 * pow(max_theta, 0.7369983254906639); // h, k, b - common constants, depending originally on star brightness
-        float k = 38581.577272697796 * pow(max_theta, 2.368787717957141);
-        v_max_tetha_hk = vec3(max_theta, h, k);
+        v_color = color;
     }
 
     gl_PointSize = pointSize;
-    v_color = color0;
     set_vp(in_Position);
 }

src/celengine/pointstarrenderer.cpp

@@ -37,6 +37,10 @@ PointStarRenderer::PointStarRenderer() :
 {
 }
 
+//const float br_limit = 1.0f / (255.0f * 12.92f);
+const float br_limit = 1.0f / 255.0f;
+const float exposure = 1.0f;
+
 void PointStarRenderer::process(const Star& star, float distance, float appMag)
 {
     if (distance > distanceLimit)