A failed attempt at BSSRDF, but I completely can't figure it out. #1554

free7187 · 2025-04-05T12:43:42Z

free7187
Apr 5, 2025

I'm trying to implement BSSRDF in Mitsuba 3. However, for some reason, I'm getting results that I completely cannot understand.
Attempting to reproduce the sampling process for SSS based on the PBRT v3 framework, just for verification purposes, I haven't yet involved the ray tracing process. The current code is as follows:

class bssrdf_data_integrator(mi.SamplingIntegrator):
    def __init__(self, props):
        super().__init__(props)
        self.m_UnitToMM=10
        self.m_MMToUnit=1/self.m_UnitToMM
        self.r_max=1.7826
        self.r_s=17.0218

    def generate_sample_ray(self, si,rInmm,phi_dir,sampler):
        offsetRay=mi.Ray3f()
        offset_x  = rInmm * dr.cos(phi_dir) * self.m_MMToUnit
        offset_y  = rInmm * dr.sin(phi_dir)* self.m_MMToUnit
        l=dr.sqrt(self.r_max*self.r_max-rInmm*rInmm)* self.m_MMToUnit*2
        
        rand_num = sampler.next_1d()
        offsetRay.o=si.p+si.sh_frame.s*offset_x+si.sh_frame.t*offset_y +si.sh_frame.n* 0.5*l
        offsetRay.d=-si.sh_frame.n
        offsetRay.o[rand_num>0.5]=si.p+si.sh_frame.t*offset_x+si.sh_frame.n*offset_y +si.sh_frame.s* 0.5*l
        offsetRay.d[rand_num>0.5]=-si.sh_frame.s
        offsetRay.o[rand_num>0.75]=si.p+si.sh_frame.n*offset_x+si.sh_frame.s*offset_y +si.sh_frame.t* 0.5*l
        offsetRay.d[rand_num>0.75]=-si.sh_frame.t
        offsetRay.maxt=l
        return offsetRay,l

    def sample(self, 
               scene: mi.Scene, 
               sampler:mi.Sampler, 
               ray: mi.RayDifferential3f, 
               medium: mi.Medium = None,
               active: bool = True):
        
        L = mi.Spectrum(0.0)

        si = scene.ray_intersect(ray, active)        

        valid_ray = active & si.is_valid()
        active &= si.is_valid()

        rInMM, pdf_r= self.sample_position(sampler.next_1d(active)*0.99, self.r_s)
        rInMM = dr.clip(rInMM, 0.0, self.r_max)

        phi_dir  = sampler.next_1d(active) * 2 * dr.pi
        offsetRay,l = self.generate_sample_ray(si, rInMM,phi_dir,sampler)
        
        so1 = scene.ray_intersect(offsetRay, si.is_valid())
        active_so1 = si.is_valid() & so1.is_valid()
        
        ray2 = mi.Ray3f()
        ray2.o = so1.p
        ray2.d = offsetRay.d
        ray2.maxt = dr.maximum(0, l - dr.norm(so1.p - offsetRay.o))
        active2 = active_so1 & (ray2.maxt > 0)
        so2 = scene.ray_intersect(ray2, active2)
        active_so2 = active2 & so2.is_valid()
        
        ray3 = mi.Ray3f()
        ray3.o = so2.p
        ray3.d = offsetRay.d
        ray3.maxt = dr.maximum(0, l - dr.norm(so2.p - offsetRay.o))
        active3 = active_so2 & (ray3.maxt > 0)
        so3 = scene.ray_intersect(ray3, active3)
        active_so3 = active3 & so3.is_valid()
        
        ray4 = mi.Ray3f()
        ray4.o = so3.p
        ray4.d = offsetRay.d
        ray4.maxt = dr.maximum(0, l - dr.norm(so3.p - offsetRay.o))
        active4 = active_so3 & (ray4.maxt > 0)
        so4 = scene.ray_intersect(ray4, active4)
        active_so4 = active4 & so4.is_valid()
        
        found_count = mi.Int32(0.0)
        found_count[active_so1] += 1
        found_count[active_so2] += 1
        found_count[active_so3] += 1
        found_count[active_so4] += 1
       
        select_idx = dr.minimum(dr.floor(sampler.next_1d(active) * found_count), found_count - 1)
    
        so = si
        mask1 = active_so1& (select_idx< 1)
        so[mask1] = so1
        
        mask2 = active_so2 & (select_idx< 2) & (select_idx >= 1)
        so[mask2] = so2

        mask3 = active_so3 & (select_idx < 3) & (select_idx >= 2)
        so[mask3] = so3
        
        mask4 = active_so4 & (select_idx < 4) & (select_idx >= 3)
        so[mask4] = so4
        
        active_final = mask1 | mask2 | mask3 | mask4
        
        any_found = found_count > 0

        ds, emitter_val = scene.sample_emitter_direction(
            so, sampler.next_2d(active_final), True, active_final
        )
        wo = so.to_local(ds.d)
        rcpPdf = dr.select(any_found, 1.0 / ( dr.maximum(found_count, 1)), 0)

        L[active_final] += rcpPdf* emitter_val*mi.Frame3f.cos_theta(wo)

        return dr.select(True, L, 0.0), valid_ray, []
    
    def sample_position(self, u, s):
        u = 1 - u
        g = 1 + 4 * u * (2 * u + dr.sqrt(1 + 4 * u * u))
        n = dr.exp2(dr.log2(g) * (-1/ 3))
        p = g * n * n
        c = 1 + p + n
        x = 3 * dr.log(c / (4 * u))
        r = x / s
        pdf = self.Burley_Pdf(r, s)
        return r, pdf

    def Burley_Pdf(self,r, s):
        x = -s * r
        return s * 0.25* (dr.exp(x) + dr.exp(x / 3))

