IceFloeTracker.jl

IceDetectionAlgorithm

Functors to detect ice regions in an image.

Each algorithm a with parameters kwargs... can be called like:

• binarize(image, a(; kwargs...))
• or a(; kwargs...)(image).

IceDetectionBrightnessPeaksMODIS721(;
    band_7_threshold::Real,
    possible_ice_threshold::Real
)(image)
binarize(
    modis_721_image, 
    a::IceDetectionBrightnessPeaksMODIS721
)

Returns pixels for a MODIS image where (band7 < threshold AND both (band2, band_1) are are above a peak value above some threshold).

IceDetectionFirstNonZeroAlgorithm(;
    algorithms::Vector{IceDetectionAlgorithm},
)(image)
binarize(image, algorithms::IceDetectionFirstNonZeroAlgorithm)

Runs each algorithm from algorithms on the image, and returns the first which detects any ice.

IceDetectionThresholdMODIS721(;
    band_7_threshold::Real,
    band_2_threshold::Real,
    band_1_threshold::Real,
)(image)
binarize(
    modis_721_image, 
    a::IceDetectionThresholdMODIS721
)

Returns pixels for a MODIS image where (band7 < threshold AND band2 > threshold AND band_1 > threshold).

Results of a segmentation comparison

Results of a segmentation comparison

IceDetectionLopezAcosta2019(;
    band_7_threshold::Float64=Float64(5 / 255),
    band_2_threshold::Float64=Float64(230 / 255),
    band_1_threshold::Float64=Float64(240 / 255),
    band_7_threshold_relaxed::Float64=Float64(10 / 255),
    band_1_threshold_relaxed::Float64=Float64(190 / 255),
    possible_ice_threshold::Float64=Float64(75 / 255),
)

Returns the first non-zero result of two threshold-based and one brightness-peak based ice detections.

Default thresholds are defined in the published Ice Floe Tracker article: Remote Sensing of the Environment 234 (2019) 111406.

addlatlon(pairedfloesdf::DataFrame, refimage::AbstractString)

Add columns latitude, longitude, and pixel coordinates x, y to pairedfloesdf.

Arguments

• pairedfloesdf: dataframe containing floe tracking data.
• refimage: path to reference image.

Find pixels in a segmented image with non-zero labels

convertcentroid!(propdf, latlondata, colstodrop)

Convert the centroid coordinates from row and column to latitude and longitude dropping unwanted columns specified in colstodrop for the output data structure. Addionally, add columns x and y with the pixel coordinates of the centroid.

converttounits!(propdf, latlondata, colstodrop)

Convert the floe properties from pixels to kilometers and square kilometers where appropiate. Also drop the columns specified in colstodrop.

find_ice_labels(falsecolor_image, landmask; band_7_threshold, band_2_threshold, band_1_threshold, band_7_relaxed_threshold, band_1_relaxed_threshold, possible_ice_threshold)

Returns pixel indices of likely ice from false color reflectance image, using the thresholds from the Ice Floe Tracker article: Remote Sensing of the Environment 234 (2019) 111406.

Arguments

• falsecolor_image: corrected reflectance false color image - bands [7,2,1]
• landmask: bitmatrix landmask for region of interest
• band_7_threshold: threshold value used to identify ice in band 7, N0f8(RGB intensity/255)
• band_2_threshold: threshold value used to identify ice in band 2, N0f8(RGB intensity/255)
• band_1_threshold: threshold value used to identify ice in band 2, N0f8(RGB intensity/255)
• band_7_relaxed_threshold: threshold value used to identify ice in band 7 if not found on first pass, N0f8(RGB intensity/255)
• band_1_relaxed_threshold: threshold value used to identify ice in band 1 if not found on first pass, N0f8(RGB intensity/255)

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,
    k
)

Identifies potential sea ice floes using two methods: selection of a relevant k-means cluster and application of adaptive threshold binarization. For the k-means section, a series of thresholds on band 7, 2, and 1 reflectance are applied in order to find the cluster containing bright sea ice pixels.

Arguments

• falsecolor_image: MODIS False Color Bands 7-2-1.
• morph_residue: Grayscale sharpened and equalized image from preprocessing workflow.
• landmask: Binary landmask.
• tiles: Iterable with tile divisions.
• binarize::Bool=true: Whether to binarize the tiling.
• band_7_threshold=5/255: The threshold for band 7.
• band_2_threshold=230/255: The threshold for band 2.
• band_1_threshold=240/255: 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/255: The threshold for possible ice.
• k=4: The number of clusters to use for k-means segmentation.

Returns

• Binary image with likely sea ice floes = 1.

Given the edges and counts from build_histogram, identify local maxima and return the location of the largest local maximum that is bright enough that it is possibly sea ice. Locations are determined by the edges, which by default are the left bin edges. Note also that peaks defaults to the left side of plateaus. Returns Inf if there are no non-zero parts of the histogram with bins larger than the possible ice threshold, or if there are no detected peaks larger than the minimum prominence.

kmeans_binarization(gray_image, false_color_image, tiles; kwargs...)

Produce a binarized image tilewise by identifying pixels of bright ice, performing k-means clustering, and then selecting the k-means cluster containing the largest fraction of bright ice pixels. If no bright ice pixels are detected, then a blank matrix is returned.

Warning: Tilewise processing may result in discontinuities at tile boundaries.

Positional Arguments

• gray_image: Grayscale image to segment using k-means.
• falsecolor_image: MODIS 721 false color image to be supplied to the cluster selection algorithm.

Keyword arguments

• k: Number of k-means clusters. Default 4.
• maxiter: Maximum number of iterations for k-means algorithm. Default 50.
• random_seed: Seed for the random number generator. Default 45.
• cluster_selection_algorithm: Binarization function to find the k-means cluster to set to 1; all other clusters set to 0.

kmeans_binarization(gray_image, false_color_image; kwargs...)

Produce a binarized image by identifying pixels of bright ice, performing k-means clustering, and then selecting the k-means cluster containing the largest fraction of bright ice pixels. If no bright ice pixels are detected, then a blank matrix is returned.

Positional Arguments

• gray_image: Grayscale image to segment using k-means.
• falsecolor_image: MODIS 7-2-1 falsecolor image, to be sent to the specified ice_labels_algorithm. It is recommended that this image be landmasked.

Keyword arguments

• ice_labels_algorithm: Binarization function to find sea ice pixels
• k: Number of k-means clusters
• maxiter: Maximum number of iterations for k-means algorithm
• random_seed: Seed for the random number generator

kmeans_segmentation(img; k=4, maxiter=50, random_seed=45, k_offset=0)

