diff --git a/docs/examples.jl b/docs/examples.jl
index b05d660..054d487 100644
--- a/docs/examples.jl
+++ b/docs/examples.jl
@@ -21,6 +21,9 @@ begin
     viewdir = normalize(ax.scene.camera.view_direction[])
 end
 
+hitpoints, centroid = RayCaster.get_centroid(bvh, viewdir)
+
+
 begin
     @time "hitpoints" hitpoints, centroid = RayCaster.get_centroid(bvh, viewdir)
     @time "illum" illum = RayCaster.get_illumination(bvh, viewdir)
@@ -33,10 +36,11 @@ begin
     f, ax, pl = mesh(world_mesh, color=:blue)
     per_face_vf = FaceView((viewfacts), [GLTriangleFace(i) for i in 1:N])
     viewfact_mesh = GeometryBasics.mesh(world_mesh, color=per_face_vf)
-    pl = Makie.mesh(f[1, 2],
+    pl = Makie.mesh(
+        f[1, 2],
         viewfact_mesh, colormap=[:black, :red], axis=(; show_axis=false),
-        shading=false, highclip=:red, lowclip=:black, colorscale=sqrt,)
-
+        shading=false, highclip=:red, lowclip=:black,  colorscale=sqrt,
+    )
     # Centroid
     cax, pl = Makie.mesh(f[2, 1], world_mesh, color=(:blue, 0.5), axis=(; show_axis=false), transparency=true)
 
@@ -58,3 +62,113 @@ begin
 
     f
 end
+
+
+using KernelAbstractions, Atomix
+
+function random_scatter_kernel!(bvh, triangle, u, v, normal)
+    point = RayCaster.random_triangle_point(triangle)
+    o = point .+ (normal .* 0.01f0) # Offset so it doesn't self intersect
+    dir = RayCaster.random_hemisphere_uniform(normal, u, v)
+    ray = RayCaster.Ray(; o=o, d=dir)
+    hit, prim, _ = RayCaster.intersect!(bvh, ray)
+    return hit, prim
+end
+
+import GeometryBasics as GB
+
+@kernel function viewfact_ka_kernel!(result, bvh, primitives, rays_per_triangle)
+    idx = @index(Global)
+    prim_idx = ((UInt32(idx) - UInt32(1)) ÷ rays_per_triangle) + UInt32(1)
+    if prim_idx <= length(primitives)
+        triangle, u, v, normal = primitives[prim_idx]
+        hit, prim = random_scatter_kernel!(bvh, triangle, u, v, normal)
+        if hit && prim.material_idx !== triangle.material_idx
+            # weigh by angle?
+            Atomix.@atomic result[triangle.material_idx, prim.material_idx] += 1
+        end
+    end
+end
+
+function view_factors!(result, bvh, prim_info, rays_per_triangle=10000)
+
+    backend = get_backend(result)
+    workgroup = 256
+    total_rays = length(bvh.primitives) * rays_per_triangle
+    per_workgroup = total_rays ÷ workgroup
+    final_rays = per_workgroup * workgroup
+    per_triangle = final_rays ÷ length(bvh.primitives)
+
+    kernel = viewfact_ka_kernel!(backend, 256)
+    kernel(result, bvh, prim_info, UInt32(per_triangle); ndrange = final_rays)
+    return result
+end
+
+result = zeros(UInt32, length(bvh.primitives), length(bvh.primitives))
+using AMDGPU
+prim_info = map(bvh.primitives) do triangle
+    n = GB.orthogonal_vector(Vec3f, GB.Triangle(triangle.vertices...))
+    normal = normalize(Vec3f(n))
+    u, v = RayCaster.get_orthogonal_basis(normal)
+    return triangle, u, v, normal
+end
+bvh_gpu = RayCaster.to_gpu(ROCArray, bvh)
+result_gpu = ROCArray(result)
+prim_info_gpu = ROCArray(prim_info)
+@time begin
+    view_factors!(result_gpu, bvh_gpu, prim_info_gpu, 10000)
+    KernelAbstractions.synchronize(get_backend(result_gpu))
+end;
+
+
+
+@kernel function viewfact_ka_kernel2!(result, bvh, primitives, rays_per_triangle)
+    idx = @index(Global)
+    prim_idx = ((UInt32(idx) - UInt32(1)) ÷ rays_per_triangle) + UInt32(1)
+    if prim_idx <= length(primitives)
+        triangle, u, v, normal = primitives[prim_idx]
+        hit, prim = random_scatter_kernel!(bvh, triangle, u, v, normal)
+        if hit && prim.material_idx !== triangle.material_idx
+            # weigh by angle?
+            @inbounds result[idx] = UInt32(1)
+        end
+    end
+end
+
+
+function view_factors2!(result, bvh, prim_info, per_triangle)
+    backend = get_backend(result)
+    kernel = viewfact_ka_kernel2!(backend, 256)
+    kernel(result, bvh, prim_info, UInt32(per_triangle); ndrange = length(result))
+    return result
+end
+
+
+using AMDGPU
+workgroup = 256
+rays_per_triangle = 10000
+total_rays = length(bvh.primitives) * rays_per_triangle
+per_workgroup = total_rays ÷ workgroup
+final_rays = per_workgroup * workgroup
+per_triangle = final_rays ÷ length(bvh.primitives)
+result = zeros(UInt32, final_rays)
+
+final_rays / 10^6
+
+prim_info = map(bvh.primitives) do triangle
+    n = GB.orthogonal_vector(Vec3f, GB.Triangle(triangle.vertices...))
+    normal = normalize(Vec3f(n))
+    u, v = RayCaster.get_orthogonal_basis(normal)
+    return triangle, u, v, normal
+end
+
+bvh_gpu = RayCaster.to_gpu(ROCArray, bvh)
+result_gpu = ROCArray(result)
+prim_info_gpu = ROCArray(prim_info)
+@time begin
+    view_factors2!(result_gpu, bvh_gpu, prim_info_gpu, per_triangle)
+    KernelAbstractions.synchronize(get_backend(result_gpu))
+end;
+
+@time view_factors2!(result, bvh, prim_info, per_triangle)
+@code_warntype random_scatter_kernel!(bvh, prim_info[1]...)
diff --git a/src/bvh.jl b/src/bvh.jl
index f5a9aeb..9046d1b 100644
--- a/src/bvh.jl
+++ b/src/bvh.jl
@@ -239,100 +239,151 @@ end
     length(bvh.nodes) > Int32(0) ? bvh.nodes[1].bounds : Bounds3()
 end
 
