Tiled workflow feats

src/IceFloeTracker.jl

@@ -51,6 +51,36 @@ export DataFrames, DataFrame, nrow, Not, select!
 export Dates, Time, Date, DateTime, @dateformat_str
 export addlatlon!, getlatlon, convertcentroid!, converttounits!, dropcols!
 
+const sk_morphology = PyNULL()
+const sk_measure = PyNULL()
+const sk_exposure = PyNULL()
+const cv2 = PyNULL()
+const getlatlon = PyNULL()
+
+function __init__()
+    try
+        copy!(cv2, pyimport("cv2"))
+    catch e
+        pyimport("subprocess").run([
+            "python", "-m", "pip", "install", "opencv-python==4.10.0.84"
+        ])
+    end
+
+    # Import the installed modules
+
+    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.6.0")
+    pyimport_conda("rasterio", "rasterio=1.3.7")
+    pyimport_conda("jinja2", "jinja2=3.1.2")
+    pyimport_conda("pandas", "pandas=2")
+    @pyinclude(joinpath(@__DIR__, "latlon.py"))
+    copy!(getlatlon, py"getlatlon")
+    return nothing
+end
+
 include("utils.jl")
 include("persist.jl")
 include("landmask.jl")
@@ -70,31 +100,17 @@ include("branch.jl")
 include("special_strels.jl")
 include("tilingutils.jl")
 include("histogram_equalization.jl")
+include("morph_fill.jl")
+include("reconstruct.jl")
 
 
-const sk_measure = PyNULL()
-const sk_exposure = PyNULL()
-const getlatlon = PyNULL()
-
 function get_version_from_toml(pth=dirname(dirname(pathof(IceFloeTracker))))::VersionNumber
     toml = TOML.parsefile(joinpath(pth, "Project.toml"))
     return VersionNumber(toml["version"])
 end
 
 const IFTVERSION = get_version_from_toml()
 
-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"))
-    pyimport_conda("pyproj", "pyproj=3.6.0")
-    pyimport_conda("rasterio", "rasterio=1.3.7")
-    pyimport_conda("jinja2", "jinja2=3.1.2")
-    pyimport_conda("pandas", "pandas=2")
-    @pyinclude(joinpath(@__DIR__, "latlon.py"))
-    copy!(getlatlon, py"getlatlon")
-    return nothing
-end
-
 include("regionprops.jl")
 include("segmentation_a_direct.jl")
 include("segmentation_b.jl")

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/find_ice_labels.jl

@@ -22,10 +22,15 @@ function find_reflectance_peaks(
     possible_ice_threshold::Float64=Float64(75 / 255),
 )::Int64
     reflectance_channel[reflectance_channel .< possible_ice_threshold] .= 0 #75 / 255
+
     _, counts = ImageContrastAdjustment.build_histogram(reflectance_channel)
+
     locs, _ = Peaks.findmaxima(counts)