Wrapper for Clustering.kmeans which accepts a grayscale image and returns a SegmentedImage object. Optionally, one can specify the number of clusters k, the maximum number of iterationsmaxiter, and the seed for the random number generator,random_seed`. Returns a SegmentedImage object.

regionprops(label_img, ; properties, connectivity)

A wrapper of the regionprops function from the skimage python library.

See its full documentation at https://scikit-image.org/docs/stable/api/skimage.measure.html#skimage.measure.regionprops.

Arguments

• label_img: Image with the labeled objects of interest
• intensity_img: (Optional) Used for generating extra_properties, integer/float array from which (presumably) label_img was generated
• extra_properties: (Optional) not yet implemented. It will be set to nothing

See also regionprops_table

Examples

julia> Random.seed!(123);

julia> bw_img = rand([0, 1], 5, 10)
5×10 Matrix{Int64}:
 1  0  1  0  0  0  0  0  0  1
 1  0  1  1  1  0  0  0  1  1
 1  1  0  1  1  0  1  0  0  1
 0  1  0  1  0  0  0  0  1  0
 1  0  0  0  0  1  0  1  0  1

 julia> label_img = Images.label_components(bw_img, trues(3,3))
 5×10 Matrix{Int64}:
  1  0  1  0  0  0  0  0  0  4
  1  0  1  1  1  0  0  0  4  4
  1  1  0  1  1  0  3  0  0  4
  0  1  0  1  0  0  0  0  4  0
  1  0  0  0  0  2  0  4  0  4

 julia> regions = regionprops(label_img, bw_img);

 julia> for region in regions
           println(region.area,"	", region.perimeter)
        end
13      11.621320343559642
1       0.0
1       0.0
7       4.621320343559642

regionprops_table(label_img, intensity_img; properties, connectivity, extra_properties)

A wrapper of the regionprops_table function from the skimage python library.

See its full documentation at https://scikit-image.org/docs/stable/api/skimage.measure.html#regionprops-table.

Arguments

• label_img: Image with the labeled objects of interest
• intensity_img: (Optional) Used for generating extra_properties, integer/float array from which (presumably) label_img was generated
• properties: List (Vector or Tuple) of properties to be generated for each connected component in label_img
• extra_properties: (Optional) not yet implemented. It will be set to nothing

Notes

• Zero indexing has been corrected for the bbox and centroid properties
• bbox data (max_col and max_row) are inclusive
• centroid data are rounded to the nearest integer

See also regionprops

Examples

julia> using IceFloeTracker, Random, Images

julia> Random.seed!(123);

julia> bw_img = rand([0, 1], 5, 10)
5×10 Matrix{Int64}:
 1  0  1  0  0  0  0  0  0  1
 1  0  1  1  1  0  0  0  1  1
 1  1  0  1  1  0  1  0  0  1
 0  1  0  1  0  0  0  0  1  0
 1  0  0  0  0  1  0  1  0  1

julia> label_img = label_components(bw_img, trues(3,3))
5×10 Matrix{Int64}:
 1  0  1  0  0  0  0  0  0  4
 1  0  1  1  1  0  0  0  4  4
 1  1  0  1  1  0  3  0  0  4
 0  1  0  1  0  0  0  0  4  0
 1  0  0  0  0  2  0  4  0  4

julia> properties = ["area", "perimeter"]
2-element Vector{String}:
 "area"
 "perimeter"

 julia> regionprops_table(label_img, bw_img, properties = properties)
 4×2 DataFrame
  Row │ area   perimeter 
      │ Int32  Float64   
 ─────┼──────────────────
    1 │    13   11.6213
    2 │     1    0.0
    3 │     1    0.0
    4 │     7    4.62132

function segmentationcomparison( validated::SegmentedImage, measured::SegmentedImage )::@NamedTuple{recall::Real, precision::Real, Fscore::Real}

Compares two SegmentedImages and returns values describing how similar the segmentations are.

This treats the segment labeled 0 as background.

Measures:

• precision: rate at which pixels in validated segments belong to measured segments
• recall: rate at which pixels in measured segments belong to validated segments
• F_score: harmonic mean of precision and recall

stitchclusters(tiles, segmentedimage, minimumoverlap, grayscalethreshold)

Stitches clusters across tile boundaries based on neighbor with largest shared boundary. The algorithm finds all pairs of segment labels at the tile edges. Then, we count the number of times each right-hand label is paired to a left-hand label, and for pairs with at least minimum_overlap pixel overlap, the right-hand label is assigned as a candidate pair to the left-hand label. If the difference in grayscale intensity is less than grayscale_threshold, the objects are merged. The function returns an image index map.

tiledadaptivebinarization(img, tiles; minimumwindowsize=).

Applies the (AdaptiveThreshold)[https://juliaimages.org/ImageBinarization.jl/v0.1/#Adaptive-Threshold-1] binarization algorithm
to each tile in the image. Following the recommendations from ImageBinarization, the default is to use the integer window size
nearest to 1/8th the tile size if the tile is large enough. So that the window is large enough to include moderately large floes,
the default minimum window size is 100 pixels (25 km for MODIS imagery). The minimum brightness parameter masks pixels with low
grayscale intensity to prevent dark regions from getting brightened (i.e., the center of a large patch of open water).
The "threshold_percentage" parameter is passed to the the AdaptiveThreshold function (percentage parameter).

ChainedFilterFunction(filters::Vector{AbstractFloeFilterFunction})

A ChainedFilterFunction is a composite function based on a set of AbstractFloeFilterFunctions. Each is applied in sequence. Thus a filter function based on the DistanceThresholdFilter and area relative error could be made as

filter_function = ChainedFilterFunction(
    filters=[DistanceThresholdFilter(), RelativeErrorThresholdFilter(variable=:area)]
)

DistanceThresholdFilter(time_column, dist_column, threshold_function, threshold_column)
DistanceThresholdFilter(floe, candidates)

The distance threshold filter creates columns for time and distance and applies a threshold function to these columns to determine if the net travel is physically possible. The struct is initialized with names for the time and distance columns, the threshold function (a TimeDistanceFunction) and the name of the column in which to store the results.

julia> dt_test = DistanceThresholdFilter(time_colum=:Δt, dist_column=:Δx, threshold_function=LinearTimeDistanceFunction())

Now, let's assume that floe and candidates are already defined. Then

julia> dt_test(floe, candidates)

will modify candidates in place to include only rows in which the LinearTimeDistanceFunction() evaluates as true. Passing Val{:raw} as the third argument will forgo the subsetting step so that the output of the test can be examined.

Tests the travel distance and time to a linear estimate of maximum travel distance using the formula

max_dx = max_vel * dt + eps

Epsilon should be the uncertainty in position, such that if for example the positional uncertainty is 250 m, then the maximum distance includes a 250 m buffer. The default maximum velocity is 1.5 m/s.

Tests the travel distance and time in log-log space against an empirically fitted quadratic function. The function is constrained by minimum and maximum times. Times less than the minimum are subject to the maximum 1-hour travel distance, while times larger than the maximum fail automatically. See Watkins et al. 2025 for details.

LopezAcostaTimeDistanceFunction(Δx, Δt; dt, dx)

Stepwise time delta function based on Lopez-Acosta et al. 2019. The time thresholds and the input Δt must be time objects so that conversion to seconds is possible. Displacement distances are assumed to be in meters. The final dx value is the maximum displacement.

MinimumWeightMatchingFunction(columns=[:scaled_distance, :relative_error_area, ...])
MinimumWeightedMatchingFunction(candidate_pairs)

Function to identify a best matching between pairs of ice floes in the DataFrame candidate_pairs. The columns variable is instantiated by the first functor call and is used to select a list of columns in candidate_pairs to sum. The result of the sum is the weight assigned to each pairing. Then, a best set of unique pairs is found by carrying out two grouped minimizations: first grouping by the first floe, identified by the head_uuid column, and finding the floe with the smallest weight, then grouping by the second floe, identified by the uuid column, and again finding the floe with the smallest weight. Finally, only pairs that exist in both the forward and backward grouped minizations are identified as likely true matches.

The piecewise linear threshold function is defined using two (area, value) pairs. For areas below the minimum area, it is constant at minimum value; likewise for above the maximum area. The threshold function is linear in between these two points. A return value true indicates that the value is below the threshold.

PsiSCorrelationThresholdFunction(area_variable, threshold_column, threshold_function)
PsiSCorrelationThresholdFunction(floe, candidates, Val(:raw))

Compute the psi-s correlation between a floe and a dataframe of candidate floes. Adds the psi-s correlation, psi-s correlation score (1 - correlation), and the result of the threshold function to the columns of candidates.

RelativeErrorThresholdFilter(variable, area_variable, threshold_column, threshold_function)
RelativeErrorThresholdFilter(floe, candidates)
RelativeErrorThresholdFilter(floe, candidates, Var(:raw))

Compute and test (absolute) relative error for variable. The relative error between scalar variables X and Y is defined as

err = abs(X - Y)/mean(X, Y)

This function takes a scalar <variable> (which must be a named column in the candidates DataFrame) and computes the relative error. Calling the function with the variable name, area_variable, threshold_column name, and athresholdfunctioninitializes the function and saves the parameter values. Once initialized, the function takes a floe (DataFrameRow) and a DataFrame of candidate floes as arguments, and subsets the candidates to only those which evaluate as true using thethresholdfunction. Including the dummy variableVar(:raw)` returns the candidates dataframe with the test results without subsetting it.

ShapeDifferenceThresholdFilter(area_variable, scale_by, threshold_column, threshold_function)
ShapeDifferenceThresholdFilter(floe, candidates)
ShapeDifferenceThresholdFilter(floe, candidates, Val(:raw))

Compute and test the scaled shape difference between input floe and each floe in the dataframe candidates. Assumes that the shape difference test operates on the shape difference scaled by a variable scale_by and the shape difference test depends on the area.

The stepwise linear threshold function is defined using a changepoint area and two levels. If the area is less than the changepoint area, the function returns true if the value is below low_value and false otherwsie; if the area is greater than or equal to the changepoint area, then the value is tested againg high_value.

FilterFunction()

The default filter function for the FloeTracker. The function is an instance of ChainedFilterFunction, applying 7 individual AbstractFloeFilterFunctions in sequence: 1. DistanceThresholdFilter 2-5. RelativeErrorThresholdFilters for area, convexarea, majoraxislength, and minoraxis_length 6. ShapeDifferenceThresholdFilter 7. PsiSCorrelationThresholdFilter Filters 2-7 use PiecewiseLinearThresholdFunctions for thresholds, while Filter 1 uses a LinearTimeDistanceFunction. The default values and settings are derived in Watkins et al. 2026.