-@inline function intersect!(bvh::BVHAccel{P}, ray::AbstractRay) where {P}
-    hit = false
-    interaction = SurfaceInteraction()
+"""
+    _traverse_bvh(bvh::BVHAccel{P}, ray::AbstractRay, hit_callback::F) where {P, F<:Function}
+
+Internal function that traverses the BVH to find ray-primitive intersections.
+Uses a callback pattern to handle different intersection behaviors.
+
+Arguments:
+- `bvh`: The BVH acceleration structure
+- `ray`: The ray to test for intersections
+- `hit_callback`: Function called when primitive is tested. Signature:
+   hit_callback(primitive, ray) -> (continue_traversal::Bool, ray::AbstractRay, results::Any)
+
+Returns:
+- The final result from the hit_callback
+"""
+@inline function traverse_bvh(hit_callback::F, bvh::BVHAccel{P}, ray::AbstractRay) where {P, F<:Function}
+    # Early return if BVH is empty
+    if length(bvh.nodes) == 0
+        return false, ray, nothing
+    end
+
+    # Prepare ray for traversal
     ray = check_direction(ray)
     inv_dir = 1f0 ./ ray.d
     dir_is_neg = is_dir_negative(ray.d)
 
-    to_visit_offset, current_node_i = Int32(1), Int32(1)
+    # Initialize traversal stack
+    to_visit_offset = Int32(1)
+    current_node_idx = Int32(1)
     nodes_to_visit = zeros(MVector{64,Int32})
     primitives = bvh.primitives
-    @_inbounds primitive = primitives[1]
     nodes = bvh.nodes
