This commit is contained in:
Yingjie Wang 2023-11-10 11:29:28 -05:00
parent a98e40d16b
commit c2228a8019
3 changed files with 85 additions and 41 deletions

BIN
image.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 20 KiB

View File

@ -2,8 +2,9 @@ import numpy as np
from raytrace_2D import Disk from raytrace_2D import Disk
import csv import csv
import colorspace import colorspace
from rich.progress import track from rich.progress import Progress, track
from PIL import Image from PIL import Image
from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor, as_completed
class Ray: class Ray:
def __init__(self, origin, direction, intensity): def __init__(self, origin, direction, intensity):
@ -13,8 +14,8 @@ class Ray:
center = [0,0] center = [0,0]
r = 1 r = 1
n = 1.3 # n = 1.3
disk = Disk(center, r, n) # disk = Disk(center, r, n)
dx = 0.00001 dx = 0.00001
N = int(2*r /dx - 1) N = int(2*r /dx - 1)
min_intensity = 0.00001 min_intensity = 0.00001
@ -105,16 +106,7 @@ def water_refraction_index(t, wavelength):
def sun_spectral(wavelength): def sun_spectral(wavelength):
return 1e16/(np.power(wavelength,5)*(np.exp(6.62607015e6*2.99792458/(wavelength*1.380649*(5250+273.15)))-1)) return 1e16/(np.power(wavelength,5)*(np.exp(6.62607015e6*2.99792458/(wavelength*1.380649*(5250+273.15)))-1))
def rainbow(n_theta, temp): def modified_trace(wavelength, temp, center, r, dx, N, n_theta, d_theta, min_intensity, max_ray, max_angle):
angles, d_theta = np.linspace(0,np.pi/2,n_theta, retstep=True)
angle_XYZ = np.zeros((n_theta, 3))
angle_sRGB = np.zeros((n_theta, 3))
num_wavelength = len(colorspace.data_wavelength)
def modified_trace(wavelength_index):
wavelength = colorspace.data_wavelength[wavelength_index]
n = water_refraction_index(temp, wavelength) n = water_refraction_index(temp, wavelength)
disk = Disk(center, r, n) disk = Disk(center, r, n)
stack = [] stack = []
@ -122,6 +114,7 @@ def rainbow(n_theta, temp):
stack.append(Ray([x, 2*r], [0,-1], 10*np.abs(x))) stack.append(Ray([x, 2*r], [0,-1], 10*np.abs(x)))
ray_count = 0 ray_count = 0
XYZ = colorspace.wavelength2XYZ(wavelength)*sun_spectral(wavelength) XYZ = colorspace.wavelength2XYZ(wavelength)*sun_spectral(wavelength)
own_angle_XYZ = np.zeros((n_theta, 3))
while stack and ray_count < max_ray: while stack and ray_count < max_ray:
ray = stack.pop() ray = stack.pop()
ray_count += 1 ray_count += 1
@ -140,9 +133,30 @@ def rainbow(n_theta, temp):
points.append(np.concatenate((ray.origin, ray.origin+ray.direction,[ray.intensity]))) points.append(np.concatenate((ray.origin, ray.origin+ray.direction,[ray.intensity])))
if direction[1] > 0: if direction[1] > 0:
angle = np.arccos(direction[1]) angle = np.arccos(direction[1])
angle_XYZ[int(angle/d_theta)] += XYZ*ray.intensity*direction[1] if angle < max_angle:
own_angle_XYZ[int(angle/d_theta)] += XYZ*ray.intensity*direction[1]
return own_angle_XYZ
map(modified_trace, range(num_wavelength)) def rainbow(n_theta, max_theta, temp):
angles, d_theta = np.linspace(0,max_theta,n_theta, retstep=True)
angle_XYZ = np.zeros((n_theta, 3))
angle_sRGB = np.zeros((n_theta, 3))
num_wavelength = len(colorspace.data_wavelength)
with Progress() as progress:
# 创建一个进度条
task = progress.add_task("[green]Simulating...", total=num_wavelength)
# 使用 ThreadPoolExecutor 并行执行
with ProcessPoolExecutor() as executor:
futures = [executor.submit(modified_trace, wavelength, temp, center, r, dx, N, n_theta, d_theta, min_intensity, max_ray, max_theta) for wavelength in colorspace.data_wavelength]
for future in as_completed(futures):
result = future.result()
angle_XYZ += result
progress.update(task, advance=1)
for i,XYZ in enumerate(angle_XYZ): for i,XYZ in enumerate(angle_XYZ):
angle_sRGB[i] = colorspace.XYZ2RGB(XYZ) angle_sRGB[i] = colorspace.XYZ2RGB(XYZ)
@ -153,12 +167,21 @@ def rainbow(n_theta, temp):
np.clip(angle_sRGB, 0, 1) np.clip(angle_sRGB, 0, 1)
return angles,angle_sRGB return angles,angle_sRGB
def bin_find(x, list): def bin_find(x, list, start):
l,r = 0, len(list)-1 l,r = 0, len(list)-1
if x < list[l]: if x < list[l]:
return 0 return 0
if x >= list[r]: if x >= list[r]:
return r return r
if list[start] <= x < list[start + 1]:
return start
elif list[start] < x:
l = start
else:
r = start
while l <= r: while l <= r:
mid = l + (r-l)//2 mid = l + (r-l)//2
if list[mid] <= x < list[mid + 1]: if list[mid] <= x < list[mid + 1]:
@ -169,18 +192,39 @@ def bin_find(x, list):
r = mid - 1 r = mid - 1
return -1 return -1
def draw_column(i, w, h, radius, angle_sRGB):
index = 0
column = np.zeros((h, 3), dtype=np.int8)
for j in range(h):
r = np.linalg.norm(np.array([i,j]) - [w/2, 0])
index = bin_find(r, radius, index)
column[j] = angle_sRGB[index]*255
return i,column
def take_picture(angles, angle_sRGB, w, h, distance, filename="image.png"): def take_picture(angles, angle_sRGB, w, h, distance, filename="image.png"):
image_array = np.zeros((w, h, 3), dtype=np.uint8) image_array = np.zeros((w, h, 3), dtype=np.uint8)
radius = distance * np.tan(angles) radius = distance * np.tan(angles)
for i in track(range(w), description="generating picture..."):
for j in range(h): with Progress() as progress:
r = np.linalg.norm(np.array([i,j]) - [w/2, 0]) # 创建一个进度条
index = bin_find(r, radius) task2 = progress.add_task("[green]Rendering...", total=w)
image_array[i,j] = angle_sRGB[index]*255
# 使用 ThreadPoolExecutor 并行执行
with ProcessPoolExecutor() as executor:
futures2 = [executor.submit(draw_column, i, w,h,radius,angle_sRGB) for i in range(w)]
for future in as_completed(futures2):
ind,column = future.result()
image_array[ind] = column
progress.update(task2, advance=1)
image = Image.fromarray(image_array, "RGB") image = Image.fromarray(image_array, "RGB")
image.save(filename) image.save(filename)
if __name__=="__main__": if __name__=="__main__":
angles,sRGB=rainbow(10000, 10) w = 7680
h = 4320
dis = 2400
max_angle = 1.1*np.arctan(np.sqrt(1+(h*h+w*w/4)/(dis*dis)))
angles,sRGB=rainbow(100, max_angle, 10)
np.savez_compressed("saved.npz", a=angles, b=sRGB) np.savez_compressed("saved.npz", a=angles, b=sRGB)
take_picture(angles,sRGB, 7680, 4320, 2400, "image.png") take_picture(angles,sRGB, 7680, 4320, 2400, "image.png")

BIN
test.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.3 MiB