technoduck/sunset_simulator
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

fixed altitude

Browse source
This commit is contained in:
Technoduck 2024-12-11 17:22:35 -05:00
parent e2069080f6
commit d5768a4bfe
4 changed files with 26 additions and 518 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

445
Midterm_submittion.py
View file

@ -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[ ]:

52
app.py
View file

@ -16,16 +16,11 @@ def get_api_params(coords):
# create tkinter window # create tkinter window
root_tk = tkinter.Tk() 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() environment_tk = tkinter.Tk()
#Prop for a colo changing button
class EnvButton(tkinter.Button): class EnvButton(tkinter.Button):
def __init__(self,*args,**kwargs): def __init__(self,*args,**kwargs):
tkinter.Button.__init__(self, *args, **kwargs) tkinter.Button.__init__(self, *args, **kwargs)
@ -40,9 +35,10 @@ class EnvButton(tkinter.Button):
self['bg'] = self.default_bg_color self['bg'] = self.default_bg_color
# main window class of the rendered Image window
class renderedImageZoom: class renderedImageZoom:
def __init__(self,root,image_window,aerosol_window): def __init__(self,root,image_window,aerosol_window):
#setting of the default parameters
self.root = root self.root = root
self.root.geometry(f"{1000}x{700}") self.root.geometry(f"{1000}x{700}")
self.root.title("map_view_simple_example.py") self.root.title("map_view_simple_example.py")
@ -72,9 +68,10 @@ class renderedImageZoom:
self.image_window.title("Simulated Sunset") self.image_window.title("Simulated Sunset")
self.image_window.config(width=256,height=256) 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 = tkinter.Canvas(self.image_window,bg="white")
self.canvas.pack(fill=tkinter.BOTH,expand=True) self.canvas.pack(fill=tkinter.BOTH,expand=True)
#separate function for ease of creating the environment options window
self.environment_window_filler() self.environment_window_filler()
@ -102,9 +99,7 @@ class renderedImageZoom:
# #
#Antarctic - 11e-6 g = 0.784 #Antarctic - 11e-6 g = 0.784
#code for changing the color of the active button
# 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]
def set_active(index): def set_active(index):
if index != self.curr_env: if index != self.curr_env:
@ -118,7 +113,7 @@ class renderedImageZoom:
buton_cont_clean['bg'] = 'cadetblue1' buton_cont_clean['bg'] = 'cadetblue1'
self.env_setter(26e-6,0.709) 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 = 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) buton_cont_avr.grid(column=0,row=2,padx=10,pady=5)
@ -154,12 +149,11 @@ class renderedImageZoom:
self.set_temp(value) self.set_temp(value)
label.config(text=str(sorted([-270,int(value),100]))) 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) t_frame = tkinter.Frame(self.aerosol_window)
text_t = tkinter.Text(t_frame,height=1,width=10) text_t = tkinter.Text(t_frame,height=1,width=10)
submit_t_button = tkinter.Button(t_frame,text="set temperature", 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")]) 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") t_label = tkinter.Label(t_frame,text="5 C")
text_t.pack(side=tkinter.LEFT) text_t.pack(side=tkinter.LEFT)
t_label.pack(side=tkinter.RIGHT) t_label.pack(side=tkinter.RIGHT)
@ -189,7 +183,6 @@ class renderedImageZoom:
s_frame.grid(column=0,row=9,columnspan=2,padx=20) s_frame.grid(column=0,row=9,columnspan=2,padx=20)
## Altitude
a_frame = tkinter.Frame(self.aerosol_window) a_frame = tkinter.Frame(self.aerosol_window)
text_a = tkinter.Text(a_frame,height=1,width=10) 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")))) 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) r_frame.grid(column=0,row=12,columnspan=2,padx=20)
#red and blue average calculations
b_frame = tkinter.Frame(self.aerosol_window) b_frame = tkinter.Frame(self.aerosol_window)
text_b_l = tkinter.Text(b_frame,height=1,width=10) text_b_l = tkinter.Text(b_frame,height=1,width=10)
text_b_u = 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) submit_b_button.pack(side=tkinter.RIGHT)
b_frame.grid(column=0,row=13,columnspan=2,padx=20) 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(): def cam_toggle():
if cam_fish_toggle.config('text')[-1] == "cam view": if cam_fish_toggle.config('text')[-1] == "cam view":
cam_fish_toggle.config(text="fisheye") cam_fish_toggle.config(text="fisheye")
@ -242,7 +232,7 @@ class renderedImageZoom:
self.fisheye = False 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_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_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) render_button = tkinter.Button(self.image_window,text="render",command=self.render)
@ -257,6 +247,7 @@ class renderedImageZoom:
self.canvas.bind("<Button-4>",self.zoom_in) self.canvas.bind("<Button-4>",self.zoom_in)
self.canvas.bind("<Button-5>", self.zoom_out) self.canvas.bind("<Button-5>", self.zoom_out)
#setter functions
def set_temp(self,temp): def set_temp(self,temp):
if temp is None: if temp is None:
self.temperature = 0 self.temperature = 0
@ -275,15 +266,16 @@ class renderedImageZoom:
if angle is None: if angle is None:
angle_degrees:cython.int = 90 # Sun direction in degrees (0 to 90, noon to sunset) angle_degrees:cython.int = 90 # Sun direction in degrees (0 to 90, noon to sunset)
elif int(angle) > 90 or int(angle) < 0: 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: 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 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)) sunDir:skydome.Vec3 = skydome.Vec3(0, math.cos(angle_radians), -math.sin(angle_radians))
self.sunDir = sunDir self.sunDir = sunDir
#setters have reasonable (ish) limiters on the values
def set_press(self,press): def set_press(self,press):
if press is None: if press is None:
self.pressure = 720 self.pressure = 720
@ -310,10 +302,8 @@ class renderedImageZoom:
for line in slice: for line in slice:
for pixel in line: for pixel in line:
if sum(pixel) == 0: if sum(pixel) == 0:
print(pixel)
red += 0 red += 0
else: else:
print(pixel)
red += pixel[0]/sum(pixel) red += pixel[0]/sum(pixel)
red = red / (100*(high-low)) red = red / (100*(high-low))
@ -331,11 +321,8 @@ class renderedImageZoom:
else: else:
self.image = skydome.renderFromCamera(self.coords,self.betaM,self.g,self.altitude,self.temperature,self.pressure,self.sunDir) 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.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.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.tk_image = ImageTk.PhotoImage(self.img)
self.canvas.delete("all") self.canvas.delete("all")
@ -365,11 +352,7 @@ class renderedImageZoom:
# def show_image(self,coords): #map window functions for adding markers
# self.coords = coords
#image_shower = renderedImageZoom(image_window,coords)
# self.image_window.mainloop()
def add_marker_event(self,coords): def add_marker_event(self,coords):
self.coords = coords self.coords = coords
print("marker added to {}".format(self.coords)) print("marker added to {}".format(self.coords))
@ -397,7 +380,6 @@ class renderedImageZoom:
image_window = tkinter.Toplevel() image_window = tkinter.Toplevel()
#image_shower = renderedImageZoom(image_window,coords)
image_shower = renderedImageZoom(root_tk,image_window,environment_tk) image_shower = renderedImageZoom(root_tk,image_window,environment_tk)
root_tk.mainloop() root_tk.mainloop()

