Physically based star rendering

[Askaniy]

Updates:
To make the code and the issue easier to maintain, I made a separate repository: https://github.com/Askaniy/CelestiaStarRenderer
For Celestia, I recommend using the Bounded algorithm.
The information in this post is recorded as of October 2023.

What is the problem:

• Celestia star styles don't model their actual difference in brightness physically accurately.
• This became even more noticeable in the branch implementing gamma correction:
  

What were ways to solve:
The problem was already noticed and described by Chris Layrel, 2010 forum discussion is here. In a roundabout way, it led to this branch. But, as onetwothree said, "it's buggy and slow". One of the problems was the squares around bright stars (Gaussians don't have a limit on the angular size, but the shader does):

What I suggest:
I've developed a rendering algorithm that solves this problem. Unfortunately, I've virtually no knowledge of the C/C++, and have limited time. Therefore, I'm creating the issue to document the current developments, perhaps I'll return to this later.

The render below shows three columns of stars renderings of different magnitudes. The first is based on the Chris's formula given on the forum, and in fact consists of two independent Gaussians. The second is an implementation of the photopic point source function (PSF) from this paper. The third column is my attempt to optimize the formula and create a constraint on the angular size (to fix the squares bug). The render is gamma corrected and enlarged three times without interpolation.

Python functions that implement optimized eye's PSF: https://github.com/Askaniy/CelestiaStarRenderer/blob/master/algorithms.py

For example, an exposure of 1 means that Vega occupies a single pixel with no glow, and 5 means that this is true for a star 5 times dimmer than Vega, while Vega itself have a noticeable glow. A general logic looks like this:

br_limit = 1 / (255 * 12.92)
# The linear brightness limit after which a star color converted to 8-bit will be zero.
# It depends on the gamma correction definition, here sRGB is assumed.

exposure = 5 # variable
linear_br = 10**(-0.4 * star_mag) * exposure # scaled brightness measured in Vegas
if linear_br <= br_limit:
    print('Star is too dim to be displayed')
else:
    draw(img, linear_br, rgb, (x,y))

What are known problems:

• Subpixel render will not work well, star is supposed to be strictly in the center. Otherwise, I assume it will lead to unexpected brightness fluctuations unless some kind of the PSF integration by pixel is implemented, but this is expensive. This means that as the DPI increases, the size of stars will decrease. Although, the other star styles have the same problems.
• A bug with squares was fixed by adding a hard limit on the glare angular diameter of 1°. However, this also resulted in very large exposures being unpresentable compared to the original PSF function. Fixed.
• The transition between point rendering and glow is done well for normal use, but breaks down at high angular resolutions (< 0.05-0.03 °/px). Thus, when changing the FOV, the angular size should not be physical, but relative to the FOV.

[Askaniy]

What to do with the faintest apparent magnitude currently in use?

• Way 1: nothing, and calculate exposure for the new star style with faintestMag2exposure.
• Way 2: completely replace faintestMag with exposure in the code, and support other star styles with exposure2faintestMag. Perhaps introducing an exposure will help future developments.

br_limit = 1 / (255 * 12.92)

def faintestMag2exposure(faintestMag):
    return br_limit * 10**(0.4*faintestMag)

def exposure2faintestMag(exposure):
    return 2.5 * np.log10(exposure/br_limit)

[375gnu]

As exposure and faintestMag are mutually convertible it doesn't matter what variable to use. But I think that for user convenience it's better to keep faintestMag.

[Askaniy]

Outdated. See original post for up-to-date information.

---

Improved design of optimized formulas:

• Render diameter increased from one degree to three
• The border will never be sharp
• The equations are even simpler (and faster)

def PSF_fast(theta: float):
    """ Optimized formula of the human eye's point source function """
    if theta < 0.02662:
        return 1 - 53.6*theta + 760*theta*theta
    elif theta < 1.5:
        temp = 1.488 / (theta - 0.012) - 1 # 1.488+0.012=1.5, it's radius
        return 1.1e-5 * temp * temp
    else:
        return 0

However, the large diameter is both an advantage and a disadvantage: now the renderer needs to draw 9 times more pixels. If there is a question of increasing performance, then for br0 < 40 you can use PSF_fast from the first post.

[Askaniy]

Outdated. See original post for up-to-date information.

---

There is an even better solution for the last PSF design: you can calculate the square size depending on the brightness, in the range from 0.012° to 1.5°. For example, for br=1 angle is 0.25°.

def PSF_max_theta(br: float):
    return 0.012 + 1.488 / (np.sqrt(br_limit * 10e6 / (11 * br)) + 1)

and changes in draw_fast():

# glare render on < 3°×3° field
square_half_size = ceil(PSF_max_theta(color0.max()) / degree_per_px)

Why not to calculate the square size for all the angles, not just in 3°?

1. The original PSF fades out more slowly, square sizes would be much larger → less performance
2. The original PSF has very different behavior at small and large angles. In any case, we will need a function that works specifically on small angles.

[Askaniy]

In the future, it may be necessary to link the “convenient linear brightness scale” to the physical one. Let me formalize how I introduce a “convenient scale”: the green color component for a 0 magnitude star illuminates exactly one pixel (exposure=1 everywhere). Let's link the zero magnitude to the Vegan spectrum from the CALSPEC database as a standard, and also assume that the green component of the pixel completely follows the green sensitivity curve of the eye. Correctly multiplying one curve by another, I obtain that the unit of the “convenient brightness scale” corresponds to 3.844750664631218e-11 W/m² in SI.

[Askaniy]

It turned out that the function PSF_max_theta did not work as expected: it returns 1.4999140248716338° even for 1 Vega, the square isn't actually cropped. Also, I tried to fix high (>100 Vega) brightnesses, and fixed them, but with no cropping squares too. I have an idea how to fix everything. However, the basic logic remains the same, so you can try to add to Celestia what I already wrote, and then update the functions when the experiments are successful.

[Askaniy]

After doing a little research, I was able to achieve a single optimized PSF work at any brightness, and, most importantly, that zeroes at an angle, determined from the brightness. Original post updated.

[375gnu]

I have implemented the same in c++ and glsl - https://github.com/CelestiaProject/Celestia/pull/1952/files (windows artifacts at https://github.com/CelestiaProject/Celestia/actions/runs/6658496014).

Problems - square_half_size is too large, in my tests Vega has > 500px. But both the bright star and glare are much smaller. Actually Vega has less than 10px core with the faintest magnitude set to 15. 99% of stars are just 1px. And when i make faintest magnitude higher, square_half_size should always grow, but around 12 it goes down.

[Askaniy]

I used two methods to reduce the square size:

1. Formula coefficients were re-selected to reduce it by 1.75 times, but accuracy compared to the original PSF was slightly reduced too. (At the same time, the formula was further simplified and the overexposure radius was reduced)
2. Cropping pixels below the minimum render brightness also reduced it by almost half. According to OpenGL Wiki: "Implementations are allowed to round the converted integer any way it likes". Therefore, the case of maximum radius was considered, i.e. regular rounding (hopefully no OpenGL implementation uses a ceil rounding, that would be weird).
   

The code became like this:

def PSF_square(theta: float, min_theta: float, max_theta: float, k: float):
    """
    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.
    """
    if theta < min_theta:
        return 1 # overexposed
    elif theta < max_theta:
        brackets = max_theta / theta - 1
        return k * brackets * brackets
    else:
        return 0 # after max_theta function starts to grow again

# Option 2: glare square render
max_br = color0.max()
max_theta = 0.2 * np.sqrt(max_br) # glare radius
k = 3.3e-5 * max_theta**-2.5 # common constant, depending originally on star brightness
min_theta = max_theta / (k**-0.5 + 1)
square_half_size = floor(max_theta / (np.sqrt(0.5 * br_limit / (k * max_br)) + 1) / degree_per_px)
...
color = color0 * PSF_square(theta, min_theta, max_theta, k) 

[375gnu]

For too bright objects, e.g. the Sun when we are inside the Solar system, core is HUGE (radius is something around 1000).

Now glare radius is too small so the image is squary (Sirius at ~1ly with faintest mag = 15):

faintest mag = 14

13:

12:

11:

10:

9:

8:

7:

6:

5:

The star is just a fade pixel.

At 2.80 it disappears. And appears back at 2.60 - when another branch is executed (for 1px stars).

[375gnu]

some thoughts.
in magnitude range from [faintestMag; faintestMag-N], where N is 3 or 4, start is 1 px, with brightness increasing from 1/255 to 1.0.
in (faintestMag-N; faintestMag-M], where M>N, only radius grows. glare is absent (radius == star radius).
in (faintestMag-M; -inf) radius grows but glare radius grows faster.

taking into account limited point sizes and that inside solar systems suns are too bright, at some faintestMag-P, where P > M, radii should grow much slower.

[Askaniy]

It's possible, but there are many problems. The main one is it’s difficult to make the “blur” and “radius” functions visually pleasing. The transition will be sharp and noticeable. Also, for gamma correction, the function must very smoothly go to limited zero. Polynomials will not work, nor will exponential.

I think the current implementation is one of the most optimal possible for glowing within a few tens to hundreds of pixels, but the high brightness problem is real. What about, as soon as a star with a magnitude higher than P appears in the field of view, to force the magnitude limit so as never to exceed P?

[375gnu]

What about, as soon as a star with a magnitude higher than P appears in the field of view, to force the magnitude limit so as never to exceed P?

Sounds ok for me

[Askaniy]

There is another optimization that (for now) can significantly reduce the size of the square: to determine the br_limit depending on the gamma correction.

if srgb:
    br_limit = 1 / (255 * 12.92)
else:
    br_limit = 1 / 255