Merge pull request #489 from WilhelmusLab/488-impose-minima

Impose minima and related functions

src/IceFloeTracker.jl

@@ -75,20 +75,23 @@ include("imadjust.jl")
 include("morph_fill.jl")
 include("imcomplement.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()
 
+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.6.0")
     pyimport_conda("rasterio", "rasterio=1.3.7")
     pyimport_conda("jinja2", "jinja2=3.1.2")

src/utils.jl

@@ -120,6 +120,38 @@ function imextendedmin(img::AbstractArray; h::Int=2, conn::Int=2)::BitMatrix
     return Bool.(mask_minima)
 end
 
+"""
+    impose_minima(I::AbstractArray{T}, BW::AbstractArray{Bool}) where {T<:Integer}
+
+Use morphological reconstruction to enforce minima on the input image `I` at the positions where the binary mask `BW` is non-zero.
+
+It supports both integer and grayscale images using different implementations for each.
+"""
+function impose_minima(I::AbstractArray{T}, BW::AbstractArray{Bool}) where {T<:Integer}
+    marker = 255 .* BW
+    mask = imcomplement(min.(I .+ 1, 255 .- marker))
+    reconstructed = sk_morphology.reconstruction(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 .- sk_morphology.reconstruction(1 .- marker, 1 .- mask)
+end
+
+function imregionalmin(A, conn=2)
+    return ImageMorphology.local_minima(A; connectivity=conn) .> 0
+end
+
 """
     bwdist(bwimg)
 
@@ -142,11 +174,11 @@ function padnhood(img, I, nhood)
     tofill = SizedMatrix{3,3}(zeros(Int, 3, 3))
     @views if I == CartesianIndex(1, 1) # top left corner`
         tofill[2:3, 2:3] = img[nhood]
-    elseif I == CartesianIndex(maxr, 1) # bottom left corner 
+    elseif I == CartesianIndex(maxr, 1) # bottom left corner
         tofill[1:2, 2:3] = img[nhood]
-    elseif I == CartesianIndex(1, maxc) # top right corner 
+    elseif I == CartesianIndex(1, maxc) # top right corner
         tofill[2:3, 1:2] = img[nhood]
-    elseif I == CartesianIndex(maxr, maxc) # bottom right corner 
+    elseif I == CartesianIndex(maxr, maxc) # bottom right corner
         tofill[1:2, 1:2] = img[nhood]
     elseif I[1] == 1 # top edge (first row)
         tofill[2:3, 1:3] = img[nhood]
@@ -169,10 +201,10 @@ Get decimal representation of a bit vector `v` with the leading bit at its leftm
 
 Example
 ```
-julia> _bin9todec([0 0 0 0 0 0 0 0 0])    
+julia> _bin9todec([0 0 0 0 0 0 0 0 0])
 0
 
-julia> _bin9todec([1 1 1 1 1 1 1 1 1])    
+julia> _bin9todec([1 1 1 1 1 1 1 1 1])
 511
 ```
 """

test/runtests.jl

@@ -9,6 +9,7 @@ using Random
 using Test
 using TestImages
 using TiledIteration
+using ZipFile
 include("test_error_rate.jl")
 include("config.jl")
 

test/test-img-processing.jl

@@ -1,6 +1,10 @@
-using ZipFile
-using DelimitedFiles: readdlm
-using IceFloeTracker: imgradientmag, to_uint8, imbinarize, adjustgamma, get_holes, se_disk4
+using IceFloeTracker:
+    imgradientmag,
+    to_uint8,
+    imbinarize,
+    adjustgamma,
+    get_holes,
+    impose_minima
 
 @testset "misc. image processing" begin
     r = ZipFile.Reader("test_inputs/coins.zip")
@@ -24,4 +28,11 @@ using IceFloeTracker: imgradientmag, to_uint8, imbinarize, adjustgamma, get_hole
         bw = coins .> 100
         @test sum(get_holes(bw)) == 2536
     end
+
+    @testset "impose_minima" begin
+        img = readdlm("test_inputs/imposemin.csv", ',', Int)
+        marker = falses(size(img))
+        marker[65:70, 65:70] .= true
+        @test sum(impose_minima(img, marker)) == 7675653
+    end
 end