diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index ff80a898..b0cfcc8c 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -15,21 +15,23 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        version:
-          - '1.9'
-        os:
-          - ubuntu-latest
-        arch:
-          - x64
+        version: ['1.9', '1.10', '~1.11.0-0']
+        os: [ubuntu-latest, macos-latest]
+        arch: [x64, arm64]
+        exclude:
+          - os: ubuntu-latest
+            arch: arm64
+          - os: macos-latest
+            arch: x64
     steps:
       - uses: actions/checkout@v4
       - uses: actions/setup-python@v5
         with:
-          python-version: '3.10'
+          python-version: '3.11'
       - name: update pip
         run: python -m pip install -U pip
       - name: install python deps
-        run: python -m pip install -U numpy==1.23 scikit-image==0.20.0 pyproj==3.6.0 rasterio==1.3.7 requests==2.31.0 skyfield==1.45.0 pandas==2 jinja2==3.1
+        run: python -m pip install -U -r requirements.txt
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 00000000..d42bbdb8
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,8 @@
+jinja2==3.1
+numpy==1.23.2
+pandas==2.1
+pyproj==3.6.0
+rasterio==1.3.7
+scikit-image==0.20.0
+requests==2.31.0
+skyfield==1.45.0
\ No newline at end of file
diff --git a/src/IceFloeTracker.jl b/src/IceFloeTracker.jl
index 9edd0590..fb0a1ede 100644
--- a/src/IceFloeTracker.jl
+++ b/src/IceFloeTracker.jl
@@ -1,21 +1,23 @@
 module IceFloeTracker
-using Images
-using DelimitedFiles: readdlm, writedlm
+using Clustering
+using DSP
+using DataFrames
 using Dates
+using DelimitedFiles: readdlm, writedlm
 using ImageContrastAdjustment
 using ImageSegmentation
+using Images
+using Interpolations
+using OffsetArrays: centered
 using Peaks
 using Pkg
+using PyCall
 using Random
+using Serialization: deserialize, serialize
 using StatsBase
-using Interpolations
-using DataFrames
-using PyCall
-using Clustering
-using DSP
 using StaticArrays
-using OffsetArrays: centered
-using Serialization: serialize, deserialize
+using StatsBase
+using TiledIteration
 using TOML
 
 export readdlm,
@@ -66,8 +68,12 @@ include("hbreak.jl")
 include("bridge.jl")
 include("branch.jl")
 include("special_strels.jl")
+include("tilingutils.jl")
+include("histogram_equalization.jl")
+
 
 const sk_measure = PyNULL()
+const sk_exposure = PyNULL()
 const getlatlon = PyNULL()
 
 function get_version_from_toml(pth=dirname(dirname(pathof(IceFloeTracker))))::VersionNumber
@@ -77,14 +83,42 @@ end
 
 const IFTVERSION = get_version_from_toml()
 
+function parse_requirements(file_path)
+    requirements = Dict{String,String}()
+    open(file_path, "r") do f
+        for line in eachline(f)
+            if occursin("==", line)
+            pkg, version = split(line, "==")
+            elseif occursin("=", line)
+            pkg, version = split(line, "=")
+            else
+            pkg, version = line, ""
+            end
+            requirements[pkg] = version
+        end
+    end
+    return requirements
+end
+
 function __init__()
-    pyimport_conda("numpy", "numpy=1.23")
-    pyimport_conda("pyproj", "pyproj=3.6.0")
-    pyimport_conda("rasterio", "rasterio=1.3.7")
-    pyimport_conda("jinja2", "jinja2=3.1.2")
-    pyimport_conda("pandas", "pandas=2")
+
+    deps = parse_requirements(joinpath(dirname(@__DIR__), "requirements.txt"))
+
+    for (pkg, version) in deps
+        if pkg == "scikit-image"
+            _modules = ["measure", "exposure"]
+            for _module in _modules
+                imported_module = pyimport_conda("skimage.$_module", "$(pkg)=$(version)")
+                reference_module = eval(Symbol("sk_$_module"))
+                copy!(reference_module, imported_module)
+            end
+
+        else
+            pyimport_conda(pkg, "$(pkg)=$(version)")
+        end
+    end
+
     @pyinclude(joinpath(@__DIR__, "latlon.py"))
-    copy!(sk_measure, pyimport_conda("skimage.measure", "scikit-image=0.20.0"))
     copy!(getlatlon, py"getlatlon")
     return nothing
 end
diff --git a/src/anisotropic_image_diffusion.jl b/src/anisotropic_image_diffusion.jl
index 1fc8640f..27f66d4a 100644
--- a/src/anisotropic_image_diffusion.jl
+++ b/src/anisotropic_image_diffusion.jl
@@ -1,19 +1,9 @@
-# Anisotropic Image Diffusion ##
-# This script is borrowed from https://github.com/Red-Portal/UltrasoundDesignGallery.jl ##
-## MIT license with permission to use 
-
-macro swap!(a::Symbol, b::Symbol)
-    blk = quote
-        c = $(esc(a))
-        $(esc(a)) = $(esc(b))
-        $(esc(b)) = c
-    end
-    return blk
-end
+# Anisotropic Image Diffusion
+# Adapted from https://github.com/Red-Portal/UltrasoundDesignGallery.jl
+## MIT license with permission to use
 
 function pmad_kernel!(image, output, g, λ)
-    M = size(image, 1)
-    N = size(image, 2)
+    M, N = size(image)
 
     @inbounds for j in 1:N
         @simd for i in 1:M
@@ -38,34 +28,36 @@ function pmad_kernel!(image, output, g, λ)
     end
 end
 
