architect<-function(inputrac=NULL, inputtps=NULL, inputrsml=NULL, res=NULL, unitlength="px", rsml.date=NULL, rsml.connect=FALSE, vertical3d="y", fitter=FALSE){
# Find out what kind of data input we have
if (is.null(inputrac)==FALSE){
if ("character" %in% class(inputrac) | "dartToTable" %in% class(inputrac)){} else {stop("inputrac must be a character string or a dartToTable object")}}
if (is.null(inputrsml)==FALSE){
if ("character" %in% class(inputrsml) | "rsmlToTable" %in% class(inputrsml)){} else {stop("inputrsml must be a character string or a rsmlToTable object")}}
if (is.null(inputrac)==TRUE & is.null(inputrsml)==TRUE){algo<-"rawfiles"}
if (is.null(inputrac)==TRUE & "rsmlToTable" %in% class(inputrsml)){algo<-"tables"}
if (is.null(inputrac)==TRUE & mode(inputrsml)=="character"){algo<-"rawfiles"}
if ("dartToTable" %in% class(inputrac) & is.null(inputrsml)==TRUE){algo<-"tables"}
if ("dartToTable" %in% class(inputrac) & "rsmlToTable" %in% class(inputrsml)){algo<-"tables"}
if ("dartToTable" %in% class(inputrac) & mode(inputrsml)=="character"){stop("inputrsml must be a rsmlToTable object")}
if (mode(inputrac)=="character" & is.null(inputrsml)==TRUE){algo<-"rawfiles"}
if (mode(inputrac)=="character" & "rsmlToTable" %in% class(inputrsml)){stop("inputrac must be a dartToTable object")}
if (mode(inputrac)=="character" & mode(inputrsml)=="character"){algo<-"rawfiles"}
##Two possible algorithms for architect
##rawfiles if raw DART and RSML files have to be loaded (the raw files are stored in a folder)
##tables if input data are stored in rsmlToTable or dartToTable objects
#################
#First algorithm
#################
if (algo=="rawfiles"){
# Errors interception
if (is.null(inputrac)==TRUE & is.null(inputtps)==TRUE & is.null(inputrsml)==TRUE){stop("inputrac/inputps and/or inputrsml must be provided")}
if (is.null(inputrac)==FALSE) {if (mode(inputrac)!="character"){stop("mode(inputrac) must be character")}}
if (is.null(inputtps)==FALSE) {if (mode(inputtps)!="character"){stop("mode(inputtps) must be character")}}
if (is.null(inputrac)==FALSE|is.null(inputtps)==FALSE){
if (is.null(inputrac)==TRUE|is.null(inputtps)==TRUE){stop("If inputrac/inputtps is not NULL, inputtps/inputrac must be provided")}}
if (is.null(inputrsml)==FALSE) {if (mode(inputrsml)!="character"){stop("mode(inputrsml) must be character")}}
if (is.null(inputrac)==FALSE & is.null(res)==TRUE & unitlength!="px"){stop("If unitlength is not px, res must be specified")}
if (is.null(res)==FALSE){
if (mode(res)!="numeric"){stop("mode(res) must be numeric")}
if (res<=0|length(res)>1){stop("res must be a single positive value")}}
if (mode(unitlength)!="character"){stop("mode(unitlength) must be character")}
if (unitlength=="px"|unitlength=="mm"|unitlength=="cm") {} else {stop("unitlength must be either px (pixels), mm (millimeters) or cm (centimeters)")}
if (is.null(rsml.date)==FALSE){
if (is.character(rsml.date)==TRUE|is.numeric(rsml.date)==TRUE){} else {stop("If rsml.date is not NULL, rsml.date must be a character string or a positive numeric value")}
if (is.numeric(rsml.date)==TRUE){if (rsml.date<=0|length(rsml.date)>1){stop("If mode(rsml.date) is numeric, rsml.date must be a single positive value")}}}
if (mode(rsml.connect)!="logical"){stop("mode(rsml.connect) must be logical")}
# Reading of DART and rsml files
if (is.null(inputtps)==FALSE){
filenames.tps<-mixedsort(list.files(path=inputtps, pattern="\\.tps$"))
path.tps<-rep(inputtps, length.out=length(filenames.tps))
filenamestps<-sub(x=filenames.tps, pattern="\\.tps$", replacement="")
message(paste("Number of DART tps files in inputtps:", length(filenames.tps), sep=" "))}
if (is.null(inputrac)==FALSE){
filenames.rac<-mixedsort(list.files(path=inputrac, pattern="\\.rac$"))
path.rac<-rep(inputrac, length.out=length(filenames.rac))
filenamesrac<-sub(x=filenames.rac, pattern="\\.rac$", replacement="")
message(paste("Number of DART rac files in inputrac:", length(filenames.rac), sep=" "))}
if (is.null(inputrsml)==FALSE) {
filenames.rsml<-mixedsort(list.files(path=inputrsml, pattern="\\.rsml$"))
path.rsml<-rep(inputrsml, length.out=length(filenames.rsml))
filenamesrsml<-sub(x=filenames.rsml, pattern="\\.rsml$", replacement="")
message(paste("Number of rsml files in inputrsml:", length(filenames.rsml), sep=" "))}
if (is.null(inputrsml)==TRUE){
if (length(filenames.rac)==0){stop("There is no rac file in inputrac")}
if (length(filenames.tps)==0){stop("There is no tps file in inputtps")}}
else {
if (is.null(inputrac)==TRUE){if (length(filenames.rsml)==0){stop("There is no rsml file in inputrsml")}}
else{
if (length(filenames.rac)==0){stop("There is no rac file in inputrac")}
if (length(filenames.tps)==0){stop("There is no tps file in inputtps")}
if (length(filenames.rsml)==0){stop("There is no rsml file in inputrsml")}}}
if (is.null(inputrsml)==TRUE){ # Only DART files
TIME<-lapply(paste(path.tps, "/", filenames.tps, sep=""), read.table, header=TRUE)
DATA<-lapply(paste(path.rac, "/", filenames.rac, sep=""), read.table, skip=1)
for (i in 1:length(DATA)) {
colnames(DATA[[i]])<-c()
colnames(DATA[[i]])[1]<-"Root"
colnames(DATA[[i]])[2]<-"Mother"
colnames(DATA[[i]])[3]<-"Ord"
colnames(DATA[[i]])[4]<-"DBase"
colnames(DATA[[i]])[5]<-"DApp"
for (j in 6:ncol(DATA[[i]])-5) {colnames(DATA[[i]])[j+5]<-paste("Lengths", j, sep="")}}
#Unit conversion DART files
if (unitlength=="mm") {cunit<-(10*cm(1)/res)}
if (unitlength=="cm") {cunit<-(cm(1)/res)}
if (unitlength=="px") {cunit<-1}
if (length(TIME)==1) {
age<-list()
obstot<-0
for (i in 1:length(DATA)) {
age[[i]]<-TIME[[1]]$Date
obstot<-obstot+length(age[[i]])}
for (i in 1:length(DATA)) {if(length(age[[i]])!=(ncol(DATA[[i]])-5)){stop("The number of observation dates between corresponding rac et tps files must be equal")}}}
else {
if (length(TIME)!=length(DATA)) {stop("If there is more than one tps file in inputtps, the number of rac files in inputrac and tps files in inputtps must be equal")}
else {
for (i in 1:length(DATA)) {if (filenamesrac[i]!=filenamestps[i]) {stop("Input rac files and their corresponding tps files must have the same name")}}
age<-list()
obstot<-0
for (i in 1:length(TIME)) {
age[[i]]<-TIME[[i]]$Date
obstot<-obstot+length(age[[i]])}
for (i in 1:length(DATA)) {if (length(age[[i]])!=(ncol(DATA[[i]])-5)) {stop("The number of observation dates between corresponding rac et tps files must be equal")}}}}}
else {
if (is.null(inputrac)==TRUE){ # Only RSML files
DATA<-list()
age<-list()
res1<-c()
unitlength1<-c()
filenamesrac<-c()
RSML <- lapply(paste(path.rsml, "/", filenames.rsml, sep=""), rsmlToDART, final.date=rsml.date, connect=rsml.connect) # Read RSML files
obstot<-0
for (i in 1:length(RSML)){
res1<-append(res1, rep(as.numeric(RSML[[i]]$resolution), length(RSML[[i]]$lie)))
unitlength1<-append(unitlength1, rep(as.character(RSML[[i]]$length), length(RSML[[i]]$lie)))
DATA<-append(DATA, RSML[[i]]$rac)
length1<-length(RSML[[i]]$rac)
l<-length(age)
for (j in 1:length1) {
age[[l+j]]<-RSML[[i]]$tps[[j]]$Date
obstot<-obstot+length(age[[l+j]])}
if (length1>1){
num<-c(1:length1)
filenamesrac[(length(filenamesrac)+1):(length(filenamesrac)+length1)]<-paste(rep(filenamesrsml[i], length.out=length1), num, sep="")}
if (length1==1){
filenamesrac[(length(filenamesrac)+1)]<-filenamesrsml[i]}}
#Unit conversion RSML
cunit1<-vector(length=length(res1))
for (i in 1:length(res1)){
if (unitlength=="cm"){
if (unitlength1[i]=="pixel") {
cunit1[i]<-1
message(paste("Unit in ", filenamesrac[i], " is pixel. Unitlength not used and results expressed in pixels", sep=""))}
if (unitlength1[i]=="m") {cunit1[i]<-100/res1[i]}
if (unitlength1[i]=="cm") {cunit1[i]<-1/res1[i]}
if (unitlength1[i]=="mm") {cunit1[i]<-1/res1[i]/10}
if (unitlength1[i]=="um") {cunit1[i]<-1/res1[i]/10000}
if (unitlength1[i]=="nm") {cunit1[i]<-1/res1[i]/10000000}
if (unitlength1[i]=="inch") {cunit1[i]<-1/res1[i]*cm(1)}}
if (unitlength=="mm"){
if (unitlength1[i]=="pixel") {
cunit1[i]<-1
message(paste("Unit in ", filenamesrac[i], " is pixel. Unitlength not used and results expressed in pixels", sep=""))}
if (unitlength1[i]=="m") {cunit1[i]<-1/res1[i]*1000}
if (unitlength1[i]=="cm") {cunit1[i]<-1/res1[i]*10}
if (unitlength1[i]=="mm") {cunit1[i]<-1/res1[i]}
if (unitlength1[i]=="um") {cunit1[i]<-1/res1[i]/1000}
if (unitlength1[i]=="nm") {cunit1[i]<-1/res1[i]/1000000}
if (unitlength1[i]=="inch") {cunit1[i]<-1/res1[i]*cm(1)*10}}
if (unitlength=="px"){cunit1[i]<-1}}}
else { # DART and RSML files
TIME<-lapply(paste(path.tps, "/", filenames.tps, sep=""), read.table, header=TRUE) # Read TPS files
DATA<-lapply(paste(path.rac, "/", filenames.rac, sep=""), read.table, skip=1) # Read RAC files
for (i in 1:length(DATA)) {
colnames(DATA[[i]])<-c()
colnames(DATA[[i]])[1]<-"Root"
colnames(DATA[[i]])[2]<-"Mother"
colnames(DATA[[i]])[3]<-"Ord"
colnames(DATA[[i]])[4]<-"DBase"
colnames(DATA[[i]])[5]<-"DApp"
for (j in 6:ncol(DATA[[i]])-5) {colnames(DATA[[i]])[j+5]<-paste("Lengths", j, sep="")}}
res1<-rep(res, length(filenames.rac))
unitlength1<-rep(unitlength, length(filenames.rac))
if (length(TIME)==1) {
age<-list()
obstot<-0
for (i in 1:length(DATA)) {
age[[i]]<-TIME[[1]]$Date
obstot<-obstot+length(age[[i]])}
for (i in 1:length(DATA)) {if (length(age[[i]])!=(ncol(DATA[[i]])-5)) {stop("The number of observation dates between corresponding rac et tps files must be equal")}}}
else {
if (length(TIME)!=length(DATA)) {stop("If there is more than one tps file in inputtps, the number of rac files in inputrac and tps files in inputtps must be equal")}
else {
for (i in 1:length(DATA)) {if (filenamesrac[i]!=filenamestps[i]) {stop("Input rac files and their corresponding tps files must have the same name")}}
age<-list()
obstot<-0
for (i in 1:length(DATA)) {
age[[i]]<-TIME[[i]]$Date
obstot<-obstot+length(age[[i]])}
for (i in 1:length(DATA)) {if (length(age[[i]])!=(ncol(DATA[[i]])-5)) {stop("The number of observation dates between corresponding rac et tps files must be equal")}}}}
RSML <- lapply(paste(path.rsml, "/", filenames.rsml, sep=""), rsmlToDART, final.date=rsml.date, connect=rsml.connect) # Read RSML files
for (i in 1:length(RSML)){
res1<-append(res1, rep(as.numeric(RSML[[i]]$resolution), length(RSML[[i]]$lie)))
unitlength1<-append(unitlength1, rep(as.character(RSML[[i]]$length), length(RSML[[i]]$lie)))
DATA<-append(DATA, RSML[[i]]$rac)
length1<-length(RSML[[i]]$rac)
l<-length(age)
for (j in 1:length1) {
age[[l+j]]<-RSML[[i]]$tps[[j]]$Date
obstot<-obstot+length(age[[l+j]])}
if (length1>1){
num<-c(1:length1)
filenamesrac[(length(filenamesrac)+1):(length(filenamesrac)+length1)]<-paste(rep(filenamesrsml[i], length.out=length1), num, sep="")}
if (length1==1){
filenamesrac[(length(filenamesrac)+1)]<-filenamesrsml[i]}}
#Unit conversion DART and RSML
cunit1<-vector(length=length(res1))
if (unitlength=="mm") {cunit1[1:length(filenames.rac)]<-(10*cm(1)/res)}
if (unitlength=="cm") {cunit1[1:length(filenames.rac)]<-cm(1)/res}
if (unitlength=="px") {cunit1[1:length(filenames.rac)]<-1}
for (i in (length(filenames.rac)+1):length(res1)){
if (unitlength=="cm"){
if (unitlength1[i]=="pixel") {
cunit[i]<-1
message(paste("Unit in ", filenamesrac[i], " is pixel. Unitlength not used and results expressed in pixels", sep=""))}
if (unitlength1[i]=="m") {cunit1[i]<-100/res1[i]}
if (unitlength1[i]=="cm") {cunit1[i]<-1/res1[i]}
if (unitlength1[i]=="mm") {cunit1[i]<-1/res1[i]/10}
if (unitlength1[i]=="um") {cunit1[i]<-1/res1[i]/10000}
if (unitlength1[i]=="nm") {cunit1[i]<-1/res1[i]/10000000}
if (unitlength1[i]=="inch") {cunit1[i]<-1/res1[i]*cm(1)}}
if (unitlength=="mm"){
if (unitlength1[i]=="pixel") {
cunit[i]<-1
message(paste("Unit in ", filenamesrac[i], " is pixel. Unitlength not used and results expressed in pixels", sep=""))}
if (unitlength1[i]=="m") {cunit1[i]<-1/res1[i]*1000}
if (unitlength1[i]=="cm") {cunit1[i]<-1/res1[i]*10}
if (unitlength1[i]=="mm") {cunit1[i]<-1/res1[i]}
if (unitlength1[i]=="um") {cunit1[i]<-1/res1[i]/1000}
if (unitlength1[i]=="nm") {cunit1[i]<-1/res1[i]/1000000}
if (unitlength1[i]=="inch") {cunit1[i]<-1/res1[i]*cm(1)*10}}
if (unitlength=="px"){cunit1[i]<-1}}}}
# Creating vectors and matrices for root architecture parameters calculation
FileNames<-c()
Time<-c()
FirstOrderRootLength<-c()
FirstOrderRootNumber<-c()
TotalRootLength<-c()
TotalLateralRootNumber<-c()
TotalLateralRootLength<-c()
LateralRootDensity<-c()
GrowthRateFirstOrderRoot<-c()
GrowthRateTotalRoot<-c()
# Calculation of root architecture parameters
k<-0
maxord<-max(sapply(DATA, function(x){max(x[[3]])}))
if (maxord>1) {latroot<-matrix(ncol=4*(maxord-1), nrow=obstot)}
for(i in 1:length(DATA)){
LastFirstOrderRootLength<-DATA[[i]][[ncol(DATA[[i]])]][DATA[[i]]$Ord==1]
for (t in 1:(ncol(DATA[[i]])-5)){
k<-k+1
# Time
Time[k]<-age[[i]][t]
#Split data per root order
data.split<-split(DATA[[i]][[paste("Lengths",t,sep="")]],as.factor(DATA[[i]]$Ord))
# File names
FileNames[k]<-filenamesrac[i]
# Unit conversion
if (is.null(inputrsml)==FALSE){cunit<-cunit1[i]}
# Total root length
TotalRootLength[k]<-sum(DATA[[i]][[paste("Lengths",t,sep="")]])*cunit
# Growth rate of the root system
if (t==1) {GrowthRateTotalRoot[k]<-TotalRootLength[k]/age[[i]][t]}
if (t>1) {GrowthRateTotalRoot[k]<-(TotalRootLength[k]-TotalRootLength[k-1])/(age[[i]][t]-age[[i]][t-1])}
# First-order root length
FirstOrderRootLength[k]<-sum(data.split$'1')*cunit
# Growth rate of the first-order root
if (t==1) {GrowthRateFirstOrderRoot[k]<-FirstOrderRootLength[k]/age[[i]][t]}
if (t>1) {GrowthRateFirstOrderRoot[k]<-(FirstOrderRootLength[k]-FirstOrderRootLength[k-1])/(age[[i]][t]-age[[i]][t-1])}
# Total number of first-order roots
FirstOrderRootNumber[k]<-sum(data.split$'1'!=0)
# Total number of lateral roots
if (FirstOrderRootLength[k]==0) {TotalLateralRootNumber[k]<-0} else {TotalLateralRootNumber[k]<-sum(DATA[[i]][[paste("Lengths",t,sep="")]]!=0)-FirstOrderRootNumber[k]}
# Total length of lateral roots
TotalLateralRootLength[k]<-TotalRootLength[k]-FirstOrderRootLength[k]
# Number, length and growth rate of lateral roots (by branching order)
if (maxord>1){
for (l in 1:(maxord-1)){
if (l<=max(DATA[[i]]$Ord)-1) {latroot[k,l]<-sum(data.split[[l+1]]!=0)} else {latroot[k,l]<-0}
if (l<=max(DATA[[i]]$Ord)-1) {latroot[k,(l+(maxord-1))]<-sum(data.split[[l+1]])*cunit} else {latroot[k,(l+(maxord-1))]<-0}
if (latroot[k,l]==0){latroot[k,(l+2*(maxord-1))]<-0} else {latroot[k,(l+2*(maxord-1))]<-latroot[k,(l+(maxord-1))]/latroot[k,l]}
if (t==1) {latroot[k,(l+3*(maxord-1))]<-latroot[k,(l+(maxord-1))]/age[[i]][t]} else {latroot[k,(l+3*(maxord-1))]<-(latroot[k,(l+(maxord-1))]-latroot[k-1,(l+(maxord-1))])/(age[[i]][t]-age[[i]][t-1])}}}
# Density of secondary roots on the first-order root
if (FirstOrderRootLength[k]==0|maxord==1) {LateralRootDensity[k]<-0} else {LateralRootDensity[k]<-latroot[k,1]/FirstOrderRootLength[k]}}}
# Summary results in a data frame
if (maxord>1) {outputresults<-data.frame(FileNames, Time, TotalRootLength, GrowthRateTotalRoot, FirstOrderRootLength, GrowthRateFirstOrderRoot, FirstOrderRootNumber, TotalLateralRootNumber, TotalLateralRootLength, latroot, LateralRootDensity)}
else {outputresults<-data.frame(FileNames, Time, TotalRootLength, GrowthRateTotalRoot, FirstOrderRootLength, GrowthRateFirstOrderRoot, FirstOrderRootNumber, TotalLateralRootNumber, TotalLateralRootLength, LateralRootDensity)}
if (maxord>1){
LRnumberheading<-c()
LRlengthheading<-c()
LRmeanlengthheading<-c()
LRgrowthrateheading<-c()
for (h in 2:maxord){
LRnumberheading[h-1]<-paste("N", h, "LR", sep="")
LRlengthheading[h-1]<-paste("L", h, "LR", sep="")
LRmeanlengthheading[h-1]<-paste("ML", h, "LR", sep="")
LRgrowthrateheading[h-1]<-paste("GR", h, "L", sep="")}}
if (maxord>1) {colnames(outputresults)<-c("FileName", "Time", "TRL", "GRTR", "L1R", "GR1R", "TN1R", "TNLR", "TLRL", t(LRnumberheading), t(LRlengthheading), t(LRmeanlengthheading), t(LRgrowthrateheading), "D2LR")}
else {colnames(outputresults)<-c("FileName", "Time", "TRL", "GRTR", "L1R", "GR1R", "TN1R", "TNLR", "TLRL", "D2LR")}
return(outputresults)}
##################
#Second algorithm
##################
if (algo=="tables"){
if (is.null(inputrsml)==FALSE & fitter==TRUE){
if ("geodesic" %in% colnames(inputrsml)) {} else {stop("No connection found between parent and daughter roots. Consider using rsml.connect=TRUE in rsmlToTable.")}}
if (vertical3d=="x"|vertical3d=="y"|vertical3d=="z") {} else {stop("vertical3d must be x, y, or z")}
if (mode(fitter)!="logical"){stop("fitter must be logical")}
if (is.null(inputrac)==TRUE & is.null(inputrsml)==FALSE) {maxord<-max(inputrsml$order)}
if (is.null(inputrac)==FALSE & is.null(inputrsml)==TRUE) {maxord<-max(inputrac$order)}
if (is.null(inputrac)==FALSE & is.null(inputrsml)==FALSE) {maxord<-max(c(inputrac$order, inputrsml$order))}
#Processing DART files
if (is.null(inputrac)==FALSE){
n<-length(unique(paste(inputrac$file, inputrac$time, sep="_")))
if (fitter==TRUE){datadart<-data.frame(FileName=rep(NA, n), Time=rep(NA, n), TRL=rep(NA, n), GRTR=rep(NA, n), L1R=rep(NA, n), GR1R=rep(NA, n), TN1R=rep(NA, n), TNLR=rep(NA, n), TLRL=rep(NA, n), D2LR=rep(NA, n), Height=rep(NA, n), Width=rep(NA, n), ConvexhullXY=rep(NA, n), Magnitude=rep(NA, n), Altitude=rep(NA, n), ExtPathLength=rep(NA, n))}
else {datadart<-data.frame(FileName=rep(NA, n), Time=rep(NA, n), TRL=rep(NA, n), GRTR=rep(NA, n), L1R=rep(NA, n), GR1R=rep(NA, n), TN1R=rep(NA, n), TNLR=rep(NA, n), TLRL=rep(NA, n), D2LR=rep(NA, n), Height=rep(NA, n), Width=rep(NA, n), ConvexhullXY=rep(NA, n))}
if (maxord>1){latroot<-matrix(ncol=4*(maxord-1), nrow=n)}
diameter<-matrix(ncol=maxord+1, nrow=n)
surface<-matrix(ncol=maxord+1, nrow=n)
volume<-matrix(ncol=maxord+1, nrow=n)
files<-unique(inputrac$file)
k<-0
for (i in 1:length(files)){
dates<-sort(unique(inputrac$time[inputrac$file==files[i]]))
for (t in 1:length(dates)){
xt<-inputrac[inputrac$file==files[i],]
k<-k+1
#File name
datadart$FileName[k]<-xt$file[1]
#Time
datadart$Time[k]<-dates[t]
#TRL (total root length)
datadart$TRL[k]<-sum(xt$length[xt$time<=dates[t]])
#GRTR (growth rate of the root system)
GR<-aggregate(xt$growth, by=list(root=xt$root, order=xt$order, time=xt$time), max)
datadart$GRTR[k]<-sum(GR$x[GR$time==dates[t]])
#L1R (total first-order root length)
datadart$L1R[k]<-sum(xt$length[xt$time<=dates[t] & xt$order==1])
#GR1R (growth rate of first-order roots)
datadart$GR1R[k]<-sum(GR$x[GR$time==dates[t] & GR$order==1])
#TN1R (number of first-order roots)
datadart$TN1R[k]<-length(unique(xt$root[xt$time<=dates[t] & xt$order==1]))
#TNLR (number of lateral roots)
datadart$TNLR[k]<-length(unique(xt$root[xt$time<=dates[t] & xt$order>1]))
#TLRL (total lateral root length)
datadart$TLRL[k]<-sum(xt$length[xt$time<=dates[t] & xt$order>1])
#D2LR (density of second-order roots on the first-order root)
datadart$D2LR[k]<-length(unique(xt$root[xt$time<=dates[t] & xt$order==2]))/datadart$L1R[k]
#Latroot (number, length, mean length, growth rate)
if (maxord>1){
for (l in 1:(maxord-1)){
latroot[k, l]<-length(unique(xt$root[xt$time<=dates[t] & xt$order==(l+1)])) #Number of lateral roots
latroot[k, l+(maxord-1)]<-sum(xt$length[xt$time<=dates[t] & xt$order==(l+1)]) #Length of lateral roots
if (latroot[k,l]==0) {latroot[k, l+2*(maxord-1)]<-0} else {latroot[k, l+2*(maxord-1)]<-latroot[k, l+(maxord-1)]/latroot[k, l]} #Mean lateral root length
latroot[k, l+3*(maxord-1)]<-sum(GR$x[GR$time==dates[t] & GR$order==(l+1)])}}
#Height
datadart$Height[k]<-abs(max(xt$y2[xt$time<=dates[t]])-min(xt$y1[xt$time<=dates[t]]))
#Width
datadart$Width[k]<-abs(max(xt$x2[xt$time<=dates[t]])-min(xt$x2[xt$time<=dates[t]]))
#Convex hull
x<-c(xt$x1[xt$time<=dates[t]], xt$x2[xt$time<=dates[t]])
y<-c(xt$y1[xt$time<=dates[t]], xt$y2[xt$time<=dates[t]])
root.coords<-cbind(x,y)
boxhull<-chull(x=x, y=y)
if (length(boxhull)<3){datadart$ConvexhullXY[k]<-0}
else {
boxhull<-c(boxhull, boxhull[1])
boxhull.coords<-root.coords[boxhull,]
chull.poly<-Polygon(boxhull.coords, hole=F)
datadart$ConvexhullXY[k]<-chull.poly@area}
#Fitter topological indices
if (fitter==TRUE){
xt<-xt[xt$time<=dates[t],]
apex<-aggregate(xt$blength, by=list(root=xt$root), max)
apicindex<-as.vector(apply(apex, 1, function(x){which(xt$root==x["root"] & xt$blength==x["x"])}))
xt$apic[apicindex]<-"true"
branindex<-which(xt$bran=="true")
#Magnitude
datadart$Magnitude[k]<-sum(xt$bran=="true")
#Altitude and external path length
xt$pathlength<-rep(1, nrow(xt))
if (nrow(xt)>1){
for (l in 1:length(branindex)){ #For each root
root<-xt$root[branindex[l]]
testbran<-which(xt$x1==xt$x2[branindex[l]] & xt$y1==xt$y2[branindex[l]] & xt$z1==xt$z2[branindex[l]]) #Is it a crossing?
if (length(testbran)>1) {testbran<-testbran[which(xt$root[testbran]==root | xt$parentroot[testbran]==root)]} #Select segments based on root ID and parentroot ID
if (length(testbran)==0) {} else {
if (length(testbran)>=2) {
xt$pathlength[testbran]<-xt$pathlength[branindex[l]]+1
index<-which(xt$bran[testbran]=="false")
suiv<-testbran[index]}
else {
xt$pathlength[testbran]<-xt$pathlength[branindex[l]]
suiv<-testbran}
while(xt$apic[suiv]=="false"){
testbran<-which(xt$x1==xt$x2[suiv] & xt$y1==xt$y2[suiv] & xt$z1==xt$z2[suiv]) #Is it a crossing?
if (length(testbran)>1) {testbran<-testbran[which(xt$root[testbran]==root | xt$parentroot[testbran]==root)]} #Select segments based on root ID and parentroot ID
if (length(testbran)>=2) {
xt$pathlength[testbran]<-xt$pathlength[suiv]+1
index<-which(xt$bran[testbran]=="false")
suiv<-testbran[index]}
else {
xt$pathlength[testbran]<-xt$pathlength[suiv]
suiv<-testbran}}}}
datadart$Altitude[k]<-max(xt$pathlength)
datadart$ExtPathLength[k]<-sum(xt$pathlength[xt$apic=="true"])}
else{
datadart$Altitude[k]<-1
datadart$ExtPathLength[k]<-1}}}}
# Results in a dataframe
diameter<-as.data.frame(diameter)
colnames(diameter)<-c(paste(rep("MD", maxord), c(1:maxord), sep=""), "MDLR")
surface<-as.data.frame(surface)
colnames(surface)<-c(paste(rep("S", maxord), c(1:maxord), sep=""), "Stot")
volume<-as.data.frame(volume)
colnames(volume)<-c(paste(rep("V", maxord), c(1:maxord), sep=""), "Vtot")
if (maxord>1){
latroot<-as.data.frame(latroot)
for (l in 1:(maxord-1)){
colnames(latroot)[l]<-paste("N", l+1, "LR", sep="")
colnames(latroot)[l+(maxord-1)]<-paste("L", l+1, "LR", sep="")
colnames(latroot)[l+2*(maxord-1)]<-paste("ML", l+1, "LR", sep="")
colnames(latroot)[l+3*(maxord-1)]<-paste("GR", l+1, "L", sep="")}
if (fitter==TRUE){datadart<-data.frame(datadart[,1:9], as.data.frame(latroot), datadart[,10:16], as.data.frame(diameter), as.data.frame(surface), as.data.frame(volume))}
else {datadart<-data.frame(datadart[,1:9], as.data.frame(latroot), datadart[,10:13], as.data.frame(diameter), as.data.frame(surface), as.data.frame(volume))}}}
#Processing RSML files
if (is.null(inputrsml)==FALSE){
inputrsml$diameter3<-inputrsml$diameter1
inputrsml$diameter3[inputrsml$apic=="true"]<-inputrsml$diameter1[inputrsml$apic=="true"]+inputrsml$diameter2[inputrsml$apic=="true"]
n<-length(unique(paste(inputrsml$file, inputrsml$plant, inputrsml$time, sep="_")))
if (fitter==TRUE){datarsml<-data.frame(FileName=rep(NA, n), Time=rep(NA, n), TRL=rep(NA, n), GRTR=rep(NA, n), L1R=rep(NA, n), GR1R=rep(NA, n), TN1R=rep(NA, n), TNLR=rep(NA, n), TLRL=rep(NA, n), D2LR=rep(NA, n), Height=rep(NA, n), Width=rep(NA, n), ConvexhullXY=rep(NA, n), ConvexhullXZ=rep(NA, n), ConvexhullYZ=rep(NA, n), Convexhull3D=rep(NA, n), Magnitude=rep(NA, n), Altitude=rep(NA, n), ExtPathLength=rep(NA, n))}
else {datarsml<-data.frame(FileName=rep(NA, n), Time=rep(NA, n), TRL=rep(NA, n), GRTR=rep(NA, n), L1R=rep(NA, n), GR1R=rep(NA, n), TN1R=rep(NA, n), TNLR=rep(NA, n), TLRL=rep(NA, n), D2LR=rep(NA, n), Height=rep(NA, n), Width=rep(NA, n), ConvexhullXY=rep(NA, n), ConvexhullXZ=rep(NA, n), ConvexhullYZ=rep(NA, n), Convexhull3D=rep(NA, n))}
if (maxord>1){latroot<-matrix(ncol=4*(maxord-1), nrow=n)}
diameter<-matrix(ncol=maxord+1, nrow=n)
surface<-matrix(ncol=maxord+1, nrow=n)
volume<-matrix(ncol=maxord+1, nrow=n)
files<-unique(inputrsml$file)
k<-0
for (i in 1:length(files)){
plants<-unique(inputrsml$plant[inputrsml$file==files[i]])
for (j in 1:length(plants)){
dates<-sort(unique(inputrsml$time[inputrsml$file==files[i] & inputrsml$plant==plants[j]]))
for (t in 1:length(dates)){
xt<-inputrsml[inputrsml$file==files[i] & inputrsml$plant==plants[j],] #x is a subset of the rsmlToTable object
k<-k+1
#File name
datarsml$FileName[k]<-paste(xt$file[1], xt$plant[1], sep="_")
#Time
datarsml$Time[k]<-dates[t]
#TRL (total root length)
datarsml$TRL[k]<-sum(xt$length[xt$time<=dates[t]])
#GRTR (growth rate of the root system)
GR<-aggregate(xt$growth, by=list(root=xt$root, order=xt$order, time=xt$time), max)
datarsml$GRTR[k]<-sum(GR$x[GR$time==dates[t]])
#L1R (total first-order root length)
datarsml$L1R[k]<-sum(xt$length[xt$time<=dates[t] & xt$order==1])
#GR1R (growth rate of first-order roots)
datarsml$GR1R[k]<-sum(GR$x[GR$time==dates[t] & GR$order==1])
#TN1R (number of first-order roots)
datarsml$TN1R[k]<-length(unique(xt$root[xt$time<=dates[t] & xt$order==1]))
#TNLR (number of lateral roots)
datarsml$TNLR[k]<-length(unique(xt$root[xt$time<=dates[t] & xt$order>1]))
#TLRL (total lateral root length)
datarsml$TLRL[k]<-sum(xt$length[xt$time<=dates[t] & xt$order>1])
#D2LR (density of second-order roots on the first-order root)
datarsml$D2LR[k]<-length(unique(xt$root[xt$time<=dates[t] & xt$order==2]))/datarsml$L1R[k]
#Latroot (number, length, mean length, growth rate)
if (maxord>1){
for (l in 1:(maxord-1)){
latroot[k, l]<-length(unique(xt$root[xt$time<=dates[t] & xt$order==(l+1)])) #Number of lateral roots
latroot[k, l+(maxord-1)]<-sum(xt$length[xt$time<=dates[t] & xt$order==(l+1)]) #Length of lateral roots
if (latroot[k,l]==0) {latroot[k, l+2*(maxord-1)]<-0} else {latroot[k, l+2*(maxord-1)]<-latroot[k, l+(maxord-1)]/latroot[k, l]} #Mean lateral root length
latroot[k, l+3*(maxord-1)]<-sum(GR$x[GR$time==dates[t] & GR$order==(l+1)])}}
#Diameter, surface and volume
if (t==length(dates)){#Diameter, surface, and volume in rsml file is for the last observation date
maxordfile<-max(xt$order)
for (l in 1:maxordfile){ #Calculate for each root order
diameter[k,l]<-sum(xt$diameter3[xt$order==l])/(nrow(xt[xt$order==l,])+sum(xt$apic[xt$order==l]=="true"))
surface[k,l]<-sum(xt$surface[xt$order==l])
volume[k,l]<-sum(xt$volume[xt$order==l])}
if (maxordfile>1) {diameter[k, ncol(diameter)]<-sum(xt$diameter3[xt$order>1])/(nrow(xt[xt$order>1,])+sum(xt$apic[xt$order>1]=="true"))} else {diameter[k, ncol(diameter)]<-NA} #Mean diameter for all lateral roots
surface[k, ncol(surface)]<-sum(xt$surface)
volume[k, ncol(volume)]<-sum(xt$volume)}
#Height
if (vertical3d=="y") {datarsml$Height[k]<-abs(max(xt$y2[xt$time<=dates[t]])-min(xt$y1[xt$time<=dates[t]]))}
if (vertical3d=="x") {datarsml$Height[k]<-abs(max(xt$x2[xt$time<=dates[t]])-min(xt$x1[xt$time<=dates[t]]))}
if (vertical3d=="z") {datarsml$Height[k]<-abs(max(xt$z2[xt$time<=dates[t]])-min(xt$z1[xt$time<=dates[t]]))}
#Width
if (vertical3d=="y") {
widthx<-abs(max(xt$x2[xt$time<=dates[t]])-min(xt$x2[xt$time<=dates[t]]))
widthz<-abs(max(xt$z2[xt$time<=dates[t]])-min(xt$z2[xt$time<=dates[t]]))
datarsml$Width[k]<-max(c(widthx, widthz))}
if (vertical3d=="x") {
widthy<-abs(max(xt$y2[xt$time<=dates[t]])-min(xt$y2[xt$time<=dates[t]]))
widthz<-abs(max(xt$z2[xt$time<=dates[t]])-min(xt$z2[xt$time<=dates[t]]))
datarsml$Width[k]<-max(c(widthy, widthz))}
if (vertical3d=="z") {
widthy<-abs(max(xt$y2[xt$time<=dates[t]])-min(xt$y2[xt$time<=dates[t]]))
widthx<-abs(max(xt$x2[xt$time<=dates[t]])-min(xt$x2[xt$time<=dates[t]]))
datarsml$Width[k]<-max(c(widthy, widthx))}
#Convex hull
#For XY
x<-c(xt$x1[xt$time<=dates[t]], xt$x2[xt$time<=dates[t]])
y<-c(xt$y1[xt$time<=dates[t]], xt$y2[xt$time<=dates[t]])
root.coords<-cbind(x,y)
boxhull<-chull(x=x, y=y)
if (length(boxhull)<3) {datarsml$ConvexhullXY[k]<-0}
else {
boxhull<-c(boxhull, boxhull[1])
boxhull.coords<-root.coords[boxhull,]
chull.poly<-Polygon(boxhull.coords, hole=F)
datarsml$ConvexhullXY[k]<-chull.poly@area}
#For YZ
y<-c(xt$y1[xt$time<=dates[t]], xt$y2[xt$time<=dates[t]])
z<-c(xt$z1[xt$time<=dates[t]], xt$z2[xt$time<=dates[t]])
root.coords<-cbind(y,z)
boxhull<-chull(x=y, y=z)
if (length(boxhull)<3) {datarsml$ConvexhullYZ[k]<-0}
else{
boxhull<-c(boxhull, boxhull[1])
boxhull.coords<-root.coords[boxhull,]
chull.poly<-Polygon(boxhull.coords, hole=F)
datarsml$ConvexhullYZ[k]<-chull.poly@area}
#For XZ
x<-c(xt$x1[xt$time<=dates[t]], xt$x2[xt$time<=dates[t]])
z<-c(xt$z1[xt$time<=dates[t]], xt$z2[xt$time<=dates[t]])
root.coords<-cbind(x,z)
boxhull<-chull(x=x, y=z)
if (length(boxhull)<3) {datarsml$ConvexhullXZ[k]<-0}
else{
boxhull<-c(boxhull, boxhull[1])
boxhull.coords<-root.coords[boxhull,]
chull.poly<-Polygon(boxhull.coords, hole=F)
datarsml$ConvexhullXZ[k]<-chull.poly@area}
#Convex hull volume (3D)
if (is.na(datarsml$ConvexhullXY[k])==TRUE|
is.na(datarsml$ConvexhullXZ[k])==TRUE|
is.na(datarsml$ConvexhullYZ[k])==TRUE){datarsml$Convexhull3D[k]<-NA}
else{
if (datarsml$ConvexhullXY[k]==0|
datarsml$ConvexhullXZ[k]==0|
datarsml$ConvexhullYZ[k]==0){datarsml$Convexhull3D[k]<-0}
else{
x<-c(xt$x1[xt$time<=dates[t]], xt$x2[xt$time<=dates[t]])
y<-c(xt$y1[xt$time<=dates[t]], xt$y2[xt$time<=dates[t]])
z<-c(xt$z1[xt$time<=dates[t]], xt$z2[xt$time<=dates[t]])
root.coords<-unique(cbind(x,y,z))
if (nrow(root.coords)>=4){ #Need at least 4 points to construct initial simplex
datarsml$Convexhull3D[k]<-convhulln(as.matrix(root.coords), options=c("QJ", "FA"))$vol}
else {
datarsml$Convexhull3D[k]<-0}}}
#Fitter topological indices
if (fitter==TRUE){
xt<-xt[xt$time<=dates[t],]
apex<-aggregate(xt$blength, by=list(root=xt$root), max)
apicindex<-as.vector(apply(apex, 1, function(x){which(xt$root==x["root"] & xt$blength==x["x"])}))
xt$apic[apicindex]<-"true"
branindex<-which(xt$bran=="true") #Equal to the number of roots
#Magnitude
datarsml$Magnitude[k]<-sum(xt$bran=="true")
#Altitude and external path length
xt$pathlength<-rep(1, nrow(xt))
if (nrow(xt)>1){
for (l in 1:length(branindex)){#For each root
root<-xt$root[branindex[l]]
testbran<-which(xt$x1==xt$x2[branindex[l]] & xt$y1==xt$y2[branindex[l]] & xt$z1==xt$z2[branindex[l]]) #Is it a crossing?
if (length(testbran)>1) {testbran<-testbran[which(xt$root[testbran]==root | xt$parentroot[testbran]==root)]} #Select segments based on root ID and parentroot ID
if (length(testbran)==0) {} else {
if (length(testbran)>=2) {
xt$pathlength[testbran]<-xt$pathlength[branindex[l]]+1
index<-which(xt$bran[testbran]=="false")
if (length(index)!=0) {suiv<-testbran[index]} else {suiv<-branindex[l]}}
else {
xt$pathlength[testbran]<-xt$pathlength[branindex[l]]
suiv<-testbran}
if (xt$apic[suiv]=="true" & suiv!=branindex[l]){
testbran<-which(xt$x1==xt$x2[suiv] & xt$y1==xt$y2[suiv] & xt$z1==xt$z2[suiv]) #Is it a crossing?
if (length(testbran)>1) {testbran<-testbran[which(xt$parentroot[testbran]==root)]} #Select segments based on parentroot ID
xt$pathlength[testbran]<-xt$pathlength[suiv]+1}
while(xt$apic[suiv]=="false"){
testbran<-which(xt$x1==xt$x2[suiv] & xt$y1==xt$y2[suiv] & xt$z1==xt$z2[suiv]) #Is it a crossing?
if (length(testbran)>1) {testbran<-testbran[which(xt$root[testbran]==root | xt$parentroot[testbran]==root)]} #Select segments based on root ID and parentroot ID
if (length(testbran)>=2) {
xt$pathlength[testbran]<-xt$pathlength[suiv]+1
index<-which(xt$bran[testbran]=="false")
suiv<-testbran[index]}
else {
xt$pathlength[testbran]<-xt$pathlength[suiv]
suiv<-testbran}
testbran<-which(xt$x1==xt$x2[suiv] & xt$y1==xt$y2[suiv] & xt$z1==xt$z2[suiv])
if (xt$apic[suiv]=="true" & length(testbran)>0){xt$pathlength[testbran]<-xt$pathlength[suiv]+1}}}}
datarsml$Altitude[k]<-max(xt$pathlength)
datarsml$ExtPathLength[k]<-sum(xt$pathlength[xt$apic=="true"])}
else{
datarsml$Altitude[k]<-1
datarsml$ExtPathLength[k]<-1}}}}}
# Results in a dataframe
diameter<-as.data.frame(diameter)
colnames(diameter)<-c(paste(rep("MD", maxord), c(1:maxord), sep=""), "MDLR")
surface<-as.data.frame(surface)
colnames(surface)<-c(paste(rep("S", maxord), c(1:maxord), sep=""), "Stot")
volume<-as.data.frame(volume)
colnames(volume)<-c(paste(rep("V", maxord), c(1:maxord), sep=""), "Vtot")
if (maxord>1){
latroot<-as.data.frame(latroot)
for (l in 1:(maxord-1)){
colnames(latroot)[l]<-paste("N", l+1, "LR", sep="")
colnames(latroot)[l+(maxord-1)]<-paste("L", l+1, "LR", sep="")
colnames(latroot)[l+2*(maxord-1)]<-paste("ML", l+1, "LR", sep="")
colnames(latroot)[l+3*(maxord-1)]<-paste("GR", l+1, "L", sep="")}
if (fitter==TRUE) {datarsml<-data.frame(datarsml[,1:9], as.data.frame(latroot), datarsml[,10:19], as.data.frame(diameter), as.data.frame(surface), as.data.frame(volume))}
else {datarsml<-data.frame(datarsml[,1:9], as.data.frame(latroot), datarsml[,10:16], as.data.frame(diameter), as.data.frame(surface), as.data.frame(volume))}}}
if (is.null(inputrac)==TRUE & is.null(inputrsml)==FALSE) {return(datarsml)}
if (is.null(inputrsml)==TRUE & is.null(inputrac)==FALSE) {return(datadart)}
if (is.null(inputrsml)==FALSE & is.null(inputrac)==FALSE){return(smartbind(datadart, datarsml))}}}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.