From d9b8013d2e3ceb57276270f9b299278520ebdee8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jakub=20Po=C4=87wiardowski?= <80002864+jpocwiar@users.noreply.github.com> Date: Sun, 19 Jan 2025 14:26:03 +0100 Subject: [PATCH] Console for photon mapping in python (#5) * PM Python can be loaded w params * Update README.md --- README.md | 5 + code/photonmapping/main.py | 221 ------------------------------------- config.ini | 6 +- main.py | 19 +++- outputs/output.png | Bin 20619 -> 5556 bytes photon_mapping.py | 191 ++++++++++++++++++++++++++++++++ 6 files changed, 217 insertions(+), 225 deletions(-) delete mode 100644 code/photonmapping/main.py create mode 100644 photon_mapping.py diff --git a/README.md b/README.md index a011c8c8..fb3e03a2 100644 --- a/README.md +++ b/README.md @@ -43,3 +43,8 @@ python main.py --algorithm ray_tracing # Wywołanie algorytmu ray tracing ze specyfikacją sceny z folderu scenes, liczbą sampli na pixek, rozdzielczością, środowiskiem i rozmyciem środowiska python main.py --scene three_spheres --samples_per_pixel 100 --resolution 100x100 --environment lake.png --env_blur 10 ``` + +```bash +# Wywołanie algorytmu photon mapping ze specyfikacją liczby fotonów i maksymalnej głębokości +python main.py --algorithm photon_mapping --max_depth --num_photons 1000 +``` diff --git a/code/photonmapping/main.py b/code/photonmapping/main.py deleted file mode 100644 index 9aec82d6..00000000 --- a/code/photonmapping/main.py +++ /dev/null @@ -1,221 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt - -# Define basic vector operations -class Vector3: - def __init__(self, x, y, z): - self.x = x - self.y = y - self.z = z - - def __add__(self, other): - return Vector3(self.x+other.x, self.y+other.y, self.z+other.z) - - def __sub__(self, other): - return Vector3(self.x-other.x, self.y-other.y, self.z-other.z) - - def __mul__(self, scalar): - return Vector3(self.x*scalar, self.y*scalar, self.z*scalar) - - def dot(self, other): - return self.x*other.x + self.y*other.y + self.z*other.z - - def norm(self): - return np.sqrt(self.dot(self)) - - def normalize(self): - n = self.norm() - return Vector3(self.x/n, self.y/n, self.z/n) - -# Define the photon -class Photon: - def __init__(self, position, direction, power): - self.position = position - self.direction = direction - self.power = power - -# Define a simple sphere -class Sphere: - def __init__(self, center, radius, color): - self.center = center - self.radius = radius - self.color = color - - def intersect(self, ray_origin, ray_direction): - # Solve quadratic equation for intersection - oc = ray_origin - self.center - a = ray_direction.dot(ray_direction) - b = 2.0 * oc.dot(ray_direction) - c = oc.dot(oc) - self.radius*self.radius - discriminant = b*b - 4*a*c - if discriminant < 0: - return None # No intersection - else: - t = (-b - np.sqrt(discriminant)) / (2.0*a) - if t < 0: - t = (-b + np.sqrt(discriminant)) / (2.0*a) - if t < 0: - return None - hit_point = ray_origin + ray_direction * t - normal = (hit_point - self.center).normalize() - return (t, hit_point, normal) - -# Define a simple plane -class Plane: - def __init__(self, point, normal, color): - self.point = point - self.normal = normal.normalize() - self.color = color - - def intersect(self, ray_origin, ray_direction): - denom = self.normal.dot(ray_direction) - if abs(denom) > 1e-6: - t = (self.point - ray_origin).dot(self.normal) / denom - if t >= 0: - hit_point = ray_origin + ray_direction * t - return (t, hit_point, self.normal) - return None - -# Scene setup -sphere = Sphere(Vector3(0, 0, -5), 1.0, np.array([1, 0, 0])) # Red sphere -plane = Plane(Vector3(0, -1, 0), Vector3(0, 1, 0), np.array([0.5, 0.5, 0.5])) # Gray plane - -objects = [sphere, plane] - -# Light source -light_position = Vector3(-5, 5, -5) -light_power = np.array([1, 1, 1]) * 1000 # Intense white light - -# Photon map -photon_map = [] - -# Parameters -num_photons = 10000 # Number of photons to emit -max_depth = 5 # Maximum number of bounces -gather_radius = 0.5 # Radius for radiance estimation - -def emit_photons(): - for _ in range(num_photons): - # Emit photons in random directions from the light source - direction = random_unit_vector() - power = light_power / num_photons - photon = Photon(light_position, direction, power) - trace_photon(photon, 0) - -def trace_photon(photon, depth): - if depth > max_depth: - return - closest_t = np.inf - hit_object = None - hit_info = None - # Find the nearest intersection - for obj in objects: - result = obj.intersect(photon.position, photon.direction) - if result: - t, hit_point, normal = result - if t < closest_t: - closest_t = t - hit_object = obj - hit_info = (hit_point, normal) - if hit_object: - hit_point, normal = hit_info - photon_map.append((hit_point, photon.power)) - # Diffuse reflection - new_direction = random_hemisphere_direction(normal) - photon.position = hit_point - photon.direction = new_direction - # Absorb some power - photon.power = photon.power * 0.8 # Simple absorption - trace_photon(photon, depth+1) - -def random_unit_vector(): - theta = np.random.uniform(0, 2*np.pi) - z = np.random.uniform(-1, 1) - r = np.sqrt(1 - z*z) - return Vector3(r * np.cos(theta), r * np.sin(theta), z) - -def random_hemisphere_direction(normal): - dir = random_unit_vector() - if dir.dot(normal) < 0: - dir = Vector3(-dir.x, -dir.y, -dir.z) - return dir - -def render_image(width, height): - aspect_ratio = width / height - fov = np.pi / 3 # 60 degrees field of view - image = np.zeros((height, width, 3)) - for y in range(height): - for x in range(width): - # Convert pixel coordinate to camera ray - px = (2 * (x + 0.5) / width - 1) * np.tan(fov / 2) * aspect_ratio - py = (1 - 2 * (y + 0.5) / height) * np.tan(fov / 2) - ray_origin = Vector3(0, 0, 0) - ray_direction = Vector3(px, py, -1).normalize() - color = trace_ray(ray_origin, ray_direction) - image[y, x, :] = np.clip(color, 0, 1) - return image - -def trace_ray(ray_origin, ray_direction): - closest_t = np.inf - hit_object = None - hit_info = None - # Find the nearest intersection - for obj in objects: - result = obj.intersect(ray_origin, ray_direction) - if result: - t, hit_point, normal = result - if t < closest_t: - closest_t = t - hit_object = obj - hit_info = (hit_point, normal, obj.color) - if hit_object: - hit_point, normal, color = hit_info - direct_light = compute_direct_light(hit_point, normal) - indirect_light = estimate_radiance(hit_point, normal) - return color * (direct_light + indirect_light) - else: - return np.array([0, 0, 0]) # Background color - -def compute_direct_light(point, normal): - # Simple Lambertian reflection from light source - direction_to_light = (light_position - point).normalize() - # Shadow ray - shadow_origin = point + normal * 1e-5 - shadow_ray = direction_to_light - in_shadow = False - for obj in objects: - result = obj.intersect(shadow_origin, shadow_ray) - if result: - in_shadow = True - break - if in_shadow: - return np.array([0, 0, 0]) - else: - intensity = max(0, normal.dot(direction_to_light)) - return intensity * light_power / (4 * np.pi * (light_position - point).norm()**2) - -def estimate_radiance(point, normal): - # Gather photons within the gather_radius - accumulated_power = np.array([0.0, 0.0, 0.0]) - for photon_pos, photon_power in photon_map: - distance = (photon_pos - point).norm() - if distance < gather_radius: - weight = max(0, normal.dot((photon_pos - point).normalize())) - accumulated_power += photon_power * weight - area = np.pi * gather_radius ** 2 - return accumulated_power / (area * num_photons) - -# Main execution -if __name__ == '__main__': - print("Emitting photons...") - emit_photons() - - print("Rendering image...") - width = 200 - height = 100 - image = render_image(width, height) - - # Display the image - plt.imshow(image) - plt.axis('off') - plt.show() diff --git a/config.ini b/config.ini index a688a058..987bfe7b 100644 --- a/config.ini +++ b/config.ini @@ -7,8 +7,10 @@ resolution = 400x300 output = output.png [ray_tracing] -max_depth = 5 ; Params for ray tracing +max_depth = 5 samples_per_pixel = 6 [photon_mapping] -photon_count = 100000 ; Params for photon mapping \ No newline at end of file +num_photons = 10000 +max_depth = 5 +gather_radius = 0.5 \ No newline at end of file diff --git a/main.py b/main.py index 9e68b148..53d8a8db 100644 --- a/main.py +++ b/main.py @@ -1,10 +1,10 @@ import argparse -from configparser import ConfigParser from rendering import ray_trace from utils import load_config, parse_resolution import importlib import os -# from scenes.cornell_box import * +import matplotlib.pyplot as plt +from photon_mapping import render_photon_mapping def main(): # default config @@ -22,6 +22,12 @@ def main(): parser.add_argument("--output", type=str, default=config.get('DEFAULT', 'output'), help="Output file name.") parser.add_argument('--num_spheres', type=int, default=3, help='Number of spheres in the scene for Ray Tracing 0') + parser.add_argument('--num_photons', type=int, default=config.getint('photon_mapping', 'num_photons'), help='Number of photons for photon mapping') + parser.add_argument('--max_depth', type=int, default=config.getint('photon_mapping', 'max_depth'), + help='Maximum depth for photon tracing') + parser.add_argument('--gather_radius', type=float, default=config.getfloat('photon_mapping', 'gather_radius'), + help='Radius for radiance estimation in photon mapping') + args = parser.parse_args() @@ -53,6 +59,15 @@ def main(): img.save(output_path) print(f"Image saved to {output_path}") img.show() + elif args.algorithm == "photon_mapping": + print("Starting photon mapping...") + image = render_photon_mapping(width, height, args.num_photons, args.max_depth, args.gather_radius) + plt.imshow(image) + plt.axis('off') + output_path = os.path.join("outputs", args.output) + plt.savefig(output_path) + print(f"Image saved to {output_path}") + plt.show() else: print(f"Unknown algorithm: {args.algorithm}") return diff --git a/outputs/output.png b/outputs/output.png index a23f6b31d73aa1a59fb1f2c5239839ddbcc00bce..dddfe4ce8b1bb57b52e75b6e3c0433bf194b3706 100644 GIT binary patch literal 5556 zcmeHLeOS`x`v0meD^uH9Q!B8xo!03Xsauwq<~r5frn!@rZ}=Kh6GTKz6jE$%c~1G# zc3K+wQtK!ZsFeXB0d1C+CYJJ&AYwV;C?O{iNKyEG?YE7x{`>v8AJ+xn>*9Kz`?>GW z{rTMY{e1uV_wW$s6>C=j0O0)jXL}<6z`hs&>_(O?1^@ATIP?qnfamUu$&JLOwkGJmF>bGUXk=)#De1Nw%=Jy3#u{r79&+6)8K_knv zKa0f!fNRRa*RCm;ngIaszW4dw-3Rk)w2EMQ>Yci!Bg|5_y`1WJpFfrC792m9{MGx_ zYbfzyEBEeu`yEbWkwb@h&4+709N6o#YFoxC-;X^IttS6MuA_W%*E@p*9c34jB%X~X z3HsjA9u~(PO0OU24HZ|$aeKgz9Z>9HWf_(LKui^|3|RFkumo840bmb&^EO}y{Jabd zI%UZxK(UYg2H+D9yH&u3vOmnS`^@G|x5j1KJj9oF#ml#Pd$a58*Il*$f0X{9eE2^w zKVD9>2N3O$hbH5&dXu_t*{zd|S4EarRaJ%jB6;MvB8zQUo$PS}cl4Ma;gqp*pSg2b zU4IfC?5<5R?y_FGVU&KS|GMEH0B4Z}5zrszya63M6-+x%s_&r)t_iV%1j*by?}!uy z^cY+BnFUK@N|&APS@FA;Vab$1`?PTMe8_a@!(pY7cOc%2U0{+kVagQ$Ar*{$%)8>c ziT#aLeYYlg#oHUI-hCza3km_gila|8F*J0jA;QxfA(_##(SGDq^N%;8SUt#OlwUD` zzw)YR4md1hq|w?1^J)tLhq z{Rltu0xHqRJWmk8We~}X7-i`;UNxC5#5?LzCe{VeT>Jevq?1Q=Zv=ubaM#KXc#pa1 z`a)!a_-6$k6Y`txKcGAQPOg;)I~HSN?OrK+O?FEtutWMCD?75&b^2;4x12WPCGV=p z<21#4Y37V*0P74CE-D3JRJo7C9}c&gz0j7d)!O>j=YcB?6hVqV!_56l%MV)LRsuH& zn0^KRW zLlj7qMGECrytTWXe~Q;PHdvIZv_b)4iETIH)~&;+5T1y2tS1lCL-2567gR#cAE^vk ztYgRVsjp?Iy0`qtgq;&(pktj7?KLGBOL+nBJSBL;w>G0M`UvP!Y=Nips+21R>6s4x zc%L-3MjtDRk;jsOD6suJv9%u*aNMCPq(wPo4s29mJzSLI2>Q&r6F8kofu+wl_zP_) z!ICdbWG#>xmNXwFpMw2LW(n^Fs+bvj$m_gnx(u?HFKmtibFPPKBhZU%{)&9|NZ!Mi zS_6Tx^a?{S$0n_ob}xABxio#=F3g5Sr<+f29#W5<&eUPeXA%6d?SA^5R6)jwwe5G>zG%P)96oD9OJ~jc1MF>rdfgEGh)y~$ zM$MhLn5d3U{9hW1=bBsROoasnoQGj5+=O092oH9WFCyl6>pxEY>)G?s${kAX!RFy8 zzhZSf>3Qe@hZ{+zXp&g&QwCCDFv`Sk?(!rtW`-M3A{B8m&M#NJ-9>%pWnd{m22o_? z5aGVt_!+H;!#dno5K8dq#H^Vgee&c=;T?5KwUNk1R0e(^@oBt$$!OZ9!eX9AslM;a z;$fR{5~7>>CHDUcYHwZ9?6_Aelzx$(5X>b0w7F3LR3Ua$&b}IKs3QPUUpQ(!JkM5}yvH*9U>fiT2a%YWfLY zc(MA7v~e24Gd;o^HP5%6RC|(m@?Nj>+mgCtQ};@YapW(No$(r_@KWM=eJn!!#rJP= znQWN}nUwE=D$fM)HmmFQ(<(0;cx0>P10{)}vlx*0c59U&@cEoVY01lOX3TggO`*Yw z`H#OXF;)wxOiL$q2toNMTvO9I%)O2G91A#StLBgF8wc@yL4+}8P7i?=T?b9A z3$SSAwV(zpW+>;B^sH_?|H(f$xyVA;mXI^c%VtbhL~hxx^8N;9MGwuuhTdOqGx$PA zc~vXyk4%GTK(Di##E8q=v0+N2NhmfA^GIMyyLn?w881~?xoh^wg{0P-Q{5$fNKJf~ z=n_jZP(j-RAz@uZ2B!1NdI`$HXlRw*M)wyp5H{kAk3(8Id{$xw^Tt~_iR5_$Xf)bMcr>g( z4sk0e)4dW(ErpF>`7@0~?>8|oYxritgYaEBsWc*2!cVm}pq#0Y{dx%k>l|k4T!#&O z9;j0GS4-EOo$qU2n%xq8e?Nynzm&nBouD5NTPRU4O~@!9BhPHqmKfKn_0Vs*OE36H{&VpX=%4D;QFta!Ft((8l zMTtlM1e*)|Z(QDh65( zdEdC>k4!qZL#nya2njiGS>ZI-qpm0O@8Xk=J-f57u4jBvsld5x=mphCEipzgG5Fez zZ1GGUx|Xvex_0b)!+B(9fx4EkG&%PDPquo!dNc4V-e!JWy52?BeDteOi5vzVy^p7O z1ZPTIg@f!o%SbhDn55&`Ro5?`KMJk$CPyukwO#752@`I%I&XnfCAbP+_hg4_Fr@g| z-s1R0fi3_&KMUr%n$irL&=SuLL>TfSE8=58PIvp`4ME$xDuXTg#RuQ)=B%dGRfZO+ zz*N!#ZfrJHoU_73^75I2xUMz&A*3mM+O{Esd<=45L{mSbtqDqi7%JpN6v$YFSHsPBkxgmdB zZxL*$2Hzxp1j@xQyb1F*jWd?203WDLu)AX31sur&m>( zbc$fJ8VFz0HAoIDvrQ?8Tiec~o$|YlkD@{@*8D7Kgzp!yMWfHQpflb3W*DCAZ_u4t z?v({W320^iRhTl;EQ~F+ZU1%1p(wZn-30|P@GeLCoxz_P;fj|Xwjgy0qa36T8maV) z=8Z)7mCf};b@iU^uyxztx+Brw@k;4z%SQD!sp+5$6^@)W`U~Uxpn66mqj#sFoxfSF z?>$7{A+1ieMp_j<>uvFaDUI@U_}5M5;OKW=wtX9WmPl?8z{cr|$AtEY*$w1A)v-c^ zHn!sUW7pT=XB~7-&Y-pPW&YF{2Y;zQQ(g0<=4?B=lp;V%0{)hEFF+8gM z?+x*n)9khwid;sj)t>99Mn+^!qz2L$h!*O#ERS?y+&$*3WU{tRx3+Iw9M^h8|5ty{ zP8$l~1w|nsTQ4X&e0`Xt0Tx0Edalc&;lfj#^6%HEp6*O@!5b1Zv1%4`{D}M6w8aDI z=`%cN}nZ zrj)}SuZWP&Bh`HMcH!E>9!=16kVH_b%MjuaT@W)vqHcOtVu)FMiTHPHrC)KVL5mA# z^1EiSN7E82m62MRTd*b$ACT6JYy8i}x>eyn*eC-s4uNvIAu+{}3)#>e83+O$*La_H zw5&npV@SMEYg_pnxiE(uVz`aNlaHk&21bJm4=o?)fSkQqh0*45OY&y8GQp(j?&48k zp;AOS#7tbQU~~mA@bv?3qgx7EemR3vwG&!!EDE_2|MgxQ5i;!t7(@2M>|Y+fi^q6_ z6qXqxm0%OrZH9cCv!ZoRoAYdTVo1pQHsbvzJi0iwmdKNR-}*HgO0ANy@brtJlt@i= zsPOZ7XyzNp-gtAxnzg(>?7f8MEvwMN#HUHgLA@Mb^ZW(cR>3=t|UQ6CMyp{2r&*L(-pYGw2eKHJ)UnC%=kj%JLOP;QjC5SVXt!Zksp1 zsaFKB6kLfvV(uv#%nV_u1EPYYmnC)%PQ+%r*gOJ9zSt?Z|Hh3@^AE{8AZ7SOc^4vx zuJ%n?Zp*`q_i)_Mx10gsfNiVoj`YWQMQ3&NOLXeGv`PCZ_y!gDd|&w9mcJ!`{Xb0! Bv>*Tg literal 20619 zcmV)jK%u{hP)(bWB`qk2@3m0xZcjF~RF4QrQW)>C{ zW)WckVPs}6T)FD~_ud-xcPl^Dn9in;lW_9z;nTs!<&BMt=pRf*w7%}jI#kiRhXLk4 z+dn*3!KKxW&1!ydGM;z(=alIi?IZoP$mgKWm3877g~73a2#bh-$mtgr01y!Y02Tp& z=7%7_APCU>==2c+0KOl6Z0H~g77veJNn#b^-r)&5+Ns8Gk=_*lMD+REx^_A?KKbCS zz1^?A-fCS0L<9jKL{{imf;(|=Tfh11|Nj2yv70Z$I2<1xAMF*DU+mw0y1RL})R`3^oRz+`eN=-Rky~UfnvTk1PV{opV zf2u<3tQxGHx40hS$0ZEfD_6;cffn%heg30fg2yu!yoG=Qf(R^v2q2P0qh;Chpix8s z&K}U*0ss&M1QZcK1QA385kU|jp+3s!>G2JKR$kPFY(>d&iEZa!Y^7;cUs6dIMY`tB zxavOl;^qPZVNVJXRT^4fK6J2=weF6OpX#t1*h6{#_3!w6cI)M@O@_l|G2MUs+4uh9 zpKfHT)!8wsv%Py6Pyxc}ivUCrMo6M8?Z*d)k2}5fqemybR)4x!^!hiWaG15SAYEn8 z+JmLlvMMwDYmTPVX`Xbx2)ys?zRteBm^^N+Z*cwAY&Mz~8*Huukboe$M9qcAg&i=) zb52EwdKhWXKdzPRWiRW~!NFfwlVfCBS?h0~TdiT4 zrrGu9UfkOKLRE(>Tl@EaadL8ih=L#pD1ab{NC<$KWu29c^{m$ul-YEGBp}{dr%d3b z574ij-RJPkz!~yVXnf5S%e}ruKb3>_tDZE?n(*gJ4fjj~-7)kKtzhgC8At`>oOV zxR+e+t*vXrYa4_9;9@@8{Ke1ym5R!4XDiTcr9%-wlqN?YfCwVP{*&8DHehgTgB?Uf zWRzv3pDbs0i)p9ZzrHBW&E~gjxhkX)ejEnnWbv0{dmKlXO}kZ$Puz0(0y9mGoXao^FTUXCEBp@JoYf&43CR+#tAqWU4)L?TY1IXem03|7T zA@)!#%CyS?%)FTG)Q?s!Y_4&UqeLM<<=v=l_KuI7PQJSm{P}X`H27?NYjF3&pD;`U zL!I7EwS5xB-{`jcSrX^VyuWsCI=TP-KmIQ-UVib`i@(+G@19;V06Yzqm@vA2>x&=0 z^@Du+G>MYQWKWw*itAW!*7k!=@74KoqSC8rac{YNG+mB2H-DuD3H1+!e_rvgweh{B z{m2AYGx+tBqaPEu!ggm-QI=g7d&D&8Z@x%q&XRd^*kF@8n==mukyMlU0RaH%^lm|= zg;Mc0q(aBuLfBK^A62h~CZS-_4fdwR#k>kZ$hpyc=y5G zkCgs~J~$X|Z*OmH?VKE%@zEWdPwMHz7cRW;;OIf?(zQi7mid+~Rj>2Hc&6@t^0OD- zSTj0&b__KeLLnAK@r5sZ`S4_aeE6^(&+oi<@WM;aEp4~m-kdESDI&2CMTh6MZcHY_ z@e(;(DHj#cOP#grl^-Ma=ZjXptgKsO>5*?$c|n0nTVDiMDHAsMpvhKn`sAm|Cjtlv0wG8e7NI&~;Yc58H9I&ar_+d@C^N74f~C%C-!4b*7s<_#;tjr4p!FA zS7o7eR4kT3H#|L@o2w{5h-+(GhO6tB)<3#`}2oVcIKy(%tEAaQ-`t67sI3Dhr^M|vZ&SHT8tM8IzRZ~pR8?dZg1oD z8_x;h`p&EDK>&}AkE*&tLS#AJ4X3XUYFi#n4{qNbef^bJ1AxbmpHxZEiN~MaeJ6}A zE$ij+&kdiMbXBk_9ew1=K5lgTnZ6gHS~T+1h$ry^$XgB4IXuF zj-Zk76oxj7>2wbwAb3|R6Ur%!z~&1H7&VXdLb+iOP+k1AlHXlUItCsU(_bXXYdO5u z?v<{drS7wS_INp5v^xn`A1a`1tYs!nN1q zesnTjJULk0y?dvwLelBMqq{nAuyeE%W)w$?(a;w|Ce2zx~ax#KzT? zv)=vr5B@c%d0t;iT1=`C%5-qzhTsU?+k8M8)4Z7` z7Sqq-_HGpSo2TJs#WauwgjoQl^g61I0CIp&xO_Xri_7YfH!G3ZmSEt^Cri5<&^(KWM+XxN{jB?QC`gI$#D2EG*_zh()eUH9iNo5N2~o8*?T~B z%2ZDYCjd~Iw4&|g-1g4?!M(eC#f4RflkE1#Kiu8k`s*M6?85ayO^afFhkc>LaY&y% z{Pa#I-PO_7Y;=dYF6JuCP|CS2@@!)#3OnA{sqBVv7tYpjladGl1$;xh0EkFXqGWZr zcSoT_on50USLBoFEQxIpANwjRmm@_TqG?%XVk?FRs7Dr~(iR)Nn?<48L3ZpD4!&5< z@9SVU3###Gg;bUX6Ht@U;X7z@O_ZO%x^X-^{=fX=|MUOo@BYJ+#l!vmz2EuOzZ-=# z{rKefzW(M1_a1skD`$l#qrH>I4}bjR@hj)fUpXZ_XGbC+5D*z8z@44z7zoDU-0_X8 z-$1r+y|q}~y7A=TgCOV^^Y2%(u~=?(TdNkv`Q)R#eBA2XT;|6s?Qiu5w?Z3=AOI_+ zqqqyG8W`LhfCzwsjHiKXA}EOMm9_cNXDgS#Km;+~o45#R*6O6w;>MylMjj>R!NT{0 z;Fb>J@&1QFblhFJy3)E{V;Ofhs>5nL;d$|CteAs~$kOgB3TiRmD$8Z%>Q-O$)$K3c zdGhwvb6>ylH(vwLt_+rW@xvegu&(P{H*U4M+3wC-7`9m1fAhWLDC^(9^Xc7BK3Ey7 z#&K6pq1Gu@Y~GCr0f5rM?Y$}vpt63DFa+%F(3kfBTa|xf@6p}k{mJ>uUkJmiI@nL} zW?jdfu6QBRy-2HudjbHeQ!XfSO61RE7(@UdKn6Yyn`~@tsvr(jdAx7?D-pux)8n=Z zHajiSus=Lxj4t%rhl{0c>0-RZIN2|8ps|wpumX+wesAZcD7)bn2d;dy96!uDaoGQT zM~7?Y^13$q>)(8H>*-*tcfD3VO5@qoKmWq!{_(?K|JHBpjh_DT_x{yrGW^nuUmFZM zUwQFMk01O~TjviReK7gsgKvKO@3ecpMnXaW1_Y3X?lXY(0tmJM&-lT~;>p3$U^EAr zF6IYC855)cWp}mDn3efvQO=T-O8Fw%fnMfo+Wf-005ygV$4WDLT7+W2vHD4 zgn-!^0$Th<6)!mC%GyEq(*l*aM#{)eC4-QU>A7W2cCr+*3#IXrs&U^^Sn1=8+-aLm5a!ERk(uFnx#<#Tb% z;ob*v@1}u~#a6p}1><<}@Wc81{qD-Ssz9ewC?0qE*Ih_#$6Bw9i&&?q z7Q7SGPkeA*=_`xTUq_ppwOa&X=&6+?TS>RCv8>=C#qIg=T_07&ojTdx?QUFM&OS=c zeR=lqzSgF_cB7ijR{Bf5oE3X;JU%&EU5oEOJQ$AV)7gA`@Y?wc7hZeq3rgYLd-r;s z{!iZhakaqF$R-obPX>2jls)s`6Ivxj;W{uiPM|;0!1r zAc4rK@(A2SOaKtoDXngHKPY3aJom=kxBhZ8d4G^PE>C8D@9^Nid*j#szAv7d)tx{N z#ZBtGVo8_Q$7UVPcIfM}DD%8ReQ^2a>rabiyPe0KE1eXV-MrnsOjHQ&i$(ms|M>5H z^zQc;^C>$6T)q0r@0^SecP_uw>%z-7ZeF{1alTl#;!MZt`+xPrD;KW5^Wp73|K7it zPKW==Klo36|M&jAjrA?ZJT2=#{PSO2yl|66oSRQ4PvUeF<7FMiIribtEUjC*BK{JTk zEu{ki%~hQ*j}kO@-}<9-U-?}rh6&2t>&3~kJ!rN1J4L?KpS<(^>2QvyO~{Dpy@&7C zuF@tbtAVrOpM3An*9V<$yR%r1p4=O~`R3oM!?eO|Z|~9n{h$BgXu|*SfASBf7ACnn zzk6P=H~yRWnM z_Hr@l4qopJE_r9QR?bhwNj`rXXjJi4QUM~OFvF?BK%cKI0clj3vojfig;=R{<4WNq z3|6*pt=9FN`&(WSxOq9Xvm<@y=-)#9AlZI(aQVxM+SB>b_W7%}9-6I}w=Vw5^6+Qr z>Lpw3neh6--d~_xEM|vQ>GAR_D~m{h_4igb6t9<|9$xnXrr03Xl%zx~c%bo%Q`*+BFE`yc-|zw-zG!PMC@^|__#TUHC9bAj6Rpd$Xr0CAEJVM2cXJ<{ z-b6$|5pEVDvj^cSU&dJr%Gt86vhI35KCT}8s2u%bGMRt!;5g`Rtev~Ev3+@m7LOmi zjlq_ZaXEW*|H)KqHJQ!uZ2a}xl(@X^s>6h-@oPopr_D%juuw3wW%ZJqnYPyga@|M2>i zTh`WJ{q?^$Jo&h;j_PHxvh&7hT#S#O9`3yrCvdp`m{k3R*ZxirZ6KMtoP=@L&h7;5 zizZr=vsJ_(!j44{2|=BO3qDU2UAPM~2ErzxbD18GMI3eV$@}fpDkH)Rak=#p#OoxooZVk8=5c#B%?1c`iirV4=ou`dhKHS1SO7qw zs_JPyKdFjo(q7xS@x?T20}ue{kaTF0PS)zK&5qw&-}#@BQ6A~$?o0F8-O&W;w#ayg%^ZJeLYpT76@{?)5rsNECC;b0@| zukB2SlkJV|M&u=;BuPfYlas^aPu~5BpOvqC_1pj9Kl)#m`Ek3u9VO{}@(^jItRRfd zFMj;JopaakeDvYb?BR4=R;o2%=K)s<-P(T5`l3_V#650cit3b%^_RJ@Y<(p5th ztfX-sL}8-w&*Z=IRZ~oT#-JRFp z{M$*poyN)G9_g?(868Y#C#|d%1j*z3?*=-YC2JcvS$jK+V{#A2GwgizuWBw_4%LzyJFOIn7E$K;7wz#&|lbJzMW=nx)?5fmwa) zXa5q2d!4PF^A}1>j^} &T=L9qk{k?_3wGZKdD((I1Dg>aCs+12dhE{_%hDzy9)T z-~C7boBs;6TJ85um`$hS{ur&+LwQ4*>n_YTZp4GE4vv{pjqZZ0AJ6{&L0HkTJ_A{72{ouX7Tpg^G#nam#nJ^h#KDY6M zfBQ$*UVe4|(Syy+i!VO6on-0VJMaJWFaGs!{N~@gd-sF4fAHTI)#6KE`*mZ?#`pa+rw<{(CdPo!R2dKu1ca>4mmEeWjx~EQL0T57nortmD`8?f$@=j4 zsFkhGrY;Cqgx%4>M<@IDUw-YIpS=6iSHJwN)n4z_=dV=7?Emm16uT!3V zR;GqigAFum?b@2kG))S3v;-?i*a%7r@mX;e2+soBMEdlr02sXWrHzx5c^HMHYHcWA zs4AD0^_?p>ZW-$Lce{`7KU!Pgy8Xf1zxN0K;HN+RldBhA-??-l30q00HJdN%eE8q| zm;c%B?yGP7>Tj4RJAdht*5Pb1`Si~18#iCdi{etM+n1l4PY#PBch&yon_qi$_ea*PF;UWv`YS7o z1FGj|XFE;K z4ohmZ!_#%t-1F>P00HlpT)YxT z850rd^sBZ_W$sI`tad{C^kR-{vWm4fdZ*1y?*^>Q52Pn z^5xQ1l`(o__tNUx?#a=Ux~>kM-0uzgfriEGQJh`g+`YyEMky1v&fZi*vCjxllf+L8 z_#g<>EG_^<vl2x!H<9T;aflX#y5WRm0x+2wEp;=w@!{9eDcx9lqB;(=!!%siH_Tn2+nhn;rdi~YolM^I$aPXwpP8ou-n6RAYX8LU3rwE_Svbw0mWEy^%mn~aB9 z+QSQ9=+qSyxl>AXAjK3+49cK)qC@n|-cns#W+9|jYjEw_SHAo2{hcdUZp5)U+}|&% zf*iJ1SKs}`PY$0Re*B9cf9Jcu{iSbwySF;9)_wf(`vxfKuEcR_jPkBz_KUoHaQF7^ zh3nJVuvpIT-@WtQ-~JznfFeYq;dqoJDF8AHAZexQs;s>y#OZVtM`5eIsz~9fCFOHx zjUi@1cV@78_78t)<(&S~X=LPd&o$8C3~Dvn*lA4(07VE9K)kPT=TDF%TeYJK}R{_ghHr_>RE@rAFiog5WWq>J3DpX?TruK{^9rk^iRI>D{ua%|J(oktH1IM z3o;yyC)4Q%AAjJhvfo+z*02A@Xfiz7fBK`p{O`X0tAArWI_j;g9qvEPmy1xV{lmTa zVtDStRhBvkND+_-06vq|&(aHc?}6)1D@FjL@eEVpvvl*h^!E&){t_q!5CB48`8+Ru zzUV&xwNrorXQ)d6xT%0QzvsQj?slAYR1m4@u?26|c3Z=fAt~||A<|@&1Bwn+Z;&Oe zZnw8t7RBYOuYUN>Pgm03>)(8{_S{?9*jV2=8SkIlxzO*e=EZV(eDdWtzX<>=u*{1f zFc0p1R@X&Y*PU)pX^jGp9^c=8a_{`*=QcLaXIUDAAtBAC(^k8S)F|tXVf9%-`s|&Z z$)iod=d?Wid4Qh<`tw+xnIg_~=hF>$whDiV96qB(Aac6wK<3LumZVP)?qYYd-Rq@6 zs-8YP?sSqU1s1WDKi)5cKnDr44^6D&xTQi192Nx+_xCPbx%!*G`}c3YaO>wk{mV`` zSlwN}dgUex4Ts0|boQ<9{^nU>6f5f#spgfHRVj$f41iCbKAMb+kF+eRhU6iAfNE zxe+CtFZLgOqP-tY54!!eYB|@P6@|LUOGMB+7hBWWvRL?|$0bD&XUWZ%F7G|Pvslzb z66)6R(aCZ$U*8!Nb^i9xzgH|KcR#(ey0LsRe1w)CJpSzB=7rz=+kXIODMb*HK#4%p znxW&caVXW+wX>vL2F3*AlRXo4G(a~^p1<~-fCTES9C9{k6_(*LUyTx&4z5-d?-- z{7TkIgW+Issn_Z3Kl$uv_!uvJq01bb!0kPncUR4HRM`^5LAnwEavaAx@LID-aQJi* zNA2@huRM8jZ)^Qr&2<{LzWVjwc=@GQRls*YdDzcdZ~of1BNLdwz~>D&AkF1(N;3=q zP&>=a%d#kna+xpMS*Em6N+(fVRwWTdQS_`x!R!#ou_sdOy@;^)N~^lE;++No@4bKm zma0HwsIyX4^_@Ep?%(n-(9)* z=Rf)@tBSSuX3wOb?%(BL(Mq!QmCGyLYYz|KTMMqJw!@1r_S&6jIGm66op%f(6kr$! zLL$tL9f+*7x|7LL1#t7ZmjyA&qQUy^ayEYQ^y!V8FMRV`zj^1Qk1m}*_x#JRc3WA$ z)8nQvCE&a}O&X0UT|QScMPOML#srA4ST5?in#~sPfAH>ypL{S~BNq!Nez(P-5c%mIQEs%TB%rvVYIT*+v%v_gyr#QScRbqBTxuQjAd!9ui_{`P_v z((SFZpWHjLB_lw0H4QSW16`I56-#3}OQ1tNpUfZLfBeD651&3guyxMjI-SA6AL(&3Zwr=#+? z6?cnVZEooO<7t4cNQbP-x~>S;hvW9aqtCQAT3Ean_Trt(r+^^rMZ9G@E*+>e43os@ zD8zJaEnSCsc_KXbUIK&b8_~&dna`@q;nw=)=dWHlc=Tvej&N{3bG76PhX7GVo98yl zI-iUu-lHphoih@~Z52k~oFqwVOb`U^3m2cO-3*xC`Pome-F%_5(yEI3=B-Bsl)*HwOUw7eO(RTefn|b>dis?_-T#^C{*MWXdJ8sd9_1xKN}tq zrDUsmIkR=4!c>UL`BZ%Ey%ivBBFi?ftf-8F*aT4;VAhKI=lfRsy7n9ED=NgU>K6~qA=-)Tnhu}K=lXRb}xDQ^l6-8wN#676~q)I$Sh$*q=Y@P_^h2)6?se5 zylD4RFwknXvNW`%lB(+qqsSmKGsQ_<6s2Q9Bt##CQQXR2d;M~|gSQ`^RF&<;c>BTa z!j9^a8?85plR96LCaqcNK$rz0TS1n#UVY(>-fE|^jt~orp;4N^T4q$fSR8$PcQKpK z^D5%1nvIHL>74TdKw1SMh!w{om>APbVnXo3B_wkv9;c{#gLc%qaB=(A8$tmp^0wAS zBY-CMBE}f33k+M4FjJt?RN$ra)OI4%<*?U_oin6Ft1-4j(QLMBHlP(ICi$$-q`{iHx8#+uxh(!cJ0Wbsj z`q59nKRi6l=aa-+&Guw5XArG)WMW~jP#FS7EzF1z2O+Qr6eJBQ$n9iWx*!U&Bwo4o zTKB>&!NA*sRG@RL6$9+kDomF=|+uDnuv$n1vOtdZCnWKW(0xARmaX<=5k+aet41m}c zbv|E2aiUBThgYMpEdo)P5bEi4A4w;1o0+|FqEL(SO6kR7{`Ao&{cLn}?1=;c1Jz&Nd!X5?ARrUP z3ga--TC?Zex-bmc)fs95Awvji;NfvXc zo}8qUlWIQ6=gTlAMFc`C{V3{nPe!i4*=u(bZ@rWSXaEskRiIoe8k~Ra*GMVP3@p(P zKAF$<*1H>9?e%UN)y{=l12i36a#rm*F1}YZTC=(@7+%c2!#Fr3+VcL%C%1w?Z zCeo;NcXh>t$p`Pg+wMOJqqTEaf9Dzs`wHr%r3isUbQWjEgi+kSw{!iS55ISC{Fyf4 z&iR+w7wrzsY+ac(y|vvI=Ut@~19PFh1B0=hEi5mV<-!iuvZL{=D(a8k`nVXDn*&ub zpa=p12%*{vV|uGA~PXHZc?{ zBkaY23&J3dVqZJweW$mzb#CY0?GM*iw}W>4^yty*>W#>B<17N_PKNh2NkCC$d0HM9 zO#RiX06$EUuJc8f>8yQa`=WrI9X-A`8WxViFpHHkx)pgr56I4ubK)FX)lH(6%Ukv6 zgtPUMm_>0EXr+SbY=Lg+ zfihzxP?OwdfhnBd?zCKxCfQ0^=C1V1+Qz-j*$9?ZIh+)6yoOF`#|%P>bm=N)0P9@k z(lqV$V&~jqk!!7kAP7QbqEIP4o6Ul_?aJcOCm$0EGnVYRT~w(lS{rL8BiG&R zg=wmk?yR(l)M$U^oQczLYwKKDmfkZF0bqM|x8DE$c=331bv>`q#C=hrW>KURc=k?{ zSX)iRs)|vnnMFmBV>zl&Gmr{}0+W~%M)YtpUuYdUz=?#wI&E(Q${<;dqv$wzs^b*$ zxdRlXAWNdUu2JJ;HeoM1Py_~$5Xf1aPKS&TCn2yqe}31N)%}k?7;J5m@VTAsogI1S zt)G_FP46c~{$xJctE=kZ;M28@n}LbfH?Q6O_|CL7@(htLKF}p@*;?&6wxUNb9K28S6|xNocef? z7s7&sN~yB4X}YSTFfoKgP++~x3-)ZCGb&Pnj>AwXDywobAJ>(QqEJA@Y6to7@L=ow zHI=l72m252-mRAT_J#AeZf@`G9XvT&55gqt_4E0?Fo*(`N6G4-zcE|C`uNFP*Pp-9 zUtb?jrcv1HRpHZ5o(?W=g%Rn1j1gNGX*(t&fCx1%XSS|XRgo+F=wNpA;K@rjzgasO zO=gX@7X_xfu{%3?nwai;ho8NE>E`a*wekGm<0qfx<^IXM6p50kBcPFVJBTW0twh!r z2FU@6)I^dXCz^IMb7b`!mu?QaE5p+AvQeW60EC$ZoOk%`@BFG4UM_L~NRvb>trW>= z66CTht@F0F!e9cOWvL<+^;e^))8AO@cl#_{70dU3{&T6zPPdPmRJdtz*ZU+2T@)0p zWR+<}>sb7gAUs$Wvm|L3`P|h-*6w^Z`Up3BD;IXcEOhmJ{Pbzs?*(CKtIAi!{wELX z-1avI+gHy`9zDtL6<4-jkKz;onYpNHBtpc|y${>7<1ooq(s;Ag?k0&r@pb9yg{|fw zI9<*n9aQWl#j={uy_?k%2y!b)YF1m_^&m>lZ|>9-+#8;>($>9`$?W{iwADdGfzxpx z0HQ~;F?$BaIF7uyfDE%=6gdbKRbg}xn4pzrfiZ}PL{ZvZzi_G5?g7bobfPrbyngA6 zUpYBC7#$xUKYcRTfQ#3qxS)lPWs z;r@dOFXF+DienVW$J2Q+39{b)M<2D;Ha9M8v{zQMR(9|1ysOR)$D^pUQ4lDly!Xh= zEYV>1>Eqkccz@AZRZ0(+qmBMnWKt6(VKz_zy;VU3IyC@_uz1fp43UVjk-xD*niS=z zsOK)fGns8)i)-r;Ar1nql>m5G<8S`fuM=qm5C-R+t*!OeIa}AR*~i8NNfLozQIzxf z0+nv9tlfJ4B{Fbwc--!Ew9>6syR6GHud;UQ@%_&ZpB#qqRvfpC4wTYb8xUKSC(GqY zoCVf<9fSa+k-e>>v~^6=^jyz$R!tn}P>Ef5?`HcCvR1mb-f_^8V?X(1rcF}hqcmGt z>GhS;pNsMU)1&?B;k)SOtth*E@s%`6LZu9fM$iBXg#ZYF8kYlNquc``161A{Z7R=) z)9Hi5y-I*CKiBRLytq0)WXC#ePp7ju%J93t`RmqI*14*z=Zi&Md!;BaK^#Yl$e6$w zqm4UqPT;R&NJXhQ4%O}O4SGzfeP!EnFd_gy3A3mNmq`CCt21? z(~ee30F@#|gb2qE@6+TU(N1Bu(&?=YHVt7wpb1HY5C{+vRKvfJ1fceQyj+mhWo-|~ z$4`%+W=S{Oy~L~AL)bcGJ1fGh%36e~nxFKLo0twca;t+YXcx~hv} z!GLiZv8C9tvXL%xVx3{hWqfi$4w3?S$XB9x@%LVD+ z)ChypO5(w}3u8sogL~~%O_$TzVw|-%`kii~Ppi)aLfE+dP|)Y);_2|nxmpl3!O3LQ z$=d1aCaiAN)>`io#Jj?lGl0-ZZR??)&L&k|S5;M44iGdMU7Ls^rS*Kd&`~^}j+}Lf zfTZ$#*&S?KTY1kHZoFu#`TXcG)}}w$7#;5i9Oh`4gE&%Y7FPypTv`R1qib zJf9Iun0Cv$1b{%7#hhzjIG}k|=B`8nwzbkq$EM`^?t}LYn0A_lQJQ3dlrsxN+5|)_ z;926(Ab@wSv#~u>=IG(QL23>8>9nrPDd{8#lh6c4DNRU-faqK~96$6dB2^R)_KqI) zyQ{jhHcD4;wgAL93<9IskR=ZoBKk1c!0Rtuvo&i?0HTzNqcn@tFphN?RF!z;^Z7D} zB3l%L&0SZ`vb5Lj_p}O7k@xIfdGz=}7`D&^v&pQid&UHWplEJXrBq;o6{5s5FPGEl zWPiChVsTLv>o5>;>^aO5CuJPBl_oFV@oAYs0HP>z;Ebk3w5o^&hyqg0gfVTREQnf& zR8$p_(n=Yl6%l(kdT=L|%EYTZ8aWX{U-~SIJ0|KHr4^wF6M~2=i&0t594}?{`Xb2A z_9!A^9EMS(#CwIzXb`n-B^}BzO=C@(R8!tWBqZiEB6gT|J14_oJ4*s7jW(`a%qL*$ zDu}Z%j8^(xbN<3R@4Qpw5=K337_bI$%jFm;6L->Tm2s5iC4sMi8A*%xbzP)cnlJN; z#c!a6ib4Z17}qeI;=CBFoy7x{#{o7ldFI!N)7|Sw9v&rxv>?b_^DBt^){lAK&0DF%c5YPM**;e031gVgX=^=5JyD5siPwTdxyY)&@`LO z^Em?*MM(tJ!7^*NO^_-gjhel40<*%zDS4b_TR zlTuUwunLS-%2I~?)ig^4fY?j3WDyZT69FNDIBQ3XIT8|zRtlL)0Ls>5 zMGyp07zBVsr00u;U`;wiU;#~^`-D)05Tq`vgYgjvsz7T4VH8$u{dAV3>=^>90VEAD z2@=O*9kaK^_^?>a*jYqK3{sRZ20+AMtvJWE1*E16$_HUK(Aps2a+woh!+VU@jjqSu zH;5OJ;y^Eo3IVhTc;||;Tr6{ItLpkg8@!0h;))1pMx*g!G)vM|7E(%wMk|F%X=4nCSEQ`9)7fw`Ic#2p2_ioM%!voz96A z5qaktxk)P(Xd@zuNP~mfag_F#O`kf7Fjx6Bh}xvo)X9^Rhhfw~1QUds2t|a5;y8}N zaJHDXTJ1DVNhxvsxz$GmNg<&K3m_5JR)#0XCQ9NcCZ!jPQCa0i={Rmzh4JE|AdMpo zL*ty^+utwBF<~=44N+?ZC(LLx5-B1UCXjxo6PN&-rd2P1po2gu&7PHz=pYLr)a(S) zI1*+N5Mf{E!ouv&?0ZTC07+{CE)`g2k1FbP+ErbbMX9x_srJr?VW5@LTB*8ig-F47eJWBGZbYV z8dZC+)@r4__sm`!qqM5)nzgR#LRUueJPBingXL_9-d9y!&PLu=I!wmLc@$+)662S@ z`VE|mKhul<$3L54mTF@bqB9?wVKPE)owKcNE}3jGCLH%WC6|cKaT(3^9iTC^cd_Q0Bck^*_f4WRqsbK+9gbN$&=2<#b zlUf@P;4N8=v%0lP{P6ijquuE7Z}yVP(SQDZssu}x7hKSLJQ!%qkepqw7krJl6=_N4 z^o1vOyzSI^q#Lz%|KXFV%v(@Kc)qe4Ez26xa9{3In`f+y7ge#iUMrzLYJIS0RA7mT z=opbN|LpTu{fbvVcnMKtenf@Ju##zBuy%fMP`f{niwzsWpJr&swyyV6WFO+oEY=qG zeavLOk&XtS!yvRd_qhG;LY{jR7ozO2u@~1Ylr7Vu8?<3>JG=>phLF^Hx!+->tcjnj zoJ=ArhN_9`3mX{FVWpDQ{67XYTAzwAAz6o??8T>ZqF9S;RaD>hWDN>KhFE;i=6l#t zj9`yLcHg{*HzsVI&|>!3$$3yAZ>>Kf@)1dutj$f-%L6ZT{&u%JLM2_!%y!`oYtZ}u zSobgv-R!V})UCzW&XSGZ@v2ap*Pgu69W^XZ3XB!apI97L&nZj}I4A0HE^fdVH42~0 z7qllig5=Bph2ImDT$Rd)j@Sn-zW$ii^Cj@X^VQx3F;bcD@7husfic2@Vy%pWl;IDr zS=p8?})dGn>Bx~vv1((j8$;-jC=2__j@5$9u2ZDX;sP`0BHTRqce^w zE9^)}`9}Noqe2r;0lcL1g!Ru3oZns)%?b%-jVbL_4CCXUN?`#3P?z$t*j(>MUNtk6 zGq>pR3p%?~aunXR!+UcR0xMI?;|Nk(DrNbDj^8^8UYyof-Fg?&U#%Gu*Fa#L(j*d5 z^*9k+3(UD6==95&p}28YdO5;kspGiI?nMLe8LZG+^1}49;sl?H93>;0cvw%)ijyQV zcK!JH!%A=cs}6byNf93CR;i1&^XX_<&fum=3IrRT&6mE52B zV_czU!IO|as&GD;PT%oi4_+GdFe5HF9L%Cdq=b)VBy?8G4X~ES;T3kl&ccCq57jwD zhL3cmXyyXWjg+Nc#)dxVxY0L>n!(?qCJIl*sFWTyA1)dd&=wYtmAA#5b9|9@L)-%r z;Lrm=u^e7k67w~K{3bcagvmn#P~clNZD2>?xcT2j&YqR_)KA{uc4w3`tpXqM3#VXF z_hD4?Ask$gLRW5YZ%@=ao-2;$R7w{X7V^y(`cJjt* zJbg;^ct{AjMhT@3IT&XRvmBk>AdcY0H#{C(hY7j7$*jS3l>Emr$>Sliki*0%^p%zM z)sV8WXtUzA{<)DORcVixki+@COC>Akd(Vk~|H>y~v1N2vsC$rmz>A`H<@yyYkm6>U zK?zYI8HvnR#F9jGQ({NYS5b5*g*u_WzZ@5Ketuukp-+z=t9pm$VyP_o<)A50sOCu( z-n?;?B;;}U2}WU2&8|p#UT$00YhJnwxYTxZfdwnM0 zLdZWb=&bG&nF5&rT5ixhs6FO`Z>>*vB>A-y7`HCbIQA>z6vmG8J;7~(o=fJnL#R4W^Ia{+4pYldbBX?lUt%?O@yFbPpij`Hcseg66z+r=7AKk~IP@Zn8xH&%?>mw*$ zWzuyQ6WYY&$JhB@X($6Diz6sh74D)!b%C&@!;9C3cSddAvh<=fJ$Q4_QnzeiX8{HV)^Wmh(P7bcHQwq}2bGI{519C;e&+}A z;`KH!VVw7zlo6pht{Bk)Dc)4q)6tJ2mEJl1a2m2S*xrsnm8gMA-T__r2NQ0i^TYnd(~kp0c?7b? z_Mk$Z&YAtT4c%S6ia-ANWoe|k^Lc`q=w;HZQKO6f3MBgG4Zl`9#Lbfc(!d~Xn`g?9 ztR{mwM||sL&B^Z7eU?>9(3fqxrDrhX>DAqpENEQ`3PU*>SlN{C^p>;Emb8??A>H7n zS>q%2LebmqUM&-45;Ls4QqVsB>~d6?Hd%{zx3UzXmNl>gxVssE;7-3j#2kySMb=z~ z33{&F+ODxRDRwd0xT!p384F7>MdVtN1D;F!L;5`Zj5#sg2GHPA_czsz*pG_mDPm+U z5goiu$Xco%T<`fdyUOiD;`T?rY@7+2@0*XasyUNQECyjJq<9d34?p2FcT?uTOou+@ z8HcMcuN@BFL@}Zo*m5XhyuUX^)bptpr3Nsram(gxZXEUMrFuy-AhGW6|HmoPGu{4m zv9JoE4Ajr_{_TkT5(j@JX*!JsDYJFBmBBQNNt@4#SGu0*D||0Yof2J$Rq0!NAwwGk ziRCQDcYUA;zTUqrt--=ehV2oO1a!NK4st2ewnmC`L64PfVaOR%Ai>`zvXjl{@sCZV znvECb-0WdZbW~wLO=y|hEV{GIx5|=Ejb==CV^PTEA-uSJx()o%SB188WkX=>MJ>Ep z(!|y^Y+3HeYkHXH5xMkzP5r^<#?Zd{{Vn`OL`pnoA85Bs-}ZF4ps+DEJZ?zw_MtJ6 zA?4xqmVU92PU$GwQ*w>v^5%;9fpcPMl+(6)e>!GmI2e1(&c#NGT8XNg9dZ%Lk#TG7f(*MSl9H{$XTU9<%dPQ$5czpT-q+C$QPfAg_0* ztZqz%L;Eohj?x;ZDY{!nTCJK$$$`U0NYVxM9Ih$R)WXahFWK$>v!w;Tccoj;s|+Ti zf-sxn&M5tI8Mk^Q5eaMj#681QSIe?HPq{UwR*p$!SeQdIA|UVd46O%#zb&)IP&5)J z1fpmSXBzX^7KK|pj4$OXyTF`odHhbXFOo_movyx_V-h6Fr67{Tqz>0|o;|yap-j8g z?=-h7!2p-m{3+Y64_&o&EU83o_`v5`J9y5M4ANhz4||>*{M2_Sp`gwR5?uoOaEray z@h(ZzMM|IK8L`9D)$p<n2ybLB(vz3YY(qsrbp7493P;+Wc&HhIvUoF~O0qK0DPgJliCoK#;>G`E+Rh~EVUC(vC6+@; z@#{_HPo$E_RC`I0#6z17TobkO+#@B+GTVsHo;9gv-qhDf;wEH|qN49z%0i%&W?{~^ zQUjnGs4NxIV^JCoO7+5lqMd^W9n4>Cq^O~hwW*5Rw5w-SxhLp9XVU%n^eYX^7qdE5 z0kN02PazL|ud~x`ir3;|Q2)i>@ekn-R8S-|c<+~+1%{rKZtccd&gj0y%VO7SW>@3x z;r-(w_jZTG3l2zJO!o_4lgWt*-LxnQ7sxrD`IuH))IL;R&d0y{PsY|`C zwiVJhhYMR?2H1Z0?&IK92a4GozWLT`hsa_wu&)%RB$iogwU64*qavtY$bwFGf1AFhVL3aHngK^A21pVj zGI}Lm^~Gan?P!0Q7Hgu`^$b}^?$&6|;olFJLD0-Fce&f(m%8l1Q$|xC1copWHq_>C zOTD{m1jSJ2fQu^MKaGMK#SYH!rzIbgZ7~MX5u$ydXAp6!-D#>3znc1rF|(624|tH- zxrIG)g(Qm`3;=rw`-+}(X?s634ElMvmjBDlBl^}B@4as~bBO?vd!9%=zc~+q-w;A6 zX261@*)X=e$wb2=Vjn3_$h(G@x$!Eu4gk7u?;KLOPAHi*B4dFUbqvS_sVN0DA0y(>c(vrC?mB|)w{{%4?8LUKmjgyeAf%6RR%Xio}eBW$~Z%v*KS?2l4pM_YXQ)-PrY3VDzicg19VV5 gJN`-rbN~0l?uaSzz2`*x-jhMVb#I)zLkY_He^*(eFaQ7m diff --git a/photon_mapping.py b/photon_mapping.py new file mode 100644 index 00000000..befe273c --- /dev/null +++ b/photon_mapping.py @@ -0,0 +1,191 @@ +import numpy as np +import matplotlib.pyplot as plt + +# Define basic vector operations +class Vector3: + def __init__(self, x, y, z): + self.x = x + self.y = y + self.z = z + + def __add__(self, other): + return Vector3(self.x + other.x, self.y + other.y, self.z + other.z) + + def __sub__(self, other): + return Vector3(self.x - other.x, self.y - other.y, self.z - other.z) + + def __mul__(self, scalar): + return Vector3(self.x * scalar, self.y * scalar, self.z * scalar) + + def dot(self, other): + return self.x * other.x + self.y * other.y + self.z * other.z + + def norm(self): + return np.sqrt(self.dot(self)) + + def normalize(self): + n = self.norm() + return Vector3(self.x / n, self.y / n, self.z / n) + + +class Photon: + def __init__(self, position, direction, power): + self.position = position + self.direction = direction + self.power = power + + +class Sphere: + def __init__(self, center, radius, color): + self.center = center + self.radius = radius + self.color = color + + def intersect(self, ray_origin, ray_direction): + oc = ray_origin - self.center + a = ray_direction.dot(ray_direction) + b = 2.0 * oc.dot(ray_direction) + c = oc.dot(oc) - self.radius * self.radius + discriminant = b * b - 4 * a * c + if discriminant < 0: + return None + else: + t = (-b - np.sqrt(discriminant)) / (2.0 * a) + if t < 0: + t = (-b + np.sqrt(discriminant)) / (2.0 * a) + if t < 0: + return None + hit_point = ray_origin + ray_direction * t + normal = (hit_point - self.center).normalize() + return (t, hit_point, normal) + + +class Plane: + def __init__(self, point, normal, color): + self.point = point + self.normal = normal.normalize() + self.color = color + + def intersect(self, ray_origin, ray_direction): + denom = self.normal.dot(ray_direction) + if abs(denom) > 1e-6: + t = (self.point - ray_origin).dot(self.normal) / denom + if t >= 0: + hit_point = ray_origin + ray_direction * t + return (t, hit_point, self.normal) + return None + + +def render_photon_mapping(width, height, num_photons, max_depth, gather_radius): + # Photon mapping logic + photon_map = [] + sphere = Sphere(Vector3(0, 0, -5), 1.0, np.array([1, 0, 0])) # Red sphere + plane = Plane(Vector3(0, -1, 0), Vector3(0, 1, 0), np.array([0.5, 0.5, 0.5])) # Gray plane + objects = [sphere, plane] + + light_position = Vector3(-5, 5, -5) + light_power = np.array([1, 1, 1]) * 1000 + + def emit_photons(): + for _ in range(num_photons): + direction = random_unit_vector() + power = light_power / num_photons + photon = Photon(light_position, direction, power) + trace_photon(photon, 0) + + def trace_photon(photon, depth): + if depth > max_depth: + return + closest_t = np.inf + hit_object = None + hit_info = None + for obj in objects: + result = obj.intersect(photon.position, photon.direction) + if result: + t, hit_point, normal = result + if t < closest_t: + closest_t = t + hit_object = obj + hit_info = (hit_point, normal) + if hit_object: + hit_point, normal = hit_info + photon_map.append((hit_point, photon.power)) + new_direction = random_hemisphere_direction(normal) + photon.position = hit_point + photon.direction = new_direction + photon.power = photon.power * 0.8 + trace_photon(photon, depth + 1) + + def random_unit_vector(): + theta = np.random.uniform(0, 2 * np.pi) + z = np.random.uniform(-1, 1) + r = np.sqrt(1 - z * z) + return Vector3(r * np.cos(theta), r * np.sin(theta), z) + + def random_hemisphere_direction(normal): + dir = random_unit_vector() + if dir.dot(normal) < 0: + dir = Vector3(-dir.x, -dir.y, -dir.z) + return dir + + def trace_ray(ray_origin, ray_direction): + closest_t = np.inf + hit_object = None + hit_info = None + for obj in objects: + result = obj.intersect(ray_origin, ray_direction) + if result: + t, hit_point, normal = result + if t < closest_t: + closest_t = t + hit_object = obj + hit_info = (hit_point, normal, obj.color) + if hit_object: + hit_point, normal, color = hit_info + direct_light = compute_direct_light(hit_point, normal) + indirect_light = estimate_radiance(hit_point, normal) + return color * (direct_light + indirect_light) + else: + return np.array([0, 0, 0]) + + def compute_direct_light(point, normal): + direction_to_light = (light_position - point).normalize() + shadow_origin = point + normal * 1e-5 + shadow_ray = direction_to_light + in_shadow = False + for obj in objects: + result = obj.intersect(shadow_origin, shadow_ray) + if result: + in_shadow = True + break + if in_shadow: + return np.array([0, 0, 0]) + else: + intensity = max(0, normal.dot(direction_to_light)) + return intensity * light_power / (4 * np.pi * (light_position - point).norm() ** 2) + + def estimate_radiance(point, normal): + accumulated_power = np.array([0.0, 0.0, 0.0]) + for photon_pos, photon_power in photon_map: + distance = (photon_pos - point).norm() + if distance < gather_radius: + weight = max(0, normal.dot((photon_pos - point).normalize())) + accumulated_power += photon_power * weight + area = np.pi * gather_radius ** 2 + return accumulated_power / area + + emit_photons() + + aspect_ratio = width / height + fov = np.pi / 3 + image = np.zeros((height, width, 3)) + for y in range(height): + for x in range(width): + px = (2 * (x + 0.5) / width - 1) * np.tan(fov / 2) * aspect_ratio + py = (1 - 2 * (y + 0.5) / height) * np.tan(fov / 2) + ray_origin = Vector3(0, 0, 0) + ray_direction = Vector3(px, py, -1).normalize() + color = trace_ray(ray_origin, ray_direction) + image[y, x, :] = np.clip(color, 0, 1) + + return image