+
+    # State variables to hold callback results
+    continue_search = true
+    prim1 = primitives[1]
+    result = hit_callback(prim1, ray, nothing)
+
+    # Traverse BVH
     @_inbounds while true
-        ln = nodes[current_node_i]
-        if intersect_p(ln.bounds, ray, inv_dir, dir_is_neg)
-            if !ln.is_interior && ln.n_primitives > Int32(0)
-                # Intersect ray with primitives in node.
-                for i in Int32(0):ln.n_primitives - Int32(1)
-                    offset = ln.offset % Int32
-                    tmp_primitive = primitives[offset+i]
-                    tmp_hit, ray, tmp_interaction = intersect_p!(
-                        tmp_primitive, ray,
-                    )
-                    if tmp_hit
-                        hit = tmp_hit
-                        interaction = tmp_interaction
-                        primitive = tmp_primitive
+        current_node = nodes[current_node_idx]
+
+        # Test ray against current node's bounding box
+        if intersect_p(current_node.bounds, ray, inv_dir, dir_is_neg)
+            if !current_node.is_interior && current_node.n_primitives > Int32(0)
+                # Leaf node - test all primitives
+                offset = current_node.offset % Int32
+
+                for i in Int32(0):(current_node.n_primitives - Int32(1))
+                    primitive = primitives[offset + i]
+
+                    # Call the callback for this primitive
+                    continue_search, ray, result = hit_callback(primitive, ray, result)
+
+                    # Early exit if callback requests it
+                    if !continue_search
+                        return false, ray, result
                     end
                 end
-                to_visit_offset == Int32(1) && break
+
+                # Done with leaf, pop next node from stack
+                if to_visit_offset == Int32(1)
+                    break
+                end
                 to_visit_offset -= Int32(1)
-                current_node_i = nodes_to_visit[to_visit_offset]
+                current_node_idx = nodes_to_visit[to_visit_offset]
             else
-                if dir_is_neg[ln.split_axis] == Int32(2)
-                    nodes_to_visit[to_visit_offset] = current_node_i + Int32(1)
-                    current_node_i = ln.offset % Int32
+                # Interior node - push children to stack
+                if dir_is_neg[current_node.split_axis] == Int32(2)
+                    nodes_to_visit[to_visit_offset] = current_node_idx + Int32(1)
+                    current_node_idx = current_node.offset % Int32
                 else
-                    nodes_to_visit[to_visit_offset] = ln.offset % Int32
-                    current_node_i += Int32(1)
+                    nodes_to_visit[to_visit_offset] = current_node.offset % Int32
+                    current_node_idx += Int32(1)
                 end
                 to_visit_offset += Int32(1)
             end
         else
-            to_visit_offset == 1 && break
+            # Miss - pop next node from stack
+            if to_visit_offset == Int32(1)
+                break
+            end
             to_visit_offset -= Int32(1)
-            current_node_i = nodes_to_visit[to_visit_offset]
+            current_node_idx = nodes_to_visit[to_visit_offset]
         end
     end
-    return hit, primitive, interaction
+
+    # Return final state
+    return continue_search, ray, result
 end
 
-@inline function intersect_p(bvh::BVHAccel, ray::AbstractRay)
+# Initialization
+closest_hit_callback(primitive, ray, ::Nothing) = (false, primitive, Point3f(0.0))
 
-    length(bvh.nodes) == Int32(0) && return false
+function closest_hit_callback(primitive, ray, prev_result::Tuple{Bool, P, Point3f}) where {P}
+    # Test intersection and update if closer
+    tmp_hit, ray, tmp_bary = intersect_p!(primitive, ray)
+    # Always continue search to find closest
+    return true, ray, ifelse(tmp_hit,  (true, primitive, tmp_bary), prev_result)
+end
 