add_floemasks!(props::DataFrame, floeimg::FloeLabelsImage)
add_floemasks!.(props::Vector{DataFrame}, floeimgs::Vector{FloeLabelsImage})

Add a column to props called mask containing the cropped floe masks from floeimg.

add_passtimes!(props::DataFrame, passtimes::DateTime)
add_passtimes!.(props::Vector{DataFrame}, passtimes::Vector{DateTime})

Add a column passtime to each DataFrame in props containing the time of the image in which the floes were captured.

Arguments

• props: array of DataFrames containing floe properties.
• passtimes: array of DateTime objects containing the time of the image in which the floes were captured.

add_uuids!(df::DataFrame)
add_uuids!.(dfs::Vector{DataFrame})

Assign a unique ID to each floe in a (vector of) table(s) of floe properties.

add_ψs!(props_df::DataFrame})
add_ψs!.(props_dfs::Vector{DataFrame})

Add the ψ-s curves to each row of props_df.

Note: each member of props must have a mask column with a binary image representing the floe. To add floe masks see addfloemasks!.

Align images by padding so that the centroids of each image are on the edge of or within the same pixel.

buildψs(floe_mask)

Alternate method of buildψs accepting binary floe mask as input.

buildψs(XY::Matrix{<:Number};rangeout::Bool=true,
unwrap::Bool=true)

Alternate method of buildψs accepting input vectors x and y as a 2-column matrix [x y].

buildψs(x::Vector{<:Number},
           y::Vector{<:Number};
           rangeout::Bool=true,
           unwrap::Bool=true)::Tuple{Vector{Float64}, Vector{Float64}}

Builds the ψ-s curve defined by vectors x and y.

Returns a tuple of vectors with the phases ψ in the first component and the traversed arclength in the second component.

Following the convention in [1], the wrapped ψ-s curve has values in [0, 2π) by default; use rangeout to control this behavior.

See also bwtraceboundary, resample_boundary

Arguments

• x: Vector of x-coordinates
• y: corresponding vector of y-coordinates
• rangeout: true (default) for phase values in [0, 2π); false for phase values in (-π, π].
• unwrap: set to true to get "unwrapped" phases (default).

Reference

[1] McConnell, Ross, et al. "psi-s correlation and dynamic time warping: two methods for tracking ice floes in SAR images." IEEE Transactions on Geoscience and Remote sensing 29.6 (1991): 1004-1012.

Example

The example below builds a cardioid and obtains its ψ-s curve.

julia> t = range(0,2pi,201);

julia> x = @. cos(t)*(1-cos(t));

julia> y = @. sin(t)*(1-cos(t));

julia> plot(x,y) # visualize the cardioid

julia> psi, s = buildψs(x,y);

julia> [s psi] # inspect psi-s data
200×2 Matrix{Float64}:
 0.00049344  0.0314159
 0.0019736   0.0733034
 0.00444011  0.11938
 0.00789238  0.166055
 0.0123296   0.212929
 0.0177505   0.259894
 0.024154    0.306907
 0.0315383   0.35395
 0.0399017   0.401012
 0.0492421   0.448087
 ⋮
 7.96772     9.02377
 7.97511     9.07083
 7.98151     9.11787
 7.98693     9.16488
 7.99137     9.21185
 7.99482     9.25872
 7.99729     9.3054
 7.99877     9.35147
 7.99926     9.39336

 julia> plot(s, psi) # inspect psi-s curve -- should be a straight line from (0, 0) to (8, 3π)

bwtraceboundary(image::Union{Matrix{Int64},Matrix{Float64},T};
                P0::Union{Tuple{Int,Int},CartesianIndex{2},Nothing}=nothing,
                closed::Bool=true) where T<:AbstractMatrix{Bool}

Trace the boundary of objects in image

Background pixels are represented as zero. The algorithm traces the boundary counterclockwise and an initial point P0 can be specified. If more than one boundary is detected and an initial point is provided, the boundary that contains this point is returned as a vector of CartesianIndex types. Otherwise an array of vectors is returned with all the detected boundaries in image.

Arguments

• image: image, preferably binary with one single object, whose objects' boundaries are to be traced.
• P0: initial point of a target boundary.
• closed: if true (default) makes the inital point of a boundary equal to the last point.

Example

julia> A = zeros(Int, 13, 16); A[2:6, 2:6] .= 1; A[4:8, 7:10] .= 1; A[10:12,13:15] .= 1; A[10:12,3:6] .= 1;

julia> A
13×16 Matrix{Int64}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
 0  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
 0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

julia> boundary = IceFloeTracker.bwtraceboundary(A);

julia> boundary[3]
9-element Vector{CartesianIndex}:
 CartesianIndex(10, 13)
 CartesianIndex(11, 13)
 CartesianIndex(12, 13)
 CartesianIndex(12, 14)
 CartesianIndex(12, 15)
 CartesianIndex(11, 15)
 CartesianIndex(10, 15)
 CartesianIndex(10, 14)
 CartesianIndex(10, 13)

Calculate the centroid of a binary image. If 'rounded', return the nearest integer.

cropfloe(floesimg, props, i)

Crops the floe delimited by the bounding box data in props at index i from the floe image floesimg.

If the dataframe has bounding box data min_row, min_col, max_row, max_col, but no label, then returns the largest contiguous component.

If the dataframe has bounding box data min_row, min_col, max_row, max_col, and a label, then returns the component with the label. In this case, floesimg must be an Array{Int}.

If the dataframe has only a label and no bounding box data, then returns the component with the label, padded by one cell of zeroes on all sides. In this case, floesimg must be an Array{Int}.

cropfloe(floesimg, min_row, min_col, max_row, max_col)

Crops the floe delimited by min_row, min_col, max_row, max_col, from the floe image floesimg.

cropfloe(floesimg, min_row, min_col, max_row, max_col, label)

Crops the floe from floesimg with the label label, returning the region bounded by min_row, min_col, max_row, max_col, and converting to a BitMatrix.

r, lags = crosscorr(u::Vector{T},
                    v::Vector{T};
                    normalize::Bool=false,
                    padmode::Symbol=:longest)

