Merge branch 'main' into 501-get_ice_labels-with-relaxation

src/IceFloeTracker.jl

@@ -71,10 +71,19 @@ include("branch.jl")
 include("special_strels.jl")
 include("tilingutils.jl")
 include("histogram_equalization.jl")
+include("watershed.jl")
 include("brighten.jl")
 include("morph_fill.jl")
 include("imcomplement.jl")
 include("imadjust.jl")
+include("ice_masks.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"))

src/ice_masks.jl

@@ -0,0 +1,86 @@
+"""
+    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.
+
+# 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.
+  - `label`: Most frequent label in the ice mask.
+"""
+function get_ice_masks(
+    falsecolor_image::Matrix{RGB{N0f8}},
+    morph_residue::Matrix{<:Integer},
+    landmask::BitMatrix,
+    tiles::S,
+    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,S<:AbstractMatrix{Tuple{UnitRange{Int64},UnitRange{Int64}}}}
+
+    # Make canvases
+    sz = size(falsecolor_image)
+    ice_mask = BitMatrix(zeros(Bool, sz))
+    binarized_tiling = zeros(Int, sz)
+
+    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))
+
+        # TODO: handle case where get_nlabel returns missing
+        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 .> 0)
+end

src/segmentation_a_direct.jl

@@ -142,4 +142,4 @@ function imgradientmag(img)
     Gx = fetch(Gx_future)
     Gy = fetch(Gy_future)
     return hypot.(Gx, Gy)
-end
+end

src/utils.jl

@@ -69,7 +69,7 @@ end
 Load an image from `dir` with filename `fname` into a matrix of `Float64` values. Returns the loaded image.
 """
 function loadimg(; dir::String, fname::String)
-    return joinpath(dir, fname) |> load |> x-> float64.(x)
+    return (x -> float64.(x))(load(joinpath(dir, fname)))
 end
 
 """
@@ -114,10 +114,39 @@ Mimics MATLAB's imextendedmin function that computes the extended-minima transfo
 - `h`: suppress minima below this depth threshold
 - `conn`: neighborhood connectivity; in 2D 1 = 4-neighborhood and 2 = 8-neighborhood
 """
-function imextendedmin(img::AbstractArray; h::Int=2, conn::Int=2)::BitMatrix
+function imextendedmin(img::AbstractArray, h::Int=2, conn::Int=2)::BitMatrix
     mask = ImageSegmentation.hmin_transform(img, h)
     mask_minima = Images.local_minima(mask; connectivity=conn)
-    return Bool.(mask_minima)
+    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
 
 """

src/watershed.jl

@@ -0,0 +1,40 @@
+function watershed1(bw::T) where {T<:Union{BitMatrix,AbstractMatrix{Bool}}}
+    seg = -IceFloeTracker.bwdist(.!bw)
+    mask2 = imextendedmin(seg)
+    seg = impose_minima(seg, mask2)
+    cc = label_components(imregionalmin(seg), trues(3, 3))
+    w = ImageSegmentation.watershed(seg, cc)
+    lmap = labels_map(w)
+    return Images.isboundary(lmap)
+end
+
+function _reconst_watershed(morph_residue::Matrix{<:Integer}, se::Matrix{Bool}=se_disk20())
+    mr_reconst = to_uint8(IceFloeTracker.reconstruct(morph_residue, se, "erosion", false))
+    mr_reconst .= to_uint8(IceFloeTracker.reconstruct(mr_reconst, se, "dilation", true))
+    mr_reconst .= imcomplement(mr_reconst)
+    return mr_reconst
+end
+
+function watershed2(morph_residue, segment_mask, ice_mask)
+    # Task 1: Reconstruct morph_residue
+    task1 = Threads.@spawn begin
+        mr_reconst = _reconst_watershed(morph_residue)
+        mr_reconst = ImageMorphology.local_maxima(mr_reconst; connectivity=2) .> 0
+    end
+
+    # Task 2: Calculate gradient magnitude
+    task2 = Threads.@spawn begin
+        gmag = imgradientmag(histeq(morph_residue))
+    end
+
+    # Wait for both tasks to complete
+    mr_reconst = fetch(task1)
+    gmag = fetch(task2)
+
+    minimamarkers = mr_reconst .| segment_mask .| ice_mask
+    gmag .= impose_minima(gmag, minimamarkers)
+    cc = label_components(imregionalmin(gmag), trues(3, 3))
+    w = ImageSegmentation.watershed(morph_residue, cc)
+    lmap = labels_map(w)
+    return (fgm=mr_reconst, L0mask=isboundary(lmap))
+end

test/test-watershed-workflows.jl

@@ -0,0 +1,24 @@
+@testset "watershed workflows" begin
+    function build_test_image()
+        center1, center2 = -40, 40
+        radius = sqrt(8 * center1^2) * 0.7
+        lims = (floor(center1 - 1.2 * radius), ceil(center2 + 1.2 * radius))
+        x = collect(lims[1]:lims[2])
+
+        function in_circle(x, y, center, radius)
+            return sqrt((x - center)^2 + (y - center)^2) <= radius
+        end
+
+        return [
+            in_circle(xi, yi, center1, radius) || in_circle(xi, yi, center2, radius) for
+            yi in x, xi in x
+        ]
+    end
+    @testset "watershed" begin
+        println("------------------------------------------------")
+        println("------------ Create Watershed Test --------------")
+        @test sum(IceFloeTracker.watershed1(build_test_image())) == 1088
+    end
+
+    # TODO: add test for watershed2
+end