From 802c4064655ad21ec62560e4ae5eb675879e040b Mon Sep 17 00:00:00 2001 From: cpaniaguam Date: Wed, 6 Nov 2024 12:44:36 -0500 Subject: [PATCH 1/3] feat: histeq for wateshed workflow --- src/histogram_equalization.jl | 95 +++++++++++++++++------- test/test-conditional-adaptive-histeq.jl | 57 +++++++++----- 2 files changed, 106 insertions(+), 46 deletions(-) diff --git a/src/histogram_equalization.jl b/src/histogram_equalization.jl index 5ade5ab3..2b6d6c28 100644 --- a/src/histogram_equalization.jl +++ b/src/histogram_equalization.jl @@ -4,7 +4,6 @@ function to_uint8(arr::AbstractMatrix{T}) where {T<:AbstractFloat} return img end - function anisotropic_diffusion_3D(I) rgbchannels = get_rgb_channels(I) @@ -13,10 +12,11 @@ function anisotropic_diffusion_3D(I) end return rgbchannels - end -function anisotropic_diffusion_2D(I::AbstractMatrix{T}; gradient_threshold::Union{T,Nothing}=nothing, niter::Int=1) where {T} +function anisotropic_diffusion_2D( + I::AbstractMatrix{T}; gradient_threshold::Union{T,Nothing}=nothing, niter::Int=1 +) where {T} if eltype(I) <: Int I = Gray.(I ./ 255) end @@ -34,11 +34,13 @@ function anisotropic_diffusion_2D(I::AbstractMatrix{T}; gradient_threshold::Unio for _ in 1:niter # These are zero-indexed offset arrays - diff_img_north = padded_img[0:end-1, 1:end-1] .- padded_img[1:end, 1:end-1] - diff_img_east = padded_img[1:end-1, 1:end] .- padded_img[1:end-1, 0:end-1] - diff_img_nw = padded_img[0:end-2, 0:end-2] .- I - diff_img_ne = padded_img[0:end-2, 2:end] .- I - diff_img_sw = padded_img[2:end, 0:end-2] .- I + diff_img_north = + padded_img[0:(end - 1), 1:(end - 1)] .- padded_img[1:end, 1:(end - 1)] + diff_img_east = + padded_img[1:(end - 1), 1:end] .- padded_img[1:(end - 1), 0:(end - 1)] + diff_img_nw = padded_img[0:(end - 2), 0:(end - 2)] .- I + diff_img_ne = padded_img[0:(end - 2), 2:end] .- I + diff_img_sw = padded_img[2:end, 0:(end - 2)] .- I diff_img_se = padded_img[2:end, 2:end] .- I # Exponential conduction coefficients @@ -49,7 +51,6 @@ function anisotropic_diffusion_2D(I::AbstractMatrix{T}; gradient_threshold::Unio conduct_coeff_sw = exp.(-(abs.(diff_img_sw) ./ gradient_threshold) .^ 2) conduct_coeff_se = exp.(-(abs.(diff_img_se) ./ gradient_threshold) .^ 2) - # Flux calculations flux_north = conduct_coeff_north .* diff_img_north flux_east = conduct_coeff_east .* diff_img_east @@ -59,28 +60,24 @@ function anisotropic_diffusion_2D(I::AbstractMatrix{T}; gradient_threshold::Unio flux_se = conduct_coeff_se .* diff_img_se # Back to regular 1-indexed arrays - flux_north_diff = flux_north[1:end-1, :] .- flux_north[2:end, :] - flux_east_diff = flux_east[:, 2:end] .- flux_east[:, 1:end-1] + flux_north_diff = flux_north[1:(end - 1), :] .- flux_north[2:end, :] + flux_east_diff = flux_east[:, 2:end] .- flux_east[:, 1:(end - 1)] # Discrete PDE solution sum_ = (1 / (dd^2)) .* (flux_nw .+ flux_ne .+ flux_sw .+ flux_se) I = I .+ diffusion_rate .* (flux_north_diff .- flux_north_diff .+ sum_) - end return I end - function imshow(img) if typeof(img) <: BitMatrix return Gray.(img) end - Gray.(img ./ 255) + return Gray.(img ./ 255) end - - function adapthisteq(img::Matrix{T}, nbins=256, clip=0.01) where {T} # Step 1: Normalize the image to [0, 1] based on its own min and max image_min, image_max = minimum(img), maximum(img) @@ -88,13 +85,15 @@ function adapthisteq(img::Matrix{T}, nbins=256, clip=0.01) where {T} # Step 2: Apply adaptive histogram equalization. equalize_adapthist handles the tiling to 1/8 of the image size (equivalent to 8x8 blocks in MATLAB) equalized_image = sk_exposure.equalize_adapthist( - normalized_image, + normalized_image; clip_limit=clip, # Equivalent to MATLAB's 'ClipLimit' - nbins=nbins # Number of histogram bins. 255 is used to match the default in MATLAB script + nbins=nbins, # Number of histogram bins. 255 is used to match the default in MATLAB script ) # Step 3: Rescale the image back to the original range [image_min, image_max] - final_image = sk_exposure.rescale_intensity(equalized_image, in_range="image", out_range=(image_min, image_max)) + final_image = sk_exposure.rescale_intensity( + equalized_image; in_range="image", out_range=(image_min, image_max) + ) # Convert back to the original data type if necessary final_image = to_uint8(final_image) @@ -120,10 +119,9 @@ function get_rgb_channels(img) greenc = green.(img) * 255 bluec = blue.(img) * 255 - return cat(redc, greenc, bluec, dims=3) + return cat(redc, greenc, bluec; dims=3) end - function _process_image_tiles( true_color_image, clouds_red, @@ -156,7 +154,6 @@ function _process_image_tiles( return rgbchannels end - """ conditional_histeq( true_color_image, @@ -192,7 +189,7 @@ function conditional_histeq( white_threshold::AbstractFloat=25.5, white_fraction_threshold::AbstractFloat=0.4, ) - tiles = get_tiles(true_color_image, rblocks=rblocks, cblocks=cblocks) + tiles = get_tiles(true_color_image; rblocks=rblocks, cblocks=cblocks) rgbchannels_equalized = _process_image_tiles( true_color_image, clouds_red, @@ -203,7 +200,6 @@ function conditional_histeq( ) return rgbchannels_equalized - end """ @@ -227,7 +223,6 @@ function conditional_histeq( white_threshold::AbstractFloat=25.5, white_fraction_threshold::AbstractFloat=0.4, ) - side_length = IceFloeTracker.get_optimal_tile_size(side_length, size(true_color_image)) tiles = IceFloeTracker.get_tiles(true_color_image, side_length) @@ -242,16 +237,16 @@ function conditional_histeq( ) return rgbchannels_equalized - end -function _get_false_color_cloudmasked(; false_color_image, +function _get_false_color_cloudmasked(; + false_color_image, prelim_threshold=110.0, band_7_threshold=200.0, band_2_threshold=190.0, ) mask_cloud_ice, clouds_view = IceFloeTracker._get_masks( - false_color_image, + false_color_image; prelim_threshold=prelim_threshold, band_7_threshold=band_7_threshold, band_2_threshold=band_2_threshold, @@ -279,3 +274,47 @@ Convert an RGB image to grayscale in the range [0, 255]. function rgb2gray(img::Matrix{RGB{Float64}}) return round.(Int, Gray.(img) * 255) end + +function _imhist(img, rng) + d = Dict(k => 0 for k in rng) + for i in img + d[i] = d[i] + 1 + end + k, heights = collect.([keys(d), values(d)]) + order = sortperm(k) + k, heights = k[order], heights[order] + return k, heights +end + +function imhist(img, imgtype::AbstractString="uint8") + + # TODO: add validation for arr: either uint8 0:255 or grayscale 0:1 + rng = imgtype == "uint8" ? range(0, 255) : range(0; stop=1, length=256) + # use range(0, stop=1, length=256) for grayscale images + + return _imhist(img, rng) +end + +""" + imhistp(img) +Parsimonious version of `imhist` that uses the maximum value of the image to determine the number of bins. +""" +function imhistp(img) + mx = maximum(img) + nbins = 2^Int(ceil(log2(mx))) + return _imhist(img, 0:(nbins - 1)) +end + +""" + histeq(img) +Histogram equalization of `img` according to [1]. +[1] R. C. Gonzalez and R. E. Woods. Digital Image Processing (3rd Edition). Upper Saddle River, NJ, USA: Prentice-Hall, 2006. +""" +function histeq(img::S)::S where {S<:AbstractArray{<:Integer}} + k, heights = imhistp(img) + cdf = last(k) * cumsum(heights) / sum(heights) + s = round.(Int, cdf, RoundNearestTiesAway) + mapping = Dict(zip(k, s)) + T(r) = mapping[r] + return T.(img) +end diff --git a/test/test-conditional-adaptive-histeq.jl b/test/test-conditional-adaptive-histeq.jl index d571ab6b..e58172a9 100644 --- a/test/test-conditional-adaptive-histeq.jl +++ b/test/test-conditional-adaptive-histeq.jl @@ -1,9 +1,20 @@ -using IceFloeTracker: _get_false_color_cloudmasked, convert_to_255_matrix, adapthisteq, conditional_histeq, rgb2gray, to_uint8 +using IceFloeTracker: + _get_false_color_cloudmasked, + convert_to_255_matrix, + adapthisteq, + conditional_histeq, + rgb2gray, + to_uint8, + histeq begin datadir = joinpath(@__DIR__, "test_inputs/") - path_true_color_image = joinpath(datadir, "NE_Greenland_truecolor.2020162.aqua.250m.tiff") - path_false_color_image = joinpath(datadir, "NE_Greenland_reflectance.2020162.aqua.250m.tiff") + path_true_color_image = joinpath( + datadir, "NE_Greenland_truecolor.2020162.aqua.250m.tiff" + ) + path_false_color_image = joinpath( + datadir, "NE_Greenland_reflectance.2020162.aqua.250m.tiff" + ) true_color_image = float64.(load(path_true_color_image)) false_color_image = float64.(load(path_false_color_image)) dilated_landmask = BitMatrix(load(joinpath(datadir, "matlab_landmask.png"))) @@ -11,7 +22,7 @@ end function test_cloud_image_workflow() @testset "Prereq cloud image" begin - false_color_cloudmasked = _get_false_color_cloudmasked( + false_color_cloudmasked = _get_false_color_cloudmasked(; false_color_image=false_color_image, prelim_threshold=110.0, band_7_threshold=200.0, @@ -22,7 +33,6 @@ function test_cloud_image_workflow() end end - function test_adaphisteq() @testset "Adaptive histogram equalization" begin img = convert_to_255_matrix(testimage("cameraman")) @@ -31,10 +41,9 @@ function test_adaphisteq() end end - function test_conditional_adaptivehisteq() @testset "Conditional adaptivehisteq" begin - clouds = _get_false_color_cloudmasked( + clouds = _get_false_color_cloudmasked(; false_color_image=false_color_image, prelim_threshold=110.0, band_7_threshold=200.0, @@ -47,29 +56,20 @@ function test_conditional_adaptivehisteq() @test sum(clouds_red) == 1_320_925_065 # Using rblocks = 8, cblocks = 6 - true_color_eq = conditional_histeq( - true_color_image, - clouds_red, - 8, - 6) + true_color_eq = conditional_histeq(true_color_image, clouds_red, 8, 6) # This differs from MATLAB script due to disparity in the implementations # of the adaptive histogram equalization / diffusion functions # For the moment testing for regression @test sum(to_uint8(true_color_eq[:, :, 1])) == 6_372_159_606 - # Use custom tile size side_length = size(true_color_eq, 1) รท 8 - true_color_eq = conditional_histeq( - true_color_image, - clouds_red, - side_length) + true_color_eq = conditional_histeq(true_color_image, clouds_red, side_length) @test sum(to_uint8(true_color_eq[:, :, 1])) == 6_328_796_398 end end - function test_rgb2gray() @testset "RGB to grayscale" begin g = rgb2gray(true_color_image) @@ -77,9 +77,30 @@ function test_rgb2gray() end end +function test_histeq() + @testset "histeq" begin + img = [ + 4 4 4 4 4 + 3 4 5 4 3 + 3 5 5 5 3 + 3 4 5 4 3 + 4 4 4 4 4 + ] + expected = [ + 6 6 6 6 6 + 2 6 7 6 2 + 2 7 7 7 2 + 2 6 7 6 2 + 6 6 6 6 6 + ] + @test histeq(img) == expected + end +end + @testset "Conditional adaptivehisteq" begin test_cloud_image_workflow() test_adaphisteq() test_conditional_adaptivehisteq() test_rgb2gray() + test_histeq() end From 011b80cc1f26e92660ac0fb68c981821050a3fa7 Mon Sep 17 00:00:00 2001 From: Carlos Paniagua Date: Fri, 8 Nov 2024 16:04:22 -0500 Subject: [PATCH 2/3] Update histogram_equalization.jl fix: bump to next power of 2 even for powers of 2 --- src/histogram_equalization.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/histogram_equalization.jl b/src/histogram_equalization.jl index 2b6d6c28..64f59bee 100644 --- a/src/histogram_equalization.jl +++ b/src/histogram_equalization.jl @@ -300,11 +300,11 @@ end Parsimonious version of `imhist` that uses the maximum value of the image to determine the number of bins. """ function imhistp(img) - mx = maximum(img) - nbins = 2^Int(ceil(log2(mx))) - return _imhist(img, 0:(nbins - 1)) + nbins = nextpow(2, maximum(img) + 1) + return _imhist(img, 0:(nbins-1)) end + """ histeq(img) Histogram equalization of `img` according to [1]. From 404ad1ad3b252144a7998a765cc557504c1fe01b Mon Sep 17 00:00:00 2001 From: Carlos Paniagua Date: Fri, 8 Nov 2024 16:44:50 -0500 Subject: [PATCH 3/3] test: add edge cases --- test/test-conditional-adaptive-histeq.jl | 48 ++++++++++++++++++------ 1 file changed, 36 insertions(+), 12 deletions(-) diff --git a/test/test-conditional-adaptive-histeq.jl b/test/test-conditional-adaptive-histeq.jl index e58172a9..748f2f20 100644 --- a/test/test-conditional-adaptive-histeq.jl +++ b/test/test-conditional-adaptive-histeq.jl @@ -79,21 +79,45 @@ end function test_histeq() @testset "histeq" begin - img = [ - 4 4 4 4 4 - 3 4 5 4 3 - 3 5 5 5 3 - 3 4 5 4 3 - 4 4 4 4 4 + imgs = [ + [ + 4 4 4 4 4 + 3 4 5 4 3 + 3 5 5 5 3 + 3 4 5 4 3 + 4 4 4 4 4 + ], + # Edge cases + # maximum a power of 2 + [ + 0 0 0 + 2 2 2 + ], + # maximum at 255 + [ + 0 0 0 + 255 255 255 + ], ] expected = [ - 6 6 6 6 6 - 2 6 7 6 2 - 2 7 7 7 2 - 2 6 7 6 2 - 6 6 6 6 6 + [ + 6 6 6 6 6 + 2 6 7 6 2 + 2 7 7 7 2 + 2 6 7 6 2 + 6 6 6 6 6 + ], + [ + 2 2 2 + 3 3 3 + ], + [ + 128 128 128 + 255 255 255 + ], ] - @test histeq(img) == expected + + @test all(IceFloeTracker.histeq(imgs[i]) == expected[i] for i in 1:3) end end