Revision 27fa41066bf09a3c1ca00fd333f72755e481b16f authored by H.F on 13 April 2021, 10:57 UTC, committed by H.F on 13 April 2021, 10:57 UTC
1 parent f3455fe
segmentation3d.jl
``````#using FileIO
using Images
using Statistics

"""
Version Comment
v0.1	hf, use yen-threshold select threshold for 3D stack at each time point
"""

#zstack = imfilter(zstack, Kernel.gaussian((2,2,2)))
#z_depth = 20
z_depth = size(zstack)[3]
#thresholds_z = [real(yen_threshold(zstack[:, :, i])) for i in 1:20]
local thresholds_z = zeros(Float64, z_depth)
@inbounds for z in 1:z_depth
zstack[:,:,z] = imfilter(zstack[:, :, z], Kernel.gaussian(2))
local select_pixels = zstack[81:end-81, 81:end-81,  z] .> 0.001
thresholds_z[z]
if sum(select_pixels) > 2e2 #TODO
thresholds_z[z] = real(otsu_threshold(
zstack[81:end-81,81:end-81,z][select_pixels]))
end
end
#threshold_3d = real(yen_threshold(imfilter(zstack, Kernel.gaussian((2,2,1)))))
#threshold_3d = real(yen_threshold(zstack))
threshold_3d = maximum(thresholds_z) # choose clearest one
end

function extract3dnucleus(stack)
#z_depth = 20
#t_len = size(stack)[3] ÷z_depth
z_depth,t_len = size(stack)[3:4]
local nucleus = zeros(Gray{Normed{UInt16,16}}, size(stack))
local nucleus_len = length(stack)
local thresholds = zeros(Float64, t_len)
println("Extracting nucleus")
#local z = (t-1)*20+1 : 20*t
nucleus_t = view(nucleus, :, :, :, t)
stack_t = view(stack, :, :, :, t)
#nucleus[:, :, z] = nucleus_3dmask .* stack[:, :, z]
nucleus_t[i] = stack_t[i]
end
end
end
end

"""
Just remove small regions in binary images
"""
selected_con = (0:length(con_size)-1)[con_size .> 2e4]
else
end
end
end

"""
Computes threshold for grayscale image using Otsu's method. (modify from Imagas.jl)
Parameters:
-    img         = Grayscale input image
-    bins        = Number of bins used to compute the histogram. Needed for floating-point images.
"""
function otsu_threshold2(img::AbstractArray{T, N}, min, bins::Int = 256) where {T<:Union{Gray,Real}, N}

#min, max = extrema(img)
max = maximum(img)
if max == 0 # avoid image without cell
return T(0)
end
edges, counts = imhist(img, range(gray(min), stop=gray(max), length=bins))
if sum(counts)<2e2 # don't count small region
return T(0)
end
histogram = counts./sum(counts)

ω0 = 0
μ0 = 0
μt = 0
μT = sum((1:(bins+1)).*histogram)
max_σb=0.0
thres=1

for t in 1:bins
ω0 += histogram[t]
ω1 = 1 - ω0
μt += t*histogram[t]

σb = (μT*ω0-μt)^2/(ω0*ω1)

if(σb > max_σb)
max_σb = σb
thres = t
end
end

return T((edges[thres-1]+edges[thres])/2)
end