LinearStretching enhancements (#28)

LinearStretching enhancements

* expose `src_minval` and `src_maxval` (fixes #27 )
* deprecate `minval`/`maxval` in favor of `dst_minval`/`dst_maxval`
* introduce a new convenient API: `LinearStretching((src_minval, src_maxval) => (dst_minval, dst_maxval))`
* correctly handles the clamp on output intensity
* expose `no_clamp` option to disable the clamp to `[dst_minval, dst_maxval]`: `FixedPoint` inputs are still clamped wrt `[typemin(T), typemax(T)]` to avoid ArgumentError
* performance tweak: 2x-20x performance boost

Project.toml

@@ -1,7 +1,7 @@
 name = "ImageContrastAdjustment"
 uuid = "f332f351-ec65-5f6a-b3d1-319c6670881a"
 authors = ["Dr. Zygmunt L. Szpak <[email protected]>"]
-version = "0.3.4"
+version = "0.3.5"
 
 [deps]
 ColorVectorSpace = "c3611d14-8923-5661-9e6a-0046d554d3a4"

src/algorithms/common.jl

@@ -57,7 +57,3 @@ function cdf2pdf!(pdf::AbstractArray, cdf::AbstractArray)
         pdf[i] = cdf[i] - cdf[i-1]
     end
 end
-
-function linear_stretch(x, A, B, a, b)
-    return (x-A) * ((b-a)/(B-A)) + a
-end

src/algorithms/linear_stretching.jl

@@ -2,13 +2,19 @@
 """
 ```
     LinearStretching <: AbstractHistogramAdjustmentAlgorithm
-    LinearStretching(; minval = 0, maxval = 1)
+    LinearStretching(; [src_minval], [src_maxval],
+                       dst_minval=0, dst_maxval=1,
+                       no_clamp=false)
+
+    LinearStretching((src_minval, src_maxval) => (dst_minval, dst_maxval))
+    LinearStretching((src_minval, src_maxval) => nothing)
+    LinearStretching(nothing => (dst_minval, dst_maxval))
 
     adjust_histogram([T,] img, f::LinearStretching)
     adjust_histogram!([out,] img, f::LinearStretching)
 ```
 
-Returns an image where the range of the intensities spans the interval [`minval`, `maxval`].
+Returns an image where the range of the intensities spans the interval [`dst_minval`, `dst_maxval`].
 
 # Details
 
@@ -39,44 +45,169 @@ channel are stretched to the specified range. The modified Y channel is then
 combined with the I and Q channels and the resulting image converted to the same
 type as the input.
 
-## Choices for `minval` and `maxval`
+## Choices for `dst_minval` and `dst_maxval`
+
+If destination value range `dst_minval` and `dst_maxval` are specified then intensities are
+mapped to the range [`dst_minval`, `dst_maxval`]. The default values are 0 and 1.
+
+## Choices for `src_minval` and `src_maxval`
 
-If minval and maxval are specified then intensities are mapped to the range
-[`minval`, `maxval`]. The default values are 0 and 1.
+The source value range `src_minval` and `src_maxval` specifies the intensity range of input
+image. By default, the values are `extrema(img)` (finite). If custom values are provided,
+the output intensity value will be clamped to range `[dst_minval, dst_maxval]` if it exceeds that.
+
+## `no_clamp`
+
+Setting `no_clamp=true` to disable the automatic clamp even if the output intensity value
+exceeds the range `[dst_minval, dst_maxval]`. Note that a clamp is still applied for types
+that has limited value range, for example, if the input eltype is `N0f8`, then the output will
+be clamped to `[0.0N0f8, 1.0N0f8]` even if `no_clamp==true`.
 
 # Example
 
 ```julia
-using ImageContrastAdjustment, ImageView, TestImages
+using ImageContrastAdjustment, TestImages
 
 img = testimage("mandril_gray")
-imgo = adjust_histogram(img, LinearStretching(minval = 0, maxval = 1))
+# Stretches the contrast in `img` so that it spans the unit interval.
+imgo = adjust_histogram(img, LinearStretching(dst_minval = 0, dst_maxval = 1))
+```
 
+For convenience, Constructing a `LinearStretching` object using `Pair` is also supported
+
+```julia
+# these two constructors are equivalent
+LinearStretching(src_minval=0.1, src_maxval=0.9, dst_minval=0.05, dst_maxval=0.95)
+LinearStretching((0.1, 0.9) => (0.05, 0.95))
+
+# replace the part with `nothing` to use default values, e.g.,
+# specify only destination value range
+LinearStretching(nothing => (0.05, 0.95))
+# specify only source value range and use default destination value range, i.e., (0, 1)
+LinearStretching((0.1, 0.9) => nothing)
 ```
 
 # References
 1. W. Burger and M. J. Burge. *Digital Image Processing*. Texts in Computer Science, 2016. [doi:10.1007/978-1-4471-6684-9](https://doi.org/10.1007/978-1-4471-6684-9)
 
 """
-@with_kw struct LinearStretching{T₁ <: Union{Real,AbstractGray},
-                                    T₂ <: Union{Real,AbstractGray}} <: AbstractHistogramAdjustmentAlgorithm
-    minval::T₁ = 0.0
-    maxval::T₂ = 1.0
+@with_kw struct LinearStretching{T} <: AbstractHistogramAdjustmentAlgorithm
+    src_minval::T = nothing
+    src_maxval::T = nothing
+    dst_minval::T = 0.0f0
+    dst_maxval::T = 1.0f0
+    minval::T = nothing
+    maxval::T = nothing
+    no_clamp::Bool = false
+    function LinearStretching(src_minval::T1,
+                              src_maxval::T2,
+                              dst_minval::T3,
+                              dst_maxval::T4,
+                              minval::T5=nothing,
+                              maxval::T6=nothing,
+                              no_clamp::Bool=false) where {T1 <: Union{Nothing,Real,AbstractGray},
+                                                 T2 <: Union{Nothing,Real,AbstractGray},
+                                                 T3 <: Union{Nothing,Real,AbstractGray},
+                                                 T4 <: Union{Nothing,Real,AbstractGray},
+                                                 T5 <: Union{Nothing,Real,AbstractGray},
+                                                 T6 <: Union{Nothing,Real,AbstractGray}}
+        # in order to deprecate old fields we have to introduce new fields if we still want to use @with_kw
+        # https://github.com/JuliaImages/ImageContrastAdjustment.jl/pull/28#discussion_r395751301
+        if !isnothing(minval)
+            dst_minval = minval
+            Base.depwarn("deprecated: use `dst_minval` for keyword `minval`", :LinearStretching)
+        end
+        if !isnothing(maxval)
+            dst_maxval = maxval
+            Base.depwarn("deprecated: use `dst_maxval` for keyword `maxval`", :LinearStretching)
+        end
+
+        dst_minval <= dst_maxval || throw(ArgumentError("dst_minval $dst_minval should be less than dst_maxval $dst_maxval"))
+        if !(isnothing(src_minval) || isnothing(src_maxval))
+            src_minval <= src_maxval || throw(ArgumentError("src_minval $src_minval should be less than src_maxval $src_maxval"))
+        end
+        T = promote_type(T1, T2, T3, T4, T5, T6)
+        new{T}(convert(T, src_minval), convert(T, src_maxval),
+               convert(T, dst_minval), convert(T, dst_maxval),
+               convert(T, dst_minval), convert(T, dst_maxval),
+               no_clamp)
+    end
+end
+function LinearStretching(rangemap::Pair{Tuple{T1, T2}, Tuple{T3, T4}}; no_clamp=false) where {T1, T2, T3, T4}
+    LinearStretching(rangemap.first..., rangemap.second..., nothing, nothing, no_clamp)
+end
+function LinearStretching(rangemap::Pair{Nothing, Tuple{T3, T4}}; no_clamp=false) where {T3, T4}
+    LinearStretching(nothing, nothing, rangemap.second..., nothing, nothing, no_clamp)
+end
+function LinearStretching(rangemap::Pair{Tuple{T1, T2}, Nothing}; no_clamp=false) where {T1, T2}
+    LinearStretching(src_minval=rangemap.first[1], src_maxval=rangemap.first[2], no_clamp=no_clamp)
 end
 
 function (f::LinearStretching)(out::GenericGrayImage, img::GenericGrayImage)
-    src_minval = minfinite(img)
-    src_maxval = maxfinite(img)
     T = eltype(out)
-    out .= img
-    map!(out,out) do val
+    FT = eltype(floattype(T))
+    img_min, img_max = minfinite(img), maxfinite(img)
+    # explicit annotation is needed because the ?: line mixes three value types:
+    # Nothing, T, and typeof(f.src_minval)
+    src_minval::FT = isnothing(f.src_minval) ? img_min : f.src_minval
+    src_maxval::FT = isnothing(f.src_maxval) ? img_max : f.src_maxval
+    dst_minval::FT = f.dst_minval
+    dst_maxval::FT = f.dst_maxval
+
+    # the kernel operation `r * x - o` is equivalent to `(x-A) * ((b-a)/(B-A)) + a`
+    # precalculate these and make inner loop contains only multiplication and addition
+    # to get better performance
+    r = (dst_maxval - dst_minval) / (src_maxval - src_minval)
+    o = (src_minval*dst_maxval - src_maxval*dst_minval) / (src_maxval - src_minval)
+
+    if 1 ≈ r && 0 ≈ o
+        # when image intensity is already adjusted, there's no need to do it again
+        # it's a trivial but common case in practice
+        out === img || (out .= img)
+        return out
+    end
+
+    # In most cases, we don't need to clamp the output
+    # this is only used when user specifies custom parameters
+    out_minval = r * img_min - o
+    out_maxval = r * img_max - o
+    do_clamp = f.no_clamp ? false : (out_minval < dst_minval) || (out_maxval > dst_maxval)
+    # re-adjust clamp option to avoid ArgumentError
+    if !do_clamp && (out_minval < typemin(T) || out_maxval > typemax(T))
+        do_clamp = true
+        dst_minval = convert(typeof(dst_minval), typemin(eltype(T)))
+        dst_maxval = convert(typeof(dst_maxval), typemax(eltype(T)))
+    else
+        do_clamp = do_clamp || out_minval < typemin(T) || out_maxval > typemax(T)
+        dst_minval = max(dst_minval, convert(typeof(dst_minval), typemin(eltype(T))))
+        dst_maxval = min(dst_maxval, convert(typeof(dst_maxval), typemax(eltype(T))))
+    end
+
+    # although slightly faster, clamp before linear stretching has some rounding issue, e.g,
+    # it's likely to get results like -9.313226f-10, which can't be directly converted to N0f8
+    # thus the line below is commented out and not used
+    # do_clamp && (img = clamp.(img, src_minval, src_maxval))
+
+    # tweak the performance of FixedPoint by fusing operations into one broadcast
+    # for Float32 the fallback implementation is faster
+    if eltype(T) <: FixedPoint
+        # ?: is faster than if-else
+        @. out = do_clamp ? clamp(r * img - o, dst_minval, dst_maxval) : r * img - o
+        return out
+    end
+
+    # fallback implementation
+    @inbounds @simd for p in eachindex(img)
+        val = img[p]
         if isnan(val)
-            return val
+            out[p] = val
         else
-            newval = linear_stretch(val, src_minval, src_maxval, f.minval, f.maxval)
-            return  T <: Integer ? round(Int, newval ) : newval
+            newval = r * val - o
+            do_clamp && (newval = clamp(newval, dst_minval, dst_maxval))
+            out[p] = T <: Integer ? round(Int, newval) : newval
         end
     end
+    out
 end
 
 function (f::LinearStretching)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})

