Tweaks to accommodate new imaging setup

Project.toml

@@ -10,11 +10,15 @@ Gtk4 = "9db2cae5-386f-4011-9d63-a5602296539b"
 GtkObservables = "8710efd8-4ad6-11eb-33ea-2d5ceb25a41c"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
+ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264"
 ImageSegmentation = "80713f31-8817-5129-9cf8-209ff8fb23e1"
 ImageView = "86fae568-95e7-573e-a6b2-d8a6b900c9ef"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OpenCV = "f878e3a2-a245-4720-8660-60795d644f2a"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 XLSX = "fdbf4ff8-1666-58a4-91e7-1b58723a45e0"
 
 [compat]
@@ -24,10 +28,12 @@ Gtk4 = "0.7"
 GtkObservables = "2"
 ImageCore = "0.10"
 ImageIO = "0.6"
+ImageMagick = "1.4.2"
 ImageMorphology = "0.4"
 ImageSegmentation = "1.9"
 ImageView = "0.12"
 JLD2 = "0.5"
+OpenCV = "4.6.1"
 XLSX = "0.10"
 julia = "1.10"
 

docs/src/index.md

@@ -13,7 +13,7 @@ The yellow spot corresponds to a stimulus provided by the experimenter, and the
 
 Tips on image quality:
 
-- Put the stimulus near one of the four corners
+- Put the stimulus near one of the four corners, keeping its location as consistent as possible between images
 - Ensure lighting is fairly uniform
 - Make sure that any extraneous marks (e.g., the black writing in the image above) are of a very different color from scent marks.
 - Ensure that all your images are of the same size (i.e., same number of pixels horizontally and vertically), even if there are some extra pixels on the edges of the image
@@ -99,12 +99,37 @@ Alternatively, you can supply a list of files:
 julia> gui("results_file_name", ["PictureA.png", "mouse7.png"])
 ```
 
+Additionally, you can supply a calibration image to improve segmentation by correcting for uneven illumination.
+
+![An image to be used for calibration](assets/blurred_calibration.png)
+
+```
+julia> gui("results_file_name", ["PictureA.png", "mouse7.png"]; background_path="calibration_image.png")
+```
+
+There are a few more keyword arguments that might be useful:
+  - `crop_top`, `crop_bottom`, `crop_left`, and `crop_right` specify a number of pixels to be cropped from images along each side
+  - `expectedloc` specifies the expected location of the stimulus spot in pixels (after cropping)
+
+```
+julia> gui("results_file_name", glob"Picture*.png"; 
+            background_path="calibration_image.png,
+            crop_top=93,
+            crop_bottom=107,
+            crop_left=55,
+            crop_right=45,    
+            expectedloc=[1600,3333]
+        )
+```
+
 However you launch it, you should see something like this:
 
 ![GUI](assets/gui.png)
 
 On the top is the raw image. On the bottom is the segmented image; you should visually compare the two to check whether you're pleased with the quality of the segmentation.
-(If not, click "Skip" to omit that file from analysis.)
+
+If the default segmentation doesn't look quite right, try adjusting the Color Similarity Threshold value using the buttons at the bottom of the GUI.
+(If you can't obtain a segmentation that you're happy with, click "Skip" to omit that file from analysis.)
 
 If you like the segmentation, your tasks are:
 - click on all the checkboxes with colors that correspond to urine spots. You'll notice that the stimulus spot is pre-clicked (you can correct its choice if it didn't pick correctly). Most of the time there will be only one you need to check, but you can click more than one.
@@ -124,6 +149,8 @@ julia> using CounterMarking, ImageView   # load packages (if this is a fresh ses
 
 julia> count = density_map("results_file_name.jld2");
 
+julia> dmap = count ./ maximum(count);
+
 julia> dct = imshow(dmap);
 ```
 

src/CounterMarking.jl

@@ -11,6 +11,8 @@ using Gtk4
 using GtkObservables
 using ImageView
 using Random
+using LinearAlgebra
+using Statistics
 
 export segment_image, stimulus_index, spots, Spot, upperleft
 export writexlsx, process_images, density_map

src/gui.jl