-    sort!(locs; rev=true)
-    return locs[2] # second greatest peak
+    Base.sort!(locs; rev=true)
+
+    length(locs) > 2 && return locs[2] # second greatest peak
+
+    return nothing
 end
 
 """
@@ -73,8 +78,12 @@ function find_ice_labels(
         if sum(abs.(ice_labels)) == 0
             ref_image_band_2 = @view(cv[2, :, :])
             ref_image_band_1 = @view(cv[3, :, :])
-            band_2_peak = find_reflectance_peaks(ref_image_band_2, possible_ice_threshold = possible_ice_threshold)
-            band_1_peak = find_reflectance_peaks(ref_image_band_1, possible_ice_threshold = possible_ice_threshold)
+            band_2_peak = find_reflectance_peaks(
+                ref_image_band_2; possible_ice_threshold=possible_ice_threshold
+            )
+            band_1_peak = find_reflectance_peaks(
+                ref_image_band_1; possible_ice_threshold=possible_ice_threshold
+            )
             mask_ice_band_2 = @view(cv[2, :, :]) .> band_2_peak / 255
             mask_ice_band_1 = @view(cv[3, :, :]) .> band_1_peak / 255
             ice = mask_ice_band_7 .* mask_ice_band_2 .* mask_ice_band_1

src/histogram_equalization.jl

@@ -1,9 +1,28 @@
-function to_uint8(arr::AbstractMatrix{T}) where {T<:AbstractFloat}
-    img = Int.(round.(arr, RoundNearestTiesAway))
-    img = clamp.(img, 0, 255)
+function to_uint8(arr::AbstractMatrix{T}) where {T<:Union{AbstractFloat,Int,Signed}}
+    img = to_uint8.(arr)
     return img
 end
 
+"""
+    rgb2gray(rgbchannels::Array{Float64, 3})
+
+Convert an array of RGB channel data to grayscale in the range [0, 255].
+"""
+function rgb2gray(rgbchannels::Array{Float64,3})
+    r, g, b = [to_uint8(rgbchannels[:, :, i]) for i in 1:3]
+    # Reusing the r array to store the equalized gray image
+    r .= to_uint8(0.2989 * r .+ 0.5870 * g .+ 0.1140 * b)
+    return r
+end
+
+function to_uint8(num::T) where {T<:Union{AbstractFloat,Int,Signed}}
+    num = Int(round(num, RoundNearestTiesAway))
+    return clamp(num, 0, 255)
+end
+
+function mat2uint8(arr::AbstractMatrix{T}) where {T<:Union{AbstractFloat,Int,Signed}}
+    return to_uint8.(arr)
+end
 
 function anisotropic_diffusion_3D(I)
     rgbchannels = get_rgb_channels(I)
@@ -13,10 +32,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 +54,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 +71,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 +80,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)
+    image_min, image_max = extrema(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 +139,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,
@@ -136,6 +154,7 @@ function _process_image_tiles(
     # Apply diffuse (anisotropic diffusion) to each channel of true color image
     true_color_diffused = IceFloeTracker.diffusion(float64.(true_color_image), 0.1, 75, 3)
 
+    # Get uint8 channels
     rgbchannels = get_rgb_channels(true_color_diffused)
 
     # For each tile, compute the entropy in the false color tile, and the fraction of white and black pixels
@@ -153,10 +172,9 @@ function _process_image_tiles(
         end
     end
 
-    return rgbchannels
+    return (equalized_gray=rgb2gray(rgbchannels), gammagreen=rgbchannels[:, :, 2])
 end
 
-
 """
     conditional_histeq(
     true_color_image,
@@ -192,7 +210,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 +221,6 @@ function conditional_histeq(
     )
 
     return rgbchannels_equalized
-
 end
 
 """
@@ -227,7 +244,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 +258,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 +295,31 @@ Convert an RGB image to grayscale in the range [0, 255].
 function rgb2gray(img::Matrix{RGB{Float64}})
     return round.(Int, Gray.(img) * 255)
 end
+
+# TODO: make imadjust more general to work with Gray{Float64} and RGB{Float64} types
+"""
+    imadjust(img; low, high)
+
+Adjust the contrast of an image using linear stretching. The image is normalized to [0, 1] and then stretched to the range [low, high].
+
+# Arguments
+- `img`: The input image.
+- `low`: The lower bound of the stretched image. Default is 0.01.
+- `high`: The upper bound of the stretched image. Default is 0.99.
+
+# Returns
+
+The contrast-adjusted image in the range [0, 255].
+"""
+function imadjust(
+    img::AbstractArray{<:Integer}; low::T=0.01, high::T=0.99
+)::Matrix{Int} where {T<:AbstractFloat}
+    img = img ./ 255
+    imgflat = vec(img)
+    plow = StatsBase.percentile(imgflat, low * 100)
+    phigh = StatsBase.percentile(imgflat, high * 100)
+
+    f = LinearStretching((plow, phigh) => (0.0, 1.0))
+
+    return to_uint8(adjust_histogram(img, f) * 255)
+end

src/landmask.jl

@@ -17,7 +17,7 @@ function create_landmask(
     fill_value_upper::Int=2000,
 ) where {T<:AbstractMatrix}
     landmask_binary = binarize_landmask(landmask_image)
-    dilated = IceFloeTracker.MorphSE.dilate(landmask_binary, centered(struct_elem))
+    dilated = sk_morphology.dilation(landmask_binary, centered(struct_elem))
     return (
         dilated=ImageMorphology.imfill(.!dilated, (fill_value_lower, fill_value_upper)),
         non_dilated=.!landmask_binary,
@@ -52,7 +52,7 @@ Zero out pixels in land and soft ice regions on truecolor image, return RGB imag
 
 # Arguments
 - `input_image`: truecolor RGB image
-- `landmask_binary`: binary landmask with 1=land, 0=water/ice 
+- `landmask_binary`: binary landmask with 1=land, 0=water/ice
 
 """
 function apply_landmask(input_image::AbstractMatrix, landmask_binary::BitMatrix)
@@ -69,7 +69,7 @@ end
 """
     remove_landmask(landmask, ice_mask)
 
-Find the pixel indexes that are floating ice rather than soft or land ice. Returns an array of pixel indexes. 
+Find the pixel indexes that are floating ice rather than soft or land ice. Returns an array of pixel indexes.
 
 # Arguments
 - `landmask`: bitmatrix landmask for region of interest