-function invert_color(color::RGB{Float64})
+function invert_color(color::RGB{T}) where {T<:AbstractFloat}
     return RGB(1.0 / color.r, 1.0 / color.g, 1.0 / color.b)
 end
-function invert_color(color::Gray{Float64})
+
+function invert_color(color::Gray{T}) where {T<:AbstractFloat}
     return Gray(1.0 / color.val)
 end
+
 function diffusion(
     image::Matrix{T}, λ::Float64, K::Int, niters::Int
 ) where {T<:Color{Float64}}
     #=
-        Perona, Pietro, and Jitendra Malik. 
-        "Scale-space and edge detection using anisotropic diffusion." 
+        Perona, Pietro, and Jitendra Malik.
+        "Scale-space and edge detection using anisotropic diffusion."
         IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), 1990.
     =#
-    if !(0 <= λ && λ <= 0.25)
-        error("Lambda must be between zero and 0.25")
-    end
+    !(0 <= λ <= 0.25) && throw(ArgumentError("Lambda must be between zero and 0.25"))
+    !(K > 0) && throw(ArgumentError("K must be greater than zero"))
+    !(niters > 0) && throw(ArgumentError("Number of iterations must be greater than zero"))
+
     @inline function g(norm∇I)
-        coef = (norm∇I / K)
-        denom = (T(1) .+ coef ⊙ coef)
+        coef = norm∇I / K
+        denom = T(1) .+ coef ⊙ coef
         return invert_color(denom)
     end
 
     output = deepcopy(image)
-    image = deepcopy(image)
-    for i in 1:niters
+    for _ in 1:niters
         pmad_kernel!(image, output, g, λ)
-        @swap!(image, output)
+        image, output = output, image
     end
     return output
 end
diff --git a/src/cloudmask.jl b/src/cloudmask.jl
index 9f5fd325..7a4d2abf 100644
--- a/src/cloudmask.jl
+++ b/src/cloudmask.jl
@@ -1,10 +1,58 @@
+function convert_to_255_matrix(img)::Matrix{Int}
+    img_clamped = clamp.(img, 0.0, 1.0)
+    return round.(Int, img_clamped * 255)
+end
+
+function _get_masks(
+    false_color_image::Matrix{RGB{Float64}};
+    prelim_threshold::Float64=Float64(110 / 255),
+    band_7_threshold::Float64=Float64(200 / 255),
+    band_2_threshold::Float64=Float64(190 / 255),
+    ratio_lower::Float64=0.0,
+    ratio_upper::Float64=0.75,
+    use_uint8::Bool=false,
+)::Tuple{BitMatrix,BitMatrix}
+
+    ref_view = channelview(false_color_image)
+    false_color_image_b7 = @view ref_view[1, :, :]
+    if use_uint8
+        false_color_image_b7 = convert_to_255_matrix(false_color_image_b7)
+    end
+
+    clouds_view = false_color_image_b7 .> prelim_threshold
+    mask_b7 = false_color_image_b7 .< band_7_threshold
+    mask_b2 = @view(ref_view[2, :, :])
+    if use_uint8
+        mask_b2 = convert_to_255_matrix(mask_b2)
+    end
+    mask_b2 = mask_b2 .> band_2_threshold
+
+    # First find all the pixels that meet threshold logic in band 7 (channel 1) and band 2 (channel 2)
+    # Masking clouds and discriminating cloud-ice
+
+    mask_b7b2 = mask_b7 .&& mask_b2
+
+    # Next find pixels that meet both thresholds and mask them from band 7 (channel 1) and band 2 (channel 2)
+    b7_masked = mask_b7b2 .* false_color_image_b7
+
+    _b2 = @view(ref_view[2, :, :])
+    b2_masked = use_uint8 ? convert_to_255_matrix(_b2) : _b2
+    b2_masked = mask_b7b2 .* b2_masked
+
+    cloud_ice = Float64.(b7_masked) ./ Float64.(b2_masked)
+    mask_cloud_ice = @. (cloud_ice >= ratio_lower) && (cloud_ice < ratio_upper)
+
+    return mask_cloud_ice, clouds_view
+end
+
+
 """
-    create_cloudmask(falsecolor_image; prelim_threshold, band_7_threshold, band_2_threshold, ratio_lower, ratio_upper)
+    create_cloudmask(false_color_image; prelim_threshold, band_7_threshold, band_2_threshold, ratio_lower, ratio_upper)
 
 Convert a 3-channel false color reflectance image to a 1-channel binary matrix; clouds = 0, else = 1. Default thresholds are defined in the published Ice Floe Tracker article: Remote Sensing of the Environment 234 (2019) 111406.
 
 # Arguments
-- `falsecolor_image`: corrected reflectance false color image - bands [7,2,1]
+- `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)
@@ -13,47 +61,41 @@ Convert a 3-channel false color reflectance image to a 1-channel binary matrix;
 
 """
 function create_cloudmask(
-    ref_image::Matrix{RGB{Float64}};
+    false_color_image::Matrix{RGB{Float64}};
     prelim_threshold::Float64=Float64(110 / 255),
     band_7_threshold::Float64=Float64(200 / 255),
     band_2_threshold::Float64=Float64(190 / 255),
     ratio_lower::Float64=0.0,
     ratio_upper::Float64=0.75,
 )::BitMatrix
-    # Setting thresholds
-    ref_view = channelview(ref_image)
-    ref_image_b7 = @view ref_view[1, :, :]
-    clouds_view = ref_image_b7 .> prelim_threshold
-    mask_b7 = ref_image_b7 .< band_7_threshold
-    mask_b2 = @view(ref_view[2, :, :]) .> band_2_threshold
-    # First find all the pixels that meet threshold logic in band 7 (channel 1) and band 2 (channel 2)
-    # Masking clouds and discriminating cloud-ice
+    mask_cloud_ice, clouds_view = _get_masks(
+        false_color_image,
+        prelim_threshold=prelim_threshold,
+        band_7_threshold=band_7_threshold,
+        band_2_threshold=band_2_threshold,
+        ratio_lower=ratio_lower,
+        ratio_upper=ratio_upper,
+    )
 