@@ -20,15 +20,31 @@ The JLD2 file can be used by [`density_map`](@ref).
 """
 function gui(
         outbase::AbstractString, files;
-        colors=distinguishable_colors(15, [RGB(1, 1, 1)]; dropseed=true),
+        colors=distinguishable_colors(15, [RGB(1, 1, 1)]; dropseed=false),
         btnclick = Condition(),         # used for testing
         whichbutton = Ref{Symbol}(),    # used for testing
         preclick::Union{Int,Nothing} = nothing,  # used for testing
+        background_path = nothing,  # used to correct for non-uniform illumination (see joinpath(pkgdir(@__MODULE__),"docs","src","assets","blurred_calibration.png"))
+        crop_top::Int = 0,          # crop this many pixels off of each side
+        crop_bottom::Int = 0,
+        crop_left::Int = 0,
+        crop_right::Int = 0,
+        colorproj = RGB{Float32}(1, 1, -2), # used for identifying the stimulus
+        expectedloc = nothing,              # ""
     )
     channelpct(x) = string(round(Int, x * 100)) * '%'
 
     outbase, _ = splitext(outbase)
 
+    # Initialize segmented image and color similarity threshold
+    img = nothing
+    rescaledimg = nothing
+    bkgimg = isnothing(background_path) ? nothing : Float32.(Gray.(load(background_path)[crop_top+1:end-crop_bottom, crop_left+1:end-crop_right]))
+    bkgmean = isnothing(bkgimg) ? nothing : Float32(mean(bkgimg))
+    seg = nothing
+    threshold = 0.15
+    labels2idx = Dict{Int,Int}()
+    idx2labels = Dict{Int,Int}()
     # Set up basic properties of the window
     winsize = round.(Int, 0.8 .* screen_size())
     win = GtkWindow("CounterMarking", winsize...)
@@ -86,6 +102,64 @@ function gui(
         seggrid[col, row] = cb.widget
         push!(cbs, cb)
     end
+
+    thrshbx = GtkBox(:h)
+    thrshlbl = GtkLabel("Color Similarity Threshold: $(threshold)")
+    push!(thrshbx, thrshlbl)
+    thrshbtnbx = GtkBox(:v)
+    push!(thrshbx, thrshbtnbx)
+    lgincbtn = GtkButton("+0.05")
+    smincbtn = GtkButton("+0.01")
+    smdecbtn = GtkButton("-0.01")
+    lgdecbtn = GtkButton("-0.05")
+    push!(thrshbtnbx, lgincbtn)
+    push!(thrshbtnbx, smincbtn)
+    push!(thrshbtnbx, smdecbtn)
+    push!(thrshbtnbx, lgdecbtn)
+    push!(guibx, thrshbx)
+
+    update_threshold = change -> begin
+        newvalue = round(threshold + change; digits=2)
+        try
+            seg = segment_image(rescaledimg; threshold = newvalue)
+            nsegs = length(segment_labels(seg))
+            nsegs > length(colors) && @warn "More than $(length(colors)) segments ($(nsegs)). Excluded ones will be displayed in white and will not be selectable"
+            empty!(labels2idx)
+            empty!(idx2labels)
+            for (i,l) in enumerate(sort(seg.segment_labels; rev=true, by=(l -> segment_pixel_count(seg,l))))
+                idx = i == 1 ? 2 : i == 2 ? 1 : i
+                idx = i > 15 ? 1 : idx  # show remaining segments in white
+                push!(labels2idx, l=>idx)
+                push!(idx2labels, idx=>l)
+            end
+            centroidsacc, _ = get_centroidsacc(seg.image_indexmap)
+            istim = labels2idx[stimulus_index(seg, centroidsacc; colorproj = colorproj, expectedloc = expectedloc)]            
+            for (j, cb) in enumerate(cbs)
+                # set_gtk_property!(cb, "active", j <= nsegs)
+                cb[] = (j == istim || j == preclick)
+            end
+            imshow(canvases[1, 1], img)
+            imshow(canvases[2, 1], map(i->colors[labels2idx[i]], labels_map(seg)))
+            threshold = newvalue
+            thrshlbl.label = "Color Similarity Threshold: $threshold"
+        catch e
+            @show e
+            # display?
+        end
+    end
+    signal_connect(lgincbtn, "clicked") do w, others...
+        update_threshold(0.05)
+    end
+    signal_connect(smincbtn, "clicked") do w, others...
+        update_threshold(0.01)
+    end
+    signal_connect(smdecbtn, "clicked") do w, others...
+        update_threshold(-0.01)
+    end
+    signal_connect(lgdecbtn, "clicked") do w, others...
+        update_threshold(-0.05)
+    end
+
     # Add "Done & Next" and "Skip" buttons
     donebtn = button("Done & Next")
     skipbtn = button("Skip")
@@ -102,31 +176,44 @@ function gui(
 
     results = []
     for (i, file) in enumerate(files)
-        img = color.(load(file))
-        seg = segment_image(img)
+        threshold = 0.15
+        thrshlbl.label = "Color Similarity Threshold: $threshold"
+        img = color.(load(file))[crop_top+1:end-crop_bottom, crop_left+1:end-crop_right]
+        rescaledimg = isnothing(bkgimg) ? img : (img ./ bkgimg .* bkgmean)
+
+        seg = segment_image(rescaledimg; threshold = threshold)
         nsegs = length(segment_labels(seg))
-        @assert nsegs < length(colors) "Too many segments for colors"
-        istim = stimulus_index(seg)
+        nsegs > length(colors) && @warn "More than $(length(colors)) segments ($(nsegs)). Excluded ones will be displayed in white and will not be selectable"
+        empty!(labels2idx)
+        empty!(idx2labels)
+        for (i,l) in enumerate(sort(seg.segment_labels; rev=true, by=(l -> segment_pixel_count(seg,l))))
+            idx = i == 1 ? 2 : i == 2 ? 1 : i
+            idx = i > 15 ? 1 : idx  # show remaining segments in white
+            push!(labels2idx, l=>idx)
+            push!(idx2labels, idx=>l)
+        end
+        centroidsacc, _ = get_centroidsacc(seg.image_indexmap)
+        istim = labels2idx[stimulus_index(seg, centroidsacc; colorproj = colorproj, expectedloc = expectedloc)]            
         for (j, cb) in enumerate(cbs)
             # set_gtk_property!(cb, "active", j <= nsegs)
             cb[] = (j == istim || j == preclick)
-        end
+        end    
         imshow(canvases[1, 1], img)
-        imshow(canvases[2, 1], map(i->colors[i], labels_map(seg)))
+        imshow(canvases[2, 1], map(l->colors[labels2idx[l]], labels_map(seg)))
 
         wait(btnclick)
         whichbutton[] == :skip && continue
 
         keep = Int[]
         for (j, cb) in enumerate(cbs)
             if cb[]
-                push!(keep, j)
+                push!(keep, idx2labels[j])
             end
         end
         pixelskeep = map(i -> i ∈ keep, labels_map(seg))
         L = label_components(pixelskeep)
         newseg = SegmentedImage(img, L)
-        spotdict, stimulus = spots(newseg)
+        spotdict, stimulus = spots(newseg; colorproj = colorproj, expectedloc = expectedloc)
         push!(results, (file, spotdict, stimulus, newseg))
     end
 

src/segment.jl

@@ -8,9 +8,9 @@ Larger `threshold` values will result in fewer segments.
 """
 function segment_image(
         img::AbstractMatrix{<:Color};
-        threshold::Real = 0.2,      # threshold for color similarity in region growing
+        threshold::Real = 0.15,     # threshold for color similarity in region growing
         prune::Bool = true,         # prune small segments
-        min_size::Int = 50,         # minimum size of segments to keep
+        min_size::Int = 500,        # minimum size of segments to keep
     )
     seg = unseeded_region_growing(img, threshold)
     if prune
