get rid of offset allocation

Author: SimonDanisch

docs/code/caustic_glass.jl

@@ -1,98 +1,96 @@
 using GeometryBasics
 using LinearAlgebra
-using Trace
+using Trace, FileIO, MeshIO
 
-function render()
-    glass = Trace.GlassMaterial(
-        Trace.ConstantTexture(Trace.RGBSpectrum(1f0)),
-        Trace.ConstantTexture(Trace.RGBSpectrum(1f0)),
-        Trace.ConstantTexture(0f0),
-        Trace.ConstantTexture(0f0),
-        Trace.ConstantTexture(1.25f0),
-        true,
-    )
-    plastic = Trace.PlasticMaterial(
-        Trace.ConstantTexture(Trace.RGBSpectrum(0.6399999857f0, 0.6399999857f0, 0.6399999857f0)),
-        Trace.ConstantTexture(Trace.RGBSpectrum(0.1000000015f0, 0.1000000015f0, 0.1000000015f0)),
-        Trace.ConstantTexture(0.010408001f0),
-        true,
-    )
+glass = Trace.GlassMaterial(
+    Trace.ConstantTexture(Trace.RGBSpectrum(1f0)),
+    Trace.ConstantTexture(Trace.RGBSpectrum(1f0)),
+    Trace.ConstantTexture(0f0),
+    Trace.ConstantTexture(0f0),
+    Trace.ConstantTexture(1.25f0),
+    true,
+)
+plastic = Trace.PlasticMaterial(
+    Trace.ConstantTexture(Trace.RGBSpectrum(0.6399999857f0, 0.6399999857f0, 0.6399999857f0)),
+    Trace.ConstantTexture(Trace.RGBSpectrum(0.1000000015f0, 0.1000000015f0, 0.1000000015f0)),
+    Trace.ConstantTexture(0.010408001f0),
+    true,
+)
 
-    model = "./scenes/models/caustic-glass.ply"
-    triangle_meshes, triangles = Trace.load_triangle_mesh(
-        model, Trace.ShapeCore(Trace.translate(Vec3f(5, -1.49, -100)), false),
-    )
-    floor_triangles = Trace.create_triangle_mesh(
-        Trace.ShapeCore(Trace.translate(Vec3f(-10, 0, -87)), false),
-        2, UInt32[1, 2, 3, 1, 4, 3],
-        4,
-        [
-            Point3f(0, 0, 0), Point3f(0, 0, -30),
-            Point3f(30, 0, -30), Point3f(30, 0, 0),
-        ],
-        [
-            Trace.Normal3f(0, 1, 0), Trace.Normal3f(0, 1, 0),
-            Trace.Normal3f(0, 1, 0), Trace.Normal3f(0, 1, 0),
-        ],
-    )
+model = load(joinpath(@__DIR__, "..", "src", "assets",  "models", "caustic-glass.ply"))
 
-    primitives = Vector{Trace.GeometricPrimitive}(undef, 0)
-    for t in triangles
-        push!(primitives, Trace.GeometricPrimitive(t, glass))
-    end
-    for t in floor_triangles
-        push!(primitives, Trace.GeometricPrimitive(t, plastic))
-    end
 
-    bvh = Trace.BVHAccel(primitives, 1)
+triangles = Trace.create_triangle_mesh(
+    model, Trace.ShapeCore(Trace.translate(Vec3f(5, -1.49, -100)), false),
+)
 
-    from, to = Point3f(0, 2, 0), Point3f(-5, 0, 5)
-    cone_angle, cone_δ_angle = 30f0, 10f0
-    dir = normalize(Vec3f(to - from))
-    dir, du, dv = Trace.coordinate_system(dir)
+floor_triangles = Trace.create_triangle_mesh(
+    Trace.ShapeCore(Trace.translate(Vec3f(-10, 0, -87)), false),
+    UInt32[1, 2, 3, 1, 4, 3],
+    [
+        Point3f(0, 0, 0), Point3f(0, 0, -30),
+        Point3f(30, 0, -30), Point3f(30, 0, 0),
+    ],
+    [
+        Trace.Normal3f(0, 1, 0), Trace.Normal3f(0, 1, 0),
+        Trace.Normal3f(0, 1, 0), Trace.Normal3f(0, 1, 0),
+    ],
+)
 
