build a family tree for all traps
H QIn
Jane 18, 2018
setup the systems
rm(list=ls())
debug = 1;
read null distribution from single cell data
moving_buffer = read.csv("data/moving_buffer.csv")
p_dist = function(x, input_buffer) {
input_buffer = input_buffer[ !is.na(input_buffer)];
lower = input_buffer[ input_buffer < x]; lower = lower[!is.na(lower)];
upper = input_buffer[input_buffer > x]; upper = upper[!is.na(upper)];
LowerP = (length(lower) / length(input_buffer))
UpperP = (length(upper) / length(input_buffer))
return(list(LowerP,UpperP, mean(input_buffer, na.rm=T), sd(input_buffer, na.rm=T)))
}
my_dist = function(x1, y1, x2, y2) {sqrt( (x1 - x2)^2 + (y1-y2)^2 ) };
Load the segmented data
DataTrap <- read.csv("data/DataSet_v6.csv",header=TRUE,stringsAsFactors = FALSE)
Take one trap for the time beling to develop family tree
tb.alltraps = DataTrap;
tb.alltraps$predecessorID = NA;
tb.alltraps$successorID = NA;
traps = sort(unique(tb.alltraps$trap_num));
Work from bottom to top to build a family tree using trace-back.
for( cur_trap in 1:length(traps)){
# cur_trap = 10;
#for( cur_trap in c(20,30,50)){
tb = tb.alltraps[tb.alltraps$trap_num==traps[cur_trap], ];
timpoints = sort(unique(tb$Image_num))
#plot( timpoints, 1:length(timpoints), type='l') #visual check timepoints.
for (i in length(timpoints):2){
if(debug > 1) { print(paste("i=",i)); }
current_frame = tb[ tb$Image_num == timpoints[i], ];
previous_frame = tb[ tb$Image_num == timpoints[i-1], ];
d.tb = p.d.tb = data.frame(matrix(nrow=length(current_frame[,1]), ncol=length(previous_frame[,1]) ));
for ( i.cur in 1:length(current_frame[,1]) ) {
for ( j.pre in 1:length(previous_frame[,1]) ) {
d.tb[i.cur, j.pre] = my_dist(current_frame$cell_X[i.cur], current_frame$cell_Y[i.cur],
previous_frame$cell_X[j.pre], previous_frame$cell_Y[j.pre]);
#my_dist = function(x1, y1, x2, y2) {sqrt( (x1 - x2)^2 + (y1-y2)^2 ) };
p.d.tb[i.cur, j.pre] = p_dist(d.tb[i.cur, j.pre], moving_buffer$d )[[2]];
}
current_frame$predecessorID[i.cur] = which.max( p.d.tb[i.cur, ]); #pick max, greedy algorithm.
}
tb$predecessorID[tb$Image_num == timpoints[i]] = current_frame$predecessorID; #assume order unchanged.
}
tb.alltraps$predecessorID[tb.alltraps$trap_num==traps[cur_trap]] = tb$predecessorID;
} tb.alltraps$predecessorID[tb.alltraps$trap_num==traps[30]]; ## [1] NA 1 1 1 1 1 1 2 1 2 1 2 1 2 2 1 1 2 1 2 2 1 1
## [24] 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2
## [47] 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 1 1 2 1
## [70] 2 1 2 1 2 1 2 1 2 1 2 1 2 1 1 1 1 1 2 1 2 1 2
## [93] 1 2 1 2 2 1 2 1 1 1 1 1 2 1 2 1 2 1 2 1 1 1 1
## [116] 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 1
## [139] 2 2 1 1 1 2 1 2 1 2 1 2 1 2 2 1 2 1 1 2 1 2 1
## [162] 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 1
## [185] 1 1 2 1 2 1 2 2 1 1 2 2 1 1 1 2 1 2 1 2 1 2 1
## [208] 1 1 1 1 2 1 2 1 2 2 1 2 1 1 1 2 1 2 1 2 1 2 1
## [231] 2 1 2 1 1 1 1 2 1 2 1 2 1 2 1 1 1 1 1 2 1 2 1
## [254] 2 1 2 1 2 1 1 1 1 2 1 2 1 2 1 2 1 1 1 1 1 2 1
## [277] 2 1 2 1 1 1 1 1 1 2 1 2 1 2 1 2 1 1 1 1 1 2 1
## [300] 2 1 2 1 2 1 1 1 1 1 1 2 1 2 1 2 1 2 1 2 1 2 1
## [323] 1 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 1 1 1 1 2 1
## [346] 2 1 2 1 2 1 2 1 1 1 1 1 1 2 1 2 1 2 1 2 1 2 1
## [369] 2 1 2 1 2 1 1 1 1 1 1 1 1 2 1 2 1 2 1 2 1 2 1
## [392] 2 1 2 1 2 1 1 1 1 1 1 1 1 2 1 2 1 2 1 2 1 2 1
## [415] 2 1 2 1 2 1 1 1 1 1 1 2 1 2 3 1 2 3 1 2 3 1 2
## [438] 3 1 2 3 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [461] 1 1 1 1 1 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## [484] 1 2 3 1 2 1 2 1 1 2 1 2 1 1 2 3 1 2 1 2 1 2 1
## [507] 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 1 2 2 1
## [530] 1 1 2 1 3 4 1 2 3 4 2 1 3 4 2 1 3 4 1 2 3 4 1
## [553] 3 2 4 1 3 2 2 1 3 1 2 3 2 1 3 1 2 3 1 2 3 1 2
## [576] 3 1 1 2 3 4 2 1 3 4 2 1 3 4 1 2 3 4 1 2 3 4 1
## [599] 2 3 1 4 2 1 3 4 2 1 3 4 1 2 3 4 1 2 3 1 2 3 1
## [622] 1 2 3 1 2 3 1 2 3 1 1 2 3 4 1 2 3 4 2 1 3 4 2
## [645] 1 3 2 1 3 2 1 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## [668] 1 2 3 1 2 3 1 2 3 1 1 2 3 4 1 2 3 4 1 2 3 4 1
## [691] 2 3 3 4 1 2 1 3 4 6 5 2 1 4 6 1 2 3 3 4 1 2 4
## [714] 5 1 2 3 3 4 1 2 4 5 1 2 1 3 4 1 2 3 4 5 1 2 3
## [737] 4 5 1 2 4 5 1 2 3 4 1 2 3 4 1 2 3 3 4 1 2 4 5
## [760] 1 2 3 4 1 2 3 1 4 1 2 3 5 1 2 3 3 4 1 2 4 5 1
## [783] 2 1 3 4 1 2 3 4 5 1 2 3 4 5 1 2 4 5 1 2 3 4 1
## [806] 2 3 4 1 2 3 3 4 1 2 4 5 1 2 3 4 1 2 3 1 4 1 2
## [829] 3 5 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## [852] 1 2 3 1 1 1 1 1 2 3 4 5 6 7 1 2 3 4 4 3 5 7 1
## [875] 1 2 3 4 5 6 8 9 7 1 2 3 4 6 7 9 8 1 2 3 4 4 7
## [898] 6 1 2 3 5 6 7 1 2 3 4 4 5 6 1 2 3 5 4 6 7 1 2
## [921] 4 3 5 7 6 1 2 3 4 5 7 6 1 2 5 4 4 3 6 1 7 1 2
## [944] 3 4 6 7 9 1 2 3 4 4 1 6 7 1 2 3 4 3 7 8 1 2 3
## [967] 4 5 6 7 1 2 5 4 3 6 7 1 2 5 4 3 6 7 1 2 3 4 5
## [990] 6 7 1 2 4 3 5 4 7 1 6 1 2 3 4 5 6 7 9 8write.csv(tb.alltraps, "data/alltraps_withtree.csv", quote=F)
No comments:
Post a Comment