-    mask_b7b2 = mask_b7 .&& mask_b2
-    # Next find pixels that meet both thresholds and mask them from band 7 (channel 1) and band 2 (channel 2)
-    b7_masked = mask_b7b2 .* ref_image_b7
-    b2_masked = mask_b7b2 .* @view(ref_view[2, :, :])
-    cloud_ice = Float64.(b7_masked) ./ Float64.(b2_masked)
-    mask_cloud_ice = @. cloud_ice >= ratio_lower && cloud_ice < ratio_upper
     # Creating final cloudmask
     cloudmask = mask_cloud_ice .|| .!clouds_view
     return cloudmask
 end
 
 """
-    apply_cloudmask(ref_image, cloudmask)
+    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
-- `ref_image`: reference image, e.g. corrected reflectance false color image bands [7,2,1] or grayscale
+- `false_color_image`: reference image, e.g. corrected reflectance false color image bands [7,2,1] or grayscale
 - `cloudmask`: binary cloudmask with clouds = 0, else = 1
 
 """
 function apply_cloudmask(
-    ref_image::Matrix{RGB{Float64}}, cloudmask::AbstractArray{Bool}
+    false_color_image::Matrix{RGB{Float64}}, cloudmask::AbstractArray{Bool}
 )::Matrix{RGB{Float64}}
