diff --git a/Midterm_submittion.py b/Midterm_submittion.py new file mode 100644 index 0000000..dfbf88d --- /dev/null +++ b/Midterm_submittion.py @@ -0,0 +1,445 @@ +#!/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 4fae26b..8446f19 100644 --- a/app.py +++ b/app.py @@ -13,89 +13,181 @@ def get_api_params(coords): # create tkinter window root_tk = tkinter.Tk() -root_tk.geometry(f"{1000}x{700}") -root_tk.title("map_view_simple_example.py") - -# create map widget -map_widget = tkintermapview.TkinterMapView(root_tk, width=1000, height=700, corner_radius=0) -map_widget.pack(fill="both", expand=True) # 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_tk = tkinter.Tk() + class renderedImageZoom: - def __init__(self,root,coords): - self.image_window = root + def __init__(self,root,image_window,aerosol_window): + self.root = root + self.root.geometry(f"{1000}x{700}") + self.root.title("map_view_simple_example.py") + + self.aerosol_window = aerosol_window +# create map widget + self.map_widget = tkintermapview.TkinterMapView(self.root, width=1000, height=700, corner_radius=0) + self.map_widget.pack(fill="both", expand=True) + self.map_widget.add_right_click_menu_command(label="Set as observer location", + command=self.add_marker_event, + pass_coords=True) + + self.map_widget.add_right_click_menu_command(label="get altitude", + command=self.get_altitude, + pass_coords=True) + + #self.map_widget.add_right_click_menu_command(label="open_window", + # command=self.show_image, + # pass_coords=True) + + + self.aerosol_window.title("options selector") + + + + environemnt_info_text = tkinter.Label(self.aerosol_window, text="Please choose in which environment the sunset is occuring") + environemnt_info_text.pack(side='top',pady=20,padx=10) + + +#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 + + buton_cont_clean = tkinter.Button(self.aerosol_window,text="Continental Clean",command=lambda: self.env_setter(26e-6,0.709)); + buton_cont_clean.pack(side='top',padx=10,pady=5) + + buton_cont_avr = tkinter.Button(self.aerosol_window,text="Continental Average ",command=lambda: self.env_setter(75e-6,0.793)); + buton_cont_avr.pack(side='top',padx=10,pady=5) + + buton_cont_poll = tkinter.Button(self.aerosol_window,text="Continental Polluted ",command=lambda: self.env_setter(175e-6,0.698)); + buton_cont_poll.pack(side='top',padx=10,pady=5) + + buton_urban = tkinter.Button(self.aerosol_window,text="Urban ",command=lambda: self.env_setter(353e-6,0.689)); + buton_urban.pack(side='top',padx=10,pady=5) + + buton_desert = tkinter.Button(self.aerosol_window,text="Desert ",command=lambda: self.env_setter(145e-6,0.729)); + buton_desert.pack(side='top',padx=10,pady=5) + + buton_mar_clean = tkinter.Button(self.aerosol_window,text="Maritime Clean ",command=lambda: self.env_setter(90e-6,0.772)); + buton_mar_clean.pack(side='top',padx=10,pady=5) + + buton_mar_poll = tkinter.Button(self.aerosol_window,text="Maritime Polluted ",command=lambda: self.env_setter(115e-6,0.756)); + buton_mar_poll.pack(side='top',padx=10,pady=5) + + buton_mar_tro = tkinter.Button(self.aerosol_window,text="Maritime Tropic ",command=lambda: self.env_setter(43e-6,0.774)); + buton_mar_tro.pack(side='top',padx=10,pady=5) + + buton_arctic = tkinter.Button(self.aerosol_window,text="Arctic ",command=lambda: self.env_setter(23e-6,0.721)); + buton_arctic.pack(side='top',padx=10,pady=5) + + buton_antarctic = tkinter.Button(self.aerosol_window,text="Antarctic ",command=lambda: self.env_setter(11e-6,0.784)); + buton_antarctic.pack(side='top',padx=10,pady=5) + + + self.image_window = image_window self.image_window.title("Simulated Sunset") self.image_window.config(width=256,height=256) - image = skydome.renderFromCamera(coords) + #image = skydome.renderFromCamera(coords) self.canvas = tkinter.Canvas(self.image_window,bg="white") self.canvas.pack(fill=tkinter.BOTH,expand=True) + #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")) + + #self.canvas.create_image(0,0, anchor="center", image=self.tk_image) + 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) + zoom_in_button.pack(side=tkinter.LEFT) + zoom_out_button.pack(side=tkinter.LEFT) + render_button.pack(side=tkinter.RIGHT) + + self.canvas.bind("",self.zoom_in) + self.canvas.bind("", self.zoom_out) + + def env_setter(self,environment_constant,g): + self.betaM = skydome.Vec3(environment_constant,environment_constant,environment_constant) + self.g = g + + def render(self): + image = skydome.renderFromCamera(self.coords,self.betaM,self.g,self.altitude) 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")) self.canvas.create_image(0,0, anchor="nw", image=self.tk_image) - 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) - zoom_in_button.pack(side=tkinter.LEFT) - zoom_out_button.pack(side=tkinter.LEFT) - self.canvas.bind("",self.zoom_in) - self.canvas.bind("", self.zoom_out) + def zoom_in(self, event=None): # Increase the image size by a factor (e.g., 1.2) self.img = self.img.resize((int(self.img.width * 1.2), int(self.img.height * 1.2))) self.tk_image = ImageTk.PhotoImage(self.img) self.canvas.delete("all") - self.canvas.create_image(0, 0, anchor=tkinter.NW, image=self.tk_image) + self.canvas.create_image(0, 0, anchor="nw", image=self.tk_image) def zoom_out(self, event=None): # Decrease the image size by a factor (e.g., 0.8) self.img = self.img.resize((int(self.img.width * 0.8), int(self.img.height * 0.8)))#,resampling) self.tk_image = ImageTk.PhotoImage(self.img) self.canvas.delete("all") - self.canvas.create_image(0, 0, anchor=tkinter.NW, image=self.tk_image) + self.canvas.create_image(0, 0, anchor="nw", image=self.tk_image) -def show_image(coords): - image_window = tkinter.Toplevel() +# def show_image(self,coords): +# self.coords = coords + #image_shower = renderedImageZoom(image_window,coords) +# self.image_window.mainloop() - image_shower = renderedImageZoom(image_window,coords) - image_window.mainloop() + def add_marker_event(self,coords): + self.coords = coords + print("marker added to {}".format(self.coords)) + self.map_widget.delete_all_marker() + new_marker = self.map_widget.set_marker(self.coords[0], self.coords[1], text="sunset observer") + self.get_altitude(coords) -def marker_click(marker): - print(f"marker clicked - text: {marker.text} position: {marker.position}") + def marker_click(self,marker): + print(f"marker clicked - text: {marker.text} position: {marker.position}") -def add_marker_event(coords): - print("marker added to {}".format(coords)) - map_widget.delete_all_marker() - new_marker = map_widget.set_marker(coords[0], coords[1], text="sunset observer") -def get_altitude(coords): - print("getting altitude for {}\n".format(coords)) - params = get_api_params(coords) - r = requests.get(URL,params) - data = r.json() - print(data['results'][0]['elevation']) + def get_altitude(self,coords): + self.coords = coords + print("getting altitude for {}\n".format(self.coords)) + params = get_api_params(self.coords) + r = requests.get(URL,params) + data = r.json() + print(data['results'][0]['elevation']) + self.altitude = data['results'][0]['elevation'] -map_widget.add_right_click_menu_command(label="Set as observer location", - command=add_marker_event, - pass_coords=True) -map_widget.add_right_click_menu_command(label="get altitude", - command=get_altitude, - pass_coords=True) -map_widget.add_right_click_menu_command(label="open_window", - command=show_image, - pass_coords=True) + +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/skydome.py b/skydome.py index 471c721..e665ef8 100644 --- a/skydome.py +++ b/skydome.py @@ -222,7 +222,7 @@ 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 +#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 @@ -235,12 +235,12 @@ Hm:cython.int = 1200 # The direction will change for each pixel @cython.cfunc -def computeIncidentLight(direction:Vec3) -> tuple[float,float,float]: +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, earthRadius + 1, 0) + orig:Vec3 = Vec3(0, observerEarthRadius + 1, 0) t0:cython.float t1:cython.float t0, t1 = raySphereIntersect(orig, direction, atmosphereRadius) @@ -262,13 +262,13 @@ def computeIncidentLight(direction:Vec3) -> tuple[float,float,float]: 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 + #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() - earthRadius + height = samplePosition.length() - observerEarthRadius hr = math.exp(-height / Hr) * segmentLength hm = math.exp(-height / Hm) * segmentLength opticalDepthR += hr @@ -288,7 +288,7 @@ def computeIncidentLight(direction:Vec3) -> tuple[float,float,float]: j:cython.int for j in range(numSamplesLight): samplePositionLight = samplePosition + (sunDir * (tCurrentLight + segmentLengthLight * 0.5)) - heightLight = samplePositionLight.length() - earthRadius + heightLight = samplePositionLight.length() - observerEarthRadius if heightLight < 0: break opticalDepthLightR += (math.exp(-heightLight / Hr) * segmentLengthLight) @@ -389,8 +389,9 @@ def renderSkydome(filename): -def renderFromCamera(filename): +def renderFromCamera(filename,betaM,g,altitude): + observerEarthRadius = earthRadius + altitude width:cython.int height:cython.int width, height = 100, 100 @@ -420,7 +421,7 @@ def renderFromCamera(filename): #This changes for each pixel, it is the direction we are looking in direction = Vec3(rayx, rayy, -1).normalize() - color = computeIncidentLight(direction) + color = computeIncidentLight(direction,betaM,g,observerEarthRadius) #image[y, x] = np.array(color) image[y][x] = np.clip(color, 0, 1)