Merge branch 'main' into deps-bump-pyproj-3_7

src/IceFloeTracker.jl

@@ -38,6 +38,7 @@ export readdlm,
     imsharpen,
     label_components,
     regionprops_table,
+    cropfloe,
     loadimg,
     matchcorr,
     centered,
@@ -70,11 +71,10 @@ include("branch.jl")
 include("special_strels.jl")
 include("tilingutils.jl")
 include("histogram_equalization.jl")
-
-
-const sk_measure = PyNULL()
-const sk_exposure = PyNULL()
-const getlatlon = PyNULL()
+include("brighten.jl")
+include("morph_fill.jl")
+include("imcomplement.jl")
+include("imadjust.jl")
 
 function get_version_from_toml(pth=dirname(dirname(pathof(IceFloeTracker))))::VersionNumber
     toml = TOML.parsefile(joinpath(pth, "Project.toml"))
@@ -83,9 +83,16 @@ end
 
 const IFTVERSION = get_version_from_toml()
 
+const sk_measure = PyNULL()
+const sk_morphology = PyNULL()
+const sk_exposure = PyNULL()
+const getlatlon = PyNULL()
+
 function __init__()
-    copy!(sk_measure, pyimport_conda("skimage.measure", "scikit-image=0.24.0"))
-    copy!(sk_exposure, pyimport_conda("skimage.exposure", "scikit-image=0.24.0"))
+    skimage = "scikit-image=0.24.0"
+    copy!(sk_measure, pyimport_conda("skimage.measure", skimage))
+    copy!(sk_exposure, pyimport_conda("skimage.exposure", skimage))
+    copy!(sk_morphology, pyimport_conda("skimage.morphology", skimage))
     pyimport_conda("pyproj", "pyproj=3.7.0")
     pyimport_conda("rasterio", "rasterio=1.4.2")
     pyimport_conda("jinja2", "jinja2=3.1.2")
