#!/usr/bin/env python # coding: utf-8 # https://journals.ametsoc.org/view/journals/bams/79/5/1520-0477_1998_079_0831_opoaac_2_0_co_2.xml?tab_body=pdf - for aerosols and clouds # # https://iopscience.iop.org/article/10.1088/1757-899X/1269/1/012011/pdf - for Rayleigh extinction # # https://sci-hub.se/10.1364/josa.48.1018_1 - temperature dependance Rayleigh scattering # # https://sci-hub.se/https://doi.org/10.1364/JOSA.47.000176 - Rayleigh coefficients (index of refraction) - Penndorf # # We are going to use the real rayleigh scattering coefficient: # # $\beta_R(h, \lambda) = \frac{8\pi(n^2 - 1)^2}{2N\lambda^4} e^{-\frac{h}{H_R}}$ # # n = index of refraction, these are calculated in the 4th paper, temperature dependance # # N = molecular density, dependent on humidity # # $\lambda$ = wavelength, but we only care for RGB (440, 550, 680 nm) # # From Penndorf, the dispersion equation is: # # $(n_s - 1) \times 10^{-8} = 6432.8 + \frac{2949810}{146 - v^2} + \frac{25540}{41 - v^2}$ # # It is important to note that $v=1/f$ which is measured in um # # Once we have calculated our dispersion index of refraction, we can see what our new value for n will be (it is dependant on temperature and pressure) # # $$ # (n - 1) = (n_s - 1) \left( \frac{(1 + \alpha t_s)}{(1 + \alpha t)} \right) \frac{p}{p_s} # $$ # # Here, we are using $n_s$ from the dispersion equation, $t_s$ is 15 degrees Celsius, and $p_s$ is 760 mmHg. $t$ represents the actual temperature, and $p$ represents the actual pressure. $\alpha$ is the coefficient of linear expansion of air = 0.00367 in units of inverse celcius. We see that the influence of air pressure really is critical # # # We also need to find N, the molecular density (depends on temperature, humidity, pressure): We are going to use the value at the altitude where our person is standing and then we scale it using the scale factor. Penndorf uses N_s = 2.54743e19, but i am not sure if it is applicable to our scenario. # # # # # Ok, so things work... but they are not super realistic. I think we should be updating the dispersion equation because Penndorf's paper is super old. We want the value of the index of refraction at the location of the observer - this is going to depend on the users inputs. # # https://refractiveindex.info/?shelf=other&book=air&page=Ciddor (1996) # # New calculation for refractive index (where lambda is in microns) # # This applies for : Standard air: dry air at 15 °C, 101.325 kPa and with 450 ppm CO2 content. # # $n - 1 = \frac{0.05792105}{238.0185 - \lambda^{-2}} + \frac{0.00167917}{57.362 - \lambda^{-2}}$ # # # # Something really cool here is that we can use the color of the sunset to predict the weather: "Red at night, sailors delight. Red in morning, sailors take warning" # When there is a bright red sunset at night, it means there is likely a high pressure cell off to the west. With west to east weather patterns, this means a high pressure cell is about to move in, high pressure = pleasant weather, "sailors delight". When there is a red sunrise in the east, it means the high pressure cell has already passed and a low pressure cell is likely to replace it. Low pressure = poor weather, "sailors take warning". # # Low pressure: 720mmHg # High pressure: 780 mmHg import numpy as np import math import imageio import time import random import cython """ This class allows us to manipulate vectors Not all of the functions are used yet - but they might be later when we start adding things """ @cython.cclass class Vec3: v = cython.declare(cython.float[3],visibility="public") def __cinit__(self, x:cython.float =0, y:cython.float=0, z:cython.float=0): self.v:cython.float[3] = [x, y, z] def __init__(self, x:cython.float =0, y:cython.float=0, z:cython.float=0): self.v:cython.float[3] = [x, y, z] def __add__(self, other ): return Vec3(self.v[0] + other.v[0],self.v[1] + other.v[1],self.v[2] + other.v[2]) def __sub__(self, other): return Vec3(self.v[0] - other.v[0],self.v[1] - other.v[1],self.v[2] - other.v[2]) def __mul__(self, other): if isinstance(other, Vec3): return Vec3(self.v[0] * other.v[0],self.v[1] * other.v[1],self.v[2] * other.v[2]) #elif isinstance(other, (int, float)): return Vec3(self.v[0] * other,self.v[1] * other,self.v[2] * other) #raise TypeError("Can only multiply by a Vec3 or scalar (int or float).") def __truediv__(self, scalar): return Vec3(self.v[0] / scalar,self.v[1] / scalar,self.v[2] / scalar) def length(self)->cython.float: return math.sqrt(self.v[0]**2 + self.v[1]**2 + self.v[2]**2) def normalize(self): length = self.length() return self / length if length > 0 else Vec3() def to_tuple(self): return (self.v[0],self.v[1],self.v[2]) def __repr__(self): return f"Vec3({self.v[0]}, {self.v[1]}, {self.v[2]})" @cython.cfunc def dot(self, other): return self.v[0]* other.v[0]+self.v[1]* other.v[1]+self.v[2]* other.v[2] """ Here we will be calculating what the value for n is based on T, P, but we will also be taking dispersion into account - remember, we are considering 440 nm, 550 nm, 680 nm These equations match the tabulated results made by Penndorf """ wavelengths = Vec3(0.68, 0.55, 0.44) # These are the wavelengths in um # The n_s values gives us the different refractive indices for the different wavelenghts # The n values give us the new refractive indices depending on the temperature and pressure @cython.cfunc def refraction_calculator(T:cython.int, P:cython.int) -> tuple[Vec3, Vec3]: alpha:cython.float = 0.00367 # in inverse Celsius t_s:cython.int = 15 # degrees Celsius p_s:cython.int = 760 # mmHg n_s_values = [] # for debugging / comparing to the internet values n_values = [] n_s:cython.double n:cython.double i:cython.int for i in [wavelengths.v[0], wavelengths.v[1], wavelengths.v[2]]: n_s = 1 + (0.05792105 / (238.0185 - i**-2) + 0.00167917 / (57.362 - i**-2)) n = ((n_s - 1) * ((1 + alpha * t_s) / (1 + alpha * T)) * (P / p_s)) + 1 n_s_values.append(n_s) n_values.append(n) return Vec3(*n_s_values), Vec3(*n_values) @cython.cfunc def beta(n_values:Vec3): N_s:cython.double = 2.54743e7 # um^-3 beta_values = [] for i, n in enumerate(n_values.v): wavelength = wavelengths.v[i] # Corresponding wavelength value beta_1 = ((8 * math.pi * (n**2 - 1)**2) / ((2 * N_s) * (wavelength)**4))*1e6 # The whole thing should be in m^-1 to match with the rest of the code beta = (beta_1*10) beta_values.append(beta) return Vec3(beta_values[0],beta_values[1],beta_values[2]) # Example usage, normal is 25, 750 T:cython.int = 5 # Temperature in degrees Celsius P:cython.int = 760 # Pressure in mmHg n_s_values:Vec3 n_values:cython.double n_s_values, n_values = refraction_calculator(T, P) beta_values = beta(n_values) #print("n_s values:", refraction_calculator_result.v) #print("n values:", n_values.v) #print("Beta values:", beta_values.v) # We will have a scattering enhancement factor for our scattering values kInfinity = float('inf') angle_degrees:cython.int = 90 # Sun direction in degrees (0 to 90, noon to sunset) angle_radians:cython.double = math.radians(angle_degrees) # Convert to radians sunDir:Vec3 = Vec3(0, math.cos(angle_radians), -math.sin(angle_radians)) """ The values for betaR depend on many things For one, they scale proportionally to barometric pressure It also varies according to temperature """ betaR = Vec3(*beta_values.v) # Keeps it as a Vec3 #betaR = Vec3(3.8e-6, 13.5e-6, 33.1e-6) # R, G, B (the smaller the value, the less it will scatter during the day) in m-1 # 3.8 e-6, 13.5 e-6, 33.1 e-6 - default # The order should be 680 nm, 550 nm, 440 nm """ There are many "categories" of aerosols and each has their own extinction coefficient (I removed the factor of 1.1) Continental Clean - 26e-6 Continental Average - 75e-6 Continental Polluted - 175e-6 Urban - 353e-6 Desert - 145e-6 Maritime Clean - 90e-6 Maritime Polluted - 115e-6 Maritime Tropic - 43e-6 Arctic - 23e-6 Antarctic - 11e-6 """ #betaM = Vec3(21e-6, 21e-6, 21e-6) # We need to have the Mie scattering be the same in all the directions # The greater the value of beta, the smaller the Mie scattering point (responsable for the halo around the sun) # If there is more pollution, we get a larger halo and the colors of the sunset become desaturated (more hazy) # Default is 21e-6 earthRadius:cython.int = 6360e3 #m atmosphereRadius:cython.int = 6420e3 #m (60 km higher than earth) Hr:cython.int = 7994 Hm:cython.int = 1200 # The direction will change for each pixel @cython.cfunc def computeIncidentLight(direction:Vec3,betaM,g,observerEarthRadius) -> tuple[float,float,float]: tmin:cython.float=0 tmax:cython.float=kInfinity # We can change the origin position if we want, but for now it is 1 meter above the surface orig:Vec3 = Vec3(0, observerEarthRadius + 1, 0) t0:cython.float t1:cython.float t0, t1 = raySphereIntersect(orig, direction, atmosphereRadius) if t1 < 0: return Vec3(0, 0, 0) if t0 > tmin and t0 > 0: tmin = t0 if t1 < tmax: tmax = t1 numSamples:int = 16 numSamplesLight:int = 8 segmentLength = (tmax - tmin) / numSamples tCurrent = tmin sumR = Vec3() sumM = Vec3() opticalDepthR = 0 opticalDepthM = 0 mu:cython.float = direction.dot(sunDir) # This is the cosine of the angle between the direction vector (V) and the sun Direction phaseR = (3 * (1 + (mu * mu))) / (16 * math.pi) #g = 0.76 phaseM = 3 / (8 * math.pi) * ((1 - g * g) * (1 + mu * mu)) / ((2 + g * g) * ((1 + g * g - 2 * g * mu) ** 1.5)) i:cython.int for i in range(numSamples): samplePosition = orig + direction * (tCurrent + segmentLength * 0.5) height = samplePosition.length() - observerEarthRadius hr = math.exp(-height / Hr) * segmentLength hm = math.exp(-height / Hm) * segmentLength opticalDepthR += hr opticalDepthM += hm # Sample position is the start, but it should be in the direction of the sun t0Light, t1Light = raySphereIntersect(samplePosition, sunDir, atmosphereRadius) if t1Light < 0: continue segmentLengthLight = (t1Light - t0Light) / numSamplesLight tCurrentLight = 0 opticalDepthLightR = 0 opticalDepthLightM = 0 j:cython.int for j in range(numSamplesLight): samplePositionLight = samplePosition + (sunDir * (tCurrentLight + segmentLengthLight * 0.5)) heightLight = samplePositionLight.length() - observerEarthRadius if heightLight < 0: break opticalDepthLightR += (math.exp(-heightLight / Hr) * segmentLengthLight) opticalDepthLightM += (math.exp(-heightLight / Hm) * segmentLengthLight) tCurrentLight += segmentLengthLight if j == numSamplesLight - 1: tau = (betaR * (opticalDepthR + opticalDepthLightR)) + (betaM * (opticalDepthM + opticalDepthLightM)) attenuation = Vec3(math.exp(-tau.v[0]), math.exp(-tau.v[1]), math.exp(-tau.v[2])) sumR += (attenuation * hr) sumM += (attenuation * hm) tCurrent += (segmentLength) # Changing the * 20 number just changes the intensity os the light, it does not much change the colors themselves # Well the colors kinda change, but more they just get super pale or super not pale. 20 should be fine I guess final_color = (sumR * betaR * phaseR + sumM * betaM * phaseM) * 20 return final_color.to_tuple() @cython.cfunc def raySphereIntersect(orig:Vec3, direction:Vec3, radius:cython.double) -> tuple[cython.float,cython.float]: A:cython.float = direction.v[0]**2 + direction.v[1]**2 +direction.v[2]**2 B:cython.float = 2 * (direction.v[0]*orig.v[0] +direction.v[1]*orig.v[1] +direction.v[2]*orig.v[2] ) C:cython.float = orig.v[0]**2 + orig.v[1]**2 +orig.v[2]**2 - radius ** 2 t0: cython.float t1: cython.float t0, t1 = solveQuadratic(A, B, C) if t1 < 0: return (kInfinity, kInfinity) #if t0 > t1: # t0, t1 = t1, t0 return (t0, t1) @cython.cfunc def solveQuadratic(a:cython.float, b:cython.float, c:cython.float) -> tuple[cython.float,cython.float]: if b == 0: if a == 0: return (0, 0) x1 = 0 x2 = math.sqrt(-c / a) return (x1, x2) discr = b ** 2 - 4 * a * c if discr < 0: return (kInfinity, kInfinity) discr:cython.float q:cython.float discrs = math.sqrt(discr) q = (-0.5 * (b - discrs)) if b < 0 else (-0.5 * (b + discrs)) return (q/a, c/q) @cython.cfunc def renderSkydome(filename): width:cython.int height:cython.int width, height = 256,256 image = np.zeros((height, width, 3), dtype=float) start_time = time.time() # Start timing j:cython.int i:cython.int x:cython.float y:cython.float z2:cython.float phi:float theta:float direction:Vec3 for j in range(height): for i in range(width): x = 2 * (i + 0.5) / (width - 1) - 1 y = 2 * (j + 0.5) / (height - 1) - 1 z2 = x * x + y * y if z2 <= 1: phi = math.atan2(y, x) theta = math.acos(1 - z2) # This changes for each pixel direction = Vec3(math.sin(theta) * math.cos(phi), math.cos(theta), math.sin(theta) * math.sin(phi)) color = computeIncidentLight(direction) #color = computeIncidentLight(direction) # Assign the clipped color directly to the image array image[j][i][0] = np.clip(color, 0, 1)[0] image[j][i][1] = np.clip(color, 0, 1)[1] image[j][i][2] = np.clip(color, 0, 1)[2] # Print elapsed time after each row elapsed_time = time.time() - start_time print(f"Rendering row {j + 1}/{height}, elapsed time: {elapsed_time:.2f} seconds") #print(f"Rendering row {j + 1}/{height}") # Save result to a PNG image image = np.clip(image, 0, 1) * 255 # change 255 imageio.imwrite(filename, image.astype(np.uint8)) def renderFromCamera(filename,betaM,g,altitude): observerEarthRadius = earthRadius + altitude width:cython.int height:cython.int width, height = 100, 100 x:cython.int y:cython.int rayx:cython.float rayy:cython.float image = np.zeros((height, width, 3), dtype=np.float32) aspectRatio = width / float(height) fov = 65 # field of view angle = math.tan(fov * np.pi / 180 * 0.5) # How much the camera can see start_time = time.time() # Start timing numPixelSamples:cython.int = 4 for y in range(height): for x in range(width): # for m in range(numPixelSamples): # for n in range(numPixelSamples): # rayx = (2 * (x + (m + random.uniform(0, 1)) / numPixelSamples) / float(width) - 1) * aspectRatio * angle # rayy = (1 - (y + (n + random.uniform(0, 1)) / numPixelSamples) / float(height) * 2) * angle rayx = (2 * x / float(width) - 1) * aspectRatio * angle rayy = (1 - y / float(height) * 2) * angle #This changes for each pixel, it is the direction we are looking in direction = Vec3(rayx, rayy, -1).normalize() color = computeIncidentLight(direction,betaM,g,observerEarthRadius) #image[y, x] = np.array(color) image[y][x] = np.clip(color, 0, 1) #image[y, x] /= numPixelSamples * numPixelSamples # Average color elapsed_time = time.time() - start_time print("what") print(f"Renderinggg row {y + 1}/{height}, elapsed time: {elapsed_time:.2f} seconds") image = np.clip(image, 0, 1) * 255 return image.astype(np.uint8) #imageio.imwrite(filename, image.astype(np.uint8)) #renderFromCamera("camera_render.png") #renderSkydome("highpress_15C.png") #renderFromCamera("highpress_camera.png")