-    dir_to_z = Trace.Transformation(transpose(Mat4f(
-        du[1], du[2], du[3], 0f0,
-        dv[1], dv[2], dv[3], 0f0,
-        dir[1], dir[2], dir[3], 0f0,
-        0f0, 0f0, 0f0, 1f0,
-    )))
-    light_to_world = (
-        Trace.translate(Vec3f(4.5, 0, -101))
-        * Trace.translate(Vec3f(from))
-        * inv(dir_to_z)
-    )
+primitives = Trace.GeometricPrimitive[]
+for t in triangles
+    push!(primitives, Trace.GeometricPrimitive(t, glass))
+end
+for t in floor_triangles
+    push!(primitives, Trace.GeometricPrimitive(t, plastic))
+end
 
-    lights = [
-        Trace.SpotLight(
-            light_to_world, Trace.RGBSpectrum(60f0),
-            cone_angle, cone_angle - cone_δ_angle,
-        ),
-    ]
+bvh = Trace.BVHAccel(primitives, 1);
 
-    scene = Trace.Scene(lights, bvh)
+from, to = Point3f(0, 2, 0), Point3f(-5, 0, 5)
+cone_angle, cone_δ_angle = 30f0, 10f0
+dir = normalize(Vec3f(to - from))
+dir, du, dv = Trace.coordinate_system(dir)
 
-    resolution = Point2f(256)
-    n_samples = 8
-    ray_depth = 5
+dir_to_z = Trace.Transformation(transpose(Mat4f(
+    du[1], du[2], du[3], 0f0,
+    dv[1], dv[2], dv[3], 0f0,
+    dir[1], dir[2], dir[3], 0f0,
+    0f0, 0f0, 0f0, 1f0,
+)))
+light_to_world = (
+    Trace.translate(Vec3f(4.5, 0, -101))
+    * Trace.translate(Vec3f(from))
+    * inv(dir_to_z)
+)
 
-    look_point = Point3f(-3, 0, -91)
-    screen = Trace.Bounds2(Point2f(-1f0), Point2f(1f0))
-    filter = Trace.LanczosSincFilter(Point2f(1f0), 3f0)
+lights = [
+    Trace.SpotLight(
+        light_to_world, Trace.RGBSpectrum(60f0),
+        cone_angle, cone_angle - cone_δ_angle,
+    ),
+]
 
-    ir = Int64.(resolution)
-    film = Trace.Film(
-        resolution, Trace.Bounds2(Point2f(0), Point2f(1)),
-        filter, 1f0, 1f0,
-        "./scenes/caustics-sppm-$(ir[1])x$(ir[2]).png",
-    )
-    camera = Trace.PerspectiveCamera(
-        Trace.look_at(Point3f(0, 150, 150), look_point, Vec3f(0, 1, 0)),
-        screen, 0f0, 1f0, 0f0, 1f6, 90f0, film,
-    )
+scene = Trace.Scene(lights, bvh)
 
-    integrator = Trace.SPPMIntegrator(camera, 0.075f0, ray_depth, 100, -1)
-    integrator(scene)
-end
+resolution = Point2f(1024)
+n_samples = 8
+ray_depth = 5
+
+look_point = Point3f(-3, 0, -91)
+screen = Trace.Bounds2(Point2f(-1f0), Point2f(1f0))
+filter = Trace.LanczosSincFilter(Point2f(1f0), 3f0)
+
+ir = Int64.(resolution)
+film = Trace.Film(
+    resolution, Trace.Bounds2(Point2f(0), Point2f(1)),
+    filter, 1f0, 1f0,
+    "./scenes/caustics-sppm-$(ir[1])x$(ir[2]).png",
+)
+camera = Trace.PerspectiveCamera(
+    Trace.look_at(Point3f(0, 150, 150), look_point, Vec3f(0, 1, 0)),
+    screen, 0f0, 1f0, 0f0, 1f6, 90f0, film,
+)
 