test/linear_stretching.jl

@@ -1,5 +1,28 @@
 @testset "Linear Stretching" begin
 
+    @testset "constructors" begin
+        @test LinearStretching() === LinearStretching(nothing, nothing, 0.0f0, 1.0f0, nothing, nothing, false)
+        @test LinearStretching(src_minval=0.1f0, src_maxval=0.9f0, dst_minval=0.05f0, dst_maxval=0.95f0, no_clamp=true) ===
+              LinearStretching(0.1f0, 0.9f0, 0.05f0, 0.95f0, nothing, nothing, true)
+
+        @test LinearStretching((0.1f0, 0.9f0)=>(0.2f0, 0.8f0)) === LinearStretching(0.1f0, 0.9f0, 0.2f0, 0.8f0)
+        @test LinearStretching((0.1f0, 0.9f0)=>(0.2f0, 0.8f0), no_clamp=true) ===
+              LinearStretching(0.1f0, 0.9f0, 0.2f0, 0.8f0, nothing, nothing, true)
+
+        @test LinearStretching(nothing=>(0.2f0, 0.8f0)) === LinearStretching((nothing, nothing)=>(0.2f0, 0.8f0))
+        @test LinearStretching(nothing=>(0.2f0, 0.8f0), no_clamp=true) ===
+              LinearStretching((nothing, nothing)=>(0.2f0, 0.8f0), no_clamp=true)
+
+        @test LinearStretching((0.1f0, 0.9f0)=>nothing, no_clamp=true) ===
+              LinearStretching(0.1f0, 0.9f0, 0.0f0, 1.0f0, nothing, nothing, true)
+
+        @test_throws MethodError LinearStretching(0.1f0, 0.9f0)
+        @test_throws MethodError LinearStretching((0.1f0, 0.9f0), (0.0f0, 1.0f0))
+        @test_throws ArgumentError LinearStretching((0.9f0, 0.1f0)=>nothing)
+        @test_throws ArgumentError LinearStretching(nothing=>(0.9f0, 0.1f0))
+        @test_throws ArgumentError LinearStretching((0.9f0, 0.1f0)=>(0.9f0, 0.1f0))
+    end
+
     for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
         #=
         Stretching an image consisting of a linear ramp should not change the image
