diff --git a/Project.toml b/Project.toml index 6e62cd0..ed85c0e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,25 +1,32 @@ 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" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] ImageCore = "0.9.3" 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/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..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 @@ -48,3 +54,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..8b1a267 100644 --- a/src/ImageContrastAdjustment.jl +++ b/src/ImageContrastAdjustment.jl @@ -4,7 +4,11 @@ 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 +using StatsBase: percentile +@reexport using IntervalSets # TODO: port HistogramAdjustmentAPI to ImagesAPI include("HistogramAdjustmentAPI/HistogramAdjustmentAPI.jl") @@ -14,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") @@ -24,6 +45,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 +57,9 @@ export GammaCorrection, 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 new file mode 100644 index 0000000..5c6a42f --- /dev/null +++ b/src/algorithms/piecewise_linear_stretching.jl @@ -0,0 +1,351 @@ +""" +``` + PiecewiseLinearStretching <: AbstractHistogramAdjustmentAlgorithm + 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)) + 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) +``` + +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 + +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. + +## 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. + +## 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: + +```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)) +``` + +## 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 +# 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 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`. + +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. + +When listing knots, the length of `src_knots` must match the length of `dst_knots`. + + +## Saturating pixels + +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 + +```julia +using ImageContrastAdjustment, TestImages + +img = testimage("mandril_gray") +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) +""" +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 + new{T₁}(intervals) + end +end + +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 + 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::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; saturate = 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) + 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("The 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. + 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 + +function (f::PiecewiseLinearStretching)(out::GenericGrayImage, img::GenericGrayImage) + # Verify that the specified destination intervals fall inside the permissible range + # of the output type. + T = eltype(out) + for (_, dst) in f.intervals + (a,b) = endpoints(dst) + 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.intervals) + rs = zeros(Float64, N) + os = zeros(Float64, N) + for (n, (src, dst)) in enumerate(f.intervals) + (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 + + # 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, (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 + break + end + end + if !is_inside_src_intervals + out[p] = val + 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..b8913ef --- /dev/null +++ b/test/piecewise_linear_stretching.jl @@ -0,0 +1,238 @@ +@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) + @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}) + 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((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 + @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((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) + @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 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 ret ≈ img + 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((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)) + 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} + 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((minval..maxval => 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((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((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 + # 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((0.3, 0.8), (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 + + # 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)) + ret = adjust_histogram(img, f) + @test eltype(ret) == eltype(img) + @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((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 + @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((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) + @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