Fixes numerical issues with N0f8 types in CLAHE

Another issue that arose is that we try to simultaneously support canonical JuliaImages image types, as well as the ability to pass in a matrix of raw numbers. Using the function `gray` on a raw number throws an error. I've refactored some key lines so that gray is only utilised on JuliaImages types.

Author: zygmuntszpak

Project.toml

@@ -1,7 +1,7 @@
 name = "ImageContrastAdjustment"
 uuid = "f332f351-ec65-5f6a-b3d1-319c6670881a"
+version = "0.3.13"
 authors = ["Dr. Zygmunt L. Szpak <zygmunt.szpak@gmail.com>"]
-version = "0.3.12"
 
 [deps]
 ImageBase = "c817782e-172a-44cc-b673-b171935fbb9e"
@@ -10,19 +10,27 @@ ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 
 [compat]
+Aqua = "0.8"
 ImageBase = "0.1"
 ImageCore = "0.9.3, 0.10"
+ImageFiltering = "0.7"
+ImageMagick = "1"
 ImageTransformations = "0.8.1, 0.9, 0.10"
+LinearAlgebra = "1"
 Parameters = "0.12"
+StableRNGs = "1.0.4"
+Test = "1"
+TestImages = "1"
 julia = "1"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [targets]
-test = ["Test", "Aqua", "ImageFiltering", "TestImages", "ImageMagick", "LinearAlgebra"]
+test = ["Test", "Aqua", "ImageFiltering", "TestImages", "ImageMagick", "LinearAlgebra", "StableRNGs"]

src/ImageContrastAdjustment.jl

@@ -10,6 +10,9 @@ using Parameters: @with_kw # Same as Base.@kwdef but works on Julia 1.0
 # TODO Relax this to all image color types
 const GenericGrayImage = AbstractArray{<:Union{Number, AbstractGray}}
 
+@inline intensity(x::AbstractGray) = float(gray(x))
+@inline intensity(x::Real) = float(x)
+
 # TODO: port HistogramAdjustmentAPI to ImagesAPI
 include("HistogramAdjustmentAPI/HistogramAdjustmentAPI.jl")
 import .HistogramAdjustmentAPI: AbstractHistogramAdjustmentAlgorithm,

src/algorithms/adaptive_equalization.jl

@@ -417,8 +417,9 @@ function perform_iterative_redistribution!(histogram::AbstractArray, limit::Numb
 end
 
 
+
 function apply_cdf_transform(val::Union{Real,AbstractGray}, minval::Union{Real,AbstractGray}, maxval::Union{Real,AbstractGray}, edges::AbstractArray, cdf::AbstractArray)
-    val, minval, maxval = gray(val), gray(minval), gray(maxval)
+    val, minval, maxval = intensity(val), intensity(minval), intensity(maxval)
 
     first_edge = first(edges)
     inv_step_size = 1 / step(edges)
@@ -494,13 +495,25 @@ function transform_image!(out, img, block_centroid_r, block_centroid_c, block_wi
                       intensity_range, block_cdf)
 end
 
-function transform_interior!(out, img, bounds, block_centroids, block_dimensions, intensity_range, block_cdf)
+@inline integer_store(::Type{T}, x, lo, hi) where {T} = convert(T, clamp(ceil(x), lo, hi))
+@inline clamp_store(::Type{T}, x, lo, hi)   where {T} = convert(T, clamp(x, lo, hi))
+
+storefn(::Type{T}) where {T<:Integer} = integer_store
+storefn(::Type{T}) where {T} = clamp_store
+
+function transform_interior!(out::AbstractArray{T}, img, bounds, block_centroids, block_dimensions, intensity_range, block_cdf) where {T}
     rows, cols = bounds
     block_centroid_r, block_centroid_c = block_centroids
     block_width, block_height = block_dimensions
     minval, maxval = intensity_range
+
+    lo = intensity(minval)
+    hi = intensity(maxval)
+    store = storefn(T)
+
     inv_block_height = 1 / block_height
     inv_block_width = 1 / block_width
+
     for r in rows
         for c in cols
             rᵢ = round(Int, r * inv_block_height)
@@ -512,53 +525,70 @@ function transform_interior!(out, img, bounds, block_centroids, block_dimensions
             T₃ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ + 1, cᵢ + 1]...)
             T₄ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ + 1]...)
             interpolated_val = (1 - t)*(1 - u)*T₁ + t*(1 - u)*T₂  + t*u*T₃ + (1 - t)*u*T₄
-            out[r,c] = eltype(img) <: Integer ? ceil(interpolated_val) : interpolated_val
+            out[r,c] = store(T, interpolated_val, lo, hi)
         end
     end
 end
 
 
-function transform_vertical_strip!(out, img, bounds, block_centroid_r, cᵢ, block_height, intensity_range, block_cdf)
+function transform_vertical_strip!(out::AbstractArray{T}, img, bounds, block_centroid_r, cᵢ, block_height, intensity_range, block_cdf) where {T}
     rows, cols = bounds
     minval, maxval = intensity_range
+
+    lo = intensity(minval)
+    hi = intensity(maxval)
+    store = storefn(T)
+
     inv_block_height = 1 / block_height