@@ -22,16 +22,30 @@ end
 segment_image(img::AbstractMatrix{<:Colorant}; kwargs...) = segment_image(color.(img); kwargs...)
 
 """
-    idx = stimulus_index(seg::SegmentedImage, colorproj = RGB(1, 1, -2))
+    idx = stimulus_index(seg::SegmentedImage, centroidsacc; colorproj = RGB(1, 1, -2), expectedloc = nothing)
 
 Given a segmented image `seg`, return the index of the segment that scores
 highest on the product of (1) projection (dot product) with `colorproj` and (2)
 number of pixels.
+
+Optionally, if images were taken with a fixed location for the stimulus, a segment's score
+is divided by the squared distance of its centroid (via `centroidsacc`) from the position given by `expectedloc`.
 """
-function stimulus_index(seg::SegmentedImage, colorproj = RGB{Float32}(1, 1, -2))
-    proj = [l => (colorproj ⋅ segment_mean(seg, l)) * segment_pixel_count(seg, l) for l in segment_labels(seg)]
-    (i, _) = argmax(last, proj)
-    return i
+function stimulus_index(seg::SegmentedImage, centroidsacc; colorproj = RGB{Float32}(1, 1, -2), expectedloc = nothing)
+    if !isnothing(expectedloc)
+        proj = map(segment_labels(seg)) do l
+            l == 0 && return 0
+            val = centroidsacc[l]
+            centroid = [round(Int, val[1] / val[3]), round(Int, val[2] / val[3])]
+            return l => (colorproj ⋅ segment_mean(seg, l) * segment_pixel_count(seg, l) / max(1, sum(abs2, centroid .- expectedloc)))
+        end
+        (i, _) = argmax(last, proj)
+        return i
+    else 
+        proj = [l => (colorproj ⋅ segment_mean(seg, l)) * segment_pixel_count(seg, l) for l in segment_labels(seg)]
+        (i, _) = argmax(last, proj)
+        return i
+    end
 end
 
 # function contiguous(seg::SegmentedImage, img::AbstractMatrix{<:Color}; min_size::Int = 50)
