diff --git a/NEWS.md b/NEWS.md
index 24d6b6024f8..3053736b4c9 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,7 @@
 
 ## master
 
+- Fixed DataInspector interaction with transformations [#3002](https://github.com/MakieOrg/Makie.jl/pull/3002)
 - Fix incomplete stroke with some Bezier markers in CairoMakie and blurry strokes in GLMakie [#2961](https://github.com/MakieOrg/Makie.jl/pull/2961)
 - Added the ability to use custom triangulations from DelaunayTriangulation.jl [#2896](https://github.com/MakieOrg/Makie.jl/pull/2896).
 - Adjusted scaling of scatter/text stroke, glow and anti-aliasing width under non-uniform 2D scaling (Vec2f markersize/fontsize) in GLMakie [#2950](https://github.com/MakieOrg/Makie.jl/pull/2950).
diff --git a/src/Makie.jl b/src/Makie.jl
index 64a19ed55e8..cb3e0f56d59 100644
--- a/src/Makie.jl
+++ b/src/Makie.jl
@@ -176,6 +176,7 @@ include("stats/hexbin.jl")
 # Interactiveness
 include("interaction/events.jl")
 include("interaction/interactive_api.jl")
+include("interaction/ray_casting.jl")
 include("interaction/inspector.jl")
 
 # documentation and help functions
diff --git a/src/interaction/inspector.jl b/src/interaction/inspector.jl
index 414e182f7fd..849c43f8392 100644
--- a/src/interaction/inspector.jl
+++ b/src/interaction/inspector.jl
@@ -73,53 +73,6 @@ function closest_point_on_line(A::Point2f, B::Point2f, P::Point2f)
     A .+ AB * dot(AP, AB) / dot(AB, AB)
 end
 
-function view_ray(scene)
-    inv_projview = inv(camera(scene).projectionview[])
-    view_ray(inv_projview, events(scene).mouseposition[], pixelarea(scene)[])
-end
-function view_ray(inv_view_proj, mpos, area::Rect2)
-    # This figures out the camera view direction from the projectionview matrix (?)
-    # and computes a ray from a near and a far point.
-    # Based on ComputeCameraRay from ImGuizmo
-    mp = 2f0 .* (mpos .- minimum(area)) ./ widths(area) .- 1f0
-    v = inv_view_proj * Vec4f(0, 0, -10, 1)
-    reversed = v[3] < v[4]
-    near = reversed ? 1f0 - 1e-6 : 0f0
-    far = reversed ? 0f0 : 1f0 - 1e-6
-
-    origin = inv_view_proj * Vec4f(mp[1], mp[2], near, 1f0)
-    origin = origin[Vec(1, 2, 3)] ./ origin[4]
-
-    p = inv_view_proj * Vec4f(mp[1], mp[2], far, 1f0)
-    p = p[Vec(1, 2, 3)] ./ p[4]
-
-    dir = normalize(p .- origin)
-    return origin, dir
-end
-
-
-# These work in 2D and 3D
-function closest_point_on_line(A, B, origin, dir)
-    closest_point_on_line(
-        to_ndim(Point3f, A, 0),
-        to_ndim(Point3f, B, 0),
-        to_ndim(Point3f, origin, 0),
-        to_ndim(Vec3f, dir, 0)
-    )
-end
-function closest_point_on_line(A::Point3f, B::Point3f, origin::Point3f, dir::Vec3f)
-    # See:
-    # https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection
-    AB_norm = norm(B .- A)
-    u_AB = (B .- A) / AB_norm
-    u_dir = normalize(dir)
-    u_perp = normalize(cross(u_dir, u_AB))
-    # e_RD, e_perp defines a plane with normal n
-    n = normalize(cross(u_dir, u_perp))
-    t = dot(origin .- A, n) / dot(u_AB, n)
-    A .+ clamp(t, 0.0, AB_norm) * u_AB
-end
-
 function point_in_triangle(A::Point2, B::Point2, C::Point2, P::Point2, ϵ = 1e-6)
     # adjusted from ray_triangle_intersection
     AO = A .- P
@@ -132,39 +85,6 @@ function point_in_triangle(A::Point2, B::Point2, C::Point2, P::Point2, ϵ = 1e-6
     return (A1 > -ϵ && A2 > -ϵ && A3 > -ϵ) || (A1 < ϵ && A2 < ϵ && A3 < ϵ)
 end
 
-function ray_triangle_intersection(A, B, C, origin, dir, ϵ = 1e-6)
-    # See: https://www.iue.tuwien.ac.at/phd/ertl/node114.html
-    AO = A .- origin
-    BO = B .- origin
-    CO = C .- origin
-    A1 = 0.5 * dot(cross(BO, CO), dir)
-    A2 = 0.5 * dot(cross(CO, AO), dir)
-    A3 = 0.5 * dot(cross(AO, BO), dir)
-
-    e = 1e-3
-    if (A1 > -ϵ && A2 > -ϵ && A3 > -ϵ) || (A1 < ϵ && A2 < ϵ && A3 < ϵ)
-        Point3f((A1 * A .+ A2 * B .+ A3 * C) / (A1 + A2 + A3))
-    else
-        Point3f(NaN)
-    end
-end
-
-
-### Surface positions
-########################################
-
-surface_x(xs::ClosedInterval, i, j, N) = minimum(xs) + (maximum(xs) - minimum(xs)) * (i-1) / (N-1)
-surface_x(xs, i, j, N) = xs[i]
-surface_x(xs::AbstractMatrix, i, j, N) = xs[i, j]
-
-surface_y(ys::ClosedInterval, i, j, N) = minimum(ys) + (maximum(ys) - minimum(ys)) * (j-1) / (N-1)
-surface_y(ys, i, j, N) = ys[j]
-surface_y(ys::AbstractMatrix, i, j, N) = ys[i, j]
-
-function surface_pos(xs, ys, zs, i, j)
-    N, M = size(zs)
-    Point3f(surface_x(xs, i, j, N), surface_y(ys, i, j, M), zs[i, j])
-end
 
 
 ### Mapping mesh vertex indices to Vector{Polygon} index
@@ -216,20 +136,20 @@ function point_in_quad_parameter(
 
     # Our initial guess is that P is in the center of the quad (in terms of AB and DC)
     f = 0.5
-
-    for i in 0:iterations
+    AB = B - A
+    DC = C - D
+    for _ in 0:iterations
         # vector between top and bottom point of the current line
         dir = (D + f * (C - D)) - (A + f * (B - A))
-        DC = C - D
-        AB = B - A
         # solves P + _ * dir = A + f1 * (B - A) (intersection point of ray & line)
         f1, _ = inv(Mat2f(AB..., dir...)) * (P - A)
         f2, _ = inv(Mat2f(DC..., dir...)) * (P - D)
 
         # next fraction estimate should be between f1 and f2
         # adding 2f to this helps avoid jumping between low and high values
+        old_f = f
         f = 0.25 * (2f + f1 + f2)
-        if abs(f2 - f1) < epsilon
+        if abs(old_f - f) < epsilon
             return f
         end
     end
@@ -241,11 +161,13 @@ end
 ## Shifted projection
 ########################################
 
-function shift_project(scene, plot, pos)
+@deprecate shift_project(scene, plot, pos) shift_project(scene, pos) false
+
+function shift_project(scene, pos)
     project(
         camera(scene).projectionview[],
         Vec2f(widths(pixelarea(scene)[])),
-        apply_transform(transform_func(plot), pos, to_value(get(plot, :space, :data)))
+        pos
     ) .+ Vec2f(origin(pixelarea(scene)[]))
 end
 
@@ -417,8 +339,8 @@ function show_data_recursion(inspector, plot, idx)
             show_data(inspector, plot, idx)
         end
 
-        if processed
-            inspector.selection = plot
+        if processed && inspector.selection != plot
+            clear_temporary_plots!(inspector, plot)
         end
 
         return processed
@@ -439,8 +361,8 @@ function show_data_recursion(inspector, plot::AbstractPlot, idx, source)
             show_data(inspector, plot, idx, source)
         end
 
-        if processed
-            inspector.selection = plot
+        if processed && inspector.selection != plot
+            clear_temporary_plots!(inspector, plot)
         end
 
         return processed
@@ -500,13 +422,14 @@ function show_data(inspector::DataInspector, plot::Scatter, idx)
     tt = inspector.plot
     scene = parent_scene(plot)
 
-    proj_pos = shift_project(scene, plot, to_ndim(Point3f, plot[1][][idx], 0))
+    pos = position_on_plot(plot, idx)
+    proj_pos = shift_project(scene, pos)
     update_tooltip_alignment!(inspector, proj_pos)
 
     if haskey(plot, :inspector_label)
-        tt.text[] = plot[:inspector_label][](plot, idx, plot[1][][idx])
+        tt.text[] = plot[:inspector_label][](plot, idx, pos)
     else
-        tt.text[] = position2string(plot[1][][idx])
+        tt.text[] = position2string(pos)
     end
     tt.offset[] = ifelse(
         a.apply_tooltip_offset[],
@@ -526,11 +449,9 @@ function show_data(inspector::DataInspector, plot::MeshScatter, idx)
     scene = parent_scene(plot)
 
     if a.enable_indicators[]
-        T = transformationmatrix(
-            plot[1][][idx],
-            _to_scale(plot.markersize[], idx),
-            _to_rotation(plot.rotations[], idx)
-        )
+        translation = apply_transform_and_model(plot, plot[1][][idx])
+        rotation = _to_rotation(plot.rotations[], idx)
+        scale = _to_scale(plot.markersize[], idx)
 
         if inspector.selection != plot
             clear_temporary_plots!(inspector, plot)
@@ -545,9 +466,12 @@ function show_data(inspector::DataInspector, plot::MeshScatter, idx)
             bbox = Rect{3, Float32}(convert_attribute(
                 plot.marker[], Key{:marker}(), Key{Makie.plotkey(plot)}()
             ))
+            T = Transformation(
+                identity; translation = translation, rotation = rotation, scale = scale
+            )
 
             p = wireframe!(
-                scene, bbox, model = T, color = a.indicator_color,
+                scene, bbox, transformation = T, color = a.indicator_color,
                 linewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
                 visible = a.indicator_visible, inspectable = false
             )
@@ -558,21 +482,21 @@ function show_data(inspector::DataInspector, plot::MeshScatter, idx)
 
         elseif !isempty(inspector.temp_plots)
             p = inspector.temp_plots[1]
-            p.model[] = T
-
+            transform!(p, translation = translation, scale = scale, rotation = rotation)
         end
 
 
         a.indicator_visible[] = true
     end
 
-    proj_pos = shift_project(scene, plot, to_ndim(Point3f, plot[1][][idx], 0))
+    pos = position_on_plot(plot, idx)
+    proj_pos = shift_project(scene, pos)
     update_tooltip_alignment!(inspector, proj_pos)
 
     if haskey(plot, :inspector_label)
-        tt.text[] = plot[:inspector_label][](plot, idx, plot[1][][idx])
+        tt.text[] = plot[:inspector_label][](plot, idx, pos)
     else
-        tt.text[] = position2string(plot[1][][idx])
+        tt.text[] = position2string(pos)
     end
     tt.visible[] = true
 
@@ -586,11 +510,9 @@ function show_data(inspector::DataInspector, plot::Union{Lines, LineSegments}, i
     scene = parent_scene(plot)
 
     # cast ray from cursor into screen, find closest point to line
-    p0, p1 = plot[1][][idx-1:idx]
-    origin, dir = view_ray(scene)
-    pos = closest_point_on_line(p0, p1, origin, dir)
+    pos = position_on_plot(plot, idx)
 
-    proj_pos = shift_project(scene, plot, to_ndim(Point3f, pos, 0))
+    proj_pos = shift_project(scene, pos)
     update_tooltip_alignment!(inspector, proj_pos)
 
     tt.offset[] = ifelse(
@@ -600,9 +522,9 @@ function show_data(inspector::DataInspector, plot::Union{Lines, LineSegments}, i
     )
 
     if haskey(plot, :inspector_label)
-        tt.text[] = plot[:inspector_label][](plot, idx, typeof(p0)(pos))
+        tt.text[] = plot[:inspector_label][](plot, idx, eltype(plot[1][])(pos))
     else
-        tt.text[] = position2string(typeof(p0)(pos))
+        tt.text[] = position2string(eltype(plot[1][])(pos))
     end
     tt.visible[] = true
     a.indicator_visible[] && (a.indicator_visible[] = false)
@@ -616,7 +538,14 @@ function show_data(inspector::DataInspector, plot::Mesh, idx)
     tt = inspector.plot
     scene = parent_scene(plot)
 
-    bbox = boundingbox(plot)
+    # Manual boundingbox including transfunc
+    bbox = let
+        points = point_iterator(plot)
+        trans_func = transform_func(plot)
+        model = plot.model[]
+        iter = iterate_transformed(points, model, to_value(get(plot, :space, :data)), trans_func)
+        limits_from_transformed_points(iter)
+    end
     proj_pos = Point2f(mouseposition_px(inspector.root))
     update_tooltip_alignment!(inspector, proj_pos)
 
@@ -632,7 +561,8 @@ function show_data(inspector::DataInspector, plot::Mesh, idx)
             end
 
             p = wireframe!(
-                scene, bbox, color = a.indicator_color,
+                scene, bbox, color = a.indicator_color, 
+                transformation = Transformation(),
                 linewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
                 visible = a.indicator_visible, inspectable = false
             )
@@ -663,43 +593,11 @@ end
 function show_data(inspector::DataInspector, plot::Surface, idx)
     a = inspector.attributes
     tt = inspector.plot
-    scene = parent_scene(plot)
 
     proj_pos = Point2f(mouseposition_px(inspector.root))
     update_tooltip_alignment!(inspector, proj_pos)
 
-    xs = plot[1][]
-    ys = plot[2][]
-    zs = plot[3][]
-    w, h = size(zs)
-    _i = mod1(idx, w); _j = div(idx-1, w)
-
-    # This isn't the most accurate so we include some neighboring faces
-    origin, dir = view_ray(scene)
-    pos = Point3f(NaN)
-    for i in _i-1:_i+1, j in _j-1:_j+1
-        (1 <= i <= w) && (1 <= j < h) || continue
-
-        if i - 1 > 0
-            pos = ray_triangle_intersection(
-                surface_pos(xs, ys, zs, i, j),
-                surface_pos(xs, ys, zs, i-1, j),
-                surface_pos(xs, ys, zs, i, j+1),
-                origin, dir
-            )
-        end
-
-        if i + 1 <= w && isnan(pos)
-            pos = ray_triangle_intersection(
-                surface_pos(xs, ys, zs, i, j),
-                surface_pos(xs, ys, zs, i, j+1),
-                surface_pos(xs, ys, zs, i+1, j+1),
-                origin, dir
-            )
-        end
-
-        isnan(pos) || break
-    end
+    pos = position_on_plot(plot, idx)
 
     if !isnan(pos)
         tt[1][] = proj_pos
@@ -730,20 +628,22 @@ function show_imagelike(inspector, plot, name, edge_based)
     a = inspector.attributes
     tt = inspector.plot
     scene = parent_scene(plot)
-    mpos = mouseposition(scene)
 
-    if plot.interpolate[]
-        i, j, z = _interpolated_getindex(plot[1][], plot[2][], plot[3][], mpos)
-        x, y = mpos
-    else
-        i, j, z = _pixelated_getindex(plot[1][], plot[2][], plot[3][], mpos, edge_based)
-        x = i; y = j
+    pos = position_on_plot(plot, -1, apply_transform = false)[Vec(1, 2)] # index irrelevant
+
+    # Not on image/heatmap
+    if isnan(pos)
+        a.indicator_visible[] = false
+        tt.visible[] = false
+        return true
     end
 
-    if haskey(plot, :inspector_label)
-        tt.text[] = plot[:inspector_label][](plot, (i, j), Point3f(mpos[1], mpos[2], z))
+    if plot.interpolate[]
+        i, j, z = _interpolated_getindex(plot[1][], plot[2][], plot[3][], pos)
+        x, y = pos
     else
-        tt.text[] = color2text(name, x, y, z)
+        i, j, z = _pixelated_getindex(plot[1][], plot[2][], plot[3][], pos, edge_based)
+        x = i; y = j
     end
 
     # in case we hover over NaN values
@@ -753,6 +653,12 @@ function show_imagelike(inspector, plot, name, edge_based)
         return true
     end
 
+    if haskey(plot, :inspector_label)
+        tt.text[] = plot[:inspector_label][](plot, (i, j), Point3f(pos[1], pos[2], z))
+    else
+        tt.text[] = color2text(name, x, y, z)
+    end
+
     a._color[] = if z isa AbstractFloat
         interpolated_getindex(
             to_colormap(plot.colormap[]), z,
@@ -762,41 +668,44 @@ function show_imagelike(inspector, plot, name, edge_based)
         z
     end
 
-    position = to_ndim(Point3f, mpos, 0)
     proj_pos = Point2f(mouseposition_px(inspector.root))
     update_tooltip_alignment!(inspector, proj_pos)
 
     if a.enable_indicators[]
         if plot.interpolate[]
-            if inspector.selection != plot
+            if inspector.selection != plot || (length(inspector.temp_plots) != 1) || 
+                    !(inspector.temp_plots[1] isa Scatter)
                 clear_temporary_plots!(inspector, plot)
                 p = scatter!(
-                    scene, position, color = a._color,
+                    scene, pos, color = a._color,
                     visible = a.indicator_visible,
-                    inspectable = false,
-                    marker=:rect, markersize = map(r -> 3r, a.range),
+                    inspectable = false, model = plot.model,
+                    # TODO switch to Rect with 2r-1 or 2r-2 markersize to have 
+                    # just enough space to always detect the underlying image
+                    marker=:rect, markersize = map(r -> 2r, a.range), 
                     strokecolor = a.indicator_color,
-                    strokewidth = a.indicator_linewidth
+                    strokewidth = a.indicator_linewidth,
+                    depth_shift = -1f-3
                 )
-                translate!(p, Vec3f(0, 0, a.depth[]-1))
                 push!(inspector.temp_plots, p)
-            elseif !isempty(inspector.temp_plots)
+            else
                 p = inspector.temp_plots[1]
-                p[1].val[1] = position
+                p[1].val[1] = pos
                 notify(p[1])
             end
         else
             bbox = _pixelated_image_bbox(plot[1][], plot[2][], plot[3][], i, j, edge_based)
-            if inspector.selection != plot
+            if inspector.selection != plot || (length(inspector.temp_plots) != 1) || 
+                    !(inspector.temp_plots[1] isa Wireframe)
                 clear_temporary_plots!(inspector, plot)
                 p = wireframe!(
-                    scene, bbox, color = a.indicator_color,
+                    scene, bbox, color = a.indicator_color, model = plot.model,
                     strokewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
-                    visible = a.indicator_visible, inspectable = false
+                    visible = a.indicator_visible, inspectable = false,
+                    depth_shift = -1f-3
                 )
-                translate!(p, Vec3f(0, 0, a.depth[]-1))
                 push!(inspector.temp_plots, p)
-            elseif !isempty(inspector.temp_plots)
+            else
                 p = inspector.temp_plots[1]
                 p[1][] = bbox
             end
@@ -902,8 +811,8 @@ function show_data(inspector::DataInspector, plot::BarPlot, idx)
     tt = inspector.plot
     scene = parent_scene(plot)
 
-    pos = plot[1][][idx]
-    proj_pos = shift_project(scene, plot, to_ndim(Point3f, pos, 0))
+    pos = apply_transform_and_model(plot, plot[1][][idx])
+    proj_pos = shift_project(scene, to_ndim(Point3f, pos, 0))
     update_tooltip_alignment!(inspector, proj_pos)
 
     if a.enable_indicators[]
@@ -942,9 +851,10 @@ function show_data(inspector::DataInspector, plot::Arrows, idx, ::LineSegments)
     return show_data(inspector, plot, div(idx+1, 2), nothing)
 end
 function show_data(inspector::DataInspector, plot::Arrows, idx, source)
-    a = inspector.plot.attributes
+    a = inspector.attributes
     tt = inspector.plot
-    pos = plot[1][][idx]
+    pos = apply_transform_and_model(plot, plot[1][][idx])
+
     mpos = Point2f(mouseposition_px(inspector.root))
     update_tooltip_alignment!(inspector, mpos)
 
@@ -953,7 +863,7 @@ function show_data(inspector::DataInspector, plot::Arrows, idx, source)
 
     tt[1][] = mpos
     if haskey(plot, :inspector_label)
-        tt.text[] = plot[:inspector_label][](plot, idx, mpos)
+        tt.text[] = plot[:inspector_label][](plot, idx, pos)
     else
         tt.text[] = "Position:\n  $p\nDirection:\n  $v"
     end
@@ -967,7 +877,7 @@ end
 # backend handle picking colors from a colormap
 function show_data(inspector::DataInspector, plot::Contourf, idx, source::Mesh)
     tt = inspector.plot
-    idx = show_poly(inspector, plot.plots[1], idx, source)
+    idx = show_poly(inspector, plot, plot.plots[1], idx, source)
     level = plot.plots[1].color[][idx]
 
     mpos = Point2f(mouseposition_px(inspector.root))
@@ -993,14 +903,13 @@ end
 #     return true
 # end
 
-function show_poly(inspector, plot, idx, source)
+function show_poly(inspector, plot, poly, idx, source)
     a = inspector.attributes
-    idx = vertexindex2poly(plot[1][], idx)
+    idx = vertexindex2poly(poly[1][], idx)
 
     if a.enable_indicators[]
-
-        line_collection = copy(convert_arguments(PointBased(), plot[1][][idx].exterior)[1])
-        for int in plot[1][][idx].interiors
+        line_collection = copy(convert_arguments(PointBased(), poly[1][][idx].exterior)[1])
+        for int in poly[1][][idx].interiors
             push!(line_collection, Point2f(NaN))
             append!(line_collection, convert_arguments(PointBased(), int)[1])
         end
@@ -1010,11 +919,11 @@ function show_poly(inspector, plot, idx, source)
             clear_temporary_plots!(inspector, plot)
 
             p = lines!(
-                scene, line_collection, color = a.indicator_color,
+                scene, line_collection, color = a.indicator_color, 
+                transformation = Transformation(source),
                 strokewidth = a.indicator_linewidth, linestyle = a.indicator_linestyle,
-                visible = a.indicator_visible, inspectable = false
+                visible = a.indicator_visible, inspectable = false, depth_shift = -1f-3
             )
-            translate!(p, Vec3f(0, 0, a.depth[]-1))
             push!(inspector.temp_plots, p)
 
         elseif !isempty(inspector.temp_plots)
@@ -1030,68 +939,46 @@ end
 function show_data(inspector::DataInspector, plot::VolumeSlices, idx, child::Heatmap)
     a = inspector.attributes
     tt = inspector.plot
-    scene = parent_scene(plot)
 
-    proj_pos = Point2f(mouseposition_px(inspector.root))
-    update_tooltip_alignment!(inspector, proj_pos)
+    pos = position_on_plot(child, -1, apply_transform = false)[Vec(1, 2)] # index irrelevant
 
-    qs = extrema(child[1][])
-    ps = extrema(child[2][])
-    data = child[3][]
-    T = child.transformation.model[]
-
-    vs = [ # clockwise
-        Point3f(T * Point4f(qs[1], ps[1], 0, 1)),
-        Point3f(T * Point4f(qs[1], ps[2], 0, 1)),
-        Point3f(T * Point4f(qs[2], ps[2], 0, 1)),
-        Point3f(T * Point4f(qs[2], ps[1], 0, 1))
-    ]
-
-    origin, dir = view_ray(scene)
-    pos = Point3f(NaN)
-    pos = ray_triangle_intersection(vs[1], vs[2], vs[3], origin, dir)
+    # Not on heatmap
     if isnan(pos)
-        pos = ray_triangle_intersection(vs[3], vs[4], vs[1], origin, dir)
+        a.indicator_visible[] && (a.indicator_visible[] = false)
+        tt.visible[] = false
+        return true
     end
 
-    if !isnan(pos)
-        child_idx = findfirst(isequal(child), plot.plots)
-        if child_idx == 2
-            x = pos[2]; y = pos[3]
-        elseif child_idx == 3
-            x = pos[1]; y = pos[3]
-        else
-            x = pos[1]; y = pos[2]
-        end
-        i = clamp(round(Int, (x - qs[1]) / (qs[2] - qs[1]) * size(data, 1) + 0.5), 1, size(data, 1))
-        j = clamp(round(Int, (y - ps[1]) / (ps[2] - ps[1]) * size(data, 2) + 0.5), 1, size(data, 2))
-        val = data[i, j]
+    i, j, val = _pixelated_getindex(child[1][], child[2][], child[3][], pos, true)
 
-        tt[1][] = proj_pos
-        if haskey(plot, :inspector_label)
-            tt.text[] = plot[:inspector_label][](plot, (i, j), pos)
-        else
-            tt.text[] = @sprintf(
-                "x: %0.6f\ny: %0.6f\nz: %0.6f\n%0.6f0",
-                pos[1], pos[2], pos[3], val
-            )
-        end
-        tt.visible[] = true
+    proj_pos = Point2f(mouseposition_px(inspector.root))
+    update_tooltip_alignment!(inspector, proj_pos)
+    tt[1][] = proj_pos
+    
+    world_pos = apply_transform_and_model(child, pos)
+
+    if haskey(plot, :inspector_label)
+        tt.text[] = plot[:inspector_label][](plot, (i, j), world_pos)
     else
-        tt.visible[] = false
+        tt.text[] = @sprintf(
+            "x: %0.6f\ny: %0.6f\nz: %0.6f\n%0.6f0",
+            world_pos[1], world_pos[2], world_pos[3], val
+        )
     end
+
+    tt.visible[] = true
     a.indicator_visible[] && (a.indicator_visible[] = false)
 
     return true
 end
 
 
-function show_data(inspector::DataInspector, plot::Band, ::Integer, ::Mesh)
+function show_data(inspector::DataInspector, plot::Band, idx::Integer, mesh::Mesh)
     scene = parent_scene(plot)
     tt = inspector.plot
     a = inspector.attributes
 
-    pos = Point2f(mouseposition(scene))
+    pos = Point2f(position_on_plot(mesh, idx, apply_transform = false)) #Point2f(mouseposition(scene))
     ps1 = plot.converted[1][]
     ps2 = plot.converted[2][]
 
@@ -1112,22 +999,20 @@ function show_data(inspector::DataInspector, plot::Band, ::Integer, ::Mesh)
 
         # Draw the line
         if a.enable_indicators[]
-            model = plot.model[]
-
-            if inspector.selection != plot || isempty(inspector.temp_plots)
+            # Why does this sometimes create 2+ plots
+            if inspector.selection != plot || (length(inspector.temp_plots) != 1)
                 clear_temporary_plots!(inspector, plot)
                 p = lines!(
-                    scene, [P1, P2], model = model, 
+                    scene, [P1, P2], transformation = Transformation(plot.transformation),
                     color = a.indicator_color, strokewidth = a.indicator_linewidth, 
                     linestyle = a.indicator_linestyle,
-                    visible = a.indicator_visible, inspectable = false
+                    visible = a.indicator_visible, inspectable = false,
+                    depth_shift = -1f-3
                 )
-                translate!(p, Vec3f(0, 0, a.depth[]))
                 push!(inspector.temp_plots, p)
             elseif !isempty(inspector.temp_plots)
                 p = inspector.temp_plots[1]
                 p[1][] = [P1, P2]
-                p.model[] = model
             end
 
             a.indicator_visible[] = true
@@ -1139,6 +1024,8 @@ function show_data(inspector::DataInspector, plot::Band, ::Integer, ::Mesh)
         if haskey(plot, :inspector_label)
             tt.text[] = plot[:inspector_label][](plot, right, (P1, P2))
         else
+            P1 = apply_transform_and_model(mesh, P1, Point2f)
+            P2 = apply_transform_and_model(mesh, P2, Point2f)
             tt.text[] = @sprintf("(%0.3f, %0.3f) .. (%0.3f, %0.3f)", P1[1], P1[2], P2[1], P2[2])
         end
         tt.visible[] = true
diff --git a/src/interaction/ray_casting.jl b/src/interaction/ray_casting.jl
new file mode 100644
index 00000000000..f87ffe4bae3
--- /dev/null
+++ b/src/interaction/ray_casting.jl
@@ -0,0 +1,390 @@
+################################################################################
+### Ray Generation
+################################################################################
+
+struct Ray
+    origin::Point3f
+    direction::Vec3f
+end
+
+"""
+    ray_at_cursor(fig/ax/scene)
+
+Returns a Ray into the scene starting at the current cursor position.
+"""
+ray_at_cursor(x) = ray_at_cursor(get_scene(x))
+function ray_at_cursor(scene::Scene)
+    return Ray(scene, mouseposition_px(scene))
+end
+
+"""
+    Ray(scene[, cam = cameracontrols(scene)], xy)
+
+Returns a `Ray` into the given `scene` passing through pixel position `xy`. Note
+that the pixel position should be relative to the origin of the scene, as it is 
+when calling `mouseposition_px(scene)`.
+"""
+Ray(scene::Scene, xy::VecTypes{2}) = Ray(scene, cameracontrols(scene), xy)
+
+
+function Ray(scene::Scene, cam::Camera3D, xy::VecTypes{2})
+    lookat = cam.lookat[]
+    eyepos = cam.eyeposition[]
+    viewdir = lookat - eyepos
+
+    u_z = normalize(viewdir)
+    u_x = normalize(cross(u_z, cam.upvector[]))
+    u_y = normalize(cross(u_x, u_z))
+    
+    px_width, px_height = widths(scene.px_area[])
+    aspect = px_width / px_height
+    rel_pos = 2 .* xy ./ (px_width, px_height) .- 1
+
+    if cam.attributes.projectiontype[] === Perspective
+        dir = (rel_pos[1] * aspect * u_x + rel_pos[2] * u_y) * tand(0.5 * cam.fov[]) + u_z
+        return Ray(cam.eyeposition[], normalize(dir))
+    else
+        # Orthographic has consistent direction, but not starting point
+        origin = norm(viewdir) * (rel_pos[1] * aspect * u_x + rel_pos[2] * u_y)
+        return Ray(origin, normalize(viewdir))
+    end
+end
+
+function Ray(scene::Scene, cam::Camera2D, xy::VecTypes{2})
+    rel_pos = xy ./ widths(scene.px_area[])
+    pv = scene.camera.projectionview[]
+    m = Vec2f(pv[1, 1], pv[2, 2])
+    b = Vec2f(pv[1, 4], pv[2, 4])
+    origin = (2 * rel_pos .- 1 - b) ./ m
+    return Ray(to_ndim(Point3f, origin, 10_000f0), Vec3f(0,0,-1))
+end
+
+function Ray(::Scene, ::PixelCamera, xy::VecTypes{2})
+    return Ray(to_ndim(Point3f, xy, 10_000f0), Vec3f(0,0,-1))
+end
+
+function Ray(scene::Scene, ::RelativeCamera, xy::VecTypes{2})
+    origin = xy ./ widths(scene.px_area[])
+    return Ray(to_ndim(Point3f, origin, 10_000f0), Vec3f(0,0,-1))
+end
+
+Ray(scene::Scene, cam, xy::VecTypes{2}) = ray_from_projectionview(scene, xy)
+
+# This method should always work 
+function ray_from_projectionview(scene::Scene, xy::VecTypes{2})
+    inv_view_proj = inv(camera(scene).projectionview[])
+    area = pixelarea(scene)[]
+
+    # This figures out the camera view direction from the projectionview matrix
+    # and computes a ray from a near and a far point.
+    # Based on ComputeCameraRay from ImGuizmo
+    mp = 2f0 .* xy ./ widths(area) .- 1f0
+    v = inv_view_proj * Vec4f(0, 0, -10, 1)
+    reversed = v[3] < v[4]
+    near = reversed ? 1f0 - 1e-6 : 0f0
+    far = reversed ? 0f0 : 1f0 - 1e-6
+
+    origin = inv_view_proj * Vec4f(mp[1], mp[2], near, 1f0)
+    origin = origin[Vec(1, 2, 3)] ./ origin[4]
+
+    p = inv_view_proj * Vec4f(mp[1], mp[2], far, 1f0)
+    p = p[Vec(1, 2, 3)] ./ p[4]
+
+    dir = normalize(p .- origin)
+
+    return Ray(origin, dir)
+end
+
+
+function transform(M::Mat4f, ray::Ray)
+    p4d = M * to_ndim(Point4f, ray.origin, 1f0)
+    dir = normalize(M[Vec(1,2,3), Vec(1,2,3)] * ray.direction)
+    return Ray(p4d[Vec(1,2,3)] / p4d[4], dir)
+end
+
+
+################################################################################
+### Ray - object intersections
+################################################################################
+
+
+# These work in 2D and 3D
+function closest_point_on_line(A::VecTypes, B::VecTypes, ray::Ray)
+    return closest_point_on_line(to_ndim(Point3f, A, 0), to_ndim(Point3f, B, 0), ray)
+end
+function closest_point_on_line(A::Point3f, B::Point3f, ray::Ray)
+    # See:
+    # https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection
+    AB_norm = norm(B .- A)
+    u_AB = (B .- A) / AB_norm
+    u_perp = normalize(cross(ray.direction, u_AB))
+    # e_RD, e_perp defines a plane with normal n
+    n = normalize(cross(ray.direction, u_perp))
+    t = dot(ray.origin .- A, n) / dot(u_AB, n)
+    return A .+ clamp(t, 0.0, AB_norm) * u_AB
+end
+
+
+function ray_triangle_intersection(A::VecTypes{3}, B::VecTypes{3}, C::VecTypes{3}, ray::Ray, ϵ = 1e-6)
+    # See: https://www.iue.tuwien.ac.at/phd/ertl/node114.html
+    # Alternative: https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
+    AO = A .- ray.origin
+    BO = B .- ray.origin
+    CO = C .- ray.origin
+    A1 = 0.5 * dot(cross(BO, CO), ray.direction)
+    A2 = 0.5 * dot(cross(CO, AO), ray.direction)
+    A3 = 0.5 * dot(cross(AO, BO), ray.direction)
+
+    # all positive or all negative
+    if (A1 > -ϵ && A2 > -ϵ && A3 > -ϵ) || (A1 < ϵ && A2 < ϵ && A3 < ϵ)
+        return Point3f((A1 * A .+ A2 * B .+ A3 * C) / (A1 + A2 + A3))
+    else
+        return Point3f(NaN)
+    end
+end
+
+function ray_rect_intersection(rect::Rect2f, ray::Ray)
+    possible_hit = ray.origin - ray.origin[3] / ray.direction[3] * ray.direction
+    min = minimum(rect); max = maximum(rect)
+    if all(min <= possible_hit[Vec(1,2)] <= max)
+        return possible_hit
+    end
+    return Point3f(NaN)
+end
+
+
+function ray_rect_intersection(rect::Rect3f, ray::Ray)
+    mins = (minimum(rect) - ray.origin) ./ ray.direction
+    maxs = (maximum(rect) - ray.origin) ./ ray.direction
+    x, y, z = min.(mins, maxs)
+    possible_hit = max(x, y, z)
+    if possible_hit < minimum(max.(mins, maxs))
+        return ray.origin + possible_hit * ray.direction
+    end
+    return Point3f(NaN)
+end
+
+function is_point_on_ray(p::Point3f, ray::Ray)
+    diff = ray.origin - p
+    return abs(dot(diff, ray.direction)) ≈ abs(norm(diff))
+end
+
+
+################################################################################
+### Ray casting (positions from ray-plot intersections)
+################################################################################
+
+
+"""
+    ray_assisted_pick(fig/ax/scene[, xy = events(fig/ax/scene).mouseposition[], apply_transform = true])
+
+This function performs a `pick` at the given pixel position `xy` and returns the
+picked `plot`, `index` and world or input space `position::Point3f`. It is equivalent to
+```
+plot, idx = pick(fig/ax/scene, xy)
+ray = Ray(parent_scene(plot), xy .- minimum(pixelarea(parent_scene(plot))[]))
+position = position_on_plot(plot, idx, ray, apply_transform = true)
+```
+See [`position_on_plot`](@ref) for more information.
+"""
+function ray_assisted_pick(obj, xy = events(obj).mouseposition[]; apply_transform = true)
+    plot, idx = pick(get_scene(obj), xy)
+    isnothing(plot) && return (plot, idx, Point3f(NaN))
+    scene = parent_scene(plot)
+    ray = Ray(scene, xy .- minimum(pixelarea(scene)[]))
+    pos = position_on_plot(plot, idx, ray, apply_transform = apply_transform)
+    return (plot, idx, pos)
+end
+
+
+"""
+    position_on_plot(plot, index[, ray::Ray; apply_transform = true])
+
+This function calculates the world or input space position of a ray - plot 
+intersection with the result `plot, idx = pick(...)` and a ray cast from the 
+picked position. If there is no intersection `Point3f(NaN)` will be returned.
+
+This should be called as
+```
+plot, idx = pick(ax, px_pos)
+pos_in_ax = position_on_plot(plot, idx, Ray(ax, px_pos .- minimum(pixelarea(ax.scene)[])))
+```
+or more simply `plot, idx, pos_in_ax = ray_assisted_pick(ax, px_pos)`.
+
+You can switch between getting a position in world space (after applying 
+transformations like `log`, `translate!()`, `rotate!()` and `scale!()`) and 
+input space (the raw position data of the plot) by adjusting `apply_transform`.
+
+Note that `position_on_plot` is only implemented for primitive plot types, i.e. 
+the  possible return types of `pick`. Depending on the plot type the calculation 
+differs:
+- `scatter` and `meshscatter` return the position of the picked marker/mesh
+- `text` is excluded, always returning `Point3f(NaN)`
+- `volume` calculates the ray - rect intersection for its bounding box
+- `lines` and `linesegments` return the closest point on the line to the ray
+- `mesh` and `surface` check for ray-triangle intersections for every triangle containing the picked vertex
+- `image` and `heatmap` check for ray-rect intersection
+"""
+function position_on_plot(plot::AbstractPlot, idx::Integer; apply_transform = true)
+    return position_on_plot(
+        plot, idx, ray_at_cursor(parent_scene(plot)); 
+        apply_transform = apply_transform
+    )
+end
+
+
+function position_on_plot(plot::Union{Scatter, MeshScatter}, idx, ray::Ray; apply_transform = true)
+    pos = to_ndim(Point3f, plot[1][][idx], 0f0)
+    if apply_transform && !isnan(pos)
+        return apply_transform_and_model(plot, pos)
+    else
+        return pos
+    end
+end
+
+function position_on_plot(plot::Union{Lines, LineSegments}, idx, ray::Ray; apply_transform = true)
+    p0, p1 = apply_transform_and_model(plot, plot[1][][idx-1:idx])
+
+    pos = closest_point_on_line(p0, p1, ray)
+    
+    if apply_transform
+        return pos
+    else
+        p4d = inv(plot.model[]) * to_ndim(Point4f, pos, 1f0)
+        p3d = p4d[Vec(1, 2, 3)] / p4d[4]
+        itf = inverse_transform(transform_func(plot))
+        return Makie.apply_transform(itf, p3d, get(plot, :space, :data))
+    end
+end
+
+function position_on_plot(plot::Union{Heatmap, Image}, idx, ray::Ray; apply_transform = true)
+    # Heatmap and Image are always a Rect2f. The transform function is currently 
+    # not allowed to change this, so applying it should be fine. Applying the 
+    # model matrix may add a z component to the Rect2f, which we can't represent.
+    # So we instead inverse-transform the ray
+    space = to_value(get(plot, :space, :data))
+    p0, p1 = map(Point2f.(extrema(plot.x[]), extrema(plot.y[]))) do p
+        return Makie.apply_transform(transform_func(plot), p, space)
+    end
+    ray = transform(inv(plot.model[]), ray)
+    pos = ray_rect_intersection(Rect2f(p0, p1 - p0), ray)
+    
+    if apply_transform
+        p4d = plot.model[] * to_ndim(Point4f, to_ndim(Point3f, pos, 0), 1)
+        return p4d[Vec(1, 2, 3)] / p4d[4]
+    else
+        pos = Makie.apply_transform(inverse_transform(transform_func(plot)), pos, space)
+        return to_ndim(Point3f, pos, 0)
+    end
+end
+
+function position_on_plot(plot::Mesh, idx, ray::Ray; apply_transform = true)
+    positions = coordinates(plot.mesh[])
+    ray = transform(inv(plot.model[]), ray)
+    tf = transform_func(plot)
+    space = to_value(get(plot, :space, :data))
+
+    for f in faces(plot.mesh[])
+        if idx in f
+            p1, p2, p3 = positions[f]
+            p1, p2, p3 = Makie.apply_transform.(tf, (p1, p2, p3), space)
+            pos = ray_triangle_intersection(p1, p2, p3, ray)
+            if pos !== Point3f(NaN)
+                if apply_transform
+                    p4d = plot.model[] * to_ndim(Point4f, pos, 1)
+                    return Point3f(p4d) / p4d[4]
+                else
+                    return Makie.apply_transform(inverse_transform(tf), pos, space)
+                end
+            end
+        end
+    end
+
+    @info "Did not find $idx"
+
+    return Point3f(NaN)
+end
+
+# Handling indexing into different surface input types
+surface_x(xs::ClosedInterval, i, j, N) = minimum(xs) + (maximum(xs) - minimum(xs)) * (i-1) / (N-1)
+surface_x(xs, i, j, N) = xs[i]
+surface_x(xs::AbstractMatrix, i, j, N) = xs[i, j]
+
+surface_y(ys::ClosedInterval, i, j, N) = minimum(ys) + (maximum(ys) - minimum(ys)) * (j-1) / (N-1)
+surface_y(ys, i, j, N) = ys[j]
+surface_y(ys::AbstractMatrix, i, j, N) = ys[i, j]
+
+function surface_pos(xs, ys, zs, i, j)
+    N, M = size(zs)
+    return Point3f(surface_x(xs, i, j, N), surface_y(ys, i, j, M), zs[i, j])
+end
+
+function position_on_plot(plot::Surface, idx, ray::Ray; apply_transform = true)
+    xs = plot[1][]
+    ys = plot[2][]
+    zs = plot[3][]
+    w, h = size(zs)
+    _i = mod1(idx, w); _j = div(idx-1, w)
+
+    ray = transform(inv(plot.model[]), ray)
+    tf = transform_func(plot)
+    space = to_value(get(plot, :space, :data))
+    
+    # This isn't the most accurate so we include some neighboring faces
+    pos = Point3f(NaN)
+    for i in _i-1:_i+1, j in _j-1:_j+1
+        (1 <= i <= w) && (1 <= j < h) || continue
+
+        if i - 1 > 0
+            # transforms only apply to x and y coordinates of surfaces
+            A = surface_pos(xs, ys, zs, i,   j)
+            B = surface_pos(xs, ys, zs, i-1, j)
+            C = surface_pos(xs, ys, zs, i, j+1)
+            A, B, C = map((A, B, C)) do p
+                xy = Makie.apply_transform(tf, Point2f(p), space)
+                Point3f(xy[1], xy[2], p[3])
+            end
+            pos = ray_triangle_intersection(A, B, C, ray)
+        end
+
+        if i + 1 <= w && isnan(pos)
+            A = surface_pos(xs, ys, zs, i,   j)
+            B = surface_pos(xs, ys, zs, i,   j+1)
+            C = surface_pos(xs, ys, zs, i+1, j+1)
+            A, B, C = map((A, B, C)) do p
+                xy = Makie.apply_transform(tf, Point2f(p), space)
+                Point3f(xy[1], xy[2], p[3])
+            end
+            pos = ray_triangle_intersection(A, B, C, ray)
+        end
+
+        isnan(pos) || break
+    end
+
+    if apply_transform
+        p4d = plot.model[] * to_ndim(Point4f, pos, 1)
+        return p4d[Vec(1, 2, 3)] / p4d[4]
+    else
+        xy = Makie.apply_transform(inverse_transform(tf), Point2f(pos), space)
+        return Point3f(xy[1], xy[2], pos[3])
+    end
+end
+
+function position_on_plot(plot::Volume, idx, ray::Ray; apply_transform = true)
+    min, max = Point3f.(extrema(plot.x[]), extrema(plot.y[]), extrema(plot.z[]))
+
+    if apply_transform
+        min = apply_transform_and_model(plot, min)
+        max = apply_transform_and_model(plot, max)
+        return ray_rect_intersection(Rect3f(min, max .- min), ray)
+    else
+        min = Makie.apply_transform(transform_func(plot), min, get(plot, :space, :data))
+        max = Makie.apply_transform(transform_func(plot), max, get(plot, :space, :data))
+        ray = transform(inv(plot.model[]), ray)
+        pos = ray_rect_intersection(Rect3f(min, max .- min), ray)
+        return Makie.apply_transform(inverse_transform(plot), pos, get(plot, :space, :data))
+    end
+end
+
+position_on_plot(plot::Text, args...; kwargs...) = Point3f(NaN)
+position_on_plot(plot::Nothing, args...; kwargs...) = Point3f(NaN)
\ No newline at end of file
diff --git a/src/layouting/transformation.jl b/src/layouting/transformation.jl
index 9454f8b50fb..6ddc7cd53d0 100644
--- a/src/layouting/transformation.jl
+++ b/src/layouting/transformation.jl
@@ -205,6 +205,34 @@ transformation(x::Attributes) = x.transformation[]
 transform_func(x) = transform_func_obs(x)[]
 transform_func_obs(x) = transformation(x).transform_func
 
+"""
+    apply_transform_and_model(plot, pos, output_type = Point3f)
+    apply_transform_and_model(model, transfrom_func, pos, output_type = Point3f)
+
+
+Applies the transform function and model matrix (i.e. transformations from 
+`translate!`, `rotate!` and `scale!`) to the given input
+"""
+function apply_transform_and_model(plot::AbstractPlot, pos, output_type = Point3f)
+    return apply_transform_and_model(
+        plot.model[], transform_func(plot), pos, 
+        to_value(get(plot, :space, :data)), 
+        output_type
+    )
+end
+function apply_transform_and_model(model::Mat4f, f, pos::VecTypes, space = :data, output_type = Point3f)
+    transformed = apply_transform(f, pos, space)
+    p4d = to_ndim(Point4f, to_ndim(Point3f, transformed, 0), 1)
+    p4d = model * p4d
+    p4d = p4d ./ p4d[4]
+    return to_ndim(output_type, p4d, NaN)
+end
+function apply_transform_and_model(model::Mat4f, f, positions::Vector, space = :data, output_type = Point3f)
+    return map(positions) do pos
+        apply_transform_and_model(model, f, pos, space, output_type)
+    end
+end
+
 """
     apply_transform(f, data, space)
 Apply the data transform func to the data if the space matches one
@@ -268,7 +296,7 @@ function apply_transform(f::PointTrans{N1}, point::Point{N2}) where {N1, N2}
 end
 
 function apply_transform(f, data::AbstractArray)
-    map(point-> apply_transform(f, point), data)
+    map(point -> apply_transform(f, point), data)
 end
 
 function apply_transform(f::Tuple{Any, Any}, point::VecTypes{2})
diff --git a/test/ray_casting.jl b/test/ray_casting.jl
new file mode 100644
index 00000000000..ec88508ba18
--- /dev/null
+++ b/test/ray_casting.jl
@@ -0,0 +1,163 @@
+@testset "Ray Casting" begin
+    @testset "View Rays" begin
+        scene = Scene()
+        xy = 0.5 * widths(pixelarea(scene)[])
+
+        orthographic_cam3d!(x) = cam3d!(x, perspectiveprojection = Makie.Orthographic)
+
+        for set_cam! in (cam2d!, cam_relative!, campixel!, cam3d!, orthographic_cam3d!)
+            @testset "$set_cam!" begin
+                set_cam!(scene)
+                ray = Makie.Ray(scene, xy)
+                ref_ray = Makie.ray_from_projectionview(scene, xy)
+                # Direction matches and is normalized
+                @test ref_ray.direction ≈ ray.direction
+                @test norm(ray.direction) ≈ 1f0
+                # origins are on the same ray
+                @test Makie.is_point_on_ray(ray.origin, ref_ray)
+            end
+        end
+    end
+
+
+    # transform() is used to apply a translation-rotation-scale matrix to rays 
+    # instead of point like data
+    # Generate random point + transform
+    rot = Makie.rotation_between(rand(Vec3f), rand(Vec3f))
+    model = Makie.transformationmatrix(rand(Vec3f), rand(Vec3f), rot)
+    point = Point3f(1) + rand(Point3f)
+    
+    # Generate rate that passes through transformed point
+    transformed = Point3f(model * Point4f(point..., 1))
+    direction = (1 + 10*rand()) * rand(Vec3f)
+    ray = Makie.Ray(transformed + direction, normalize(direction))
+
+    @test Makie.is_point_on_ray(transformed, ray)
+    transformed_ray = Makie.transform(inv(model), ray)
+    @test Makie.is_point_on_ray(point, transformed_ray)
+
+
+    @testset "Intersections" begin
+        p = rand(Point3f)
+        v = rand(Vec3f)
+        ray = Makie.Ray(p + 10v, normalize(v))
+
+        # ray - line
+        w = cross(v, rand(Vec3f))
+        A = p - 5w
+        B = p + 5w
+        result = Makie.closest_point_on_line(A, B, ray)
+        @test result ≈ p
+
+        # ray - triangle
+        w2 = cross(v, w)
+        A = p - 5w - 5w2
+        B = p + 5w
+        C = p + 5w2
+        result = Makie.ray_triangle_intersection(A, B, C, ray)
+        @test result ≈ p
+
+        # ray - rect3
+        rect = Rect(Vec(A), 10w + 10w2 + 10v)
+        result = Makie.ray_rect_intersection(rect, ray)
+        @test Makie.is_point_on_ray(result, ray)
+
+        # ray - rect2
+        p2 = Point2f(ray.origin - ray.origin[3] / ray.direction[3] * ray.direction)
+        w = rand(Vec2f)
+        rect = Rect2f(p2 - 5w, 10w)
+        result = Makie.ray_rect_intersection(rect, ray)
+        @test result ≈ Point3f(p2..., 0)
+    end
+
+
+    # Note that these tests depend on the exact placement of plots and may 
+    # error when cameras are adjusted
+    @testset "position_on_plot()" begin
+    
+        # Lines (2D) & Linesegments (3D)
+        ps = [exp(-0.01phi) * Point2f(cos(phi), sin(phi)) for phi in range(0, 20pi, length = 501)]
+        scene = Scene(resolution = (400, 400))
+        p = lines!(scene, ps)
+        cam2d!(scene)
+        ray = Makie.Ray(scene, (325.0, 313.0))
+        pos = Makie.position_on_plot(p, 157, ray)
+        @test pos ≈ Point3f(0.6087957666683925, 0.5513198993583837, 0.0)
+    
+        scene = Scene(resolution = (400, 400))
+        p = linesegments!(scene, ps)
+        cam3d!(scene)
+        ray = Makie.Ray(scene, (238.0, 233.0))
+        pos = Makie.position_on_plot(p, 178, ray)
+        @test pos ≈ Point3f(-0.7850463447725504, -0.15125213957100314, 0.0)
+    
+    
+        # Heatmap (2D) & Image (3D)
+        scene = Scene(resolution = (400, 400))
+        p = heatmap!(scene, 0..1, -1..1, rand(10, 10))
+        cam2d!(scene)
+        ray = Makie.Ray(scene, (228.0, 91.0))
+        pos = Makie.position_on_plot(p, 0, ray)
+        @test pos ≈ Point3f(0.13999999, -0.54499996, 0.0)
+    
+        scene = Scene(resolution = (400, 400))
+        p = image!(scene, -1..1, -1..1, rand(10, 10))
+        cam3d!(scene)
+        ray = Makie.Ray(scene, (309.0, 197.0))
+        pos = Makie.position_on_plot(p, 3, ray)
+        @test pos ≈ Point3f(-0.7830243, 0.8614166, 0.0)
+    
+    
+        # Mesh (3D)
+        scene = Scene(resolution = (400, 400))
+        p = mesh!(scene, Rect3f(Point3f(0), Vec3f(1)))
+        cam3d!(scene)
+        ray = Makie.Ray(scene, (201.0, 283.0))
+        pos = Makie.position_on_plot(p, 15, ray)
+        @test pos ≈ Point3f(0.029754717, 0.043159597, 1.0)
+    
+        # Surface (3D)
+        scene = Scene(resolution = (400, 400))
+        p = surface!(scene, -2..2, -2..2, [sin(x) * cos(y) for x in -10:10, y in -10:10])
+        cam3d!(scene)
+        ray = Makie.Ray(scene, (52.0, 238.0))
+        pos = Makie.position_on_plot(p, 57, ray)
+        @test pos ≈ Point3f(0.80910987, -1.6090667, 0.137722)
+    
+        # Volume (3D)
+        scene = Scene(resolution = (400, 400))
+        p = volume!(scene, rand(10, 10, 10))
+        cam3d!(scene)
+        center!(scene)
+        ray = Makie.Ray(scene, (16.0, 306.0))
+        pos = Makie.position_on_plot(p, 0, ray)
+        @test pos ≈ Point3f(10.0, 0.18444633, 9.989262)
+    end
+
+    # For recreating the above:
+    #= 
+    # Scene setup from tests:
+    scene = Scene(resolution = (400, 400))
+    p = surface!(scene, -2..2, -2..2, [sin(x) * cos(y) for x in -10:10, y in -10:10])
+    cam3d!(scene)
+    
+    pos = Observable(Point3f(0.5))
+    on(events(scene).mousebutton, priority = 100) do event
+        if event.button == Mouse.left && event.action == Mouse.press
+            mp = events(scene).mouseposition[]
+            _p, idx = pick(scene, mp, 10)
+            pos[] = Makie.position_on_plot(p, idx)
+            println(_p == p)
+            println("ray = Makie.Ray(scene, $mp)")
+            println("pos = Makie.position_on_plot(p, $idx, ray)")
+            println("@test pos ≈ Point3f(", pos[][1], ", ", pos[][2], ", ", pos[][3], ")")
+        end
+    end
+
+    # Optional - show selected positon
+    # This may change the camera, so don't use it for test values
+    # scatter!(scene, pos)
+
+    scene
+    =#
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 5d46b19a92e..1e2e257d552 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -32,4 +32,5 @@ using Makie: volume
     include("events.jl")
     include("text.jl")
     include("boundingboxes.jl")
+    include("ray_casting.jl")
 end