+
     for r in rows
         for c in cols
             rᵢ = round(Int, r * inv_block_height)
             t = (r - block_centroid_r[rᵢ]) / block_height
             T₁ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ]...)
             T₂ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ + 1, cᵢ]...)
             interpolated_val = (1 - t)*T₁ + t*T₂
-            out[r,c] = eltype(img) <: Integer ? ceil(interpolated_val) : interpolated_val
+            out[r,c] = store(T, interpolated_val, lo, hi)
         end
     end
 end
 
 
-function transform_horizontal_strip!(out, img, bounds, block_centroid_c, rᵢ, block_width, intensity_range, block_cdf)
+function transform_horizontal_strip!(out::AbstractArray{T}, img, bounds, block_centroid_c, rᵢ, block_width, intensity_range, block_cdf) where {T}
     rows, cols = bounds
     minval, maxval = intensity_range
+
+    lo = intensity(minval)
+    hi = intensity(maxval)
+    store = storefn(T)
+
     inv_block_width = 1 / block_width
+
     for r in rows
         for c in cols
             cᵢ = round(Int, c * inv_block_width)
             u = (c - block_centroid_c[cᵢ]) / block_width
             T₁ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ]...)
             T₂ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ + 1]...)
             interpolated_val = (1 - u)*T₁ + u*T₂
-            out[r,c] = eltype(img) <: Integer ? ceil(interpolated_val) : interpolated_val
+            out[r,c] = store(T, interpolated_val, lo, hi)
         end
     end
 end
 
 
-function transform_corner!(out, img, bounds, rᵢ, cᵢ, intensity_range, block_cdf)
+function transform_corner!(out::AbstractArray{T}, img, bounds, rᵢ, cᵢ, intensity_range, block_cdf) where {T}
     rows, cols = bounds
     minval, maxval = intensity_range
+
+    lo = intensity(minval)
+    hi = intensity(maxval)
+    store = storefn(T)
+
     for r in rows
         for c in cols
             val = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ]...)
-            out[r,c] = eltype(img) <: Integer ? ceil(val) : val
+            out[r,c] = store(T, val, lo, hi)
         end
     end
 end

src/algorithms/common.jl

@@ -3,7 +3,7 @@ function transform_density!(out::GenericGrayImage, img::GenericGrayImage, edges:
     first_newval, last_newval = first(newvals), last(newvals)
     inv_step_size = 1/step(edges)
     function transform(val)
-        val = gray(val)
+        val = intensity(val)
         if val >= last_edge
             return last_newval
         elseif val < first_edge

src/algorithms/matching.jl

@@ -144,7 +144,7 @@ function match_pdf!(img::GenericGrayImage, edges::AbstractArray, pdf::AbstractAr
     first_edge = first(edges)
     last_edge = last(edges)
     map!(img, img) do val
-        val = gray(val)
+        val = intensity(val)
         if isnan(val)
             return val
         else

src/build_histogram.jl

@@ -200,7 +200,7 @@ function build_histogram(img::GenericGrayImage, edges::AbstractRange)
             elseif val < first_edge
                 counts[lb] += 1
             else
-                index = floor(Int, gray((val-first_edge)*inv_step_size)) + 1
+                index = floor(Int, intensity((val-first_edge)*inv_step_size)) + 1
                 counts[index] += 1
             end
         end

test/adaptive_equalization.jl

@@ -122,3 +122,22 @@
     imgeq₂ = adjust_histogram(imgg₂, algo)
     @test norm(imgeq₁ .- imgeq₂) ≈ 0.0
 end
+
+@testset "CLAHE regression: Avoid failure on Gray{N0f8} images " begin 
+    # https://github.com/JuliaImages/ImageContrastAdjustment.jl/issues/64
+    rng = StableRNG(123)
+    img = Gray{N0f8}.([only(rand(rng,1)) for r = 1:600, c = 1:600])
+    algo = AdaptiveEqualization( nbins = 256, minval = 0, maxval = 1, rblocks = 4, cblocks = 4, clip = 0.2)
+    imgeq = adjust_histogram(img, algo)
+    @test minimum(imgeq) >= Gray{N0f8}(0)
+    @test maximum(imgeq) <= Gray{N0f8}(1)
+end 
+
+@testset "CLAHE supports raw UInt8 arrays without conversion failure" begin
+    img = UInt8.([round(Int, 240 + 15 * sin(r / 8)) for r = 1:64, c = 1:64])
+    algo = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 255, rblocks = 8, cblocks = 8, clip = 0.2)
+    imgeq = adjust_histogram(img, algo)
+    @test eltype(imgeq) == UInt8
+    @test minimum(imgeq) >= 0
+    @test maximum(imgeq) <= 255    
+end

test/runtests.jl

@@ -1,5 +1,5 @@
 using ImageContrastAdjustment
-using Test, ImageCore, ImageFiltering, TestImages, LinearAlgebra
+using Test, ImageCore, ImageFiltering, TestImages, LinearAlgebra, StableRNGs
 using Aqua
 
 if Base.VERSION >= v"1.6"