From 4feea1b7d86e121d21e7d555d68a86638dd0c2bb Mon Sep 17 00:00:00 2001 From: Zygmunt Szpak Date: Mon, 13 Sep 2021 11:04:56 +0930 Subject: [PATCH 1/4] Adds PiecewiseLinearStretching --- docs/Project.toml | 1 + docs/src/reference.md | 5 + src/ImageContrastAdjustment.jl | 2 + src/algorithms/piecewise_linear_stretching.jl | 179 ++++++++++++++++++ test/piecewise_linear_stretching.jl | 170 +++++++++++++++++ test/runtests.jl | 1 + 6 files changed, 358 insertions(+) create mode 100644 src/algorithms/piecewise_linear_stretching.jl create mode 100644 test/piecewise_linear_stretching.jl diff --git a/docs/Project.toml b/docs/Project.toml index 8a682ee..c6c3abf 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,7 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a" ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31" diff --git a/docs/src/reference.md b/docs/src/reference.md index 72d8af4..a0eb2d2 100644 --- a/docs/src/reference.md +++ b/docs/src/reference.md @@ -48,3 +48,8 @@ Matching ```@docs MidwayEqualization ``` + +### PiecewiseLinearStretching +```@docs +PiecewiseLinearStretching +``` \ No newline at end of file diff --git a/src/ImageContrastAdjustment.jl b/src/ImageContrastAdjustment.jl index 101d2b0..97c6cf3 100644 --- a/src/ImageContrastAdjustment.jl +++ b/src/ImageContrastAdjustment.jl @@ -24,6 +24,7 @@ include("algorithms/contrast_stretching.jl") include("algorithms/gamma_correction.jl") include("algorithms/matching.jl") include("algorithms/midway_equalization.jl") +include("algorithms/piecewise_linear_stretching.jl") include("compat.jl") export @@ -35,6 +36,7 @@ export GammaCorrection, LinearStretching, ContrastStretching, + PiecewiseLinearStretching, build_histogram, adjust_histogram, adjust_histogram! diff --git a/src/algorithms/piecewise_linear_stretching.jl b/src/algorithms/piecewise_linear_stretching.jl new file mode 100644 index 0000000..2f439fa --- /dev/null +++ b/src/algorithms/piecewise_linear_stretching.jl @@ -0,0 +1,179 @@ +""" +``` + PiecewiseLinearStretching <: AbstractHistogramAdjustmentAlgorithm + PiecewiseLinearStretching(; src_intervals = (0.0 => 0.1, 0.1 => 0.9, 0.9 => 1.0), + dst_intervals = (0.0 => 0.2, 0.2 => 1.0, 1.0 => 1.0)) + + adjust_histogram([T,] img, f::PiecewiseLinearStretching) + adjust_histogram!([out,] img, f::PiecewiseLinearStretching) +``` + +Returns an image for which the intensity range was partitioned into subintervals (specified +by `src_intervals`) and mapped linearly to a new set of subintervals +(specified by `dst_intervals`). + +# Details + +Piecewise linear stretching is a contrast enhancing transformation that is used to modify +the dynamic range of an image. The concept involves partitioning the intensity +range of an image, ``I = [A,B]`` into a sequence of ``m`` disjoint intervals + +```math +I_1 = [A_1,A_2], I_2 = (A_2, A_3], \\ldots, I_{m} = (A_m,A_{m+1}] +``` +which encompass the entire range of intensity values such that ``\\sum_j |I_j | = | I |``. +One subsequently constructs a concomitant new set of disjoint intervals +```math +I_{1}^{′} = [a_1,a_2], I_{2}^{′} = (a_2, a_3], \\ldots, I_{m}^{′} = (a_m,a_{m+1}] +``` +such that ``\\sum_j |I_{j}^{′} | = | I^{′} |``, where ``I^{′} = [a,b]``. Finally, one introduces a +sequence of mappings ``f_j : I_j \\to I_{j}^{′}`` ``(j = 1, \\ldots, m)`` given by +```math +f_j(x) = (x-A_{j}) \\frac{a_{j+1}-a_j}{A_{j+1}-A_{j}} + a_j. +``` +which linearly map intensities from the interval ``I_j`` to the interval ``I_{j}^{′}`` . + +# Options + +Various options for the parameters of the `adjust_histogram` and +`PiecewiseLinearStretching` type are described in more detail below. + +## Choices for `img` + +The function can handle a variety of input types. The returned image depends on the input +type. + +For colored images, the input is converted to the [YIQ](https://en.wikipedia.org/wiki/YIQ) +type and the intensities of the Y 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 `src_intervals` + +The `src_intervals` must be specified as an `ntuple` of pairs. For example, +`src_intervals = (0.0 => 0.2, 0.2 => 0.8, 0.8 => 1.0)` partitions the unit +interval into three subintervals, [0.0, 0.2], (0.2, 0.8] and (0.8, 1.0]. + +To specify a single interval you must include a trailing comma, e.g. +`src_intervals = (0.0 => 1.0,)`. + +The start and end of a subinterval in `src_intervals` cannot be equal. For example, +`src_intervals = (0.0 => 0.5, 0.5 => 0.5, 0.5 => 1.0)` will result in an `ArgumentError` +because of the pair `0.5 => 0.5`. + +## Choices for `dst_intervals` + +The `dst_intervals` must be specified as an `ntuple` of pairs and there needs to be a +corresponding pair for each pair in `src_intervals`. + +Unlike `src_intervals`, the start and end of an interval in `src_intervals` can be equal. +For example, the combination `src_intervals = (0.0 => 0.5, 0.5 => 1.0)` and `dst_intervals = +(0.0 => 0.0 , 1.0 => 1.0)` will produce a binary image. + +Care must be taken to ensure that the subintervals specified in `dst_intervals` are +representable within the type of image on which the piecewise linear stretching will be +applied. For example, if you specify a negative output interval like `dst_intervals = (1.0 +=> -1.0),` and try to apply the transformation on an image with type `Gray{N0f8}`, you will +receive a DomainError. + +## Saturating pixels + +If the specified `src_intervals` don't span the entire range of intensities in an image, +then intensities that fall outside `src_intervals` are automatically set to the boundary +values of the `dst_intervals`. For example, given the combination `src_intervals = (0.1 => +0.9, )` and `dst_intervals = (0.0 => 1.0)`, intensities below `0.1` are set to zero, and +intensities above `0.9` are set to one. + +# Example + +```julia +using ImageContrastAdjustment, TestImages + +img = testimage("mandril_gray") +# Stretches the contrast in `img` so that it spans the unit interval. +f = PiecewiseLinearStretching(src_intervals = (0.0 => 0.2, 0.2 => 0.8, 0.8 => 1.0), + dst_intervals = (0.0 => 0.4, 0.4 => 0.9, 0.9 => 1.0)) +imgo = adjust_histogram(img, f) +``` + +# 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 PiecewiseLinearStretching{T} <: AbstractHistogramAdjustmentAlgorithm where T <: NTuple{N, Pair{<: Union{Real, AbstractGray}, <: Union{Real, AbstractGray}}} where {N} + src_intervals::T = (0.0 => 0.1, 0.1 => 0.9, 0.9 => 1.0) + dst_intervals::T = (0.0 => 0.2, 0.2 => 1.0, 1.0 => 1.0) + function PiecewiseLinearStretching(src_intervals::T₁, dst_intervals::T₂) where {T₁ <: NTuple{N₁, Pair{<: Union{Real, AbstractGray}, <: Union{Real, AbstractGray}}} where {N₁}, T₂ <: NTuple{N₂, Pair{<: Union{Real, AbstractGray}, <: Union{Real, AbstractGray}}} where {N₂}} + length(src_intervals) == length(dst_intervals) || throw(ArgumentError("The dst_intervals must define an interval pair for each consecutive src_interval pair, i.e. length(src_intervals) == length(dst_intervals).")) + for src_pair in src_intervals + if first(src_pair) == last(src_pair) + throw(ArgumentError("The start and end of an interval in src_intervals cannot be equal.")) + end + end + T = promote_type(T₁, T₂) + new{T}(src_intervals, dst_intervals) + end +end + +function (f::PiecewiseLinearStretching)(out::GenericGrayImage, img::GenericGrayImage) + # Verify that the specified dst_intervals fall inside the permissible range + # of the output type. + T = eltype(out) + for (a,b) in f.dst_intervals + p = typemin(T) + r = typemax(T) + if !(p <= a <= r) || !(p <= b <= r) + throw(DomainError("The specified interval [$a, $b] falls outside the permissible range [$p, $r] associated with $T")) + end + end + # A linear mapping of a value x in the interval [A, B] to the interval [a,b] is given by + # the formula f(x) = (x-A) * ((b-a)/(B-A)) + a. The operation r * x - o is equivalent to + # (x-A) * ((b-a)/(B-A)) + a. Precalculate the r and o so that later on an inner loop + # contains only multiplication and addition to get better performance. + N = length(f.src_intervals) + rs = zeros(Float64, N) + os = zeros(Float64, N) + for (src, dst) in zip(pairs(f.src_intervals), pairs(f.dst_intervals)) + n, (A,B) = src + _, (a,b) = dst + rs[n] = (b - a) / (B - A) + os[n] = (A*b - B*a) / (B - A) + end + + # Determine which interval the source pixel fall into, and apply the accompanying linear + # transformation taking into account NaN values. + @inbounds for p in eachindex(img) + val = img[p] + if isnan(val) + out[p] = val + else + is_inside_src_intervals = false + for (n, (A,B)) in pairs(f.src_intervals) + if (A <= val < B) || (A <= val <= B && n == N) + newval = rs[n] * val - os[n] + out[p] = T <: Integer ? round(Int, newval) : newval + is_inside_src_intervals = true + break + end + end + if !is_inside_src_intervals + # Saturate pixels that fall outside the specified src_intervals by setting + # them equal to the start/end of the edge dst_intervals. + a = first(first(f.dst_intervals)) + b = last(last(f.dst_intervals)) + A = first(first(f.src_intervals)) + newval = val < A ? a : b + out[p] = T <: Integer ? round(Int, newval) : newval + end + end + end + return out +end + +function (f::PiecewiseLinearStretching)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3}) + T = eltype(out) + yiq = convert.(YIQ, img) + yiq_view = channelview(yiq) + adjust_histogram!(view(yiq_view,1,:,:), f) + out .= convert.(T, yiq) +end + diff --git a/test/piecewise_linear_stretching.jl b/test/piecewise_linear_stretching.jl new file mode 100644 index 0000000..ba9aa14 --- /dev/null +++ b/test/piecewise_linear_stretching.jl @@ -0,0 +1,170 @@ +@testset "Piecewise Linear Stretching" begin + + @testset "constructors" begin + @test_throws ArgumentError PiecewiseLinearStretching( + src_intervals = ((0.0 => 0.2),(0.2 => 0.2), (0.8 => 1.0)), + dst_intervals = ((0.0 => 0.2),(-1.0 => 0.0),(1 => -2))) + + @test_throws ArgumentError PiecewiseLinearStretching( + src_intervals = ((0.0 => 0.2),(0.2 => 0.8), (0.8 => 1.0)), + dst_intervals = ((0.0 => 0.2), (1 => -2))) + end + @testset "miscellaneous" begin + # We should be able to invert the image. + f = PiecewiseLinearStretching(src_intervals = (0.0 => 1.0,), dst_intervals = (1.0 => 0.0,)) + for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) + img = T.(testimage("mandril_gray")) + ret = adjust_histogram(img, f) + @test ret ≈ (1 .- img) + end + + # We should be able to binarize the image. + f = PiecewiseLinearStretching( + src_intervals = ((0.0 => 0.5), (0.5 => 1.0)), + dst_intervals = ((0.0 => 0.0),(1.0 => 1.0))) + for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) + img = T.(testimage("mandril_gray")) + + ret = adjust_histogram(img, f) + vals = sort(unique(ret)) + @test first(vals) == 0.0 + @test last(vals) == 1.0 + end + + # Setting the src_intervals equal to the dst_intervals should result in an identity + # mapping (no change to the image) provided that the src_intervals span the entire + # permissible contrast range, i.e. the unit interval. + f = PiecewiseLinearStretching( + src_intervals = ((0.0 => 0.3), (0.3 => 0.8), (0.8 => 1.0)), + dst_intervals = ((0.0 => 0.3), (0.3 => 0.8), (0.8 => 1.0))) + for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) + img = T.(testimage("mandril_gray")) + ret = adjust_histogram(img, f) + @test ret ≈ img + end + + # If the src_intervals equal the dst_intervals, but the src_intervals don't span the + # entire permissible contrast range, i.e. the unit interval, then the resulting + # image will not match the original image because pixels that fall outside the + # specified src_intervals are satuared to the edge value of src_intervals. + f = PiecewiseLinearStretching( + src_intervals = ((0.3 => 0.8),), + dst_intervals = ((0.3 => 0.8),)) + for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) + img = T.(testimage("mandril_gray")) + ret = adjust_histogram(img, f) + @test !isapprox(ret, 1 .- img, atol=1e-1) + 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 + if the specified minimum and maximum values match the minimum and maximum of + the image. + =# + img = T.(collect(reshape(1/100:1/100:1, 10, 10))) + minval = minimum(img) + maxval = maximum(img) + f = PiecewiseLinearStretching(src_intervals = (minval => maxval,), dst_intervals = (minval => maxval,)) + ret = adjust_histogram(img, f) + if T <: Gray{Float32} || T <: Gray{Float64} + @test all(map((i, r) -> isapprox(i, r), img, ret)) + else + @test ret == img + end + + # Verify that NaN is also handled correctly. + if T <: Gray{Float32} || T <: Gray{Float64} + img[10] = NaN + ret = adjust_histogram(img, f) + @test isapprox(first(img), first(ret)) + @test isapprox(last(img), last(ret)) + @test isnan(ret[10]) + end + + # Verify that the smallest and largest values match the specified minval and maxval. + img = T.(collect(reshape(1/100:1/100:1, 10, 10))) + minval = minimum(img) + maxval = maximum(img) + f = PiecewiseLinearStretching(src_intervals = (minval => maxval,), dst_intervals = (0.0 => 1.0,)) + ret = adjust_histogram(img, f) + @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")) + minval = minimum(img) + maxval = maximum(img) + f = PiecewiseLinearStretching(src_intervals = (minval => maxval,), dst_intervals = (0.0 => 1.0,)) + ret = adjust_histogram(img, f) + @test eltype(ret) == eltype(img) + @test isapprox(0, minimum(ret)) + @test isapprox(1, maximum(ret)) + + f = PiecewiseLinearStretching(src_intervals = (0.0 => 1.0,), dst_intervals = (0.2 => 0.8,)) + ret = adjust_histogram(img, f) + @test eltype(ret) == eltype(img) + # We are mapping the interval [0, 1] to [0.2, 0.8] but since the + # smallest intensity in the mandril image need not be zero, and the largest + # need to be one, smallest and largest observed values in the + # image need not be 0.2 and 0.8, but should still be within those bounds. + @test minimum(ret) >= 0.2 + @test maximum(ret) <= 0.8 + + # Verify that results are correctly clamped to [0.2, 0.9] if it exceeds the range + f = PiecewiseLinearStretching(src_intervals = (0.3 => 0.8,), + dst_intervals = (0.2 => 0.9,)) + ret = adjust_histogram(img, f) + @test eltype(ret) == eltype(img) + @test minimum(ret) == T(0.2) + @test maximum(ret) == T(0.9) + end + + img = Float32.(collect(reshape(1/100:1/100:1, 10, 10))) + + f = PiecewiseLinearStretching(src_intervals = ((0.0 => 1.0),), dst_intervals = ((-1.0 => 0.0),)) + ret = adjust_histogram(img, f) + @test eltype(ret) == eltype(img) + # @test minimum(ret) ≈ -1.0 + # @test maximum(ret) ≈ 0.0 + @test isapprox(minimum(ret), -1.0, atol=1e-1) + @test isapprox(maximum(ret), 0.0, atol=1e-1) + @test_throws DomainError adjust_histogram(Gray{N0f8}.(img), f) + + img = Gray{Float32}.(collect(reshape(1/100:1/100:1, 10, 10))) + f = PiecewiseLinearStretching(src_intervals = ((0.0 => 0.2), (0.2 => 0.8), (0.8 => 1.0)), dst_intervals = ((0.0 => 0.2), (-1.0 => 0.0), (1 => -2))) + ret = adjust_histogram(img, f) + @test eltype(ret) == eltype(img) + @test minimum(ret) ≈ -2.0 + @test maximum(ret) ≈ 1.0 + + for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64}) + #= + Create a color image that spans a narrow graylevel range. Then + quantize the 256 bins down to 32 and determine how many bins have + non-zero counts. + =# + imgg = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10]) + img = colorview(RGB,imgg,imgg,imgg) + img = T.(img) + _, counts_before = build_histogram(img,32, minval = 0, maxval = 1) + nonzero_before = sum(counts_before .!= 0) + + #= + Stretch the histogram. Then quantize the 256 bins down to 32 and + verify that all 32 bins have non-zero counts. This will confirm that + the dynamic range of the original image has been increased. + =# + f = PiecewiseLinearStretching(src_intervals = (0.3 => 0.5,), + dst_intervals = (0.0 => 1.0,)) + ret = adjust_histogram(img, f) + 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 + end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 95a908d..2044b73 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,4 +11,5 @@ using Test, ImageCore, ImageFiltering, TestImages, LinearAlgebra include("gamma_adjustment.jl") include("linear_stretching.jl") include("contrast_stretching.jl") + include("piecewise_linear_stretching.jl") end From 3da3ae0de2456b9c0a492f039a1addf98ee587cb Mon Sep 17 00:00:00 2001 From: Zygmunt Szpak Date: Thu, 30 Sep 2021 14:06:21 +0930 Subject: [PATCH 2/4] Refactors PiecewiseLinearStretching --- Project.toml | 8 +- src/ImageContrastAdjustment.jl | 3 + src/algorithms/piecewise_linear_stretching.jl | 218 +++++++++++++----- test/piecewise_linear_stretching.jl | 74 +++--- 4 files changed, 201 insertions(+), 102 deletions(-) diff --git a/Project.toml b/Project.toml index 6e62cd0..1c36adb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,17 +1,21 @@ name = "ImageContrastAdjustment" uuid = "f332f351-ec65-5f6a-b3d1-319c6670881a" -authors = ["Dr. Zygmunt L. Szpak "] -version = "0.3.10" +authors = ["Dr. Zygmunt L. Szpak and contributors"] +version = "0.4.0" [deps] ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795" +IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [compat] ImageCore = "0.9.3" ImageTransformations = "0.8.1, 0.9" +IntervalSets = "0.5.3" Parameters = "0.12" +Reexport = "1.2.2" julia = "1" [extras] diff --git a/src/ImageContrastAdjustment.jl b/src/ImageContrastAdjustment.jl index 97c6cf3..0c0365b 100644 --- a/src/ImageContrastAdjustment.jl +++ b/src/ImageContrastAdjustment.jl @@ -4,7 +4,10 @@ using ImageCore using ImageTransformations: imresize # Where possible we avoid a direct dependency to reduce the number of [compat] bounds using ImageCore.MappedArrays +using IntervalSets using Parameters: @with_kw # Same as Base.@kwdef but works on Julia 1.0 +using Reexport +@reexport using IntervalSets # TODO: port HistogramAdjustmentAPI to ImagesAPI include("HistogramAdjustmentAPI/HistogramAdjustmentAPI.jl") diff --git a/src/algorithms/piecewise_linear_stretching.jl b/src/algorithms/piecewise_linear_stretching.jl index 2f439fa..1fdac6c 100644 --- a/src/algorithms/piecewise_linear_stretching.jl +++ b/src/algorithms/piecewise_linear_stretching.jl @@ -1,16 +1,19 @@ """ ``` PiecewiseLinearStretching <: AbstractHistogramAdjustmentAlgorithm - PiecewiseLinearStretching(; src_intervals = (0.0 => 0.1, 0.1 => 0.9, 0.9 => 1.0), - dst_intervals = (0.0 => 0.2, 0.2 => 1.0, 1.0 => 1.0)) - + PiecewiseLinearStretching((0.0, 0.2, 0.8, 1.0), (0.0, 0.1, 0.7, 1.0) ; saturate = true) + PiecewiseLinearStretching([0.0, 0.2, 0.8, 1.0], [0.0, 0.1, 0.7, 1.0] ; saturate = true) + PiecewiseLinearStretching(; src_knots = (0.0, 0.2, 0.8, 1.0), + dst_knots = (0.0, 0.1, 0.7, 1.0), saturate = true) + PiecewiseLinearStretching(ClosedInterval(0.0, 0.2) => ClosedInterval(0.0, 0.1), + Interval{:open, :closed}(0.2, 0.8) => Interval{:open, :closed}(0.1, 0.7), + Interval{:open, :closed}(0.8, 1,0) => Interval{:open, :closed}(0.7, 1.0)) adjust_histogram([T,] img, f::PiecewiseLinearStretching) adjust_histogram!([out,] img, f::PiecewiseLinearStretching) ``` -Returns an image for which the intensity range was partitioned into subintervals (specified -by `src_intervals`) and mapped linearly to a new set of subintervals -(specified by `dst_intervals`). +When used with `adjust_histogram`, returns an image for which the intensity range was +partitioned into subintervals and mapped linearly to a new set of subintervals. # Details @@ -48,41 +51,65 @@ type and the intensities of the Y channel are stretched to the specified range. 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 `src_intervals` +## Specifying intervals explicitly +One can specify the mappings between source and destination intervals using +the interval types provided by [IntervalSets.jl](https://github.com/JuliaMath/IntervalSets.jl). +For instance, the mappings [0, 0.2] => [0, 0.1], (0.2, 0.8] => (0.1, 0.7] and (0.8, 1.0] => +(0.7, 1.0] can be specified as follows: +```julia +PiecewiseLinearStretching(ClosedInterval(0.0, 0.2) => ClosedInterval(0.0, 0.1), + Interval{:open, :closed}(0.2, 0.8) => Interval{:open, :closed}(0.1, 0.7), + Interval{:open, :closed}(0.8, 1,0) => Interval{:open, :closed}(0.7, 1.0)) +``` +While this manner of specifying intervals permits the most flexibility, it is +cumbersome for many types of mappings one would like to specify in practice. +Hence, one can instead specify the mappings by designating a sequence of knot pairs. -The `src_intervals` must be specified as an `ntuple` of pairs. For example, -`src_intervals = (0.0 => 0.2, 0.2 => 0.8, 0.8 => 1.0)` partitions the unit -interval into three subintervals, [0.0, 0.2], (0.2, 0.8] and (0.8, 1.0]. +## Specifying intervals implicitly +As an alternative to specifying the endpoints of each interval, one can list the knots of +the intervals. For instance, the mappings [0.0, 0.2] => [0, 0.1], (0.2, 0.8] => (0.1, 0.7] +and (0.8, 1.0] => (0.7, 1.0] can be listed as follows: +```julia +# Without the use of keyword arguments. +PiecewiseLinearStretching((0.0, 0.2, 0.8, 1.0), (0.0, 0.1, 0.7, 1.0)) +# With the use of keyword arguments. +PiecewiseLinearStretching(src_knots = (0.0, 0.2, 0.8, 1.0), dst_knots = (0.0, 0.1, 0.7, 1.0)) +``` + +## Constraints on intervals and knots. To specify a single interval you must include a trailing comma, e.g. -`src_intervals = (0.0 => 1.0,)`. +```julia +# Using the ellipses notation provided by IntervalSets.jl +f = PiecewiseLinearStretching((0.0..1.0 => 0.2..1.0,)) +# Using the types provided by IntervalSets.jl +f = PiecewiseLinearStretching(ClosedInterval(0.0, 1.0) => ClosedInterval(0.2, 1.0),) +``` -The start and end of a subinterval in `src_intervals` cannot be equal. For example, -`src_intervals = (0.0 => 0.5, 0.5 => 0.5, 0.5 => 1.0)` will result in an `ArgumentError` -because of the pair `0.5 => 0.5`. +The endpoints of any source interval cannot be equal. For example, `ClosedInterval(0.5, 0.5) +=> ClosedInterval(0.0, 1.0)` is not permitted because the value `0.5` cannot be +simultaneously mapped to `0.0` and `1.0`. There is no such constraint for the destination +interval. For instance, `ClosedInterval(0.0, 1.0) => ClosedInterval(0.5, 0.5)` is permitted +and corresponds to mapping all intensities in the unit interval to the value `0.5`. -## Choices for `dst_intervals` +Care must be taken to ensure that the destination subintervals are representable within the +type of image on which the piecewise linear stretching will be applied. For example, if you +specify a negative output interval like `(0.0..1.0 => -1.0..0,),` and try to apply the +transformation on an image with type `Gray{N0f8}`, you will receive a DomainError. -The `dst_intervals` must be specified as an `ntuple` of pairs and there needs to be a -corresponding pair for each pair in `src_intervals`. - -Unlike `src_intervals`, the start and end of an interval in `src_intervals` can be equal. -For example, the combination `src_intervals = (0.0 => 0.5, 0.5 => 1.0)` and `dst_intervals = -(0.0 => 0.0 , 1.0 => 1.0)` will produce a binary image. +When listing knots, the length of `src_knots` must match the length of `dst_knots`. -Care must be taken to ensure that the subintervals specified in `dst_intervals` are -representable within the type of image on which the piecewise linear stretching will be -applied. For example, if you specify a negative output interval like `dst_intervals = (1.0 -=> -1.0),` and try to apply the transformation on an image with type `Gray{N0f8}`, you will -receive a DomainError. ## Saturating pixels -If the specified `src_intervals` don't span the entire range of intensities in an image, -then intensities that fall outside `src_intervals` are automatically set to the boundary -values of the `dst_intervals`. For example, given the combination `src_intervals = (0.1 => -0.9, )` and `dst_intervals = (0.0 => 1.0)`, intensities below `0.1` are set to zero, and -intensities above `0.9` are set to one. +When listing knots, one has the option of saturating pixels that fall outside the endpoints +of the source knots to the endpoints of the destination knots (this option is enabled by default). For example, +```julia +f = PiecewiseLinearStretching(src_knots = (0.2, 0.8), dst_knots = (0.0, 1.0); saturate = true) +``` +will ensure that any image intensities below `0.2` are set to zero, and any intensities above +`0.8` are set to one, whilst also linearly stretching the interval `[0.2, 0.8]` to `[0.0. 1.0]`. + # Example @@ -90,35 +117,109 @@ intensities above `0.9` are set to one. using ImageContrastAdjustment, TestImages img = testimage("mandril_gray") -# Stretches the contrast in `img` so that it spans the unit interval. -f = PiecewiseLinearStretching(src_intervals = (0.0 => 0.2, 0.2 => 0.8, 0.8 => 1.0), - dst_intervals = (0.0 => 0.4, 0.4 => 0.9, 0.9 => 1.0)) +f = PiecewiseLinearStretching((0.0, 0.2, 0.8, 1.0), (0.0, 0.1, 0.7, 1.0)) imgo = adjust_histogram(img, f) ``` # 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 PiecewiseLinearStretching{T} <: AbstractHistogramAdjustmentAlgorithm where T <: NTuple{N, Pair{<: Union{Real, AbstractGray}, <: Union{Real, AbstractGray}}} where {N} - src_intervals::T = (0.0 => 0.1, 0.1 => 0.9, 0.9 => 1.0) - dst_intervals::T = (0.0 => 0.2, 0.2 => 1.0, 1.0 => 1.0) - function PiecewiseLinearStretching(src_intervals::T₁, dst_intervals::T₂) where {T₁ <: NTuple{N₁, Pair{<: Union{Real, AbstractGray}, <: Union{Real, AbstractGray}}} where {N₁}, T₂ <: NTuple{N₂, Pair{<: Union{Real, AbstractGray}, <: Union{Real, AbstractGray}}} where {N₂}} - length(src_intervals) == length(dst_intervals) || throw(ArgumentError("The dst_intervals must define an interval pair for each consecutive src_interval pair, i.e. length(src_intervals) == length(dst_intervals).")) - for src_pair in src_intervals - if first(src_pair) == last(src_pair) - throw(ArgumentError("The start and end of an interval in src_intervals cannot be equal.")) +struct PiecewiseLinearStretching{T} <: AbstractHistogramAdjustmentAlgorithm where T <: Tuple{Vararg{Pair}} + intervals::T + function PiecewiseLinearStretching(intervals::T₁) where T₁ <: Tuple{Vararg{Pair}} + for (src, _) in intervals + if leftendpoint(src) == rightendpoint(src) + throw(ArgumentError("The start and end of a src interval cannot be equal.")) end end - T = promote_type(T₁, T₂) - new{T}(src_intervals, dst_intervals) + new{T₁}(intervals) end end +function PiecewiseLinearStretching(img::GenericGrayImage) + T = eltype(img) + img_min, img_max = minfinite(img), maxfinite(img) + if (img_min ≈ img_max) + throw(DomainError("The image consists of a single intensity, which gives rise to a degenerate source interval.")) + end + src_knots = (img_min, img_max) + dst_knots - (zero(T), one(T)) + return PiecewiseLinearStretching(src_knots, dst_knots; saturate = true) +end + +function PiecewiseLinearStretching(; src_knots, dst_knots, saturate::Bool = true) + PiecewiseLinearStretching(src_knots, dst_knots; saturate = saturate) +end + +function PiecewiseLinearStretching(src_knots::AbstractVector, dst_knots::AbstractVector; saturate::Bool = true) + return PiecewiseLinearStretching(tuple(src_knots...), tuple(dst_knots...); saturate = saturate) +end + +function PiecewiseLinearStretching(src_knots::Tuple{Vararg{Union{Real,AbstractGray}}}, dst_knots::Tuple{Vararg{Union{Real,AbstractGray}}}; saturate::Bool = true) + if length(src_knots) != length(dst_knots) + throw(ArgumentError("You number of src and destination knots must match.")) + end + + if length(src_knots) < 2 || length(dst_knots) < 2 + throw(ArgumentError("You need to specify at least two knots (the start and end of an interval).")) + end + + # Promote the src_knots and dst_knots to a common type to address the case + # where someone might use a mixture of integers and floating point numbers + # when specifying the knots. + T₁ = promote_type(typeof.(src_knots)...) + T₂ = promote_type(typeof.(dst_knots)...) + T = promote_type(T₁, T₂) + + if saturate + intervals = convert_knots_to_intervals_and_saturate(T.(src_knots), T.(dst_knots)) + else + intervals = convert_knots_to_intervals(T.(src_knots), T.(dst_knots)) + end + + return PiecewiseLinearStretching(tuple(intervals...)) +end + +function convert_knots_to_intervals(src_knots::NTuple{N, <:Union{Real,AbstractGray}}, dst_knots::NTuple{N, <:Union{Real,AbstractGray}}) where {N} + intervals = Vector{Pair{<:AbstractInterval, <: AbstractInterval}}(undef, N-1) + for n = 1:(N-1) + iₛ = src_knots[n] + jₛ = dst_knots[n] + iₑ = src_knots[n+1] + jₑ = dst_knots[n+1] + if n == 1 + interval = ClosedInterval(iₛ, iₑ) => ClosedInterval(jₛ, jₑ) + intervals[n] = interval + else + interval = Interval{:open, :closed}(iₛ, iₑ) => Interval{:open, :closed}(jₛ, jₑ) + intervals[n] = interval + end + end + return intervals +end + +function convert_knots_to_intervals_and_saturate(src_knots::NTuple{N, <:Union{Real,AbstractGray}}, dst_knots::NTuple{N, <:Union{Real,AbstractGray}}) where {N} + intervals = convert_knots_to_intervals(src_knots, dst_knots) + lb = gray(typemin(eltype(src_knots))) + ub = gray(typemax(eltype(src_knots))) + lb = lb == -Inf ? lb = nextfloat(-Inf) : lb + ub = ub == Inf ? ub = prevfloat(Inf) : ub + + # Values that fall outside the first and last source knot get clamped to the first and + # last destination knot, respectively. + interval₁ = Interval{:closed, :open}(lb, first(src_knots)) => ClosedInterval(first(dst_knots), first(dst_knots)) + interval₂ = Interval{:open, :closed}(last(src_knots), ub) => ClosedInterval(last(dst_knots), last(dst_knots)) + push!(intervals, interval₁) + push!(intervals, interval₂) + return intervals +end + function (f::PiecewiseLinearStretching)(out::GenericGrayImage, img::GenericGrayImage) - # Verify that the specified dst_intervals fall inside the permissible range + # Verify that the specified destination intervals fall inside the permissible range # of the output type. T = eltype(out) - for (a,b) in f.dst_intervals + for (_, dst) in f.intervals + (a,b) = endpoints(dst) p = typemin(T) r = typemax(T) if !(p <= a <= r) || !(p <= b <= r) @@ -129,12 +230,12 @@ function (f::PiecewiseLinearStretching)(out::GenericGrayImage, img::GenericGrayI # the formula f(x) = (x-A) * ((b-a)/(B-A)) + a. The operation r * x - o is equivalent to # (x-A) * ((b-a)/(B-A)) + a. Precalculate the r and o so that later on an inner loop # contains only multiplication and addition to get better performance. - N = length(f.src_intervals) + N = length(f.intervals) rs = zeros(Float64, N) os = zeros(Float64, N) - for (src, dst) in zip(pairs(f.src_intervals), pairs(f.dst_intervals)) - n, (A,B) = src - _, (a,b) = dst + for (n, (src, dst)) in enumerate(f.intervals) + (A,B) = endpoints(src) + (a,b) = endpoints(dst) rs[n] = (b - a) / (B - A) os[n] = (A*b - B*a) / (B - A) end @@ -147,8 +248,8 @@ function (f::PiecewiseLinearStretching)(out::GenericGrayImage, img::GenericGrayI out[p] = val else is_inside_src_intervals = false - for (n, (A,B)) in pairs(f.src_intervals) - if (A <= val < B) || (A <= val <= B && n == N) + for (n, (src, _)) in enumerate(f.intervals) + if val ∈ src newval = rs[n] * val - os[n] out[p] = T <: Integer ? round(Int, newval) : newval is_inside_src_intervals = true @@ -156,14 +257,8 @@ function (f::PiecewiseLinearStretching)(out::GenericGrayImage, img::GenericGrayI end end if !is_inside_src_intervals - # Saturate pixels that fall outside the specified src_intervals by setting - # them equal to the start/end of the edge dst_intervals. - a = first(first(f.dst_intervals)) - b = last(last(f.dst_intervals)) - A = first(first(f.src_intervals)) - newval = val < A ? a : b - out[p] = T <: Integer ? round(Int, newval) : newval - end + out[p] = val + end end end return out @@ -175,5 +270,4 @@ function (f::PiecewiseLinearStretching)(out::AbstractArray{<:Color3}, img::Abstr yiq_view = channelview(yiq) adjust_histogram!(view(yiq_view,1,:,:), f) out .= convert.(T, yiq) -end - +end diff --git a/test/piecewise_linear_stretching.jl b/test/piecewise_linear_stretching.jl index ba9aa14..9ef7cb6 100644 --- a/test/piecewise_linear_stretching.jl +++ b/test/piecewise_linear_stretching.jl @@ -1,30 +1,23 @@ @testset "Piecewise Linear Stretching" begin @testset "constructors" begin - @test_throws ArgumentError PiecewiseLinearStretching( - src_intervals = ((0.0 => 0.2),(0.2 => 0.2), (0.8 => 1.0)), - dst_intervals = ((0.0 => 0.2),(-1.0 => 0.0),(1 => -2))) - - @test_throws ArgumentError PiecewiseLinearStretching( - src_intervals = ((0.0 => 0.2),(0.2 => 0.8), (0.8 => 1.0)), - dst_intervals = ((0.0 => 0.2), (1 => -2))) + @test_throws ArgumentError PiecewiseLinearStretching((0..0.2 => 0..0.4, 0.2..0.2 => 0.4..1.0)) + @test_throws ArgumentError PiecewiseLinearStretching((0, 0.2, 0.2, 1), (0.0, 0, 0, 1)) + @test_throws ArgumentError PiecewiseLinearStretching(src_knots = (0, 0.2, 1), dst_knots = (0.0, 0, 0, 1)) end + @testset "miscellaneous" begin # We should be able to invert the image. - f = PiecewiseLinearStretching(src_intervals = (0.0 => 1.0,), dst_intervals = (1.0 => 0.0,)) - for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) - img = T.(testimage("mandril_gray")) - ret = adjust_histogram(img, f) - @test ret ≈ (1 .- img) - end - + f = PiecewiseLinearStretching((0, 1), (1, 0)) + for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) + img = T.(testimage("mandril_gray")) + ret = adjust_histogram(img, f) + @test ret ≈ (1 .- img) + end # We should be able to binarize the image. - f = PiecewiseLinearStretching( - src_intervals = ((0.0 => 0.5), (0.5 => 1.0)), - dst_intervals = ((0.0 => 0.0),(1.0 => 1.0))) + f = PiecewiseLinearStretching((0.0..0.5 => 0..0, 0.5..1.0 => 1..1)) for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) img = T.(testimage("mandril_gray")) - ret = adjust_histogram(img, f) vals = sort(unique(ret)) @test first(vals) == 0.0 @@ -34,9 +27,7 @@ # Setting the src_intervals equal to the dst_intervals should result in an identity # mapping (no change to the image) provided that the src_intervals span the entire # permissible contrast range, i.e. the unit interval. - f = PiecewiseLinearStretching( - src_intervals = ((0.0 => 0.3), (0.3 => 0.8), (0.8 => 1.0)), - dst_intervals = ((0.0 => 0.3), (0.3 => 0.8), (0.8 => 1.0))) + f = PiecewiseLinearStretching((0, 0.3, 0.8 , 1.0), (0, 0.3, 0.8, 1.0)) for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) img = T.(testimage("mandril_gray")) ret = adjust_histogram(img, f) @@ -46,14 +37,25 @@ # If the src_intervals equal the dst_intervals, but the src_intervals don't span the # entire permissible contrast range, i.e. the unit interval, then the resulting # image will not match the original image because pixels that fall outside the - # specified src_intervals are satuared to the edge value of src_intervals. - f = PiecewiseLinearStretching( - src_intervals = ((0.3 => 0.8),), - dst_intervals = ((0.3 => 0.8),)) + # specified src_intervals are saturated to the edge value of src_intervals if + # saturate = true. + f = PiecewiseLinearStretching([0.3, 0.8], [0.3, 0.8]; saturate = true ) + for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) + img = T.(testimage("mandril_gray")) + ret = adjust_histogram(img, f) + @test !isapprox(ret, img, atol=1e-1) + end + + # If the src_intervals equal the dst_intervals, but the src_intervals don't span the + # entire permissible contrast range, i.e. the unit interval, then the resulting + # image will match the original image because pixels that fall outside the + # specified src_intervals are not saturated to the edge value of src_intervals if + # saturate = false + f = PiecewiseLinearStretching([0.3, 0.8], [0.3, 0.8]; saturate = false ) for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) img = T.(testimage("mandril_gray")) ret = adjust_histogram(img, f) - @test !isapprox(ret, 1 .- img, atol=1e-1) + @test ret ≈ img end for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) @@ -65,7 +67,7 @@ img = T.(collect(reshape(1/100:1/100:1, 10, 10))) minval = minimum(img) maxval = maximum(img) - f = PiecewiseLinearStretching(src_intervals = (minval => maxval,), dst_intervals = (minval => maxval,)) + f = PiecewiseLinearStretching((minval..maxval => minval..maxval,)) ret = adjust_histogram(img, f) if T <: Gray{Float32} || T <: Gray{Float64} @test all(map((i, r) -> isapprox(i, r), img, ret)) @@ -86,7 +88,7 @@ img = T.(collect(reshape(1/100:1/100:1, 10, 10))) minval = minimum(img) maxval = maximum(img) - f = PiecewiseLinearStretching(src_intervals = (minval => maxval,), dst_intervals = (0.0 => 1.0,)) + f = PiecewiseLinearStretching((minval..maxval => 0.0..1.0,)) ret = adjust_histogram(img, f) @test isapprox(0, first(ret)) @test isapprox(1, last(ret)) @@ -97,13 +99,13 @@ img = T.(testimage("mandril_gray")) minval = minimum(img) maxval = maximum(img) - f = PiecewiseLinearStretching(src_intervals = (minval => maxval,), dst_intervals = (0.0 => 1.0,)) + f = PiecewiseLinearStretching((minval..maxval => 0.0..1.0,)) ret = adjust_histogram(img, f) @test eltype(ret) == eltype(img) @test isapprox(0, minimum(ret)) @test isapprox(1, maximum(ret)) - f = PiecewiseLinearStretching(src_intervals = (0.0 => 1.0,), dst_intervals = (0.2 => 0.8,)) + f = PiecewiseLinearStretching((0.0..1.0 => 0.2..0.8,)) ret = adjust_histogram(img, f) @test eltype(ret) == eltype(img) # We are mapping the interval [0, 1] to [0.2, 0.8] but since the @@ -114,8 +116,7 @@ @test maximum(ret) <= 0.8 # Verify that results are correctly clamped to [0.2, 0.9] if it exceeds the range - f = PiecewiseLinearStretching(src_intervals = (0.3 => 0.8,), - dst_intervals = (0.2 => 0.9,)) + f = PiecewiseLinearStretching((0.3, 0.8), (0.2, 0.9)) ret = adjust_histogram(img, f) @test eltype(ret) == eltype(img) @test minimum(ret) == T(0.2) @@ -124,17 +125,15 @@ img = Float32.(collect(reshape(1/100:1/100:1, 10, 10))) - f = PiecewiseLinearStretching(src_intervals = ((0.0 => 1.0),), dst_intervals = ((-1.0 => 0.0),)) + f = PiecewiseLinearStretching((0.0,1.0), (-1.0, 0.0)) ret = adjust_histogram(img, f) @test eltype(ret) == eltype(img) - # @test minimum(ret) ≈ -1.0 - # @test maximum(ret) ≈ 0.0 @test isapprox(minimum(ret), -1.0, atol=1e-1) @test isapprox(maximum(ret), 0.0, atol=1e-1) @test_throws DomainError adjust_histogram(Gray{N0f8}.(img), f) img = Gray{Float32}.(collect(reshape(1/100:1/100:1, 10, 10))) - f = PiecewiseLinearStretching(src_intervals = ((0.0 => 0.2), (0.2 => 0.8), (0.8 => 1.0)), dst_intervals = ((0.0 => 0.2), (-1.0 => 0.0), (1 => -2))) + f = PiecewiseLinearStretching((0.0..0.2 => 0.0..0.2, 0.2..0.8 => -1.0..0.0, 0.8..1.0 => 1.0..(-2.0))) ret = adjust_histogram(img, f) @test eltype(ret) == eltype(img) @test minimum(ret) ≈ -2.0 @@ -157,8 +156,7 @@ verify that all 32 bins have non-zero counts. This will confirm that the dynamic range of the original image has been increased. =# - f = PiecewiseLinearStretching(src_intervals = (0.3 => 0.5,), - dst_intervals = (0.0 => 1.0,)) + f = PiecewiseLinearStretching((0.3, 0.5,), (0.0, 1.0)) ret = adjust_histogram(img, f) edges, counts_after = build_histogram(ret, 32, minval = 0, maxval = 1) nonzero_after = sum(counts_after .!= 0) From 02669a9ab8b33bc7e59b7ed62b5919421b143471 Mon Sep 17 00:00:00 2001 From: Zygmunt Szpak Date: Tue, 5 Oct 2021 22:41:29 +1030 Subject: [PATCH 3/4] Adds Percentiles and MinMax dispatch options --- Project.toml | 5 +- docs/src/reference.md | 6 ++ src/ImageContrastAdjustment.jl | 20 +++++ src/algorithms/piecewise_linear_stretching.jl | 90 +++++++++++++++++-- test/piecewise_linear_stretching.jl | 28 ++++++ 5 files changed, 139 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index 1c36adb..ed85c0e 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795" IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] ImageCore = "0.9.3" @@ -16,14 +17,16 @@ ImageTransformations = "0.8.1, 0.9" IntervalSets = "0.5.3" Parameters = "0.12" Reexport = "1.2.2" +StatsBase = "0.33" julia = "1" [extras] ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" +ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" [targets] -test = ["Test", "ImageFiltering", "TestImages", "ImageMagick", "LinearAlgebra"] +test = ["Test", "ImageFiltering", "TestImages", "ImageMagick", "ImageIO", "LinearAlgebra"] diff --git a/docs/src/reference.md b/docs/src/reference.md index a0eb2d2..81772df 100644 --- a/docs/src/reference.md +++ b/docs/src/reference.md @@ -13,6 +13,12 @@ adjust_histogram! build_histogram ``` +## Adjunt Types +```@docs +Percentiles +MinMax +``` + ## Algorithms ```@docs diff --git a/src/ImageContrastAdjustment.jl b/src/ImageContrastAdjustment.jl index 0c0365b..8b1a267 100644 --- a/src/ImageContrastAdjustment.jl +++ b/src/ImageContrastAdjustment.jl @@ -7,6 +7,7 @@ using ImageCore.MappedArrays using IntervalSets using Parameters: @with_kw # Same as Base.@kwdef but works on Julia 1.0 using Reexport +using StatsBase: percentile @reexport using IntervalSets # TODO: port HistogramAdjustmentAPI to ImagesAPI @@ -17,6 +18,23 @@ import .HistogramAdjustmentAPI: AbstractHistogramAdjustmentAlgorithm, # TODO Relax this to all image color types const GenericGrayImage = AbstractArray{<:Union{Number, AbstractGray}} +# Temporary definition which will be moved to ImageCore. +""" + Percentiles(p₁, p₂, p₃, ..., pₙ) +Specifies a length-n list of [percentiles](https://en.wikipedia.org/wiki/Percentile). +""" +struct Percentiles{T} <: Real + p::T +end +Percentiles(args...) = Percentiles(tuple(args...)) +""" + MinMax() +A type that can be used to signal to [`PiecewiseLinearStretching`](@ref) that it should +use the mininum and maximum value in an image as it source or destination interval. +""" +struct MinMax end + + include("core.jl") include("build_histogram.jl") include("algorithms/common.jl") @@ -40,6 +58,8 @@ export LinearStretching, ContrastStretching, PiecewiseLinearStretching, + Percentiles, + MinMax, build_histogram, adjust_histogram, adjust_histogram! diff --git a/src/algorithms/piecewise_linear_stretching.jl b/src/algorithms/piecewise_linear_stretching.jl index 1fdac6c..76752a7 100644 --- a/src/algorithms/piecewise_linear_stretching.jl +++ b/src/algorithms/piecewise_linear_stretching.jl @@ -8,6 +8,9 @@ PiecewiseLinearStretching(ClosedInterval(0.0, 0.2) => ClosedInterval(0.0, 0.1), Interval{:open, :closed}(0.2, 0.8) => Interval{:open, :closed}(0.1, 0.7), Interval{:open, :closed}(0.8, 1,0) => Interval{:open, :closed}(0.7, 1.0)) + PiecewiseLinearStretching(Percentiles(1,99), (0,1), img::GenericGrayImage) + PiecewiseLinearStretching(Percentiles(1,50,99), (0, 0.5, 1), img::GenericGrayImage) + PiecewiseLinearStretching(MinMax(), (0,1), img::GenericGrayImage) adjust_histogram([T,] img, f::PiecewiseLinearStretching) adjust_histogram!([out,] img, f::PiecewiseLinearStretching) ``` @@ -65,7 +68,7 @@ While this manner of specifying intervals permits the most flexibility, it is cumbersome for many types of mappings one would like to specify in practice. Hence, one can instead specify the mappings by designating a sequence of knot pairs. -## Specifying intervals implicitly +## Specifying intervals implicitly as knots As an alternative to specifying the endpoints of each interval, one can list the knots of the intervals. For instance, the mappings [0.0, 0.2] => [0, 0.1], (0.2, 0.8] => (0.1, 0.7] and (0.8, 1.0] => (0.7, 1.0] can be listed as follows: @@ -77,6 +80,28 @@ PiecewiseLinearStretching((0.0, 0.2, 0.8, 1.0), (0.0, 0.1, 0.7, 1.0)) PiecewiseLinearStretching(src_knots = (0.0, 0.2, 0.8, 1.0), dst_knots = (0.0, 0.1, 0.7, 1.0)) ``` +## Specifying intervals implicitly as [`Percentiles`](@ref). +Instead of directly listing knots, one can set them implicitly as percentiles of the +image pixels. For example, to may intensities between the 1st and 99th percentile to the +unit interval one could use the following code: +```julia +img = testimage("mandril_gray") +# Note that you have to pass in the img as an argument to obtain the percentiles. +f = PiecewiseLinearStretching(Percentile((1,99)), (0,1), img) +img_adjusted = adjust_histogram(img, f) +``` + +## Specifing intervals implicitly as [`MinMax`](@ref). +One can also specify and interval implicitly as the minimum and maximum intensity in an +image. For example, one can map intensities within the minimum/maximum intensity in an image +to the unit interval: +```julia +img = testimage("mandril_gray") +# Note that you have to pass in the img as an argument to obtain the minimum and maximum intensities. +f = PiecewiseLinearStretching(MinMax(), (0,1), img) +img_adjusted = adjust_histogram(img, f) +``` + ## Constraints on intervals and knots. To specify a single interval you must include a trailing comma, e.g. ```julia @@ -136,19 +161,66 @@ struct PiecewiseLinearStretching{T} <: AbstractHistogramAdjustmentAlgorithm wher end end -function PiecewiseLinearStretching(img::GenericGrayImage) - T = eltype(img) +function PiecewiseLinearStretching(; src_knots, dst_knots, saturate::Bool = true) + return PiecewiseLinearStretching(src_knots, dst_knots; saturate = saturate) +end + +function PiecewiseLinearStretching(img::GenericGrayImage; src_knots, dst_knots, saturate::Bool = true) + return PiecewiseLinearStretching(src_knots, dst_knots, img; saturate = saturate) +end + +function PiecewiseLinearStretching(src_knots::MinMax, dst_knots::T, img::GenericGrayImage; saturate::Bool = true) where T <: Union{Tuple{Vararg{Union{Real,AbstractGray}}}, AbstractVector, Percentiles} img_min, img_max = minfinite(img), maxfinite(img) if (img_min ≈ img_max) throw(DomainError("The image consists of a single intensity, which gives rise to a degenerate source interval.")) end - src_knots = (img_min, img_max) - dst_knots - (zero(T), one(T)) - return PiecewiseLinearStretching(src_knots, dst_knots; saturate = true) + if T <: Percentiles + return PiecewiseLinearStretching((img_min, img_max), dst_knots, img; saturate = saturate) + else + return PiecewiseLinearStretching((img_min, img_max), dst_knots; saturate = saturate) + end end -function PiecewiseLinearStretching(; src_knots, dst_knots, saturate::Bool = true) - PiecewiseLinearStretching(src_knots, dst_knots; saturate = saturate) +function PiecewiseLinearStretching(src_knots::T, dst_knots::MinMax, img::GenericGrayImage; saturate::Bool = true) where T <: Union{Tuple{Vararg{Union{Real,AbstractGray}}}, AbstractVector, Percentiles} + img_min, img_max = minfinite(img), maxfinite(img) + if T <: Percentiles + return PiecewiseLinearStretching(src_knots, (img_min, img_max), img ; saturate = saturate) + else + return PiecewiseLinearStretching(src_knots, (img_min, img_max); saturate = saturate) + end +end + +function PiecewiseLinearStretching(src_knots::AbstractVector, dst_knots::Tuple{Vararg{Union{Real,AbstractGray}}}; saturate::Bool = true) + return PiecewiseLinearStretching(tuple(src_knots...), dst_knots; satureate = saturate) +end + +function PiecewiseLinearStretching(src_knots::Tuple{Vararg{Union{Real,AbstractGray}}}, dst_knots::AbstractVector; saturate::Bool = true) + return PiecewiseLinearStretching(src_knots, tuple(dst_knots...); saturate = saturate) +end + +function PiecewiseLinearStretching(src_percentile::Percentiles, dst_percentile::Percentiles, img::GenericGrayImage; saturate::Bool = true) + if length(src_percentile.p) != length(dst_percentile.p) + throw(ArgumentError("The number of src percentiles and dst percentiles must match.")) + end + src_knots = percentile(vec(img), collect(src_percentile.p)) + dst_knots = percentile(vec(img), collect(dst_percentile.p)) + return PiecewiseLinearStretching(tuple(src_knots...), tuple(dst_knots...); saturate = saturate) +end + +function PiecewiseLinearStretching(src_percentile::Percentiles, dst_knots::T, img::GenericGrayImage; saturate::Bool = true) where T <: Union{Tuple{Vararg{Union{Real,AbstractGray}}}, AbstractVector} + if length(src_percentile.p) != length(dst_knots) + throw(ArgumentError("The number of src percentiles and destination knots must match.")) + end + src_knots = percentile(vec(img), collect(src_percentile.p)) + return PiecewiseLinearStretching(tuple(src_knots...), dst_knots; saturate = saturate) +end + +function PiecewiseLinearStretching(src_knots::T, dst_percentile::Percentiles, img::GenericGrayImage; saturate::Bool = true) where T <: Union{Tuple{Vararg{Union{Real,AbstractGray}}}, AbstractVector} + if length(src_knots) != length(dst_percentile.p) + throw(ArgumentError("The number of src knots and destination percentiles must match.")) + end + dst_knots = percentile(vec(img), collect(dst_percentile.p)) + return PiecewiseLinearStretching(src_knots, tuple(dst_knots...); saturate = saturate) end function PiecewiseLinearStretching(src_knots::AbstractVector, dst_knots::AbstractVector; saturate::Bool = true) @@ -157,7 +229,7 @@ end function PiecewiseLinearStretching(src_knots::Tuple{Vararg{Union{Real,AbstractGray}}}, dst_knots::Tuple{Vararg{Union{Real,AbstractGray}}}; saturate::Bool = true) if length(src_knots) != length(dst_knots) - throw(ArgumentError("You number of src and destination knots must match.")) + throw(ArgumentError("The number of src and destination knots must match.")) end if length(src_knots) < 2 || length(dst_knots) < 2 diff --git a/test/piecewise_linear_stretching.jl b/test/piecewise_linear_stretching.jl index 9ef7cb6..e02cec8 100644 --- a/test/piecewise_linear_stretching.jl +++ b/test/piecewise_linear_stretching.jl @@ -1,9 +1,12 @@ @testset "Piecewise Linear Stretching" begin @testset "constructors" begin + img = testimage("mandril_gray") @test_throws ArgumentError PiecewiseLinearStretching((0..0.2 => 0..0.4, 0.2..0.2 => 0.4..1.0)) @test_throws ArgumentError PiecewiseLinearStretching((0, 0.2, 0.2, 1), (0.0, 0, 0, 1)) @test_throws ArgumentError PiecewiseLinearStretching(src_knots = (0, 0.2, 1), dst_knots = (0.0, 0, 0, 1)) + @test_throws ArgumentError PiecewiseLinearStretching(src_knots = (0,), dst_knots = (1,)) + @test_throws ArgumentError PiecewiseLinearStretching(Percentiles(10,50, 100,150), (0, 0.5, 1), img) end @testset "miscellaneous" begin @@ -123,6 +126,31 @@ @test maximum(ret) == T(0.9) end + # Verify that Percentiles and MinMax works correctly + for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) + img = T.(collect(reshape(1/100:1/100:1, 10, 10))) + f = PiecewiseLinearStretching(Percentiles((10,90)), MinMax(), img) + ret = adjust_histogram(img, f) + @test all(isapprox.(ret[:,1], 0.01, atol = 1e-2)) + @test all(isapprox.(ret[:,end], 1.0, atol = 1e-2)) + + # Same as above, but with keyword arguments. + f = PiecewiseLinearStretching(img; src_knots = Percentiles((10,90)), dst_knots = MinMax()) + ret = adjust_histogram(img, f) + @test all(isapprox.(ret[:,1], 0.01, atol=1e-2)) + @test all(isapprox.(ret[:,end], 1.0, atol=1e-2)) + + f = PiecewiseLinearStretching(MinMax(), Percentiles((10,90)), img) + ret = adjust_histogram(img, f) + @test isapprox(ret[1,1], 0.109, atol=1e-2) + @test isapprox(ret[end,end], 0.901, atol=1e-2) + + f = PiecewiseLinearStretching(Percentiles(10,90), Percentiles((10,90)), img) + ret = adjust_histogram(img, f) + @test all(isapprox.(ret[:,1], 0.109, atol=1e-2)) + @test all(isapprox.(ret[:,end], 0.901, atol=1e-2)) + end + img = Float32.(collect(reshape(1/100:1/100:1, 10, 10))) f = PiecewiseLinearStretching((0.0,1.0), (-1.0, 0.0)) From 87bee8a8a2976f43fb6ac971a2128d4ed0c5676c Mon Sep 17 00:00:00 2001 From: Zygmunt Szpak Date: Sat, 9 Oct 2021 22:40:23 +1030 Subject: [PATCH 4/4] Expands test coverage --- src/algorithms/piecewise_linear_stretching.jl | 20 ++++--- test/piecewise_linear_stretching.jl | 54 ++++++++++++++++--- 2 files changed, 61 insertions(+), 13 deletions(-) diff --git a/src/algorithms/piecewise_linear_stretching.jl b/src/algorithms/piecewise_linear_stretching.jl index 76752a7..5c6a42f 100644 --- a/src/algorithms/piecewise_linear_stretching.jl +++ b/src/algorithms/piecewise_linear_stretching.jl @@ -191,7 +191,7 @@ function PiecewiseLinearStretching(src_knots::T, dst_knots::MinMax, img::Generic end function PiecewiseLinearStretching(src_knots::AbstractVector, dst_knots::Tuple{Vararg{Union{Real,AbstractGray}}}; saturate::Bool = true) - return PiecewiseLinearStretching(tuple(src_knots...), dst_knots; satureate = saturate) + return PiecewiseLinearStretching(tuple(src_knots...), dst_knots; saturate = saturate) end function PiecewiseLinearStretching(src_knots::Tuple{Vararg{Union{Real,AbstractGray}}}, dst_knots::AbstractVector; saturate::Bool = true) @@ -279,10 +279,16 @@ function convert_knots_to_intervals_and_saturate(src_knots::NTuple{N, <:Union{Re # Values that fall outside the first and last source knot get clamped to the first and # last destination knot, respectively. - interval₁ = Interval{:closed, :open}(lb, first(src_knots)) => ClosedInterval(first(dst_knots), first(dst_knots)) - interval₂ = Interval{:open, :closed}(last(src_knots), ub) => ClosedInterval(last(dst_knots), last(dst_knots)) - push!(intervals, interval₁) - push!(intervals, interval₂) + if lb != first(src_knots) + interval₁ = Interval{:closed, :open}(lb, first(src_knots)) => ClosedInterval(first(dst_knots), first(dst_knots)) + push!(intervals, interval₁) + end + + if ub != last(src_knots) + interval₂ = Interval{:open, :closed}(last(src_knots), ub) => ClosedInterval(last(dst_knots), last(dst_knots)) + push!(intervals, interval₂) + end + return intervals end @@ -306,8 +312,8 @@ function (f::PiecewiseLinearStretching)(out::GenericGrayImage, img::GenericGrayI rs = zeros(Float64, N) os = zeros(Float64, N) for (n, (src, dst)) in enumerate(f.intervals) - (A,B) = endpoints(src) - (a,b) = endpoints(dst) + (A,B) = Float64.(endpoints(src)) + (a,b) = Float64.(endpoints(dst)) rs[n] = (b - a) / (B - A) os[n] = (A*b - B*a) / (B - A) end diff --git a/test/piecewise_linear_stretching.jl b/test/piecewise_linear_stretching.jl index e02cec8..b8913ef 100644 --- a/test/piecewise_linear_stretching.jl +++ b/test/piecewise_linear_stretching.jl @@ -6,10 +6,16 @@ @test_throws ArgumentError PiecewiseLinearStretching((0, 0.2, 0.2, 1), (0.0, 0, 0, 1)) @test_throws ArgumentError PiecewiseLinearStretching(src_knots = (0, 0.2, 1), dst_knots = (0.0, 0, 0, 1)) @test_throws ArgumentError PiecewiseLinearStretching(src_knots = (0,), dst_knots = (1,)) - @test_throws ArgumentError PiecewiseLinearStretching(Percentiles(10,50, 100,150), (0, 0.5, 1), img) + @test_throws ArgumentError PiecewiseLinearStretching(Percentiles(10, 50, 100, 150), (0, 0.5, 1), img) + @test_throws ArgumentError PiecewiseLinearStretching((0, 1), Percentiles(0, 0.5, 1), img) + @test_throws ArgumentError PiecewiseLinearStretching(Percentiles(10, 50, 100, 150), Percentiles(0, 100), img) end @testset "miscellaneous" begin + # An image that consists of a single intensity gives rise to a degenerate source + # interval. + @test_throws DomainError PiecewiseLinearStretching(MinMax(), (0, 1), ones(10,10)) + # We should be able to invert the image. f = PiecewiseLinearStretching((0, 1), (1, 0)) for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) @@ -29,8 +35,8 @@ # Setting the src_intervals equal to the dst_intervals should result in an identity # mapping (no change to the image) provided that the src_intervals span the entire - # permissible contrast range, i.e. the unit interval. - f = PiecewiseLinearStretching((0, 0.3, 0.8 , 1.0), (0, 0.3, 0.8, 1.0)) + # permissible contrast range, i.e. the unit interval. + f = PiecewiseLinearStretching((0, 0.3, 0.8 , 1.0), (0, 0.3, 0.8, 1.0)) for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) img = T.(testimage("mandril_gray")) ret = adjust_histogram(img, f) @@ -53,7 +59,7 @@ # entire permissible contrast range, i.e. the unit interval, then the resulting # image will match the original image because pixels that fall outside the # specified src_intervals are not saturated to the edge value of src_intervals if - # saturate = false + # saturate = false. f = PiecewiseLinearStretching([0.3, 0.8], [0.3, 0.8]; saturate = false ) for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) img = T.(testimage("mandril_gray")) @@ -77,6 +83,42 @@ else @test ret == img end + #MinMax & Tuple. + f = PiecewiseLinearStretching(MinMax(), (minval, maxval), img) + ret = adjust_histogram(img, f) + @test ret ≈ img + if T <: Gray{Float32} || T <: Gray{Float64} + @test all(map((i, r) -> isapprox(i, r), img, ret)) + else + @test ret == img + end + # Tuple & MinMax. + f = PiecewiseLinearStretching((minval, maxval), MinMax(), img) + ret = adjust_histogram(img, f) + @test ret ≈ img + if T <: Gray{Float32} || T <: Gray{Float64} + @test all(map((i, r) -> isapprox(i, r), img, ret)) + else + @test ret == img + end + # Tuple & Vector. + f = PiecewiseLinearStretching((minval, maxval), [minval, maxval]) + ret = adjust_histogram(img, f) + @test ret ≈ img + if T <: Gray{Float32} || T <: Gray{Float64} + @test all(map((i, r) -> isapprox(i, r), img, ret)) + else + @test ret == img + end + # Vector & Tuple. + f = PiecewiseLinearStretching([minval, maxval], (minval, maxval)) + ret = adjust_histogram(img, f) + @test ret ≈ img + if T <: Gray{Float32} || T <: Gray{Float64} + @test all(map((i, r) -> isapprox(i, r), img, ret)) + else + @test ret == img + end # Verify that NaN is also handled correctly. if T <: Gray{Float32} || T <: Gray{Float64} @@ -118,7 +160,7 @@ @test minimum(ret) >= 0.2 @test maximum(ret) <= 0.8 - # Verify that results are correctly clamped to [0.2, 0.9] if it exceeds the range + # Verify that results are correctly clamped to [0.2, 0.9] if it exceeds the range. f = PiecewiseLinearStretching((0.3, 0.8), (0.2, 0.9)) ret = adjust_histogram(img, f) @test eltype(ret) == eltype(img) @@ -126,7 +168,7 @@ @test maximum(ret) == T(0.9) end - # Verify that Percentiles and MinMax works correctly + # Verify that Percentiles and MinMax works correctly. for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) img = T.(collect(reshape(1/100:1/100:1, 10, 10))) f = PiecewiseLinearStretching(Percentiles((10,90)), MinMax(), img)