-    masked_image = cloudmask .* ref_image
+    masked_image = cloudmask .* false_color_image
     image_view = channelview(masked_image)
     cloudmasked_view = StackedView(
         zeroarray, @view(image_view[2, :, :]), @view(image_view[3, :, :])
@@ -63,13 +105,14 @@ function apply_cloudmask(
 end
 
 function apply_cloudmask(
-    ref_image::Matrix{Gray{Float64}}, cloudmask::AbstractArray{Bool}
+    false_color_image::Matrix{Gray{Float64}}, cloudmask::AbstractArray{Bool}
 )::Matrix{Gray{Float64}}
-    return Gray.(cloudmask .* ref_image)
+    return Gray.(cloudmask .* false_color_image)
 end
 
 function create_clouds_channel(
-    cloudmask::AbstractArray{Bool}, ref_image::Matrix{RGB{Float64}}
+    cloudmask::AbstractArray{Bool}, false_color_image::Matrix{RGB{Float64}}
 )::Matrix{Gray{Float64}}
-    return Gray.(@view(channelview(cloudmask .* ref_image)[1, :, :]))
+    return Gray.(@view(channelview(cloudmask .* false_color_image)[1, :, :]))
 end
+
diff --git a/src/histogram_equalization.jl b/src/histogram_equalization.jl
new file mode 100644
index 00000000..ae0b5740
--- /dev/null
+++ b/src/histogram_equalization.jl
@@ -0,0 +1,287 @@
+function to_uint8(img::AbstractMatrix{T}) where {T<:AbstractFloat}
+    img = UInt8.(round.(img))
+    img = clamp.(img, 0, 255)
+    return img
+end
+
+
+function anisotropic_diffusion_3D(I)
+    rgbchannels = get_rgb_channels(I)
+
+    for i in 1:3
+        rgbchannels[:, :, i] .= anisotropic_diffusion_2D(rgbchannels[:, :, i])
+    end
+
+    return rgbchannels
+
+end
+
+function anisotropic_diffusion_2D(I::AbstractMatrix{T}; gradient_threshold::Union{T,Nothing}=nothing, niter::Int=1) where {T}
+    if eltype(I) <: Int
+        I = Gray.(I ./ 255)
+    end
+
+    # Determine the gradient threshold if not provided
+    if gradient_threshold === nothing
+        dynamic_range = maximum(I) - minimum(I)
+        gradient_threshold = 0.1 * dynamic_range
+    end
+
+    # Padding the image (corrected)
+    padded_img = padarray(I, Pad(:replicate, (1, 1)))
+    dd = sqrt(2)
+    diffusion_rate = 1 / 8  # Fixed for maximal connectivity (8 neighbors)
+
+    for _ in 1:niter
+        # These are zero-indexed offset arrays
+        diff_img_north = padded_img[0:end-1, 1:end-1] .- padded_img[1:end, 1:end-1]
+        diff_img_east = padded_img[1:end-1, 1:end] .- padded_img[1:end-1, 0:end-1]
+        diff_img_nw = padded_img[0:end-2, 0:end-2] .- I
+        diff_img_ne = padded_img[0:end-2, 2:end] .- I
+        diff_img_sw = padded_img[2:end, 0:end-2] .- I
+        diff_img_se = padded_img[2:end, 2:end] .- I
+
+        # Exponential conduction coefficients
+        conduct_coeff_north = exp.(-(abs.(diff_img_north) ./ gradient_threshold) .^ 2)
+        conduct_coeff_east = exp.(-(abs.(diff_img_east) ./ gradient_threshold) .^ 2)
+        conduct_coeff_nw = exp.(-(abs.(diff_img_nw) ./ gradient_threshold) .^ 2)
+        conduct_coeff_ne = exp.(-(abs.(diff_img_ne) ./ gradient_threshold) .^ 2)
+        conduct_coeff_sw = exp.(-(abs.(diff_img_sw) ./ gradient_threshold) .^ 2)
+        conduct_coeff_se = exp.(-(abs.(diff_img_se) ./ gradient_threshold) .^ 2)
+
+
+        # Flux calculations
+        flux_north = conduct_coeff_north .* diff_img_north
+        flux_east = conduct_coeff_east .* diff_img_east
+        flux_nw = conduct_coeff_nw .* diff_img_nw
+        flux_ne = conduct_coeff_ne .* diff_img_ne
+        flux_sw = conduct_coeff_sw .* diff_img_sw
+        flux_se = conduct_coeff_se .* diff_img_se
+
+        # Back to regular 1-indexed arrays
+        flux_north_diff = flux_north[1:end-1, :] .- flux_north[2:end, :]
+        flux_east_diff = flux_east[:, 2:end] .- flux_east[:, 1:end-1]
+
+        # Discrete PDE solution
+        sum_ = (1 / (dd^2)) .* (flux_nw .+ flux_ne .+ flux_sw .+ flux_se)
+        I = I .+ diffusion_rate .* (flux_north_diff .- flux_north_diff .+ sum_)
+
+    end
+
+    return I
+end
+
+
+function imshow(img)
+    if typeof(img) <: BitMatrix
+        return Gray.(img)
+    end
+    Gray.(img ./ 255)
+end
+
+
+
+function adapthisteq(img::Matrix{T}, nbins=256, clip=0.01) where {T}
+    # Step 1: Normalize the image to [0, 1] based on its own min and max
+    image_min, image_max = minimum(img), maximum(img)
+    normalized_image = (img .- image_min) / (image_max - image_min)
+
+    # Step 2: Apply adaptive histogram equalization. equalize_adapthist handles the tiling to 1/8 of the image size (equivalent to 8x8 blocks in MATLAB)
+    equalized_image = sk_exposure.equalize_adapthist(
+        normalized_image,
+        clip_limit=clip,  # Equivalent to MATLAB's 'ClipLimit'
+        nbins=nbins         # Number of histogram bins. 255 is used to match the default in MATLAB script
+    )
+
+    # Step 3: Rescale the image back to the original range [image_min, image_max]
+    final_image = sk_exposure.rescale_intensity(equalized_image, in_range="image", out_range=(image_min, image_max))
+
+    # Convert back to the original data type if necessary
+    final_image = to_uint8(final_image)
+
+    return final_image
+end
+
+"""
+    get_rgb_channels(img)
+
+Get the RBC (Red, Blue, and Green) channels of an image.
+
+# Arguments
+- `img`: The input image.
+
+# Returns
+An m x n x 3 array the Red, Blue, and Green channels of the input image.
+
+"""
+function get_rgb_channels(img)
+    # TODO: might be able to use channelview instead
+    redc = red.(img) * 255
+    greenc = green.(img) * 255
+    bluec = blue.(img) * 255
+
+    return cat(redc, greenc, bluec, dims=3)
+end
+
+function get_tiles(false_color_image; rblocks, cblocks)
+    rtile, ctile = size(false_color_image)
+    tile_size = Tuple{Int,Int}((rtile / rblocks, ctile / cblocks))
+    return TileIterator(axes(false_color_image), tile_size)
+end
+
+function _process_image_tiles(true_color_image, false_color_image, landmask, tiles, white_threshold, entropy_threshold, white_fraction_threshold, prelim_threshold, band_7_threshold, band_2_threshold)
+    # Get the red channel for cloud detection
+    clouds_red = _get_red_channel_cloud_cae(
+        false_color_image=false_color_image, landmask=landmask,
+        prelim_threshold=prelim_threshold,
+        band_7_threshold=band_7_threshold,
+        band_2_threshold=band_2_threshold,
+    )
+
+    # Apply diffuse (anisotropic diffusion) to each channel of true color image
+    true_color_diffused = IceFloeTracker.diffusion(float64.(true_color_image), 0.1, 75, 3)
+
+    rgbchannels = get_rgb_channels(true_color_diffused)
+
+    # For each tile, compute the entropy in the false color tile, and the fraction of white and black pixels
+    for tile in tiles
+        clouds_tile = clouds_red[tile...]
+        entropy = Images.entropy(clouds_tile)
+        whitefraction = sum(clouds_tile .> white_threshold) / length(clouds_tile)
+
+        # If the entropy is above a threshold, and the fraction of white pixels is above a threshold, then apply histogram equalization to the tiles of each channel of the true color image. Otherwise, keep the original tiles.
+        if entropy > entropy_threshold && whitefraction > white_fraction_threshold
+            for i in 1:3
+                eqhist = adapthisteq(rgbchannels[:, :, i][tile...])
+                @view(rgbchannels[:, :, i])[tile...] .= eqhist
+            end
+        end
+    end
+
+    return rgbchannels
+end
+
+
+"""
+    conditional_histeq(true_color_image, clouds, landmask, rblocks::Int=8, cblocks::Int=6, entropy_threshold::AbstractFloat=4.0, white_threshold::AbstractFloat=25.5, white_fraction_threshold::AbstractFloat=0.4)
+
+Performs conditional histogram equalization on a true color image.
+
+# Arguments
+- `true_color_image`: The true color image to be equalized.
+- `false_color_image`: The false color image used to determine the regions to equalize.
+- `landmask`: The land mask indicating the land regions in the image.
+- `rblocks`: The number of row-blocks to divide the image into for histogram equalization. Default is 8.
+- `cblocks`: The number of column-blocks to divide the image into for histogram equalization. Default is 6.
+- `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.
+
+"""
+function conditional_histeq(
+    true_color_image,
+    false_color_image,
+    landmask,
+    rblocks::Int,
+    cblocks::Int,
+    entropy_threshold::AbstractFloat=4.0,
+    white_threshold::AbstractFloat=25.5,
+    white_fraction_threshold::AbstractFloat=0.4,
+    prelim_threshold=110.0,
+    band_7_threshold=200.0,
+    band_2_threshold=190.0
+)
+
+    tiles = get_tiles(false_color_image, rblocks=rblocks, cblocks=cblocks)
+
+    rgbchannels_equalized = _process_image_tiles(
+        true_color_image,
+        false_color_image,
+        landmask,
+        tiles,
+        white_threshold,
+        entropy_threshold,
+        white_fraction_threshold,
+        prelim_threshold,
+        band_7_threshold,
+        band_2_threshold
+    )
+
+    return rgbchannels_equalized
+
+end
+
+"""
+    conditional_histeq(
+    true_color_image,
+    false_color_image,
+    landmask,
+    side_length::Int,
+    entropy_threshold::AbstractFloat=4.0,
+    white_threshold::AbstractFloat=25.5,
+    white_fraction_threshold::AbstractFloat=0.4,
+    prelim_threshold=110.0,
+    band_7_threshold=200.0,
+    band_2_threshold=190.0
+)
+
+Performs conditional histogram equalization on a true color image using tiles of approximately sidelength size `side_length`. If a perfect tiling is not possible, the tiling on the egde of the image is adjusted to ensure that the tiles are as close to `side_length` as possible. See `get_tiles(array, side_length)` for more details.
+"""
+function conditional_histeq(
+    true_color_image,
+    false_color_image,
+    landmask,
+    side_length::Int,
+    entropy_threshold::AbstractFloat=4.0,
+    white_threshold::AbstractFloat=25.5,
+    white_fraction_threshold::AbstractFloat=0.4,
+    prelim_threshold=110.0,
+    band_7_threshold=200.0,
+    band_2_threshold=190.0
+)
+
+    side_length = IceFloeTracker.get_optimal_tile_size(side_length, size(false_color_image))
+
+    tiles = IceFloeTracker.get_tiles(false_color_image, side_length)
+
+    rgbchannels_equalized = _process_image_tiles(
+        true_color_image,
+        false_color_image,
+        landmask,
+        tiles,
+        white_threshold,
+        entropy_threshold,
+        white_fraction_threshold,
+        prelim_threshold,
+        band_7_threshold,
+        band_2_threshold
+    )
+
+    return rgbchannels_equalized
+
+end
+
+function _get_red_channel_cloud_cae(; false_color_image, landmask,
+    prelim_threshold=110.0,
+    band_7_threshold=200.0,
+    band_2_threshold=190.0,
+)
+    mask_cloud_ice, clouds_view = IceFloeTracker._get_masks(
+        false_color_image,
+        prelim_threshold=prelim_threshold,
+        band_7_threshold=band_7_threshold,
+        band_2_threshold=band_2_threshold,
+        use_uint8=true,
+    )
+
+    clouds_view[mask_cloud_ice] .= 0
+
+    # remove clouds and land from red channel
+    red_channel = Int.(red.(false_color_image) * 255)
+    red_channel[clouds_view.|landmask] .= 0
+
+    return red_channel
+end
diff --git a/src/normalization.jl b/src/normalization.jl
index 256b6d39..e16619e7 100644
--- a/src/normalization.jl
+++ b/src/normalization.jl
@@ -74,7 +74,7 @@ Sharpen `truecolor_image`.
 - `nbins`: number of bins during histogram equalization
 - `rblocks`: number of row blocks to divide input image during equalization
 - `cblocks`: number of column blocks to divide input image during equalization
-- `clip`: Thresholds for clipping histogram bins (0–1); values closer to one minimize contrast enhancement, values closer to zero maximize contrast enhancement 
+- `clip`: Thresholds for clipping histogram bins (0–1); values closer to one minimize contrast enhancement, values closer to zero maximize contrast enhancement
 - `smoothing_param`: pixel radius for gaussian blurring (1–10)
 - `intensity`: amount of sharpening to perform
 """
@@ -85,14 +85,16 @@ function imsharpen(
     kappa::Real=75,
     niters::Int64=3,
     nbins::Int64=255,
-    rblocks::Int64=10,
-    cblocks::Int64=10,
-    clip::Float64=0.86,
+    rblocks::Int64=10, # matlab default is 8 CP
+    cblocks::Int64=10, # matlab default is 8 CP
+    clip::Float64=0.86, # matlab default is 0.01 CP
     smoothing_param::Int64=10,
     intensity::Float64=2.0,
 )::Matrix{Float64}
     input_image = IceFloeTracker.apply_landmask(truecolor_image, landmask_no_dilate)
+
     input_image .= IceFloeTracker.diffusion(input_image, lambda, kappa, niters)
+
     masked_view = Float64.(channelview(input_image))
 
     eq = [
@@ -115,7 +117,7 @@ end
     imsharpen_gray(imgsharpened, landmask)
 
 Apply landmask and return Gray type image in colorview for normalization.
-    
+
 """
 function imsharpen_gray(
     imgsharpened::Matrix{Float64}, landmask::AbstractArray{Bool}
@@ -123,3 +125,5 @@ function imsharpen_gray(
     image_sharpened_landmasked = apply_landmask(imgsharpened, landmask)
     return colorview(Gray, image_sharpened_landmasked)
 end
+
+
diff --git a/src/tilingutils.jl b/src/tilingutils.jl
new file mode 100644
index 00000000..babcd993
--- /dev/null
+++ b/src/tilingutils.jl
@@ -0,0 +1,112 @@
+function getfit(dims::Tuple{Int,Int}, side_length::Int)::Tuple{Int,Int}
+    return dims .÷ side_length
+end
+
+function get_area_missed(side_length::Int, dims::Tuple{Int,Int}, area::Int)::Float64
+    return 1 - prod(getfit(dims, side_length)) * side_length^2 / area
+end
+
+"""
+    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 during tiling. In case of a tie, it prefers the larger tile size.
+
+# Example
+```
+julia> get_optimal_tile_size(3, (10, 7))
+2
+```
+"""
+function get_optimal_tile_size(l0::Int, dims::Tuple{Int,Int})::Int
+    l0 < 2 && error("l0 must be at least 2")
+    any(l0 .> dims) && error("l0 = $l0 is too large for the given dimensions $dims")
+
+    area = prod(dims)
+    minimal_shift = l0 == 2 ? 0 : 1
+    candidates = [l0 + i for i in -minimal_shift:1]
+
+    minl, M = 0, Inf
+    for side_length in candidates
+        missedarea = get_area_missed(side_length, dims, area)
+        if missedarea <= M # prefer larger side_length in case of tie
+            M, minl = missedarea, side_length
+        end
+    end
+    return minl
+end
+
+"""
+    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`.
+"""
+function get_tile_meta(tile)
+    a, c = first.(tile)
+    b, d = last.(tile)
+    return [a, b, c, d]
+end
+
+
+function bump_tile(tile, dims)
+    extrarows, extracols = dims
+    a, b, c, d = get_tile_meta(tile)
+    b += extrarows
+    d += extracols
+    return (a:b, c:d)
+end
+
+function get_tile_dims(tile)
+    a, b, c, d = get_tile_meta(tile)
+    width, height = d - c + 1, b - a + 1
+    return (width, height)
+end
+
+
+"""
+    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`.
+"""
+function get_tiles(array, side_length)
+
+    tiles = TileIterator(axes(array), (side_length, side_length)) |> collect
+
+    bottombump, rightbump = mod.(size(array), side_length)
+
+    if bottombump == 0 && rightbump == 0
+        return tiles
+    end
+
+    crop_height, crop_width = 0, 0
+
+    # Adjust bottom edge if necessary
+    if bottombump <= side_length ÷ 2
+        bottom_edge = tiles[end-1, :]
+        tiles[end-1, :] .= bump_tile.(bottom_edge, Ref((bottombump, 0)))
+        crop_height += 1
+    end
+
+    # Adjust right edge if necessary
+    if rightbump <= side_length ÷ 2
+        right_edge = tiles[:, end-1]
+        tiles[:, end-1] .= bump_tile.(right_edge, Ref((0, rightbump)))
+        crop_width += 1
+    end
+
+    return tiles[1:end-crop_height, 1:end-crop_width]
+end
diff --git a/test/Project.toml b/test/Project.toml
index 2baeb7e4..53499230 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,3 +10,5 @@ Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
+TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
diff --git a/test/runtests.jl b/test/runtests.jl
index aca5652b..5a77ed59 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,50 +1,52 @@
+using ArgParse: @add_arg_table!, ArgParseSettings, add_arg_group!, parse_args
+using DataFrames
+using Dates
+using DelimitedFiles
 using IceFloeTracker
+using ImageTransformations: imrotate
 using Images
-using Test
-using DelimitedFiles
-using Dates
-using DataFrames
 using Random
-using ImageTransformations: imrotate
-using ArgParse: ArgParseSettings, @add_arg_table!, add_arg_group!, parse_args
+using Test
+using TestImages
+using TiledIteration
 include("test_error_rate.jl")
 include("config.jl")
 
 # Setting things up (see config.jl)
 
-## Get all test files filenames "test-*" in test folder and their corresponding names/label 
+## Get all test files filenames "test-*" in test folder and their corresponding names/label
 alltests = [f for f in readdir() if startswith(f, "test-")]
-testnames = [n[6:(end - 3)] for n in alltests]
+testnames = [n[6:(end-3)] for n in alltests]
 
 ## Put the filenames to test below
 
-to_test = alltests # uncomment this line to run all tests or add individual files below 
+to_test = alltests # uncomment this line to run all tests or add individual files below
 [
-    # "test-create-landmask.jl",
-    # "test-create-cloudmask.jl",
-    # "test-normalize-image.jl",
-    # "test-persist.jl",
-    # "test-utils-padding.jl",
-    # "test-discrim-ice-water.jl",
-    # "test-find-ice-labels.jl",
-    # "test-segmentation-a.jl",
-    # "test-segmentation-b.jl",
-    # "test-segmentation-watershed.jl",
-    # "test-segmentation-f.jl",
-    # "test-bwtraceboundary.jl",
-    # "test-resample-boundary.jl",
-    # "test-regionprops.jl",
-    # "test-psi-s.jl",
-    # "test-crosscorr.jl"
-    # "test-bwperim.jl",
-    # "test-bwareamaxfilt.jl"
-    # "test-register-mismatch.jl",
-    # "test-utils-imextendedmin.jl",
-    # "test-morphSE.jl",
-    # "test-hbreak.jl",
-    # "test-bridge.jl",
-    # "test-branch.jl"
-    # "test-pipeline.jl"
+# "test-create-landmask.jl",
+# "test-create-cloudmask.jl",
+# "test-normalize-image.jl",
+# "test-persist.jl",
+# "test-utils-padding.jl",
+# "test-discrim-ice-water.jl",
+# "test-find-ice-labels.jl",
+# "test-segmentation-a.jl",
+# "test-segmentation-b.jl",
+# "test-segmentation-watershed.jl",
+# "test-segmentation-f.jl",
+# "test-bwtraceboundary.jl",
+# "test-resample-boundary.jl",
+# "test-regionprops.jl",
+# "test-psi-s.jl",
+# "test-crosscorr.jl"
+# "test-bwperim.jl",
+# "test-bwareamaxfilt.jl"
+# "test-register-mismatch.jl",
+# "test-utils-imextendedmin.jl",
+# "test-morphSE.jl",
+# "test-hbreak.jl",
+# "test-bridge.jl",
+# "test-branch.jl"
+# "test-pipeline.jl"
 ]
 
 # Run the tests
diff --git a/test/test-conditional-adaptive-histeq.jl b/test/test-conditional-adaptive-histeq.jl
new file mode 100644
index 00000000..48fa8871
--- /dev/null
+++ b/test/test-conditional-adaptive-histeq.jl
@@ -0,0 +1,65 @@
+begin
+    datadir = joinpath(@__DIR__, "test_inputs/")
+    path_true_color_image = joinpath(datadir, "NE_Greenland_truecolor.2020162.aqua.250m.tiff")
+    path_false_color_image = joinpath(datadir, "NE_Greenland_reflectance.2020162.aqua.250m.tiff")
+    true_color_image = float64.(load(path_true_color_image))
+    false_color_image = float64.(load(path_false_color_image))
+    dilated_landmask = BitMatrix(load(joinpath(datadir, "matlab_landmask.png")))
+end
+
+function test_cloud_image_workflow()
+    @testset "Prereq cloud image" begin
+        redchannel_cahe = IceFloeTracker._get_red_channel_cloud_cae(
+            false_color_image=false_color_image,
+            landmask=dilated_landmask,
+            prelim_threshold=110.0,
+            band_7_threshold=200.0,
+            band_2_threshold=190.0,
+        )
+
+        @test sum(redchannel_cahe) == 1_320_925_065
+    end
+end
+
+
+function test_adaphisteq()
+    @testset "Adaptive histogram equalization" begin
+        img = IceFloeTracker.convert_to_255_matrix(testimage("cameraman"))
+        img_eq = IceFloeTracker.adapthisteq(img)
+        @test sum(img_eq) == 32_387_397
+    end
+end
+
+function test_conditional_adaptivehisteq()
+    @testset "Conditional adaptivehisteq" begin
+
+        # Using rblocks = 8, cblocks = 6
+        true_color_eq = IceFloeTracker.conditional_histeq(
+            true_color_image,
+            false_color_image,
+            dilated_landmask,
+            8,
+            6)
+
+        # This differs from MATLAB script due to disparity in the implementations
+        # of the adaptive histogram equalization / diffusion functions
+        # For the moment testing for regression
+        @test sum(IceFloeTracker.to_uint8(true_color_eq[:, :, 1])) == 6_372_159_606
+
+
+        # Use custom tile size
+        side_length = size(true_color_eq, 1) ÷ 8
+        true_color_eq = IceFloeTracker.conditional_histeq(
+            true_color_image,
+            false_color_image,
+            dilated_landmask,
+            side_length)
+        @test sum(IceFloeTracker.to_uint8(true_color_eq[:, :, 1])) == 6_328_796_398
+    end
+end
+
+@testset "Conditional adaptivehisteq" begin
+    test_cloud_image_workflow()
+    test_adaphisteq()
+    test_conditional_adaptivehisteq()
+end
diff --git a/test/test-create-cloudmask.jl b/test/test-create-cloudmask.jl
index be066606..a44595f2 100644
--- a/test/test-create-cloudmask.jl
+++ b/test/test-create-cloudmask.jl
@@ -12,12 +12,12 @@
     @time masked_image = IceFloeTracker.apply_cloudmask(ref_image, cloudmask)
 
     # test for percent difference in cloudmask images
-    @test (@test_approx_eq_sigma_eps masked_image matlab_cloudmask [0, 0] 0.005) == nothing
+    @test (@test_approx_eq_sigma_eps masked_image matlab_cloudmask [0, 0] 0.005) === nothing
 
     # test for create_clouds_channel
     clouds_channel_expected = load(clouds_channel_test_file)
     clds_channel = IceFloeTracker.create_clouds_channel(cloudmask, ref_image)
-    @test (@test_approx_eq_sigma_eps (clds_channel) (clouds_channel_expected) [0, 0] 0.005) ==
+    @test (@test_approx_eq_sigma_eps (clds_channel) (clouds_channel_expected) [0, 0] 0.005) ===
         nothing
 
     # Persist output images
diff --git a/test/test-discrim-ice-water.jl b/test/test-discrim-ice-water.jl
index fb1c065f..a4c2f981 100644
--- a/test/test-discrim-ice-water.jl
+++ b/test/test-discrim-ice-water.jl
@@ -17,7 +17,7 @@
     ice_water_discrim = IceFloeTracker.discriminate_ice_water(
         falsecolor_image, normalized_image, landmask, cloudmask
     )
-    @test (@test_approx_eq_sigma_eps ice_water_discrim matlab_ice_water_discrim [0, 0] 0.065) ==
+    @test (@test_approx_eq_sigma_eps ice_water_discrim matlab_ice_water_discrim [0, 0] 0.065) ===
         nothing
 
     # persist generated image
diff --git a/test/test-normalize-image.jl b/test/test-normalize-image.jl
index 909b512e..d0e0d84c 100644
--- a/test/test-normalize-image.jl
+++ b/test/test-normalize-image.jl
@@ -20,12 +20,12 @@
 
     @time image_diffused = IceFloeTracker.diffusion(input_landmasked, 0.1, 75, 3)
 
-    @test (@test_approx_eq_sigma_eps image_diffused matlab_diffused [0, 0] 0.0054) ==
+    @test (@test_approx_eq_sigma_eps image_diffused matlab_diffused [0, 0] 0.0054) ===
         nothing
 
-    @test (@test_approx_eq_sigma_eps input_landmasked image_diffused [0, 0] 0.004) ==
+    @test (@test_approx_eq_sigma_eps input_landmasked image_diffused [0, 0] 0.004) ===
         nothing
-    @test (@test_approx_eq_sigma_eps input_landmasked matlab_diffused [0, 0] 0.007) ==
+    @test (@test_approx_eq_sigma_eps input_landmasked matlab_diffused [0, 0] 0.007) ===
         nothing
 
     diffused_image_filename =
@@ -43,7 +43,7 @@
         i in 1:3
     ]
     image_equalized = colorview(RGB, eq...)
-    @test (@test_approx_eq_sigma_eps image_equalized matlab_equalized [0, 0] 0.051) ==
+    @test (@test_approx_eq_sigma_eps image_equalized matlab_equalized [0, 0] 0.051) ===
         nothing
 
     equalized_image_filename =
@@ -59,7 +59,7 @@
     @time image_sharpened_gray = IceFloeTracker.imsharpen_gray(
         sharpenedimg, landmask_bitmatrix
     )
-    @test (@test_approx_eq_sigma_eps image_sharpened_gray matlab_sharpened [0, 0] 0.046) ==
+    @test (@test_approx_eq_sigma_eps image_sharpened_gray matlab_sharpened [0, 0] 0.046) ===
         nothing
 
     sharpened_image_filename =
@@ -76,7 +76,7 @@
     )
 
     #test for percent difference in normalized images
-    @test (@test_approx_eq_sigma_eps normalized_image matlab_norm_image [0, 0] 0.045) ==
+    @test (@test_approx_eq_sigma_eps normalized_image matlab_norm_image [0, 0] 0.045) ===
         nothing
 
     normalized_image_filename =
diff --git a/test/test-segmentation-b.jl b/test/test-segmentation-b.jl
index bc593b3b..04356c9e 100644
--- a/test/test-segmentation-b.jl
+++ b/test/test-segmentation-b.jl
@@ -21,7 +21,7 @@
     IceFloeTracker.@persist matlab_ice_intersect "./test_outputs/matlab_ice_intersect.png" true
 
     @test typeof(segB.not_ice) == typeof(matlab_not_ice_mask)
-    @test (@test_approx_eq_sigma_eps segB.not_ice matlab_not_ice_mask [0, 0] 0.001) ==
+    @test (@test_approx_eq_sigma_eps segB.not_ice matlab_not_ice_mask [0, 0] 0.001) ===
         nothing
 
     @test typeof(segB.not_ice_bit) == typeof(matlab_not_ice_mask .> 0.499)
diff --git a/test/test-tilingutils.jl b/test/test-tilingutils.jl
new file mode 100644
index 00000000..40c152dc
--- /dev/null
+++ b/test/test-tilingutils.jl
@@ -0,0 +1,74 @@
+using IceFloeTracker: get_optimal_tile_size, get_tile_meta, bump_tile, get_tiles
+gots = get_optimal_tile_size
+
+@testset "Tiling utils" begin
+    @testset "get_optimal_tile_size" begin
+        @test gots(2, (10, 10)) == 2 # disregard tiles of 1 pixel
+        @test gots(3, (15, 15)) == 3
+        @test gots(4, (20, 20)) == 5 # prefer larger tile size
+
+        # Test with edge case for minimum l0
+        @test gots(3, (5, 5)) == 4
+        @test gots(5, (5, 5)) == 5
+
+        # Test with non-square dimensions
+        @test gots(3, (10, 20)) == 2
+        @test gots(3, (10, 7)) == 2
+
+        # Test error handling for invalid l0
+        @test_throws ErrorException gots(1, (10, 10))
+        @test_throws ErrorException gots(7, (5, 5)) == 5
+    end
+
+    tile = (1:2, 3:4)
+
+    @testset "get_tile_meta" begin
+        @test get_tile_meta(tile) == [1, 2, 3, 4]
+    end
+
+    @testset "bump_tile" begin
+        extrarows, extracols = rand(1:100, 2)
+        bumpby = (extrarows, extracols)
+        @test bump_tile(tile, bumpby) == (1:2+extrarows, 3:4+extracols)
+    end
+
+    @testset "get_tiles" begin
+        # unadjusted tiles
+        _get_tiles(array, side_length) = TileIterator(axes(array), (side_length, side_length)) |> collect
+
+        array = rand(40, 20)
+
+        # Total coverage, no adjustment needed
+        side_length = 5
+        tiles = _get_tiles(array, side_length)
+        @test tiles == get_tiles(array, side_length)
+
+        # Leftovers greater than side_length ÷ 2, no adjustment needed
+        side_length = 7 # bumpby = (5, 6)
+        tiles = _get_tiles(array, side_length)
+        adjusted_tiles = get_tiles(array, side_length)
+        @test tiles == get_tiles(array, side_length)
+
+        # adjust right edge tiles
+        side_length = 6 # bumpby = (4, 2) => crop 1 column and bump by 2
+        tiles = _get_tiles(array, side_length)
+        adjusted_tiles = get_tiles(array, side_length)
+        @test all(size(tiles) .- size(adjusted_tiles) .== (0, 1))
+
+        # general case with both edges adjusted
+        expected_tiles = [(1:5, 1:5) (1:5, 6:10) (1:5, 11:16);
+            (6:12, 1:5) (6:12, 6:10) (6:12, 11:16)]
+
+        array = rand(12, 16)
+        side_length = 5
+        newtiles = get_tiles(array, side_length)
+        @test expected_tiles == newtiles
+
+        lowerleft_tile = newtiles[end, 1]
+        _, b, _, _ = get_tile_meta(lowerleft_tile)
+
+        lowerright_tile = newtiles[end, end]
+        _, _, _, d = get_tile_meta(lowerright_tile)
+        @test (b, d) == size(array)
+    end
+end