diff --git a/Project.toml b/Project.toml index df1d32cf..ea1c53af 100644 --- a/Project.toml +++ b/Project.toml @@ -38,7 +38,7 @@ Graphics = "0.4, 1.0" ImageAxes = "0.6.9" ImageBase = "0.1.5" ImageBinarization = "0.2, 0.3" -ImageContrastAdjustment = "0.3.3" +ImageContrastAdjustment = "0.3.12" ImageCore = "0.9.3, 0.10" ImageCorners = "0.1.2" ImageDistances = "0.2.5" @@ -61,6 +61,7 @@ TiledIteration = "0.2, 0.3, 0.4, 0.5" julia = "1.6" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -69,4 +70,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" [targets] -test = ["FFTW", "LinearAlgebra", "Random", "Suppressor", "Test", "TestImages"] +test = ["Aqua", "FFTW", "LinearAlgebra", "Random", "Suppressor", "Test", "TestImages"] diff --git a/src/Images.jl b/src/Images.jl index 93a33213..80f19b52 100644 --- a/src/Images.jl +++ b/src/Images.jl @@ -48,33 +48,17 @@ import FileIO: metadata import Graphics: Point include("compat.jl") -include("labeledarrays.jl") include("algorithms.jl") include("deprecations.jl") include("edge.jl") export - # types - ColorizedArray, - # macros @test_approx_eq_sigma_eps, # algorithms imedge, # TODO: deprecate? imgaussiannoise, - otsu_threshold, - yen_threshold, - - #Exposure - imhist, - histeq, - adjust_gamma, - histmatch, - clahe, - imadjustintensity, - imstretch, - cliphist, magnitude, magnitude_phase, @@ -85,10 +69,8 @@ export thin_edges_nonmaxsup, thin_edges_nonmaxsup_subpix, canny, - gaussian_pyramid, + gaussian_pyramid - # phantoms - shepp_logan """ Images.jl is an "umbrella package" that exports a set of packages which are useful for diff --git a/src/algorithms.jl b/src/algorithms.jl index 6de80bbb..45e2d2c6 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -120,85 +120,3 @@ function pyramid_scale(img::OffsetArray, downsample) off = (.-ceil.(Int,(.-iterate.(map(x->UnitRange(x).-1,axes(img)))[1])./downsample)) OffsetArray(imresize(img, sz_next), off) end - -""" -``` -thres = otsu_threshold(img) -thres = otsu_threshold(img, bins) -``` - -Computes threshold for grayscale image using Otsu's method. - -Parameters: -- img = Grayscale input image -- bins = Number of bins used to compute the histogram. Needed for floating-point images. - -""" -function otsu_threshold(img::AbstractArray{T, N}, bins::Int = 256) where {T<:Union{Gray,Real}, N} - - min, max = extrema(img) - edges, counts = imhist(img, range(gray(min), stop=gray(max), length=bins)) - histogram = counts./sum(counts) - - ω0 = 0 - μ0 = 0 - μt = 0 - μT = sum((1:(bins+1)).*histogram) - max_σb=0.0 - thres=1 - - for t in 1:bins - ω0 += histogram[t] - ω1 = 1 - ω0 - μt += t*histogram[t] - - σb = (μT*ω0-μt)^2/(ω0*ω1) - - if(σb > max_σb) - max_σb = σb - thres = t - end - end - - return T((edges[thres-1]+edges[thres])/2) -end - -""" -``` -thres = yen_threshold(img) -thres = yen_threshold(img, bins) -``` - -Computes threshold for grayscale image using Yen's maximum correlation criterion for bilevel thresholding - -Parameters: -- img = Grayscale input image -- bins = Number of bins used to compute the histogram. Needed for floating-point images. - - -#Citation -Yen J.C., Chang F.J., and Chang S. (1995) “A New Criterion for Automatic Multilevel Thresholding” IEEE Trans. on Image Processing, 4(3): 370-378. DOI:10.1109/83.366472 -""" -function yen_threshold(img::AbstractArray{T, N}, bins::Int = 256) where {T<:Union{Gray, Real}, N} - - min, max = extrema(img) - if(min == max) - return T(min) - end - - edges, counts = imhist(img, range(gray(min), stop=gray(max), length=bins)) - - prob_mass_function = counts./sum(counts) - clamp!(prob_mass_function,eps(),Inf) - prob_mass_function_sq = prob_mass_function.^2 - cum_pmf = cumsum(prob_mass_function) - cum_pmf_sq_1 = cumsum(prob_mass_function_sq) - cum_pmf_sq_2 = reverse!(cumsum(reverse!(prob_mass_function_sq))) - - #Equation (4) cited in the paper. - criterion = log.(((cum_pmf[1:end-1].*(1.0 .- cum_pmf[1:end-1])).^2) ./ (cum_pmf_sq_1[1:end-1].*cum_pmf_sq_2[2:end])) - - thres = edges[findmax(criterion)[2]] - return T(thres) - -end diff --git a/src/deprecations.jl b/src/deprecations.jl index b879250f..ab1afaa5 100644 --- a/src/deprecations.jl +++ b/src/deprecations.jl @@ -1,1154 +1,2 @@ -function build_histogram(img::AbstractArray, nbins::Integer, minval::Union{Real,AbstractGray}, maxval::Union{Real,AbstractGray}) - Base.depwarn("`build_histogram(img, nbins, minval, maxval)` is deprecated, use build_histogram(img, nbins; minval = minval, maxval = maxval)) instead.", :build_histogram) - ImageContrastAdjustment.build_histogram(img, nbins; minval = minval, maxval = maxval) -end - -function adjust_histogram(operation::Equalization, img::AbstractArray, nbins::Integer, minval::Union{Real,AbstractGray} = 0, maxval::Union{Real,AbstractGray} = 1) - Base.depwarn("`adjust_histogram(Equalization(),img)` is deprecated, use adjust_histogram(img, Equalization(; nbins, minval, maxval)) instead.", :adjust_histogram) - ImageContrastAdjustment.adjust_histogram(img, Equalization(nbins = nbins, minval = minval, maxval = maxval)) -end - -function adjust_histogram(operation::Equalization, img::AbstractArray{T}, nbins::Integer, minval::Union{Real,AbstractGray} = 0, maxval::Union{Real,AbstractGray} = 1) where {T<:Color3} - Base.depwarn("`adjust_histogram(Equalization(),img)` is deprecated, use adjust_histogram(img, Equalization(; nbins, minval, maxval)) instead.", :adjust_histogram) - ImageContrastAdjustment.adjust_histogram(img, Equalization(nbins = nbins, minval = minval, maxval = maxval)) -end - -function adjust_histogram!(operation::Equalization, img::AbstractArray, nbins::Integer, minval::Union{Real,AbstractGray} = 0, maxval::Union{Real,AbstractGray} = 1) - Base.depwarn("`adjust_histogram!(Equalization(),img)` is deprecated, use adjust_histogram!(img,img, Equalization(; nbins, minval, maxval)) instead.", :adjust_histogram!) - ImageContrastAdjustment.adjust_histogram!(img, img, Equalization(nbins = nbins, minval = minval, maxval = maxval)) -end - -function adjust_histogram(operation::Matching, img::AbstractArray, targetimg::AbstractArray, nbins::Integer = 256) - Base.depwarn("`adjust_histogram(Matching(),img, targetimg, nbins)` is deprecated, use adjust_histogram(img, Matching(; targetimg = targetimg, nbins = nbins)) instead.", :adjust_histogram) - ImageContrastAdjustment.adjust_histogram(img, Matching(targetimg = targetimg, nbins = nbins)) -end - -function adjust_histogram(operation::Matching, img::AbstractArray, targetimg::AbstractArray, edges::AbstractRange) - Base.depwarn("`adjust_histogram(Matching(),img, targetimg, edges)` is deprecated, use adjust_histogram(img, Matching(; targetimg = targetimg, edges = edges)) instead.", :adjust_histogram) - ImageContrastAdjustment.adjust_histogram(img, Matching(targetimg = targetimg, edges = edges)) -end - -function adjust_histogram!(operation::Matching, img::AbstractArray{T}, targetimg::AbstractArray{T}, edges::AbstractRange ) where T <: Color3 - Base.depwarn("`adjust_histogram!(Matching(),img, targetimg, edges)` is deprecated, use adjust_histogram!(img, Matching(; targetimg = targetimg, edges = edges)) instead.", :adjust_histogram!) - ImageContrastAdjustment.adjust_histogram!(img, img, Matching(targetimg = targetimg, edges = edges)) -end - -function adjust_histogram!(operation::Matching, img::AbstractArray{T}, targetimg::AbstractArray{T}, nbins::Integer = 256 ) where T <: Color3 - Base.depwarn("`adjust_histogram!(Matching(),img, targetimg, nbins)` is deprecated, use adjust_histogram!(img, img, Matching(; targetimg = targetimg, nbins = nbins)) instead.", :adjust_histogram!) - ImageContrastAdjustment.adjust_histogram!(img, img, Matching(targetimg = targetimg, nbins = nbins)) -end - -function adjust_histogram!(operation::Matching, img::AbstractArray, targetimg::AbstractArray, edges::AbstractRange) - Base.depwarn("`adjust_histogram!(Matching(),img, targetimg, edges)` is deprecated, use adjust_histogram!(img, Matching(; targetimg = targetimg, edges = edges)) instead.", :adjust_histogram!) - ImageContrastAdjustment.adjust_histogram!(img, img, Matching(targetimg = targetimg, edges = edges)) -end - -function adjust_histogram!(operation::Matching, img::AbstractArray, targetimg::AbstractArray, nbins::Integer = 256 ) - Base.depwarn("`adjust_histogram!(Matching(),img, targetimg, nbins)` is deprecated, use adjust_histogram!(img, Matching(; targetimg = targetimg, nbins = nbins)) instead.", :adjust_histogram!) - ImageContrastAdjustment.adjust_histogram!(img, img, Matching(targetimg = targetimg, nbins = nbins)) -end - -""" - `imadjustintensity(img [, (minval,maxval)]) -> Image` - - Map intensities over the interval `(minval,maxval)` to the interval - `[0,1]`. This is equivalent to `map(ScaleMinMax(eltype(img), minval, - maxval), img)`. (minval,maxval) defaults to `extrema(img)`. -""" -function imadjustintensity(img::AbstractArray{T}, range::Tuple{Any,Any}) where {T} - Base.depwarn("`imadjustintensity` will be removed in a future release, please use `adjust_histogram(img, LinearStretching())` instead.", :imadjustintensity) - return map(scaleminmax(T, range...), img) -end -function imadjustintensity(img::AbstractArray, range::AbstractArray) - Base.depwarn("`imadjustintensity` will be removed in a future release, please use `adjust_histogram(img, LinearStretching())` instead.", :imadjustintensity) - return imadjustintensity(img, (range...,)) -end -function imadjustintensity(img::AbstractArray) - Base.depwarn("`imadjustintensity` will be removed in a future release, please use `adjust_histogram(img, LinearStretching())` instead.", :imadjustintensity) - return map(takemap(scaleminmax, img), img) -end - -_imstretch(img::AbstractArray{T}, m::Number, slope::Number) where {T} = map(i -> 1 / (1 + (m / (i + eps(T))) ^ slope), img) - -""" -`imgs = imstretch(img, m, slope)` enhances or reduces (for -slope > 1 or < 1, respectively) the contrast near saturation (0 and 1). This is -essentially a symmetric gamma-correction. For a pixel of brightness `p`, the new -intensity is `1/(1+(m/(p+eps))^slope)`. - -This assumes the input `img` has intensities between 0 and 1. -""" -function imstretch(img::AbstractArray, m::Number, slope::Number) - Base.depwarn("`imstretch` will be removed in a future release, please use `adjust_histogram(img, ContrastStretching())` instead.", :imstretch) - return _imstretch(float(img), m, slope) -end -function imstretch(img::ImageMeta, m::Number, slope::Number) - Base.depwarn("`imstretch` will be removed in a future release, please use `adjust_histogram(img, ContrastStretching())` instead.", :imstretch) - return shareproperties(img, imstretch(arraydata(img), m, slope)) -end - -function imhist(img::AbstractArray{T}, nbins::Integer = 400) where {T<:Colorant} - Base.depwarn("`imhist` will be removed in a future release, please use `build_histogram` instead.", :imhist) - return imhist(convert(Array{Gray}, img), nbins) -end - -function imhist(img::AbstractArray{T}, nbins::Integer = 400) where T<:NumberLike - Base.depwarn("`imhist` will be removed in a future release, please use `build_histogram` instead.", :imhist) - minval = minfinite(img) - maxval = maxfinite(img) - imhist(img, nbins, minval, maxval) -end - - -""" -``` -edges, count = imhist(img, nbins) -edges, count = imhist(img, nbins, minval, maxval) -edges, count = imhist(img, edges) -``` - -Generates a histogram for the image over nbins spread between `(minval, maxval]`. -Color images are automatically converted to grayscale. - -# Output -Returns `edges` which is a [`range`](@ref) type that specifies how the interval -`(minval, maxval]` is divided into bins, and an array `count` which records the -concomitant bin frequencies. In particular, `count` has the following -properties: - -* `count[i+1]` is the number of values `x` that satisfy `edges[i] <= x < edges[i+1]`. -* `count[1]` is the number satisfying `x < edges[1]`, and -* `count[end]` is the number satisfying `x >= edges[end]`. -* `length(count) == length(edges)+1`. - -# Details - -One can consider a histogram as a piecewise-constant model of a probability -density function ``f`` [1]. Suppose that ``f`` has support on some interval ``I = -[a,b]``. Let ``m`` be an integer and ``a = a_1 < a_2 < \\ldots < a_m < a_{m+1} = -b`` a sequence of real numbers. Construct a sequence of intervals - -```math -I_1 = [a_1,a_2], I_2 = (a_2, a_3], \\ldots, I_{m} = (a_m,a_{m+1}] -``` - -which partition ``I`` into subsets ``I_j`` ``(j = 1, \\ldots, m)`` on which -``f`` is constant. These subsets satisfy ``I_i \\cap I_j = \\emptyset, \\forall -i \\neq j``, and are commonly referred to as *bins*. Together they encompass the -entire range of data values such that ``\\sum_j |I_j | = | I |``. Each bin has -width ``w_j = |I_j| = a_{j+1} - a_j`` and height ``h_j`` which is the constant -probability density over the region of the bin. Integrating the constant -probability density over the width of the bin ``w_j`` yields a probability mass -of ``\\pi_j = h_j w_j`` for the bin. - -For a sample ``x_1, x_2, \\ldots, x_N``, let - - -```math -n_j = \\sum_{n = 1}^{N}\\mathbf{1}_{(I_j)}(x_n), -\\quad \\text{where} \\quad -\\mathbf{1}_{(I_j)}(x) = -\\begin{cases} - 1 & \\text{if} x \\in I_j,\\\\ - 0 & \\text{otherwise}, -\\end{cases}, -``` -represents the number of samples falling into the interval ``I_j``. An estimate -for the probability mass of the ``j``th bin is given by the relative frequency -``\\hat{\\pi} = \\frac{n_j}{N}``, and the histogram estimator of the probability -density function is defined as -```math -\\begin{aligned} -\\hat{f}_n(x) & = \\sum_{j = 1}^{m}\\frac{n_j}{Nw_j} \\mathbf{1}_{(I_j)}(x) \\\\ -& = \\sum_{j = 1}^{m}\\frac{\\hat{\\pi}_j}{w_j} \\mathbf{1}_{(I_j)}(x) \\\\ -& = \\sum_{j = 1}^{m}\\hat{h}_j \\mathbf{1}_{(I_j)}(x). -\\end{aligned} -``` - -The function ``\\hat{f}_n(x)`` is a genuine density estimator because ``\\hat{f}_n(x) \\ge 0`` and -```math -\\begin{aligned} -\\int_{-\\infty}^{\\infty}\\hat{f}_n(x) \\operatorname{d}x & = \\sum_{j=1}^{m} \\frac{n_j}{Nw_j} w_j \\\\ -& = 1. -\\end{aligned} -``` - -# Options -Various options for the parameters of this function are described in more detail -below. - -## Choices for `nbins` -You can specify the number of discrete bins for the histogram. - -## Choices for `minval` -You have the option to specify the lower bound of the interval over which the -histogram will be computed. If `minval` is not specified then the minimum -value present in the image is taken as the lower bound. - -## Choices for `maxval` -You have the option to specify the upper bound of the interval over which the -histogram will be computed. If `maxval` is not specified then the maximum -value present in the image is taken as the upper bound. - -## Choices for `edges` -If you do not designate the number of bins, nor the lower or upper bound of the -interval, then you have the option to directly stipulate how the intervals will -be divided by specifying a [`range`](@ref) type. - -# Example - -Compute the histogram of a grayscale image. -```julia - -using TestImages, FileIO, ImageView - -img = testimage("mandril_gray"); -edges, counts = imhist(img,256); -``` - -Given a color image, compute the hisogram of the red channel. -```julia -img = testimage("mandrill") -r = red(img) -edges, counts = imhist(r,256); -``` - -# References -[1] E. Herrholz, "Parsimonious Histograms," Ph.D. dissertation, Inst. of Math. and Comp. Sci., University of Greifswald, Greifswald, Germany, 2011. -""" -function imhist(img::AbstractArray, nbins::Integer, minval::NumberLike, maxval::NumberLike) - Base.depwarn("`imhist` will be removed in a future release, please use `build_histogram` instead.", :imhist) - edges = StatsBase.histrange([Float64(minval), Float64(maxval)], nbins, :left) - imhist(img, edges) -end - -function imhist(img::AbstractArray, edges::AbstractRange) - Base.depwarn("`imhist` will be removed in a future release, please use `build_histogram` instead.", :imhist) - histogram = zeros(Int, length(edges) + 1) - o = Base.Order.Forward - G = _graytype(eltype(img)) - for v in img - val = real(convert(G, v)) - if val >= edges[end] - histogram[end] += 1 - continue - end - index = searchsortedlast(edges, val, o) - histogram[index + 1] += 1 - end - edges, histogram -end - - -function _histeq_pixel_rescale(pixel::T, cdf, minval, maxval) where T<:NumberLike - n = length(cdf) - bin_pixel = clamp(ceil(Int, gray((pixel - minval) * length(cdf) / (maxval - minval))), 1, n) - rescaled_pixel = minval + ((cdf[bin_pixel] - cdf[1]) * (maxval - minval) / (cdf[end] - cdf[1])) - convert(T, rescaled_pixel) -end -function _histeq_pixel_rescale(pixel::C, cdf, minval, maxval) where C<:Color - yiq = convert(YIQ, pixel) - y = _histeq_pixel_rescale(yiq.y, cdf, minval, maxval) - convert(C, YIQ(y, yiq.i, yiq.q)) -end -function _histeq_pixel_rescale(pixel::C, cdf, minval, maxval) where C<:TransparentColor - base_colorant_type(C)(_histeq_pixel_rescale(color(pixel), cdf, minval, maxval), alpha(pixel)) -end - -""" -``` -hist_equalised_img = histeq(img, nbins) -hist_equalised_img = histeq(img, nbins, minval, maxval) -``` - -Returns a histogram equalised image with a granularity of approximately `nbins` -number of bins. - -# Details - -Histogram equalisation was initially conceived to improve the contrast in a -single-channel grayscale image. The method transforms the -distribution of the intensities in an image so that they are as uniform as -possible [1]. The natural justification for uniformity -is that the image has better contrast if the intensity levels of an image span -a wide range on the intensity scale. As it turns out, the necessary -transformation is a mapping based on the cumulative histogram. - -One can consider an ``L``-bit single-channel ``I \\times J`` image with gray -values in the set ``\\{0,1,\\ldots,L-1 \\}``, as a collection of independent and -identically distributed random variables. Specifically, let the sample space -``\\Omega`` be the set of all ``IJ``-tuples ``\\omega -=(\\omega_{11},\\omega_{12},\\ldots,\\omega_{1J},\\omega_{21},\\omega_{22},\\ldots,\\omega_{2J},\\omega_{I1},\\omega_{I2},\\ldots,\\omega_{IJ})``, -where each ``\\omega_{ij} \\in \\{0,1,\\ldots, L-1 \\}``. Furthermore, impose a -probability measure on ``\\Omega`` such that the functions ``\\Omega \\ni -\\omega \\to \\omega_{ij} \\in \\{0,1,\\ldots,L-1\\}`` are independent and -identically distributed. - -One can then regard an image as a matrix of random variables ``\\mathbf{G} = -[G_{i,j}(\\omega)]``, where each function ``G_{i,j}: \\Omega \\to \\mathbb{R}`` -is defined by -```math -G_{i,j}(\\omega) = \\frac{\\omega_{ij}}{L-1}, -``` -and each ``G_{i,j}`` is distributed according to some unknown density ``f_{G}``. -While ``f_{G}`` is unknown, one can approximate it with a normalised histogram -of gray levels, - -```math -\\hat{f}_{G}(v)= \\frac{n_v}{IJ}, -``` -where -```math -n_v = \\left | \\left\\{(i,j)\\, |\\, G_{i,j}(\\omega) = v \\right \\} \\right | -``` -represents the number of times a gray level with intensity ``v`` occurs in -``\\mathbf{G}``. To transforming the distribution of the intensities so that -they are as uniform as possible one needs to find a mapping ``T(\\cdot)`` such -that ``T(G_{i,j}) \\thicksim U ``. The required mapping turns out to be the -cumulative distribution function (CDF) of the empirical density -``\\hat{f}_{G}``, -```math - T(G_{i,j}) = \\int_0^{G_{i,j}}\\hat{f}_{G}(w)\\mathrm{d} w. -``` - -# Options - -Various options for the parameters of this function are described in more detail -below. - -## Choices for `img` - -The `histeq` function can handle a variety of input types. The returned image -depends on the input type. If the input is an `Image` then the resulting image -is of the same type and has the same properties. - -For coloured images, the input is converted to -[YIQ](https://en.wikipedia.org/wiki/YIQ) type and the Y channel is equalised. -This is the combined with the I and Q channels and the resulting image converted -to the same type as the input. - -## Choices for `nbins` - -You can specify the total number of bins in the histogram. - -## Choices for `minval` and `maxval` - -If minval and maxval are specified then intensities are equalized to the range -(minval, maxval). The default values are 0 and 1. - -# Example - -```julia - -using TestImages, FileIO, ImageView - -img = testimage("mandril_gray"); -imgeq = histeq(img,256); - -imshow(img) -imshow(imgeq) -``` - -# References -1. R. C. Gonzalez and R. E. Woods. *Digital Image Processing (3rd Edition)*. Upper Saddle River, NJ, USA: Prentice-Hall, 2006. - -See also: [histmatch](@ref),[clahe](@ref), [imhist](@ref) and [adjust_gamma](@ref). - -""" -function histeq(img::AbstractArray, nbins::Integer, minval::NumberLike, maxval::NumberLike) - Base.depwarn("`histeq` will be removed in a future release, please use `adjust_histogram(img, Equalization())` instead.", :histeq) - bins, histogram = imhist(img, nbins, minval, maxval) - cdf = cumsum(histogram[2:end-1]) - img_shape = size(img) - minval == maxval && return map(identity, img) - # Would like to use `map` here, but see https://github.com/timholy/Images.jl/pull/523#issuecomment-235236460 - hist_equalised_img = similar(img) - for I in eachindex(img) - hist_equalised_img[I] = _histeq_pixel_rescale(img[I], cdf, minval, maxval) - end - hist_equalised_img -end - -function histeq(img::AbstractArray, nbins::Integer) - Base.depwarn("`histeq` will be removed in a future release, please use `adjust_histogram(img, Equalization())` instead.", :histeq) - T = _graytype(eltype(img)) - histeq(img, nbins, zero(T), oneunit(T)) -end - -function histeq(img::ImageMeta, nbins::Integer, minval::NumberLike, maxval::NumberLike) - Base.depwarn("`histeq` will be removed in a future release, please use `adjust_histogram(img, Equalization())` instead.", :histeq) - newimg = histeq(arraydata(img), nbins, minval, maxval) - shareproperties(img, newimg) -end - -function histeq(img::ImageMeta, nbins::Integer) - Base.depwarn("`histeq` will be removed in a future release, please use `adjust_histogram(img, Equalization())` instead.", :histeq) - return shareproperties(img, histeq(arraydata(img), nbins)) -end - -function adjust_gamma(img::ImageMeta, gamma::Number) - Base.depwarn("`adjust_gamma` will be removed in a future release, please use `adjust_histogram(img, GammaCorrection())` instead.", :adjust_gamma) - return shareproperties(img, adjust_gamma(arraydata(img), gamma)) -end - -_gamma_pixel_rescale(pixel::T, gamma::Number) where {T<:NumberLike} = pixel ^ gamma - -function _gamma_pixel_rescale(pixel::C, gamma::Number) where C<:Color - yiq = convert(YIQ, pixel) - y = _gamma_pixel_rescale(yiq.y, gamma) - convert(C, YIQ(y, yiq.i, yiq.q)) -end - -function _gamma_pixel_rescale(pixel::C, gamma::Number) where C<:TransparentColor - base_colorant_type(C)(_gamma_pixel_rescale(color(pixel), gamma), alpha(pixel)) -end - -function _gamma_pixel_rescale(original_val::Number, gamma::Number, minval::Number, maxval::Number) - Float32(minval + (maxval - minval) * ((original_val - minval) / (maxval - minval)) ^ gamma) -end - -""" -``` -gamma_corrected_img = adjust_gamma(img, gamma) -``` - -Returns a gamma corrected image. - -# Details - - -Gamma correction is a non-linear transformation given by the relation -```math -f(x) = x^\\gamma \\quad \\text{for} \\; x \\in \\mathbb{R}, \\gamma > 0. -``` -It is called a *power law* transformation because one quantity varies as a power -of another quantity. - -Gamma correction has historically been used to preprocess -an image to compensate for the fact that the intensity of light generated by a -physical device is not usually a linear function of the applied signal but -instead follows a power law [1]. For example, for many Cathode Ray Tubes (CRTs) the -emitted light intensity on the display is approximately equal to the voltage -raised to the power of γ, where γ ∈ [1.8, 2.8]. Hence preprocessing a raw image with -an exponent of 1/γ would have ensured a linear response to brightness. - -Research in psychophysics has also established an [empirical power law -](https://en.wikipedia.org/wiki/Stevens%27s_power_law) between light intensity and perceptual -brightness. Hence, gamma correction often serves as a useful image enhancement -tool. - - -# Options - -Various options for the parameters of this function are described in more detail -below. - -## Choices for `img` - -The `adjust_gamma` function can handle a variety of input types. The returned -image depends on the input type. If the input is an `Image` then the resulting -image is of the same type and has the same properties. - -For coloured images, the input is converted to YIQ type and the Y channel is -gamma corrected. This is the combined with the I and Q channels and the -resulting image converted to the same type as the input. - -## Choice for `gamma` - -The `gamma` value must be a non-zero positive number. - -# Example - -```julia -using Images, ImageView - -# Create an example image consisting of a linear ramp of intensities. -n = 32 -intensities = 0.0:(1.0/n):1.0 -img = repeat(intensities, inner=(20,20))' - -# Brighten the dark tones. -imgadj = adjust_gamma(img,1/2) - -# Display the original and adjusted image. -imshow(img) -imshow(imgadj) -``` - -# 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) - - -See also: [histmatch](@ref),[clahe](@ref), and [imhist](@ref). - - -""" -function adjust_gamma(img::AbstractArray{Gray{T}}, gamma::Number) where T<:Normed - Base.depwarn("`adjust_gamma` will be removed in a future release, please use `adjust_histogram(img, GammaCorrection())` instead.", :adjust_gamma) - raw_type = FixedPointNumbers.rawtype(T) - gamma_inv = 1.0 / gamma - table = zeros(T, typemax(raw_type) + 1) - for i in zero(raw_type):typemax(raw_type) - table[i + 1] = T((i / typemax(raw_type)) ^ gamma_inv) - end - gamma_corrected_img = similar(img) - for I in eachindex(img) - gamma_corrected_img[I] = Gray(table[convert(base_colorant_type(typeof(img[I])){T}, img[I]).val.i + 1]) - end - gamma_corrected_img -end - -function adjust_gamma(img::AbstractArray{T}, gamma::Number) where {T<:Number} - Base.depwarn("`adjust_gamma` will be removed in a future release, please use `adjust_histogram(img, GammaCorrection())` instead.", :adjust_gamma) - return _adjust_gamma(img, gamma, Float64) -end -function adjust_gamma(img::AbstractArray{T}, gamma::Number) where {T<:Colorant} - Base.depwarn("`adjust_gamma` will be removed in a future release, please use `adjust_histogram(img, GammaCorrection())` instead.", :adjust_gamma) - return _adjust_gamma(img, gamma, T) -end - -function _adjust_gamma(img::AbstractArray, gamma::Number, C::Type) - gamma_corrected_img = _fill(oneunit(C), axes(img)) - for I in eachindex(img) - gamma_corrected_img[I] = _gamma_pixel_rescale(img[I], gamma) - end - gamma_corrected_img -end - -function adjust_gamma(img::AbstractArray{T}, gamma::Number, minval::Number, maxval::Number) where T<:Number - Base.depwarn("`adjust_gamma` will be removed in a future release, please use `adjust_histogram(img, GammaCorrection())` instead.", :adjust_gamma) - gamma_corrected_img = _fill(oneunit(T), axes(img)) - for I in eachindex(img) - gamma_corrected_img[I] = _gamma_pixel_rescale(img[I], gamma, minval, maxval) - end - gamma_corrected_img -end - -_fill(val, dim) = fill(val, dim) # fallback -_fill(val, dim::NTuple{N,Base.OneTo}) where {N} = fill(val, map(length, dim)) - -""" -``` -hist_matched_img = histmatch(img, oimg, nbins) -``` - -Returns a histogram matched image with a granularity of `nbins` number -of bins. The first argument `img` is the image to be matched, and the second -argument `oimg` is the image having the desired histogram to be matched to. - -# Details -The purpose of histogram matching is to transform the intensities in a source -image so that the intensities distribute according to the histogram of a -specified target image. If one interprets histograms as piecewise-constant -models of probability density functions (see [imhist](@ref)), then the histogram -matching task can be modelled as the problem of transforming one probability -distribution into another [1]. It turns out that the solution to this -transformation problem involves the cumulative and inverse cumulative -distribution functions of the source and target probability density functions. - -In particular, let the random variables ``x \\thicksim p_{x} `` and ``z -\\thicksim p_{z}`` represent an intensity in the source and target image -respectively, and let - -```math - S(x) = \\int_0^{x}p_{x}(w)\\mathrm{d} w \\quad \\text{and} \\quad - T(z) = \\int_0^{z}p_{z}(w)\\mathrm{d} w -``` -represent their concomitant cumulative disitribution functions. Then the -sought-after mapping ``Q(\\cdot)`` such that ``Q(x) \\thicksim p_{z} `` is given -by - -```math -Q(x) = T^{-1}\\left( S(x) \\right), -``` - -where ``T^{-1}(y) = \\operatorname{min} \\{ x \\in \\mathbb{R} : y \\leq T(x) -\\}`` is the inverse cumulative distribution function of ``T(x)``. - -The mapping suggests that one can conceptualise histogram matching as performing -histogram equalisation on the source and target image and relating the two -equalised histograms. Refer to [histeq](@ref) for more details on histogram -equalisation. - -# Options - -Various options for the parameters of this function are described in more detail -below. - -## Choices for `img` and `oimg` - -The `histmatch` function can handle a variety of input types. The returned -image depends on the input type. If the input is an `Image` then the resulting -image is of the same type and has the same properties. - -For coloured images, the input is converted to -[YIQ](https://en.wikipedia.org/wiki/YIQ) type and the Y channel is gamma -corrected. This is then combined with the I and Q channels and the resulting -image converted to the same type as the input. - -## Choices for `nbins` - -You can specify the total number of bins in the histogram. - -# Example - -```julia -using Images, TestImages, ImageView - -img_source = testimage("mandril_gray") -img_target = adjust_gamma(img_source,1/2) -img_transformed = histmatch(img_source, img_target) -#= - A visual inspection confirms that img_transformed resembles img_target - much more closely than img_source. -=# -imshow(img_source) -imshow(img_target) -imshow(img_transformed) -``` - -# 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) - - -See also: [histeq](@ref),[clahe](@ref), and [imhist](@ref). - -""" -function histmatch(img::ImageMeta, oimg::AbstractArray, nbins::Integer = 400) - Base.depwarn("`histmatch` will be removed in a future release, please use `adjust_histogram(img, Matching())` instead.", :histmatch) - return shareproperties(img, histmatch(arraydata(img), oimg, nbins)) -end - -_hist_match_pixel(pixel::T, bins, lookup_table) where {T<:NumberLike} = T(bins[lookup_table[searchsortedlast(bins, pixel)]]) - -function _hist_match_pixel(pixel::T, bins, lookup_table) where T<:Color - yiq = convert(YIQ, pixel) - y = _hist_match_pixel(yiq.y, bins, lookup_table) - convert(T, YIQ(y, yiq.i, yiq.q)) -end - -_hist_match_pixel(pixel::T, bins, lookup_table) where {T<:TransparentColor} = base_colorant_type(T)(_hist_match_pixel(color(pixel), bins, lookup_table), alpha(pixel)) - -function histmatch(img::AbstractArray{T}, oimg::AbstractArray, nbins::Integer = 400) where T<:Colorant - Base.depwarn("`histmatch` will be removed in a future release, please use `adjust_histogram(img, Matching())` instead.", :histmatch) - el_gray = _graytype(eltype(img)) - oedges, ohist = imhist(oimg, nbins, zero(el_gray), oneunit(el_gray)) - _histmatch(img, oedges, ohist) -end - -function _histmatch(img::AbstractArray, oedges::AbstractRange, ohist::AbstractArray{Int}) - bins, histogram = imhist(img, oedges) - ohist[1] = zero(eltype(ohist)) - ohist[end] = zero(eltype(ohist)) - histogram[1] = zero(eltype(histogram)) - histogram[end] = zero(eltype(histogram)) - cdf = cumsum(histogram) - norm_cdf = cdf / cdf[end] - ocdf = cumsum(ohist) - norm_ocdf = ocdf / ocdf[end] - lookup_table = zeros(Int, length(norm_cdf)) - for I in eachindex(cdf) - lookup_table[I] = argmin(abs.(norm_ocdf .- norm_cdf[I])) - end - hist_matched_img = similar(img) - for I in eachindex(img) - hist_matched_img[I] = _hist_match_pixel(img[I], bins, lookup_table) - end - hist_matched_img -end - -_graytype(::Type{T}) where {T<:Number} = T -_graytype(::Type{C}) where {C<:AbstractGray} = C -_graytype(::Type{C}) where {C<:Colorant} = Gray{eltype(C)} - -""" -``` -hist_equalised_img = clahe(img, nbins, xblocks = 8, yblocks = 8, clip = 3) - -``` - -Performs Contrast Limited Adaptive Histogram Equalisation (CLAHE) on the input -image. It differs from ordinary histogram equalization in the respect that the -adaptive method computes several histograms, each corresponding to a distinct -section of the image, and uses them to redistribute the lightness values of the -image. It is therefore suitable for improving the local contrast and enhancing -the definitions of edges in each region of an image. - -# Details - -Histogram equalisation was initially conceived to improve the contrast in a -single-channel grayscale image. The method transforms the -distribution of the intensities in an image so that they are as uniform as -possible [1]. The natural justification for uniformity -is that the image has better contrast if the intensity levels of an image span -a wide range on the intensity scale. As it turns out, the necessary -transformation is a mapping based on the cumulative histogram---see [histeq](@ref) -for more details. - -A natural extension of histogram equalisation is to apply the contrast -enhancement locally rather than globally [2]. Conceptually, one can imagine that -the process involves partitioning the image into a grid of rectangular regions -and applying histogram equalisation based on the local CDF of each contextual -region. However, to smooth the transition of the pixels from one contextual -region to another, the mapping of a pixel is not done soley based on the local -CDF of its contextual region. Rather, the mapping of a pixel is a bilinear blend -based on the CDF of its contextual region, and the CDFs of the immediate -neighbouring regions. - -In adaptive histogram equalisation the image ``\\mathbf{G}`` is partitioned into -``P \\times Q`` equisized submatrices, -```math -\\mathbf{G} = \\begin{bmatrix} -\\mathbf{G}_{11} & \\mathbf{G}_{12} & \\ldots & \\mathbf{G}_{1C} \\\\ -\\mathbf{G}_{21} & \\mathbf{G}_{22} & \\ldots & \\mathbf{G}_{2C} \\\\ -\\vdots & \\vdots & \\ldots & \\vdots \\\\ -\\mathbf{G}_{R1} & \\mathbf{G}_{R2} & \\ldots & \\mathbf{G}_{RC} \\\\ -\\end{bmatrix}. -``` - -For each submatrix ``\\mathbf{G}_{rc}``, one computes a concomitant CDF, which we -shall denote by ``T_{rc}(G_{i,j})``. In order to determine which CDFs will be -used in the bilinear interpolation step, it is useful to introduce the function - -```math -\\Phi(\\mathbf{G}_{rc}) = \\left( \\phi_{rc}, \\phi'_{rc}\\right) \\triangleq \\left(\\frac{rP}{2}, \\frac{cQ}{2} \\right) -``` - -and to form the sequences ``\\left(\\phi_{11}, \\phi_{21}, \\ldots, \\phi_{R1} \\right)`` -and ``\\left(\\phi'_{11}, \\phi'_{12}, \\ldots, \\phi'_{1C} \\right)``. -For a given pixel ``G_{i,j}(\\omega)``, values of ``r`` and ``c`` are implicitly -defined by the solution to the inequalities -```math -\\phi_{r1} \\le i \\le \\phi_{(r+1)1} \\quad \\text{and} \\quad \\phi'_{1c} \\le j \\le \\phi'_{1(c+1)}. -``` -With ``r`` and ``c`` appropriately defined, the requisite CDFs are given by - -```math -\\begin{aligned} -T_1(v) & \\triangleq T_{rc}(G_{i,j}) \\\\ -T_2(v) & \\triangleq T_{(r+1)c}(G_{i,j}) \\\\ -T_3(v) & \\triangleq T_{(r+1)(c+1)}(G_{i,j}) \\\\ -T_4(v) & \\triangleq T_{r(c+1)}(G_{i,j}). -\\end{aligned} -``` - -Finally, with - -```math -\\begin{aligned} -t & \\triangleq \\frac{i - \\phi_{r1}}{\\phi_{(r+1)1} - \\phi_{r1} } \\\\ -u & \\triangleq \\frac{j - \\phi'_{1c}}{\\phi'_{1(c+1)} - \\phi'_{1c} }, -\\end{aligned} -``` - -the bilinear interpolated transformation that maps an intensity ``v`` at location ``(i,j)`` in the image -to an intensity ``v'`` is given by [3] - -```math -v' \\triangleq \\bar{T}(v) = (1-t) (1-u)T_1(G_{i,j}) + t(1-u)T_2(G_{i,j}) + tuT_3(G_{i,j}) +(1-t)uT_4(G_{i,j}). -``` - -An unfortunate side-effect of contrast enhancement is that it has a tendency to -amplify the level of noise in an image, especially when the magnitude of the -contrast enhancement is very high. The magnitude of contrast enhancement is -associated with the gradient of ``T(\\cdot)``, because the gradient determines the -extent to which consecutive input intensities are stretched across the -grey-level spectrum. One can diminish the level of noise amplification by -limiting the magnitude of the contrast enhancement, that is, by limiting the -magnitude of the gradient. - -Since the derivative of ``T(\\cdot)`` is the empirical density ``\\hat{f}_{G}``, -the slope of the mapping function at any input intensity is proportional to the -height of the histogram ``\\hat{f}_{G}`` at that intensity. Therefore, -limiting the slope of the local mapping function is equivalent to clipping the -height of the histogram. A detailed description of the implementation details -of the clipping process can be found in [2]. - -# Options - -Various options for the parameters of this function are described in more detail -below. - -## Choices for `img` - -The `clahe` function can handle a variety of input types. The returned image -depends on the input type. If the input is an `Image` then the resulting image -is of the same type and has the same properties. - -For coloured images, the input is converted to -[YIQ](https://en.wikipedia.org/wiki/YIQ) type and the Y channel is equalised. -This is the combined with the I and Q channels and the resulting image converted -to the same type as the input. - -## Choices for `nbins` - -You can specify the total number of bins in the histogram of each local region. - -## Choices for `xblocks` and `yblocks` - -The `xblocks` and `yblocks` specify the number of blocks to divide the input -image into in each direction. By default both values are set to eight. - -## Choices for `clip` - -The `clip` parameter specifies the value at which the histogram is clipped. The -default value is three. The excess in the histogram bins with value exceeding -`clip` is redistributed among the other bins. - -# Example - -```julia - -using Images, TestImages, ImageView - -img = testimage("mandril_gray") -imgeq = clahe(img,256, xblocks = 50, yblocks = 50) - -imshow(img) -imshow(imgeq) -``` - -# References -1. R. C. Gonzalez and R. E. Woods. *Digital Image Processing (3rd Edition)*. Upper Saddle River, NJ, USA: Prentice-Hall, 2006. -2. S. M. Pizer, E. P. Amburn, J. D. Austin, R. Cromartie, A. Geselowitz, T. Greer, B. ter Haar Romeny, J. B. Zimmerman and K. Zuiderveld “Adaptive histogram equalization and its variations,” *Computer Vision, Graphics, and Image Processing*, vol. 38, no. 1, p. 99, Apr. 1987. [10.1016/S0734-189X(87)80186-X](https://doi.org/10.1016/s0734-189x(87)80156-1) -3. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. *Numerical Recipes: The Art of Scientific Computing (3rd Edition)*. New York, NY, USA: Cambridge University Press, 2007. - -See also: [histmatch](@ref),[histeq](@ref), [imhist](@ref) and [adjust_gamma](@ref). - - -""" -function clahe(img::AbstractArray{C, 2}, nbins::Integer = 100; xblocks::Integer = 8, yblocks::Integer = 8, clip::Number = 3) where C - Base.depwarn("`clahe` will be removed in a future release, please use `adjust_histogram(img, AdaptiveEqualization())` instead.", :clahe) - h, w = size(img) - y_padded = ceil(Int, h / (2 * yblocks)) * 2 * yblocks - x_padded = ceil(Int, w / (2 * xblocks)) * 2 * xblocks - - img_padded = imresize(img, (y_padded, x_padded)) - - hist_equalised_img = _clahe(img_padded, nbins, xblocks, yblocks, clip) - out = similar(img) - ImageTransformations.imresize!(out, hist_equalised_img) -end - -function clahe(img::ImageMeta, nbins::Integer = 100; xblocks::Integer = 8, yblocks::Integer = 8, clip::Number = 3) - Base.depwarn("`clahe` will be removed in a future release, please use `adjust_histogram(img, AdaptiveEqualization())` instead.", :clahe) - shareproperties(clahe(arraydata(img), nbins, xblocks = xblocks, yblocks = yblocks, clip = clip), img) -end - -function _clahe(img::AbstractArray{C, 2}, nbins::Integer = 100, xblocks::Integer = 8, yblocks::Integer = 8, clip::Number = 3) where C - h, w = size(img) - xb = 0:1:xblocks - 1 - yb = 0:1:yblocks - 1 - blockw = Int(w / xblocks) - blockh = Int(h / yblocks) - temp_cdf = Array{Float64, 1}[] - T = _graytype(eltype(img)) - edges = StatsBase.histrange([Float64(zero(T)), Float64(oneunit(T))], nbins, :left) - - for i in xb - for j in yb - temp_block = img[j * blockh + 1 : (j + 1) * blockh, i * blockw + 1 : (i + 1) * blockw] - _, histogram = imhist(temp_block, edges) - clipped_hist = cliphist(histogram[2:end - 1], clip) - cdf = cumsum(clipped_hist) - n_cdf = cdf / max(convert(eltype(cdf), 1), cdf[end]) - push!(temp_cdf, n_cdf) - end - end - - norm_cdf = reshape(temp_cdf, yblocks, xblocks) - res_img = zeros(C, size(img)) - - #Interpolations - xb = 1:2:xblocks * 2 - 2 - yb = 1:2:yblocks * 2 - 2 - - for j in yb - for i in xb - p_block = img[j * Int(blockh / 2) + 1 : (j + 2) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 2) * Int(blockw / 2)] - bnum_y = floor(Int, (j + 1) / 2) - bnum_x = floor(Int, (i + 1) / 2) - top_left = norm_cdf[bnum_y, bnum_x] - top_right = norm_cdf[bnum_y, bnum_x + 1] - bot_left = norm_cdf[bnum_y + 1, bnum_x] - bot_right = norm_cdf[bnum_y + 1, bnum_x + 1] - temp_block = zeros(C, size(p_block)) - for l in 1:blockw - for m in 1:blockh - temp_block[m, l] = _clahe_pixel_rescale(p_block[m, l], top_left, top_right, bot_left, bot_right, edges, l, m, blockw, blockh) - end - end - res_img[j * Int(blockh / 2) + 1 : (j + 2) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 2) * Int(blockw / 2)] = temp_block - end - end - - #Corners - - block_tl = img[1 : Int(blockh / 2), 1 : Int(blockw / 2)] - corner_tl = map(i -> _clahe_pixel_rescale(i, norm_cdf[1, 1], edges), block_tl) - res_img[1 : Int(blockh / 2), 1 : Int(blockw / 2)] = corner_tl - - block_tr = img[1 : Int(blockh / 2), (xblocks * 2 - 1) * Int(blockw / 2) + 1 : (xblocks * 2) * Int(blockw / 2)] - corner_tr = map(i -> _clahe_pixel_rescale(i, norm_cdf[1, xblocks], edges), block_tr) - res_img[1 : Int(blockh / 2), (xblocks * 2 - 1) * Int(blockw / 2) + 1 : (xblocks * 2) * Int(blockw / 2)] = corner_tr - - block_bl = img[(yblocks * 2 - 1) * Int(blockh / 2) + 1 : (yblocks * 2) * Int(blockh / 2), 1 : Int(blockw / 2)] - corner_bl = map(i -> _clahe_pixel_rescale(i, norm_cdf[yblocks, 1], edges), block_bl) - res_img[(yblocks * 2 - 1) * Int(blockh / 2) + 1 : (yblocks * 2) * Int(blockh / 2), 1 : Int(blockw / 2)] = corner_bl - - block_br = img[(yblocks * 2 - 1) * Int(blockh / 2) + 1 : (yblocks * 2) * Int(blockh / 2), (xblocks * 2 - 1) * Int(blockw / 2) + 1 : (xblocks * 2) * Int(blockw / 2)] - corner_br = map(i -> _clahe_pixel_rescale(i, norm_cdf[yblocks, xblocks], edges), block_br) - res_img[(yblocks * 2 - 1) * Int(blockh / 2) + 1 : (yblocks * 2) * Int(blockh / 2), (xblocks * 2 - 1) * Int(blockw / 2) + 1 : (xblocks * 2) * Int(blockw / 2)] = corner_br - - #Horizontal Borders - - for j in [0, yblocks * 2 - 1] - for i in xb - p_block = img[j * Int(blockh / 2) + 1 : (j + 1) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 2) * Int(blockw / 2)] - bnum_x = floor(Int, (i + 1) / 2) - left = norm_cdf[ceil(Int, (j + 1) / 2), bnum_x] - right = norm_cdf[ceil(Int, (j + 1) / 2), bnum_x + 1] - temp_block = zeros(C, size(p_block)) - for l in 1:blockw - for m in 1:Int(blockh / 2) - temp_block[m, l] = _clahe_pixel_rescale(p_block[m, l], left, right, edges, l, blockw) - end - end - res_img[j * Int(blockh / 2) + 1 : (j + 1) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 2) * Int(blockw / 2)] = temp_block - end - end - - #Vertical Borders - - for i in [0, xblocks * 2 - 1] - for j in yb - p_block = img[j * Int(blockh / 2) + 1 : (j + 2) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 1) * Int(blockw / 2)] - bnum_y = floor(Int, (j + 1) / 2) - top = norm_cdf[bnum_y, ceil(Int, (i + 1) / 2)] - bot = norm_cdf[bnum_y + 1, ceil(Int, (i + 1) / 2)] - temp_block = zeros(C, size(p_block)) - for l in 1:Int(blockw / 2) - for m in 1:blockh - temp_block[m, l] = _clahe_pixel_rescale(p_block[m, l], top, bot, edges, m, blockh) - end - end - res_img[j * Int(blockh / 2) + 1 : (j + 2) * Int(blockh / 2), i * Int(blockw / 2) + 1 : (i + 1) * Int(blockw / 2)] = temp_block - end - end - res_img -end - -_clahe_pixel_rescale(pixel::T, cdf, edges) where {T<:NumberLike} = cdf[searchsortedlast(edges, pixel, Base.Order.Forward)] - -function _clahe_pixel_rescale(pixel::T, first, second, edges, pos, length) where T<:NumberLike - id = searchsortedlast(edges, pixel, Base.Order.Forward) - f = first[id] - s = second[id] - T(((length - pos) * f + (pos - 1) * s) / (length - 1)) -end - -function _clahe_pixel_rescale(pixel::T, top_left, top_right, bot_left, bot_right, edges, i, j, w, h) where T<:NumberLike - id = searchsortedlast(edges, pixel, Base.Order.Forward) - tl = top_left[id] - tr = top_right[id] - bl = bot_left[id] - br = bot_right[id] - r1 = ((w - i) * tl + (i - 1) * tr) / (w - 1) - r2 = ((w - i) * bl + (i - 1) * br) / (w - 1) - T(((h - j) * r1 + (j - 1) * r2) / (h - 1)) -end - -function _clahe_pixel_rescale(pixel::C, args...) where C<:Color - yiq = convert(YIQ, pixel) - y = _clahe_pixel_rescale(yiq.y, args...) - convert(C, YIQ(y, yiq.i, yiq.q)) -end - -function _clahe_pixel_rescale(pixel::C, args...) where C<:TransparentColor - base_colorant_type(C)(_clahe_pixel_rescale(color(pixel), args...), alpha(pixel)) -end - -""" -``` -clipped_hist = cliphist(hist, clip) -``` - -Clips the histogram above a certain value `clip`. The excess left in the bins -exceeding `clip` is redistributed among the remaining bins. -""" -function cliphist(hist::AbstractArray{T, 1}, clip::Number) where T - Base.depwarn("`cliphist` will be removed in a future release.", :cliphist) - hist_length = length(hist) - excess = sum(map(i -> i > clip ? i - clip : zero(i - clip), hist)) - increase = excess / hist_length - clipped_hist = zeros(Float64, hist_length) - removed_sum = zero(Float64) - - for i in 1:hist_length - if hist[i] > clip - increase - clipped_hist[i] = Float64(clip) - if hist[i] < clip - removed_sum += hist[i] - (clip - increase) - end - else - clipped_hist[i] = hist[i] + Float64(increase) - removed_sum += increase - end - end - - leftover = excess - removed_sum - - while true - oleftover = leftover - for i in 1:hist_length - leftover <= 0 && break - step = ceil(Int, max(1, hist_length / leftover)) - for h in i:step:hist_length - leftover <= 0 && break - diff = clip - clipped_hist[h] - if diff > 1 - clipped_hist[h] += 1 - leftover -= 1 - elseif diff > 0 - clipped_hist[h] = clip - leftover -= diff - end - end - end - (leftover <= 0 || leftover >= oleftover) && break - end - clipped_hist -end - -# phantom images - -""" -``` -phantom = shepp_logan(N,[M]; highContrast=true) -``` - -output the NxM Shepp-Logan phantom, which is a standard test image usually used -for comparing image reconstruction algorithms in the field of computed -tomography (CT) and magnetic resonance imaging (MRI). If the argument M is -omitted, the phantom is of size NxN. When setting the keyword argument -`highConstrast` to false, the CT version of the phantom is created. Otherwise, -the high contrast MRI version is calculated. -""" -function shepp_logan(M,N; highContrast=true) - # Initially proposed in Shepp, Larry; B. F. Logan (1974). - # "The Fourier Reconstruction of a Head Section". IEEE Transactions on Nuclear Science. NS-21. - - Base.depwarn("`shepp_logan` will be removed in a future release, please use `TestImages.shepp_logan` instead.", :shepp_logan) - - P = zeros(M,N) - - x = range(-1, stop=1, length=M)' - y = range(1, stop=-1, length=N) - - centerX = [0, 0, 0.22, -0.22, 0, 0, 0, -0.08, 0, 0.06] - centerY = [0, -0.0184, 0, 0, 0.35, 0.1, -0.1, -0.605, -0.605, -0.605] - majorAxis = [0.69, 0.6624, 0.11, 0.16, 0.21, 0.046, 0.046, 0.046, 0.023, 0.023] - minorAxis = [0.92, 0.874, 0.31, 0.41, 0.25, 0.046, 0.046, 0.023, 0.023, 0.046] - theta = [0, 0, -18.0, 18.0, 0, 0, 0, 0, 0, 0] - - # original (CT) version of the phantom - grayLevel = [2, -0.98, -0.02, -0.02, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01] - - if(highContrast) - # high contrast (MRI) version of the phantom - grayLevel = [1, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] - end - - for l=1:length(theta) - P += grayLevel[l] * ( - ( (cos(theta[l] / 360*2*pi) * (x .- centerX[l]) .+ - sin(theta[l] / 360*2*pi) * (y .- centerY[l])) / majorAxis[l] ).^2 .+ - ( (sin(theta[l] / 360*2*pi) * (x .- centerX[l]) .- - cos(theta[l] / 360*2*pi) * (y .- centerY[l])) / minorAxis[l] ).^2 .< 1 ) - end - - return P -end - -shepp_logan(N;highContrast=true) = shepp_logan(N,N;highContrast=highContrast) - -# `complement` is now moved to ColorVectorSpace but we still want to keep backward compatibility -# until Images 1.0 -function ColorVectorSpace.complement(x::AbstractArray) - Base.depwarn("`complement(img)` is deprecated, please use the broadcasting version `complement.(img)`", :complement) - complement.(x) -end - -@deprecate forwarddiffx(X) ImageBase.FiniteDiff.fdiff(X, dims=2, boundary=:zero) -@deprecate forwarddiffy(X) ImageBase.FiniteDiff.fdiff(X, dims=1, boundary=:zero) -@deprecate backdiffx(X) ImageBase.FiniteDiff.fdiff(X, dims=2, rev=true, boundary=:zero) -@deprecate backdiffy(X) ImageBase.FiniteDiff.fdiff(X, dims=1, rev=true, boundary=:zero) -@deprecate div(A::AbstractArray{<:Any,3}) ImageBase.FiniteDiff.fdiv(view(A, :, :, 1), view(A, :, :, 2)) false -@deprecate imROF(img::AbstractMatrix, λ::Number, iterations::Integer) ImageFiltering.Models.solve_ROF_PD(img, λ, iterations) - - -# This is now replaced by ImageTransformations and Interpolations -using ImageTransformations.Interpolations: BSpline, Linear, interpolate -function bilinear_interpolation(img::AbstractArray{T,N}, xs::Vararg{Number, N}) where {T,N} - Base.depwarn("`bilinear_interpolation` is deprecated, please use `warp`, `imresize` from ImageTransformations or `extrapolate`, `interpolate` from Interpolations", :bilinear_interpolation) - pad_ax = map(axes(img), xs) do ax, x - min(floor(Int, x), first(ax)):max(ceil(Int, x), last(ax)) - end - padded_img = PaddedView(zero(T), img, pad_ax) - itp = interpolate(padded_img, BSpline(Linear())) - return itp(xs...) -end -export bilinear_interpolation - -import ImageFiltering: findlocalmaxima, findlocalminima, blob_LoG -dims2window(img, dims) = ntuple(d -> d ∈ dims ? 3 : 1, ndims(img)) -function find_depwarn(sym, dims, window, edges) - Base.depwarn("`$sym(img, $dims, $edges)` is deprecated, please use `$sym(img; window=$window, edges=$edges)`.\nSee the documentation for details about `window`.", sym) -end -function findlocalmaxima(img::AbstractArray, region::Union{Tuple{Int,Vararg{Int}},Vector{Int},UnitRange{Int},Int}, edges=true) - window = dims2window(img, region) - find_depwarn(:findlocalmaxima, region, window, edges) - findlocalmaxima(img; window=window, edges=edges) -end -function findlocalminima(img::AbstractArray, region::Union{Tuple{Int,Vararg{Int}},Vector{Int},UnitRange{Int},Int}, edges=true) - window = dims2window(img, region) - find_depwarn(:findlocalminima, region, window, edges) - findlocalminima(img; window=window, edges=edges) -end -@deprecate blob_LoG(img::AbstractArray{T,N}, σscales::Union{AbstractVector,Tuple}, - edges::Union{Bool,Tuple{Vararg{Bool}}}, σshape=ntuple(d->1, Val(N))) where {T,N} blob_LoG(img, σscales; edges=edges, σshape=(σshape...,), rthresh=0) -@deprecate blob_LoG(img::AbstractArray{T,N}, σscales::Union{AbstractVector,Tuple}, σshape::Union{AbstractVector,NTuple{N,Real}}) where {T,N} blob_LoG(img, σscales; edges=(true, ntuple(d->false,Val(N))...), σshape=(σshape...,), rthresh=0) - - -import Statistics: std, var -@deprecate var(A::AbstractArray{C}; kwargs...) where C<:ColorVectorSpace.MathTypes varmult(⊙, A; kwargs...) -@deprecate std(A::AbstractArray{C}; kwargs...) where C<:ColorVectorSpace.MathTypes stdmult(⊙, A; kwargs...) -# Doing stats on non-MathTypes is a bad idea -function var(A::AbstractArray{C}; dims=nothing, kwargs...) where C<:Colorant - dims === nothing || error("dims was never supported, but is now for MathTypes") - newdims = 2:ndims(A)+1 - Cbase = base_colorant_type(C) - forced_depwarn("`var(A::AbstractArray{<:Colorant})` is deprecated (and not recommended, use a MathTypes colorspace instead), please use `colorview($Cbase, var(channelview(A); dims=$newdims))[]` instead.", :var) - return colorview(Cbase, var(channelview(A); dims=newdims, kwargs...))[] -end -function std(A::AbstractArray{C}; dims=nothing, kwargs...) where C<:Colorant - dims === nothing || error("dims was never supported, but is now for MathTypes") - newdims = 2:ndims(A)+1 - Cbase = base_colorant_type(C) - forced_depwarn("`std(A::AbstractArray{<:Colorant})` is deprecated (and not recommended, use a MathTypes colorspace instead), please use `colorview($Cbase, std(channelview(A); dims=$newdims))[]` instead.", :var) - return colorview(Cbase, std(channelview(A); dims=newdims, kwargs...))[] -end - -# Integral arrays -@deprecate integral_image(A) IntegralArray(A) -@deprecate boxdiff(Ai::IntegralArray{T,2}, y::UnitRange, x::UnitRange) where T Ai[first(y)..last(y), first(x)..last(x)] -@deprecate boxdiff(Ai::IntegralArray{T,2}, tl::CartesianIndex, br::CartesianIndex) where T Ai[tl[1]..br[1], tl[2]..br[2]] -@deprecate boxdiff(Ai::IntegralArray{T,2}, tl_y::Integer, tl_x::Integer, br_y::Integer, br_x::Integer) where T Ai[tl_y..br_y, tl_x..br_x] - - -function imaverage(filter_size=(3,3)) - Base.depwarn("imaverage(m, n) is deprecated, use `Kernel.box` or an `IntegralImage`.", :imaverage) - Kernel.box(filter_size) -end - -@deprecate imlaplacian(alpha::Number) Kernel.laplacian2d(alpha) false +@deprecate otsu_threshold(img::AbstractArray{T}, nbins::Int = 256) where {T<:Union{Gray,Real}} T(find_threshold(img, Otsu(); nbins)) +@deprecate yen_threshold(img::AbstractArray{T}, nbins::Int = 256) where {T<:Union{Gray,Real}} T(find_threshold(img, Yen(); nbins)) diff --git a/src/labeledarrays.jl b/src/labeledarrays.jl deleted file mode 100644 index 779c1a58..00000000 --- a/src/labeledarrays.jl +++ /dev/null @@ -1,48 +0,0 @@ -struct ColorizedArray{C<:Colorant,N,A<:AbstractArray,L<:AbstractArray} <: AbstractArray{C,N} - intensity::A - label::L -end - -""" - ColorizedArray(intensity, label::IndirectArray) -> A - -Create an array, combining a `label` array (where each pixel is -assigned one of a list of discrete colors) and an `intensity` array -(where each pixel has a scalar value). `A` satisfies - - A[i,j,...] = intensity[i,j,...] * label[i,j,...] - -The label array "tinges" the grayscale intensity with the color -associated with that point's label. - -This computation is performed lazily, as to be suitable even for large arrays. -""" -function ColorizedArray(intensity, label::IndirectArray{C,N}) where {C<:Colorant,N} - Base.depwarn("`ColorizedArray(intensity, label)` is deprecated, use `mappedarray(*, intensity, label)` from `MappedArrays` instead", :ColorizedArray) - - axes(intensity) == axes(label) || throw(DimensionMismatch("intensity and label must have the same axes, got $(axes(intensity)) and $(axes(label))")) - CI = typeof(zero(C)*zero(eltype(intensity))) - ColorizedArray{CI,N,typeof(intensity),typeof(label)}(intensity, label) -end - -# TODO: an implementation involving AxisArray that matches the shared -# axes, and therefore allows `label` to be of lower dimensionality -# than `intensity`. - -intensitytype(::Type{ColorizedArray{C,N,A,L}}) where {C<:Colorant,N,A<:AbstractArray,L<:AbstractArray} = A -labeltype(::Type{ColorizedArray{C,N,A,L}}) where {C<:Colorant,N,A<:AbstractArray,L<:AbstractArray} = L - -Base.size(A::ColorizedArray) = size(A.intensity) -Base.axes(A::ColorizedArray) = axes(A.intensity) -Base.IndexStyle(::Type{CA}) where {CA<:ColorizedArray} = IndexStyle(IndexStyle(intensitytype(CA)), IndexStyle(labeltype(CA))) - -@inline function Base.getindex(A::ColorizedArray, i::Integer) - @boundscheck checkbounds(A, i) - @inbounds ret = A.intensity[i]*A.label[i] - ret -end -@inline function Base.getindex(A::ColorizedArray{C,N}, I::Vararg{Int,N}) where {C,N} - @boundscheck checkbounds(A, I...) - @inbounds ret = A.intensity[I...]*A.label[I...] - ret -end diff --git a/test/algorithms.jl b/test/algorithms.jl index 69867650..c5f33012 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -9,7 +9,7 @@ using ImageBase.FiniteDiff: fdiff for T in (N0f8, Float32) A = rand(RGB{T}, 5, 4) Ac = channelview(A) - s = @suppress_err std(A) + s = stdmult(⊙, A) @test red(s) ≈ std(Ac[1,:,:]) @test green(s) ≈ std(Ac[2,:,:]) @test blue(s) ≈ std(Ac[3,:,:]) @@ -99,40 +99,15 @@ using ImageBase.FiniteDiff: fdiff @test img2 ≈ A end - - # deprecated - @suppress_err @testset "Phantoms" begin - P = [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; - 0.0 0.0 1.0 0.2 0.2 1.0 0.0 0.0; - 0.0 0.0 0.2 0.3 0.3 0.2 0.0 0.0; - 0.0 0.0 0.2 0.0 0.2 0.2 0.0 0.0; - 0.0 0.0 0.2 0.0 0.0 0.2 0.0 0.0; - 0.0 0.0 0.2 0.2 0.2 0.2 0.0 0.0; - 0.0 0.0 1.0 0.2 0.2 1.0 0.0 0.0; - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ] - Q = Images.shepp_logan(8) - @test norm((P-Q)[:]) < 1e-10 - P = [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; - 0.0 0.0 2.0 1.02 1.02 2.0 0.0 0.0; - 0.0 0.0 1.02 1.03 1.03 1.02 0.0 0.0; - 0.0 0.0 1.02 1.0 1.02 1.02 0.0 0.0; - 0.0 0.0 1.02 1.0 1.0 1.02 0.0 0.0; - 0.0 0.0 1.02 1.02 1.02 1.02 0.0 0.0; - 0.0 0.0 2.0 1.02 1.02 2.0 0.0 0.0; - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ] - Q = Images.shepp_logan(8, highContrast=false) - @test norm((P-Q)[:]) < 1e-10 - end - # functionality moved to ImageTransformations # tests are here as well to make sure everything # is exported properly. @testset "Image resize" begin img = zeros(10,10) - img2 = Images.imresize(img, (5,5)) + img2 = imresize(img, (5,5)) @test length(img2) == 25 img = rand(RGB{Float32}, 10, 10) - img2 = Images.imresize(img, (6,7)) + img2 = imresize(img, (6,7)) @test size(img2) == (6,7) @test eltype(img2) == RGB{Float32} end @@ -174,69 +149,6 @@ using ImageBase.FiniteDiff: fdiff @test sort(B)==sort(C) end - # deprecated - @suppress_err @testset "Thresholding" begin - - #otsu_threshold - img = testimage("cameraman") - thres = otsu_threshold(img) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), convert(N0f8, 87/256), atol=eps(N0f8)) - thres = otsu_threshold(img, 512) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), convert(N0f8, 87/256), atol=eps(N0f8)) - - img = map(x->convert(Gray{Float64}, x), img) - thres = otsu_threshold(img) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), 87/256, atol=0.01) - thres = otsu_threshold(img, 512) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), 87/256, atol=0.01) - - img = map(x->convert(Float64, x), img) - thres = otsu_threshold(img) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), 87/256, atol=0.01) - thres = otsu_threshold(img, 512) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), 87/256, atol=0.01) - - #test for multidimension arrays - img = rand(Float64, 10, 10, 3) - @test otsu_threshold(img) == otsu_threshold(cat(img[:,:,1], img[:,:,2], img[:,:,3], dims=1)) - - #yen_threshold - img = testimage("cameraman") - thres = yen_threshold(img) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), convert(N0f8, 199/256), atol=eps(N0f8)) - thres = yen_threshold(img, 512) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), convert(N0f8, 199/256), atol=eps(N0f8)) - - img = map(x->convert(Gray{Float64}, x), img) - thres = yen_threshold(img) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), 199/256, atol=0.01) - thres = yen_threshold(img, 512) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), 199/256, atol=0.01) - - img = map(x->convert(Float64, x), img) - thres = yen_threshold(img) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), 199/256, atol=0.01) - thres = yen_threshold(img, 512) - @test typeof(thres) == eltype(img) - @test ≈(gray(thres), 199/256, atol=0.01) - - img = rand(Float64, 10, 10, 3) - @test yen_threshold(img) == yen_threshold(cat(img[:,:,1], img[:,:,2], img[:,:,3], dims=1)) - - img = zeros(Gray{Float64},10,10,3) - @test yen_threshold(img) == 0 - end end nothing diff --git a/test/arrays.jl b/test/arrays.jl index 03873b34..6dbd7565 100644 --- a/test/arrays.jl +++ b/test/arrays.jl @@ -4,29 +4,26 @@ using Images.MappedArrays: mappedarray @testset "ColorizedArray" begin intensity = [0.1 0.3; 0.2 0.4] labels = IndirectArray([1 2; 2 1], [RGB(1,0,0), RGB(0,0,1)]) - colorized_A = @suppress_err ColorizedArray(intensity, labels) - mapped_A = mappedarray(*, intensity, labels) + A = mappedarray(*, intensity, labels) target = intensity .* labels - for A in (colorized_A, mapped_A) - @test eltype(A) == RGB{Float64} - @test size(A) == (2,2) - @test axes(A) == (Base.OneTo(2), Base.OneTo(2)) - for i = 1:4 - @test A[i] === target[i] - end - for j = 1:2, i = 1:2 - @test A[i,j] === target[i,j] - end - for (a,t) in zip(A, target) - @test a === t - end + @test eltype(A) == RGB{Float64} + @test size(A) == (2,2) + @test axes(A) == (Base.OneTo(2), Base.OneTo(2)) + for i = 1:4 + @test A[i] === target[i] + end + for j = 1:2, i = 1:2 + @test A[i,j] === target[i,j] + end + for (a,t) in zip(A, target) + @test a === t end - intensity1 = [0.1 0.3 0.6; 0.2 0.4 0.1] + intensity1 = [0.1 0.3 0.6; 0.2 0.4 0.1] labels1 = IndirectArray([1 2 3; 2 1 3], [RGB(1,0,0), RGB(0,0,1),RGB(0,1,0)]) - A1 = @suppress_err ColorizedArray(intensity1, labels1) + A1 = mappedarray(*, intensity1, labels1) target1 = intensity1 .* labels1 @test eltype(A1) == RGB{Float64} @test size(A1) == (2,3) diff --git a/test/deprecated.jl b/test/deprecated.jl index a5615f11..8ec5582f 100644 --- a/test/deprecated.jl +++ b/test/deprecated.jl @@ -1,115 +1,67 @@ # When we deprecate old APIs, we copy the old tests here to make sure they still work @testset "deprecations" begin + # deprecated + @suppress_err @testset "Thresholding" begin - @testset "Complement" begin - # deprecated (#690) - img = Gray{N0f16}.([0.01164 0.01118; 0.01036 0.01187]) - @test all(complement(img) .== 1 .- img) - end - - @testset "Stats" begin - a = [0.1, 0.2, 0.1] - @test var(Gray.(a)) == Gray(var(a)) - @test std(Gray.(a)) == Gray(std(a)) - an0f8 = N0f8.(a) - @test var(Gray.(an0f8)) == Gray(var(an0f8)) - @test std(Gray.(an0f8)) == Gray(std(an0f8)) - c = [RGB(1, 0, 0), RGB(0, 0, 1)] - @test var(c) == varmult(⊙, c) - @test std(c) ≈ mapc(sqrt, varmult(⊙, c)) - # We don't want to support this next one at all: - chsv = HSV.(c) - @test var(chsv) == HSV{Float32}(28800.0f0, 0.0f0, 0.0f0) - @test std(chsv) == mapc(sqrt, HSV{Float32}(28800.0f0, 0.0f0, 0.0f0)) - end + #otsu_threshold + img = testimage("cameraman") + thres = otsu_threshold(img) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), convert(N0f8, 87/256), atol=eps(N0f8)) + thres = otsu_threshold(img, 512) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), convert(N0f8, 87/256), atol=eps(N0f8)) - @testset "Restrict with vector dimensions" begin - imgcol = colorview(RGB, rand(3,5,6)) - imgmeta = ImageMeta(imgcol, myprop=1) - @test isa(restrict(imgmeta, [1, 2]), ImageMeta) - end - - @testset "Interpolations" begin - img = zeros(Float64, 5, 5) - @test bilinear_interpolation(img, 4.5, 5.5) == 0.0 - @test bilinear_interpolation(img, 4.5, 3.5) == 0.0 + img = map(x->convert(Gray{Float64}, x), img) + thres = otsu_threshold(img) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), 87/256, atol=0.01) + thres = otsu_threshold(img, 512) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), 87/256, atol=0.01) - for i in [1.0, 2.0, 5.0, 7.0, 9.0] - img = ones(Float64, 5, 5) * i - @test (bilinear_interpolation(img, 3.5, 4.5) == i) - @test (bilinear_interpolation(img, 3.2, 4) == i) # X_MAX == X_MIN - @test (bilinear_interpolation(img, 3.2, 4) == i) # Y_MAX == Y_MIN - @test (bilinear_interpolation(img, 3.2, 4) == i) # BOTH EQUAL - @test (bilinear_interpolation(img, 2.8, 1.9) == i) - # One dim out of bounds - @test isapprox(bilinear_interpolation(img, 0.5, 1.5), 0.5 * i) - @test isapprox(bilinear_interpolation(img, 0.5, 1.6), 0.5 * i) - @test isapprox(bilinear_interpolation(img, 0.5, 1.8), 0.5 * i) - # Both out of bounds (corner) - @test isapprox(bilinear_interpolation(img, 0.5, 0.5), 0.25 * i) - end + img = map(x->convert(Float64, x), img) + thres = otsu_threshold(img) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), 87/256, atol=0.01) + thres = otsu_threshold(img, 512) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), 87/256, atol=0.01) - img = reshape(1.0:1.0:25.0, 5, 5) - @test bilinear_interpolation(img, 1.5, 2) == 6.5 - @test bilinear_interpolation(img, 2, 1.5) == 4.5 - @test bilinear_interpolation(img, 2, 1) == 2.0 - @test bilinear_interpolation(img, 1.5, 2.5) == 9.0 - @test bilinear_interpolation(img, 1.5, 3.5) == 14.0 - @test bilinear_interpolation(img, 1.5, 4.5) == 19.0 - @test bilinear_interpolation(img, 1.5, 5.5) == 10.75 - end + #test for multidimension arrays + img = rand(Float64, 10, 10, 3) + @test otsu_threshold(img) == otsu_threshold(cat(img[:,:,1], img[:,:,2], img[:,:,3], dims=1)) - @testset "boxdiff" begin - a = zeros(10, 10) - int_img = integral_image(a) - @test all(int_img == a) + #yen_threshold + img = testimage("cameraman") + thres = yen_threshold(img) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), convert(N0f8, 199/256), atol=eps(N0f8)) + thres = yen_threshold(img, 512) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), convert(N0f8, 199/256), atol=eps(N0f8)) - a = ones(10,10) - int_img = integral_image(a) - chk = Array(1:10) - @test all([vec(int_img[i, :]) == chk * i for i in 1:10]) + img = map(x->convert(Gray{Float64}, x), img) + thres = yen_threshold(img) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), 199/256, atol=0.01) + thres = yen_threshold(img, 512) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), 199/256, atol=0.01) - int_sum = boxdiff(int_img, 1, 1, 5, 2) - @test int_sum == 10.0 - int_sum = boxdiff(int_img, 1:5, 1:2) - @test int_sum == 10.0 - int_sum = boxdiff(int_img, CartesianIndex((1, 1)), CartesianIndex((5, 2))) - @test int_sum == 10.0 - int_sum = boxdiff(int_img, 1, 1, 2, 5) - @test int_sum == 10.0 - int_sum = boxdiff(int_img, 1:2, 1:5) - @test int_sum == 10.0 - int_sum = boxdiff(int_img, CartesianIndex((1, 1)), CartesianIndex((2, 5))) - @test int_sum == 10.0 - int_sum = boxdiff(int_img, 4, 4, 8, 8) - @test int_sum == 25.0 - int_sum = boxdiff(int_img, 4:8, 4:8) - @test int_sum == 25.0 - int_sum = boxdiff(int_img, CartesianIndex((4, 4)), CartesianIndex((8, 8))) - @test int_sum == 25.0 + img = map(x->convert(Float64, x), img) + thres = yen_threshold(img) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), 199/256, atol=0.01) + thres = yen_threshold(img, 512) + @test typeof(thres) == eltype(img) + @test ≈(gray(thres), 199/256, atol=0.01) - a = reshape(1:100, 10, 10) - int_img = integral_image(a) - @test int_img[diagind(int_img)] == Array([1, 26, 108, 280, 575, 1026, 1666, 2528, 3645, 5050]) + img = rand(Float64, 10, 10, 3) + @test yen_threshold(img) == yen_threshold(cat(img[:,:,1], img[:,:,2], img[:,:,3], dims=1)) - int_sum = boxdiff(int_img, 1, 1, 3, 3) - @test int_sum == 108 - int_sum = boxdiff(int_img, 1:3, 1:3) - @test int_sum == 108 - int_sum = boxdiff(int_img, CartesianIndex((1, 1)), CartesianIndex((3, 3))) - @test int_sum == 108 - int_sum = boxdiff(int_img, 1, 1, 5, 2) - @test int_sum == 80 - int_sum = boxdiff(int_img, 1:5, 1:2) - @test int_sum == 80 - int_sum = boxdiff(int_img, CartesianIndex((1, 1)), CartesianIndex((5, 2))) - @test int_sum == 80 - int_sum = boxdiff(int_img, 4, 4, 8, 8) - @test int_sum == 1400 - int_sum = boxdiff(int_img, 4:8, 4:8) - @test int_sum == 1400 - int_sum = boxdiff(int_img, CartesianIndex((4, 4)), CartesianIndex((8, 8))) - @test int_sum == 1400 + img = zeros(Gray{Float64},10,10,3) + @test yen_threshold(img) == 0 end end diff --git a/test/exposure.jl b/test/exposure.jl index f6d7cb7d..6efef3c3 100644 --- a/test/exposure.jl +++ b/test/exposure.jl @@ -12,32 +12,32 @@ eye(m,n) = Matrix{Float64}(I,m,n) # should have pixel values between 0 and 1. Still, it doesn't # hurt to get the integer case right too. img = 1:10 - bins, hist = imhist(img, 10) + bins, hist = build_histogram(img, 10) @test length(hist) == length(bins)+1 - @test bins == 1.0:1.0:11.0 - @test hist == [0,1, 1, 1, 1, 1, 1, 1, 1, 1, 1,0] - bins, hist = imhist(img, 5, 2, 6) + @test_broken bins == 1.0:1.0:11.0 + @test_broken hist == [0,1, 1, 1, 1, 1, 1, 1, 1, 1, 1,0] + bins, hist = build_histogram(img, 5; minval=2, maxval=6) @test length(hist) == length(bins)+1 - @test bins == 2.0:1.0:7.0 - @test hist == [1, 1, 1, 1, 1, 1, 4] + @test_broken bins == 2.0:1.0:7.0 + @test_broken hist == [1, 1, 1, 1, 1, 1, 4] img = reshape(0:99, 10, 10) - bins, hist = imhist(img, 10) + bins, hist = build_histogram(img, 10) @test length(hist) == length(bins)+1 - @test bins == 0.0:10.0:100.0 - @test hist == [0, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 0] - bins, hist = imhist(img, 7, 25, 59) + @test_broken bins == 0.0:10.0:100.0 + @test_broken hist == [0, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 0] + bins, hist = build_histogram(img, 7; minval=25, maxval=59) @test length(hist) == length(bins)+1 - @test bins == 25.0:5.0:60.0 - @test hist == [25, 5, 5, 5, 5, 5, 5, 5, 40] + @test_broken bins == 25.0:5.0:60.0 + @test_broken hist == [25, 5, 5, 5, 5, 5, 5, 5, 40] # Test the more typical case img = reinterpret(Gray{N0f8}, [0x20,0x40,0x80,0xd0]) - @test imhist(img, 5) == (0.0:0.2:1.0,[0,1,1,1,0,1,0]) + @test_broken build_histogram(img, 5) == (0.0:0.2:1.0,[0,1,1,1,0,1,0]) img = reinterpret(Gray{N0f8}, [0x00,0x40,0x80,0xd0]) - @test imhist(img, 5) == (0.0:0.2:1.0,[0,1,1,1,0,1,0]) + @test_broken build_histogram(img, 5) == (0.0:0.2:1.0,[0,1,1,1,0,1,0]) img = reinterpret(Gray{N0f8}, [0x00,0x40,0x80,0xff]) - @test all((≈).(imhist(img, 6), (0.0:0.2:1.2,[0,1,1,1,0,0,1,0]))) + @test_broken all((≈).(build_histogram(img, 6), (0.0:0.2:1.2,[0,1,1,1,0,0,1,0]))) # Consider an image where each intensity occurs only once and vary the number # of bins used in the histogram in powers of two. With the exception of the @@ -46,7 +46,7 @@ eye(m,n) = Matrix{Float64}(I,m,n) bins = [2^i for i = 8:-1:1] for i = 1:length(bins) for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) - edges, counts = build_histogram(T.(collect(0:1/255:1)),bins[i],0,1) + edges, counts = build_histogram(T.(collect(0:1/255:1)), bins[i]; minval=0, maxval=1) @test length(edges) == length(counts) - 1 @test all(counts[1:end] .== expected_counts[i]) && counts[0] == 0 @test axes(counts) == (0:length(edges),) @@ -56,14 +56,14 @@ eye(m,n) = Matrix{Float64}(I,m,n) for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64}) imgg = collect(0:1/255:1) img = colorview(RGB,imgg,imgg,imgg) - edges, counts = build_histogram(T.(img),bins[i],0,1) + edges, counts = build_histogram(T.(img), bins[i]; minval=0, maxval=1) @test length(edges) == length(counts) - 1 @test all(counts[1:end] .== expected_counts[i]) && counts[0] == 0 @test axes(counts) == (0:length(edges),) end # Consider also integer-valued images. - edges, counts = build_histogram(0:1:255,bins[i],0,255) + edges, counts = build_histogram(0:1:255, bins[i]; minval=0, maxval=255) @test length(edges) == length(counts) - 1 @test all(counts[1:end] .== expected_counts[i]) && counts[0] == 0 @test axes(counts) == (0:length(edges),) @@ -79,59 +79,63 @@ eye(m,n) = Matrix{Float64}(I,m,n) #DataTypes img = oneunits(Gray{Float64}, 10, 10) - ret = histeq(img, 100) + ret = adjust_histogram(img, Equalization(;nbins=100)) @test img == ret @test eltype(ret) == eltype(img) img = oneunits(Gray{N0f8}, 10, 10) - ret = histeq(img, 100) + ret = adjust_histogram(img, Equalization(;nbins=100)) @test img == ret @test eltype(ret) == eltype(img) img = oneunits(Gray{N0f16}, 10, 10) - ret = histeq(img, 100) + ret = adjust_histogram(img, Equalization(;nbins=100)) @test eltype(ret) == eltype(img) - img = oneunits(AGray{N0f8}, 10, 10) - ret = histeq(img, 100) - @test img == ret - @test eltype(ret) == eltype(img) + @info "Test skippped due to lack of support for `AGray` in adjust_histogram" + # img = oneunits(AGray{N0f8}, 10, 10) + # ret = adjust_histogram(img, Equalization(;nbins=100)) + # @test img == ret + # @test eltype(ret) == eltype(img) img = oneunits(RGB{N0f8}, 10, 10) - ret = histeq(img, 100) + ret = adjust_histogram(img, Equalization(;nbins=100)) @test img == ret @test eltype(ret) == eltype(img) img = oneunits(RGB{N0f16}, 10, 10) - ret = histeq(img, 100) + ret = adjust_histogram(img, Equalization(;nbins=100)) @test img == ret @test eltype(ret) == eltype(img) img = oneunits(RGB{Float64}, 10, 10) - ret = histeq(img, 100) + ret = adjust_histogram(img, Equalization(;nbins=100)) @test all(map((i, r) -> isapprox(i, r), img, ret)) @test eltype(ret) == eltype(img) - img = oneunits(ARGB{N0f8}, 10, 10) - ret = histeq(img, 100) - @test img == ret - @test eltype(ret) == eltype(img) - - #Working + @info "Test skipped due to lack of support for `ARGB` in adjust_histogram" + # img = oneunits(ARGB{N0f8}, 10, 10) + # ret = adjust_histogram(img, Equalization(;nbins=100)) + # @test img == ret + # @test eltype(ret) == eltype(img) img = zeros(10, 10) for i in 1:10 img[i, :] .= 10 * (i - 1) end - @test img == histeq(img, 10, 0, 90) + @test_broken img == adjust_histogram(img, Equalization(;nbins=10, minval=0, maxval=90)) - ret = histeq(img, 2, 0, 90) + ret = adjust_histogram(img, Equalization(;nbins=2, minval=0, maxval=90)) @test all(ret[1:5, :] .== 0) @test all(ret[6:10, :] .== 90) - ret = histeq(img, 5, 0, 90) + ret = adjust_histogram(img, Equalization(;nbins=5, minval=0, maxval=90)) for i in 1:2:10 - @test all(ret[i:i+1, :] .== 22.5 * floor(i / 2)) + if i ∈ (1, 9) + @test all(ret[i:i+1, :] .== 22.5 * floor(i / 2)) + else + @test_broken all(ret[i:i+1, :] .== 22.5 * floor(i / 2)) + end end img = [0.0, 21.0, 29.0, 38.0, 38.0, 51.0, 66.0, 79.0, 79.0, 91.0, @@ -146,11 +150,11 @@ eye(m,n) = Matrix{Float64}(I,m,n) 15.0, 27.0, 36.0, 36.0, 44.0, 69.0, 74.0, 74.0, 86.0, 92.0] img = reshape(img, 10, 10)' - ret = histeq(img, 10, 0, 99) - cdf = cumsum(imhist(img, 10)[2][2:end-1]) + ret = adjust_histogram(img, Equalization(; nbins=10, minval=0, maxval=99)) + cdf = cumsum(build_histogram(img, 10)[2][2:end-1]) @test all(ret[1:cdf[1]] .== 0.0) for i in 1:(size(cdf)[1]-1) - @test all(ret[cdf[i] + 1 : cdf[i + 1]] .== (cdf[i + 1] - cdf[1]) * 99.0 / (cdf[end] - cdf[1])) + @test_broken all(ret[cdf[i] + 1 : cdf[i + 1]] .== (cdf[i + 1] - cdf[1]) * 99.0 / (cdf[end] - cdf[1])) end for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64}) @@ -162,7 +166,7 @@ eye(m,n) = Matrix{Float64}(I,m,n) img = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10]) img = T.(img) - _, counts_before = build_histogram(img,32,0,1) + _, counts_before = build_histogram(img,32; minval=0, maxval=1) nonzero_before = sum(counts_before .!= 0) #= @@ -170,8 +174,8 @@ eye(m,n) = Matrix{Float64}(I,m,n) and verify that all 32 bins have non-zero counts. This will confirm that the dynamic range of the original image has been increased. =# - imgeq = adjust_histogram(Equalization(),img,256,0,1) - edges, counts_after = build_histogram(imgeq,32,0,1) + imgeq = adjust_histogram(img, Equalization(256,0,1)) + edges, counts_after = build_histogram(imgeq,32; minval=0, maxval=1) nonzero_after = sum(counts_after .!= 0) @test nonzero_before < nonzero_after @test nonzero_after == 32 @@ -188,7 +192,7 @@ eye(m,n) = Matrix{Float64}(I,m,n) 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,0,1) + _, counts_before = build_histogram(img,32; minval=0, maxval=1) nonzero_before = sum(counts_before .!= 0) #= @@ -196,8 +200,8 @@ eye(m,n) = Matrix{Float64}(I,m,n) verify that all 32 bins have non-zero counts. This will confirm that the dynamic range of the original image has been increased. =# - imgeq = adjust_histogram(Equalization(),img,256,0,1) - edges, counts_after = build_histogram(imgeq,32,0,1) + imgeq = adjust_histogram(img, Equalization(256,0,1)) + edges, counts_after = build_histogram(imgeq,32; minval=0, maxval=1) nonzero_after = sum(counts_after .!= 0) @test nonzero_before < nonzero_after @test nonzero_after == 32 @@ -206,11 +210,11 @@ eye(m,n) = Matrix{Float64}(I,m,n) # Verify that the minimum and maximum values of the equalised image match the # specified minimum and maximum values, i.e. that the intensities of the equalised # image are in the interval [minvalue, maxvalue]. - imgeq = adjust_histogram(Equalization(),collect(0:1:255),256,64,128) + imgeq = adjust_histogram(collect(0:1:255), Equalization(256,64,128)) @test all(imgeq[1:65] .== 64) @test all(imgeq[128+1:end] .== 128) - imgeq = adjust_histogram(Equalization(),collect(0:1/255:1),256,64/255,128/255) + imgeq = adjust_histogram(collect(0:1/255:1), Equalization(256,64/255,128/255)) @test all(imgeq[1:65] .== 64/255) @test all(imgeq[128+1:end] .== 128/255) @@ -224,110 +228,97 @@ eye(m,n) = Matrix{Float64}(I,m,n) @testset "Gamma Correction" begin img = oneunits(Gray{Float64}, 10, 10) - ret = adjust_gamma(img, 1) + ret = adjust_histogram(img, GammaCorrection(; gamma=1)) @test img == ret @test eltype(ret) == eltype(img) imgp = padarray(img, Fill(0, (2,2))) - retp = adjust_gamma(imgp, 1) + retp = adjust_histogram(imgp, GammaCorrection(; gamma=1)) @test imgp == retp @test eltype(retp) == eltype(imgp) img = oneunits(Gray{N0f8}, 10, 10) - ret = adjust_gamma(img, 1) + ret = adjust_histogram(img, GammaCorrection(; gamma=1)) @test img == ret @test eltype(ret) == eltype(img) imgp = padarray(img, Fill(0, (2,2))) - retp = adjust_gamma(imgp, 1) + retp = adjust_histogram(imgp, GammaCorrection(; gamma=1)) @test imgp == retp @test eltype(retp) == eltype(imgp) img = oneunits(Gray{N0f16}, 10, 10) - ret = adjust_gamma(img, 1) + ret = adjust_histogram(img, GammaCorrection(; gamma=1)) @test eltype(ret) == eltype(img) imgp = padarray(img, Fill(0, (2,2))) - retp = adjust_gamma(imgp, 1) + retp = adjust_histogram(imgp, GammaCorrection(; gamma=1)) @test imgp == retp @test eltype(retp) == eltype(imgp) - img = oneunits(AGray{N0f8}, 10, 10) - ret = adjust_gamma(img, 1) - @test img == ret - @test eltype(ret) == eltype(img) + @info "Test skipped due to lack of support for `AGray` in adjust_histogram" + # img = oneunits(AGray{N0f8}, 10, 10) + # ret = adjust_histogram(img, GammaCorrection(; gamma=1)) + # @test img == ret + # @test eltype(ret) == eltype(img) imgp = padarray(img, Fill(0, (2,2))) - retp = adjust_gamma(imgp, 1) + retp = adjust_histogram(imgp, GammaCorrection(; gamma=1)) @test imgp == retp @test eltype(retp) == eltype(imgp) img = oneunits(RGB{N0f8}, 10, 10) - ret = adjust_gamma(img, 1) + ret = adjust_histogram(img, GammaCorrection(; gamma=1)) @test img == ret @test eltype(ret) == eltype(img) imgp = padarray(img, Fill(zero(eltype(img)), (2,2))) - retp = adjust_gamma(imgp, 1) + retp = adjust_histogram(imgp, GammaCorrection(; gamma=1)) @test imgp == retp @test eltype(retp) == eltype(imgp) img = oneunits(RGB{N0f16}, 10, 10) - ret = adjust_gamma(img, 1) + ret = adjust_histogram(img, GammaCorrection(; gamma=1)) @test img == ret @test eltype(ret) == eltype(img) imgp = padarray(img, Fill(zero(eltype(img)), (2,2))) - retp = adjust_gamma(imgp, 1) + retp = adjust_histogram(imgp, GammaCorrection(; gamma=1)) @test imgp == retp @test eltype(retp) == eltype(imgp) img = oneunits(RGB{Float64}, 10, 10) - ret = adjust_gamma(img, 1) + ret = adjust_histogram(img, GammaCorrection(; gamma=1)) @test all(map((i, r) -> isapprox(i, r), img, ret)) @test eltype(ret) == eltype(img) imgp = padarray(img, Fill(zero(eltype(img)), (2,2))) - retp = adjust_gamma(imgp, 1) + retp = adjust_histogram(imgp, GammaCorrection(; gamma=1)) @test all(map((i, r) -> isapprox(i, r), imgp, retp)) @test eltype(retp) == eltype(imgp) - img = oneunits(ARGB{N0f8}, 10, 10) - ret = adjust_gamma(img, 1) - @test img == ret + @info "Test skipped due to lack of support for `ARGB` in adjust_histogram" + # img = oneunits(ARGB{N0f8}, 10, 10) + # ret = adjust_histogram(img, GammaCorrection(; gamma=1)) + # @test img == ret imgp = padarray(img, Fill(zero(eltype(img)), (2,2))) - retp = adjust_gamma(imgp, 1) - @test imgp == retp - - - #Working + retp = adjust_histogram(imgp, GammaCorrection(; gamma=1)) + @test_broken imgp == retp img = reshape(1:1:100, 10, 10) - ret = adjust_gamma(img, 2) + ret = adjust_histogram(img, GammaCorrection(; gamma=2)) @test ret == img .^ 2 - ret = adjust_gamma(img, 0.3) - @test ret == img .^ 0.3 - ret = adjust_gamma(img, 1, 1, 100) - @test ret == img + ret = adjust_histogram(img, GammaCorrection(; gamma=0.3)) + @test_broken ret == img .^ 0.3 + @info "Test skipped due to lack of support for `minval` and `maxval` in `GammaCorrection`" + # ret = adjust_histogram(img, GammaCorrection(; gamma=1, minval=1, maxval=100)) + # @test ret == img img = zeros(10, 10) - ret = adjust_gamma(img, 2) + ret = adjust_histogram(img, GammaCorrection(; gamma=2)) @test all(ret .== 0) - a = ARGB(0.2, 0.3, 0.4, 0.5) - r = Images._gamma_pixel_rescale(a, 1) - @test isapprox(a, r, rtol = 0.0001) - r = Images._gamma_pixel_rescale(a, 2) - @test alpha(r) == alpha(a) - - b = AGray(0.2, 0.6) - r = Images._gamma_pixel_rescale(b, 1) - @test b == r - r = Images._gamma_pixel_rescale(b, 2) - @test alpha(r) == alpha(b) - @test isapprox(r.val, b.val ^ 2) - # Verify exported ImageContrastAdjustment symbol img = Gray{N0f8}.(rand(10, 10)); @inferred adjust_histogram(img, GammaCorrection(0.1)) @@ -339,47 +330,49 @@ eye(m,n) = Matrix{Float64}(I,m,n) #DataTypes img = oneunits(Gray{Float64}, 10, 10) - ret = histmatch(img, img) - @test all(ret .== zero(eltype(img))) + ret = adjust_histogram(img, Matching(targetimg=img)) + @test_broken all(ret .== zero(eltype(img))) @test eltype(ret) == eltype(img) img = oneunits(Gray{N0f8}, 10, 10) - ret = histmatch(img, img) - @test all(ret .== zero(eltype(img))) + ret = adjust_histogram(img, Matching(targetimg=img)) + @test_broken all(ret .== zero(eltype(img))) @test eltype(ret) == eltype(img) img = oneunits(Gray{N0f16}, 10, 10) - ret = histmatch(img, img) + ret = adjust_histogram(img, Matching(targetimg=img)) @test eltype(ret) == eltype(img) - img = oneunits(AGray{N0f8}, 10, 10) - ret = histmatch(img, img) - @test eltype(ret) == eltype(img) + @info "Test skipped due to lack of support for `AGray` in adjust_histogram" + # img = oneunits(AGray{N0f8}, 10, 10) + # ret = adjust_histogram(img, Matching(targetimg=img)) + # @test eltype(ret) == eltype(img) img = oneunits(RGB{N0f8}, 10, 10) - ret = histmatch(img, img) + ret = adjust_histogram(img, Matching(targetimg=img)) @test eltype(ret) == eltype(img) img = oneunits(RGB{N0f16}, 10, 10) - ret = histmatch(img, img) + ret = adjust_histogram(img, Matching(targetimg=img)) @test eltype(ret) == eltype(img) img = oneunits(RGB{Float64}, 10, 10) - ret = histmatch(img, img) - @test all(map((i, r) -> isapprox(zero(RGB), r, atol = 0.001), img, ret)) + ret = adjust_histogram(img, Matching(targetimg=img)) + @test_broken all(map((i, r) -> isapprox(zero(RGB), r, atol = 0.001), img, ret)) @test eltype(ret) == eltype(img) - img = oneunits(ARGB{N0f8}, 10, 10) - ret = histmatch(img, img) - @test eltype(ret) == eltype(img) + @info "Test skipped due to lack of support for `ARGB` in adjust_histogram" + # img = oneunits(ARGB{N0f8}, 10, 10) + # ret = adjust_histogram(img, Matching(targetimg=img)) + # @test eltype(ret) == eltype(img) img = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] - edges, hist = imhist(img, 2) - himg = Images._histmatch(img, edges, hist) - @test himg == [0, 0, 0, 0, 0, 5, 5, 5, 5, 5] - edges, hist = imhist(img, 5) - himg = Images._histmatch(img, edges, hist) - @test himg == [0, 0, 2, 2, 4, 4, 6, 6, 8, 8] + edges, hist = build_histogram(img, 2) + himg = adjust_histogram(img, Matching(; targetimg=img, edges)) + @test_broken himg == [0, 0, 0, 0, 0, 5, 5, 5, 5, 5] + edges, hist = build_histogram(img, 5) + himg = adjust_histogram(img, Matching(; targetimg=img, edges)) + @test_broken himg == [0, 0, 2, 2, 4, 4, 6, 6, 8, 8] @test_throws ErrorException Matching() @@ -389,114 +382,55 @@ eye(m,n) = Matrix{Float64}(I,m,n) @inferred adjust_histogram!(img, Matching(targetimg = img)) end - @testset "CLAHE" begin + @testset "AdaptiveEqualization" begin #DataTypes img = oneunits(Gray{Float64}, 10, 10) - ret = clahe(img, 100) + ret = adjust_histogram(img, AdaptiveEqualization(nbins=100)) @test size(ret) == size(img) @test eltype(ret) == eltype(img) img = oneunits(Gray{N0f8}, 10, 10) - ret = clahe(img, 100) + ret = adjust_histogram(img, AdaptiveEqualization(nbins=100)) @test size(ret) == size(img) @test eltype(ret) == eltype(img) img = oneunits(Gray{N0f16}, 10, 10) - ret = clahe(img, 100) + ret = adjust_histogram(img, AdaptiveEqualization(nbins=100)) @test size(ret) == size(img) @test eltype(ret) == eltype(img) - img = oneunits(AGray{N0f8}, 10, 10) - ret = clahe(img, 100) - @test size(ret) == size(img) - @test eltype(ret) == eltype(img) + @info "Test skipped for lack of support for `AGray` in adjust_histogram" + # img = oneunits(AGray{N0f8}, 10, 10) + # ret = adjust_histogram(img, AdaptiveEqualization(nbins=100)) + # @test size(ret) == size(img) + # @test eltype(ret) == eltype(img) img = oneunits(RGB{N0f8}, 10, 10) - ret = clahe(img, 100) + ret = adjust_histogram(img, AdaptiveEqualization(nbins=100)) @test size(ret) == size(img) @test eltype(ret) == eltype(img) img = oneunits(RGB{N0f16}, 10, 10) - ret = clahe(img, 100) + ret = adjust_histogram(img, AdaptiveEqualization(nbins=100)) @test size(ret) == size(img) @test eltype(ret) == eltype(img) img = oneunits(RGB{Float64}, 10, 10) - ret = clahe(img, 100) - @test size(ret) == size(img) - @test eltype(ret) == eltype(img) - - img = oneunits(ARGB{N0f8}, 10, 10) - ret = clahe(img, 100) + ret = adjust_histogram(img, AdaptiveEqualization(nbins=100)) @test size(ret) == size(img) @test eltype(ret) == eltype(img) - #Working + @info "Test skipped for lack of support for `ARGB` in adjust_histogram" + # img = oneunits(ARGB{N0f8}, 10, 10) + # ret = adjust_histogram(img, AdaptiveEqualization(nbins=100)) + # @test size(ret) == size(img) + # @test eltype(ret) == eltype(img) - cdf = Array(0.1:0.1:1.0) - res_pix = Images._clahe_pixel_rescale(0.5654, cdf, cdf, cdf, cdf, 0.1:0.1:1.0, 1, 4, 10, 10) - @test isapprox(res_pix, 0.5) - res_pix = Images._clahe_pixel_rescale(0.2344, cdf, cdf, cdf, cdf, 0.1:0.1:1.0, 1, 4, 10, 10) - @test isapprox(res_pix, 0.2) - res_pix = Images._clahe_pixel_rescale(0.8123, cdf, cdf, cdf, cdf, 0.1:0.1:1.0, 1, 4, 10, 10) - @test isapprox(res_pix, 0.8) - cdf2 = [0.0, 0.0, 0.1, 0.1, 0.2, 0.3, 0.5, 0.6, 0.8, 1.0] - - res_pix = Images._clahe_pixel_rescale(0.3, cdf2, cdf, 0.1:0.1:1, 3, 10) - @test isapprox(res_pix, 0.144444444) - - res_pix = Images._clahe_pixel_rescale(0.3, cdf2, cdf, 0.1:0.1:1, 6, 10) - @test isapprox(res_pix, 0.211111111) - - res_pix = Images._clahe_pixel_rescale(0.5654, cdf2, 0.1:0.1:1.0) - @test isapprox(res_pix, 0.2) - - cdf3 = [0.0, 0.0, 0.1, 0.3, 0.5, 0.5, 0.6, 0.9, 1.0, 1.0] - cdf4 = [0.0, 0.1, 0.3, 0.4, 0.5, 0.6, 0.6, 0.6, 0.8, 1.0] - - res_pix = Images._clahe_pixel_rescale(0.8123, cdf, cdf2, cdf3, cdf4, 0.1:0.1:1.0, 3, 4, 10, 10) - @test isapprox(res_pix, 0.781481481) - res_pix = Images._clahe_pixel_rescale(0.3, cdf, cdf2, cdf3, cdf4, 0.1:0.1:1.0, 5, 4, 10, 10) - @test isapprox(res_pix, 0.2037037037) - - img = [zeros(4,4) 0.1*ones(4,4); 0.2*ones(4,4) 0.3*ones(4,4)] - hist_eq_img = clahe(img, xblocks = 2, yblocks = 2, clip = 5) - # TODO: it would be good to understand where these numbers come from - @test all(x->x ≈ 0.2960140679953115, hist_eq_img[1:3, 1:3]) - @test all(x->x ≈ 0.4196951934349368, hist_eq_img[1:3, end-2:end]) - @test all(x->x ≈ 0.4841735052754997, hist_eq_img[end-2:end, 1:3]) - @test all(x->x ≈ 0.5486518171160624, hist_eq_img[end-2:end, end-2:end]) # TODO: the "boundary regions" img[4:5, :] and img[:, 4:5] are not # "symmetric," is this a problem? - # The following test is sensitive to any change in the - # algorithm. Is there way to make this more generic? - # img = [ 0.907976 0.0601377 0.747873 0.938683 0.935189 0.49492 0.508545 0.394573 0.991136 0.210218 - # 0.219737 0.727806 0.606995 0.25272 0.279777 0.584529 0.708244 0.958645 0.979217 0.418678 - # 0.793384 0.314803 0.324222 0.840858 0.49438 0.102003 0.724399 0.133663 0.312085 0.662211 - # 0.905336 0.423511 0.564616 0.692544 0.367656 0.630386 0.133439 0.529039 0.0175596 0.644065 - # 0.255092 0.00242822 0.723724 0.673323 0.244508 0.518068 0.204353 0.34222 0.106092 0.0331161 - # 0.432963 0.383491 0.903465 0.579605 0.719874 0.571533 0.728544 0.1864 0.950187 0.470226 - # 0.877475 0.554114 0.987133 0.947148 0.710115 0.131948 0.711611 0.0221843 0.470008 0.952806 - # 0.231911 0.177463 0.742054 0.0333307 0.481319 0.716638 0.332057 0.978177 0.0610481 0.439462 - # 0.935457 0.159602 0.178357 0.163585 0.275052 0.0557963 0.066368 0.199349 0.9238 0.85161 - # 0.692148 0.503563 0.0918827 0.0206237 0.702344 0.546088 0.04163 0.174103 0.45499 0.90019 - # ] - # hist_eq_img = clahe(img, xblocks = 2, yblocks = 2, clip = 5) - # ret_expected = [ 0.972222 0.154261 0.672734 0.962918 0.960518 0.603715 0.559332 0.684868 0.919632 0.11213 - # 0.26227 0.577172 0.615628 0.420623 0.356921 0.581208 0.845865 0.888908 0.908469 0.518953 - # 0.79315 0.398521 0.454306 0.692093 0.404984 0.287649 0.619472 0.47354 0.442781 0.826681 - # 0.875951 0.325109 0.561585 0.746554 0.355356 0.542615 0.308523 0.303799 0.12818 0.662404 - # 0.218908 0.145377 0.627466 0.703423 0.255746 0.543485 0.270196 0.240487 0.119202 0.0762326 - # 0.443491 0.419296 0.857736 0.778583 0.78256 0.660936 0.66722 0.420796 0.793904 0.50695 - # 0.765105 0.583873 0.870947 0.825773 0.755839 0.346564 0.598017 0.329962 0.564487 0.897893 - # 0.512372 0.328346 0.50314 0.339795 0.493125 0.476082 0.45704 0.597256 0.385295 0.814294 - # 0.866771 0.256846 0.140199 0.118129 0.287474 0.175723 0.115673 0.343809 0.735832 0.928661 - # 0.833503 0.481259 0.135795 0.116606 0.737954 0.57691 0.103692 0.132731 0.539984 0.999525 - # ] - # @test all(map((i, j) -> isapprox(i, j, rtol = 0.00001), hist_eq_img, ret_expected)) # Verify exported ImageContrastAdjustment symbol img = Gray{N0f8}.(rand(10, 10)); @inferred adjust_histogram(img, AdaptiveEqualization()) @@ -506,35 +440,17 @@ eye(m,n) = Matrix{Float64}(I,m,n) @testset "Other" begin # Issue #282 - img = Gray{N0f8}.(eye(2,2)) - imgs = imstretch(img, 0.3, 0.4) - @test imgs ≈ 1 ./ (1 .+ (0.3 ./ (eye(2,2) .+ eps(float(N0f8)))).^0.4) + rawimg = [1 0; 0 1] + img = Gray{N0f8}.(rawimg) + imgs = adjust_histogram(img, ContrastStretching(t=0.3, slope=0.4, ϵ=eps(Float32))) + @test imgs ≈ N0f8.(1 ./ (1 .+ (0.3 ./ (rawimg .+ eps(float(N0f8)))).^0.4)) img = Gray{N0f16}.([0.01164 0.01118; 0.01036 0.01187]) - @test imadjustintensity(img,[0.0103761, 0.0252166])[2,1] == 0.0 - @test eltype(imadjustintensity(img)) == Gray{N0f16} + @test adjust_histogram(img, LinearStretching(; src_minval=0.0103761, src_maxval=0.0252166))[2,1] == 0.0 + @test eltype(adjust_histogram(img, LinearStretching())) == Gray{N0f16} hist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] - clipped_hist = cliphist(hist, 2) - @test all(clipped_hist .== 2.0) - clipped_hist = cliphist(hist, 3) - @test all(clipped_hist .== 3.0) - clipped_hist = cliphist(hist, 4) - @test all(clipped_hist .== 4.0) - clipped_hist = cliphist(hist, 5) - @test all(clipped_hist .== 5.0) - clipped_hist = cliphist(hist, 6) - @test clipped_hist == vec([3, 4, 6, 6, 6, 6, 6, 6, 6, 6]) - clipped_hist = cliphist(hist, 7) - @test clipped_hist == vec([2.6, 3.6, 3.6, 4.6, 5.6, 7, 7, 7, 7, 7]) - clipped_hist = cliphist(hist, 8) - @test clipped_hist == vec([2.3, 2.3, 3.3, 4.3, 5.3, 6.3, 7.3, 8, 8, 8]) - clipped_hist = cliphist(hist, 9) - @test clipped_hist == vec([2.1, 2.1, 3.1, 4.1, 5.1, 6.1, 7.1, 8.1, 9, 9]) - clipped_hist = cliphist(hist, 10) - @test clipped_hist == hist - # Verify exported ImageContrastAdjustment symbol img = Gray{N0f8}.(rand(10, 10)); @inferred adjust_histogram(img, ContrastStretching()) diff --git a/test/runtests.jl b/test/runtests.jl index 07a57037..1ecd57b9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,16 +4,29 @@ using Test using ImageBase using ImageBase.OffsetArrays using Suppressor +using Aqua -@testset "Images" begin -include("arrays.jl") -include("algorithms.jl") -@suppress_err include("exposure.jl") # deprecated -include("edge.jl") -include("corner.jl") -include("writemime.jl") +@testset "Aqua" begin + Aqua.test_all(Images; + ambiguities=false, # TODO? fix ambiguities? (may not be easy/possible) Currently 14 ambiguities + undefined_exports=false, # TODO: remove `data` export from ImageMetadata + stale_deps=false, # ImageMagick and ImageIO are not loaded + piracy=false, # TODO: fix piracy of `zero(Point)` in edges.jl + ) +end -@suppress_err include("legacy.jl") -@suppress_err include("deprecated.jl") +@testset "Images" begin + include("arrays.jl") + include("algorithms.jl") + @suppress_err include("exposure.jl") # deprecated + include("edge.jl") + include("corner.jl") + include("writemime.jl") + # @suppress_err include("legacy.jl") + if Base.JLOptions().depwarn < 2 + @suppress_err include("deprecated.jl") + else + @info "Skipping deprecated tests because of --depwarn=error" + end end