a few improvements

Author: SimonDanisch

docs/code/TraceMakie.jl

@@ -42,7 +42,7 @@ function extract_material(plot::Plot, color_obs::Observable)
             return
         end
     elseif color isa Colorant || color isa Union{String,Symbol}
-        tex = Trace.ConstantTexture(to_spectrum(color))
+        tex = Trace.ConstantTexture(to_spectrum(to_color(color)))
     elseif color isa Nothing
         # ignore!
         nothing
@@ -193,9 +193,6 @@ function render_gpu(mscene::Makie.Scene, ArrayType; samples_per_pixel=8, max_dep
     return Array(gpu_img)
 end
 
-function render(w::Whitten5, scene)
-
-end
 
 function render_interactive(mscene::Makie.Scene, ArrayType; max_depth=5)
     scene, camera, film = convert_scene(mscene)

src/film.jl

@@ -4,6 +4,45 @@ struct Pixel
     splat_xyz::Point3f
 end
 Pixel() = Pixel(Point3f(0.0f0), 0.0f0, Point3f(0.0f0))
+
+
+function filter_offset(x, discrete_point, inv_filter_radius, filter_table_width)
+    fx = abs((x - discrete_point) * inv_filter_radius * filter_table_width)
+    return clamp(u_int32(ceil(fx)), Int32(1), Int32(filter_table_width))  # TODO is clipping ok?
+end
+
+
+function generate_filter_table(filter)
+    filter_table_width = 16
+    filter_table = Matrix{Float32}(undef, filter_table_width, filter_table_width)
+    r = filter.radius ./ filter_table_width
+    for y in 0:filter_table_width-1, x in 0:filter_table_width-1
+        p = Point2f((x + 0.5f0) * r[1], (y + 0.5f0) * r[2])
+        filter_table[y+1, x+1] = filter(p)
+    end
+
+    point = Point2f(filter_table_width)
+    # Compute sample's raster bounds.
+    discrete_point = point .- 0.5f0
+    # Compute sample radius around point
+    p0 = ceil.(Int, discrete_point .- filter.radius)
+    p1 = floor.(Int, discrete_point .+ filter.radius) .+ 1
+    # Make sure we're inbounds
+    inv_radius = 1.0f0 ./ filter.radius
+    # Precompute x & y filter offsets.
+    offsets_x = filter_offsets(p0[1], p1[1], discrete_point[1], inv_radius[1], filter_table_width)
+    offsets_y = filter_offsets(p0[2], p1[2], discrete_point[2], inv_radius[2], filter_table_width)
+    # Loop over filter support & add sample to pixel array.
+    xrange = p0[1]:p1[1]
+    yrange = p0[2]:p1[2]
+    weights = zeros(Float32, length(xrange), length(yrange))
+    for i in 1:length(xrange), j in 1:length(yrange)
+        w = filter_table[offsets_y[j], offsets_x[i]]
+        weights[i, j] = w
+    end
+    return weights
+end
+
 struct Film{Pixels<:AbstractMatrix{Pixel}}
     resolution::Point2f
     """
@@ -92,153 +131,68 @@ struct FilmTilePixel{S<:Spectrum}
 end
 FilmTilePixel() = FilmTilePixel(RGBSpectrum(), 0f0)
 
-struct FilmTile{Pixels<:AbstractMatrix{<:FilmTilePixel}}
-    """
-    Bounds should start from 1 not 0.
-    """
-    bounds::Bounds2
-    filter_radius::Point2f
-    inv_filter_radius::Point2f
-    filter_table_width::Int32
-    pixels::Pixels
-end
 
-function FilmTile(
-        bounds::Bounds2, filter_radius::Point2f,
-        filter_table_width::Int32,
-    )
-    # Include some padding for over rounding (since we re-use the tiles)
-    tile_res = (Int32.(inclusive_sides(bounds))) .+ 2
-    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))
-    FilmTile{typeof(pixels)}(
-        bounds, filter_radius, 1.0f0 ./ filter_radius,
-        filter_table_width,
-        pixels,
-    )
-end
 
 """
-Bounds should start from 1 not 0.
+Point in (x, y) format.
 """
-function FilmTile(f::Film, sample_bounds::Bounds2, radius)
-    p0 = ceil.(sample_bounds.p_min .- 0.5f0 .- radius)
-    p1 = floor.(sample_bounds.p_max .- 0.5f0 .+ radius) .+ 1f0
-    tile_bounds = Bounds2(p0, p1) ∩ f.crop_bounds
-    FilmTile(tile_bounds, radius, f.filter_table_width)
-end
-
-function reset!(tile::FilmTile)
-
-    tile.pixels.contrib_sum .= (RGBSpectrum(0f0),)
-    tile.pixels.filter_weight_sum .= 0f0
-end
-
-function update_bounds!(f::Film, tile::FilmTile, sample_bounds::Bounds2)
-    reset!(tile)
-    radius = tile.filter_radius
-    p0 = ceil.(sample_bounds.p_min .- 0.5f0 .- radius)
-    p1 = floor.(sample_bounds.p_max .- 0.5f0 .+ radius) .+ 1.0f0
-    bounds = Bounds2(p0, p1) ∩ f.crop_bounds
-    tile_res = (Int32.(inclusive_sides(bounds)))
-    @assert all(reverse(tile_res) .<= size(tile.pixels)) "$(reverse(tile_res)) != $(size(tile.pixels)) $(sample_bounds)"
-    FilmTile(bounds, radius, tile.inv_filter_radius, tile.filter_table_width, tile.pixels)
-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) * inv_filter_radius * filter_table_width)
-        offsets[i] = clamp(ceil(fx), 1, filter_table_width)  # TODO is clipping ok?
+@inline function get_pixel_index(crop_bounds, p::Point2)
+    i1, i2 = u_int32.((p .- crop_bounds.p_min .+ 1.0f0))
+    return CartesianIndex(i1, i2)
+end
+
+@inline function merge_film_tile!(f::AbstractMatrix{Pixel}, crop_bounds::Bounds2, ft::AbstractMatrix{FilmTilePixel}, tile::Bounds2, tile_col::Int32)
+    ft_contrib_sum = ft.contrib_sum
+    ft_filter_weight_sum = ft.filter_weight_sum
+    f_xyz = f.xyz
+    f_filter_weight_sum = f.filter_weight_sum
+    linear = Int32(1)
+    @inbounds for pixel in tile
+        f_idx = get_pixel_index(crop_bounds, pixel)
+        f_xyz[f_idx] += to_XYZ(ft_contrib_sum[linear, tile_col])
+        f_filter_weight_sum[f_idx] += ft_filter_weight_sum[linear, tile_col]
+        linear += Int32(1)
     end
+    return
 end
 
-@inline function filter_offset(x, discrete_point, inv_filter_radius, filter_table_width)
-    fx = abs((x - discrete_point) * inv_filter_radius * filter_table_width)
-    return clamp(u_int32(ceil(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
+@inline function get_tile_index(bounds::Bounds2, p::Point2)
+    j, i = u_int32.((p .- bounds.p_min .+ 1.0f0))
+    ncols = u_int32(inclusive_sides(bounds)[1])
+    return (i - Int32(1)) * ncols + j
 end
 
-"""
-Add sample contribution to the film tile.
-
-- `point::Point2f`:
-    should start from 1 not 0.
-    And is relative to the film, not the film tile.
-"""
-function add_sample!(
-        film::Film, t::FilmTile, point::Point2f, spectrum::S,
-        sample_weight::Float32 = 1f0,
-    ) where S<:Spectrum
-
+@inline function add_sample!(
+    tiles::AbstractMatrix{FilmTilePixel}, tile::Bounds2, tile_column::Int32, point::Point2f, spectrum::RGBSpectrum,
+    filter_table, filter_radius::Point2f, sample_weight::Float32=1.0f0,
+)
     # 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) .+ 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 = filter_offsets(p0[1], p1[1], discrete_point[1], t.inv_filter_radius[1], t.filter_table_width)
-    offsets_y = filter_offsets(p0[2], p1[2], discrete_point[2], t.inv_filter_radius[2], t.filter_table_width)
+    # Compute sample radius around point
+    p0 = u_int32.(ceil.(discrete_point .- filter_radius))
+    p1 = u_int32.(floor.(discrete_point .+ filter_radius)) .+ Int32(1)
+    # Make sure we're inbounds
+    pmin = u_int32.(tile.p_min)
+    pmax = u_int32.(tile.p_max)
+    p0 = max.(p0, max.(pmin, Point2{Int32}(1)))::Point2{Int32}
+    p1 = min.(p1, pmax)::Point2{Int32}
     # Loop over filter support & add sample to pixel array.
-    pixels = t.pixels
-    contrib_sum = pixels.contrib_sum
-    filter_weight_sum = pixels.filter_weight_sum
-    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))
-        contrib_sum[idx] += spectrum * sample_weight * w
-        filter_weight_sum[idx] += w
+    contrib_sum = tiles.contrib_sum
+    filter_weight_sum = tiles.filter_weight_sum
+    xrange = p0[1]:p1[1]
+    yrange = p0[2]:p1[2]
+    xn = length(xrange) % Int32
+    yn = length(yrange) % Int32
+    @inbounds for i in Int32(1):xn, j in Int32(1):yn
+        x = xrange[i]
+        y = yrange[j]
+        w = filter_table[i, j]
+        idx = get_tile_index(tile, Point2(x, y))
+        contrib_sum[idx, tile_column] += spectrum * sample_weight * w
+        filter_weight_sum[idx, tile_column] += w
     end
 end
 
-"""
-Point in (x, y) format.
-"""
-@inline function get_pixel_index(t::FilmTile, p::Point2)
-    i1, i2 = u_int32.((p .- t.bounds.p_min .+ 1f0))
-    return CartesianIndex(i2, i1)
-end
-
-"""
-Point in (x, y) format.
-"""
-@inline function get_pixel_index(f::Film, p::Point2)
-    i1, i2 = u_int32.((p .- f.crop_bounds.p_min .+ 1f0))
-    return CartesianIndex(i2, i1)
-end
-
-function merge_film_tile!(f::Film, ft::FilmTile)
-    x_range = Int(ft.bounds.p_min[1]):Int(ft.bounds.p_max[1])
-    y_range = Int(ft.bounds.p_min[2]):Int(ft.bounds.p_max[2])
-    ft_contrib_sum = ft.pixels.contrib_sum
-    ft_filter_weight_sum = ft.pixels.filter_weight_sum
-    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{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])
-        f_filter_weight_sum[f_idx] += ft_filter_weight_sum[ft_idx]
-    end
-end
-
-
 function set_image!(f::Film, spectrum::Matrix{S}) where {S<:Spectrum}
     @real_assert size(f.pixels) == size(spectrum)
     f.pixels.xyz .= to_XYZ.(spectrum)

src/integrators/sampler.jl

@@ -155,9 +155,9 @@ function li(
 end
 
 @inline function specular_reflect(
-    bsdf, sampler, max_depth, ray::RayDifferentials,
-    si::SurfaceInteraction, scene::Scene, depth::Int32,
-)
+        bsdf, sampler, max_depth, ray::RayDifferentials,
+        si::SurfaceInteraction, scene::Scene, depth::Int32,
+    )
 
     # Compute specular reflection direction `wi` and BSDF value.
     wo = si.core.wo
@@ -408,29 +408,15 @@ end
 end
 
 
-struct Tile
-    tile_indx::NTuple{2, Int32}
-    width::Int32
+@inline function trace_pixel(camera, scene, pixel, sampler, max_depth)
+    camera_sample = get_camera_sample(sampler, pixel)
+    ray, ω = generate_ray_differential(camera, camera_sample)
+    if ω > 0.0f0
+        return li_iterative(sampler, max_depth, ray, scene)
+    end
+    return RGBSpectrum(0.0f0)
 end
 
-# function sample_kernel_inner(i::A, scene::B, t_sampler::C, film::D, film_tile::E, camera::F, pixel::G, spp_sqr::H) where {A,B,C,D,E,F,G,H}
-#     for _ in 1:t_sampler.samples_per_pixel
-#         camera_sample = get_camera_sample(t_sampler, pixel)
-#         ray, ω = generate_ray_differential(camera, camera_sample)
-#         ray = scale_differentials(ray, spp_sqr)
-#         l = RGBSpectrum(0.0f0)
-#         if ω > 0.0f0
-#             # l = li(t_sampler, Int32(i.max_depth), ray, scene, Int32(1))
-#             l = li_iterative(t_sampler, Int32(i.max_depth), ray, scene)
-#         end
-#         # TODO check l for invalid values
-#         if isnan(l)
-#             l = RGBSpectrum(0.0f0)
-#         end
-#         add_sample!(film, film_tile, camera_sample.film, l, ω)
-#     end
-# end
-
 @noinline function sample_tile(sampler, camera, scene, film, film_tile, tile_bounds, max_depth)
     spp_sqr = 1.0f0 / √Float32(sampler.samples_per_pixel)
     for pixel in tile_bounds

src/kernel-abstractions.jl

@@ -63,29 +63,15 @@ function to_gpu(ArrayType, scene::Trace.Scene; preserve=[])
     return Trace.Scene(scene.lights, bvh, scene.bound)
 end
 
-
-@inline function trace_pixel(camera, scene, pixel, samples_per_pixel, max_depth, l)
-    s = UniformSampler(samples_per_pixel)
-    camera_sample = get_camera_sample(s, pixel)
-    ray, ω = generate_ray_differential(camera, camera_sample)
-    if ω > 0.0f0
-        l += li_iterative(s, max_depth, ray, scene)
-    end
-    return l
-end
-
-@kernel function ka_trace_image!(img, camera, scene, samples_per_pixel, max_depth, niter)
+@kernel function ka_trace_image!(img, camera, scene, sampler, max_depth)
     _idx = @index(Global)
     idx = _idx % Int32
     @inbounds if checkbounds(Bool, img, idx)
         cols = size(img, 2) % Int32
         row = (idx - Int32(1)) ÷ cols + Int32(1)
         col = (idx - Int32(1)) % cols + Int32(1)
-        rgb = img[idx]
-        ninv = 0.5f0
-        l = RGBSpectrum(ninv * rgb.r, ninv * rgb.g, ninv * rgb.b)
         pixel = Point2f((row, cols - col))
-        l = trace_pixel(camera, scene, pixel, samples_per_pixel, max_depth, l)
+        l = trace_pixel(camera, scene, pixel, sampler, max_depth)
         img[idx] = RGB{Float32}(( l.c)...)
     end
     nothing
@@ -94,7 +80,8 @@ end
 function launch_trace_image!(img, camera, scene, samples_per_pixel::Int32, max_depth::Int32, niter::Int32)
     backend = KA.get_backend(img)
     kernel! = ka_trace_image!(backend)
-    kernel!(img, camera, scene, samples_per_pixel, max_depth, niter, ndrange=size(img), workgroupsize=(16, 16))
+    sampler = UniformSampler(samples_per_pixel)
+    kernel!(img, camera, scene, sampler, max_depth, niter, ndrange=size(img), workgroupsize=(16, 16))
     KA.synchronize(backend)
     return img
 end