LinearStretching enhancements

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,17 @@
 """
 ```
     LinearStretching <: AbstractHistogramAdjustmentAlgorithm
-    LinearStretching(; minval = 0, maxval = 1)
+    LinearStretching(; [src_minval], [src_maxval], dest_minval = 0, dest_maxval = 1)
+
+    LinearStretching((src_minval, src_maxval) => (dest_minval, dest_maxval))
+    LinearStretching((src_minval, src_maxval))
+    LinearStretching(nothing => (dest_minval, dest_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 [`dest_minval`, `dest_maxval`].
 
 # Details
 
@@ -39,44 +43,131 @@ 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 `dest_minval` and `dest_maxval`
+
+If dest_minval and dest_maxval are specified then intensities are mapped to the range
+[`dest_minval`, `dest_maxval`]. The default values are 0 and 1.
 
-If minval and maxval are specified then intensities are mapped to the range
-[`minval`, `maxval`]. The default values are 0 and 1.
+## Choices for `src_minval` and `src_maxval`
+
+`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 `(dest_minval, dest_maxval)` if it exceeds that.
 
 # Example
 
 ```julia
-using ImageContrastAdjustment, ImageView, TestImages
+using ImageContrastAdjustment, TestImages
 
 img = testimage("mandril_gray")
-imgo = adjust_histogram(img, LinearStretching(minval = 0, maxval = 1))
+imgo = adjust_histogram(img, LinearStretching(nothing=>(0, 1)))
 
 ```
 
 # 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
+    dest_minval::T = 0.0f0
+    dest_maxval::T = 1.0f0
+    minval::T = nothing
+    maxval::T = nothing
+    function LinearStretching(src_minval::T1,
+                              src_maxval::T2,
+                              dest_minval::T3,
+                              dest_maxval::T4,
+                              minval::T5=nothing,
+                              maxval::T6=nothing) 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)
+            dest_minval = minval
+            Base.depwarn("deprecated: use `dest_minval` for keyword `minval`", :LinearStretching)
+        end
+        if !isnothing(maxval)
+            dest_maxval = maxval
+            Base.depwarn("deprecated: use `dest_maxval` for keyword `maxval`", :LinearStretching)
+        end
+
+        dest_minval <= dest_maxval || throw(ArgumentError("dest_minval $dest_minval should be less than dest_maxval $dest_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, dest_minval), convert(T, dest_maxval),
+               convert(T, dest_minval), convert(T, dest_maxval))
+    end
+end
+function LinearStretching(rangemap::Pair{Tuple{T1, T2}, Tuple{T3, T4}}) where {T1, T2, T3, T4}
+    LinearStretching(rangemap.first..., rangemap.second...)
+end
+function LinearStretching(rangemap::Pair{Nothing, Tuple{T3, T4}}) where {T3, T4}
+    LinearStretching(nothing, nothing, rangemap.second...)
+end
+function LinearStretching(src_range::Tuple{T1, T2}) where {T1, T2}
+    LinearStretching(src_minval=src_range[1], src_maxval=src_range[2])
 end
 
 function (f::LinearStretching)(out::GenericGrayImage, img::GenericGrayImage)
-    src_minval = minfinite(img)
-    src_maxval = maxfinite(img)
+    img_min, img_max = minfinite(img), maxfinite(img)
+    src_minval = isnothing(f.src_minval) ? img_min : f.src_minval
+    src_maxval = isnothing(f.src_maxval) ? img_max : f.src_maxval
+    dest_minval = f.dest_minval
+    dest_maxval = f.dest_maxval
     T = eltype(out)
-    out .= img
-    map!(out,out) do val
+
+    # 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 = convert(floattype(T), (dest_maxval - dest_minval) / (src_maxval - src_minval))
+    o = convert(floattype(T), (src_minval*dest_maxval - src_maxval*dest_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 `(src_minval, src_maxval)`
+    out_minval = r * img_min - o
+    out_maxval = r * img_max - o
+    do_clamp = (out_minval < dest_minval) || (out_maxval > dest_maxval)
+
+    # strip the colorant type to hit faster clamp version
+    dest_minval = convert(floattype(eltype(typeof(out_minval))), dest_minval)
+    dest_maxval = convert(floattype(eltype(typeof(out_maxval))), dest_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, dest_minval, dest_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, dest_minval, dest_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,14 @@
 @testset "Linear Stretching" begin
 
+    @test LinearStretching() === LinearStretching(nothing, nothing, 0.0f0, 1.0f0)
+    @test LinearStretching(src_minval=0.1f0, src_maxval=0.9f0, dest_minval=0.0f0, dest_maxval=1.0f0) ===
+          LinearStretching(0.1f0, 0.9f0, 0.0f0, 1.0f0)
+    @test LinearStretching((0.1f0, 0.9f0)=>(0.2f0, 0.8f0)) === LinearStretching(0.1f0, 0.9f0, 0.2f0, 0.8f0)
+    @test LinearStretching(nothing=>(0.2f0, 0.8f0)) === LinearStretching((nothing, nothing)=>(0.2f0, 0.8f0))
+    @test LinearStretching((0.1f0, 0.9f0)) === LinearStretching(0.1f0, 0.9f0, 0.0f0, 1.0f0)
+    @test_throws MethodError LinearStretching(0.1f0, 0.9f0)
+    @test_throws MethodError LinearStretching((0.1f0, 0.9f0), (0.0f0, 1.0f0))
+
     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 +18,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 +28,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,23 +38,29 @@
         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
 
     for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
@@ -65,12 +80,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(dest_minval = 0.1)
+        @test LinearStretching(maxval = 0.9) === LinearStretching(dest_maxval = 0.9)
+        @test LinearStretching(minval = 0.1, maxval = 0.9) === LinearStretching(dest_minval = 0.1, dest_maxval = 0.9)
+    end
 end