diff --git a/src/ImageSegmentation.jl b/src/ImageSegmentation.jl
index 67115ff..b398205 100644
--- a/src/ImageSegmentation.jl
+++ b/src/ImageSegmentation.jl
@@ -4,6 +4,7 @@ import Base: show
 
 using LinearAlgebra, Statistics
 using DataStructures, StaticArrays, ImageCore, ImageFiltering, ImageMorphology, LightGraphs, SimpleWeightedGraphs, RegionTrees, Distances, StaticArrays, Clustering, MetaGraphs
+using ImageCore.MappedArrays: of_eltype
 import Clustering: kmeans, fuzzy_cmeans
 
 include("compat.jl")
@@ -41,7 +42,7 @@ export
     kmeans,
     fuzzy_cmeans,
     merge_segments,
-    
+
     # types
     SegmentedImage,
     ImageEdge
diff --git a/src/clustering.jl b/src/clustering.jl
index 7133f7d..1ad5532 100644
--- a/src/clustering.jl
+++ b/src/clustering.jl
@@ -1,13 +1,11 @@
 function img_to_data(img::AbstractArray{T,N}) where T<:Colorant where N
-    AT = accum_type(T)
-    aimg = AT.(img)
-    pimg = parent(aimg)
-    ch = channelview(pimg)
-    data = reshape(ch, :, *(size(pimg)...))
+    aimg = of_eltype(accum_type(T), img)
+    ch = channelview(aimg)
+    return reshape(ch, :, *(size(img)...))
 end
 
 kmeans(img::AbstractArray{T,N}, args...; kwargs...) where {T<:Colorant,N} =
-    kmeans(img_to_data(img), args...; kwargs...)
+    SegmentedImage(kmeans(img_to_data(img), args...; kwargs...), img)
 
 fuzzy_cmeans(img::AbstractArray{T,N}, args...; kwargs...) where {T<:Colorant,N} =
     fuzzy_cmeans(img_to_data(img), args...; kwargs...)
diff --git a/src/core.jl b/src/core.jl
index 166f1fa..2932638 100644
--- a/src/core.jl
+++ b/src/core.jl
@@ -14,6 +14,13 @@ struct SegmentedImage{T<:AbstractArray,U<:Union{Colorant,Real}}
     segment_pixel_count::Dict{Int,Int}
 end
 
+SegmentedImage(r::Clustering.KmeansResult, img) =
+    SegmentedImage(reshape(r.assignments, axes(img)),
+                   collect(1:length(r.counts)),
+                   Dict(zip(1:length(r.counts), colorview(ImageCore.ColorTypes.base_colorant_type(eltype(img)), r.centers))),
+                   Dict(zip(1:length(r.counts), r.counts))
+    )
+
 """
     edge = ImageEdge(index1, index2, weight)
 
diff --git a/test/clustering.jl b/test/clustering.jl
new file mode 100644
index 0000000..857a6e9
--- /dev/null
+++ b/test/clustering.jl
@@ -0,0 +1,24 @@
+@testset "Clustering" begin
+    # Test the example in the docs
+    path = download("https://github.com/JuliaImages/juliaimages.github.io/raw/source/docs/src/pkgs/segmentation/assets/flower.jpg")
+    img = load(path)
+    r = fuzzy_cmeans(img, 3, 2)
+    @test size(r.centers) == (3,3)
+    @test size(r.weights, 1) == length(img)
+    @test all(≈(1), sum(r.weights; dims=2))
+    cmin, cmax = extrema(sum(r.weights; dims=1))
+    @test cmax < 3*cmin
+    # Make sure we support OffsetArrays
+    imgo = OffsetArray(img, (-1, -1))
+    ro = fuzzy_cmeans(imgo, 3, 2)
+    @test size(ro.centers) == (3,3)
+
+    # Also with kmeans
+    r = kmeans(img, 3)
+    nc = last.(sort([pr for pr in segment_pixel_count(r)]; by=first))
+    @test sum(nc) == length(img)
+    @test 3*minimum(nc) > maximum(nc)
+    ro = kmeans(imgo, 3)
+    nc = last.(sort([pr for pr in segment_pixel_count(ro)]; by=first))
+    @test sum(nc) == length(img)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index a8c257c..59e11b8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 using ImageSegmentation, Images, Test, SimpleWeightedGraphs, LightGraphs, StaticArrays, DataStructures
-using RegionTrees: isleaf, Cell, split! 
+using Images.OffsetArrays
+using RegionTrees: isleaf, Cell, split!
 using MetaGraphs: MetaGraph, clear_props!, get_prop, has_prop, set_prop!, props, vertices
 
 using Documenter
@@ -12,4 +13,5 @@ include("fast_scanning.jl")
 include("watershed.jl")
 include("region_merging.jl")
 include("meanshift.jl")
+include("clustering.jl")
 include("merge_segments.jl")