From d5768a4bfe8f3162be7b38dac538297900e1705b Mon Sep 17 00:00:00 2001 From: Technoduck Date: Wed, 11 Dec 2024 17:22:35 -0500 Subject: [PATCH] fixed altitude --- Midterm_submittion.py | 445 ------------------------------------------ app.py | 52 ++--- astropy_test.py | 21 -- skydome.py | 26 +-- 4 files changed, 26 insertions(+), 518 deletions(-) delete mode 100644 Midterm_submittion.py delete mode 100644 astropy_test.py diff --git a/Midterm_submittion.py b/Midterm_submittion.py deleted file mode 100644 index dfbf88d..0000000 --- a/Midterm_submittion.py +++ /dev/null @@ -1,445 +0,0 @@ -#!/usr/bin/env python -# coding: utf-8 - - -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.double =0, y:cython.double=0, z:cython.double=0): - self.v:cython.float[3] = [x, y, z] - - def __init__(self, x:cython.double =0, y:cython.double=0, z:cython.double=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).") - - @cython.cfunc - def __truediv__(self, scalar): - if scalar == 0: - raise ValueError("Cannot divide by zero.") - return Vec3(*(self.v[i] / scalar for i in range(3))) - - def length(self)->cython.float: - return math.sqrt(self.v[0]**2 + self.v[1]**2 + self.v[2]**2) - - @cython.cfunc - 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) - - -# In[7]: - - -""" -Here we are able to change the sun's direction - once the GUI is in place, this will be done from there -""" - -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 betaR values come from the previous cell -""" - -#betaR = Vec3(3.8e-6, 13.5e-6, 33.1e-6), this is what the Nishita paper used -betaR = Vec3(*beta_values.v) # Keeps it as a Vec3 - - -""" -There are many "categories" of aerosols and each has their own extinction coefficient (I removed the factor of 1.1) - -Continental Clean - 26e-6 g = 0.709 -Continental Average - 75e-6 g = 0.703 -Continental Polluted - 175e-6 g = 0.698 - -Urban - 353e-6 g = 0.689 - -Desert - 145e-6 g = 0.729 - -Maritime Clean - 90e-6 g = 0.772 -Maritime Polluted - 115e-6 g = 0.756 -Maritime Tropic - 43e-6 g = 0.774 - -Arctic - 23e-6 g = 0.721 - -Antarctic - 11e-6 g = 0.784 - - -""" -betaM = Vec3(11e-6, 11e-6, 11e-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 - -# For the value of betaM, we also need to change the anisotropy factor. This is in the same paper. -g = 0.784 - -""" -The atmosphere model we use is as follows, with the radii and the scale heights -""" -earthRadius:cython.int = 6360e3 #m -atmosphereRadius:cython.int = 6420e3 #m (60 km higher than earth) -Hr:cython.int = 7994 -Hm:cython.int = 1200 - - -# In[8]: - - -""" -This is the compute incident light funciton from Nishita's paper -We translated it to Python and then are using cython to make it go faster - -This is all explained in more detail in the submitted PDF -""" - -@cython.cfunc -def computeIncidentLight(direction:Vec3) -> 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, earthRadius + 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 - - """ - Phase functions - the anisotropy factor depends on the Mie scattering we have chosen, it is defined in the previous cell - """ - phaseR = (3 * (1 + (mu * mu))) / (16 * math.pi) - 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() - earthRadius - 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() - earthRadius - 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, this is the "magic number" Nishita used - final_color = (sumR * betaR * phaseR + sumM * betaM * phaseM) * 20 - return final_color.to_tuple() - - -# In[9]: - - -""" -These are complimentary functions that are used to solve the intensity of light. Again, from nishita's paper -""" - -@cython.cfunc -def raySphereIntersect(orig:Vec3, direction:Vec3, radius:cython.double) -> tuple[cython.float,cython.float]: - A:cython.double = direction.v[0]**2 + direction.v[1]**2 +direction.v[2]**2 - B:cython.double = 2 * (direction.v[0]*orig.v[0] +direction.v[1]*orig.v[1] +direction.v[2]*orig.v[2] ) - C:cython.double = 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.double, b:cython.double, c:cython.double) -> 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) - - q = (-0.5 * (b - np.sqrt(discr))) if b < 0 else (-0.5 * (b + np.sqrt(discr))) - x1 = q / a - x2 = c / q - return (x1, x2) - - -# In[10]: - - -""" -Rendering Skydome image -""" - -@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)) - - - -# In[11]: - - -""" -Rendering camera image -""" - - -def renderFromCamera(filename: str): - width: cython.int = 252 - height: cython.int = 252 - image = np.zeros((height, width, 3), dtype=np.float32) - - aspectRatio: cython.float = width / float(height) - fov: cython.float = 65.0 # Field of view - angle: cython.float = np.tan(fov * np.pi / 180 * 0.5) # Camera's view angle - - start_time = time.time() # Start timing - - numPixelSamples: cython.int = 4 - y: cython.int - x: cython.int - m: cython.int - n: cython.int - rayx: cython.float - rayy: cython.float - direction: Vec3 - color: tuple # Assuming computeIncidentLight returns a tuple - - for y in range(height): - for x in range(width): - sum_color = np.zeros(3, dtype=np.float32) # Use numpy array for accumulation - - for m in range(numPixelSamples): - for n in range(numPixelSamples): - # Compute ray direction with jitter - 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 - - # Create the direction vector and normalize it - direction = Vec3(rayx, rayy, -1).normalize() - color = computeIncidentLight(direction) # Assuming this returns a tuple - - # Accumulate color values - sum_color[0] += color[0] - sum_color[1] += color[1] - sum_color[2] += color[2] - - # Average the accumulated color and clip to [0, 1] - image[y, x] = np.clip(sum_color / (numPixelSamples * numPixelSamples), 0, 1) - - # Print elapsed time after each row - elapsed_time = time.time() - start_time - print(f"Rendering row {y + 1}/{height}, elapsed time: {elapsed_time:.2f} seconds") - - # Save the result to a PNG image - image = np.clip(image, 0, 1) * 255 - imageio.imwrite(filename, image.astype(np.uint8)) - - - - -# In[12]: - - -""" -Testing -""" - -renderSkydome("test.png") - - -# In[ ]: - - - - diff --git a/app.py b/app.py index 7b55582..8b478eb 100644 --- a/app.py +++ b/app.py @@ -16,16 +16,11 @@ def get_api_params(coords): # create tkinter window root_tk = tkinter.Tk() -# set other tile server (standard is OpenStreetMap) -# map_widget.set_tile_server("https://mt0.google.com/vt/lyrs=m&hl=en&x={x}&y={y}&z={z}&s=Ga", max_zoom=22) # google normal -# map_widget.set_tile_server("https://mt0.google.com/vt/lyrs=s&hl=en&x={x}&y={y}&z={z}&s=Ga", max_zoom=22) # google satellite - -def printer(): - print("button") +#environment setting window definition environment_tk = tkinter.Tk() - +#Prop for a colo changing button class EnvButton(tkinter.Button): def __init__(self,*args,**kwargs): tkinter.Button.__init__(self, *args, **kwargs) @@ -40,9 +35,10 @@ class EnvButton(tkinter.Button): self['bg'] = self.default_bg_color - +# main window class of the rendered Image window class renderedImageZoom: def __init__(self,root,image_window,aerosol_window): + #setting of the default parameters self.root = root self.root.geometry(f"{1000}x{700}") self.root.title("map_view_simple_example.py") @@ -72,9 +68,10 @@ class renderedImageZoom: self.image_window.title("Simulated Sunset") self.image_window.config(width=256,height=256) - #image = skydome.renderFromCamera(coords) + #creating the image canvas for displaying the rendered image self.canvas = tkinter.Canvas(self.image_window,bg="white") self.canvas.pack(fill=tkinter.BOTH,expand=True) + #separate function for ease of creating the environment options window self.environment_window_filler() @@ -102,9 +99,7 @@ class renderedImageZoom: # #Antarctic - 11e-6 g = 0.784 - -# button_list = [buton_cont_clean,buton_cont_avr,buton_cont_poll,buton_urban,buton_desert,buton_mar_clean,buton_mar_poll,button_mar_tro,button_arctic,buton_antarctic] - +#code for changing the color of the active button def set_active(index): if index != self.curr_env: @@ -118,7 +113,7 @@ class renderedImageZoom: buton_cont_clean['bg'] = 'cadetblue1' self.env_setter(26e-6,0.709) - + #setup of all the buttons for the environment buton_cont_avr = EnvButton(self.aerosol_window,text="Continental Average ",command=lambda: [self.env_setter(75e-6,0.793),set_active(1)]); buton_cont_avr.grid(column=0,row=2,padx=10,pady=5) @@ -154,12 +149,11 @@ class renderedImageZoom: self.set_temp(value) label.config(text=str(sorted([-270,int(value),100]))) - + #frames for placement of temperature, pressure, and sun angle buttons and info t_frame = tkinter.Frame(self.aerosol_window) text_t = tkinter.Text(t_frame,height=1,width=10) submit_t_button = tkinter.Button(t_frame,text="set temperature", command=lambda: [self.set_temp(text_t.get("1.0","end-1c")),t_label.config(text=str(sorted([-270,int(text_t.get("1.0","end-1c")),100])[1])+" C")]) - #submit_t_button = tkinter.Button(t_frame,text="set temperature",command=lambda: [print("what"),t_label.config(text="bruh")]) t_label = tkinter.Label(t_frame,text="5 C") text_t.pack(side=tkinter.LEFT) t_label.pack(side=tkinter.RIGHT) @@ -189,7 +183,6 @@ class renderedImageZoom: s_frame.grid(column=0,row=9,columnspan=2,padx=20) -## Altitude a_frame = tkinter.Frame(self.aerosol_window) text_a = tkinter.Text(a_frame,height=1,width=10) submit_a_button = tkinter.Button(a_frame,text="Override altitude",command=lambda: self.set_altitute(int(text_a.get("1.0","end-1c")))) @@ -217,6 +210,7 @@ class renderedImageZoom: r_frame.grid(column=0,row=12,columnspan=2,padx=20) + #red and blue average calculations b_frame = tkinter.Frame(self.aerosol_window) text_b_l = tkinter.Text(b_frame,height=1,width=10) text_b_u = tkinter.Text(b_frame,height=1,width=10) @@ -229,10 +223,6 @@ class renderedImageZoom: submit_b_button.pack(side=tkinter.RIGHT) b_frame.grid(column=0,row=13,columnspan=2,padx=20) - #self.img = Image.fromarray(image,mode="RGB") - #self.tk_image = ImageTk.PhotoImage(width=256,height=256,image=self.img) - #img = ImageTk.PhotoImage(Image.open("highpress_camera.png")) - def cam_toggle(): if cam_fish_toggle.config('text')[-1] == "cam view": cam_fish_toggle.config(text="fisheye") @@ -242,7 +232,7 @@ class renderedImageZoom: self.fisheye = False - #self.canvas.create_image(0,0, anchor="center", image=self.tk_image) + #rendered image window buttons for zoom in/out, save image, switch render mode, and render zoom_in_button = tkinter.Button(self.image_window, text="Zoom In", command=self.zoom_in) zoom_out_button = tkinter.Button(self.image_window, text="Zoom Out", command=self.zoom_out) render_button = tkinter.Button(self.image_window,text="render",command=self.render) @@ -257,6 +247,7 @@ class renderedImageZoom: self.canvas.bind("",self.zoom_in) self.canvas.bind("", self.zoom_out) + #setter functions def set_temp(self,temp): if temp is None: self.temperature = 0 @@ -275,15 +266,16 @@ class renderedImageZoom: if angle is None: angle_degrees:cython.int = 90 # Sun direction in degrees (0 to 90, noon to sunset) elif int(angle) > 90 or int(angle) < 0: - angle_degrees:cython.int = 90 # Sun direction in degrees (0 to 90, noon to sunset) + angle_degrees:cython.int = 90 else: - angle_degrees:cython.int = int(angle) # Sun direction in degrees (0 to 90, noon to sunset) + angle_degrees:cython.int = int(angle) angle_radians:cython.double = math.radians(angle_degrees) # Convert to radians sunDir:skydome.Vec3 = skydome.Vec3(0, math.cos(angle_radians), -math.sin(angle_radians)) self.sunDir = sunDir + #setters have reasonable (ish) limiters on the values def set_press(self,press): if press is None: self.pressure = 720 @@ -310,10 +302,8 @@ class renderedImageZoom: for line in slice: for pixel in line: if sum(pixel) == 0: - print(pixel) red += 0 else: - print(pixel) red += pixel[0]/sum(pixel) red = red / (100*(high-low)) @@ -331,11 +321,8 @@ class renderedImageZoom: else: self.image = skydome.renderFromCamera(self.coords,self.betaM,self.g,self.altitude,self.temperature,self.pressure,self.sunDir) self.img = Image.fromarray(self.image,mode="RGB") - #self.tk_image = ImageTk.PhotoImage(width=256,height=256,image=self.img) - #img = ImageTk.PhotoImage(Image.open("highpress_camera.png")) - - #self.canvas.create_image(0,0, anchor="nw", image=self.tk_image) + #scaling function so every rendered image is the same scale of zoomed in as the previous one self.img = self.img.resize((int(self.img.width * self.zoom_factor), int(self.img.height * self.zoom_factor))) self.tk_image = ImageTk.PhotoImage(self.img) self.canvas.delete("all") @@ -365,11 +352,7 @@ class renderedImageZoom: -# def show_image(self,coords): -# self.coords = coords - #image_shower = renderedImageZoom(image_window,coords) -# self.image_window.mainloop() - + #map window functions for adding markers def add_marker_event(self,coords): self.coords = coords print("marker added to {}".format(self.coords)) @@ -397,7 +380,6 @@ class renderedImageZoom: image_window = tkinter.Toplevel() -#image_shower = renderedImageZoom(image_window,coords) image_shower = renderedImageZoom(root_tk,image_window,environment_tk) root_tk.mainloop() diff --git a/astropy_test.py b/astropy_test.py deleted file mode 100644 index 60e56ba..0000000 --- a/astropy_test.py +++ /dev/null @@ -1,21 +0,0 @@ -from astropy.time import Time -from astropy.coordinates import solar_system_ephemeris, EarthLocation , AltAz -from astropy.coordinates import get_body_barycentric, get_body - -def get_location(long,lat): - return EarthLocation.from_geodetic(long,lat) - -loc = get_location(45.508,-73.561) - -t = Time("2024-10-07 18:16",location=loc) - -with solar_system_ephemeris.set('builtin'): - sun = get_body('sun', t, loc) - moon = get_body('moon', t, loc) - saturn = get_body('saturn', t, loc) - -new_sun = sun.transform_to(AltAz(obstime=t,location=loc)) -new_moon = moon.transform_to(AltAz(obstime=t,location=loc)) -new_saturn = saturn.transform_to(AltAz(obstime=t,location=loc)) - -print(saturn) diff --git a/skydome.py b/skydome.py index 93463ff..b0600d7 100644 --- a/skydome.py +++ b/skydome.py @@ -18,8 +18,10 @@ Not all of the functions are used yet - but they might be later when we start ad """ +#typed Vec3 class with appropriate functions, so we dont have to mess with ndarray[3] all the time @cython.cclass class Vec3: + #cython specific declaration of the init and the v vactor 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] @@ -36,9 +38,7 @@ class Vec3: 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) @@ -76,8 +76,6 @@ wavelengths = Vec3(0.68, 0.55, 0.44) # These are the wavelengths in um @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 = T # degrees Celsius - #p_s:cython.int = P # mmHg t_s:cython.int = 25 # Temperature in degrees Celsius p_s:cython.int = 760 # Pressure in mmHg @@ -215,8 +213,9 @@ def computeIncidentLight(direction:Vec3,betaM,betaR,g,observerEarthRadius,sunDir i:cython.int for i in range(numSamples): + #first ray in the ray walking algorithm samplePosition = orig + direction * (tCurrent + segmentLength * 0.5) - height = samplePosition.length() - observerEarthRadius + height = samplePosition.length() - earthRadius hr = math.exp(-height / Hr) * segmentLength hm = math.exp(-height / Hm) * segmentLength opticalDepthR += hr @@ -235,8 +234,9 @@ def computeIncidentLight(direction:Vec3,betaM,betaR,g,observerEarthRadius,sunDir j:cython.int for j in range(numSamplesLight): + #second ray in the ray walking samplePositionLight = samplePosition + (sunDir * (tCurrentLight + segmentLengthLight * 0.5)) - heightLight = samplePositionLight.length() - observerEarthRadius + heightLight = samplePositionLight.length() - earthRadius if heightLight < 0: break opticalDepthLightR += (math.exp(-heightLight / Hr) * segmentLengthLight) @@ -246,6 +246,7 @@ def computeIncidentLight(direction:Vec3,betaM,betaR,g,observerEarthRadius,sunDir 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])) + #separate contributions of Rayleigh and Mie contributions, with attentuation dependant on the BetaM/R, and altitude and pressure (through tau) sumR += (attenuation * hr) sumM += (attenuation * hm) tCurrent += (segmentLength) @@ -256,6 +257,7 @@ def computeIncidentLight(direction:Vec3,betaM,betaR,g,observerEarthRadius,sunDir final_color = (sumR * betaR * phaseR + sumM * betaM * phaseM) * 20 return final_color.to_tuple() +#intersection of a ray and a sphere @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 @@ -267,10 +269,9 @@ def raySphereIntersect(orig:Vec3, direction:Vec3, radius:cython.double) -> tuple t0, t1 = solveQuadratic(A, B, C) if t1 < 0: return (kInfinity, kInfinity) - #if t0 > t1: - # t0, t1 = t1, t0 return (t0, t1) +#quadratic equation solved. @cython.cfunc def solveQuadratic(a:cython.float, b:cython.float, c:cython.float) -> tuple[cython.float,cython.float]: if b == 0: @@ -337,10 +338,8 @@ def renderSkydome(filename,betaM,g,altitude,T,P,sunDir): 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 return image.astype(np.uint8) - #imageio.imwrite(filename, image.astype(np.uint8)) @@ -383,7 +382,6 @@ def renderFromCamera(filename,betaM,g,altitude,T,P,sunDir): #This changes for each pixel, it is the direction we are looking in direction = Vec3(rayx, rayy, -1).normalize() color = computeIncidentLight(direction,betaM,betaR,g,observerEarthRadius,sunDir) - #image[y, x] = np.array(color) image[y][x] = np.clip(color, 0, 1) elapsed_time = time.time() - start_time @@ -391,12 +389,6 @@ def renderFromCamera(filename,betaM,g,altitude,T,P,sunDir): 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") -