@@ -9,7 +32,7 @@
         img = T.(collect(reshape(1/100:1/100:1, 10, 10)))
         minval = minimum(img)
         maxval = maximum(img)
-        ret =  adjust_histogram(img, LinearStretching(minval = minval, maxval = maxval))
+        ret =  adjust_histogram(img, LinearStretching(nothing=>(minval, maxval)))
         if T <: Gray{Float32} || T <: Gray{Float64}
             @test all(map((i, r) -> isapprox(i, r), img, ret))
         else
@@ -19,7 +42,7 @@
         # Verify that NaN is also handled correctly.
         if T <: Gray{Float32} || T <: Gray{Float64}
             img[10] = NaN
-            ret = adjust_histogram(img, LinearStretching(minval = minval, maxval = maxval))
+            ret = adjust_histogram(img, LinearStretching(nothing=>(minval, maxval)))
             @test isapprox(first(img), first(ret))
             @test isapprox(last(img), last(ret))
             @test isnan(ret[10])
@@ -29,25 +52,53 @@
         img = T.(collect(reshape(1/100:1/100:1, 10, 10)))
         minval = minimum(img)
         maxval = maximum(img)
-        ret =  adjust_histogram(img, LinearStretching(minval = 0, maxval = 1))
+        ret =  adjust_histogram(img, LinearStretching(nothing=>(0, 1)))
         @test isapprox(0, first(ret))
         @test isapprox(1, last(ret))
         @test isapprox(0, minimum(ret[.!isnan.(ret)]))
         @test isapprox(1, maximum(ret[.!isnan.(ret)]))
 
         # Verify that the return type matches the input type.
         img = T.(testimage("mandril_gray"))
