# setwd("~/Desktop/laura/")
setwd("D:/Projects/Pattern_matter2.0/2D")

####Read in data ####
df <- read.csv("all_root_data_2D.csv") #single root data
summary(df)
gt<-read.csv("ground_truth_data_2D.csv") # ground truth per image
meas<-read.table("all_analyses_2D.txt",header = T, #the image analyses data
                 sep="\t", 
                 strip.white = TRUE
) 


df$dia_length<-df$diameter*df$length #needed for wm diameter
meas$dpi<-as.numeric(substr(meas$RHIZO.2013e,1,4)) #
meas$image<-substr(meas$RHIZO.2013e,6,35)
summary(meas)
head(meas)

Suma<-aggregate(length ~ image, df, sum) #summarize the data by image
names(Suma)[2]<-"Sum_length"
Suma$diam<-aggregate(dia_length~image,df,sum)[,2]/Suma$Sum_length#length weighted mean diameter
Suma$N<-aggregate(surface ~ image, df, length)[,2] #Number of root segments for sd
df<-merge(df,Suma,by="image") #add the total length and length weighted mean diameter per image to each root segment
Suma$diam_sd<-sqrt(aggregate((diameter-diam)^2*length~image,df,sum)[,2] 
                   /(((Suma$N-1)*Suma$Sum_length)/Suma$N))#length weighted sd of diameter

Suma$vol<-aggregate(volume ~ image, df, sum)[,2]
Suma$sufa<-aggregate(surface ~ image, df, sum)[,2]

gt_o<-gt[,c(1,20)] #overlap

all_D1<-merge(meas,Suma,by="image")
all_D<-merge(all_D1,gt_o,by="image") # add overlap
names(all_D)

# correct for different units, SufArea not corrected because not used
# diamter mm, length+vol cm
all_D$IJ_Eq_Vol<-all_D$IJ_Eq_Vol/1000 #mm3 to cm3
all_D$vol<-all_D$vol/1000
all_D$IJ_Kimura_length<-all_D$IJ_Kimura_length/10 #mm to cm
all_D$Sum_length<-all_D$Sum_length/10
all_D$CV_diam<-all_D$diam_sd/all_D$diam #coefficient of variation
all_D$density<-all_D$Sum_length/all_D$AnalysedRegionArea.cm2. #root length per image size

# Create groups

all_D$group <- "fibrous"
all_D$group[grepl("dicot", all_D$image)] <- "taproot"
all_D<-all_D[order(all_D$overlap) , ]

#### Create Subset resolution ####
all_D<-subset(all_D,dpi==400)

#### Compute rel errors for Tab 1 ####
# relative differences #
rel_diff_start<-function(Var,Ref,name) { #funktion to create summary of errors txt
  err<-(Var-Ref)/Ref
  # summary(abs(err)*100)
  # summary((err)*100)#error % 
  # 
  model <-name #the variable
  min <- as.numeric(summary((err)*100)[1]) #standarderror
  mean <- as.numeric(summary((err)*100)[4]) #standarderror
  max<- as.numeric(summary((err)*100)[6]) #adjR2
  abs<-  as.numeric(summary(abs(err)*100)[4]) #adjR2
  
  sink(file = "400errors.txt", append = T)
  d<-data.frame(model=model,min=min, mean=mean, max=max, abs=abs)
  print(d)
  sink()
  
  
}

rel_diff<-function(Var,Ref,name) { #funktion to create summary of errors txt
  err<-(Var-Ref)/Ref
  # summary(abs(err)*100)
  # summary((err)*100)#error % 
  # 
  model <-name #the variable
  min <- as.numeric(summary((err)*100)[1]) #standarderror
  mean <- as.numeric(summary((err)*100)[4]) #standarderror
  max<- as.numeric(summary((err)*100)[6]) #adjR2
  abs<-  as.numeric(summary(abs(err)*100)[4]) #adjR2
  
  sink(file = "400errors.txt", append = T)
  d<-data.frame(model=model,min=min, mean=mean, max=max, abs=abs)
  names(d)<-NULL
  print(d)
  sink()
  
  
}
# Diamter #
#WinRhizo average
all_D$d_d_WR<-(all_D$WR_AvgDiam-all_D$diam)/all_D$diam
summary(abs(all_D$d_d_WR)*100)
summary((all_D$d_d_WR)*100)#error %
rel_diff_start(all_D$WR_AvgDiam,all_D$diam,"WR_AvgDiam") #write in txt

#WinRhizo weighted mean
all_D$d_d_WR_wm<-(all_D$WR_Weighted_mean_diam-all_D$diam)/all_D$diam
summary(abs(all_D$d_d_WR_wm)*100)
summary((all_D$d_d_WR_wm)*100)
rel_diff(all_D$WR_Weighted_mean_diam,all_D$diam,"WR_Weighted_mean_diam") #write in txt

#IJ
all_D$d_d_IJ<-(all_D$IJ_Mean_Dia-all_D$diam)/all_D$diam
summary(abs(all_D$d_d_IJ)*100)
summary((all_D$d_d_IJ)*100)
rel_diff(all_D$IJ_Mean_Dia,all_D$diam,"IJ_Mean_Dia") #write in txt

#IJC
all_D$d_d_IJC<-(all_D$IJ_Mean_Dia_C-all_D$diam)/all_D$diam
summary(abs(all_D$d_d_IJC)*100)
summary((all_D$d_d_IJC)*100)
rel_diff(all_D$IJ_Mean_Dia_C,all_D$diam,"IJ_Mean_Dia_C") #write in txt