@@ -47,6 +61,35 @@ end
 # contiguous(seg::SegmentedImage, img::AbstractMatrix{<:Colorant}; kwargs...) =
 #     contiguous(seg, color.(img); kwargs...)
 
+"""
+    centroidsacc, nadj = get_centroidsacc(seg::SegmentedImage)
+
+Given a the index map `indexmap` of a segmented image, return an accumulator for each segment's centroid
+as well as the number of times two segments are adjacent.
+"""
+function get_centroidsacc(indexmap::Matrix{Int64})
+    keypair(i, j) = i < j ? (i, j) : (j, i)
+    R = CartesianIndices(indexmap)
+    Ibegin, Iend = extrema(R)
+    I1 = oneunit(Ibegin)
+    centroidsacc = Dict{Int, Tuple{Int, Int, Int}}()   # accumulator for centroids
+    nadj = Dict{Tuple{Int, Int}, Int}()             # number of times two segments are adjacent
+    for idx in R
+        l = indexmap[idx]
+        l == 0 && continue
+        acc = get(centroidsacc, l, (0, 0, 0))
+        centroidsacc[l] = (acc[1] + idx[1], acc[2] + idx[2], acc[3] + 1)
+        for j in max(Ibegin, idx - I1):min(Iend, idx + I1)
+            lj = indexmap[j]
+            if lj != l && lj != 0
+                k = keypair(l, lj)
+                nadj[k] = get(nadj, k, 0) + 1
+            end
+        end
+    end
+    return centroidsacc, nadj
+end
+
 struct Spot
     npixels::Int
     centroid::Tuple{Int, Int}
@@ -67,40 +110,21 @@ Spots larger than `max_size_frac * npixels` (default: 10% of the image) are igno
 function spots(
         seg::SegmentedImage;
         max_size_frac=0.1,            # no spot is bigger than max_size_frac * npixels
+        kwargs...
     )
-    keypair(i, j) = i < j ? (i, j) : (j, i)
-
-    istim = stimulus_index(seg)
+    centroidsacc, nadj = get_centroidsacc(seg.image_indexmap)
+    istim = stimulus_index(seg, centroidsacc; kwargs...)
 
-    label = seg.image_indexmap
-    R = CartesianIndices(label)
-    Ibegin, Iend = extrema(R)
-    I1 = oneunit(Ibegin)
-    centroidsacc = Dict{Int, Tuple{Int, Int, Int}}()   # accumulator for centroids
-    nadj = Dict{Tuple{Int, Int}, Int}()             # number of times two segments are adjacent
-    for idx in R
-        l = label[idx]
-        l == 0 && continue
-        acc = get(centroidsacc, l, (0, 0, 0))
-        centroidsacc[l] = (acc[1] + idx[1], acc[2] + idx[2], acc[3] + 1)
-        for j in max(Ibegin, idx - I1):min(Iend, idx + I1)
-            lj = label[j]
-            if lj != l && lj != 0
-                k = keypair(l, lj)
-                nadj[k] = get(nadj, k, 0) + 1
-            end
-        end
-    end
     stimulus = Ref{Pair{Int,Spot}}()
     filter!(centroidsacc) do (key, val)
         if key == istim
             stimulus[] = key => Spot(val[3], (round(Int, val[1] / val[3]), round(Int, val[2] / val[3])))
             return false
         end
-        val[3] <= max_size_frac * length(label) || return false
+        val[3] <= max_size_frac * length(seg.image_indexmap) || return false
         # # is the centroid within the segment?
         # x, y = round(Int, val[1] / val[3]), round(Int, val[2] / val[3])
-        # l = label[x, y]
+        # l = seg.image_indexmap[x, y]
         # @show l
         # l == key || return false
         # is the segment lighter than most of its neighbors?