Because I'm not sure about the execution mechanism of multiple loops under cuda_ad_rgb and issues with custom eval() functions, I've implemented multiple collision tests along rays in this unusual way.
The XML is as follows

<scene version="3.0.0">
<!-- Defaults, these can be set via the command line: -Darg=value -->
	<default name="spp" value="32"/>
	<default name="resx" value="960"/>
	<default name="resy" value="540"/>
<!-- Camera and Rendering Parameters -->
	<integrator type="My_integrator"/>
	<sensor type="perspective">
		<string name="fov_axis" value="x"/>
		<float name="fov" value="39.597755"/>
		<float name="principal_point_offset_x" value="0.000000"/>
		<float name="principal_point_offset_y" value="-0.000000"/>
		<float name="near_clip" value="0.100000"/>
		<float name="far_clip" value="100.000000"/>
		<transform name="to_world">
			<rotate x="1" angle="-90.00000250447816"/>
			<rotate z="1" angle="90.00000250447816"/>
			<translate value="4.000000 0.000000 0.000000"/>
		</transform>
		<sampler type="independent">
			<integer name="sample_count" value="$spp"/>
		</sampler>
		<film type="hdrfilm">
			<integer name="width" value="$resx"/>
			<integer name="height" value="$resy"/>
			<rfilter type="tent"/>
            <string name="pixel_format" value="rgb"/>
            <string name="component_format" value="float32"/>
		</film>
	</sensor>
    <emitter type="directional">
        <vector name="direction" value="-1.0, 0.5, 0.2"/>
        <rgb name="irradiance" value="3.0"/>
    </emitter>
	<shape type="cube">
		<transform name="to_world">
			<scale x="0.4" y="0.4" z="0.4"/>
			<rotate angle="-10" x="1"/>
			<rotate angle="20" z="1"/>
		</transform>
	</shape>
</scene>

but strangely, the rendering result looks like this.

I haven't implemented the complete BSSRDF process yet, so I wasn't expecting good results, but this still doesn't explain the strange lines on the cube surface. Why is this happening? Is there any problem with my sampling setup? Or are there any features in DrJit that I haven't noticed?