# length #
#WinRhizo
all_D$d_l_WR<-(all_D$WR_Length_Sum-all_D$Sum_length)/all_D$Sum_length
summary(abs(all_D$d_l_WR)*100)
summary((all_D$d_l_WR)*100)
rel_diff(all_D$WR_Length_Sum,all_D$Sum_length,"WR_Length_Sum") #write in txt
#IJ
all_D$d_l_IJ<-(all_D$IJ_Kimura_length-all_D$Sum_length)/all_D$Sum_length
summary(abs(all_D$d_l_IJ)*100)
summary((all_D$d_l_IJ)*100)
rel_diff(all_D$IJ_Kimura_length,all_D$Sum_length,"IJ_Kimura_length") #write in txt

# volume #
#WinRhizo correct
all_D$d_Vc_WR<-(all_D$WR_V_sum-all_D$vol)/all_D$vol
summary(abs(all_D$d_Vc_WR)*100)
summary((all_D$d_Vc_WR)*100)
rel_diff(all_D$WR_V_sum,all_D$vol,"WR_V_sum") #write in txt
#WinRhizo wrong
all_D$d_V_WR<-(all_D$WR_RootVolume-all_D$vol)/all_D$vol
summary(abs(all_D$d_V_WR)*100)
summary((all_D$d_V_WR)*100)
rel_diff(all_D$WR_RootVolume,all_D$vol,"WR_RootVolume") #write in txt
#IJ
all_D$d_V_IJ<-(all_D$IJ_Eq_Vol-all_D$vol)/all_D$vol
summary(abs(all_D$d_V_IJ)*100)
summary((all_D$d_V_IJ)*100)
rel_diff(all_D$IJ_Eq_Vol,all_D$vol,"IJ_Eq_Vol") #write in txt


#### Analyse Relationships to ground truth ####
# plots relationships to ground truth#

outp_models_start<-function(Vari,x) { #with columnames
  model <-Vari #model
  int <- as.numeric(summary(x)$coefficients[1]) #intercept
  stderr_int <- as.numeric(summary(x)$coefficients[1,2]) #standarderror
  p_int<- as.numeric(summary(x)$coefficients[1,4]) #p_int
  est <- as.numeric(summary(x)$coefficients[2]) #estimates
  stderr_est <- as.numeric(summary(x)$coefficients[2,2]) #standarderror
  p_est <- as.numeric(summary(x)$coefficients[2,4]) #standarderror
  r2<- as.numeric(summary(x)$adj.r.squared) #adjR2
  
  sink(file = "400models.txt", append = T)
  d<-data.frame(model=model,int =int, stderr_int=stderr_int, p_int=p_int, est=est, stderr_est=stderr_est, p_est=p_est, r2=r2)
  print(d)
  sink()
}
outp_models<-function(Vari,x) { #without columnames
model <-Vari #model
int <- as.numeric(summary(x)$coefficients[1]) #intercept
stderr_int <- as.numeric(summary(x)$coefficients[1,2]) #standarderror
p_int<- as.numeric(summary(x)$coefficients[1,4]) #p_int
est <- as.numeric(summary(x)$coefficients[2]) #estimates
stderr_est <- as.numeric(summary(x)$coefficients[2,2]) #standarderror
p_est <- as.numeric(summary(x)$coefficients[2,4]) #standarderror
r2<- as.numeric(summary(x)$adj.r.squared) #adjR2

sink(file = "400models.txt", append = T)
d<-data.frame(model=model,int =int, stderr_int=stderr_int, p_int=p_int, est=est, stderr_est=stderr_est, p_est=p_est, r2=r2)
names(d)<-NULL
print(d)
sink()
}

# diameter WinRhizo average
plot(all_D$diam,all_D$WR_AvgDiam)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$WR_AvgDiam~all_D$diam)
summary(r1)

outp_models_start("WR_AvgDiam",r1)

# diameter WinRhizo weighted mean
plot(all_D$diam,all_D$WR_Weighted_mean_diam)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$WR_Weighted_mean_diam~all_D$diam)
summary(r1)

outp_models("WR_Weighted_mean_diam",r1)

# diameter IJ
plot(all_D$diam,all_D$IJ_Mean_Dia)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$IJ_Mean_Dia~all_D$diam)
summary(r1)
outp_models("IJ_Mean_Dia",r1)

# diameter IJC
plot(all_D$diam,all_D$IJ_Mean_Dia_C)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$IJ_Mean_Dia_C~all_D$diam)
summary(r1)
outp_models("IJ_Mean_Dia_C",r1)

# length WinRhizo
plot(all_D$Sum_length,all_D$WR_Length_Sum)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$WR_Length_Sum~all_D$Sum_length)
summary(r1)
outp_models("WR_Length_Sum",r1)

# length IJ
plot(all_D$Sum_length,all_D$IJ_Kimura_length)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$IJ_Kimura_length~all_D$Sum_length)
summary(r1)
outp_models("IJ_Kimura_length",r1)

# volume WinRhizo wrong
plot(all_D$vol,all_D$WR_RootVolume)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$WR_RootVolume~all_D$vol) 
abline(lm(all_D$WR_RootVolume~all_D$vol))
summary(r1)
outp_models("WR_RootVolume",r1)

# volume WinRhizo correct
plot(all_D$vol,all_D$WR_V_sum)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$WR_V_sum~all_D$vol) 
summary(r1)
outp_models("WR_V_sum",r1)

# volume IJ
plot(all_D$vol,all_D$IJ_Eq_Vol)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$IJ_Eq_Vol~all_D$vol) 
summary(r1)
outp_models("IJ_Eq_Vol",r1)

# volume IJ vs. volume WinRhizo wrong
plot(all_D$IJ_Eq_Vol,all_D$WR_RootVolume)
abline(0,1, col="grey",lty="dashed")
r1<-lm(all_D$WR_RootVolume~all_D$IJ_Eq_Vol)
summary(r1)