feat: preprocess with tiling

src/IceFloeTracker.jl

@@ -182,6 +182,8 @@ module MorphSE
     include("morphSE/fill_holes.jl")
 end
 
+include("preprocess_tiling.jl")
+
 module Register
     include("Register/CenterIndexedArrays.jl-0.2.0/CenterIndexedArrays.jl")
     include("Register/RegisterCore.jl-0.2.4/src/RegisterCore.jl")

src/find_ice_labels.jl

@@ -84,171 +84,3 @@ function find_ice_labels(
     # @info "Done with ice labels" # to uncomment when logger is added
     return ice_labels
 end
-
-function get_image_peaks(arr, imgtype="uint8")
-    _, heights = imhist(arr, imgtype)
-
-    locs, heights, _ = Peaks.findmaxima(heights)
-
-    # TODO: make this conditional on input args
-    order = sortperm(heights; rev=true)
-    locs, heights = locs[order], heights[order]
-
-    return (locs=locs, heights=heights)
-end
-
-
-function get_ice_labels_mask(ref_img::Matrix{RGB{N0f8}}, thresholds, factor=1)
-    cv = channelview(ref_img)
-    cv = [float64.(cv[i, :, :]) .* factor for i in 1:3]
-    mask_ice_band_7 = cv[1] .< thresholds[1]
-    mask_ice_band_2 = cv[2] .> thresholds[2]
-    mask_ice_band_1 = cv[3] .> thresholds[3]
-    mask = mask_ice_band_7 .* mask_ice_band_2 .* mask_ice_band_1
-    @debug "Found $(sum(mask)) ice pixels"
-    return mask
-end
-
-function get_nlabel(
-    ref_img,
-    morph_residue_labels,
-    factor;
-    band_7_threshold::T=5,
-    band_2_threshold::T=230,
-    band_1_threshold::T=240,
-    band_7_threshold_relaxed::T=10,
-    band_1_threshold_relaxed::T=190,
-    possible_ice_threshold::T=75,
-) where {T<:Integer}
-    _getnlabel(morphr, mask) = StatsBase.mode(morphr[mask])
-
-    # Initial attempt to get ice labels
-    thresholds = (band_7_threshold, band_2_threshold, band_1_threshold)
-    ice_labels_mask = get_ice_labels_mask(ref_img, thresholds, 255)
-    sum(ice_labels_mask) > 1 && return _getnlabel(morph_residue_labels, ice_labels_mask)
-
-    # First relaxation
-    thresholds = (band_7_threshold_relaxed, band_2_threshold, band_1_threshold_relaxed)
-    ice_labels_mask = get_ice_labels_mask(ref_img, thresholds, 255)
-    sum(ice_labels_mask) > 0 && return _getnlabel(morph_residue_labels, ice_labels_mask)
-
-    # Second/Third relaxation
-    return get_nlabel_relaxation(
-        ref_img,
-        morph_residue_labels,
-        factor,
-        possible_ice_threshold,
-        band_7_threshold_relaxed,
-        band_2_threshold,
-    )
-end
-
-function get_nlabel_relaxation(
-    ref_img,
-    morph_residue_labels,
-    factor,
-    possible_ice_threshold,
-    band_7_threshold_relaxed,
-    band_2_threshold,
-)
-    # filter b/c channels (landmasked channels 2 and 3) and compute peaks
-    b, c = [float64.(channelview(ref_img)[i, :, :]) .* factor for i in 2:3]
-    b[b .< possible_ice_threshold] .= 0
-    c[c .< possible_ice_threshold] .= 0
-    pksb, pksc = get_image_peaks.([b, c])
-
-    # return early if no peaks are found
-    !all(length.([pksb.locs, pksc.locs]) .> 2) && return 1
-
-    relaxed_thresholds = [band_7_threshold_relaxed, pksb.locs[2], pksc.locs[2]]
-    ice_labels = get_ice_labels_mask(ref_img, relaxed_thresholds, factor)
-
-    sum(ice_labels) > 0 && return StatsBase.mode(morph_residue_labels[ice_labels])
-
-    # Final relaxation
-    mask_b = b .> band_2_threshold
-    sum(mask_b) > 0 && return StatsBase.mode(morph_residue_labels[mask_b])
-
-    # No mode found
-    return missing
-end
-
-"""
-    get_ice_masks(
-    falsecolor_image,
-    morph_residue,
-    landmask,
-    tiles,
-    binarize;
-    band_7_threshold,
-    band_2_threshold,
-    band_1_threshold,
-    band_7_threshold_relaxed,
-    band_1_threshold_relaxed,
-    possible_ice_threshold,
-    factor,
-)
-Get the ice masks from the falsecolor image and morphological residue given a particualr tiling configuration.
-# Arguments
-- `falsecolor_image`: The falsecolor image.
-- `morph_residue`: The morphological residue image.
-- `landmask`: The landmask.
-- `tiles`: The tiles.
-- `binarize::Bool=true`: Whether to binarize the tiling.
-- `band_7_threshold=5`: The threshold for band 7.
-- `band_2_threshold=230`: The threshold for band 2.
-- `band_1_threshold=240`: The threshold for band 1.
-- `band_7_threshold_relaxed=10`: The relaxed threshold for band 7.
-- `band_1_threshold_relaxed=190`: The relaxed threshold for band 1.
-- `possible_ice_threshold=75`: The threshold for possible ice.
-- `factor=255`: normalization factor to convert images to uint8.
-# Returns
-- A named tuple `(icemask, bin)` where:
-  - `icemask`: The ice mask.
-  - `bin`: The binarized tiling.
-"""
-function get_ice_masks(
-    falsecolor_image,
-    morph_residue,
-    landmask::BitMatrix,
-    tiles,
-    binarize::Bool=true;
-    band_7_threshold::T=5,
-    band_2_threshold::T=230,
-    band_1_threshold::T=240,
-    band_7_threshold_relaxed::T=10,
-    band_1_threshold_relaxed::T=190,
-    possible_ice_threshold::T=75,
-    factor::T=255,
-) where {T<:Integer}
-
-    # Make canvases
-    ice_mask = BitMatrix(zeros(Bool, size(falsecolor_image)))
-    binarized_tiling = zeros(Int, size(falsecolor_image))
-
-    fc_landmasked = apply_landmask(falsecolor_image, landmask)
-
-    for tile in tiles
-        #  Conditionally update binarized_tiling as it's not used in some workflows
-        if binarize
-            binarized_tiling[tile...] .= imbinarize(morph_residue[tile...])
-        end
-
-        morph_residue_seglabels = kmeans_segmentation(Gray.(morph_residue[tile...] / 255))
-
-        floes_label = get_nlabel(
-            fc_landmasked[tile...],
-            morph_residue_seglabels,
-            factor;
-            band_7_threshold=band_7_threshold,
-            band_2_threshold=band_2_threshold,
-            band_1_threshold=band_1_threshold,
-            band_7_threshold_relaxed=band_7_threshold_relaxed,
-            band_1_threshold_relaxed=band_1_threshold_relaxed,
-            possible_ice_threshold=possible_ice_threshold,
-        )
-
-        ice_mask[tile...] .= (morph_residue_seglabels .== floes_label)
-    end
-    return (icemask=ice_mask, bin=binarized_tiling)
-end

src/histogram_equalization.jl

@@ -9,6 +9,11 @@ function to_uint8(arr::AbstractMatrix{T}) where {T<:Integer}
     return img
 end
 
+function to_uint8(num::T) where {T<:Union{AbstractFloat,Int,Signed}}
+    num = Int(round(num, RoundNearestTiesAway))
+    return clamp(num, 0, 255)
+end
+
 function anisotropic_diffusion_3D(I)
     rgbchannels = get_rgb_channels(I)
 
@@ -127,6 +132,18 @@ function get_rgb_channels(img)
     return cat(redc, greenc, bluec; dims=3)
 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 _process_image_tiles(
     true_color_image,
     clouds_red,
@@ -287,7 +304,7 @@ end
 Histogram equalization of `img` using `nbins` bins.
 """
 function histeq(img::S; nbins=64)::S where {S<:AbstractArray{<:Integer}}
-    return to_uint8(sk_exposure.equalize_hist(img, nbins=nbins) * 255)
+    return to_uint8(sk_exposure.equalize_hist(img; nbins=nbins) * 255)
 end
 
 function _imhist(img, rng)

src/ice_masks.jl

@@ -1,18 +1,107 @@
-"""
-    get_ice_masks(
-    falsecolor_image,
-    morph_residue,
-    landmask,
-    tiles,
-    binarize;
-    band_7_threshold,
-    band_2_threshold,
-    band_1_threshold,
-    band_7_threshold_relaxed,
-    band_1_threshold_relaxed,
-    possible_ice_threshold,
+function get_image_peaks(arr, imgtype="uint8")
+    _, heights = imhist(arr, imgtype)
+
+    locs, heights, _ = Peaks.findmaxima(heights)
+
+    # TODO: make this conditional on input args
+    order = sortperm(heights; rev=true)
+    locs, heights = locs[order], heights[order]
+
+    return (locs=locs, heights=heights)
+end
+
+function get_ice_labels_mask(ref_img::Matrix{RGB{N0f8}}, thresholds, factor=255)
+    cv = channelview(ref_img)
+    cv = [float64.(cv[i, :, :]) .* factor for i in 1:3]
+    mask_ice_band_7 = cv[1] .< thresholds[1]
+    mask_ice_band_2 = cv[2] .> thresholds[2]
+    mask_ice_band_1 = cv[3] .> thresholds[3]
+    mask = mask_ice_band_7 .* mask_ice_band_2 .* mask_ice_band_1
+    @debug "Found $(sum(mask)) ice pixels"
+    return mask
+end
+
+function get_nlabel(
+    ref_img,
+    morph_residue_labels,
+    factor;
+    band_7_threshold::T=5,
+    band_2_threshold::T=230,
+    band_1_threshold::T=240,
+    band_7_threshold_relaxed::T=10,
+    band_1_threshold_relaxed::T=190,
+    possible_ice_threshold::T=75,
+) where {T<:Integer}
+    _getnlabel(morphr, mask) = StatsBase.mode(morphr[mask])
+
+    # Initial attempt to get ice labels
+    thresholds = (band_7_threshold, band_2_threshold, band_1_threshold)
+    ice_labels_mask = get_ice_labels_mask(ref_img, thresholds)
+    sum(ice_labels_mask) > 0 && return _getnlabel(morph_residue_labels, ice_labels_mask)
+    @debug "Trying first relaxation."
+
+    # First relaxation
+    thresholds = (band_7_threshold_relaxed, band_2_threshold, band_1_threshold_relaxed)
+    ice_labels_mask = get_ice_labels_mask(ref_img, thresholds)
+    sum(ice_labels_mask) > 0 && return _getnlabel(morph_residue_labels, ice_labels_mask)
+
+    @debug "Trying second/third relaxation."
+    # Second/Third relaxation
+    return get_nlabel_relaxation(
+        ref_img,
+        morph_residue_labels,
+        factor,
+        possible_ice_threshold,
+        band_7_threshold_relaxed,
+        band_2_threshold,
+    )
+end
+
+function get_nlabel_relaxation(
+    ref_img,
+    morph_residue_labels,
     factor,
+    possible_ice_threshold,
+    band_7_threshold_relaxed,
+    band_2_threshold,
 )
+    # filter b/c channels (landmasked channels 2 and 3) and compute peaks
+    b, c = [float64.(channelview(ref_img)[i, :, :]) .* factor for i in 2:3]
+    b[b .< possible_ice_threshold] .= 0
+    c[c .< possible_ice_threshold] .= 0
+    pksb, pksc = get_image_peaks.([b, c])
+
+    # return early if no peaks are found
+    !all(length.([pksb.locs, pksc.locs]) .> 2) && return 1
+
+    relaxed_thresholds = [band_7_threshold_relaxed, pksb.locs[2], pksc.locs[2]]
+    ice_labels = get_ice_labels_mask(ref_img, relaxed_thresholds, factor)
+
+    sum(ice_labels) > 0 && return StatsBase.mode(morph_residue_labels[ice_labels])
+
+    # Final relaxation
+    mask_b = b .> band_2_threshold
+    sum(mask_b) > 0 && return StatsBase.mode(morph_residue_labels[mask_b])
+
+    # No mode found
+    return 1
+end
+
+"""
+    get_ice_masks(
+        falsecolor_image,
+        morph_residue,
+        landmask,
+        tiles,
+        binarize;
+        band_7_threshold,
+        band_2_threshold,
+        band_1_threshold,
+        band_7_threshold_relaxed,
+        band_1_threshold_relaxed,
+        possible_ice_threshold,
+        factor
+    )
 
 Get the ice masks from the falsecolor image and morphological residue given a particular tiling configuration.
 
@@ -29,12 +118,10 @@ Get the ice masks from the falsecolor image and morphological residue given a pa
 - `band_1_threshold_relaxed=190`: The relaxed threshold for band 1.
 - `possible_ice_threshold=75`: The threshold for possible ice.
 - `factor=255`: normalization factor to convert images to uint8.
-
 # Returns
 - A named tuple `(icemask, bin)` where:
   - `icemask`: The ice mask.
   - `bin`: The binarized tiling.
-  - `label`: Most frequent label in the ice mask.
 """
 function get_ice_masks(
     falsecolor_image::Matrix{RGB{N0f8}},
@@ -48,6 +135,7 @@ function get_ice_masks(
     band_7_threshold_relaxed::T=10,
     band_1_threshold_relaxed::T=190,
     possible_ice_threshold::T=75,
+    k::T=3,
     factor::T=255,
 ) where {T<:Integer,S<:AbstractMatrix{Tuple{UnitRange{Int64},UnitRange{Int64}}}}
 
@@ -56,15 +144,13 @@ function get_ice_masks(
     ice_mask = BitMatrix(zeros(Bool, sz))
     binarized_tiling = zeros(Int, sz)
 
-    fc_landmasked = apply_landmask(falsecolor_image, landmask)
+    fc_landmasked = apply_landmask(falsecolor_image, .!landmask)
 
-    Threads.@threads for tile in tiles
-        #  Conditionally update binarized_tiling as its not used in some workflows
-        if binarize
-            binarized_tiling[tile...] .= imbinarize(morph_residue[tile...])
-        end
-
-        morph_residue_seglabels = kmeans_segmentation(Gray.(morph_residue[tile...] / 255))
+    # Threads.@threads
+    for tile in tiles
+        @debug "Processing tile: $tile"
+        mrt = morph_residue[tile...]
+        morph_residue_seglabels = kmeans_segmentation(Gray.(mrt / 255); k=k)
 
         # TODO: handle case where get_nlabel returns missing
         floes_label = get_nlabel(
@@ -80,7 +166,10 @@ function get_ice_masks(
         )
 
         ice_mask[tile...] .= (morph_residue_seglabels .== floes_label)
+
+        #  Conditionally update binarized_tiling as its not used in some workflows
+        binarize && (binarized_tiling[tile...] .= imbinarize(mrt))
     end
-       
+
     return (icemask=ice_mask, bin=binarized_tiling .> 0)
 end