21
astropy_test.py
View file

@ -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)

26
skydome.py
View file

@ -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 @cython.cclass
class Vec3: class Vec3:
#cython specific declaration of the init and the v vactor
v = cython.declare(cython.float[3],visibility="public") v = cython.declare(cython.float[3],visibility="public")
def __cinit__(self, x:cython.float =0, y:cython.float=0, z:cython.float=0): def __cinit__(self, x:cython.float =0, y:cython.float=0, z:cython.float=0):
self.v:cython.float[3] = [x, y, z] self.v:cython.float[3] = [x, y, z]
@ -36,9 +38,7 @@ class Vec3:
def __mul__(self, other): def __mul__(self, other):
if isinstance(other, Vec3): if isinstance(other, Vec3):
return Vec3(self.v[0] * other.v[0],self.v[1] * other.v[1],self.v[2] * other.v[2]) 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) 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): def __truediv__(self, scalar):
return Vec3(self.v[0] / scalar,self.v[1] / scalar,self.v[2] / 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 @cython.cfunc
def refraction_calculator(T:cython.int, P:cython.int) -> tuple[Vec3, Vec3]: def refraction_calculator(T:cython.int, P:cython.int) -> tuple[Vec3, Vec3]:
alpha:cython.float = 0.00367 # in inverse Celsius 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 t_s:cython.int = 25 # Temperature in degrees Celsius
p_s:cython.int = 760 # Pressure in mmHg p_s:cython.int = 760 # Pressure in mmHg
@ -215,8 +213,9 @@ def computeIncidentLight(direction:Vec3,betaM,betaR,g,observerEarthRadius,sunDir
i:cython.int i:cython.int
for i in range(numSamples): for i in range(numSamples):
#first ray in the ray walking algorithm
samplePosition = orig + direction * (tCurrent + segmentLength * 0.5) samplePosition = orig + direction * (tCurrent + segmentLength * 0.5)
height = samplePosition.length() - observerEarthRadius height = samplePosition.length() - earthRadius
hr = math.exp(-height / Hr) * segmentLength hr = math.exp(-height / Hr) * segmentLength
hm = math.exp(-height / Hm) * segmentLength hm = math.exp(-height / Hm) * segmentLength
opticalDepthR += hr opticalDepthR += hr
@ -235,8 +234,9 @@ def computeIncidentLight(direction:Vec3,betaM,betaR,g,observerEarthRadius,sunDir
j:cython.int j:cython.int
for j in range(numSamplesLight): for j in range(numSamplesLight):
#second ray in the ray walking
samplePositionLight = samplePosition + (sunDir * (tCurrentLight + segmentLengthLight * 0.5)) samplePositionLight = samplePosition + (sunDir * (tCurrentLight + segmentLengthLight * 0.5))
heightLight = samplePositionLight.length() - observerEarthRadius heightLight = samplePositionLight.length() - earthRadius
if heightLight < 0: if heightLight < 0:
break break
opticalDepthLightR += (math.exp(-heightLight / Hr) * segmentLengthLight) opticalDepthLightR += (math.exp(-heightLight / Hr) * segmentLengthLight)
@ -246,6 +246,7 @@ def computeIncidentLight(direction:Vec3,betaM,betaR,g,observerEarthRadius,sunDir
if j == numSamplesLight - 1: if j == numSamplesLight - 1:
tau = (betaR * (opticalDepthR + opticalDepthLightR)) + (betaM * (opticalDepthM + opticalDepthLightM)) tau = (betaR * (opticalDepthR + opticalDepthLightR)) + (betaM * (opticalDepthM + opticalDepthLightM))
attenuation = Vec3(math.exp(-tau.v[0]), math.exp(-tau.v[1]), math.exp(-tau.v[2])) 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) sumR += (attenuation * hr)
sumM += (attenuation * hm) sumM += (attenuation * hm)
tCurrent += (segmentLength) tCurrent += (segmentLength)
@ -256,6 +257,7 @@ def computeIncidentLight(direction:Vec3,betaM,betaR,g,observerEarthRadius,sunDir
final_color = (sumR * betaR * phaseR + sumM * betaM * phaseM) * 20 final_color = (sumR * betaR * phaseR + sumM * betaM * phaseM) * 20
return final_color.to_tuple() return final_color.to_tuple()
#intersection of a ray and a sphere
@cython.cfunc @cython.cfunc
def raySphereIntersect(orig:Vec3, direction:Vec3, radius:cython.double) -> tuple[cython.float,cython.float]: 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 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) t0, t1 = solveQuadratic(A, B, C)
if t1 < 0: if t1 < 0:
return (kInfinity, kInfinity) return (kInfinity, kInfinity)
#if t0 > t1:
# t0, t1 = t1, t0
return (t0, t1) return (t0, t1)
#quadratic equation solved.
@cython.cfunc @cython.cfunc
def solveQuadratic(a:cython.float, b:cython.float, c:cython.float) -> tuple[cython.float,cython.float]: def solveQuadratic(a:cython.float, b:cython.float, c:cython.float) -> tuple[cython.float,cython.float]:
if b == 0: 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}, elapsed time: {elapsed_time:.2f} seconds")
#print(f"Rendering row {j + 1}/{height}") #print(f"Rendering row {j + 1}/{height}")
# Save result to a PNG image
image = np.clip(image, 0, 1) * 255 # change 255 image = np.clip(image, 0, 1) * 255 # change 255
return image.astype(np.uint8) 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 #This changes for each pixel, it is the direction we are looking in
direction = Vec3(rayx, rayy, -1).normalize() direction = Vec3(rayx, rayy, -1).normalize()
color = computeIncidentLight(direction,betaM,betaR,g,observerEarthRadius,sunDir) color = computeIncidentLight(direction,betaM,betaR,g,observerEarthRadius,sunDir)
#image[y, x] = np.array(color)
image[y][x] = np.clip(color, 0, 1) image[y][x] = np.clip(color, 0, 1)
elapsed_time = time.time() - start_time 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 image = np.clip(image, 0, 1) * 255
return image.astype(np.uint8) return image.astype(np.uint8)
#imageio.imwrite(filename, image.astype(np.uint8))
#renderFromCamera("camera_render.png")
#renderSkydome("highpress_15C.png")
#renderFromCamera("highpress_camera.png")