Merge pull request #35 from tlnagy/tn/add-encapsulation-support

Improve encapsulation support

.github/workflows/CI.yml

@@ -16,7 +16,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.4'
+          - '1.6'
           - '1'
           - 'nightly'
         python-version: ['3']

Project.toml

@@ -1,7 +1,7 @@
 name = "FluorescenceExclusion"
 uuid = "b1480678-572a-4d60-ab1d-db2f300f29c9"
 authors = ["Tamas Nagy <[email protected]>"]
-version = "0.2.4"
+version = "0.2.5"
 
 [deps]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"

docs/src/index.md

@@ -21,7 +21,7 @@ img = TiffImages.load(path) # load a FxM image
 
 fimg = float.(img) # convert to floating point
 
-correct!(fimg)
+correct!(fimg, fimg)
 
 hcat(img, fimg)
 ```

src/correct.jl

@@ -35,20 +35,22 @@ julia> img = TiffImages.load(path);
 
 julia> fimg = float.(img); #convert to Gray{Float32}
 
-julia> correct!(fimg); # correct FxM channel
+julia> correct!(fimg, fimg); # correct FxM channel, we don't have a denoise image so we re-use the same one
 ```
 """
-function correct!(data::A; mask=falses(size(data)...)) where {A <: AbstractMatrix{<: AbstractGray{<: AbstractFloat}}}
-    pillar_mask, pillar_centers = identify(Pillars(), data)
-    cell_mask = identify(Cells(), data, pillar_mask) .| mask
+function correct!(data::A, dn_img) where {A <: AbstractMatrix{<: AbstractGray{<: AbstractFloat}}}
+    pillar_mask, pillar_centers = identify(Pillars(), dn_img)
+    cell_mask = identify(Cells(), dn_img, pillar_mask)
 
     bkg = darkfield(data, pillar_centers)
 
-    flatfield = compute_flatfield(data, cell_mask .| pillar_mask, len=20) .- bkg;
+    data .-= bkg
+
+    flatfield = compute_flatfield(data, cell_mask .| pillar_mask, len=30);
 
     labeled = label_components(cell_mask)
 
-    data .= @. clamp((data - bkg) / flatfield, Gray(-0.05), Gray(1.05))
+    data .= @. clamp(data / flatfield, Gray(-0.1), Gray(1.1))
 
     labeled
 end

src/particles.jl

@@ -10,7 +10,8 @@ distance away in pixels from each cell to include in its local background
 calculation.
 """
 function build_tp_df(img::AxisArray{T1, 4},
-                     components::AxisArray{T2, 3}; dist=(2, 10)) where {T1, T2 <: Integer}
+                     components::AxisArray{T2, 3}; dist=(2, 10), 
+                     fxmchannel = :EPI_DeepRed, min_pillar_dist = 10) where {T1, T2 <: Integer}
 
     particles = DataFrames.DataFrame[]
 