-render()
+integrator = Trace.SPPMIntegrator(camera, 0.075f0, ray_depth, 10)
+integrator(scene)

src/film.jl

@@ -99,21 +99,20 @@ struct FilmTile{Pixels<:AbstractMatrix{<:FilmTilePixel}}
     bounds::Bounds2
     filter_radius::Point2f
     inv_filter_radius::Point2f
-    filter_table::Matrix{Float32}
     filter_table_width::Int32
     pixels::Pixels
 
     function FilmTile(
             bounds::Bounds2, filter_radius::Point2f,
-            filter_table::Matrix{Float32}, filter_table_width::Int32,
+            filter_table_width::Int32,
         )
         tile_res = (Int32.(inclusive_sides(bounds)))
-        pixels = StructArray{FilmTilePixel{RGBSpectrum}}(undef, tile_res[2], tile_res[1])
-        pixels.contrib_sum .= (RGBSpectrum(),)
-        pixels.filter_weight_sum .= 0.0f0
+        contrib_sum = fill(RGBSpectrum(), tile_res[2], tile_res[1])
+        filter_weight_sum = fill(0.0f0, tile_res[2], tile_res[1])
+        pixels = StructArray{FilmTilePixel{RGBSpectrum}}((contrib_sum, filter_weight_sum))
         new{typeof(pixels)}(
             bounds, filter_radius, 1f0 ./ filter_radius,
-            filter_table, filter_table_width,
+            filter_table_width,
             pixels,
         )
     end
@@ -126,7 +125,30 @@ function FilmTile(f::Film, sample_bounds::Bounds2)
     p0 = ceil.(sample_bounds.p_min .- 0.5f0 .- f.filter.radius)
     p1 = floor.(sample_bounds.p_max .- 0.5f0 .+ f.filter.radius) .+ 1f0
     tile_bounds = Bounds2(p0, p1) ∩ f.crop_bounds
-    FilmTile(tile_bounds, f.filter.radius, f.filter_table, f.filter_table_width)
+    FilmTile(tile_bounds, f.filter.radius, f.filter_table_width)
+end
+
+function filter_offset!(offsets, start, stop, discrete_point, inv_filter_radius, filter_table_width)
+    @inbounds for (i, x) in enumerate(Int(start):Int(stop))
+        fx = abs((x - discrete_point[1]) * inv_filter_radius * filter_table_width)
+        offsets[i] = clamp(ceil(fx), 1, filter_table_width)  # TODO is clipping ok?
+    end
+end
+
+@inline function filter_offset(x, discrete_point, inv_filter_radius, filter_table_width)
+    fx = abs((x - discrete_point[1]) * inv_filter_radius * filter_table_width)
+    return clamp(ceil(Int32, fx), Int32(1), Int32(filter_table_width))  # TODO is clipping ok?
+end
+
+function filter_offsets(start, stop, discrete_point, inv_filter_radius, filter_table_width)::NTuple{8, Int32}
+    range = Int32(start):Int32(stop)
+    return ntuple(8) do i
+        if i <= length(range)
+            filter_offset(range[i], discrete_point, inv_filter_radius, filter_table_width)::Int32
+        else
+            Int32(0)
+        end
+    end
 end
 
 """
@@ -137,33 +159,26 @@ Add sample contribution to the film tile.
     And is relative to the film, not the film tile.
 """
 function add_sample!(
-        t::FilmTile, point::Point2f, spectrum::S,
+        film::Film, t::FilmTile, point::Point2f, spectrum::S,
         sample_weight::Float32 = 1f0,
     ) where S<:Spectrum
 
     # Compute sample's raster bounds.
     discrete_point = point .- 0.5f0
     p0 = ceil.(Int32, discrete_point .- t.filter_radius)
