From 0cfa1208c5fd4daa654286d06953c2705f43ea0f Mon Sep 17 00:00:00 2001 From: Johnny Chen Date: Fri, 27 Mar 2020 08:29:10 +0800 Subject: [PATCH] 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 | 2 +- src/algorithms/common.jl | 4 - src/algorithms/linear_stretching.jl | 167 +++++++++++++++++++++++++--- test/linear_stretching.jl | 69 +++++++++++- 4 files changed, 213 insertions(+), 29 deletions(-) diff --git a/Project.toml b/Project.toml index f68632b..2f438ab 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ImageContrastAdjustment" uuid = "f332f351-ec65-5f6a-b3d1-319c6670881a" authors = ["Dr. Zygmunt L. Szpak "] -version = "0.3.4" +version = "0.3.5" [deps] ColorVectorSpace = "c3611d14-8923-5661-9e6a-0046d554d3a4" diff --git a/src/algorithms/common.jl b/src/algorithms/common.jl index 6a64d54..146d287 100644 --- a/src/algorithms/common.jl +++ b/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 diff --git a/src/algorithms/linear_stretching.jl b/src/algorithms/linear_stretching.jl index f70baef..dc482b0 100644 --- a/src/algorithms/linear_stretching.jl +++ b/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}) diff --git a/test/linear_stretching.jl b/test/linear_stretching.jl index 562b57b..4470a51 100644 --- a/test/linear_stretching.jl +++ b/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,7 +52,7 @@ 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)])) @@ -37,17 +60,45 @@ # 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,7 +116,7 @@ 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 @@ -73,4 +124,10 @@ @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