-    ray = check_direction(ray)
-    inv_dir = 1f0 ./ ray.d
-    dir_is_neg = is_dir_negative(ray.d)
+"""
+    intersect!(bvh::BVHAccel{P}, ray::AbstractRay) where {P}
 
-    to_visit_offset, current_node_i = Int32(1), Int32(1)
-    nodes_to_visit = zeros(MVector{64,Int32})
-    primitives = bvh.primitives
-    @_inbounds while true
-        ln = bvh.nodes[current_node_i]
-        if intersect_p(ln.bounds, ray, inv_dir, dir_is_neg)
-            if !ln.is_interior && ln.n_primitives > Int32(0)
-                for i in Int32(0):ln.n_primitives-Int32(1)
-                    offset = ln.offset % Int32
-                    intersect_p(
-                        primitives[offset + i], ray,
-                    ) && return true
-                end
-                to_visit_offset == 1 && break
-                to_visit_offset -= Int32(1)
-                current_node_i = nodes_to_visit[to_visit_offset]
-            else
-                if dir_is_neg[ln.split_axis] == Int32(2)
-                    # @setindex 64 nodes_to_visit[to_visit_offset] = Int32(current_node_i + 1)
-                    nodes_to_visit[to_visit_offset] = current_node_i + Int32(1)
-                    current_node_i = ln.offset % Int32
-                else
-                    # @setindex 64 nodes_to_visit[to_visit_offset] = Int32(ln.offset)
-                    nodes_to_visit[to_visit_offset] = ln.offset % Int32
-                    current_node_i += Int32(1)
-                end
-                to_visit_offset += Int32(1)
-            end
-        else
-            to_visit_offset == Int32(1) && break
-            to_visit_offset -= Int32(1)
-            current_node_i = Int32(nodes_to_visit[to_visit_offset])
-        end
+Find the closest intersection between a ray and the primitives stored in a BVH.
+
+Returns:
+- `hit_found`: Boolean indicating if an intersection was found
+- `hit_primitive`: The primitive that was hit (if any)
+- `barycentric_coords`: Barycentric coordinates of the hit point
+"""
+@inline function intersect!(bvh::BVHAccel{P}, ray::AbstractRay) where {P}
+    # Traverse BVH with closest-hit callback
+    _, _, result = traverse_bvh(closest_hit_callback, bvh, ray)
+    return result::Tuple{Bool, Triangle, Point3f}
+end
+
+
+any_hit_callback(primitive, current_ray, result::Nothing) = ()
+
+# Define any-hit callback
+function any_hit_callback(primitive, current_ray, ::Tuple{})
+    # Test for intersection
+    if intersect_p(primitive, current_ray)
+        # Stop traversal on first hit
+        return false, current_ray, true
     end
-    false
+    # Continue search if no hit
+    return true, current_ray, false
+end
+
+"""
+    intersect_p(bvh::BVHAccel, ray::AbstractRay)
+
+Test if a ray intersects any primitive in the BVH (without finding the closest hit).
+
+Returns:
+- `hit_found`: Boolean indicating if any intersection was found
+"""
+@inline function intersect_p(bvh::BVHAccel, ray::AbstractRay)
+    # Traverse BVH with any-hit callback
+    continue_search, _, result = traverse_bvh(any_hit_callback, bvh, ray)
+    # If traversal completed without finding a hit, return false
+    # Otherwise return the hit result (true)
+    return !continue_search ? result : false
 end
 
 function calculate_ray_grid_bounds(bounds::GeometryBasics.Rect, ray_direction::Vec3f)
diff --git a/src/kernel-abstractions.jl b/src/kernel-abstractions.jl
index 6943c9e..f346a1c 100644
--- a/src/kernel-abstractions.jl
+++ b/src/kernel-abstractions.jl
@@ -19,6 +19,5 @@ end
 function to_gpu(ArrayType, bvh::RayCaster.BVHAccel; preserve=[])
     primitives = to_gpu(ArrayType, bvh.primitives; preserve=preserve)
     nodes = to_gpu(ArrayType, bvh.nodes; preserve=preserve)
-    materials = to_gpu(ArrayType, to_gpu.((ArrayType,), bvh.materials; preserve=preserve); preserve=preserve)
-    return RayCaster.BVHAccel(primitives, materials, bvh.max_node_primitives, nodes)
+    return RayCaster.BVHAccel(primitives, bvh.max_node_primitives, nodes)
 end