Answered by free7187

Apr 8, 2025

It seems I've solved the problems:
1.si.sh_frame.n is not equivalent to si.n, and because of this issue, diagonal lines appeared on the cube surface when rendering.
2.

ray2 = mi.Ray3f(offsetRay); 
     ray2.o = so1.p; 
     ray2.maxt = dr.maximum(0, l - dr.norm(so1.p - offsetRay.o))

was incorrect as it didn't account for ray self-intersection; a small offset should be added.

I still have another question: since I don't fully understand how to implement nested loops in DrJit under cuda_ad_rgb, especially how while loops should be implemented, I'm currently implementing transmittance through multiple collisions. Could anyone explain how to implement while() loop logic in DrJit?

View full answer

free7187 · 2025-04-07T07:32:15Z

free7187
Apr 7, 2025
Author

Although I have completed the BSSRDF implementation, there are still strange diagonal lines visible on the cube in the final rendering results. Does anyone know what might be causing this issue?

0 replies

free7187 · 2025-04-08T07:38:40Z

free7187
Apr 8, 2025
Author

I've tried to further clarify my problem:
Testing has proven that si.sh_frame.n and si.n are not the same, here is my code now:

class bssrdf_data_integrator(mi.SamplingIntegrator):
    def __init__(self, props):
        super().__init__(props)
        self.m_UnitToMM=10
        self.m_MMToUnit=1/self.m_UnitToMM
        self.r_max=1.7826
        self.r_s=17.0218
        self.effctive_albedo=0.05
        self.value=new_config.value
    def generate_sample_ray(self, si,rInmm,phi_dir,sampler):
        offsetRay=mi.RayDifferential3f()
        offset_x  = rInmm * dr.cos(phi_dir) * self.m_MMToUnit
        offset_y  = rInmm * dr.sin(phi_dir)* self.m_MMToUnit
        l=dr.sqrt(self.r_max*self.r_max-rInmm*rInmm)* self.m_MMToUnit*2
        
        rand_num = sampler.next_1d()

        N=dr.normalize(si.n)
        T=dr.normalize(dr.cross(si.n, si.sh_frame.s))
        S=dr.normalize(dr.cross(T,si.n ))

        # offsetRay.o=si.p+si.sh_frame.s*offset_x+si.sh_frame.t*offset_y +si.sh_frame.n* 0.5*l
        # offsetRay.d=-si.sh_frame.n
        # offsetRay.o[rand_num>0.5]=si.p+si.sh_frame.t*offset_x+si.sh_frame.n*offset_y +si.sh_frame.s* 0.5*l
        # offsetRay.d[rand_num>0.5]=-si.sh_frame.s
        # offsetRay.o[rand_num>0.75]=si.p+si.sh_frame.n*offset_x+si.sh_frame.s*offset_y +si.sh_frame.t* 0.5*l
        # offsetRay.d[rand_num>0.75]=-si.sh_frame.t

        offsetRay.o=si.p+S*offset_x+T*offset_y -N* 0.5*l
        offsetRay.d=N
        offsetRay.o[rand_num>0.5]=si.p+T*offset_x+N*offset_y -S* 0.5*l
        offsetRay.d[rand_num>0.5]=S
        offsetRay.o[rand_num>0.75]=si.p+N*offset_x+S*offset_y -T* 0.5*l
        offsetRay.d[rand_num>0.75]=T
        offsetRay.maxt=l
        return offsetRay,l,S,T,N

    def sample(self, 
               scene: mi.Scene, 
               sampler:mi.Sampler, 
               ray: mi.RayDifferential3f, 
               medium: mi.Medium = None,
               active: bool = True):
        
        L = mi.Spectrum(0.0)
        
        si = scene.ray_intersect(ray, active)        
        valid_ray = active & si.is_valid()
        active &= si.is_valid()

        rInMM, pdf_r= self.sample_position(sampler.next_1d(active)*0.99, self.r_s)
        rInMM = dr.clip(rInMM, 0.0, self.r_max)
        phi_dir  = sampler.next_1d(active) * 2 * dr.pi
        offsetRay,l,S,T,N = self.generate_sample_ray(si, rInMM,phi_dir,sampler)
        
        so1 = scene.ray_intersect(offsetRay, si.is_valid())
        active_so1 = si.is_valid() & so1.is_valid()
        
        ray2 = mi.Ray3f(offsetRay)
        ray2.o = so1.p
        ray2.maxt = dr.maximum(0, l - dr.norm(so1.p - offsetRay.o))
        so2 = scene.ray_intersect(ray2, active)
        active_so2 = active_so1 & so2.is_valid()
        
        ray3 = mi.Ray3f(offsetRay)
        ray3.o = so2.p
        ray3.maxt = dr.maximum(0, l - dr.norm(so2.p - offsetRay.o))
        so3 = scene.ray_intersect(ray3, active_so2)
        active_so3 = active_so2 & so3.is_valid()
        
        ray4 = mi.Ray3f(offsetRay)
        ray4.o = so3.p
        ray4.maxt = dr.maximum(0, l - dr.norm(so3.p - offsetRay.o))
        so4 = scene.ray_intersect(ray4, active_so3)
        active_so4 = active_so3 & so4.is_valid()
        
        found_count = mi.UInt32(0.0)
        found_count_1 = mi.Spectrum(0.0)
        found_count[active_so1] += 1
        found_count_1[active_so1] += 1.0
        found_count[active_so2] += 1
        found_count_1[active_so2] += 1.0
        found_count[active_so3] += 1
        found_count_1[active_so3] += 1.0
        found_count[active_so4] += 1
        found_count_1[active_so4] += 1.0
    
        select_idx = dr.minimum(dr.floor(sampler.next_1d(active) * found_count), found_count - 1)

        so = so1
        mask1 = active_so1& (select_idx< 1)
        so = dr.select(mask1,so1,so) 
        
        mask2 = active_so2 & (select_idx< 2) & (select_idx >= 1)
        so = dr.select(mask2,so2,so) 

        mask3 = active_so3 & (select_idx < 3) & (select_idx >= 2)
        so = dr.select(mask3,so3,so) 
        
        mask4 = active_so4 & (select_idx < 4) & (select_idx >= 3)
        so = dr.select(mask4,so4,so) 

        active_final = mask1 | mask2 | mask3 | mask4
        
        any_found = found_count > 0
        pdf_r=self.pdf_sp(si, so,S,T,N)/found_count_1
        # pdf_r=self.pdf_sp(si, so)
        rcpPdf = dr.select((pdf_r > 0) & any_found, 1.0 / (pdf_r), 0)
        ds, emitter_val = scene.sample_emitter_direction(
            so, sampler.next_2d(active_final), True, active_final
        )
        wo = so.to_local(ds.d)
        cosThetaO = mi.Frame3f.cos_theta(wo)

        real_r=dr.norm(so.p-si.p)*self.m_UnitToMM

        R_1=self.Burley_bssrdf(real_r, self.r_s)
        L[active_final] +=  R_1*rcpPdf * emitter_val*cosThetaO
        return dr.select(True, L, 0.0), valid_ray, []
    
    def sample_position(self, u, s):
        u = 1 - u
        g = 1 + 4 * u * (2 * u + dr.sqrt(1 + 4 * u * u))
        n = dr.exp2(dr.log2(g) * (-1/ 3))
        p = g * n * n
        c = 1 + p + n
        x = 3 * dr.log(c / (4 * u))
        r = x / s
        pdf = self.Burley_Pdf(r, s)
        return r, pdf
    
    def pdf_sp(self, si,so,S,T,N):
        d = si.p - so.p  
        
        dLocal = mi.Vector3f(
            dr.dot(S, d),
            dr.dot(T, d),
            dr.dot(N, d)
        )
        
        nLocal = mi.Vector3f(
            dr.dot(S, so.n),
            dr.dot(T, so.n),
            dr.dot(N, so.n)
        )
        
        rProj_0 = dr.sqrt(dLocal.y * dLocal.y + dLocal.z * dLocal.z)
        rProj_1 = dr.sqrt(dLocal.z * dLocal.z + dLocal.x * dLocal.x)
        rProj_2 = dr.sqrt(dLocal.x * dLocal.x + dLocal.y * dLocal.y)
        
        axisProb = mi.Vector3f(0.25, 0.25, 0.5)
        
        pdf = mi.Spectrum(0.0)
        
        sr_pdf = self.Burley_Pdf(rProj_0*self.m_UnitToMM, self.r_s)
        pdf += sr_pdf * dr.abs(nLocal.x) * axisProb.x
        sr_pdf = self.Burley_Pdf(rProj_1*self.m_UnitToMM, self.r_s)
        pdf += sr_pdf * dr.abs(nLocal.y) * axisProb.y
        sr_pdf = self.Burley_Pdf(rProj_2*self.m_UnitToMM, self.r_s)
        pdf += sr_pdf * dr.abs(nLocal.z) * axisProb.z        
        return pdf    
     
    def Burley_bssrdf(self,r, s):
        x = -s * r
        return s * 0.25* (dr.exp(x) + dr.exp(x / 3))*self.effctive_albedo
    def Burley_Pdf(self,r, s):
        x = -s * r
        return s * 0.25* (dr.exp(x) + dr.exp(x / 3))

