La/pow issue

docs/Project.toml

@@ -1,10 +1,15 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
+ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
 Trace = "afc56b53-c9a9-482a-a956-d1d800e05558"
 
 [compat]

src/Trace.jl

@@ -52,7 +52,7 @@ function concentric_sample_disk(u::Point2f)::Point2f
     # Map uniform random numbers to [-1, 1].
     offset = 2f0 * u - Vec2f(1f0)
     # Handle degeneracy at the origin.
-    offset[1] ≈ 0 && offset[2] ≈ 0 && return Point2f(0)
+    offset[1] ≈ 0f0 && offset[2] ≈ 0f0 && return Point2f(0)
     if abs(offset[1]) > abs(offset[2])
         r = offset[1]
         θ = (offset[2] / offset[1]) * π / 4f0

src/materials/material.jl

@@ -59,7 +59,7 @@ Base.Base.@propagate_inbounds function glass_material(g::UberMaterial, si::Surfa
     t_black = is_black(t)
     r_black && t_black && return BSDF(si, η)
 
-    is_specular = u_roughness ≈ 0 && v_roughness ≈ 0
+    is_specular = u_roughness ≈ 0f0 && v_roughness ≈ 0f0
     if is_specular && allow_multiple_lobes
         return BSDF(si, η, FresnelSpecular(true, r, t, 1.0f0, η, transport))
     end

src/materials/tonemap.jl

@@ -1,17 +1,17 @@
 using Colors, Statistics
 
-luminosity(c::RGB{T}) where {T} = (max(c.r, c.g, c.b) + min(c.r, c.g, c.b)) / 2.0
+luminosity(c::RGB{T}) where {T} = (max(c.r, c.g, c.b) + min(c.r, c.g, c.b)) / 2.0f0
 
 function lum_max(rgb_m)
-    lum_max = 0.0
+    lum_max = 0.0f0
     for pix in rgb_m
         (lum_max > luminosity(pix)) || (lum_max = luminosity(pix))
     end
     lum_max
 end
 
-function avg_lum(rgb_m, δ::Number=1e-10)
-    cumsum = 0.0
+function avg_lum(rgb_m, δ::Number=1f-10)
+    cumsum = 0.0f0
     for pix in rgb_m
         cumsum += log10(δ + luminosity(pix))
     end
@@ -20,13 +20,13 @@ end
 
 function normalize_image(
         rgb_m,
-        a::Float64=0.18,
+        a::Float32=0.18f0,
         lum::Union{Number,Nothing}=nothing,
-        δ::Number=1e-10
+        δ::Number=1f-10
     )
 
-    (isnothing(lum) || lum ≈ 0.0) && (lum = avg_lum(rgb_m, δ))
-    return rgb_m .* a .* (1.0 / lum)
+    (isnothing(lum) || lum ≈ 0.0f0) && (lum = avg_lum(rgb_m, δ))
+    return rgb_m .* a .* (1.0f0 / lum)
 end
 
 function clamp_image(img::AbstractMatrix{T}) where {T}
@@ -35,7 +35,7 @@ function clamp_image(img::AbstractMatrix{T}) where {T}
     end
 end
 
-function γ_correction(img::AbstractMatrix{T}, γ::Float64=1.0, k::Float64=1.0) where T
+function γ_correction(img::AbstractMatrix{T}, γ::Float32=1.0f0, k::Float32=1.0f0) where T
     return map(img) do c
         return T(
             floor(255 * c.r^(1 / γ)),
@@ -46,8 +46,8 @@ function γ_correction(img::AbstractMatrix{T}, γ::Float64=1.0, k::Float64=1.0)
 end
 
 function tone_mapping(img;
-        a::Float64=0.18,
-        γ::Float64=1.0,
+        a::Float32=0.18f0,
+        γ::Float32=1.0f0,
         lum::Union{Number,Nothing}=nothing
     )
     img = normalize_image(img, a, lum)

src/reflection/microfacet.jl

@@ -95,8 +95,11 @@ with the surface normal `w`.
 function D(trd::TrowbridgeReitzDistribution, w::Vec3f)::Float32
     tan_θ² = tan_θ(w)^2
     isinf(tan_θ²) && return 0f0
-
-    cos_θ⁴ = cos_θ(w)^4
+
+    # Calculate cos_θ⁴ without using ^4
+    cos_θ² = cos_θ(w) * cos_θ(w)
+    cos_θ⁴ = cos_θ² * cos_θ²
+
     e = (cos_ϕ(w)^2 / (trd.α_x^2) + sin_ϕ(w)^2 / (trd.α_y^2)) * tan_θ²
     1f0 / (π * trd.α_x * trd.α_y * cos_θ⁴ * (1f0 + e)^2)
 end
@@ -210,7 +213,7 @@ function distribution_microfacet_reflection(m::UberBxDF{S}, wo::Vec3f, wi::Vec3f
     wh = wi + wo
     # Degenerate cases for microfacet reflection.
 
-    (cosθi ≈ 0 || cosθo ≈ 0) && return S(0f0)
+    (cosθi ≈ 0f0 || cosθo ≈ 0f0) && return S(0f0)
     wh ≈ Vec3f(0) && return S(0f0)
     wh = normalize(wh)
     f = m.fresnel(wi ⋅ face_forward(wh, Vec3f(0, 0, 1)))
@@ -222,12 +225,12 @@ end
         m::UberBxDF{S}, wo::Vec3f, u::Point2f,
     )::Tuple{Vec3f,Float32,RGBSpectrum,UInt8} where {S<:Spectrum}
 
-    wo[3] ≈ 0 && return Vec3f(0.0f0), 0.0f0, S(0.0f0), UInt8(0)
+    wo[3] ≈ 0f0 && return Vec3f(0.0f0), 0.0f0, S(0.0f0), UInt8(0)
 
     # Sample microfacet orientation `wh` and reflected direction `wi`.
 
     wh = sample_wh(m.distribution, wo, u)
-    (wo ⋅ wh) < 0 && return Vec3f(0.0f0), 0.0f0, S(0.0f0), UInt8(0)
+    (wo ⋅ wh) < 0f0 && return Vec3f(0.0f0), 0.0f0, S(0.0f0), UInt8(0)
 
     wi = reflect(wo, wh)
     !same_hemisphere(wo, wi) && return Vec3f(0.0f0), 0.0f0, S(0.0f0), UInt8(0)
@@ -252,7 +255,7 @@ function distribution_microfacet_transmission(m::UberBxDF{S}, wo::Vec3f, wi::Vec
     same_hemisphere(wo, wi) && return S(0f0)
 
     cosθo, cosθi = cos_θ(wo), cos_θ(wi)
-    (cosθo ≈ 0 || cosθi ≈ 0) && return S(0f0)
+    (cosθo ≈ 0f0 || cosθi ≈ 0f0) && return S(0f0)
     # Compute `wh` from `wo` & `wi` for microfacet transmission.
     η = cos_θ(wo) > 0f0 ? (m.η_b / m.η_a) : (m.η_a / m.η_b)
     wh = normalize((wo + wi * η))
@@ -275,7 +278,7 @@ end
 
 @inline function sample_microfacet_transmission(m::UberBxDF{S}, wo::Vec3f, u::Point2f) where {S<:Spectrum}
 
-    wo[3] ≈ 0 && return Vec3f(0f0), 0f0, S(0f0), UInt8(0)
+    wo[3] ≈ 0f0 && return Vec3f(0f0), 0f0, S(0f0), UInt8(0)
     wh = sample_wh(m.distribution, wo, u)
     (wo ⋅ wh) < 0 && return Vec3f(0f0), 0f0, S(0f0), UInt8(0)
 

src/shapes/sphere.jl

@@ -55,7 +55,7 @@ end
 
 function refine_intersection(p::Point, s::Sphere)
     p *= s.radius ./ distance(Point3f(0), p)
-    p[1] ≈ 0 && p[2] ≈ 0 && (p = Point3f(1f-6 * s.radius, p[2], p[3]))
+    p[1] ≈ 0f0 && p[2] ≈ 0f0 && (p = Point3f(1f-6 * s.radius, p[2], p[3]))
     p
 end
 

src/shapes/triangle_mesh.jl

@@ -181,7 +181,7 @@ end
     δuv_13, δuv_23 = uv[1] - uv[3], uv[2] - uv[3]
     δp_13, δp_23 = Vec3f(vs[1] - vs[3]), Vec3f(vs[2] - vs[3])
     det = δuv_13[1] * δuv_23[2] - δuv_13[2] * δuv_23[1]
-    if det ≈ 0
+    if det ≈ 0f0
         v = normalize((vs[3] - vs[1]) × (vs[2] - vs[1]))
         _, ∂p∂u, ∂p∂v = coordinate_system(Vec3f(v))
         return ∂p∂u, ∂p∂v, δp_13, δp_23
@@ -201,7 +201,7 @@ end
     δuv_13, δuv_23 = uv[1] - uv[3], uv[2] - uv[3]
     δn_13, δn_23 = t_normals[1] - t_normals[3], t_normals[2] - t_normals[3]
     det = δuv_13[1] * δuv_23[2] - δuv_13[2] * δuv_23[1]
-    det ≈ 0 && return Normal3f(0), Normal3f(0)
+    det ≈ 0f0 && return Normal3f(0), Normal3f(0)
 
     inv_det = 1f0 / det
     ∂n∂u = (δuv_23[2] * δn_13 - δuv_13[2] * δn_23) * inv_det

test/Project.toml

@@ -1,8 +1,15 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
+ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Trace = "afc56b53-c9a9-482a-a956-d1d800e05558"

test/gpu-threading-benchmarks.jl

@@ -100,6 +100,13 @@ function launch_trace_image!(img, camera, scene)
     KA.synchronize(backend)
     return img
 end
+# using AMDGPU
+# ArrayType = ROCArray
+# using CUDA
+# ArrayType = CuArray
+
+# using Metal
+# ArrayType = MtlArray
 
 preserve = []
 gpu_scene = to_gpu(ArrayType, scene; preserve=preserve);