diff --git a/src/kernels.jl b/src/kernels.jl
index 734be41..2e7faed 100644
--- a/src/kernels.jl
+++ b/src/kernels.jl
@@ -12,8 +12,9 @@ function hits_from_grid(bvh, viewdir; grid_size=32)
     Threads.@threads for idx in CartesianIndices(ray_origins)
         o = ray_origins[idx]
         ray = RayCaster.Ray(; o=o, d=ray_direction)
-        hit, prim, si = RayCaster.intersect!(bvh, ray)
-        @inbounds result[idx] = RayHit(hit, si.core.p, prim.material_idx)
+        hit, prim, bary = RayCaster.intersect!(bvh, ray)
+        hitpoint = sum_mul(bary, prim.vertices)
+        @inbounds result[idx] = RayHit(hit, hitpoint, prim.material_idx)
     end
     return result
 end
@@ -34,7 +35,7 @@ function view_factors!(result, bvh, rays_per_triangle=10000)
                 point_on_triangle = random_triangle_point(triangle)
                 o = point_on_triangle .+ (normal .* 0.01f0) # Offset so it doesn't self intersect
                 ray = Ray(; o=o, d=random_hemisphere_uniform(normal, u, v))
-                hit, prim, si = intersect!(bvh, ray)
+                hit, prim, _ = intersect!(bvh, ray)
                 if hit && prim.material_idx != triangle.material_idx
                     # weigh by angle?
                     result[triangle.material_idx, prim.material_idx] += 1
diff --git a/src/ray.jl b/src/ray.jl
index ada7826..93739bf 100644
--- a/src/ray.jl
+++ b/src/ray.jl
@@ -61,11 +61,9 @@ end
 increase_hit(ray::Ray, t_hit) = Ray(ray; t_max=t_hit)
 increase_hit(ray::RayDifferentials, t_hit) = RayDifferentials(ray; t_max=t_hit)
 
-@inline function intersect_p!(
-        shape::AbstractShape, ray::R,
-    )::Tuple{Bool, R, SurfaceInteraction} where {R<:AbstractRay}
-    intersects, t_hit, interaction = intersect(shape, ray)
-    !intersects && return false, ray, interaction
+@inline function intersect_p!(shape::AbstractShape, ray::R) where {R<:AbstractRay}
+    intersects, t_hit, barycentric = intersect(shape, ray)
+    !intersects && return false, ray, barycentric
     ray = increase_hit(ray, t_hit)
-    return true, ray, interaction
+    return true, ray, barycentric
 end
diff --git a/src/shapes/sphere.jl b/src/shapes/sphere.jl
index 7bc4bbd..b9b4371 100644
--- a/src/shapes/sphere.jl
+++ b/src/shapes/sphere.jl
@@ -85,7 +85,7 @@ end
 """
 Compute partial derivatives of intersection point in parametric form.
 """