@@ -38,7 +39,9 @@ function build_tp_df(img::AxisArray{T1, 4},
         # select the current time slice and enforce storage order to match
         # components
         slice = view(img, Axis{:y}(:), Axis{:x}(:), Axis{:channel}(:), Axis{:time}(timepoint))
-        localities = identify(Locality(), component_slice, ids, dist=dist)
+        pillars, _ = identify(Pillars(), view(slice, Axis{:channel}(fxmchannel)))
+        localities = identify(Locality(), component_slice, pillars;
+                              dist=dist, min_pillar_dist=min_pillar_dist)
         # dictionary of ids to areas
         data = OrderedDict(:x=>xs,
                            :y=>ys,
@@ -84,7 +87,8 @@ function build_tp_df(img::AxisArray{T1, 3},
     build_tp_df(AxisArray(reshape(img, size(img)..., 1),
                           AxisArrays.axes(img)..., Axis{:channel}([:FxM])),
                 thresholds;
-                dist=dist
+                dist=dist,
+                fxmchannel=:FxM
                )
 end
 
@@ -103,7 +107,8 @@ function build_tp_df(img::AxisArray{T1, 3},
     build_tp_df(AxisArray(reshape(img, size(img)..., 1),
                           AxisArrays.axes(img)..., Axis{:channel}([:FxM])),
                 thresholds;
-                dist=dist
+                dist=dist,
+                fxmchannel=:FxM
                )
 end
 

src/segment.jl

@@ -62,14 +62,19 @@ thresholds using the distribution of spatial gradient sizes across cells.
 """
 function identify(::Cells,
                  img::AbstractArray{T, 2},
-                 pillar_masks::AbstractArray{Bool, 2}) where {T}
+                 pillar_masks::AbstractArray{Bool, 2};
+                 min_pixel = 50,
+                 min_pillar_dist = 30) where {T}
 
-    grad_y, grad_x = imgradients(img, KernelFactors.scharr);
+    grad_y, grad_x = imgradients(img, KernelFactors.scharr, "reflect");
     mag = hypot.(grad_x, grad_y)
 
     # remove the pillars and their vicinities from the calculations since they
     # can skew the results
-    foreground = distance_transform(feature_transform((pillar_masks))) .> 30
+    foreground = distance_transform(feature_transform((pillar_masks))) .> min_pillar_dist
+    nx, ny = size(img)
+    foreground[1:ny, [1, 2, ny-1, ny]] .= false
+    foreground[[1, 2, nx-1, nx], 1:nx] .= false
     mag .*= foreground
 
     flattened = filter(x->x > 0.0, reshape(mag, :))
@@ -82,7 +87,7 @@ function identify(::Cells,
     normfit = fit_mle(LogNormal, view(flattened, lo .< flattened .< hi))
 
     thres = opening(.~imfill(mag .< exp(normfit.μ + 1*normfit.σ), (0, 500)))
-    imfill(thres, (0, 50))
+    imfill(thres, (0, min_pixel))
 end
 
 function identify(::Cells,
@@ -117,28 +122,28 @@ function get_locality(foreground::AbstractArray{Bool, 2}, cellmask::AbstractArra
 end
 
 """
-    identify(Locality(), labels, ids; dist) -> Dict
+    identify(Locality(), labels, pillar_masks; min_pillar_dist, dist) -> Dict
 
-Gets the local background areas given the output of `label_components` and a
-list of objects whose areas should be found defined by `ids`. The local
+Gets the local background areas given the output of `label_components`. The local
 background is defined as a ring around the object that is at
 minimum `mindist` away from every object (in pixels) and a maximum of `maxdist`
-away from the target object. The result is a dictionary mapping the object id to
-all indices in its local background.
+away from the target object. Pillars, as defined by true values in
+`pillar_masks`, are ignored as are areas that are `min_pillar_dist` away from a
+pillar. The result is a dictionary mapping the object id to all indices in its
+local background. 
+
 !!! warning
     It's really important that `labels` has all foreground objects labeled
     because otherwise they might inadvertently be included in an object's
     locality.
 """
 function identify(::Locality, 
-                  labels::AbstractArray{Int, 2}, 
-                  ids::Vector{Int};
+                  labels::AbstractMatrix{Int}, 
+                  pillar_masks::AbstractMatrix{Bool};
+                  min_pillar_dist = 30,
                   dist=(2, 10))
 
     boxes = component_boxes(labels)
-    # boxes always contains the background label as well, while there's no guarantee
-    # that ids does as well
-    (!(0 in ids)) && deleteat!(boxes, 1)
 
     # if we detect gaps in the numbering it's probably because some labels were
     # removed, which can make our foreground detection incorrect
@@ -147,26 +152,33 @@ function identify(::Locality,
         @warn "Do not remove values from `labels`, only from `ids`. Calc might be wrong!"
     end
 
+    pillar_area = distance_transform(feature_transform((pillar_masks))) .<= min_pillar_dist
+    nx, ny = size(labels)
+    pillar_area[1:nx, [1, 2, ny-1, ny]] .= true
+    pillar_area[[1, 2, nx-1, nx], 1:ny] .= true
+
     nx, ny = size(labels)
     localities = OrderedDict{Int, Vector{CartesianIndex{2}}}()
 
     allobjects = labels .> 0
 
     δ = dist[1] + dist[2]
 
-    for id in ids
+    for id in computed_ids
         # unpack bounding box
-        (minx, miny), (maxx, maxy) = boxes[id]
+        (minx, miny), (maxx, maxy) = boxes[id+1]
 
         # extend window around bounding box by the max distance
         xrange = max(minx-δ, 1):min(maxx+δ, nx)
         yrange = max(miny-δ, 1):min(maxy+δ, ny)
 
         local_allobjects = view(allobjects, xrange, yrange)
         local_cellmask = view(labels, xrange, yrange) .== id
-        @assert sum(local_cellmask) > 0 "No object of $id found"
+        @assert sum(local_cellmask) > 0 "No object of id=$id found"
 
         locality = get_locality(local_allobjects, local_cellmask; dist=dist)
+        localpillar = view(pillar_area, xrange, yrange)
+        locality[localpillar] .= false
 
         # use offset arrays to return CartesianIndices in the original image coordinates
         localities[id] = findall(OffsetArray(locality, xrange, yrange))
@@ -223,7 +235,7 @@ function isencapsulated(locality::Vector{CartesianIndex{2}})
     # extremely thin localities where the area of the convex hull is very close
     # to the area of the hole
     hulltmp = dilate(tmp)
-    (sum(hulltmp) < 3) && return false
+    (count(tmp) <= 3) && return false
     hullc = convexhull(hulltmp)
     # area contained within the convex hull, we need this to set the 
     # maximum size of the flood fill algorithm
@@ -239,5 +251,5 @@ function isencapsulated(locality::Vector{CartesianIndex{2}})
     # flood fill any regions that are up to the area of the convex hull of
     # the locality and then check if anything has changed. If it has, that
     # means the locality fully encapsulated the cell.
-    any(imfill(tmp2, 1:hullarea) .⊻ tmp2)
+    any(imfill(tmp2, (1, hullarea)) .⊻ tmp2)
 end

test/Project.toml

@@ -3,6 +3,7 @@ ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a"
 ImageDraw = "4381153b-2b60-58ae-a1ba-fd683676385f"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"

test/segment.jl

@@ -1,4 +1,5 @@
 using ImageDraw
+using FixedPointNumbers
 using FluorescenceExclusion: isencapsulated
 
 @testset "Segmenting pillars" begin
@@ -58,6 +59,30 @@ end
 
 end
 
+@testset "Segment locality" begin
+    labels = Gray.(zeros(UInt8, 100, 200))
+
+    # draw two cells
+    draw!(labels, Ellipse(CirclePointRadius(50, 50, 20)), Gray(reinterpret(N0f8, UInt8(1))))
+    draw!(labels, Ellipse(CirclePointRadius(120, 50, 20)), Gray(reinterpret(N0f8, UInt8(2))))
+
+    # draw pillar
+    pillar_mask = Gray.(falses(size(labels)))
+    draw!(pillar_mask, Ellipse(CirclePointRadius(200, 50, 50)))
+
+    localities = identify(FluorescenceExclusion.Locality(), Int.(reinterpret(UInt8, labels)), Bool.(pillar_mask))
+    @test length(localities[1]) == 1244
+    @test length(localities[2]) == 681 # cell 2 is close to the pillar and thus its locality is smaller
+
+    # if we change it so there's no padding then the two localities should be same
+    localities = identify(FluorescenceExclusion.Locality(), Int.(reinterpret(UInt8, labels)), Bool.(pillar_mask), min_pillar_dist = 0)
+    @test length(localities[1]) == length(localities[2])
+
+    # test that the function complains if a label is missing
+    draw!(labels, Ellipse(CirclePointRadius(120, 50, 20)), Gray(reinterpret(N0f8, UInt8(4))))    
+    @test_logs (:warn, "Do not remove values from `labels`, only from `ids`. Calc might be wrong!") identify(FluorescenceExclusion.Locality(), Int.(reinterpret(UInt8, labels)), Bool.(pillar_mask))
+end
+
 @testset "Encapsulation" begin
     img = Gray.(falses(100, 100))
 
@@ -98,6 +123,8 @@ end
     # if there are <3 points in the locality than the convex hull calculation
     # fails, lets make sure we handle this properly
     img .= false
+    img[10, 10] = true
+    img[10, 20] = true
     @test !isencapsulated(findall(imgr))