Merge pull request #502 from WilhelmusLab/501-get_ice_labels-with-rel…

…axation

501 get ice labels with relaxation

src/find_ice_labels.jl

@@ -84,3 +84,171 @@ 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

@@ -289,3 +289,46 @@ 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)
 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
+
+"""
+    imhist(img, imgtype::AbstractString="uint8")
+
+Compute the histogram of an image where each possible value is represented in the histogram. The function returns a tuple with the bins and counts of each bin.
+
+# Example
+```jldoctest; setup = :(using IceFloeTracker)
+julia> 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
+]
+
+julia> bins, heights = imhist(img);
+
+julia> [bins[heights .> 0] heights[heights .>0]] # display only non-zero bins and heights
+3×2 Matrix{Int64}:
+ 3   6
+ 4  14
+ 5   5
+
+"""
+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)
+
+    return _imhist(img, rng)
+end

src/utils.jl

@@ -120,35 +120,6 @@ function imextendedmin(img::AbstractArray, h::Int=2, conn::Int=2)::BitMatrix
     return mask_minima .> 0
 end
 
-function impose_minima(I::AbstractArray{T}, BW::AbstractArray{Bool}) where {T<:Integer}
-    marker = 255 .* BW
-    mask = imcomplement(min.(I .+ 1, 255 .- marker))
-    reconstructed = IceFloeTracker.MorphSE.mreconstruct(
-        IceFloeTracker.MorphSE.dilate, marker, mask
-    )
-    return IceFloeTracker.imcomplement(Int.(reconstructed))
-end
-
-function impose_minima(
-    I::AbstractArray{T}, BW::AbstractMatrix{Bool}
-) where {T<:AbstractFloat}
-    # compute shift
-    a, b = extrema(I)
-    rng = b - a
-    h = rng == 0 ? 0.1 : rng / 1000
-
-    marker = -Inf * BW .+ (Inf * .!BW)
-    mask = min.(I .+ h, marker)
-
-    return 1 .- IceFloeTracker.MorphSE.mreconstruct(
-        IceFloeTracker.MorphSE.dilate, 1 .- marker, 1 .- mask
-    )
-end
-
-function imregionalmin(A, conn=2)
-    return ImageMorphology.local_minima(A; connectivity=conn) .> 0
-end
-
 """
     impose_minima(I::AbstractArray{T}, BW::AbstractArray{Bool}) where {T<:Integer}
 
@@ -177,8 +148,19 @@ function impose_minima(
     return 1 .- sk_morphology.reconstruction(1 .- marker, 1 .- mask)
 end
 
-function imregionalmin(A, conn=2)
-    return ImageMorphology.local_minima(A; connectivity=conn) .> 0
+"""
+    imregionalmin(img, conn=2)
+
+Compute the regional minima of the image `img` using the connectivity `conn`.
+
+Returns a bitmatrix of the same size as `img` with the regional minima.
+
+# Arguments
+- `img`: Image object
+- `conn`: Neighborhood connectivity; in 2D, 1 = 4-neighborhood and 2 = 8-neighborhood
+"""
+function imregionalmin(img, conn=2)
+    return ImageMorphology.local_minima(img; connectivity=conn) .> 0
 end
 
 """

test/test-conditional-adaptive-histeq.jl

@@ -99,10 +99,12 @@ function test_histeq()
                 255 255 255
             ],
         ]
+
         _exp = [
             128 128 128
             255 255 255
         ]
+
         expected = [
             [
                 204 204 204 204 204
@@ -115,7 +117,7 @@ function test_histeq()
             _exp,
         ]
 
-        @test all(IceFloeTracker.histeq(imgs[i]) == expected[i] for i in 1:3)
+        @test all(histeq(imgs[i]) == expected[i] for i in 1:3)
     end
 end
 

test/test-tiled-ice-labels.jl

@@ -0,0 +1,30 @@
+using IceFloeTracker: get_image_peaks, get_tiles, get_ice_labels_mask, get_nlabel_relaxation
+
+@testset "tiled ice labels" begin
+    @testset "get_image_peaks" begin
+        Random.seed!(123)
+        img = rand(0:255, 10, 10)
+        l, h = get_image_peaks(img)
+        @test sum(l[1:5]) == 324
+        @test sum(h[1:5]) == 11
+    end
+
+    ref_img = load(falsecolor_test_image_file)
+    tiles = get_tiles(ref_img; rblocks=8, cblocks=6)
+    tile = tiles[1]
+    factor = 255
+    thresholds = [10, 118, 120]
+    morph_residue = readdlm("test_inputs/morph_residue_tile.csv", ',', Int)
+
+    @testset "get_ice_labels" begin
+        # regular use case applies landmask
+        @test sum(get_ice_labels_mask(ref_img[tile...], thresholds, 255)) == 6515
+    end
+
+    @testset "get_nlabel_relaxation" begin
+        # regular use case applies landmask
+        @test get_nlabel_relaxation(
+            ref_img[tile...], morph_residue[tile...], factor, 75, 10, 230
+        ) == 1
+    end
+end