Wrapper of DSP.xcorr with normalization (see https://docs.juliadsp.org/stable/convolutions/#DSP.xcorr)

Returns the pair (r, lags) with the cross correlation scores r and corresponding lags according to padmode.

Arguments

• u,v: real vectors which could have unequal length.
• normalize: return normalized correlation scores (false by default).
• padmode: either :longest (default) or :none to control padding of shorter vector with zeros.

Examples

The example below builds two vectors, one a shifted version of the other, and computes various cross correlation scores.

julia> n = 1:5;

julia> x = 0.48.^n;

julia> y = circshift(x,3);

julia> r, lags = crosscorr(x,y,normalize=true);

julia> [r lags]
9×2 Matrix{Float64}:
0.369648    -4.0
0.947531    -3.0
0.495695    -2.0
0.3231      -1.0
0.332519     0.0
0.15019      1.0
0.052469     2.0
0.0241435    3.0
0.00941878   4.0

julia> r, lags = crosscorr(x,y,normalize=true,padmode=:none);

julia> [r lags]
9×2 Matrix{Float64}:
0.369648    1.0
0.947531    2.0
0.495695    3.0
0.3231      4.0
0.332519    5.0
0.15019     6.0
0.052469    7.0
0.0241435   8.0
0.00941878  9.0

This final example builds two vectors of different length and computes some cross correlation scores.

julia> n = 1:5; m = 1:3;

julia> x = 0.48.^n; y = 0.48.^m;

julia> r, lags = crosscorr(x,y,normalize=true);

julia> [r lags]
9×2 Matrix{Float64}:
0.0          -4.0
4.14728e-17  -3.0
0.178468     -2.0
0.457473     -1.0
0.994189      0.0
0.477211      1.0
0.229061      2.0
0.105402      3.0
0.0411191     4.0

julia> r, lags = crosscorr(x,y,normalize=true,padmode=:none);

julia> [r lags]
7×2 Matrix{Float64}:
0.178468   1.0
0.457473   2.0
0.994189   3.0
0.477211   4.0
0.229061   5.0
0.105402   6.0
0.0411191  7.0

floe_tracker(props; filter_function, matching_function, minimum_floe_size, maximum_floe_size, maximum_time_step)

Track ice floes over multiple observations.

Trajectories are built as follows:

• Assume the floes detected in observation 1 are trajectories of length 1.
• For each subsequent observation at time t`:
  • Determine the latest observation for each trajectory – these are the "current trajectory heads".
  • Select the subset of trajectory heads observed within the window maximum_time_step, t
  • Apply the filter function in order to determine possible floe pairings
  • Apply the matching function to produce unique pairs of floes
  • Update the trajectories to include the newly paired floes
  • Add all unmatched floes as heads for new trajectories.

Arguments

• props::Vector{DataFrame}: A vector of DataFrames, each containing ice floe properties for a single observation time. Each DataFrame must have the following columns:
  • "area"
  • "min_row"
  • "min_col"
  • "max_row"
  • "max_col"
  • "row_centroid"
  • "col_centroid"
  • "convex_area"
  • "majoraxislength"
  • "minoraxislength"
  • "orientation"
  • "perimeter"
  • "mask": 2D boolean floe image cropped to the floe location
  • "passtime": A timestamp for the floe
  • "psi": the psi-s curve for the floe
  • "uuid": a universally unique identifier for each segmented floe
• filter_function: A function that accepts a floe::DataFrameRow and a candidates::DataFrame argument, and subsets the candidates dataframe to those rows that are possible matches for floe.
• matching_function: A function that takes the dataframe of candidate pairs and resolves conflicts to find at most one match for each floe.

Returns

A DataFrame with the above columns, plus extra columns:

• columns added by the filter function, such as similarity measures
• head_uuid, the floe which was best matched by this floe.
• Trajectories are identified by:
  • a unique identifier ID and the
  • UUID of the trajectory, trajectory_uuid.

Note: the props dataframes are modified in place.

Calculate the angle and rotation rate between two observations in DataFrameRows row1 and row2. image_column and time_column specify which columns to use from the DataFrameRows. registration_function is used to compare the two images and should return an angle. Returns a NamedTuple with the angle theta_rad, time difference dt_sec and rotation rate omega_rad_per_sec, and the two input rows.

Calculate the angle and rotation rate between a measurement in a DataFrameRow measurement, and all the other rows in DataFrame df.

• image_column is the column with the image to compare,
• time_column is the column with the timepoint of each observation,
• registration_function is used to compare the two images and should return an angle.

Returns a vector of NamedTuples with one entry for each comparison, with the angle theta_rad, time difference dt_sec and rotation rate omega_rad_per_sec, and the two input rows for each comparison row1 and row2.

Calculate the angle and rotation rate between observations in DataFrame df.

• id_column is the column with the ID of the image over several observations, e.g. the floe ID.
• image_column is the column with the image to compare,
• time_column is the column with the timepoint of each observation,
• registration_function is used to compare the two images and should return an angle.

Each row is compared to each other row in df which are:

• for the same object ID,
• strictly older,
• not older than the previous day.

Returns a DataFrame with one row for each comparison, with the angle theta_rad, time difference dt_sec and rotation rate omega_rad_per_sec, and all the other values from df with the column name suffix 1 for the first observation and 2 for the second.

Calculate the angle and rotation rate between two images image1 and image2 at times time1 and time2. Returns a NamedTuple with the angle theta_rad, time difference dt_sec and rotation rate omega_rad_per_sec. registration_function is used to compare the two images and should return an angle.

dx, dy = grad(A::Matrix{<:Number})

Make gradient vector field for the set of points with coordinates in the rows of the matrix A with x-coordinates down column 1 and y-coordinates down column 2. Return a tuple with dx and dy in that order.

dx, dy = grad(x::Vector{<:Number}, y::Vector{<:Number})

Make gradient vector field for the set of points with coordinates in vectors x and y. Return a tuple with dx and dy in that order.

mismatch(
    fixed::AbstractArray,
    moving::AbstractArray,
    mxrot::Real,
    step::Real,
)

Estimate a rotation that minimizes the 'mismatch' of aligning moving with fixed.

Returns a pair with the mismatch score mm and the associated registration angle rot.

Arguments

• fixed,moving: images to align via a rigid transformation
• mxrot: maximum rotation angle in degrees
• step: rotation angle step size in degrees

The default registration angles are evenly distributed in steps of 5º around a full rotation, ensuring that no angles are repeated (since -180º == +180º).

Angles are ordered so that smaller absolute angles which are positive will be returned in the event of a tie in the shape difference. ```

mismatch(
    fixed::AbstractArray,
    moving::AbstractArray,
    test_angles::AbstractArray,
)

Estimate a rotation that minimizes the 'mismatch' of aligning moving with fixed.

Returns a pair with the mismatch score mm and the associated registration angle rot.

Arguments

• fixed,moving: images to align via a rigid transformation
• test_angles: candidate angles to check for rotations by, in degrees. In the case of a tie in the shape difference, the earlier angle from this array will be returned.

```

norm(v)

Get the euclidean norm of the vector v.

normalized_cross_corr(f1,f2)

Return the normalized cross-correlation between the psi-s curves p1 and p2.

Finds the image rotation angle in test_angles which minimizes the shape difference between im_reference and im_target. The default test angles are shown in register_default_angles_rad. Use imrotate_function=imrotate_bin_<clockwise|counterclockwise>_<radians|degrees> to get angles <clockwise|counterclockwise> in <radians|degrees>.

resample_boundary(bd_points::Vector{<:CartesianIndex}, reduc_factor::Int64=2, bd::String="natural")

Get a uniform set of resampled boundary points from bd_points using cubic splines with specified boundary conditions

The resampled set of points is obtained using parametric interpolation of the points in bd_points. It is assumed that the separation between a pair of adjacent points is 1.

Arguments

• bd_points: Sequetial set of boundary points for the object of interest
• reduc_factor: Factor by which to reduce the number of points in bd_points (2 by default)

-bd: Boundary condition, either 'natural' (default) or 'periodic'

See also bwtraceboundary

Example

```jldoctest; setup = :(using IceFloeTracker) julia> A = zeros(Int, 13, 16); A[2:6, 2:6] .= 1; A[4:8, 7:10] .= 1; A[10:12,13:15] .= 1; A[10:12,3:6] .= 1; julia> A 13×16 Matrix{Int64}: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 1 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

julia> boundary = bwtraceboundary(A);

julia> boundary[3] 9-element Vector{CartesianIndex}: CartesianIndex(10, 13) CartesianIndex(11, 13) CartesianIndex(12, 13) CartesianIndex(12, 14) CartesianIndex(12, 15) CartesianIndex(11, 15) CartesianIndex(10, 15) CartesianIndex(10, 14) CartesianIndex(10, 13)

julia> resample_boundary(boundary[3]) 4×2 Matrix{Float64}: 10.0 13.0 12.0357 13.5859 10.5859 15.0357 10.0 13.0

Computes the shape difference between imreference and imtarget for each angle in testangles. The reference image is held constant, while the target image is rotated. The testangles are interpreted as the angle of rotation from target to reference, so to find the best match, we rotate the reverse direction. A perfect match at angle A would imply imtarget is the same shape as if imreference was rotated by A. Use imrotate_function=imrotate_bin_<clockwise|counterclockwise>_<radians|degrees> to get angles <clockwise|counterclockwise> in <radians|degrees>.

LopezAcostaCloudMask(prelim_threshold, band_7_threshold, band_2_threshold, ratio_lower, ratio_offset, ratio_upper)

AbstractCloudMaskAlgorithm implementation of the cloud mask from Lopez-Acosta et al. 2019. Cloud masks algorithms are initialized with a set of parameters, then can be supplied to create_cloudmask as an argument. The Lopez-Acosta et al. cloudmask creates a piecewise linear bifurcation of band 2 and band 7 brightness in a MODIS 7-2-1 false color image using a sequence of thresholds on band 2 and band 7 and on the ratio of band 7 to band 2 brightness.

Example:

using IceFloeTracker
using IceFloeTracker: Watkins2026Dataset

dataset = Watkins2026Dataset(; ref="v0.1")
case = first(filter(c -> (c.case_number == 6 && c.satellite == "terra"), dataset))
cm_algo = LopezAcostaCloudMask()
cloud_mask = create_cloudmask(modis_falsecolor(case), cm_algo)

# show image:
Gray.(cloud_mask)

Watkins2025CloudMask(prelimthreshold, band7threshold, band2threshold, ratiolower, ratiooffset, ratioupper, markerstrel, openingstrel)

Extension of the Lopez-Acosta et al. 2019 with parameters calibrated to the Ice Floe Validation Dataset. The Lopez-Acosta et al. cloudmask creates a piecewise linear bifurcation of band 2 and band 7 brightness in a MODIS 7-2-1 false color image using a sequence of thresholds on band 2 and band 7 and on the ratio of band 7 to band 2 brightness. This extension first creates a cloud mask using the LopezAcostaCloudMask, then applies morphological operations to remove speckle and smooth boundaries.

Example:

using IceFloeTracker
using IceFloeTracker: Watkins2026Dataset

dataset = Watkins2026Dataset(; ref="v0.1")
case = first(filter(c -> (c.case_number == 6 && c.satellite == "terra"), dataset))
cm_algo = Watkins2025CloudMask()
cloud_mask = create_cloudmask(modis_falsecolor(case), cm_algo)

# show image:
Gray.(cloud_mask)

apply_cloudmask(false_color_image, cloudmask)

Zero out pixels containing clouds where clouds and ice are not discernable. Arguments should be of the same size.

Arguments

• img: RGB, RGBA, or Gray image to be masked
• cloudmask: binary cloudmask with clouds = 1, else = 0
• modify_channel_1: optional keyword argument for RGB images. If true, set the first channel to 0 in the returned image.

apply_landmask(input_image, landmask_binary)

Zero out pixels in all channels of the input image using the binary landmask.

Arguments

• input_image: truecolor RGB image
• landmask_binary: binary landmask with 1=land, 0=water/ice

apply_landmask(img, landmask; as_indices::Bool=false)

Apply the landmask to the input image, optionally returning the indices of non-masked (ocean/ice) pixels.

Arguments

• img: input image (e.g., ice mask or RGB image)
• landmask: binary landmask (1=ocean/ice, 0=land)
• as_indices: if true, return indices of non-masked pixels; otherwise, return masked image

create_cloudmask(img, f::AbstractCloudMaskAlgorithm)

Cloud masks in the IFT are BitMatrix objects such that for an image I and cloudmask C, cloudy pixels can be selected by I[C], and clear-sky pixels can be selected with I[.!C]. Construction of a cloud mask uses the syntax

f = CloudMaskAlgorithm(parameters)
C = create_cloudmask(img; CloudMaskAlgorithm)

By default, create_cloudmask uses the algorithm found in [1]. This algorithm converts a 3-channel MODIS 7-2-1 false color image into a 1-channel binary matrix in which clouds = 1 and anything else = 0. The algorithm aims to identify patches of opaque cloud while allowing thin and transparent cloud to remain. This algorithm is instantiated using

f = LopezAcostaCloudMask()

In this case, the default values are applied. It can also called using a set of customized parameters. These values must be real numbers between 0 and 1. To reproduce the default parameters, you may call

f = LopezAcostaCloudMask(prelim_threshold=110/255, band_7_threshold=200/255, band_2_threshold=190/255, ratio_lower=0.0, ratio_upper=0.75).

A stricter cloud mask was defined in [2], covering more cloudy pixels while minimally impacting the masking of cloud-covered ice pixels.

f = LopezAcostaCloudMask(prelim_threshold=53/255, band_7_threshold=130/255, band_2_threshold=169/255, ratio_lower=0.0, ratio_upper=0.53).

These parameters together define a piecewise linear partition of pixels based on their Band 7 and Band 2 callibrated reflectance. Pixels with intensity above prelim_threshold are considered as potential cloudy pixels. Then, pixels with Band 7 reflectance less than band_7_threshold, Band 2 reflectance greater than band_2_threshold, and Band 7 to Band 2 ratios between ratio_lower and ratio_upper are removed from the cloud mask (i.e., set to cloud-free).

Arguments

• false_color_image: corrected reflectance false color image - bands [7,2,1]
• prelim_threshold: threshold value used to identify clouds in band 7, N0f8(RGB intensity/255)
• band_7_threshold: threshold value used to identify cloud-ice in band 7, N0f8(RGB intensity/255)
• band_2_threshold: threshold value used to identify cloud-ice in band 2, N0f8(RGB intensity/255)
• ratio_lower: threshold value used to set lower ratio of cloud-ice in bands 7 and 2
• ratio_upper: threshold value used to set upper ratio of cloud-ice in bands 7 and 2
• ratio_offset: offset value used to adjust the upper ratio of cloud-ice in bands 7 and 2

1. Lopez-Acosta, R., Schodlok, M. P., & Wilhelmus, M. M. (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), 1–15. (https://doi.org/10.1016/j.rse.2019.111406)[https://doi.org/10.1016/j.rse.2019.111406]
2. Watkins, D.M., Kim, M., Paniagua, C., Divoll, T., Holland, J.G., Hatcher, S., Hutchings, J.K., and Wilhelmus, M.M. (in prep). Calibration and validation of the Ice Floe Tracker algorithm.

create_landmask(landmask_image, struct_elem, fill_value_lower, fill_value_upper)

Convert a land mask image to a 1-channel binary matrix, and use a structuring element to extend a buffer to mask complex coastal features, and fill holes in the dilated image. In the resulting mask, land = 0 and ocean = 1. Returns a named tuple with the dilated and non-dilated landmask.

Arguments

• landmask_image: RGB land mask image from fetchdata
• struct_elem: structuring element for dilation (optional)
• fill_value_lower: fill holes having at least these many pixels (optional)
• fill_value_upper: fill holes having at most these many pixels (optional)

make_landmask_se()

Create a structuring element for dilating the landmask.

branch(img::AbstractArray{Bool})

Find branch points in skeletonized image img according to Definition 3 of [1].

[1] Arcelli, Carlo, and Gabriella Sanniti di Baja. "Skeletons of planar patterns." Machine Intelligence and Pattern Recognition. Vol. 19. North-Holland, 1996. 99-143.

bridge(bw)

Set 0-valued pixels to 1 if they have two nonzero neighbors that are not connected. Note the following exceptions:

0 0 0 0 0 0 1 0 1 becomes 1 1 1 0 0 0 0 0 0

1 0 1 1 1 1 0 0 0 becomes 0 0 0 0 0 0 0 0 0

The same applies to all their corresponding rotations.

Examples

julia> bw = [0 0 0; 0 0 0; 1 0 1]
3×3 Matrix{Int64}:
 0  0  0
 0  0  0
 1  0  1

julia> bridge(bw)
3×3 BitMatrix:
 0  0  0
 0  0  0
 1  1  1

julia> bw = [1 0 0; 1 0 1; 0 0 1]
3×3 Matrix{Int64}:
 1  0  0
 1  0  1
 0  0  1

julia> bridge(bw)
3×3 BitMatrix:
 1  1  0
 1  1  1
 0  1  1

 julia> bw = [1 0 1; 0 0 0; 1 0 1]
3×3 Matrix{Int64}:
 1  0  1
 0  0  0
 1  0  1

julia> bridge(bw)
3×3 BitMatrix:
 1  1  1
 1  1  1
 1  1  1

bwareamaxfilt(bwimg::AbstractArray{Bool}, conn)

Filter the smaller (by area) connected components in bwimg keeping the (assumed unique) largest.

Uses 8-pixel connectivity by default (conn=8). Use conn=4 for 4-pixel connectivity.

bwareamaxfilt!(bwimg::AbstractArray)

In-place version of bwareamaxfilt.

See also bwareamaxfilt

bwperim(bwimg)

Locate the pixels at the boundary of objects in an binary image bwimg using 8-pixel connectivity.

Arguments

• bwimg: Binary (black/white – 1/0) image

Examples

julia> A = zeros(Bool, 13, 16); A[2:6, 2:6] .= 1; A[4:8, 7:10] .= 1; A[10:12,13:15] .= 1; A[10:12,3:6] .= 1;

julia> A
13×16 Matrix{Bool}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
 0  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
 0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

 julia> bwperim(A)
13×16 Matrix{Bool}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  1  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  1  1  1  1  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  1  0  0  0  0  0  0
 0  1  1  1  1  1  0  0  0  1  0  0  0  0  0  0
 0  0  0  0  0  0  1  0  0  1  0  0  0  0  0  0
 0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

filt_except_label(labeled_arr::Array{Int64, 2}, label::Int64)

Make 0 all values in labeled_arr that are not equal to label.

See also filt_except_label!

get_areas(labeled_arr::Array{T, 2})::Dict{T, Int} where T

Get the "areas" (count of pixels of a given label) of the connected components in labeled_arr.

Return a dictionary with the frequency distribution: label => countoflabel.

get_max_label(d::Dict{Int, Int})

Get the label k in dictionary d for which d[k] is maximal.

hbreak!(img::AbstractArray{Bool})

Inplace version of hbreak. See also hbreak.

hbreak(img::AbstractArray{Bool})

Remove H-connected pixels in the binary image img. See also hbreak! for an inplace version of this function.

Examples

julia> h1 = trues(3,3); h1[[1 3], 2] .= false; h1     
3×3 BitMatrix:    
 1  0  1
 1  1  1
 1  0  1

julia> h2 = trues(3,3); h2[2, [1 3]] .= false; h2     
3×3 BitMatrix:    
 1  1  1
 0  1  0
 1  1  1

julia> hbreak!(h1); h1 # modify h1 inplace
3×3 BitMatrix:
 1  0  1
 1  0  1
 1  0  1

julia> hbreak(h2) 
3×3 BitMatrix:
 1  1  1
 0  0  0
 1  1  1

imextendedmin(img)

Mimics MATLAB's imextendedmin function that computes the extended-minima transform, which is the regional minima of the H-minima transform. Regional minima are connected components of pixels with a constant intensity value. This function returns a transformed bitmatrix.

Arguments

• img: image object
• h: suppress minima below this depth threshold
• conn: neighborhood connectivity; in 2D 1 = 4-neighborhood and 2 = 8-neighborhood

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.

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

make_hbreak_dict()

Build dict with the two versions of an H-connected 3x3 neighboorhood.

h1 = [1 0 1 1 1 1 1 0 1]

h2 = [1 1 1 0 1 0 1 1 1]

morph_fill(bw::T)::T where {T<:AbstractArray{Bool}}

Fill holes in binary image bw by setting 0-valued pixels to 1 if they are surrounded by 1-valued pixels.

Examples

julia> bw = Bool[
        0 0 0 0 0
        0 1 1 1 0
        0 1 0 1 0
        0 1 1 1 0
        0 0 0 0 0
    ];

julia> morph_fill(bw)
5×5 Matrix{Bool}:
 0  0  0  0  0
 0  1  1  1  0
 0  1  1  1  0
 0  1  1  1  0
 0  0  0  0  0

reconstruct(img, se, type, invert)

Perform closing/opening by reconstruction on img.

Arguments

• img::AbstractArray: The input image.
• se::AbstractArray: The structuring element.
• type::String: The type of morphological operation to perform. Must be either "dilation" (close by reconstruction) or "erosion" (open by reconstruction).
• invert::Bool=true: Invert marker and mask before reconstruction.

conditional_histeq(
    image,
    clouds_red,
    tiles;
    entropy_threshold::Real=4.0,
    white_threshold::Real=25.5,
    white_fraction_threshold::Real=0.4,
)

Performs conditional histogram equalization on a true color image.

Arguments

• image: The true color image to be equalized.
• clouds_red: The land/cloud masked red channel of the false color image.
• tiles: the output from get_tiles(image) specifying the tiling to use on the image.
• entropy_threshold: The entropy threshold used to determine if a block should be equalized. Default is 4.0.
• white_threshold: The white threshold used to determine if a pixel should be considered white. Default is 25.5.
• white_fraction_threshold: The white fraction threshold used to determine if a block should be equalized. Default is 0.4.

Returns

The equalized true color image.

histeq(img)
histeq(img; nbins=64)

Histogram equalization of img using nbins bins.

imadjust(img; low, high)

Adjust the contrast of an image using linear stretching. The image is normalized to [0, 1] and then stretched to the range [low, high].

Arguments

• img: The input image.
• low: The lower bound of the stretched image. Default is 0.01.
• high: The upper bound of the stretched image. Default is 0.99.

Returns

The contrast-adjusted image in the range [0, 255].

imgradientmag(img)

Compute the gradient magnitude of an image using the Sobel operator.

Perform nonlinear diffusion on an input image. By default, use the Perona-Malik method.

rgb2gray(rgbchannels::Array{Float64, 3})

Convert an array of RGB channel data to grayscale in the range [0, 255].

Identical to MATLAB rgb2gray (https://www.mathworks.com/help/matlab/ref/rgb2gray.html).

rgb2gray(img::Matrix{RGB{Float64}})

Convert an RGB image to grayscale in the range [0, 255].

unsharp_mask(img, radius, amount, threshold)

Enhance image sharpness by weighted differencing of the image and a Gaussian blurred image.
If ``B`` is the blurred version of image ``I``, then an unsharp mask sharpened image is obtained by
``S = I + (I - B)*A``
The amount of sharpening is determined by the factor A. An option threshold can be supplied such
that the sharpening is only applied where ``I - B`` is greater than some factor.

# Arguments
img: input image
radius: standard deviation of the Gaussian blur
amount: multiplicative factor
threshold: minimum difference for applying the sharpening

# Returns
Sharpened image

unsharp_mask(image_gray, smoothing_param, intensity, clampmax)

Apply unsharp masking on (equalized) grayscale ([0, clampmax]) image to enhance its sharpness.

Arguments

• image_gray: The input grayscale image, typically already equalized.
• smoothing_param::Int: The pixel radius for Gaussian blurring (typically between 1 and 10).
• intensity: The amount of sharpening to apply. Higher values result in more pronounced sharpening.
• clampmax: upper limit of intensity values in the returned image.`

Returns

The sharpened grayscale image with values clipped between 0 and clapmax.

Module for loading validated ice floe data.

Watkins2026Dataset()
Watkins2026Dataset(; [ref, url, dataset_metadata_path, cache_dir])

Validated ice floe data from the Watkins 2026 Ice Floe Validation Dataset.

The dataset is initialized with a specific git tag, branch or commit ID from which to load the data.

julia> dataset = Watkins2026Dataset(; ref="v0.1")

Watkins2026Dataset fields:

• ref (optional): git tag, commit-id or branch from which to load the data
• cache_dir (optional): local path where the data will be stored, which defaults to /tmp/Watkins2026/.
• url (optional): URL of the GitHub repository with the dataset
• dataset_metadata_path (optional): path within the repository to a CSV file describing the data

The dataset can be filtered using Base.filter or DataFrames.subset. This checks each case's info, and if the function returns true, the case is included in the returned dataset.

julia> dataset = filter(c -> (
                        c.visible_floes == "yes" &&
                        c.cloud_category_manual == "none" &&
                        c.artifacts == "no"
                    ), dataset);

Equivalently:

julia> dataset = subset(dataset, 
                        :visible_floes => c -> c .== "yes",
                        :cloud_category_manual => c -> c .== "none",
                        :artifacts => c -> c .== "no",
                    );

The returned dataset (a Dataset) has an info accessor which returns a DataFrame of the cases which passed the filter:

julia> info(dataset)
8×28 DataFrame
    Row │ case_number  region        start_date  center_lon  center_lat  center_x  center_y  month  sea_ice_fr ⋯
        │ Int64        String31      Date        Float64     Float64     Int64     Int64     Int64  Float64    ⋯
   ─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────
      1 │          11  baffin_bay    2011-07-02    -70.7347     72.3303   -837500  -1737500      7             ⋯
      2 │          14  baffin_bay    2022-07-06    -69.0755     72.3157   -787500  -1762500      7
      3 │          48  beaufort_sea  2021-04-27   -140.612      70.1346  -2162500    212500      4
      4 │          48  beaufort_sea  2021-04-27   -140.612      70.1346  -2162500    212500      4
      5 │          54  beaufort_sea  2015-05-16   -136.675      70.4441  -2137500     62500      5             ⋯
      6 │          54  beaufort_sea  2015-05-16   -136.675      70.4441  -2137500     62500      5
      7 │         128  hudson_bay    2019-04-15    -91.9847     57.853   -2612500  -2437500      4
      8 │         166  laptev_sea    2016-09-04    136.931      79.7507    -37500   1112500      9
                                                                                        20 columns omitted

The dataset contains Case objects. Each Case has info fields including:

• name: name of the case
• info: metadata for the case, corresponding to a row in the dataset.info DataFrame

Each Case has functions to access its contents:

• modis_truecolor: MODIS true color image
• modis_falsecolor: MODIS false color image
• modis_landmask: MODIS landmask image
• modis_cloudfraction: MODIS cloud fraction image
• masie_landmask: MASIE landmask image
• masie_seaice: MASIE sea ice image
• validated_binary_floes: binary image of validated floes
• validated_labeled_floes: labeled image of validated floes
• validated_floe_properties: CSV file of validated floe properties

The dataset can be iterated over to get each Case: Example:

julia> for case in dataset
           println(name(case) * 
                   ": sea ice fraction: " * string(info(case).sea_ice_fraction) *
                   ", true color image size: " * string(size(modis_truecolor(case))))
       end
011-baffin_bay-100km-20110702-aqua-250m: sea ice fraction: 0.8, true color image size: (400, 400)
014-baffin_bay-100km-20220706-terra-250m: sea ice fraction: 1.0, true color image size: (400, 400)
048-beaufort_sea-100km-20210427-aqua-250m: sea ice fraction: 1.0, true color image size: (400, 400)
048-beaufort_sea-100km-20210427-terra-250m: sea ice fraction: 1.0, true color image size: (400, 400)
054-beaufort_sea-100km-20150516-aqua-250m: sea ice fraction: 1.0, true color image size: (400, 400)
054-beaufort_sea-100km-20150516-terra-250m: sea ice fraction: 1.0, true color image size: (400, 400)
128-hudson_bay-100km-20190415-aqua-250m: sea ice fraction: 1.0, true color image size: (400, 400)
166-laptev_sea-100km-20160904-aqua-250m: sea ice fraction: 1.0, true color image size: (400, 400)

To get the first case in the dataset, you can use first(...):

julia> first(dataset)
Case(GitHubLoader("https://github.com/danielmwatkins/ice_floe_validation_dataset/", "v0.1", "/tmp/Watkins2026"), DataFrameRow
 Row │        case_number  region      start_date  center_lon  center_lat  center_x  center_y  month  sea_ice_fraction  mean_sea_ice_concentration  init_case_number  satellite  visible_sea_ice  visible_la ⋯
     │ Int64  Int64        String      Dates.Date  Float64     Float64     Int64     Int64     Int64  Float64           Float64                     Int64             String     String           String     ⋯
─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │    20           11  baffin_bay  2011-07-02    -70.7347     72.3303   -837500  -1737500      7               0.8                       0.388                -1  aqua       yes              no         ⋯
                                                                                                                                                                                            16 columns omitted)

julia> validated_labeled_floes(first(dataset))
Segmented Image with:
labels map: 400×400 Matrix{Int64}
number of labels: 105

julia> modis_truecolor(first(dataset))
400×400 Array{RGBA{N0f8},2} with eltype ColorTypes.RGBA{FixedPointNumbers.N0f8}:
 RGBA{N0f8}(0.094,0.133,0.169,1.0)  RGBA{N0f8}(0.051,0.094,0.118,1.0) ...

```

latlon(imgpath::AbstractString)

Reads the GeoTiff located at <imgpath>, extracts the coordinate reference system, and produces a lookup table with for the column and row values in the same projection as the GeoTiff and a 2D array for latitude and longitude.

Module for general image utilities.

bump_tile(tile::Tuple{UnitRange{Int64}, UnitRange{Int64}}, dims::Tuple{Int,Int})::Tuple{UnitRange{Int}, UnitRange{Int}}

Adjust the tile dimensions by adding extra rows and columns.

Arguments

• tile::Tuple{Int,Int,Int,Int}: A tuple representing the tile dimensions (a, b, c, d).
• dims::Tuple{Int,Int}: A tuple representing the extra rows and columns to add (extrarows, extracols).

Returns

• Tuple{UnitRange{Int}, UnitRange{Int}}: A tuple of ranges representing the new tile dimensions.

Examples

```julia julia> bump_tile((1:3, 1:4), (1, 1)) (1:4, 1:5)

get_area_missed(side_length::Int, dims::Tuple{Int,Int})::Float64

Calculate the proportion of the area that is not covered by tiles of a given side length.

Arguments

• side_length::Int: The side length of the tile.
• dims::Tuple{Int,Int}: A tuple representing the dimensions (width, height).

Returns

• Float64: The proportion of the area that is not covered by the tiles.

Examples

``` julia> getareamissed(5, (10, 20)) 0.0

julia> getareamissed(7, (10, 20)) 0.51

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 equalizedgrayreconstructedimg - gammagreen clamped between 0 and 255.

get_optimal_tile_size(l0::Int, dims::Tuple{Int,Int}) -> Int

Calculate the optimal tile size in the range [l0-1, l0+1] for the given size l0 and image dimensions dims.

Description

This function computes the optimal tile size for tiling an area with given dimensions. It ensures that the initial tile size l0 is at least 2 and not larger than any of the given dimensions. The function evaluates candidate tile sizes and selects the one that minimizes the area missed by its corresponding tiling. In case of a tie, it prefers the larger tile size.

Example

julia> get_optimal_tile_size(3, (10, 7))
2

get_tile_dims(tile)

Calculate the dimensions of a tile.

Arguments

• tile::Tuple{UnitRange{Int},UnitRange{Int}}: A tuple representing the tile dimensions.

Returns

• Tuple{Int,Int}: A tuple representing the width and height of the tile.

Examples

```julia julia> gettiledims((1:3, 1:4)) (4, 3)

get_tile_meta(tile)

Extracts metadata from a given tile.

Arguments

• tile: A collection of tuples, where each tuple represents a coordinate pair.

Returns

• A tuple (a, b, c, d) where:
  • a: The first element of the first tuple in tile.
  • b: The last element of the first tuple in tile.
  • c: The first element of the last tuple in tile.
  • d: The last element of the last tuple in tile.

get_tiles(array, side_length)

Generate a collection of tiles from an array.

Unlike TileIterator, the function adjusts the bottom and right edges of the tile matrix if they are smaller than half the tile size side_length.

get_tiles(array, t::Tuple{Int,Int})

Generate a collection of tiles from an array.

The function adjusts the bottom and right edges of the tile matrix if they are smaller than half the tile sizes in t.

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.

masker(mask::AbstractArray, img::AbstractArray{<:Colorant})
masker(mask::AbstractArray)

Returns a version of img with masked pixels made transparent. If img has an alpha channel, it is combined with the mask.

masker(mask) returns a function which can be used to apply the mask.

Examples

With a BitMatrix type of mask, truthy values in the mask are transparent in the output.

julia> using Images

julia> hide = true;

julia> pass = false;

julia> bit_mask = [hide hide pass; pass pass hide]
2×3 Matrix{Bool}:
 1  1  0
 0  0  1

 julia> img = parse.(Colorant, ["red" "green" "blue"; "cyan" "magenta" "yellow"])
2×3 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(1.0,0.0,0.0)  RGB{N0f8}(0.0,0.502,0.0)  RGB{N0f8}(0.0,0.0,1.0)
 RGB{N0f8}(0.0,1.0,1.0)  RGB{N0f8}(1.0,0.0,1.0)    RGB{N0f8}(1.0,1.0,0.0)

 julia> masker(bit_mask, img)
2×3 Array{ARGB{N0f8},2} with eltype ARGB{N0f8}:
 ARGB{N0f8}(1.0,0.0,0.0,0.0)  ARGB{N0f8}(0.0,0.502,0.0,0.0)  ARGB{N0f8}(0.0,0.0,1.0,1.0)
 ARGB{N0f8}(0.0,1.0,1.0,1.0)  ARGB{N0f8}(1.0,0.0,1.0,1.0)    ARGB{N0f8}(1.0,1.0,0.0,0.0)

Using the masker in a pipeline is also possible:

julia> img |> masker(bit_mask)
2×3 Array{ARGB{N0f8},2} with eltype ARGB{N0f8}:
 ARGB{N0f8}(1.0,0.0,0.0,0.0)  ARGB{N0f8}(0.0,0.502,0.0,0.0)  ARGB{N0f8}(0.0,0.0,1.0,1.0)
 ARGB{N0f8}(0.0,1.0,1.0,1.0)  ARGB{N0f8}(1.0,0.0,1.0,1.0)    ARGB{N0f8}(1.0,1.0,0.0,0.0)

Where the mask is itself an image with transparency, areas which are opaque in the mask are transparent in the output. This corresopnds to overlaying the mask over the image, and hiding in the output those areas which were masked.

julia> hide = AGray(0.5, 1);

julia> pass = AGray(1, 0);

julia> agray_mask = [hide hide pass; pass pass hide];

julia> masker(agray_mask, img)
2×3 Array{ARGB{N0f8},2} with eltype ARGB{N0f8}:
 ARGB{N0f8}(1.0,0.0,0.0,0.0)  ARGB{N0f8}(0.0,0.502,0.0,0.0)  ARGB{N0f8}(0.0,0.0,1.0,1.0)
 ARGB{N0f8}(0.0,1.0,1.0,1.0)  ARGB{N0f8}(1.0,0.0,1.0,1.0)    ARGB{N0f8}(1.0,1.0,0.0,0.0)

Any pixels which are partially opaque in the mask, will be partially obscured in the output:

julia> part = AGray(0.75, 0.5)
AGray{Float64}(0.75,0.5)

julia> partial_mask = [hide part pass; pass part hide]
2×3 Array{AGray{Float64},2} with eltype AGray{Float64}:
 AGray{Float64}(0.5,1.0)  AGray{Float64}(0.75,0.5)  AGray{Float64}(1.0,0.0)
 AGray{Float64}(1.0,0.0)  AGray{Float64}(0.75,0.5)  AGray{Float64}(0.5,1.0)

julia> masker(partial_mask, img)
2×3 Array{ARGB{N0f8},2} with eltype ARGB{N0f8}:
 ARGB{N0f8}(1.0,0.0,0.0,0.0)  ARGB{N0f8}(0.0,0.502,0.0,0.502)  ARGB{N0f8}(0.0,0.0,1.0,1.0)
 ARGB{N0f8}(0.0,1.0,1.0,1.0)  ARGB{N0f8}(1.0,0.0,1.0,0.502)    ARGB{N0f8}(1.0,1.0,0.0,0.0)

If the image already has transparency, this is combined with the mask.

julia> imga = RGBA.(parse.(Colorant, ["red" "transparent" "blue"; "cyan" "transparent" "yellow"]))

2×3 Array{RGBA{N0f8},2} with eltype RGBA{N0f8}:
 RGBA{N0f8}(1.0,0.0,0.0,1.0)  RGBA{N0f8}(0.0,0.0,0.0,0.0)  RGBA{N0f8}(0.0,0.0,1.0,1.0)
 RGBA{N0f8}(0.0,1.0,1.0,1.0)  RGBA{N0f8}(0.0,0.0,0.0,0.0)  RGBA{N0f8}(1.0,1.0,0.0,1.0)

julia> masker(agray_mask, imga)
2×3 Array{ARGB{N0f8},2} with eltype ARGB{N0f8}:
 ARGB{N0f8}(1.0,0.0,0.0,0.0)  ARGB{N0f8}(0.0,0.0,0.0,0.0)  ARGB{N0f8}(0.0,0.0,1.0,1.0)
 ARGB{N0f8}(0.0,1.0,1.0,1.0)  ARGB{N0f8}(0.0,0.0,0.0,0.0)  ARGB{N0f8}(1.0,1.0,0.0,0.0)

Where the mask is an image without transparency, any non-zero pixels are masked:

julia> gray_mask = Gray.([0.5 0.1 0.0; 0.0 0.0 1.0])
2×3 Array{Gray{Float64},2} with eltype Gray{Float64}:

 Gray{Float64}(0.5)  Gray{Float64}(0.1)  Gray{Float64}(0.0)
 Gray{Float64}(0.0)  Gray{Float64}(0.0)  Gray{Float64}(1.0)

julia> masker(gray_mask, img)
2×3 Array{ARGB{N0f8},2} with eltype ARGB{N0f8}:
 ARGB{N0f8}(1.0,0.0,0.0,0.0)  ARGB{N0f8}(0.0,0.502,0.0,0.0)  ARGB{N0f8}(0.0,0.0,1.0,1.0)
 ARGB{N0f8}(0.0,1.0,1.0,1.0)  ARGB{N0f8}(1.0,0.0,1.0,1.0)    ARGB{N0f8}(1.0,1.0,0.0,0.0)

julia> rgb_mask = parse.(Colorant, ["red" "red" "black"; "black" "purple" "green"])
2×3 Array{RGB{N0f8},2} with eltype RGB{N0f8}:
 RGB{N0f8}(1.0,0.0,0.0)  RGB{N0f8}(1.0,0.0,0.0)      RGB{N0f8}(0.0,0.0,0.0)
 RGB{N0f8}(0.0,0.0,0.0)  RGB{N0f8}(0.502,0.0,0.502)  RGB{N0f8}(0.0,0.502,0.0)

julia> masker(rgb_mask, img)
2×3 Array{ARGB{N0f8},2} with eltype ARGB{N0f8}:
 ARGB{N0f8}(1.0,0.0,0.0,0.0)  ARGB{N0f8}(0.0,0.502,0.0,0.0)  ARGB{N0f8}(0.0,0.0,1.0,1.0)
 ARGB{N0f8}(0.0,1.0,1.0,1.0)  ARGB{N0f8}(1.0,0.0,1.0,0.0)    ARGB{N0f8}(1.0,1.0,0.0,0.0)

Where the mask is an array of Real numbers, 0-pixels will be completely unmasked, 1-pixels completely masked, and values between partially masked:

julia> real_mask = [1.0 0.5 0.1; 0.1 0.2 1.0]

2×3 Matrix{Float64}:
 1.0  0.5  0.1
 0.1  0.2  1.0

julia> masker(real_mask, img)
2×3 Array{ARGB{N0f8},2} with eltype ARGB{N0f8}:
 ARGB{N0f8}(1.0,0.0,0.0,0.0)    ARGB{N0f8}(0.0,0.502,0.0,0.502)  ARGB{N0f8}(0.0,0.0,1.0,0.902)
 ARGB{N0f8}(0.0,1.0,1.0,0.902)  ARGB{N0f8}(1.0,0.0,1.0,0.8)      ARGB{N0f8}(1.0,1.0,0.0,0.0)```

Where values are outside of the range [0, 1], 
they are clamped to whichever of 0 and 1 is nearer:

julia-repl julia> outofrange_mask = [5 2 0.75; -1 -2 1] 2×3 Matrix{Float64}: 5.0 2.0 0.75 -1.0 -2.0 1.0

julia> masker(outofrange_mask, img) 2×3 Array{ARGB{N0f8},2} with eltype ARGB{N0f8}: ARGB{N0f8}(1.0,0.0,0.0,0.0) ARGB{N0f8}(0.0,0.502,0.0,0.0) ARGB{N0f8}(0.0,0.0,1.0,0.251) ARGB{N0f8}(0.0,1.0,1.0,1.0) ARGB{N0f8}(1.0,0.0,1.0,1.0) ARGB{N0f8}(1.0,1.0,0.0,0.0)```

callable_store()

Create a store and a callback function to add key-value pairs to the store.

Returns a store::Dict and a callback::Function which stores any kwargs passed to it in the store.

Examples

Basic usage is to store values using the callback function

julia> store, callback = callable_store()
julia> store
Dict{Any, Any}()

julia> callback(;foo="bar")  # echoes the updated store
Dict{Any, Any} with 1 entry:
  :foo => "bar"

julia> store  # values are available from the store object
Dict{Any, Any} with 1 entry:
  :foo => "bar"

A real-world use case is to extract data from a segmentation algorithm run:

julia> intermediate_results, intermediate_results_callback = callable_store()
julia> data = first(Watkins2026Dataset(; ref="v0.1")());
julia> segments = LopezAcosta2019Tiling.Segment()(
    modis_truecolor(data),
    modis_falsecolor(data),
    modis_landmask(data);
    intermediate_results_callback,
)
Segmented Image with:
  labels map: 400×400 Matrix{Int64}
  number of labels: 12

julia> intermediate_data
Dict{Any, Any} with 16 entries:
  :binarized_tiling                       => Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0]
  :icemask                                => Bool[1 1 … 1 1; 1 1 … 1 1; … ; 0 0 … 1 1; 0 0 … 1 1]
  :equalized_gray                         => [0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0]
  :morphed_residue                        => [0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0]
  :L0mask                                 => Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0]
  :segmented                              => Segmented Image with:…
  :prelim_icemask2                        => [255 255 … 255 255; 255 255 … 255 255; … ; 255 255 … 255 255; 255 255 … 255 255]
  :equalized_gray_sharpened_reconstructed => [0 0 … 0 0; 0 0 … 0 0; … ; 255 255 … 255 255; 255 255 … 255 255]
  :gammagreen                             => [190.35 190.23 … 182.93 185.03; 191.68 190.6 … 185.04 192.08; … ; 163.87 173.33 … 108.02 108.18; 166.14 173.3 … 112.35 112.32]
  :segment_mask                           => Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0]
  :ref_img_cloudmasked                    => RGB{N0f8}[RGB{N0f8}(0.0,0.0,0.0) RGB{N0f8}(0.0,0.0,0.0) … RGB{N0f8}(0.008,0.706,0.761) RGB{N0f8}(0.0,0.722,0.769); RGB{N0f8}(0.0,0.0,0.0) RGB{N0f8}(0.0,0.0,0.0) … RGB{N0f8}(0.039,0.702,0.784) RGB{N0f8}(0.075,0.784,0.859); … ; RGB{…
  :prelim_icemask                         => Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0]
  :equalized_gray_reconstructed           => [0 0 … 0 0; 0 0 … 0 0; … ; 255 255 … 255 255; 255 255 … 255 255]
  :final                                  => Bool[0 0 … 0 0; 0 1 … 1 0; … ; 0 0 … 1 0; 0 0 … 0 0]
  :local_maxima_mask                      => [255 255 … 255 255; 255 255 … 255 255; … ; 255 255 … 255 255; 255 255 … 255 255]
  :labeled                                => [0 0 … 0 0; 0 1 … 1 0; … ; 0 0 … 9 0; 0 0 … 0 0]

@persist img fname
@persist(img,fname)
@persist img
@persist(img)
@persist img fname ts
@persist(img, fname, ts)

Given a reference to an image object img, the macro persists (saves to a file) img to the current working directory using fname as filename. Returns img.

Arguments

• img: Symbol expression representing an image object loaded in memory.
• fname: Optional filename for the persisted image.
• ts: Optional boolean to attach timestamp to fname.