-function ∂p(s::Sphere, p::Point3f, θ::Float32, sin_ϕ::Float32, cos_ϕ::Float32)
+function partial_derivatives(s::Sphere, p::Point3f, θ::Float32, sin_ϕ::Float32, cos_ϕ::Float32)
     ∂p∂u = Vec3f(-s.ϕ_max * p[2], s.ϕ_max * p[1], 0f0)
     ∂p∂v = (s.θ_max - s.θ_min) * Vec3f(
         p[3] * cos_ϕ, p[3] * sin_ϕ, -s.radius * sin(θ),
@@ -93,7 +93,7 @@ function ∂p(s::Sphere, p::Point3f, θ::Float32, sin_ϕ::Float32, cos_ϕ::Float
     ∂p∂u, ∂p∂v, sin_ϕ, cos_ϕ
 end
 
-function ∂n(
+function normal_derivatives(
     s::Sphere, p::Point3f,
     sin_ϕ::Float32, cos_ϕ::Float32,
     ∂p∂u::Vec3f, ∂p∂v::Vec3f,
@@ -154,8 +154,8 @@ function intersect(
     v = (θ - s.θ_min) / (s.θ_max - s.θ_min)
 
     sin_ϕ, cos_ϕ = precompute_ϕ(hit_point)
-    ∂p∂u, ∂p∂v = ∂p(s, hit_point, θ, sin_ϕ, cos_ϕ)
-    ∂n∂u, ∂n∂v = ∂n(s, hit_point, sin_ϕ, cos_ϕ, ∂p∂u, ∂p∂v)
+    ∂p∂u, ∂p∂v = partial_derivatives(s, hit_point, θ, sin_ϕ, cos_ϕ)
+    ∂n∂u, ∂n∂v = normal_derivatives(s, hit_point, sin_ϕ, cos_ϕ, ∂p∂u, ∂p∂v)
     reverse_normal = (s.core.reverse_orientation ⊻ s.core.transform_swaps_handedness)
     si = SurfaceInteraction(
         hit_point, ray.time, -ray.d, Point2f(u, v),
diff --git a/src/shapes/triangle_mesh.jl b/src/shapes/triangle_mesh.jl
index 4da875f..21cb904 100644
--- a/src/shapes/triangle_mesh.jl
+++ b/src/shapes/triangle_mesh.jl
@@ -178,7 +178,7 @@ end
     return SVector{3, Point3f}(tvs), shear
 end
 
-@inline function ∂p(
+@inline function partial_derivatives(
         ::Triangle, vs::AbstractVector{Point3f}, uv::AbstractVector{Point2f},
     )::Tuple{Vec3f,Vec3f,Vec3f,Vec3f}
 
@@ -202,7 +202,7 @@ end
 end
 
 
-@inline function ∂n(
+@inline function normal_derivatives(
         t::Triangle, uv::AbstractVector{Point2f},
     )::Tuple{Normal3f,Normal3f}
     t_normals = normals(t)
@@ -220,97 +220,88 @@ end
 end
 
 
-@inline function _init_triangle_shading_geometry(
-        t::Triangle, si::SurfaceInteraction,
-        barycentric::Point3f, uv::AbstractVector{Point2f},
+@inline function init_triangle_shading_geometry(
+        triangle::Triangle, surf_interact::SurfaceInteraction,
+        bary_coords::Point3f, tex_coords::AbstractVector{Point2f},
     )
-    has_normals = _all(x -> _all(isfinite, x), t.normals)
-    has_tangents = _all(x -> _all(isfinite, x), t.tangents)
-    !has_normals && !has_tangents && return si
-    # Initialize triangle shading geometry.
-    # Compute shading normal, tangent & bitangent.
-    ns = si.core.n
+    # Check if the triangle has valid normal and tangent vectors
+    has_normals = _all(x -> _all(isfinite, x), triangle.normals)
+    has_tangents = _all(x -> _all(isfinite, x), triangle.tangents)
+
+    # If no valid shading geometry exists, return the original surface interaction
+    !has_normals && !has_tangents && return surf_interact
+
+    # Initialize triangle shading geometry by computing shading normal, tangent & bitangent
+    shading_normal = surf_interact.core.n  # Start with geometric normal
+
+    # If we have valid normals, interpolate them using barycentric coordinates
     if has_normals
-        ns = normalize(sum_mul(barycentric, t.normals))
+        shading_normal = normalize(sum_mul(bary_coords, triangle.normals))
     end
+
+    # Calculate shading tangent - either from triangle tangents or from position derivatives
+    shading_tangent = Vector3f(0)
     if has_tangents
-        ss = normalize(sum_mul(barycentric, t.tangents))
-    else
-        ss = normalize(si.∂p∂u)
-    end
-    ts = ns × ss
-    if (ts ⋅ ts) > 0f0
-        ts = Vec3f(normalize(ts))
-        ss = Vec3f(ts × ns)
+        shading_tangent = normalize(sum_mul(bary_coords, triangle.tangents))
     else
-        _, ss, ts = coordinate_system(Vec3f(ns))
+        shading_tangent = normalize(surf_interact.pos_deriv_u)  # Assuming ∂p∂u was renamed to pos_deriv_u
     end
-    ∂n∂u, ∂n∂v = ∂n(t, uv)
-    return set_shading_geometry(si, ss, ts, ∂n∂u, ∂n∂v, true)
-end
 
+    # Calculate shading bitangent from normal and tangent
+    shading_bitangent = shading_normal × shading_tangent
 
-@inline function intersect(
-        t::Triangle, ray::AbstractRay, ::Bool = false,
-    )::Tuple{Bool,Float32,SurfaceInteraction}
+    # Check if bitangent is valid, otherwise create a new coordinate system
+    if (shading_bitangent ⋅ shading_bitangent) > 0f0
+        shading_bitangent = Vec3f(normalize(shading_bitangent))
+        shading_tangent = Vec3f(shading_bitangent × shading_normal)  # Ensure orthogonality
+    else
+        # Create a new coordinate system if the vectors are nearly parallel
+        _, shading_tangent, shading_bitangent = coordinate_system(Vec3f(shading_normal))
+    end
 
-    vs = vertices(t)
-    hit, t_hit, barycentric = intersect_triangle(vs, ray)
-    !hit && return false, t_hit, SurfaceInteraction()
-    # TODO check that t_hit > 0
-    uv = uvs(t)
-    ∂p∂u, ∂p∂v, δp_13, δp_23 = ∂p(t, vs, uv)
-    # Interpolate (u, v) paramteric coordinates and hit point.
-    hit_point = sum_mul(barycentric, vs)
-    uv_hit = sum_mul(barycentric, uv)
-    normal = normalize(δp_13 × δp_23)
-
-    si = SurfaceInteraction(
-        normal, hit_point, ray.time, -ray.d, uv_hit,
-        ∂p∂u, ∂p∂v, Normal3f(0f0), Normal3f(0f0)
+    # Calculate normal derivatives
+    nd_u, nd_v = normal_derivatives(triangle, tex_coords)
+
+    # Set the shading geometry on the surface interaction
+    return set_shading_geometry(
+        surf_interact,
+        shading_tangent,
+        shading_bitangent,
+        nd_u,
+        nd_v,
+        true
     )
-    si = _init_triangle_shading_geometry(t, si, barycentric, uv)
-    # TODO test against alpha texture if present.
-    return true, t_hit, si
 end
 
-function intersect_triangle(
-    ray_origin::Point3f,
-    ray_direction::Vec3f,
-    v0::Point3f,
-    v1::Point3f,
-    v2::Point3f)
 
-    edge1 = v1 - v0
-    edge2 = v2 - v0
+function surface_interaction(triangle, ray, bary_coords)
 
-    h = cross(ray_direction, edge2)
-    determinant = dot(edge1, h)
+    verts = vertices(triangle)
+    tex_coords = uvs(triangle)  # Get texture coordinates
 
-    if abs(determinant) < 1e-6
-        return (false, 0.0f0, 0f0)
-    end
-
-    inv_determinant = 1f0 / determinant
-    s = ray_origin - v0
-    u = inv_determinant * dot(s, h)
-
-    if u < 0.0f0  || u > 1.0f0
-        return (false, 0.0f0, 0.0f0)
-    end
-
-    edge1_cross_s = cross(s, edge1)
-    barycentric_v = inv_determinant * dot(ray_direction, edge1_cross_s)
+    # Calculate position derivatives and triangle edges
+    pos_deriv_u, pos_deriv_v, edge1, edge2 = partial_derivatives(triangle, verts, tex_coords)
 
-    if barycentric_v < 0.0f0 || u + barycentric_v > 1.0f0
-        return (false, 0.0f0, 0.0f0)
-    end
+    # Interpolate hit point and texture coordinates using barycentric coordinates
+    hit_point = sum_mul(bary_coords, verts)
+    hit_uv = sum_mul(bary_coords, tex_coords)
 
-    intersection_distance = inv_determinant * dot(edge2, edge1_cross_s)
+    # Calculate surface normal from triangle edges
+    normal = normalize(edge1 × edge2)
 
-    return (true, intersection_distance, inv_determinant)
+    # Create surface interaction data at hit point
+    surf_interact = SurfaceInteraction(
+        normal, hit_point, ray.time, -ray.d, hit_uv,
+        pos_deriv_u, pos_deriv_v, Normal3f(0f0), Normal3f(0f0)
+    )
+    # TODO test against alpha texture if present.
+    return init_triangle_shading_geometry(triangle, surf_interact, bary_coords, tex_coords)
 end
 
+@inline function intersect(triangle::Triangle, ray::AbstractRay)::Tuple{Bool,Float32,Point3f}
+    verts = vertices(triangle)  # Get triangle vertices
+    return intersect_triangle(verts, ray)  # Check if ray hits triangle
+end
 
 @inline function intersect_p(t::Triangle, ray::Union{Ray,RayDifferentials}, ::Bool=false)
     intersect_triangle(t.vertices, ray)[1]
@@ -331,6 +322,7 @@ end
     # Perform triangle edge & determinant tests.
     # Point is inside a triangle if all edges have the same sign.
     any(edges .< 0.0f0) && any(edges .> 0.0f0) && return false, t_hit, barycentric
+
     det = sum(edges)
     det ≈ 0f0 && return false, t_hit, barycentric
     # Compute scaled hit distance to triangle.
diff --git a/test/runtests.jl b/test/runtests.jl
index 3647059..fef5a3c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -38,11 +38,3 @@ end
     @test sb[1] == Point3f(-1f0)
     @test sb[2] == Point3f(1f0)
 end
-
-# LanczosSincFilter functionality removed in refactored RayCaster
-
-# Film functionality removed in refactored RayCaster
-
-# FilmTile functionality removed in refactored RayCaster
-
-# Camera functionality removed in refactored RayCaster
diff --git a/test/test_intersection.jl b/test/test_intersection.jl
index d4cb914..1d6a5a3 100644
--- a/test/test_intersection.jl
+++ b/test/test_intersection.jl
@@ -104,28 +104,23 @@ end
     # In the refactored API, object_bound returns world bounds since transformation is applied during creation
     @test RayCaster.object_bound(triangle) ≈ target_wb
 
-    # Test ray intersection.
+    # Test ray intersection - API has changed: intersect now returns (Bool, Float32, Point3f) with barycentric coords
     ray = RayCaster.Ray(o = Point3f(0, 0, -2), d = Vec3f(0, 0, 1))
     intersects_p = RayCaster.intersect_p(triangle, ray)
-    intersects, t_hit, interaction = RayCaster.intersect(triangle, ray)
+    intersects, t_hit, bary_coords = RayCaster.intersect(triangle, ray)
     @test intersects_p == intersects == true
     @test t_hit ≈ 4f0
-    @test RayCaster.apply(ray, t_hit) ≈ interaction.core.p ≈ Point3f(0, 0, 2)
-    @test interaction.uv ≈ Point2f(0)
-    # Normal is flipped in the refactored API
-    @test interaction.core.n ≈ RayCaster.Normal3f(0, 0, 1)
-    @test interaction.core.wo ≈ -ray.d
-    # Test ray intersection (lower-left corner).
-    ray = RayCaster.Ray(o = Point3f(1, 0.5, 0), d = Vec3f(0, 0, 1))
+    @test RayCaster.apply(ray, t_hit) ≈ Point3f(0, 0, 2)
+    # Barycentric coordinates for vertex 0 (corner hit)
+    @test bary_coords ≈ Point3f(1, 0, 0)
+    
+    # Test ray intersection (different point).
+    ray = RayCaster.Ray(o = Point3f(0.5, 0.25, 0), d = Vec3f(0, 0, 1))
     intersects_p = RayCaster.intersect_p(triangle, ray)
-    intersects, t_hit, interaction = RayCaster.intersect(triangle, ray)
+    intersects, t_hit, bary_coords = RayCaster.intersect(triangle, ray)
     @test intersects_p == intersects == true
     @test t_hit ≈ 2f0
-    @test RayCaster.apply(ray, t_hit) ≈ interaction.core.p ≈ Point3f(1, 0.5, 2)
-    @test interaction.uv ≈ Point2f(1, 0.5)
-    # Normal is flipped in the refactored API
-    @test interaction.core.n ≈ RayCaster.Normal3f(0, 0, 1)
-    @test interaction.core.wo ≈ -ray.d
+    @test RayCaster.apply(ray, t_hit) ≈ Point3f(0.5, 0.25, 2)
 end
 
 # BVH tests with spheres removed - refactored RayCaster only supports triangle meshes in BVH