Merge branch 'main' into 490-add-watershed-workflows

docs/src/index.md

@@ -1,6 +1,19 @@
-# IceFloeTracker
+# IceFloeTracker.jl
 
-Documentation for IceFloeTracker.jl package.
+## Overview
+IceFloeTracker.jl is a collection of routines and tools for processing remote sensing imagery, identifying sea ice floes, and tracking the displacement and rotation of ice floes across multiple images. It can be used either standalone to create custom processing pathways or with the [Ice Floe Tracker Pipeline](https://github.com/WilhelmusLab/ice-floe-tracker-pipeline).
+
+## Algorithm components
+The Ice Floe Tracker (IFT) package includes the core functions for the three main steps of the algorithm. These functions can be used independently and can be customized for specific use cases. 
+
+### Preprocessing
+IFT operates on optical satellite imagery. The main functions are designed with "true color" and "false color" imagery in mind, and have thus far primarily been tested on imagery from the Moderate Resolution Imaging Spectroradiometer (MODIS) from the NASA _Aqua_ and _Terra_ satellites. The preprocessing routines mask land and cloud features, and aim to adjust and sharpen the remainder of the images to amplify the contrast along the edges of sea ice floes. (TBD: Link to main preprocessing page)
+
+### Segmentation
+The IFT segmentation functions include functions for semantic segmentation (pixel-by-pixel assignment into predefined categories) and object-based segmentation (groupings of pixels into distinct objects). The semantic segmentation steps use $k$-means to group pixels into water and ice regions. A combination of watershed functions, morphological operations, and further applications of $k$-means are used to identify candidate ice floes. (TBD: Link to main segmentation page)
+
+### Tracking
+Ice floe tracking is carried out by comparing the shapes produced in the segmentation step. Shapes with similar area are rotated until the difference in surface area is minimized, and then the edge shapes are compared using a Ѱ-s curve. If thresholds for correlation and area differences are met, then the floe with the best correlation and smallest area differences are considered matches and the objects are assigned the same label. In the end, trajectories for individual floes are recorded in a dataframe.
 
 ```@contents
 ```
@@ -13,3 +26,18 @@ Order   = [:function, :macro, :type]
 ## Index
 ```@index
 ```
+
+## Developers
+IceFloeTracker.jl is a product of the [Wilhelmus Lab](https://www.wilhelmuslab.me) at Brown University, led by Monica M. Wilhelmus. The original algorithm was developed by Rosalinda Lopez-Acosta during her PhD work at University of California Riverside, advised by Dr. Wilhelmus. The translation of the original Matlab code into the current modular, open source Julia package has been carried out in conjunction with the Center for Computing and Visualization at Brown University. Contributors include Daniel Watkins, Maria Isabel Restrepo, Carlos Paniagua, Tim DiVoll, John Holland, and Bradford Roarr.
+
+## Citing
+
+If you use IceFloeTracker.jl in research, teaching, or elsewhere, please mention the IceFloeTracker package and cite our journal article outlining the algorithm:
+
+Lopez-Acosta et al., (2019). Ice Floe Tracker: An algorithm to automatically retrieve Lagrangian trajectories via feature matching from moderate-resolution visual imagery. _Remote Sensing of Environment_, **234(111406)**, doi:[10.1016/j.rse.2019.111406](https://doi.org/10.1016/j.rse.2019.111406).
+
+## Papers using Ice Floe Tracker
+1. Manucharyan, Lopez-Acosta, and Wilhelmus (2022)\*. Spinning ice floes reveal intensification of mesoscale eddies in the western Arctic Ocean. _Scientific Reports_, **12(7070)**, doi:[10.1038/s41598-022-10712-z](https://doi.org/10.1038/s41598-022-10712-z)
+2. Watkins, Bliss, Hutchings, and Wilhelmus (2023)\*. Evidence of Abrupt Transitions Between Sea Ice Dynamical Regimes in the East Greenland Marginal Ice Zone. _Geophysical Research Letters_, **50(e2023GL103558)**, pp. 1-10, doi:[10.1029/2023GL103558](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023GL103558)
+
+\* Papers using data from the Matlab implementation of Ice Floe Tracker.

src/IceFloeTracker.jl

@@ -72,13 +72,10 @@ include("special_strels.jl")
 include("tilingutils.jl")
 include("histogram_equalization.jl")
 include("watershed.jl")
-include("imadjust.jl")
+include("brighten.jl")
 include("morph_fill.jl")
 include("imcomplement.jl")
-
-const sk_measure = PyNULL()
-const sk_exposure = PyNULL()
-const getlatlon = PyNULL()
+include("imadjust.jl")
 
 function get_version_from_toml(pth=dirname(dirname(pathof(IceFloeTracker))))::VersionNumber
     toml = TOML.parsefile(joinpath(pth, "Project.toml"))
@@ -87,9 +84,16 @@ 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/bridge.jl

@@ -1,15 +1,15 @@
 include("./lut/lutbridge.jl")
 
+const LUTBRIDGE = make_lutbridge()
+
 function _bridge_operator_lut(
     I::CartesianIndex{2},
     img::AbstractArray{Bool},
     nhood::CartesianIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}},
 )
-    lutbridge = make_lutbridge()
-    return _operator_lut(I, img, nhood, lutbridge)
+    return _operator_lut(I, img, nhood, LUTBRIDGE)
 end
 
-# 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)

src/brighten.jl

@@ -0,0 +1,28 @@
+"""
+    get_brighten_mask(equalized_gray_reconstructed_img, gamma_green)
+# Arguments
+- `equalized_gray_reconstructed_img`: The equalized gray reconstructed image (uint8 in Matlab).
+- `gamma_green`: The gamma value for the green channel (also uint8).
+# Returns
+Difference equalized_gray_reconstructed_img - gamma_green clamped between 0 and 255.
+"""
+function get_brighten_mask(equalized_gray_reconstructed_img, gamma_green)
+    return to_uint8(equalized_gray_reconstructed_img - gamma_green)
+end
+
+"""
+    imbrighten(img, brighten_mask, bright_factor)
+Brighten the image using a mask and a brightening factor.
+# Arguments
+- `img`: The input image.
+- `brighten_mask`: A mask indicating the pixels to brighten.
+- `bright_factor`: The factor by which to brighten the pixels.
+# Returns
+- The brightened image.
+"""
+function imbrighten(img, brighten_mask, bright_factor)
+    img = Float64.(img)
+    brighten_mask = brighten_mask .> 0
+    img[brighten_mask] .= img[brighten_mask] * bright_factor
+    return img = to_uint8(img)
+end

src/hbreak.jl

@@ -61,7 +61,8 @@ hbreak(img) = hbreak!(copy(img))
 Inplace version of `hbreak`. See also [`hbreak`](@ref).
 """
 function hbreak!(img::T)::T where {T<:AbstractArray{Bool}}
-    f = x -> get(make_hbreak_dict(), x, false)
+    hbreak_dict = make_hbreak_dict()
+    f = x -> get(hbreak_dict, x, false)
     R = CartesianIndices(img)
     I_first, I_last = first(R), last(R)
     Δ = CartesianIndex(1, 1)

src/histogram_equalization.jl

@@ -4,6 +4,11 @@ function to_uint8(arr::AbstractMatrix{T}) where {T<:AbstractFloat}
     return img
 end
 
+function to_uint8(arr::AbstractMatrix{T}) where {T<:Integer}
+    img = clamp.(arr, 0, 255)
+    return img
+end
+
 function anisotropic_diffusion_3D(I)
     rgbchannels = get_rgb_channels(I)
 
@@ -275,46 +280,12 @@ function rgb2gray(img::Matrix{RGB{Float64}})
     return round.(Int, Gray.(img) * 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
-
-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)
-    # use range(0, stop=1, length=256) for grayscale images
-
-    return _imhist(img, rng)
-end
-
-"""
-    imhistp(img)
-Parsimonious version of `imhist` that uses the maximum value of the image to determine the number of bins.
-"""
-function imhistp(img)
-    nbins = nextpow(2, maximum(img) + 1)
-    return _imhist(img, 0:(nbins-1))
-end
-
-
 """
     histeq(img)
-Histogram equalization of `img` according to [1].
-[1] R. C. Gonzalez and R. E. Woods. Digital Image Processing (3rd Edition). Upper Saddle River, NJ, USA: Prentice-Hall, 2006.
+    histeq(img; nbins=64)
+
+Histogram equalization of `img` using `nbins` bins.
 """
-function histeq(img::S)::S where {S<:AbstractArray{<:Integer}}
-    k, heights = imhistp(img)
-    cdf = last(k) * cumsum(heights) / sum(heights)
-    s = round.(Int, cdf, RoundNearestTiesAway)
-    mapping = Dict(zip(k, s))
-    T(r) = mapping[r]
-    return T.(img)
+function histeq(img::S; nbins=64)::S where {S<:AbstractArray{<:Integer}}
+    return to_uint8(sk_exposure.equalize_hist(img, nbins=nbins) * 255)
 end

src/morph_fill.jl

@@ -1,11 +1,13 @@
 include("./lut/lutfill.jl")
 
+const LUTFILL = make_lutfill()
+
 function _fill_operator_lut(
     I::CartesianIndex{2},
     img::AbstractArray{Bool},
     nhood::CartesianIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}},
 )
-    return _operator_lut(I, img, nhood, make_lutfill())
+    return _operator_lut(I, img, nhood, LUTFILL)
 end
 
 """

src/tracker/tracker-funcs.jl

@@ -103,7 +103,7 @@ end
 
 Sort the floes in `tracked` by area in descending order.
 """
-function sort!(tracked::Tracked)
+function Base.sort!(tracked::Tracked)
     for container in tracked.data
         p = sortperm(container.props1, "area"; rev=true)
         for nm in namesof(container)

src/utils.jl

@@ -149,6 +149,38 @@ 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}
+
+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)
 

test/test-brighten.jl

@@ -0,0 +1,22 @@
+using IceFloeTracker: get_brighten_mask, imbrighten
+
+@testset "brighten tests" begin
+    @testset "get_brighten_mask" begin
+        img = rand(0:255, 5, 5)
+        bumped_img = img .+ 1
+        mask = get_brighten_mask(img, bumped_img)
+        @test all(mask .== 0)
+    end
+
+    @testset "imbrighten tests" begin
+        img = [1 2; 3 4]
+        brighten_mask = [1 0; 1 0]
+
+        test_cases = [(1.25, [1 2; 4 4]), (0.1, [0 2; 0 4]), (0.9, img)]
+
+        for (bright_factor, expected_result) in test_cases
+            result = imbrighten(img, brighten_mask, bright_factor)
+            @test result == expected_result
+        end
+    end
+end

test/test-conditional-adaptive-histeq.jl

@@ -99,22 +99,20 @@ function test_histeq()
                 255 255 255
             ],
         ]
+        _exp = [
+            128 128 128
+            255 255 255
+        ]
         expected = [
             [
-                6 6 6 6 6
-                2 6 7 6 2
-                2 7 7 7 2
-                2 6 7 6 2
-                6 6 6 6 6
-            ],
-            [
-                2 2 2
-                3 3 3
-            ],
-            [
-                128 128 128
-                255 255 255
+                204 204 204 204 204
+                61 204 255 204 61
+                61 255 255 255 61
+                61 204 255 204 61
+                204 204 204 204 204
             ],
+            _exp,
+            _exp,
         ]
 
         @test all(IceFloeTracker.histeq(imgs[i]) == expected[i] for i in 1:3)

test/test-imadjust.jl

@@ -1,5 +1,5 @@
 @testset "imadjust" begin
     Random.seed!(123)
     img = rand(0:255, 100, 100)
-    sum(IceFloeTracker.imadjust(img)) == 1291155
+    @test sum(IceFloeTracker.imadjust(img)) == 1291155
 end