-        ret = adjust_histogram(img, LinearStretching(minval = 0, maxval = 1))
+        ret = adjust_histogram(img, LinearStretching(nothing=>(0, 1)))
         @test eltype(ret) == eltype(img)
         @test isapprox(0, minimum(ret))
         @test isapprox(1, maximum(ret))
 
-        ret = adjust_histogram(img, LinearStretching(minval = 0.2, maxval = 0.8))
+        ret = adjust_histogram(img, LinearStretching(nothing=>(0.2, 0.8)))
         @test eltype(ret) == eltype(img)
         @test isapprox(0.2, minimum(ret))
         @test isapprox(0.8, maximum(ret))
+
+        # Verify that results are correctly clamped to [0.2, 0.9] if it exceeds the range
+        ret = adjust_histogram(img, LinearStretching((0.1, 0.8)=>(0.2, 0.9)))
+        @test eltype(ret) == eltype(img)
+        @test isapprox(T(0.2), minimum(ret))
+        @test isapprox(T(0.9), maximum(ret), atol=1e-2)
     end
 
+    # Verify that no_clamp option handles different input types correctly without ArgumentError
+    img = Float32.(testimage("mandril_gray"))
+    ret_clamp = adjust_histogram(img, LinearStretching((0.3, 0.8)=>(0.1, 1.1), no_clamp=true))
+    @test eltype(ret_clamp) == eltype(img)
+    @test isapprox(-0.5f0, minimum(ret_clamp))
+    @test isapprox(1.272549f0, maximum(ret_clamp))
+    ret_noclamp = adjust_histogram(img, LinearStretching((0.3, 0.8)=>(0.1, 1.1), no_clamp=false))
+    @test eltype(ret_noclamp) == eltype(img)
+    @test isapprox(0.1f0, minimum(ret_noclamp))
+    @test isapprox(1.1f0, maximum(ret_noclamp))
+    # when no_clamp==true, the output is still clamped by (typemin(T), typemax(T))
+    img = N0f8.(testimage("mandril_gray"))
+    ret_clamp = adjust_histogram(img, LinearStretching((0.3, 0.8)=>(0.1, 1.1), no_clamp=true))
+    @test eltype(ret_clamp) == eltype(img)
+    @test isapprox(0.0N0f8, minimum(ret_clamp))
+    @test isapprox(1.0N0f8, maximum(ret_clamp))
+    ret_noclamp = adjust_histogram(img, LinearStretching((0.3, 0.8)=>(0.1, 1.1), no_clamp=false))
+    @test eltype(ret_noclamp) == eltype(img)
+    @test isapprox(0.1N0f8, minimum(ret_noclamp))
+    @test isapprox(1.0N0f8, maximum(ret_noclamp))
+    @test ret_clamp != ret_noclamp
+
     for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
         #=
         Create a color image that spans a narrow graylevel range.  Then
@@ -65,12 +116,18 @@
         verify that all 32 bins have non-zero counts. This will confirm that
         the dynamic range of the original image has been increased.
         =#
-        ret = adjust_histogram(img, LinearStretching(minval = 0, maxval = 1))
+        ret = adjust_histogram(img, LinearStretching(nothing=>(0, 1)))
         edges, counts_after = build_histogram(ret, 32, minval = 0, maxval = 1)
         nonzero_after = sum(counts_after .!= 0)
         @test nonzero_before < nonzero_after
         @test nonzero_after == 32
         @test eltype(img) == eltype(ret)
     end
 
+    @testset "deprecations" begin
+        @info "four depwarns are expected"
+        @test LinearStretching(minval = 0.1) === LinearStretching(dst_minval = 0.1)
+        @test LinearStretching(maxval = 0.9) === LinearStretching(dst_maxval = 0.9)
+        @test LinearStretching(minval = 0.1, maxval = 0.9) === LinearStretching(dst_minval = 0.1, dst_maxval = 0.9)
+    end
 end