-    p1 = floor.(Int32, discrete_point .+ t.filter_radius) .+ 1
+    p1 = floor.(Int32, discrete_point .+ t.filter_radius) .+ Int32(1)
     p0 = Int32.(max.(p0, max.(t.bounds.p_min, Point2{Int32}(1))))
     p1 = Int32.(min.(p1, t.bounds.p_max))
     # Precompute x & y filter offsets.
-    offsets_x = Vector{Int32}(undef, p1[1] - p0[1] + 1)
-    offsets_y = Vector{Int32}(undef, p1[2] - p0[2] + 1)
-    @inbounds for (i, x) in enumerate(p0[1]:p1[1])
-        fx = abs((x - discrete_point[1]) * t.inv_filter_radius[1] * t.filter_table_width)
-        offsets_x[i] = clamp(ceil(fx), 1, t.filter_table_width)  # TODO is clipping ok?
-    end
-    @inbounds for (i, y) in enumerate(p0[2]:p1[2])
-        fy = abs((y - discrete_point[2]) * t.inv_filter_radius[2] * t.filter_table_width)
-        offsets_y[i] = clamp(floor(fy), 1, t.filter_table_width)
-    end
+    offsets_x = filter_offsets(p0[1], p1[1], discrete_point, t.inv_filter_radius[1], t.filter_table_width)
+    offsets_y = filter_offsets(p0[2], p1[2], discrete_point, t.inv_filter_radius[2], t.filter_table_width)
     # Loop over filter support & add sample to pixel array.
     pixels = t.pixels
     contrib_sum = pixels.contrib_sum
     filter_weight_sum = pixels.filter_weight_sum
-    @inbounds for (j, y) in enumerate(p0[2]:p1[2]), (i, x) in enumerate(p0[1]:p1[1])
-        w = t.filter_table[offsets_y[j], offsets_x[i]]
+    filter_table = film.filter_table
+    @inbounds for (j, y) in enumerate(Int(p0[2]):Int(p1[2])), (i, x) in enumerate(Int(p0[1]):Int(p1[1]))
+        w = filter_table[offsets_y[j], offsets_x[i]]
         @real_assert sample_weight <= 1
         @real_assert w <= 1
         idx = get_pixel_index(t, Point2(x, y))
@@ -196,7 +211,7 @@ function merge_film_tile!(f::Film, ft::FilmTile)
     f_xyz = f.pixels.xyz
     f_filter_weight_sum = f.pixels.filter_weight_sum
     @inbounds for y in y_range, x in x_range
-        pixel = Point2(x, y)
+        pixel = Point2{Int32}(x, y)
         ft_idx = get_pixel_index(ft, pixel)
         f_idx = get_pixel_index(f, pixel)
         f_xyz[f_idx] += to_XYZ(ft_contrib_sum[ft_idx])
@@ -212,21 +227,28 @@ function set_image!(f::Film, spectrum::Matrix{S}) where {S<:Spectrum}
 end
 
 function save(film::Film, splat_scale::Float32 = 1f0)
-    image = Array{Float32}(undef, size(film.pixels)..., 3)
-    for y in 1:size(film.pixels)[1], x in 1:size(film.pixels)[2]
-        pixel = film.pixels[y, x]
-        image[y, x, :] .= XYZ_to_RGB(pixel.xyz)
+    image = film.framebuffer
+    xyz = film.pixels.xyz
+    filter_weight_sum = film.pixels.filter_weight_sum
+    splat_xyz = film.pixels.splat_xyz
+
+    @inbounds for idx in eachindex(film.pixels)
+        rgb = XYZ_to_RGB(xyz[idx])
         # Normalize pixel with weight sum.
-        filter_weight_sum = pixel.filter_weight_sum
-        if filter_weight_sum != 0
-            inv_weight = 1f0 / filter_weight_sum
-            image[y, x, :] .= max.(0f0, image[y, x, :] .* inv_weight)
+        fws = filter_weight_sum[idx]
+        if fws != 0
+            inv_weight = 1.0f0 / fws
+            rgb = max.(0.0f0, rgb .* inv_weight)
         end
         # Add splat value at pixel & scale.
