diff --git a/CHANGELOG.md b/CHANGELOG.md
index 334df3a5727..45f82637f21 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 - Show DataInspector tooltip on NaN values if `nan_color` has been set to other than `:transparent` [#4310](https://github.com/MakieOrg/Makie.jl/pull/4310)
 - Fix `linestyle` not being used in `triplot` [#4332](https://github.com/MakieOrg/Makie.jl/pull/4332)
+- Fix voxel clipping not being based on voxel centers [#4397](https://github.com/MakieOrg/Makie.jl/pull/4397)
 
 ## [0.21.11] - 2024-09-13
 
@@ -22,6 +23,7 @@
 - `plotfunc()` and `func2type()` support functions ending with `!` [#4275](https://github.com/MakieOrg/Makie.jl/pull/4275).
 - Fixed Boundserror in clipped multicolor lines in CairoMakie [#4313](https://github.com/MakieOrg/Makie.jl/pull/4313)
 - Fix float precision based assertions error in GLMakie.volume [#4311](https://github.com/MakieOrg/Makie.jl/pull/4311)
+ - Support images with reversed axes [#4338](https://github.com/MakieOrg/Makie.jl/pull/4338)
 
 ## [0.21.9] - 2024-08-27
 
diff --git a/CairoMakie/src/primitives.jl b/CairoMakie/src/primitives.jl
index 6478d404ac1..c3a860f5a8d 100644
--- a/CairoMakie/src/primitives.jl
+++ b/CairoMakie/src/primitives.jl
@@ -714,14 +714,14 @@ function draw_atomic(scene::Scene, screen::Screen{RT}, @nospecialize(primitive::
     image = primitive[3][]
     xs, ys = primitive[1][], primitive[2][]
     if xs isa Makie.EndPoints
-        l, r = extrema(xs)
+        l, r = xs
         N = size(image, 1)
         xs = range(l, r, length = N+1)
     else
         xs = regularly_spaced_array_to_range(xs)
     end
     if ys isa Makie.EndPoints
-        l, r = extrema(ys)
+        l, r = ys
         N = size(image, 2)
         ys = range(l, r, length = N+1)
     else
@@ -1237,16 +1237,15 @@ function draw_atomic(scene::Scene, screen::Screen, @nospecialize(primitive::Maki
     pos = Makie.voxel_positions(primitive)
     scale = Makie.voxel_size(primitive)
     colors = Makie.voxel_colors(primitive)
+    marker = GeometryBasics.clear_faceviews(normal_mesh(Rect3f(Point3f(-0.5), Vec3f(1))))
+
 
-    # Face culling
     if !isempty(primitive.clip_planes[]) && Makie.is_data_space(primitive.space[])
         valid = [is_visible(primitive.clip_planes[], p) for p in pos]
         pos = pos[valid]
         colors = colors[valid]
     end
 
-    marker = GeometryBasics.clear_faceviews(normal_mesh(Rect3f(Point3f(-0.5), Vec3f(1))))
-
     # For correct z-ordering we need to be in view/camera or screen space
     model = copy(primitive.model[])::Mat4d
     view = scene.camera.view[]
diff --git a/GLMakie/assets/shader/voxel.frag b/GLMakie/assets/shader/voxel.frag
index 71f11e2c356..c9c6873210e 100644
--- a/GLMakie/assets/shader/voxel.frag
+++ b/GLMakie/assets/shader/voxel.frag
@@ -85,18 +85,16 @@ vec4 get_color(sampler2D color, Nothing color_map, int id) {
 
 bool is_clipped()
 {
-    float d1, d2;
+    float d;
+    // Center of voxel
     ivec3 size = ivec3(textureSize(voxel_id, 0).xyz);
-    vec3 xyz = vec3(ivec3(o_uvw * size));
+    vec3 xyz = vec3(ivec3(o_uvw * size)) + vec3(0.5);
     for (int i = 0; i < _num_clip_planes; i++) {
-        // distance from clip planes with negative clipped
-        d1 = dot(xyz, clip_planes[i].xyz) - clip_planes[i].w;
-        d2 = dot(xyz, clip_planes[i].xyz) - clip_planes[i].w;
+        // distance between clip plane and center
+        d = dot(xyz, clip_planes[i].xyz) - clip_planes[i].w;
 
-        // both outside - clip everything
-        if (d1 < 0.0 || d2 < 0.0) {
+        if (d < 0.0) 
             return true;
-        }
     }
 
     return false;
diff --git a/GLMakie/src/drawing_primitives.jl b/GLMakie/src/drawing_primitives.jl
index 98983b6c64e..12c05e0712b 100644
--- a/GLMakie/src/drawing_primitives.jl
+++ b/GLMakie/src/drawing_primitives.jl
@@ -691,8 +691,8 @@ end
 function draw_image(screen::Screen, scene::Scene, plot::Union{Heatmap, Image})
     return cached_robj!(screen, scene, plot) do gl_attributes
         position = lift(plot, plot[1], plot[2]) do x, y
-            xmin, xmax = extrema(x)
-            ymin, ymax = extrema(y)
+            xmin, xmax = x
+            ymin, ymax = y
             rect = Rect2(xmin, ymin, xmax - xmin, ymax - ymin)
             return decompose(Point2d, rect)
         end
diff --git a/Project.toml b/Project.toml
index d34ed86ea0b..184f59e60d3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -85,9 +85,9 @@ FreeTypeAbstraction = "0.10.3"
 GeometryBasics = "0.4.11"
 GridLayoutBase = "0.11"
 ImageBase = "0.1.7"
-ImageIO = "0.2, 0.3, 0.4, 0.5, 0.6"
+ImageIO = "0.5, 0.6"
 InteractiveUtils = "1.0, 1.6"
-Interpolations = "0.15.1"
+Interpolations = "0.14, 0.15.1"
 IntervalSets = "0.3, 0.4, 0.5, 0.6, 0.7"
 Isoband = "0.1"
 KernelDensity = "0.5, 0.6"
diff --git a/ReferenceTests/src/tests/examples3d.jl b/ReferenceTests/src/tests/examples3d.jl
index 27879658561..364541c679b 100644
--- a/ReferenceTests/src/tests/examples3d.jl
+++ b/ReferenceTests/src/tests/examples3d.jl
@@ -682,7 +682,7 @@ end
 @reference_test "Clip planes - voxel" begin
     f = Figure()
     a = LScene(f[1, 1])
-    a.scene.theme[:clip_planes][] = [Plane3f(Vec3f(-2, -1, -0.5), 0.0), Plane3f(Vec3f(-0.5, -1, -2), 0.0)]
+    a.scene.theme[:clip_planes][] = [Plane3f(Vec3f(-2, -1, -0.5), 0.1), Plane3f(Vec3f(-0.5, -1, -2), 0.1)]
     r = -10:10
     p = voxels!(a, [cos(sin(x+y)+z) for x in r, y in r, z in r])
     f
diff --git a/ReferenceTests/src/tests/primitives.jl b/ReferenceTests/src/tests/primitives.jl
index 862b03cd51a..8f8465ba58b 100644
--- a/ReferenceTests/src/tests/primitives.jl
+++ b/ReferenceTests/src/tests/primitives.jl
@@ -818,3 +818,31 @@ end
     end
     f
 end
+
+@reference_test "Reverse image, heatmap and surface axes" begin
+    img = [2 0 0 3; 0 0 0 0; 1 1 0 0; 1 1 0 4]
+
+    f = Figure(size = (600, 400))
+
+    for (i, interp) in enumerate((true, false))
+        for (j, plot_func) in enumerate((
+            (fp, x, y, cs, interp) -> image(fp, x, y, cs, colormap = :viridis, interpolate = interp), 
+            (fp, x, y, cs, interp) -> heatmap(fp, x, y, cs, colormap = :viridis, interpolate = interp), 
+            (fp, x, y, cs, interp) -> surface(fp, x, y, zeros(size(cs)), color = cs, colormap = :viridis, interpolate = interp, shading = NoShading)
+        ))
+
+            gl = GridLayout(f[i, j])
+
+            a, p = plot_func(gl[1, 1], 1:4, 1:4, img, interp)
+            hidedecorations!(a)
+            a, p = plot_func(gl[2, 1], 1:4, 4..1, img, interp)
+            hidedecorations!(a)
+            a, p = plot_func(gl[1, 2], 4:-1:1, 1:4, img, interp)
+            hidedecorations!(a)
+            a, p = plot_func(gl[2, 2], 4:-1:1, [4, 3, 2, 1], img, interp)
+            hidedecorations!(a)
+        end
+    end
+
+    f
+end
\ No newline at end of file
diff --git a/WGLMakie/assets/voxel.frag b/WGLMakie/assets/voxel.frag
index 3fe7a6dfeee..4e095306341 100644
--- a/WGLMakie/assets/voxel.frag
+++ b/WGLMakie/assets/voxel.frag
@@ -85,18 +85,15 @@ vec3 blinnphong(vec3 N, vec3 V, vec3 L, vec3 color){
 
 bool is_clipped()
 {
-    float d1, d2;
+    float d;
+    // get center pos of this voxel
     vec3 size = vec3(textureSize(voxel_id, 0).xyz);
-    vec3 xyz = vec3(ivec3(o_uvw * size));
+    vec3 xyz = vec3(ivec3(o_uvw * size)) + vec3(0.5);
     for (int i = 0; i < num_clip_planes; i++) {
-        // distance from clip planes with negative clipped
-        d1 = dot(xyz, clip_planes[i].xyz) - clip_planes[i].w;
-        d2 = dot(xyz, clip_planes[i].xyz) - clip_planes[i].w;
-
-        // both outside - clip everything
-        if (d1 < 0.0 || d2 < 0.0) {
+        // distance between clip plane and voxel center
+        d = dot(xyz, clip_planes[i].xyz) - clip_planes[i].w;
+        if (d < 0.0)
             return true;
-        }
     }
 
     return false;
diff --git a/WGLMakie/src/imagelike.jl b/WGLMakie/src/imagelike.jl
index 28952f4fda5..7734d8cf2b4 100644
--- a/WGLMakie/src/imagelike.jl
+++ b/WGLMakie/src/imagelike.jl
@@ -129,7 +129,7 @@ end
 
 
 
-xy_convert(x::Makie.EndPoints, n) = LinRange(extrema(x)..., n + 1)
+xy_convert(x::Makie.EndPoints, n) = LinRange(x..., n + 1)
 xy_convert(x::AbstractArray, n) = x
 
 # TODO, speed up GeometryBasics
@@ -166,8 +166,8 @@ function limits_to_uvmesh(plot, f32c)
     py = lift(identity, plot, py; ignore_equal_values=true)
     if px[] isa Makie.EndPoints && py[] isa Makie.EndPoints && Makie.is_identity_transform(t)
         rect = lift(plot, px, py) do x, y
-            xmin, xmax = extrema(x)
-            ymin, ymax = extrema(y)
+            xmin, xmax = x
+            ymin, ymax = y
             return Rect2f(xmin, ymin, xmax - xmin, ymax - ymin)
         end
         ps = lift(rect -> decompose(Point2f, rect), plot, rect)
diff --git a/src/conversions.jl b/src/conversions.jl
index 5cd5a80a6e2..5815e061566 100644
--- a/src/conversions.jl
+++ b/src/conversions.jl
@@ -364,12 +364,10 @@ function to_endpoints(x::Tuple{<:Real,<:Real})
     T = float_type(x...)
     return EndPoints(T.(x))
 end
-to_endpoints(x::ClosedInterval) = to_endpoints(endpoints(x))
-function to_endpoints(x::Union{Interval,AbstractVector,ClosedInterval})
-    return to_endpoints((minimum(x), maximum(x)))
-end
+to_endpoints(x::Interval) = to_endpoints(endpoints(x))
+to_endpoints(x::EndPoints) = x
+to_endpoints(x::AbstractVector) = to_endpoints((first(x), last(x)))
 function to_endpoints(x, dim)
-    # having minimum and maximum here actually invites bugs
     x isa AbstractVector && !(x isa EndPoints) && print_range_warning(dim, x)
     return to_endpoints(x)
 end
@@ -692,7 +690,8 @@ end
 #                               Helper Functions                               #
 ################################################################################
 
-to_linspace(interval, N) = range(minimum(interval), stop = maximum(interval), length = N)
+to_linspace(interval::Interval, N) = range(leftendpoint(interval), stop = rightendpoint(interval), length = N)
+to_linspace(x, N) = range(first(x), stop = last(x), length = N)
 
 """
 Converts the element array type to `T1` without making a copy if the element type matches
diff --git a/src/interaction/inspector.jl b/src/interaction/inspector.jl
index d2d548347bb..64a1ee9b45c 100644
--- a/src/interaction/inspector.jl
+++ b/src/interaction/inspector.jl
@@ -677,10 +677,12 @@ function show_imagelike(inspector, plot, name, edge_based)
         tt.text[] = plot[:inspector_label][](plot, (i, j), ins_p)
     end
 
-    a._color[] = if z isa Real
-        get(plot.calculated_colors[], z)
+    if z isa Real
+        if haskey(plot, :calculated_colors)
+            a._color[] = get(plot.calculated_colors[], z)
+        end
     else
-        z
+        a._color[] = z
     end
 
     proj_pos = Point2f(mouseposition_px(inspector.root))
diff --git a/test/boundingboxes.jl b/test/boundingboxes.jl
index 66897c8444f..88cda42f179 100644
--- a/test/boundingboxes.jl
+++ b/test/boundingboxes.jl
@@ -92,6 +92,11 @@ end
     @test bb.origin ≈ Point3f(0)
     @test bb.widths ≈ Vec3f(10.0, 10.0, 0)
 
+    fig, ax, p = image(1..0, 1:10, rand(10, 10))
+    bb = boundingbox(p)
+    @test bb.origin ≈ Point3f(0, 1, 0)
+    @test bb.widths ≈ Vec3f(1.0, 9.0, 0)
+
     # text transforms to pixel space atm (TODO)
     fig = Figure(size = (400, 400))
     ax = Axis(fig[1, 1])
diff --git a/test/conversions.jl b/test/conversions.jl
index 64959083868..2ecd39f57dc 100644
--- a/test/conversions.jl
+++ b/test/conversions.jl
@@ -342,9 +342,11 @@ end
 
     v1 = collect(1:10)
     v2 = collect(1:6)
+    v3 = reverse(v1)
 
     i1 = 1 .. 10
     i2 = 1 .. 6
+    i3 = 10 .. 1
 
     o3 = Float32.(m3)
 
@@ -354,6 +356,8 @@ end
         @test convert_arguments(Image, m3) == ((0.0f0, 10.0f0), (0.0f0, 6.0f0), o3)
         @test convert_arguments(Image, v1, r2, m3) == ((1.0f0, 10.0f0), (1.0f0, 6.0f0), o3)
         @test convert_arguments(Image, i1, v2, m3) == ((1.0f0, 10.0f0), (1.0f0, 6.0f0), o3)
+        @test convert_arguments(Image, v3, i1, m3) == ((10, 1), (1, 10), o3)
+        @test convert_arguments(Image, v1, i3, m3) == ((1, 10), (10, 1), o3)
         @test convert_arguments(Image, m1, m2, m3) === (m1, m2, m3)
         @test convert_arguments(Heatmap, m1, m2) === (m1, m2)
     end