Add meanfinite (like nanmean, except ignores all non-finite values)

doc/function_reference.md

@@ -840,6 +840,11 @@ maxabsfinite(img)
 ```
 Return the minimum and maximum value in the image, respectively, ignoring any values that are not finite (Inf or NaN).
 
+```
+meanfinite(img, region)
+```
+Calculate the mean value along the dimensions listed in `region`, ignoring any non-finite values.
+
 <br />
 ```
 ssd(img1, img2)

src/Images.jl

@@ -224,6 +224,7 @@ export # types
     label_components!,
     magnitude,
     magnitude_phase,
+    nanmean,
     ncc,
     orientation,
     padarray,

src/algorithms.jl

@@ -86,6 +86,47 @@ function sum(img::AbstractImageDirect, region::Union(AbstractVector,Tuple,Intege
     out
 end
 
+meanfinite(A::AbstractArray, region) = mean(A, region)
+function meanfinite{T<:FloatingPoint}(A::AbstractArray{T}, region)
+    sz = Base.reduced_dims(A, region)
+    K = zeros(Int, sz)
+    S = zeros(T, sz)
+    sumfinite!(S, K, A)
+    S./K
+end
+# Note that you have to zero S and K upon entry
+@ngenerate N typeof((S,K)) function sumfinite!{T,N}(S, K, A::AbstractArray{T,N})
+    isempty(A) && return S, K
+    @nexprs N d->(sizeS_d = size(S,d))
+    sizeA1 = size(A, 1)
+    if size(S, 1) == 1 && sizeA1 > 1
+        # When we are reducing along dim == 1, we can accumulate to a temporary
+        @inbounds @nloops N i d->(d>1? (1:size(A,d)) : (1:1)) d->(j_d = sizeS_d==1 ? 1 : i_d) begin
+            s = @nref(N, S, j)
+            k = @nref(N, K, j)
+            for i_1 = 1:sizeA1
+                tmp = @nref(N, A, i)
+                if isfinite(tmp)
+                    s += tmp
+                    k += 1
+                end
+            end
+            @nref(N, S, j) = s
+            @nref(N, K, j) = k
+        end
+    else
+        # Accumulate to array storage
+        @inbounds @nloops N i A d->(j_d = sizeS_d==1 ? 1 : i_d) begin
+            tmp = @nref(N, A, i)
+            if isfinite(tmp)
+                @nref(N, S, j) += tmp
+                @nref(N, K, j) += 1
+            end
+        end
+    end
+    S, K
+end
+
 # Logical operations
 (.<)(img::AbstractImageDirect, n::Number) = data(img) .< n
 (.>)(img::AbstractImageDirect, n::Number) = data(img) .> n

test/algorithms.jl

@@ -64,11 +64,18 @@ let
     img = Images.colorim(A, "RGB")
     img["limits"] = (0.0, 1.0)
     s12 = sum(img, (1,2))
-    @test colorspace(s12) == "RGB"
-    @test limits(s12) == (0.0,25.0)
+    @test Images.colorspace(s12) == "RGB"
+    @test Images.limits(s12) == (0.0,25.0)
     s3 = sum(img, (3,))
-    @test colorspace(s3) == "Unknown"
-    @test limits(s3) == (0.0,3.0)
+    @test Images.colorspace(s3) == "Unknown"
+    @test Images.limits(s3) == (0.0,3.0)
+    A = [NaN, 1, 2, 3]
+    @test_approx_eq Images.meanfinite(A, 1) [2]
+    A = [NaN 1 2 3;
+         NaN 6 5 4]
+    @test_approx_eq Images.meanfinite(A, 1) [NaN 3.5 3.5 3.5]
+    @test_approx_eq Images.meanfinite(A, 2) [2, 5]'
+    @test_approx_eq Images.meanfinite(A, (1,2)) [3.5]
 end
 
 # Array padding
@@ -137,7 +144,7 @@ Aimg = permutedims(convert(Images.Image, A), [3,1,2])
 @assert approx_equal(Images.imfilter_gaussian(ones(4,4), [5,5]), 1.0)
 
 A = zeros(Int, 9, 9); A[5, 5] = 1
-@test maximum(abs(Images.imfilter_LoG(A, [1,1]) - imlog(1.0))) < EPS
+@test maximum(abs(Images.imfilter_LoG(A, [1,1]) - Images.imlog(1.0))) < EPS
 
 # restriction
 A = reshape(uint16(1:60), 4, 5, 3)