-        splat_rgb = XYZ_to_RGB(pixel.splat_xyz)
-        image[y, x, :] .+= splat_scale .* splat_rgb
-        image[y, x, :] .*= film.scale
+        splat_rgb = XYZ_to_RGB(splat_xyz[idx])
+        rgb = rgb .+ splat_scale .* splat_rgb
+        rgb = rgb .* film.scale
+        rgb = map(rgb) do c
+            c = ifelse(isfinite(c), c, 0.0f0)
+            return clamp(c, 0.0f0, 1.0f0)
+        end
+        image[idx] = RGB(rgb...)
     end
-    clamp!(image, 0f0, 1f0) # TODO remap instead of clamping?
-    FileIO.save(film.filename, image[end:-1:begin, :, :])
+    FileIO.save(film.filename, @view image[end:-1:begin, :])
 end

src/integrators/sampler.jl

@@ -11,6 +11,7 @@ end
 Render scene.
 """
 function (i::SamplerIntegrator)(scene::Scene)
+
     sample_bounds = get_sample_bounds(get_film(i.camera))
     sample_extent = diagonal(sample_bounds)
     tile_size = 16
@@ -23,6 +24,7 @@ function (i::SamplerIntegrator)(scene::Scene)
 
     @info "Utilizing $(Threads.nthreads()) threads"
     mempools = [MemoryPool(round(Int, 2*16384)) for _ in 1:Threads.maxthreadid()]
+    film = get_film(i.camera)
     Threads.@threads for k in 0:total_tiles
         x, y = k % width, k ÷ width
         tile = Point2f(x, y)
@@ -32,7 +34,7 @@ function (i::SamplerIntegrator)(scene::Scene)
         tb_max = min.(tb_min .+ (tile_size - 1), sample_bounds.p_max)
         tile_bounds = Bounds2(tb_min, tb_max)
 
-        film_tile = FilmTile(get_film(i.camera), tile_bounds)
+        film_tile = FilmTile(film, tile_bounds)
         spp_sqr = 1f0 / √Float32(t_sampler.samples_per_pixel)
         pool = mempools[Threads.threadid()]
         for pixel in tile_bounds
@@ -46,7 +48,7 @@ function (i::SamplerIntegrator)(scene::Scene)
                 ω > 0.0f0 && (l = li(pool, i, ray, scene, 1))
                 # TODO check l for invalid values
                 isnan(l) && (l = RGBSpectrum(0f0))
-                add_sample!(film_tile, camera_sample.film, l, ω)
+                add_sample!(film, film_tile, camera_sample.film, l, ω)
                 start_next_sample!(t_sampler)
             end
         end
@@ -108,6 +110,7 @@ function specular_reflect(
 ) where I<:SamplerIntegrator
 
     # Compute specular reflection direction `wi` and BSDF value.
+
     wo = surface_intersect.core.wo
     type = BSDF_REFLECTION | BSDF_SPECULAR
     wi, f, pdf, sampled_type = sample_f(

src/materials/bsdf.jl

@@ -61,6 +61,7 @@ function BSDF(si::SurfaceInteraction, sbdfs::Vararg{UberBxDF{S}, N}) where {S<:S
 end
 
 function BSDF(si::SurfaceInteraction, η::Float32, sbdfs::Vararg{UberBxDF{S},N}) where {S<:Spectrum, N}
+
     ng = si.core.n
     ns = si.shading.n
     ss = normalize(si.shading.∂p∂u)

src/primitive.jl

@@ -31,9 +31,9 @@ function compute_scattering!(
         pool, p::GeometricPrimitive, si::SurfaceInteraction,
         allow_multiple_lobes::Bool, transport::UInt8,
     )
+    @real_assert (si.core.n ⋅ si.shading.n) ≥ 0f0
     if p.material.type !== NO_MATERIAL
         return p.material(pool, si, allow_multiple_lobes, transport)
     end
-    # @real_assert (si.core.n ⋅ si.shading.n) ≥ 0f0
     return BSDF()
 end