@@ -152,34 +159,34 @@ julia> IceFloeTracker.MorphSE.dilate(a, se)
 ```
 """
 module MorphSE
-using ImageCore
-using ColorTypes
-using LoopVectorization
-using OffsetArrays
-using TiledIteration: EdgeIterator
-using DataStructures
-include("morphSE/StructuringElements.jl")
-using .StructuringElements
-include("morphSE/extreme_filter.jl")
-include("morphSE/utils.jl")
-include("morphSE/dilate.jl")
-include("morphSE/erode.jl")
-include("morphSE/opening.jl")
-include("morphSE/closing.jl")
-include("morphSE/bothat.jl")
-include("morphSE/mreconstruct.jl")
-include("morphSE/fill_holes.jl")
+    using ImageCore
+    using ColorTypes
+    using LoopVectorization
+    using OffsetArrays
+    using TiledIteration: EdgeIterator
+    using DataStructures
+    include("morphSE/StructuringElements.jl")
+    using .StructuringElements
+    include("morphSE/extreme_filter.jl")
+    include("morphSE/utils.jl")
+    include("morphSE/dilate.jl")
+    include("morphSE/erode.jl")
+    include("morphSE/opening.jl")
+    include("morphSE/closing.jl")
+    include("morphSE/bothat.jl")
+    include("morphSE/mreconstruct.jl")
+    include("morphSE/fill_holes.jl")
 end
 
 module Register
-include("Register/CenterIndexedArrays.jl-0.2.0/CenterIndexedArrays.jl")
-include("Register/RegisterCore.jl-0.2.4/src/RegisterCore.jl")
-include("Register/RegisterMismatchCommon.jl-master/RegisterMismatchCommon.jl")
-include("Register/RegisterUtilities.jl-master/RegisterUtilities.jl")
-include("Register/RFFT.jl-master/RFFT.jl")
-include("Register/RegisterDeformation.jl-0.4.4/RegisterDeformation.jl")
-include("Register/QuadDIRECT.jl-master/QuadDIRECT.jl")
-include("Register/RegisterQD.jl-0.3.1/RegisterQD.jl")
-include("Register/RegisterMismatch.jl-0.4.0/RegisterMismatch.jl")
+    include("Register/CenterIndexedArrays.jl-0.2.0/CenterIndexedArrays.jl")
+    include("Register/RegisterCore.jl-0.2.4/src/RegisterCore.jl")
+    include("Register/RegisterMismatchCommon.jl-master/RegisterMismatchCommon.jl")
+    include("Register/RegisterUtilities.jl-master/RegisterUtilities.jl")
+    include("Register/RFFT.jl-master/RFFT.jl")
+    include("Register/RegisterDeformation.jl-0.4.4/RegisterDeformation.jl")
+    include("Register/QuadDIRECT.jl-master/QuadDIRECT.jl")
+    include("Register/RegisterQD.jl-0.3.1/RegisterQD.jl")
+    include("Register/RegisterMismatch.jl-0.4.0/RegisterMismatch.jl")
 end
 end

src/bridge.jl

@@ -9,7 +9,8 @@ function _bridge_operator_lut(
     return _operator_lut(I, img, nhood, lutbridge)
 end
 
-function _bridge_filter(img::T, operator::Function)::T where {T<:AbstractArray{Bool}}
+# TODO: see about implemting _filter using parallelization
+function _filter(img::T, operator::Function)::T where {T<:AbstractArray{Bool}}
     out = zeros(Bool, size(img))
     R = CartesianIndices(img)
     I_first, I_last = first(R), last(R)
@@ -79,5 +80,5 @@ julia> bridge(bw)
 ```
 """
 function bridge(bw::T)::T where {T<:AbstractArray{Bool}}
-    return _bridge_filter(bw, _bridge_operator_lut)
+    return _filter(bw, _bridge_operator_lut)
 end

src/brighten.jl

@@ -0,0 +1,28 @@
+"""
+    get_brighten_mask(equalized_gray_reconstructed_img, gamma_green)
+# Arguments
+- `equalized_gray_reconstructed_img`: The equalized gray reconstructed image (uint8 in Matlab).
+- `gamma_green`: The gamma value for the green channel (also uint8).
+# Returns
+Difference equalized_gray_reconstructed_img - gamma_green clamped between 0 and 255.
+"""
+function get_brighten_mask(equalized_gray_reconstructed_img, gamma_green)
+    return to_uint8(equalized_gray_reconstructed_img - gamma_green)
+end
+
+"""
+    imbrighten(img, brighten_mask, bright_factor)
+Brighten the image using a mask and a brightening factor.
+# Arguments
+- `img`: The input image.
+- `brighten_mask`: A mask indicating the pixels to brighten.
+- `bright_factor`: The factor by which to brighten the pixels.
+# Returns
+- The brightened image.
+"""
+function imbrighten(img, brighten_mask, bright_factor)
+    img = Float64.(img)
+    brighten_mask = brighten_mask .> 0
+    img[brighten_mask] .= img[brighten_mask] * bright_factor
+    return img = to_uint8(img)
+end

src/bwtraceboundary.jl

@@ -54,7 +54,7 @@ julia> boundary[3]
 ```
 """
 function bwtraceboundary(
-    image::Union{Matrix{Int64},Matrix{Float64},T};
+    image::Union{Matrix{UInt8},Matrix{Int64},Matrix{Float64},T};
     P0::Union{Tuple{Int,Int},CartesianIndex{2},Nothing}=nothing,
     closed::Bool=true,
 ) where {T<:AbstractArray{Bool,2}}

src/histogram_equalization.jl

@@ -4,6 +4,10 @@ function to_uint8(arr::AbstractMatrix{T}) where {T<:AbstractFloat}
     return img
 end
 
+function to_uint8(arr::AbstractMatrix{T}) where {T<:Integer}
+    img = clamp.(arr, 0, 255)
+    return img
+end
 
 function anisotropic_diffusion_3D(I)
     rgbchannels = get_rgb_channels(I)
@@ -13,10 +17,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 +39,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 +56,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,42 +65,40 @@ 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)
     normalized_image = (img .- image_min) / (image_max - image_min)
 
     # 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 +124,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 +159,6 @@ function _process_image_tiles(
     return rgbchannels
 end
 
-
 """
     conditional_histeq(
     true_color_image,
@@ -192,7 +194,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 +205,6 @@ function conditional_histeq(
     )
 
     return rgbchannels_equalized
-
 end
 
 """
@@ -227,7 +228,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 +242,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 +279,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)
+    nbins = nextpow(2, maximum(img) + 1)
+    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