|
| 1 | +""" |
| 2 | + seg = segment_image(img; threshold=0.1, min_size=20) |
| 3 | +
|
| 4 | +Given an image `img`, segment it into regions using a region growing algorithm. |
| 5 | +`min_size` is the minimum number of pixels per segment, and `threshold` determines |
| 6 | +how different two colors must be to be considered different segments. |
| 7 | +Larger `threshold` values will result in fewer segments. |
| 8 | +""" |
| 9 | +function segment_image( |
| 10 | + img::AbstractMatrix{<:Color}; |
| 11 | + threshold::Real = 0.2, # threshold for color similarity in region growing |
| 12 | + prune::Bool = true, # prune small segments |
| 13 | + min_size::Int = 50, # minimum size of segments to keep |
| 14 | + ) |
| 15 | + seg = unseeded_region_growing(img, threshold) |
| 16 | + L = label_components(labels_map(seg)) # insist on contiguous regions |
| 17 | + seg = SegmentedImage(img, L) |
| 18 | + if prune |
| 19 | + println("Pruning segments smaller than $min_size pixels") |
| 20 | + seg = prune_segments(seg, label -> segment_pixel_count(seg, label) < min_size, (l1, l2) -> colordiff(segment_mean(seg, l1), segment_mean(seg, l2))) |
| 21 | + end |
| 22 | + return seg |
| 23 | +end |
| 24 | +segment_image(img::AbstractMatrix{<:Colorant}; kwargs...) = segment_image(color.(img); kwargs...) |
| 25 | + |
| 26 | +""" |
| 27 | + idx = stimulus_index(seg::SegmentedImage, colorproj = RGB(1, 1, -2)) |
| 28 | +
|
| 29 | +Given a segmented image `seg`, return the index of the segment that scores |
| 30 | +highest on the product of (1) projection (dot product) with `colorproj` and (2) |
| 31 | +number of pixels. |
| 32 | +""" |
| 33 | +function stimulus_index(seg::SegmentedImage, colorproj = RGB{Float32}(1, 1, -2)) |
| 34 | + proj = [l => (colorproj ⋅ segment_mean(seg, l)) * segment_pixel_count(seg, l) for l in segment_labels(seg)] |
| 35 | + (i, _) = argmax(last, proj) |
| 36 | + return i |
| 37 | +end |
| 38 | + |
| 39 | +struct Spot |
| 40 | + npixels::Int |
| 41 | + centroid::Tuple{Int, Int} |
| 42 | +end |
| 43 | + |
| 44 | +""" |
| 45 | + spotdict, stimulus = spots(seg; max_size_frac=0.1) |
| 46 | +
|
| 47 | +Given a segmented image `seg`, return a `Dict(idx => spot)` where `idx` is the segment index |
| 48 | +and `spot` is a `Spot` object where `spot.npixels` is the number of pixels in the segment |
| 49 | +and `spot.centroid` is the centroid of the segment. |
| 50 | +
|
| 51 | +`stimulus` is a `Pair{Int, Spot}` where the first element is the index of the |
| 52 | +stimulus segment and the second element is the `Spot` object for that segment. |
| 53 | +
|
| 54 | +Spots larger than `max_size_frac * npixels` (default: 10% of the image) are ignored. |
| 55 | +""" |
| 56 | +function spots( |
| 57 | + seg; |
| 58 | + max_size_frac=0.1, # no spot is bigger than max_size_frac * npixels |
| 59 | + ) |
| 60 | + keypair(i, j) = i < j ? (i, j) : (j, i) |
| 61 | + |
| 62 | + istim = stimulus_index(seg) |
| 63 | + |
| 64 | + label = seg.image_indexmap |
| 65 | + R = CartesianIndices(label) |
| 66 | + Ibegin, Iend = extrema(R) |
| 67 | + I1 = one(Ibegin) |
| 68 | + centroidsacc = Dict{Int, Tuple{Int, Int, Int}}() # accumulator for centroids |
| 69 | + nadj = Dict{Tuple{Int, Int}, Int}() # number of times two segments are adjacent |
| 70 | + for idx in R |
| 71 | + l = label[idx] |
| 72 | + l == 0 && continue |
| 73 | + acc = get(centroidsacc, l, (0, 0, 0)) |
| 74 | + centroidsacc[l] = (acc[1] + idx[1], acc[2] + idx[2], acc[3] + 1) |
| 75 | + for j in max(Ibegin, idx - I1):min(Iend, idx + I1) |
| 76 | + lj = label[j] |
| 77 | + if lj != l && lj != 0 |
| 78 | + k = keypair(l, lj) |
| 79 | + nadj[k] = get(nadj, k, 0) + 1 |
| 80 | + end |
| 81 | + end |
| 82 | + end |
| 83 | + stimulus = Ref{Pair{Int,Spot}}() |
| 84 | + filter!(centroidsacc) do (key, val) |
| 85 | + if key == istim |
| 86 | + stimulus[] = key => Spot(val[3], (round(Int, val[1] / val[3]), round(Int, val[2] / val[3]))) |
| 87 | + return false |
| 88 | + end |
| 89 | + val[3] <= max_size_frac * length(label) || return false |
| 90 | + # # is the centroid within the segment? |
| 91 | + # x, y = round(Int, val[1] / val[3]), round(Int, val[2] / val[3]) |
| 92 | + # l = label[x, y] |
| 93 | + # @show l |
| 94 | + # l == key || return false |
| 95 | + # is the segment lighter than most of its neighbors? |
| 96 | + dcol, ncol = zero(valtype(seg.segment_means)), 0 |
| 97 | + for (k, n) in nadj |
| 98 | + if key == k[1] || key == k[2] |
| 99 | + l1, l2 = k[1], k[2] |
| 100 | + if l1 == key |
| 101 | + l1, l2 = l2, l1 |
| 102 | + end |
| 103 | + dcol += n * (segment_mean(seg, l1) - segment_mean(seg, l2)) |
| 104 | + ncol += n |
| 105 | + end |
| 106 | + end |
| 107 | + return reducec(+, dcol) < 0 |
| 108 | + end |
| 109 | + return Dict(l => Spot(val[3], (round(Int, val[1] / val[3]), round(Int, val[2] / val[3]))) for (l, val) in centroidsacc), stimulus[] |
| 110 | +end |
| 111 | + |
| 112 | +""" |
| 113 | + spotdict, stimulus = upperleft(spotdict::AbstractDict{Int, Spot}, stimulus, imgsize) |
| 114 | +
|
| 115 | +Given a `spotdict` of `Spot` objects and a `stimulus` segment, return a new |
| 116 | +`spotdict` where the centroids of the spots are flipped so that the stimlus spot |
| 117 | +is in the upper left corner. |
| 118 | +""" |
| 119 | +function upperleft(spotdict::AbstractDict{Int, Spot}, stimulus, imgsize) |
| 120 | + sidx, ss = stimulus |
| 121 | + midpoint = imgsize .÷ 2 |
| 122 | + c1, c2 = ss.centroid .< midpoint |
| 123 | + imsz1, imsz2 = imgsize |
| 124 | + |
| 125 | + function flip(spot::Spot) |
| 126 | + x1, x2 = spot.centroid |
| 127 | + return Spot(spot.npixels, (c1 * x1 + (1 - c1) * (imsz1 - x1), c2 * x2 + (1 - c2) * (imsz2 - x2))) |
| 128 | + end |
| 129 | + return Dict(k => flip(v) for (k, v) in spotdict), sidx => flip(ss) |
| 130 | +end |
0 commit comments