-
Notifications
You must be signed in to change notification settings - Fork 0
/
normalization3d.jl
129 lines (118 loc) · 4.07 KB
/
normalization3d.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
"""
Normalize brightness for each time point
Version Comment
0.1 hf, initial version
0.2 hf, don't use fitting
0.2.1 hf, medianfilter with repeat padding
"""
using Images
using Statistics
using LsqFit
using Clustering
"""
Extract feature from image
"""
function extract_feature(imgs)
z_depth, t_len = size(imgs)[3:4]
low_line = zeros(t_len)
mean_line = zeros(t_len)
@inbounds Threads.@threads for t in 1:t_len
nozeros_pixel = imgs[:,:,:,t][imgs[:,:,:,t].≠0]
if length(nozeros_pixel) == 0
#println("no cell in time point $t")
continue # skip empty stack
end
low_line[t] = minimum(nozeros_pixel)
mean_line[t] = mean(nozeros_pixel)
end
low_line, mean_line
end
medianfilter1(v, ws) = [median(v[i:(i+ws-1)]) for i=1:(length(v)-ws+1) ]
#medianfilter(v, ws) = medianfilter1(vcat(0,v,0),ws)
medianfilter(v, ws) = medianfilter1(vcat(v[1], v, v[end]),ws) # repeat
"""
Use k-means to remove outliers
Input raw data series, return index without outliers
"""
function kmeans_rm_outliers( raw_line )
# TODO: if only one cluster?
t_len = length(raw_line)
X = zeros(2, t_len)
X[1, :] = fill(0.02, t_len) # just choose 0.02 in random.
X[2, :] = raw_line; # forgive me just add one more dimension for k-means
result = kmeans(X, 2; init=:kmpp)
if result.counts[1]> result.counts[2] # chose the one have more points
index = 1
else
index = 2
end
t_nooutliers = [ i for i in 1:t_len if result.assignments[i] == index ]
end
"""
Removing outliers and smoothing feature
"""
function correct_feature( low_line, mean_line, fitted::Bool=false )
t_len = length(low_line)
ws = 3 # use small window to remove spark
low_line = medianfilter(low_line, ws)
mean_line = medianfilter(mean_line, ws)
# TODO: if only one cluster? Add additional judge by std(mean_line[1:end-1].-mean_line[2:end])
if fitted
t_nooutliers = kmeans_rm_outliers(mean_line)
mean_line_nooutliers = mean_line[t_nooutliers]
# assume mean_line and low_line share same shape
low_line_nooutliers = low_line[t_nooutliers]
@. model(x, p) = p[1]*exp(-x*p[2]) # fit with expontional function
fit_mean = curve_fit(model, t_nooutliers, mean_line_nooutliers, [0.003, 0.0003])
fit_low = curve_fit(model, t_nooutliers, low_line_nooutliers, [0.002, 0.0003])
corrected_mean_line = model(1:t_len, fit_mean.param)
corrected_low_line = model(1:t_len, fit_low.param)
# TODO: what if mean.(fit, no_outliers) ?
else
corrected_mean_line = mean_line
corrected_low_line = low_line
end
corrected_low_line, corrected_mean_line
end
"Normalize nucleus to uniform intensity"
function normalize(imgs)
#z_depth = 20
#t_len = size(imgs)[3] ÷ z_depth
#d1_len, d2_len = size(imgs)[1], size(imgs)[2]
d1_len, d2_len, z_depth, t_len = size(imgs)
low_line, mean_line = extract_feature(imgs)
corrected_low_line, corrected_mean_line = correct_feature(low_line, mean_line)
imgs_norm = zeros(N0f16, size(imgs))
imgs_len = length(imgs_norm[:, :, :, 1])
@inbounds Threads.@threads for t in 1:t_len
"""
z = (t-1)*20+1 : 20*t
imgs_norm[:,:,z] = ((imgs[:,:,z] .- corrected_low_line[t]) ./ (corrected_mean_line[t]-corrected_low_line[t]).*0.005 .+ 0.02) .* ( imgs[:, :, z] .>0)
#print("$t ")
"""
tmp = view(imgs, :, :, :, t)
tmp_norm = view(imgs_norm, :, :, :, t)
for i in 1:imgs_len
if tmp[i] > 0
tmp_norm_value = (tmp[i] - corrected_low_line[t]) /
(corrected_mean_line[t] - corrected_low_line[t])*0.005 + 0.02
if tmp_norm_value < 0
tmp_norm[i] = 0
elseif tmp_norm_value > 1
tmp_norm[i] = 1
else
tmp_norm[i]= tmp_norm_value
end
end
end
end
imgs_norm, Dict("low_line"=>low_line, "mean_line"=>mean_line,
"corrected_low_line"=>corrected_low_line, "corrected_mean_line"=>corrected_mean_line)
end
#@time s14c1 = load(File(format"TIFF", "$data_dir/$(cell_name[14])"));
#@time imgs_norm, low_line, mean_line = normalize(s14c1);
#plot(low_line)
#plot!(low_line_raw, marker=(:star))
#plot!(mean_line)
#plot!(mean_line_raw, marker=(:dot))
#plot!(ylim=(0,0.004))