My further experimental results show:

It looks much better, but there are still some strange stripes. Does anyone know why this is happening? Or is there something inadequate in my description of the problem?

0 replies

free7187 · 2025-04-08T13:03:34Z

free7187
Apr 8, 2025
Author

It seems I've solved the problems:
1.si.sh_frame.n is not equivalent to si.n, and because of this issue, diagonal lines appeared on the cube surface when rendering.
2.

ray2 = mi.Ray3f(offsetRay); 
     ray2.o = so1.p; 
     ray2.maxt = dr.maximum(0, l - dr.norm(so1.p - offsetRay.o))

was incorrect as it didn't account for ray self-intersection; a small offset should be added.

I still have another question: since I don't fully understand how to implement nested loops in DrJit under cuda_ad_rgb, especially how while loops should be implemented, I'm currently implementing transmittance through multiple collisions. Could anyone explain how to implement while() loop logic in DrJit?

2 replies

merlinND Apr 8, 2025
Collaborator

Hello @free7187,

Regarding loops in DrJit, please check out the documentation on control flow: https://drjit.readthedocs.io/en/latest/cflow.html

wjakob Apr 8, 2025
Maintainer

@free7187 For this kind of advanced coding, reading the Dr.Jit documentation back-to-back would be advised ;)

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

A failed attempt at BSSRDF, but I completely can't figure it out. #1554

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Replies: 3 comments 2 replies

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

A failed attempt at BSSRDF, but I completely can't figure it out. #1554

Uh oh!

Uh oh!

free7187 Apr 5, 2025

Replies: 3 comments · 2 replies

Uh oh!

free7187 Apr 7, 2025 Author

Uh oh!

Uh oh!

free7187 Apr 8, 2025 Author

Uh oh!

Uh oh!

free7187 Apr 8, 2025 Author

Uh oh!

merlinND Apr 8, 2025 Collaborator

Uh oh!

wjakob Apr 8, 2025 Maintainer

free7187
Apr 5, 2025

Replies: 3 comments 2 replies

free7187
Apr 7, 2025
Author

free7187
Apr 8, 2025
Author

free7187
Apr 8, 2025
Author

merlinND Apr 8, 2025
Collaborator

wjakob Apr 8, 2025
Maintainer