Initial functionality (#1)

Introduces `quantize(img, alg)` and the following algorithms:
* `UniformQuantization` 
* `KMeansQuantization` via Clustering.jl

Author: adrhill

.github/workflows/CI.yml

@@ -19,6 +19,7 @@ jobs:
       matrix:
         version:
           - '1.6'
+          - '1'
           - 'nightly'
         os:
           - ubuntu-latest

.gitignore

@@ -1,5 +1,5 @@
 *.jl.*.cov
 *.jl.cov
 *.jl.mem
-/Manifest.toml
+Manifest.toml
 /docs/build/

Project.toml

@@ -3,5 +3,15 @@ uuid = "652893fb-f6a0-4a00-a44a-7fb8fac69e01"
 authors = ["Adrian Hill <[email protected]>"]
 version = "0.1.0"
 
+[deps]
+Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+ImageBase = "c817782e-172a-44cc-b673-b171935fbb9e"
+LazyModules = "8cdb02fc-e678-4876-92c5-9defec4f444e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
 [compat]
+Clustering = "0.14.3"
+Colors = "0.12"
+LazyModules = "0.3"
 julia = "1.6"

src/ColorQuantization.jl

@@ -1,5 +1,22 @@
 module ColorQuantization
 
-# Write your package code here.
+using Colors
+using ImageBase: FixedPoint, floattype, FixedPointNumbers.rawtype
+using ImageBase: channelview, colorview, restrict
+using Random: AbstractRNG, GLOBAL_RNG
+using LazyModules: @lazy
+#! format: off
+@lazy import Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
+#! format: on
+
+abstract type AbstractColorQuantizer end
+
+include("api.jl")
+include("utils.jl")
+include("uniform.jl")
+include("clustering.jl") # lazily loaded
+
+export AbstractColorQuantizer, quantize
+export UniformQuantization, KMeansQuantization
 
 end

src/api.jl

@@ -0,0 +1,16 @@
+"""
+    quantize([T,] cs, alg)
+
+Apply color quantization algorithm `alg` to an iterable collection of Colorants,
+e.g. an image or any `AbstractArray`.
+The return type `T` can be specified and defaults to the element type of `cs`.
+"""
+function quantize(cs::AbstractArray{T}, alg::AbstractColorQuantizer) where {T<:Colorant}
+    return quantize(T, cs, alg)
+end
+
+function quantize(
+    ::Type{T}, cs::AbstractArray{<:Colorant}, alg::AbstractColorQuantizer
+) where {T}
+    return convert.(T, alg(cs)[:])
+end

src/clustering.jl

@@ -0,0 +1,69 @@
+# The following code is lazily loaded from Clustering.jl using LazyModules.jl
+# Code adapted from @cormullion's ColorSchemeTools (https://github.com/JuliaGraphics/ColorSchemeTools.jl)
+
+# The following type definition is taken from from Clustering.jl for the kwarg `init`:
+const ClusteringInitType = Union{
+    Symbol,Clustering.SeedingAlgorithm,AbstractVector{<:Integer}
+}
+
+const KMEANS_DEFAULT_COLORSPACE = RGB{Float32}
+
+"""
+    KMeansQuantization([T=RGB,] ncolors)
+
+Quantize colors by applying the K-means method, where `ncolors` corresponds to
+the amount of clusters and output colors.
+
+The colorspace `T` in which K-means are computed defaults to `RGB`.
+
+## Optional arguments
+The following keyword arguments from Clustering.jl can be specified:
+- `init`: specifies how cluster seeds are initialized
+- `maxiter`: maximum number of iterations
+- `tol`:  minimal allowed change of the objective during convergence.
+    The algorithm is considered to be converged when the change of objective value between
+    consecutive iterations drops below `tol`.
+
+The default values are carried over from are imported from Clustering.jl.
+For more details, refer to the [documentation](https://juliastats.org/Clustering.jl/stable/)
+of Clustering.jl.
+"""
+struct KMeansQuantization{T<:Colorant,I<:ClusteringInitType,R<:AbstractRNG} <:
+       AbstractColorQuantizer
+    ncolors::Int
+    maxiter::Int
+    tol::Float64
+    init::I
+    rng::R
+
+    function KMeansQuantization(
+        T::Type{<:Colorant},
+        ncolors::Integer;
+        init=Clustering._kmeans_default_init,
+        maxiter=Clustering._kmeans_default_maxiter,
+        tol=Clustering._kmeans_default_tol,
+        rng=GLOBAL_RNG,
+    )
+        ncolors ≥ 2 ||
+            throw(ArgumentError("K-means clustering requires ncolors ≥ 2, got $(ncolors)."))
+        return new{T,typeof(init),typeof(rng)}(ncolors, maxiter, tol, init, rng)
+    end
+end
+function KMeansQuantization(ncolors::Integer; kwargs...)
+    return KMeansQuantization(KMEANS_DEFAULT_COLORSPACE, ncolors; kwargs...)
+end
+
+function (alg::KMeansQuantization{T})(cs::AbstractArray{<:Colorant}) where {T}
+    # Clustering on the downsampled image already generates good enough colormap estimation.
+    # This significantly reduces the algorithmic complexity.
+    cs = _restrict_to(cs, alg.ncolors * 100)
+    return _kmeans(alg, convert.(T, cs))
+end
+
+function _kmeans(alg::KMeansQuantization, cs::AbstractArray{<:Colorant{T,N}}) where {T,N}
+    data = reshape(channelview(cs), N, :)
+    R = Clustering.kmeans(
+        data, alg.ncolors; maxiter=alg.maxiter, tol=alg.tol, init=alg.init, rng=alg.rng
+    )
+    return colorview(eltype(cs), R.centers)
+end

src/uniform.jl

@@ -0,0 +1,51 @@
+"""
+    UniformQuantization(n::Int)
+
+Quantize colors in RGB color space by dividing each dimension of the ``[0, 1]³``
+RGB color cube into `n` equidistant steps for a total of `n³` cubes of equal size.
+Each color in `cs` is then quantized to the center of the cube it is in.
+Only unique colors are returned. The amount of output colors is therefore bounded by `n³`.
+"""
+struct UniformQuantization{N} <: AbstractColorQuantizer
+    function UniformQuantization(n::Integer)
+        return n < 1 ? throw(ArgumentError("n has to be ≥ 1, got $(n).")) : new{n}()
+    end
+end
+
+# UniformQuantization is applied in RGB colorspace:
+const UQ_RGB = RGB{Float32}
+@inline (alg::UniformQuantization)(cs::AbstractArray) = alg(convert.(UQ_RGB, cs))
+@inline (alg::UniformQuantization)(c::Colorant) = alg(convert(UQ_RGB, c))
+
+@inline (alg::UniformQuantization)(cs::AbstractArray{<:RGB}) = unique(alg.(cs))
+@inline (alg::UniformQuantization)(c::RGB) = mapc(alg, c)
+
+@inline (alg::UniformQuantization)(x::Real) = _uq_round(alg, x)
+@inline (alg::UniformQuantization)(x::FixedPoint) = _uq_lookup(alg, x)
+
+# For general real numbers, use round
+@inline @generated function _uq_round(::UniformQuantization{N}, x::T) where {N,T<:Real}
+    return quote
+        x < $(T(1 / N)) && return $(T(1 / (2 * N)))
+        x ≥ $(T((N - 1) / N)) && return $(T((2 * N - 1) / (2 * N)))
+        return (round(x * $N - $(T(0.5))) + $(T(0.5))) / $N
+    end
+end
+
+# For fixed-point numbers, use the internal integer to index onto a lookup table
+@inline function _uq_lookup(alg::UniformQuantization{N}, x::T) where {N,T<:FixedPoint}
+    @inbounds _uq_lookup_table(alg, T)[x.i + 1]
+end
+# the lookup table is generated at compile time
+@generated function _uq_lookup_table(
+    ::UniformQuantization{N}, ::Type{T}
+) where {N,T<:FixedPoint}
+    RT = rawtype(T)
+    tmax = typemax(RT)
+    table = Vector{T}(undef, tmax + 1)
+    for raw_x in zero(RT):tmax
+        x = reinterpret(T, raw_x)
+        table[raw_x + 1] = _uq_round(UniformQuantization(N), x)
+    end
+    return table
+end

src/utils.jl

@@ -0,0 +1,9 @@
+# Reduce the size of img until there are at most n elements left
+function _restrict_to(img, n)
+    length(img) <= n && return img
+    out = restrict(img)
+    while length(out) > n
+        out = restrict(out)
+    end
+    return out
+end

test/Manifest.toml

@@ -1,2 +1,13 @@
 [deps]
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
+
+[compat]
+Colors = "0.12"
+ReferenceTests = "0.10"
+StableRNGs = "1"
+TestImages = "1"