introduce sumfinite, enhancement to *finite, rename minfinite/maxfinite (#17)

Features:

* (new function) add `sumfinite([f=identity], A)`
* add `meanfinite([f=identity], A)`
* add `minimum_finite([f=identity], A)`
* add `maximum_finite([f=identity], a)`

Bug fixes:

* "fixes" `varfinite` for `RGB` array inputs using dot product. This requires ColorVectorSpace v0.9.7
* support `minimum_finite(A; dims)` and `maximum_finite(A; dims)`

Deprecations:

* rename: `minfinite` -> `minimum_finite`
* rename: `maxfinite` -> `maximum_finite`
* deprecate `maxabsfinite(A)` in favor of `maximum_finite(abs, A)`

There are some type infer issues, which I leave as future work.

Author: johnnychen94

Project.toml

@@ -7,7 +7,7 @@ ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
-ImageCore = "0.8, 0.9"
+ImageCore = "0.9.3"
 Reexport = "0.2, 1"
 julia = "1"
 

src/ImageBase.jl

@@ -11,12 +11,15 @@ export
     fdiff!,
 
     # basic image statistics, from Images.jl
-    minfinite,
-    maxfinite,
-    maxabsfinite,
+    minimum_finite,
+    maximum_finite,
     meanfinite,
-    varfinite
+    varfinite,
+    sumfinite
 
+# Introduced in ColorVectorSpace v0.9.3
+# https://github.com/JuliaGraphics/ColorVectorSpace.jl/pull/172
+using ImageCore.ColorVectorSpace.Future: abs2
 
 using Reexport
 

src/deprecated.jl

@@ -2,6 +2,10 @@
 
 @deprecate restrict(A::AbstractArray, region::Vector{Int}) restrict(A, (region...,))
 
-@deprecate meanfinite(A, region) meanfinite(A; dims=region)
+@deprecate meanfinite(A::AbstractArray, region) meanfinite(A; dims=region)
+
+@deprecate minfinite(A; kwargs...) minimum_finite(A; kwargs...)
+@deprecate maxfinite(A; kwargs...) maximum_finite(A; kwargs...)
+@deprecate maxabsfinite(A; kwargs...) maximum_finite(abs, A; kwargs...)
 
 # END 0.1 deprecation

src/statistics.jl

@@ -1,48 +1,70 @@
 """
-    minfinite(A; kwargs...)
+    minimum_finite([f=identity], A; kwargs...)
 
-Calculate the minimum value in `A`, ignoring any values that are not finite (Inf or NaN).
+Calculate `minimum(f, A)` while ignoring any values that are not finite, e.g., `Inf` or
+`NaN`.
+
+If `A` is a colorant array with multiple channels (e.g., `Array{RGB}`), the `min` comparison
+is done in channel-wise sense.
 
 The supported `kwargs` are those of `minimum(f, A; kwargs...)`.
 """
-minfinite(A; kwargs...) = mapreduce(IfElse(isfinite, identity, typemax), minc, A; kwargs...)
+function minimum_finite(f, A::AbstractArray{T}; kwargs...) where T
+    # FIXME(johnnychen94): if `typeof(f(first(A))) != eltype(A)`, this function is not type-stable.
+    mapreduce(IfElse(isfinite, f, typemax), minc, A; kwargs...)
+end
+minimum_finite(A::AbstractArray; kwargs...) = minimum_finite(identity, A; kwargs...)
 
 """
-    maxfinite(A; kwargs...)
+    maximum_finite([f=identity], A; kwargs...)
+
+Calculate `maximum(f, A)` while ignoring any values that are not finite, e.g., `Inf` or
+`NaN`.
 
-Calculate the maximum value in `A`, ignoring any values that are not finite (Inf or NaN).
+If `A` is a colorant array with multiple channels (e.g., `Array{RGB}`), the `max` comparison
+is done in channel-wise sense.
 
 The supported `kwargs` are those of `maximum(f, A; kwargs...)`.
 """
-maxfinite(A; kwargs...) = mapreduce(IfElse(isfinite, identity, typemin), maxc, A; kwargs...)
+function maximum_finite(f, A::AbstractArray{T}; kwargs...) where T
+    # FIXME(johnnychen94): if `typeof(f(first(A))) != eltype(A)`, this function is not type-stable
+    mapreduce(IfElse(isfinite, f, typemin), maxc, A; kwargs...)
+end
+maximum_finite(A::AbstractArray; kwargs...) = maximum_finite(identity, A; kwargs...)
 
 """
-    maxabsfinite(A; kwargs...)
+    sumfinite([f=identity], A; kwargs...)
 
-Calculate the maximum absolute value in `A`, ignoring any values that are not finite (Inf or NaN).
+Compute `sum(f, A)` while ignoring any non-finite values.
 
-The supported `kwargs` are those of `maximum(f, A; kwargs...)`.
+The supported `kwargs` are those of `sum(f, A; kwargs...)`.
 """
-maxabsfinite(A; kwargs...) = mapreduce(IfElse(isfinite, abs, typemin), maxc, A; kwargs...)
+sumfinite(A; kwargs...) = sumfinite(identity, A; kwargs...)
+
+if Base.VERSION >= v"1.1"
+    sumfinite(f, A; kwargs...) = sum(IfElse(isfinite, f, zero), A; kwargs...)
+else
+    sumfinite(f, A; kwargs...) = sum(IfElse(isfinite, f, zero).(A); kwargs...)
+end
 
 """
-    meanfinite(A; kwargs...)
+    meanfinite([f=identity], A; kwargs...)
 
-Compute the mean value of `A`, ignoring any non-finite values.
+Compute `mean(f, A)` while ignoring any non-finite values.
 
 The supported `kwargs` are those of `sum(f, A; kwargs...)`.
 """
-function meanfinite end
+meanfinite(A; kwargs...) = meanfinite(identity, A; kwargs...)
 
 if Base.VERSION >= v"1.1"
-    function meanfinite(A; kwargs...)
-        s = sum(IfElse(isfinite, identity, zero), A; kwargs...)
+    function meanfinite(f, A; kwargs...)
+        s = sumfinite(f, A; kwargs...)
         n = sum(IfElse(isfinite, x->true, x->false), A; kwargs...)   # TODO: replace with `Returns`
         return s./n
     end
 else
-    function meanfinite(A; kwargs...)
-        s = sum(IfElse(isfinite, identity, zero).(A); kwargs...)
+    function meanfinite(f, A; kwargs...)
+        s = sumfinite(f, A; kwargs...)
         n = sum(IfElse(isfinite, x->true, x->false).(A); kwargs...)
         return s./n
     end
@@ -54,21 +76,28 @@ end
 Compute the variance of `A`, ignoring any non-finite values.
 
 The supported `kwargs` are those of `sum(f, A; kwargs...)`.
+
+!!! note
+    This function can produce a seemingly suprising result if the input array is an RGB
+    image. To make it more clear, the implementation is made so that
+    `varfinite(img) ≈ varfinite(RGB.(img))` holds for any gray-scale image. See also
+    https://github.com/JuliaGraphics/ColorVectorSpace.jl#abs-and-abs2 for more information.
+
 """
 function varfinite end
 
 if Base.VERSION >= v"1.1"
     function varfinite(A; kwargs...)
         m = meanfinite(A; kwargs...)
         n = sum(IfElse(isfinite, x->true, x->false), A; kwargs...)   # TODO: replace with `Returns`
-        s = sum(IfElse(isfinite, identity, zero), (A .- m).^2; kwargs...)
+        s = sum(IfElse(isfinite, identity, zero), abs2.(A .- m); kwargs...)
         return s ./ max.(0, (n .- 1))
     end
 else
     function varfinite(A; kwargs...)
         m = meanfinite(A; kwargs...)
         n = sum(IfElse(isfinite, x->true, x->false).(A); kwargs...)
-        s = sum(IfElse(isfinite, identity, zero).((A .- m).^2); kwargs...)
+        s = sum(IfElse(isfinite, identity, zero).(abs2.(A .- m)); kwargs...)
         return s ./ max.(0, (n .- 1))
     end
 end

src/utils.jl

@@ -13,7 +13,14 @@ struct IfElse{C,F1,F2}
 end
 (m::IfElse)(x) = m.condition(x) ? m.f1(x) : m.f2(x)
 
+# channelwise min/max
 minc(x, y) = min(x, y)
 minc(x::Color, y::Color) = mapc(min, x, y)
 maxc(x, y) = max(x, y)
 maxc(x::Color, y::Color) = mapc(max, x, y)
+
+for (f1, f2) in ((:minc, :min), (:maxc, :max))
+    @eval function Base.reducedim_init(f, ::typeof($f1), A::AbstractArray, region)
+        Base.reducedim_init(f, $f2, A, region)
+    end
+end

test/deprecated.jl

@@ -3,4 +3,15 @@
         A = rand(N0f8, 4, 5, 3)
         @test restrict(A, [1, 2]) == restrict(A, (1, 2))
     end
+
+    @testset "statistics" begin
+        A = rand(Float32, 4, 4) .- 0.5
+        @test minfinite(A, dims=1) == minimum_finite(A, dims=1)
+        @test minfinite(A) == minimum_finite(A)
+        @test maxfinite(A, dims=1) == maximum_finite(A, dims=1)
+        @test maxfinite(A) == maximum_finite(A)
+
+        @test maxabsfinite(A) == maximum_finite(abs, A)
+        @test maxabsfinite(A, dims=1) == maximum_finite(abs, A, dims=1)
+    end
 end

test/runtests.jl

@@ -2,6 +2,7 @@ using ImageBase, OffsetArrays, StackViews
 using Test, TestImages, Aqua, Documenter
 
 using OffsetArrays: IdentityUnitRange
+include("testutils.jl")
 
 @testset "ImageBase.jl" begin
 

test/statistics.jl

@@ -1,54 +1,134 @@
 using ImageBase
+using ImageBase: varmult
 using Statistics
 using Test
 
 @testset "Reductions" begin
     _abs(x::Colorant) = mapreducec(abs, +, 0, x)
 
-    A = rand(5,5,3)
-    img = colorview(RGB, PermutedDimsArray(A, (3,1,2)))
-    s12 = sum(img, dims=(1,2))
-    @test eltype(s12) <: RGB
-
-    A = [NaN, 1, 2, 3]
-    @test meanfinite(A, dims=1) ≈ [2]
-    @test varfinite(A, dims=1) ≈ [1]
-
-    A = [NaN NaN 1;
-        1 2 3]
-    vf = varfinite(A, dims=2)
-    @test isnan(vf[1])
-
-    A = [NaN 1 2 3;
-         NaN 6 5 4]
-    mf = meanfinite(A, dims=1)
-    vf = varfinite(A, dims=1)
-    @test isnan(mf[1])
-    @test mf[2:end] ≈ [3.5,3.5,3.5]
-    @test isnan(vf[1]) 
-    @test vf[2:end] ≈ [12.5,4.5,0.5]
-
-    @test meanfinite(A, dims=2) ≈ reshape([2, 5], 2, 1)
-    @test varfinite(A, dims=2) ≈ reshape([1, 1], 2, 1)
-
-    @test meanfinite(A, dims=(1,2)) ≈ [3.5]
-    @test varfinite(A, dims=(1,2)) ≈ [3.5]
-
-    @test minfinite(A) == 1
-    @test maxfinite(A) == 6
-    @test maxabsfinite(A) == 6
-    A = rand(10:20, 5, 5)
-    @test minfinite(A) == minimum(A)
-    @test maxfinite(A) == maximum(A)
-    A = reinterpret(N0f8, rand(0x00:0xff, 5, 5))
-    @test minfinite(A) == minimum(A)
-    @test maxfinite(A) == maximum(A)
-    A = rand(Float32,3,5,5)
-    img = colorview(RGB, A)
-    dc = meanfinite(img, dims=1)-reshape(reinterpretc(RGB{Float32}, mean(A, dims=2)), (1,5))
-    @test maximum(map(_abs, dc)) < 1e-6
-    dc = minfinite(img)-RGB{Float32}(minimum(A, dims=(2,3))...)
-    @test _abs(dc) < 1e-6
-    dc = maxfinite(img)-RGB{Float32}(maximum(A, dims=(2,3))...)
-    @test _abs(dc) < 1e-6
+    @testset "sumfinite, meanfinite, varfinite" begin
+        for T in generate_test_types([N0f8, Float32], [Gray, RGB])
+            A = rand(T, 5, 5)
+            s12 = sum(A, dims=(1,2))
+            @test eltype(s12) <: Union{T, float(T), float64(T)}
+
+            @test sumfinite(A) ≈ sum(A)
+            @test sumfinite(A, dims=1) ≈ sum(A, dims=1)
+            @test sumfinite(A, dims=(1, 2)) ≈ sum(A, dims=(1, 2))
+
+            @test meanfinite(A) ≈ mean(A)
+            @test meanfinite(A, dims=1) ≈ mean(A, dims=1)
+            @test meanfinite(A, dims=(1, 2)) ≈ mean(A, dims=(1, 2))
+
+            @test varfinite(A) ≈ varmult(⋅, A)
+            @test varfinite(A, dims=1) ≈ varmult(⋅, A, dims=1)
+            @test varfinite(A, dims=(1, 2)) ≈ varmult(⋅, A, dims=(1, 2))
+
+            # test NaN/Inf
+            if eltype(T) != N0f8
+                A = rand(T, 5, 5) .- 0.5 .* oneunit(T)
+                A[1] = Inf
+                @test sum(A) ≈ A[1]
+                @test sum(abs, A) ≈ A[1]
+                @test sumfinite(A) ≈ sum(A[2:end])
+                @test sumfinite(abs, A) ≈ sum(abs, A[2:end])
+                A[1] = NaN
+                @test isnan(sum(A))
+                @test isnan(sum(abs, A))
+                @test sumfinite(A) ≈ sum(A[2:end])
+                @test sumfinite(abs, A) ≈ sum(abs, A[2:end])
+
+                A = rand(T, 5, 5) .- 0.5 .* oneunit(T)
+                A[1] = Inf
+                @test mean(A) ≈ A[1]
+                @test mean(abs, A) ≈ A[1]
+                @test meanfinite(A) ≈ mean(A[2:end])
+                @test meanfinite(abs, A) ≈ mean(abs, A[2:end])
+                A[1] = NaN
+                @test isnan(mean(A))
+                @test isnan(mean(abs, A))
+                @test meanfinite(A) ≈ mean(A[2:end])
+                @test meanfinite(abs, A) ≈ mean(abs, A[2:end])
+
+                A = rand(T, 5, 5)
+                A[1] = Inf
+                @test isnan(varmult(⋅, A))
+                @test varfinite(A) ≈ varmult(⋅, A[2:end])
+                A[1] = NaN
+                @test isnan(varmult(⋅, A))
+                @test varfinite(A) ≈ varmult(⋅, A[2:end])
+            end
+        end
+
+        A = [NaN, 1, 2, 3]
+        @test meanfinite(A, dims=1) ≈ [2]
+        @test varfinite(A, dims=1) ≈ [1]
+
+        A = [NaN NaN 1;
+            1 2 3]
+        vf = varfinite(A, dims=2)
+        @test isnan(vf[1])
+
+        A = [NaN 1 2 3;
+            NaN 6 5 4]
+        mf = meanfinite(A, dims=1)
+        vf = varfinite(A, dims=1)
+        @test isnan(mf[1])
+        @test mf[2:end] ≈ [3.5,3.5,3.5]
+        @test isnan(vf[1])
+        @test vf[2:end] ≈ [12.5,4.5,0.5]
+
+        @test meanfinite(A, dims=2) ≈ reshape([2, 5], 2, 1)
+        @test varfinite(A, dims=2) ≈ reshape([1, 1], 2, 1)
+
+        @test meanfinite(A, dims=(1,2)) ≈ [3.5]
+        @test varfinite(A, dims=(1,2)) ≈ [3.5]
+
+        # Ensure we're consistant with our decision to `abs2` in ColorVectorSpace
+        # See also: https://github.com/JuliaGraphics/ColorVectorSpace.jl/blob/master/README.md#abs-and-abs2
+        A = rand(Gray, 4, 4)
+        @test varfinite(A) ≈ varfinite(RGB.(A))
+        A[1] = Inf
+        @test varfinite(A) ≈ varfinite(RGB.(A))
+        A[1] = NaN
+        @test varfinite(A) ≈ varfinite(RGB.(A))
+    end
+
+    @testset "minfinite, maxfinite, maxabsfinite" begin
+        for T in generate_test_types([N0f8, Float32], [Gray, ])
+            A = rand(T, 5, 5)
+            @test @inferred(minimum_finite(A)) == minimum(A)
+            @test @inferred(maximum_finite(A)) == maximum(A)
+            @test minimum_finite(A; dims=1) == minimum(A; dims=1)
+            @test maximum_finite(A; dims=1) == maximum(A; dims=1)
+            @test_broken @inferred maximum_finite(A; dims=1)
+            @test_broken @inferred minimum_finite(A; dims=1)
+
+            @test maximum_finite(abs2, A) == maximum(abs2, A)
+            @test_broken @inferred maximum(abs2, A)
+            @test_broken @inferred minimum(abs2, A)
+
+            if eltype(T) != N0f8
+                A = rand(T, 5, 5) .- 0.5 * rand(T, 5, 5)
+                A[1] = Inf
+
+                @test @inferred(minimum_finite(A)) == minimum(A[2:end])
+                @test @inferred(maximum_finite(A)) == maximum(A[2:end])
+                @test minimum_finite(abs, A) == minimum(abs, A[2:end])
+                @test maximum_finite(abs2, A) == maximum(abs2, A[2:end])
+            end
+        end
+
+        # minimum_finite and maximum_finite for RGB are processed per channel
+        A = rand(RGB{Float32}, 5, 5)
+        @test minimum_finite(A) == RGB(minimum(channelview(A), dims=(2, 3))...)
+        @test maximum_finite(A) == RGB(maximum(channelview(A), dims=(2, 3))...)
+
+        # Container of abstract type
+        A = Any[1, 2, Gray(0.3)]
+        @test minimum_finite(A) == 0.3
+        @test maximum_finite(A) == 2.0
+        @test Base.return_types(minimum_finite, (typeof(A),)) == Base.return_types(minimum, (typeof(A), ))
+    end
+
 end

test/testutils.jl

@@ -0,0 +1,15 @@
+# AbstractGray and Color3 tests should be generated seperately
+function generate_test_types(number_types::AbstractArray{<:DataType}, color_types::AbstractArray{<:UnionAll})
+    test_types = map(Iterators.product(number_types, color_types)) do T
+        try
+            T[2]{T[1]}
+        catch err
+            !isa(err, TypeError) && rethrow(err)
+        end
+    end
+    test_types = filter(x->x != false, test_types)
+    if isempty(filter(x->x<:Color3, test_types))
+        test_types = [number_types..., test_types...]
+    end
+    test_types
+end