BDPT debug

README.md

@@ -6,11 +6,11 @@
 
 This renderer is implemented based on **MY OWN** understanding of path tracing and other CG knowledge, therefore I **DO NOT** guarantee usability (also, I have done no verification experiments). The output results just... look decent:
 
-|         "The cornell spheres"          |         "The cornell boxes"         |
-| :------------------------------------: | :---------------------------------: |
-| ![](./assets/adapt-cornell-sphere.png) | ![](./assets/adapt-cornell-box.png) |
-|         "The cornell volume box"       |         "BDPT cbox 64 spp"         |
-| ![adapt-volume-box](https://user-images.githubusercontent.com/46109954/223172315-18888c22-3699-42ba-801d-bbd054e9246b.png) | ![cbox-64-bdpt](https://user-images.githubusercontent.com/46109954/223172423-bec7ac02-8533-432e-9bef-4f02bb4ddbb9.png) |
+|         "The cornell spheres"          |         "The cornell boxes"         | "Fresnel Blend" |
+| :------------------------------------: | :---------------------------------: | :------------------------------------: |
+| ![](./assets/adapt-cornell-sphere.png) | ![](./assets/adapt-cornell-box.png) | ![pbr-big-bdpt](https://user-images.githubusercontent.com/126778364/225679926-f75aab9f-0f47-4f45-ab4a-3ea7eaf34055.png)|
+|         "The cornell volume box"       |         "BDPT cbox 64 spp"         | "Giant mirror ball" |
+| ![pbr-cbox-bdpt](https://user-images.githubusercontent.com/126778364/225680094-8084c378-1533-4b74-871e-4524fff88f28.png)| ![cbox-64-bdpt](https://user-images.githubusercontent.com/46109954/223172423-bec7ac02-8533-432e-9bef-4f02bb4ddbb9.png) | ![pbr-single-ball-bdpt-single-ball](https://user-images.githubusercontent.com/126778364/225680022-ffeb3380-eeab-4beb-9bff-d3c631c36204.png)|
 
 
 
@@ -20,7 +20,7 @@ Here are the features I currently implemented and supports:
 - A unidirectional / bidirectional Monte-Carlo MIS path tracer: supports as many bounce times as you wish, and the rendering process is based on Taichi Lang, therefore it can be very fast (not on the first run, the first run of a scene might take a long time due to taichi function inlining, especially for BDPT). The figures displayed above can be rendered within 15-20s (with cuda-backend, GPU supported). The rendering result is displayed incrementally, or maximum iteration number can be pre-set.
 - Volumetric path tracer that supports uni/bidirectional path tracing in both bounded and unbounded condition
 - Global / indirect illumination & Ability to handle simple caustics
-- BRDFs: `Lambertian`, `Modified Phong` (Lafortune and Willems 1994), `Frensel Blend` (Ashikhmin and Shirley 2002), `Blinn-Phong`, `Mirror-specular`.
+- BRDFs: `Lambertian`, `Modified Phong` (Lafortune and Willems 1994), `Fresnel Blend` (Ashikhmin and Shirley 2002), `Blinn-Phong`, `Mirror-specular`.
 - BSDFs (with medium): deterministic refractive (glass-like)
 - mitusba-like XML scene file definition, supports mesh (from wavefront `.obj` file) and analytical sphere.
 - scene visualizer: which visualizes the scene you are going to render, helping to set some parameters like the relative position and camera pose

bxdf/brdf.py

@@ -1,6 +1,6 @@
 """
     All the BRDFs are here, note that only three kinds of simple BRDF are supported
-    Blinn-Phong / Lambertian / Mirror specular / Modified Phong / Frensel Blend
+    Blinn-Phong / Lambertian / Mirror specular / Modified Phong / Fresnel Blend
     @author: Qianyue He
     @date: 2023-1-23
 """
@@ -17,7 +17,7 @@
 from la.cam_transform import *
 from sampler.general_sampling import *
 from scene.general_parser import rgb_parse
-from renderer.constants import TRANSPORT_RAD, INV_PI
+from renderer.constants import INV_PI, ZERO_V3
 
 __all__ = ['BRDF_np', 'BRDF']
 
@@ -34,7 +34,7 @@ class BRDF_np:
     __all_specular_name     = {"specular", "k_s"}
     __all_absorption_name   = {"absorptions", "k_a"}
     # Attention: microfacet support will not be added recently
-    __type_mapping          = {"blinn-phong": 0, "lambertian": 1, "specular": 2, "microfacet": 3, "mod-phong": 4, "frensel-blend": 5}
+    __type_mapping          = {"blinn-phong": 0, "lambertian": 1, "specular": 2, "microfacet": 3, "mod-phong": 4, "fresnel-blend": 5}
 
     def __init__(self, elem: xet.Element, no_setup = False):
         self.type: str = elem.get("type")
@@ -74,9 +74,8 @@ def setup(self):
             raise NotImplementedError(f"Unknown BRDF type: {self.type}")
         self.type_id = BRDF_np.__type_mapping[self.type]
         if self.type_id == 2:
-            if self.k_g.max() < 1e-4:       # glossiness (actually means roughness) in specular BRDF being too "small"
-                self.is_delta = True
-        elif self.type_id == 5:             # precomputed coefficient for Frensel Blend BRDF
+            self.is_delta = True
+        elif self.type_id == 5:             # precomputed coefficient for Fresnel Blend BRDF
             self.k_g[2] = np.sqrt((self.k_g[0] + 1) * (self.k_g[1] + 1)) / (8. * np.pi)
 
     def export(self):
@@ -165,28 +164,29 @@ def sample_mod_phong(self, incid: vec3, normal: vec3):
         # Sample around reflected view dir (while blinn-phong samples around normal)
         return ray_out_d, spec, pdf
 
-    # ======================= Frensel-Blend =======================
+    # ======================= Fresnel-Blend =======================
     """
-        For Frensel Blend (by Ashikhmin and Shirley 2002), n_u and n_v will be stored in k_g
+        For Fresnel Blend (by Ashikhmin and Shirley 2002), n_u and n_v will be stored in k_g
         since k_g will not be used, k_d and k_s preserve their original meaning
     """
 
     @ti.func
-    def frensel_blend_dir(self, incid: vec3, half: vec3, normal: vec3, power_coeff: float):
+    def fresnel_blend_dir(self, incid: vec3, half: vec3, normal: vec3, power_coeff: float):
         reflected, dot_incid = inci_reflect_dir(incid, half)
         half_pdf = self.k_g[2] * tm.pow(tm.dot(half, normal), power_coeff)
         pdf = half_pdf / ti.max(ti.abs(dot_incid), EPS)
+
         valid_sample = tm.dot(normal, reflected) > 0.
         return reflected, pdf, valid_sample
 
     @ti.func
-    def frensel_cos2_sin2(self, half_vec: vec3, normal: vec3, R: mat3, dot_half: float):
+    def fresnel_cos2_sin2(self, half_vec: vec3, normal: vec3, R: mat3, dot_half: float):
         transed_x = (R @ vec3([1, 0, 0])).normalized()
         cos_phi2  = tm.dot(transed_x, (half_vec - dot_half * normal).normalized()) ** 2       # azimuth angle of half vector 
         return cos_phi2, 1. - cos_phi2
 
     @ti.func
-    def eval_frensel_blend(self, ray_in: vec3, ray_out: vec3, normal: vec3, R: mat3):
+    def eval_fresnel_blend(self, ray_in: vec3, ray_out: vec3, normal: vec3, R: mat3):
         # specular part, note that ray out is actually incident light in forward tracing
         half_vec = (ray_out - ray_in)
         dot_out  = tm.dot(normal, ray_out)
@@ -196,11 +196,11 @@ def eval_frensel_blend(self, ray_in: vec3, ray_out: vec3, normal: vec3, R: mat3)
             dot_in   = -tm.dot(normal, ray_in)              # incident dot should always be positive (otherwise it won't hit this point)
             dot_half = ti.abs(tm.dot(normal, half_vec))
             dot_hk   = ti.abs(tm.dot(half_vec, ray_out))
-            frensel  = schlick_frensel(self.k_s, dot_hk)
-            cos_phi2, sin_phi2 = self.frensel_cos2_sin2(half_vec, normal, R, dot_half)
+            fresnel  = schlick_fresnel(self.k_s, dot_hk)
+            cos_phi2, sin_phi2 = self.fresnel_cos2_sin2(half_vec, normal, R, dot_half)
             # k_g[2] should store sqrt((n_u + 1)(n_v + 1)) / 8pi
             denom = dot_hk * tm.max(dot_in, dot_out)
-            specular = self.k_g[2] * tm.pow(dot_half, self.k_g[0] * cos_phi2 + self.k_g[1] * sin_phi2) * frensel / denom
+            specular = self.k_g[2] * tm.pow(dot_half, self.k_g[0] * cos_phi2 + self.k_g[1] * sin_phi2) * fresnel / denom
             # diffusive part
             diffuse  = 28. / (23. * tm.pi) * self.k_d * (1. - self.k_s)
             pow5_in  = tm.pow(1. - dot_in / 2., 5)
@@ -210,13 +210,14 @@ def eval_frensel_blend(self, ray_in: vec3, ray_out: vec3, normal: vec3, R: mat3)
         return spec
 
     @ti.func
-    def sample_frensel_blend(self, incid: vec3, normal: vec3):
-        local_new_dir, power_coeff = frensel_hemisphere(self.k_g[0], self.k_g[1])
+    def sample_fresnel_blend(self, incid: vec3, normal: vec3):
+        local_new_dir, power_coeff = fresnel_hemisphere(self.k_g[0], self.k_g[1])
         ray_half, R = delocalize_rotate(normal, local_new_dir)
-        ray_out_d, pdf, is_valid = self.frensel_blend_dir(incid, ray_half, normal, power_coeff)
-        spec = vec3([0, 0, 0])
-        if is_valid:
-            spec = self.eval_frensel_blend(incid, ray_out_d, normal, R)
+        ray_out_d, pdf, is_valid = self.fresnel_blend_dir(incid, ray_half, normal, power_coeff)
+        if ti.random(float) > 0.5:
+            ray_out_d, _s, _p = self.sample_lambertian(normal)
+        pdf = 0.5 * (pdf + ti.abs(tm.dot(ray_out_d, normal)) * INV_PI)
+        spec = ti.select(is_valid, self.eval_fresnel_blend(incid, ray_out_d, normal, R), ZERO_V3)
         return ray_out_d, spec, pdf
 
     # ======================= Lambertian ========================
@@ -233,15 +234,6 @@ def sample_lambertian(self, normal: vec3):
         return ray_out_d, spec, pdf
 
     # ======================= Mirror-Specular ========================
-    @ti.func
-    def eval_specular(self, ray_in: vec3, ray_out: vec3, normal: vec3):
-        """ Attention: ray_in (in backward tracing) is actually out-going direction (in forward tracing) """
-        reflect_dir, _ = inci_reflect_dir(ray_in, normal)
-        spec = vec3([0, 0, 0])
-        if tm.dot(ray_out, reflect_dir) > 1 - 1e-4:
-            spec = self.k_d
-        return spec
-
     @ti.func
     def sample_specular(self, ray_in: vec3, normal: vec3):
         ray_out_d, _ = inci_reflect_dir(ray_in, normal)
@@ -261,15 +253,11 @@ def eval(self, incid: vec3, out: vec3, normal: vec3) -> vec3:
                 ret_spec = self.eval_blinn_phong(incid, out, normal)
             elif self._type == 1:       # Lambertian
                 ret_spec = self.eval_lambertian(out, normal)
-            elif self._type == 2:       # Specular
-                ret_spec = self.eval_specular(incid, out, normal)
             elif self._type == 4:
                 ret_spec = self.eval_mod_phong(incid, out, normal)
             elif self._type == 5:
                 R = rotation_between(vec3([0, 1, 0]), normal)
-                ret_spec = self.eval_frensel_blend(incid, out, normal, R)
-            else:
-                print(f"Warnning: unknown or unsupported BRDF type: {self._type} during evaluation.")
+                ret_spec = self.eval_fresnel_blend(incid, out, normal, R)
         return ret_spec
 
     @ti.func
@@ -290,8 +278,8 @@ def sample_new_rays(self, incid: vec3, normal: vec3):
             ret_dir, ret_spec, pdf = self.sample_specular(incid, normal)
         elif self._type == 4:       # Modified-Phong
             ret_dir, ret_spec, pdf = self.sample_mod_phong(incid, normal)
-        elif self._type == 5:       # Frensel-Blend
-            ret_dir, ret_spec, pdf = self.sample_frensel_blend(incid, normal)
+        elif self._type == 5:       # Fresnel-Blend
+            ret_dir, ret_spec, pdf = self.sample_fresnel_blend(incid, normal)
         else:
             print(f"Warnning: unknown or unsupported BRDF type: {self._type} during sampling.")
         return ret_dir, ret_spec, pdf
@@ -311,10 +299,6 @@ def get_pdf(self, outdir: vec3, normal: vec3, incid: vec3):
                 pdf = dot_outdir * INV_PI       # dot is cosine term
             elif self._type == 1:
                 pdf = dot_outdir * INV_PI
-            elif self._type == 2:
-                reflect_view, _ = inci_reflect_dir(incid, normal)
-                if tm.dot(reflect_view, outdir) > 1 - 1e-4:
-                    pdf = 1.0
             elif self._type == 4:
                 glossiness      = self.mean[2]
                 reflect_view, _ = inci_reflect_dir(incid, normal)
@@ -326,6 +310,7 @@ def get_pdf(self, outdir: vec3, normal: vec3, incid: vec3):
                 half_vec = (outdir - incid).normalized()
                 dot_half = tm.dot(half_vec, normal)
                 R = rotation_between(vec3([0, 1, 0]), normal)
-                cos_phi2, sin_phi2 = self.frensel_cos2_sin2(half_vec, normal, R, dot_half)
-                pdf = self.k_g[2] * tm.pow(dot_half, self.k_g[0] * cos_phi2 + self.k_g[1] * sin_phi2) * 4.
+                cos_phi2, sin_phi2 = self.fresnel_cos2_sin2(half_vec, normal, R, dot_half)
+                pdf = self.k_g[2] * tm.pow(dot_half, self.k_g[0] * cos_phi2 + self.k_g[1] * sin_phi2) / ti.abs(tm.dot(incid, half_vec))
+                pdf = 0.5 * (pdf + dot_outdir * INV_PI)
         return pdf

bxdf/bsdf.py

@@ -38,6 +38,8 @@ def __init__(self, elem: xet.Element):
         self.medium = Medium_np(elem.find("medium"))
         self.is_delta = False
         self.setup()
+        if self.type_id == 0:
+            self.is_delta = True
 
         # for BSDF, there will be medium defined in it
 
@@ -87,22 +89,19 @@ def sample_det_refraction(self, incid: vec3, normal: vec3, medium, mode):
         ret_dir = vec3([0, 1, 0])
         ret_int = self.k_d
         if is_total_reflection(dot_normal, ni, nr):
-            # Only sample BSDF reflection
-            ret_dir = (incid - 2 * normal * dot_normal).normalized()
+            # ret_dir = (incid - 2 * normal * dot_normal).normalized()      # we can account for total reflection... or not
+            ret_int.fill(0.)
         else:
-            refra_vec, valid_ref = snell_refraction(incid, normal, dot_normal, ni, nr)
-            if valid_ref:
-                reflect_ratio = frensel_equation(ni, nr, ti.abs(dot_normal), ti.abs(tm.dot(refra_vec, normal)))
-                if ti.random(float) > reflect_ratio:        # refraction
-                    ret_pdf = 1. - reflect_ratio
-                    ret_dir = refra_vec
-                    if mode == TRANSPORT_RAD:       # non-symmetry effect
-                        ret_int *= (ni * ni) / (nr * nr)
-                else:                                       # reflection
-                    ret_dir = (incid - 2 * normal * dot_normal).normalized()
-                    ret_pdf = reflect_ratio
-            else:
+            refra_vec, _v = snell_refraction(incid, normal, dot_normal, ni, nr)
+            reflect_ratio = fresnel_equation(ni, nr, ti.abs(dot_normal), ti.abs(tm.dot(refra_vec, normal)))
+            if ti.random(float) > reflect_ratio:        # refraction
+                ret_pdf = 1. - reflect_ratio
+                ret_dir = refra_vec
+                if mode == TRANSPORT_RAD:       # non-symmetry effect
+                    ret_int *= (ni * ni) / (nr * nr)
+            else:                                       # reflection
                 ret_dir = (incid - 2 * normal * dot_normal).normalized()
+                ret_pdf = reflect_ratio
         return ret_dir, ret_int * ret_pdf, ret_pdf
 
     @ti.func
@@ -115,32 +114,31 @@ def eval_det_refraction(self, ray_in: vec3, ray_out: vec3, normal: vec3, medium,
         ret_int = vec3([0, 0, 0])
         if is_total_reflection(dot_out, ni, nr):
             ref_dir = (ray_out - 2 * normal * dot_out).normalized()
-            if tm.dot(ref_dir, ray_in) > 1 - 1e-4:
+            if tm.dot(ref_dir, ray_in) > 1 - 5e-5:
                 ret_int = self.k_d
         else:
             # in sampling: ray_in points to the intersection, here ray_out points away from the intersection
             ref_dir = (ray_out - 2 * normal * dot_out).normalized()
             refra_vec, valid_ref = snell_refraction(ray_out, normal, dot_out, ni, nr)
             if valid_ref:
-                reflect_ratio = frensel_equation(ni, nr, ti.abs(dot_out), ti.abs(tm.dot(refra_vec, normal)))
-                if tm.dot(refra_vec, ray_in) > 1 - 1e-4:            # ray_in close to refracted dir
+                reflect_ratio = fresnel_equation(ni, nr, ti.abs(dot_out), ti.abs(tm.dot(refra_vec, normal)))
+                if tm.dot(refra_vec, ray_in) > 1 - 5e-5:            # ray_in close to refracted dir
                     ret_int = self.k_d * (1. - reflect_ratio)
                     if mode == TRANSPORT_RAD:    # consider non-symmetric effect due to different transport mode
                         ret_int *= (ni * ni) / (nr * nr)
-                elif tm.dot(ref_dir, ray_in) > 1 - 1e-4:            # ray_in close to reflected dir
+                elif tm.dot(ref_dir, ray_in) > 1 - 5e-5:            # ray_in close to reflected dir
                     ret_int = self.k_d * reflect_ratio
             else:
-                if tm.dot(ref_dir, ray_in) > 1 - 1e-4:            # ray_in close to reflected dir
+                if tm.dot(ref_dir, ray_in) > 1 - 5e-5:            # ray_in close to reflected dir
                     ret_int = self.k_d
         return ret_int
     # ========================= General operations =========================
 
     @ti.func
     def get_pdf(self, outdir: vec3, normal: vec3, incid: vec3, medium):
-        """ TODO: currently, deterministic refraction / Null transmission yields PDF = 0.0 """
         pdf = 0.
         if self._type == -1:
-            pdf = ti.select(tm.dot(incid, outdir) > 1 - 1e-4, 1., 0.)
+            pdf = ti.select(tm.dot(incid, outdir) > 1 - 5e-5, 1., 0.)
         elif self._type == 0:
             dot_out = tm.dot(outdir, normal)
             entering_this = dot_out < 0
@@ -150,13 +148,13 @@ def get_pdf(self, outdir: vec3, normal: vec3, incid: vec3, medium):
             ref_dir = (outdir - 2 * normal * dot_out).normalized()
             refra_vec, valid_refra = snell_refraction(outdir, normal, dot_out, ni, nr)
             if valid_refra:             # not total reflection, so there is not only one possible choice
-                reflect_ratio = frensel_equation(ni, nr, ti.abs(dot_out), ti.abs(tm.dot(refra_vec, normal)))
-                if tm.dot(refra_vec, incid) > 1 - 1e-4:            # ray_in close to refracted dir
+                reflect_ratio = fresnel_equation(ni, nr, ti.abs(dot_out), ti.abs(tm.dot(refra_vec, normal)))
+                if tm.dot(refra_vec, incid) > 1 - 5e-5:            # ray_in close to refracted dir
                     pdf = 1. - reflect_ratio
-                elif tm.dot(ref_dir, incid) > 1 - 1e-4:            # ray_in close to reflected dir
+                elif tm.dot(ref_dir, incid) > 1 - 5e-5:            # ray_in close to reflected dir
                     pdf = reflect_ratio
             else:
-                if tm.dot(ref_dir, incid) > 1 - 1e-4:
+                if tm.dot(ref_dir, incid) > 1 - 5e-5:
                     pdf = 1.
         return pdf
 

emitters/point.py

@@ -24,5 +24,5 @@ def __init__(self, elem: xet.Element = None):
         self.pos: np.ndarray = vec3d_parse(pos_elem)
 
     def export(self) -> TaichiSource:
-        bool_bits = 0x01 | (self.in_free_space << 4)
+        bool_bits = 0x21 + (self.in_free_space << 4)        # position delta (0x01), delta vertex (0x20)
         return TaichiSource(_type = 0, intensity = vec3(self.intensity), pos = vec3(self.pos), bool_bits = bool_bits)

inputs/csphere/balls-mono.xml

@@ -28,7 +28,7 @@
 		</film>
 	</sensor>
 
-	<brdf type="frensel-blend" id="frensel">
+	<brdf type="fresnel-blend" id="fresnel">
 		<rgb name="k_d" value="#CACACA"/>
 		<!-- <rgb name="k_g" value="1.0"/>
 		<rgb name="k_s" value="0.0"/> -->
@@ -140,7 +140,7 @@
 	<shape type="sphere">
 		<point name="center" x="4.2" y="0.5" z="4.1"/>
 		<float name="radius" value="0.5"/>
-		<ref type="material" id="frensel"/>
+		<ref type="material" id="fresnel"/>
 	</shape>
 
 	<shape type="sphere">

inputs/csphere/balls-multi.xml

@@ -28,7 +28,7 @@
 		</film>
 	</sensor>
 
-	<brdf type="frensel-blend" id="frensel">
+	<brdf type="fresnel-blend" id="fresnel">
 		<rgb name="k_d" value="#CACACA"/>
 		<!-- <rgb name="k_g" value="1.0"/>
 		<rgb name="k_s" value="0.0"/> -->
@@ -151,7 +151,7 @@
 	<shape type="sphere">
 		<point name="center" x="4.2" y="0.5" z="4.1"/>
 		<float name="radius" value="0.5"/>
-		<ref type="material" id="frensel"/>
+		<ref type="material" id="fresnel"/>
 	</shape>
 
 	<shape type="sphere">