# R/CM.processCenterline.r In AntoniusGolly/cmgo: Derive principle Channel metrics from bank points

#### Documented in CM.processCenterline

#' Process the centerline
#'
#' Derive width, slope and further principle channel metrics for the channel centerline previously created.
#'
#' \code{CM.processCenterline()} calculate the channel metrics (Fig. 9) based on the centerline previously calculated.
#' It does that by first deriving channel transects. The transects are lines perpendicular to a group of centerline points
#' where the size of that group is defined by the parameter \code{transects.span}. By default this span
#' equals three which means for each group of three centerline points a line is created through the outer points of
#' that group to which the perpendicular -- the transect -- is calculated (see Fig. 8b). In the final step the intersections of the
#' transects with the banks are calculated (Fig. 8c).
#'
#' \if{html}{\figure{08-processing.png}{options: width="800px" alt="Figure: processing"}}
#' \if{latex}{\figure{08-processing.pdf}{options: width=9cm}}
#' \emph{Figure 8: A visualization of the transect calculation. a) the channel centerline, b) for each
#' group of \code{n} points (in the example n=3) a line is fitted through the outer points (black line).
#' The perpendicular to this line is the transect. c) the vectors from each centerline point to the left
#' and to the right bank are calculated and stored separately.}
#'
#' When the transects cross the banks multiple times, the minimum distance is taken. In addition to the width,
#' the distances of the centerline points to the banks is stored sepearately for the left and the right bank. This
#' is particularly of importance when using multiple time series of channel profiles. See paragraph in
#' the (see \link[=cmgo]{paragraph "Time series analyses" of package documentation}.
#'
#' The function returns the global data object extended by the following variables (length equals number
#' of points of the reference centerline): \preformatted{
#' $metrics$tr       # linear equations of the transects
#' $metrics$cp.r     # coordinates of crossing points transects / right bank
#' $metrics$cp.l     # coordinates of crossing points transects / left bank
#' $metrics$d.r      # distance of reference centerline point / right bank
#' $metrics$d.l      # distance of reference centerline point / left bank
#' $metrics$w        # channel width
#' $metrics$r.r      # direction value: -1 for right, +1 for left to the centerline
#' $metrics$r.l      # direction value: -1 for right, +1 for left to the centerline
#' $metrics$diff.r   # difference between right bank point of actual time series and right bank point of reference series
#' $metrics$diff.l   # difference between left bank point of actual time series and left bank point of reference series
#'}
#' If you calculate channel metrics for one channel survey, the right bank is always to the right side of the centerline and
#' the left bank always left to the centerline. Thus, the values of \code{r.r} are all -1 and the values of \code{r.l} all +1. However,
#' when using a reference centerline to compare different channel surveys, these values can be required. If for example
#' you examine the bank metrics for a reference centerline that is very different to the actual survey, the centerline
#' is not necessarily between the two banks. This can happen when a massive shift of the channel bed occurred (see Fig. 9).
#' For further calculations consider then the sign of the \code{r.r} and \code{r.l} values!
#'
#'
#'
#' @template param_global_data_object
#' @param set an optional argument for processing a specific data set, if \code{NULL} all available data sets are used
#' @return returns the global data object extended by the centerline data \code{$metrics} for the respective data set #' @author Antonius Golly #' @examples #' # get demo data (find instructions on how to use own data in the documentation of CM.ini()) #' cmgo.obj = CM.ini("demo1") #' #' # calculate channel metrics from centerline #' cmgo.obj = CM.processCenterline(cmgo.obj) #' #' @export CM.processCenterline CM.processCenterline <- function(cmgo.obj, set=NULL){ par = cmgo.obj$par
data = cmgo.obj$data sets = if(is.null(set)) names(data) else set notice = function(x,prim=FALSE){cat(paste((if(prim) "\n--> " else " "), x, sep=""), sep="\n")} error = function(x){stop(x, call.=FALSE)} warn = function(x){warning(x, call.=FALSE)} notice("process centerline (calculate channel metrics)", TRUE) CM.calculateIntersections = function(tr, cb){ ## get special cases horizontals = which(cb[,"Ny"] - cb[,"Py"] == 0) verticals = which(cb[,"Nx"] - cb[,"Px"] == 0) # create all crossing points cp.Px = (cb[,"n"] - tr["n"]) / (tr["m"] - cb[,"m"]) cp.Px[verticals] = cb[verticals, "Px"] cp.Py = tr["m"] * cp.Px + tr["n"] # get scalar factors to filter to those which cross bank segments scalar = (cp.Py - cb[,"Py"]) / (cb[,"Ny"] - cb[,"Py"]) # scalar calculations fail for segments that horizontal, fix it: scalar[horizontals] = (cp.Px[horizontals] - cb[horizontals,"Px"]) / (cb[horizontals,"Nx"] - cb[horizontals,"Px"]) # get all crossing coordinates X = cp.Px[scalar >= 0 & scalar <= 1] Y = cp.Py[scalar >= 0 & scalar <= 1] # check length of coordinate vectors if(length(X) != length(Y)) stop("inconsistent number of crossing points (X,Y)") # find shortest section when more than one available d.ix = which.min(((X - tr["Px"])^2 + (Y - tr["Py"])^2 )^(1/2)) if(length(d.ix) == 0){ warn(paste("no intersection found for the transect at x=", tr[["Px"]], ", y=", tr[["Py"]])) return(c(NA, NA, NA)) } # get final crossing coordinate X = X[d.ix] Y = Y[d.ix] # get direction r = (X - tr["Px"]) / tr["m"] #r = X - tr["Px"] ## bug: consider not only plain X-coordinate!! r = if(r < 0) -1 else if(r > 0) 1 else 0 # return return(c(X, Y, r)) } CM.linearEquationsBanks = function(x){ if(!is.matrix(x)) return(F) # first and last elements represent averages over moving window x1 = x[1,1] x2 = x[2,1] y1 = x[1,2] y2 = x[2,2] # calculate slope (m) and center point (P) m = (( y2 - y1 ) / ( x2 - x1 )) Px = x1 Py = y1 n = Py - ( m * Px ) Nx = x2 Ny = y2 # return return(data.frame(m=m,n=n,Px=Px,Py=Py,Nx=Nx,Ny=Ny)) } CM.linearEquationsTransects = function(x){ # x has to be matrix representing x-y-pairs of centerline points in the number of transects.span if(!is.matrix(x)) return(F) # first and last elements represent averages over moving window x1 = head(as.matrix(x)[,1], n=1) x2 = tail(as.matrix(x)[,1], n=1) y1 = head(as.matrix(x)[,2], n=1) y2 = tail(as.matrix(x)[,2], n=1) # calculate orientation (m) and center point (P) of transect vectors m = -1 / ((( y2 - y1 ) / ( x2 - x1 ))) Px = (( x2 - x1 ) / 2 ) + x1 Py = (( y2 - y1 ) / 2 ) + y1 n = Py - ( m * Px ) # return return( data.frame(m=m,n=n,Px=Px,Py=Py) ) } CM.calculateWidthChange = function(cl_len, w, width_range){ width_change = list() for(w_range in width_range){ # calculates the local change of the width for a given point dw.fit = apply(as.array(cl_len), 1, function(x){ # select indices within the matching range ixs = which(cl_len >= x & cl_len < (x + w_range)) # fit a line through values (reverse indices to reflect water direction) w_fit = if(!all(is.na(w[ixs]))) lm(w[ixs] ~ rev(cl_len[ixs])) else return(NA) # return slope of fit return (w_fit$coefficients[2])

})

# calculates the local change of the width for a given point
dw.avg = apply(as.array(cl_len), 1, function(x){

# select indices within the matching range
ixs = which(cl_len >= x & cl_len < (x + w_range))

# return slope of average
return ((w[head(ixs, n=1)] - w[tail(ixs, n=1)]) / (cl_len[tail(ixs, n=1)] - cl_len[head(ixs, n=1)]))

})

# calculates the local difference of the width for a given point
dw.d = apply(as.array(cl_len), 1, function(x){

# select indices within the matching range
ixs = which(cl_len >= x & cl_len < (x + w_range))

# return difference
return (w[head(ixs, n=1)] - w[tail(ixs, n=1)])

})

width_change[[paste0("dw_fit_", w_range)]] = dw.fit
width_change[[paste0("dw_avg_", w_range)]] = dw.avg
width_change[[paste0("dw_d_",   w_range)]] = dw.d

}

return(width_change)

}

set = "set1" # ignored, only for debugging

for(set in sets){

if(par$calculate.metrics){ if(is.null(data[[set]]$metrics) || par$force.calc.metrics || !identical(data[[set]]$metrics$tr.span, par$transects.span)){

cl.type = "smoothed" # it is highly recommended to do the processing on the "smoothed" version of the centerline
set.ref  = if(par$centerline.use.reference) par$centerline.reference else set

notice(paste("process data set '", set, "' with centerline reference '", set.ref,"' now...", sep=""), TRUE)
reasons = c()
if(is.null(data[[set]]$metrics)) reasons = append(reasons, "data does not exis") if(par$force.calc.metrics)           reasons = append(reasons, "calculation is forced")
if(!identical(data[[set]]$metrics$tr.span, par$transects.span)) reasons = append(reasons, "tr spans differ") notice(paste("reason for calculation:", paste(reasons, collapse = " + "))) ### calculate linear equations of transsects ############## tr = rollapply(as.data.frame(data[[set.ref]]$cl[[cl.type]][c("x","y")]), par$transects.span, CM.linearEquationsTransects, by.column=FALSE, partial=TRUE) if(par$transects.span==2) tr[nrow(tr),] = tr[nrow(tr)-1,] # when span==2 and partial=TRUE the last element will be 0,0,0,0 which will be substituted
notice(paste("linear equations for transects calculated with span = ", par$transects.span,"!", sep="")) ### calculate linear equations of bank segments ########### cb.r = rollapply(data.frame(x = data[[set]]$channel$x[data[[set]]$ix.r], y = data[[set]]$channel$y[data[[set]]$ix.r]), 2, CM.linearEquationsBanks, by.column=FALSE) cb.l = rollapply(data.frame(x = data[[set]]$channel$x[data[[set]]$ix.l], y = data[[set]]$channel$y[data[[set]]$ix.l]), 2, CM.linearEquationsBanks, by.column=FALSE) notice("linear equations for bank segments calculated!") ### calculate intersections of transsects and bank segments cp.r = t(apply(tr, 1, function(x){ CM.calculateIntersections(x, cb.r) })) cp.l = t(apply(tr, 1, function(x){ CM.calculateIntersections(x, cb.l) })) notice("intersections calculated!") ### calculate distances ################################### d.r = apply(cbind(cp.r[,c(1,2)], tr[,c(3,4)]), 1, function(x){ ( (x[1]-x[3])^2 + (x[2]-x[4])^2 )^(1/2) }) d.l = apply(cbind(cp.l[,c(1,2)], tr[,c(3,4)]), 1, function(x){ ( (x[1]-x[3])^2 + (x[2]-x[4])^2 )^(1/2) }) ### calculate width ####################################### w = apply(cbind(cp.r[,c(1,2)], cp.l[,c(1,2)]), 1, function(x){ ( (x[1]-x[3])^2 + (x[2]-x[4])^2 )^(1/2) }) ### calculate width change ################################ dw = CM.calculateWidthChange(data[[set.ref]]$cl[[cl.type]]$cum_dist_2d, w, par$metrics.local.change.range)

notice("width calculated!")

# check
#if(length(data[[set]]$cl[[cl.type]]$x) != length(cw)) error("length of CL and CW differ!")
if(length(data[[set.ref]]$cl[[cl.type]]$x) != length(d.r)) error("length of CL and CW differ!")
if(length(data[[set.ref]]$cl[[cl.type]]$x) != length(d.l)) error("length of CL and CW differ!")

# store
data[[set]]$metrics$cl.ref      = set.ref
data[[set]]$metrics$cl.type     = cl.type
data[[set]]$metrics$tr.span     = par$transects.span data[[set]]$metrics$tr = tr data[[set]]$metrics$cp.r = cp.r[,c(1,2)] data[[set]]$metrics$cp.l = cp.l[,c(1,2)] data[[set]]$metrics$d.r = d.r data[[set]]$metrics$d.l = d.l data[[set]]$metrics$w = w data[[set]]$metrics$r.r = cp.r[,3] data[[set]]$metrics$r.l = cp.l[,3] data[[set]]$metrics$diff.r = if(set == set.ref) rep(0, length(w)) else -(data[[set.ref]]$metrics$d.r * data[[set.ref]]$metrics$r.r - d.r * cp.r[,3]) data[[set]]$metrics$diff.l = if(set == set.ref) rep(0, length(w)) else data[[set.ref]]$metrics$d.l * data[[set.ref]]$metrics$r.l - d.l * cp.l[,3] data[[set]]$metrics             = modifyList(data[[set]]$metrics, dw) ### project features onto centerline if(typeof(data[[set]]$features) == "list"){

notice("feature list found!")

cl = data[[set]]$cl nrp = length(cl$smoothed$x) ## create ppp feature centerline_points = ppp( cl$smoothed$x, cl$smoothed$y, window=owin(range(cl$smoothed$x), range(cl$smoothed$y)) ) for(feature in names(data[[set]]$features)){

notice(paste("processing feature '", feature, "'...", sep=""))

if(!is.null(data[[set]]$features[[feature]]$x) && !is.null(data[[set]]$features[[feature]]$y)){

## create ppp feature
feature_points = ppp(
data[[set]]$features[[feature]]$x,
data[[set]]$features[[feature]]$y,
window=owin(range(data[[set]]$features[[feature]]$x), range(data[[set]]$features[[feature]]$y))
)

# project centerline to feature and store in $feature prj_to_feature = nncross(feature_points, centerline_points) data[[set]]$features[[feature]]$cl_index = prj_to_feature$which

# project feature to centerline and store in $centerline prj_to_cl = nncross(centerline_points, feature_points) data[[set]]$metrics[[paste(feature, "_which", sep="")]] = prj_to_cl$which data[[set]]$metrics[[paste(feature, "_dist",  sep="")]] = prj_to_cl$dist for(variable in names(data[[set]]$features[[feature]][- which(names(data[[set]]$features[[feature]]) %in% c("x", "y", "cl_index"))])){ data[[set]]$metrics[[paste(feature, "_", variable,  sep="")]] = data[[set]]$features[[feature]][[variable]][prj_to_cl$which]
}

# calculate 2d dist
diffs = diff(as.matrix(data[[set]]$features[[feature]][which(names(data[[set]]$features[[feature]]) %in% c("x", "y"))]))
data[[set]]$features[[feature]]$seg_dist_2d = c(0, sqrt(diffs[,"x"]^2 + diffs[,"y"]^2))
data[[set]]$features[[feature]]$cum_dist_2d = cumsum(data[[set]]$features[[feature]]$seg_dist_2d)

# calculate 3d dist
if(!is.null(data[[set]]$features[[feature]]$z)){
diffs = diff(as.matrix(data[[set]]$features[[feature]][which(names(data[[set]]$features[[feature]]) %in% c("x", "y", "z"))]))
data[[set]]$features[[feature]]$seg_dist_3d = c(0, sqrt(diffs[,"x"]^2 + diffs[,"y"]^2 + diffs[,"z"]^2))
data[[set]]$features[[feature]]$cum_dist_3d = cumsum(data[[set]]$features[[feature]]$seg_dist_3d)
}

notice(paste("feature '", feature, "' projected to centerline", sep=""))

} else { warning("feature does not contain proper x,y coordinates lists!") }

}

}

notice("features processed successfully!")

# # #
if(par$steps.identify){ notice("identify steps in long profile...", TRUE) if(!is.null(data[[set]]$features$lp)){ # set parameters sl_min = par$steps.minimum.step.length    * par$steps.bank.full.width / 100 sl_max = par$steps.maximum.step.length    * par$steps.bank.full.width / 100 pl_min = par$steps.minimum.pool.length    * par$steps.bank.full.width / 100 rd_min = par$steps.minimum.residual.depth * par$steps.bank.full.width / 100 dh_min = par$steps.minimum.drop.height    * par$steps.bank.full.width / 100 ssl_min = if(par$steps.average.slope.fix) par$steps.average.slope + 10 else { # calculate minimum slope based on elevation and length of long profile ( atan(diff(range(data[[set]]$features$lp$z)) / tail(data[[set]]$features$lp$cum_dist_2d, n=1)) * 360 / (2*pi) ) + 10 } ssg_min = tan(ssl_min/360 * (2*pi)) * 100 notice("****** parameters **************************************************************") notice(paste0("min step length (sl_min = ", par$steps.minimum.step.length,    "% of bank full width): ",   sl_min, "[m]"));
notice(paste0("max step length (sl_max = ",    par$steps.maximumg.step.length, "% of bank full width): ", sl_max, "[m]")); notice(paste0("min pool length (pl_min = ", par$steps.minimum.pool.length,    "% of bank full width): ",   pl_min, "[m]"));
notice(paste0("min residual depth (rd_min = ", par$steps.minimum.residual.depth, "% of bank full width): ", rd_min, "[m]")); notice(paste0("min drop height (dh_min = ", par$steps.minimum.drop.height,    "% of bank full width): ",   dh_min, "[m]"));
notice(paste0("min step slope (ssl_min = ",                                      "mean slope + 10 degree): ", ssl_min, "[deg]"));
notice("********************************************************************************")

# prepare variables
lp           = data[[set]]$features$lp
steps        = data.frame()
verbose      = par$steps.verbose step.i = 0 nr.of.checks = 4 crash.i = 0 ix = 1 notice("identified steps:") checks = "start iteration:" while(ix < length(lp$id)){

crash.i = crash.i + 1; if(crash.i > 50000){ stop("program crashed!"); break;}

#if(Parameters$step.max) if(nrow(steps) >= Parameters$step.max.count) break;

if(verbose) notice(paste("index: ", ix, " / dist_cum: ", round(lp$cum_dist_2d[ix], digits=2), sep="")) # reset checks checks = rep(FALSE, nr.of.checks) # set indices ixl = ix + 1 # start looking from the next point on ixr = ix + 40 # end after 40 points if(ixr > length(lp$id)) ixr = length(lp$id) # actual point lpx = lp$x[ix]
lpy  = lp$y[ix] lpz = lp$z[ix]
lpd  = lp$cum_dist_2d[ix] # i-points to check lpxi = lp$x[c(ixl:ixr)]                               # all points upstream of ix (end of pool) excluding ix
lpyi = lp$y[c(ixl:ixr)] # all points upstream of ix (end of pool) excluding ix lpzi = lp$z[c(ixl:ixr)]                               # all points upstream of ix (end of pool) excluding ix
lpdi = lp$cum_dist_2d[c(ixl:ixr)] # cumulative distances lpfi = ( ( lpxi - lpx )^2 + ( lpyi - lpy )^2 )^(1/2) # real distances from end of pool upstream # get distances of points dist3d = ( ( lpxi - lpx )^2 + ( lpyi - lpy )^2 + ( lpzi - lpz )^2 )^(1/2) dist2d = ( ( lpxi - lpx )^2 + ( lpyi - lpy )^2 )^(1/2) ### check 1: higher than the next point? if(lpz > lpzi[1]) checks[1] = TRUE else{ ix = ix+1; next } ### check 2: highest point in elevation among all points upstream that fall within a minimum pool length if(all(lpzi < lpz)){ notice("end of lp reached"); break;} first_higher_point = min(which((lpz > lpzi) == FALSE)) last_lower_point = first_higher_point - 1 # pl_min.pts = which(dist2d > pl_min) if(dist2d[last_lower_point]>pl_min) checks[2] = TRUE else { ix = ix+1; next } # generate interpolated point lp_llp = c(lpxi[last_lower_point], lpyi[last_lower_point], lpzi[last_lower_point]) # LongProfile_LastLowerPoint lp_fhp = c(lpxi[first_higher_point], lpyi[first_higher_point], lpzi[first_higher_point]) # LongProfile_FirstHigherPoint lp_h = lp_fhp - lp_llp # LongProfile_Helpingpoint lambda = (lpz -lp_llp[3]) / lp_h[3] lp_sp = c(lp_llp[1] + (lp_h[1] * lambda), lp_llp[2] + (lp_h[2] * lambda), lpz) # LongProfile_StartPool (interpolated point) # calculate cum. dist of interpolated start of pool point (for plotting, cum. distances are required) lpd_sp = lpdi[last_lower_point] - lpd # LongProfileDistance_StartPool if(par$steps.thalweg.dist == "2d") lpd_sp = lpd_sp + ( ( lp_sp[1] - lpxi[last_lower_point] )^2 + ( lp_sp[2] - lpyi[last_lower_point] )^2 ) ^(1/2)
if(par$steps.thalweg.dist == "3d") lpd_sp = lpd_sp + ( ( lp_sp[1] - lpxi[last_lower_point] )^2 + ( lp_sp[2] - lpyi[last_lower_point] )^2 + ( lp_sp[3] - lpzi[last_lower_point] )^2 ) ^(1/2) # calculate real pool length pool.length = ( ( lp_sp[1] - lpx )^2 + ( lp_sp[2] - lpy )^2 ) ^(1/2) # calculate residual depth pool.depth.ix = which.min(lp$z[ix + seq(1, last_lower_point)])
pool.depth.z  = min( lp$z[ix + seq(1, last_lower_point)] ) pool.depth.r = lpz - pool.depth.z ### check 3: min residual depth if(pool.depth.r >= rd_min) checks[3] = TRUE else { ix = ix + first_higher_point; next } # re-calculate real distances step.lengths = ( ( lpxi[first_higher_point:length(lpxi)] - lp_sp[1] )^2 + ( lpyi[first_higher_point:length(lpyi)] - lp_sp[2] )^2 )^(1/2) # real distances from end of pool upstream # find start of step (first point above lp_sp for which min. step drop height and min step length applies) if(!any( (lpzi[first_higher_point:length(lpzi)] - lpz) > dh_min)){ notice("end of lp reached"); break;} first_top_point = first_higher_point - 1 + max( min(which( (lpzi[first_higher_point:length(lpzi)] - lpz) > dh_min)), # minimum drop height min(which( step.lengths > sl_min )) # minimum step length (not pool length), step length == length from start of pool (end of step) to first_top_point ) # point cloud calculations ############################################################ gr.ix = 0 r.squared = 0; ssl = 0 group = list() pointcloud = TRUE while(pointcloud){ if(verbose) notice(paste("_________________ ix", ix, "_________________")) crash.i = crash.i + 1; if(crash.i > 5000){ error("program crashed!"); break;} if(first_top_point + gr.ix > length(lpxi)) {notice("end of point cloud!"); break;} # point cloud of points to fit slx = c(lp_sp[1], lpxi[first_higher_point : ( first_top_point + gr.ix )]) sly = c(lp_sp[2], lpyi[first_higher_point : ( first_top_point + gr.ix )]) slz = c(lp_sp[3], lpzi[first_higher_point : ( first_top_point + gr.ix )]) slz0 = rep(min(slz), length(slz)) # re-calculate step length step.length = ( ( slx[1] - slx[length(slx)] )^2 + ( sly[1] - sly[length(sly)] )^2 )^(1/2) # real distances from end of pool upstream # check min and max step length if(step.length < sl_min){ if(verbose) notice("minimum step length is not reached!") error("program aborted") pointcloud = FALSE } if(step.length > sl_max){ if(verbose) notice("maximum step length reached!") if(verbose) notice("No step detected. Discard and go to next...") ix.add = first_higher_point - 1; pointcloud = FALSE break; } # gradient index gr.ix = gr.ix + 1 # 3 dimensional fit fit.3d = TRUE if(fit.3d){ m = 0 r.squared = 0 if(length(slx) < 2){ notice("point cloud only contains 1 point!") error("program aborted!") } else if(length(slx) == 2){ if(verbose) notice("point cloud contains two points."); if(verbose) notice("set r2 to 0.5") r.squared = 0.5; m = (slz[2]-slz[1]) / (( (slx[2]-slx[1])^2 + (sly[2]-sly[1]) ^2) ^(1/2)) } else { if(verbose) notice(paste("point cloud size:", length(slx))) if(verbose) notice(paste("step length [m]:", round(step.length,1))) ### PCA for bottom line ######################## xyz <- data.frame(x = slx, y = sly, z = slz) xyz0 <- data.frame(x = slx, y = sly, z = slz0) N <- nrow(xyz) mean_xyz <- apply(xyz, 2, mean) mean_xyz0 <- apply(xyz0, 2, mean) xyz_pca <- princomp(xyz) xyz_pca0 <- princomp(xyz0) dirVector <- xyz_pca$loadings[, 1]    # PC1
dirVector0 <- xyz_pca0$loadings[, 1] # PC1 xyz_fit <- matrix(rep(mean_xyz, each = N), ncol=3) + xyz_pca$score[, 1] %*% t(dirVector)
xyz_fit0   <- matrix(rep(mean_xyz0, each = N), ncol=3) + xyz_pca0$score[, 1] %*% t(dirVector0) t_ends <- c(min(xyz_pca$score[,1]) - 0.2, max(xyz_pca$score[,1]) + 0.2) # for both ends of line t_ends0 <- c(min(xyz_pca0$score[,1]) - 0.2, max(xyz_pca0$score[,1]) + 0.2) # for both ends of line endpts <- rbind(mean_xyz0 + t_ends0[1]*dirVector0, mean_xyz0 + t_ends0[2]*dirVector0) ### interactive plot ############################ #requires library(rgl) if(verbose){ clear3d() #Sys.sleep(2) plot3d(xyz, type="p", size=10) #rgl.bbox(xlen = 0, ylen = 0, zlen = 0, color = c('grey100')) #Wait() # plot fits abclines3d(mean_xyz, a = dirVector, col="orange", lwd=2) abclines3d(mean_xyz0, a = dirVector0, col="blue", lwd=2) # mean + t * direction_vector # plot connections for(i in 1:N) segments3d(rbind(xyz[i,], xyz_fit[i,]), lty=2,col="green2") for(i in 1:N) segments3d(rbind(xyz[i,], xyz_fit0[i,]), col="gray") # plot secondary coordinate system for(i in 1:N) points3d(rbind(c(xyz_fit0[i,c(1,2)], slz[i])), size=8, col="red") for(i in 1:N) points3d(rbind(xyz_fit0[i,]), size=4, col="red") } origin = xyz_fit0[1,] # fit lm() in 2d sub-space x2d = ( (xyz_fit0[,1] - origin[1])^2 + (xyz_fit0[,2] - origin[2])^2 ) ^ (1/2) y2d = slz final_fit = lm(y2d ~ x2d) # plot this final fit in 2d space if(verbose) plot(y2d ~ x2d) if(verbose) abline(final_fit, col="red") # extract fit parameters n = final_fit$coefficients[1]
m = final_fit$coefficients[2] r.squared = summary(final_fit)$r.squared

### plot this final fit in 3d space
dirVectorFit = c(dirVector0[c(1,2)], m)
intercept = c(xyz_fit[1,c(1,2)], n)
if(verbose) abclines3d(intercept, a = dirVectorFit, col="red")

#stop()

} # if(length(slx) < 2) {} else if(length(slx) == 2) {} else

#if(gr.ix==2) stop("end")

# eval slope
ssl = atan( m ) / (2*pi) * 360
if(verbose) notice(paste("step slope [deg]", round(ssl,2)))

### check 4: min step slope
if(ssl >= ssl_min){
group = rbind(group, data.frame(
"ix"  = gr.ix,
"n"   = length(slx),
"ssl" = ssl,
"r2"  = r.squared,
"sl"  = step.length,
"px"  = slx[length(slx)],
"py"  = sly[length(sly)],
"pz"  = slz[length(slz)],
"dummy" = TRUE
))
checks[4] = TRUE
} else {
if(verbose) notice("ssl_min not reached!");
}

} else {} # 2 dimensional fit to come

#Wait();

} # while(pointcloud)

### step identified ###################################################################

if(all(checks)){

### find start of step

stepgroup = group
group = stepgroup
bestfit = which.max(stepgroup$r2) start.of.step = first_top_point - 1 + group$ix[bestfit]
step.slope    = group$ssl[bestfit] step.i = step.i + 1 cat(paste("#",step.i, " ", sep="")) if(step.i %% 20 == 0) cat("\n") ### store step data steps = rbind(steps, data.frame( "id" = step.i, "ix" = ix, "pool.end.x" = lpx, "pool.end.y" = lpy, "pool.end.z" = lpz, "pool.end.d" = lp$cum_dist_2d[ix],

"pool.start.x"    = lp_sp[1],
"pool.start.y"    = lp_sp[2],
"pool.start.z"    = lp_sp[3],
"pool.start.d"    = lp$cum_dist_2d[ix] + lpd_sp, "pool.length" = pool.length, "pool.dist" = lpd_sp, "pool.depth.x" = lp$x[ix + pool.depth.ix],
"pool.depth.y"    = lp$y[ix + pool.depth.ix], "pool.depth.z" = pool.depth.z, "pool.depth.r" = pool.depth.r, "pool.depth.ix" = ix + pool.depth.ix, "pool.depth.d" = lp$cum_dist_2d[ix + pool.depth.ix],

"step.height"     = lpzi[start.of.step] - pool.depth.z,

"start.of.step.x" = lpxi[start.of.step],
"start.of.step.y" = lpyi[start.of.step],
"start.of.step.z" = lpzi[start.of.step],
"start.of.step.d" = lpdi[start.of.step],

"top.first.x"     = lpxi[first_top_point],
"top.first.y"     = lpyi[first_top_point],
"top.first.z"     = lpzi[first_top_point],
"top.first.d"     = lpdi[first_top_point],

"tx.id"           = as.character(step.i),
"tx.rd"           = as.character(round(pool.depth.r, digits=2)),
"tx.pl"           = as.character(round(pool.length, digits=2)),

"dh_min"          = lpz + dh_min,

stringsAsFactors  = FALSE

))

ix = ix + first_higher_point

} else {

# proceed with next ix
ix = ix + 1 + ix.add

}

} # while(ix < ixs[2])

# store
data[[set]]$steps = steps } else { # if(!is.null(data[[set]]$features$lp) warning("feature list does not contain long profile information to identify steps!") } } # if(par$steps.identify) # # #

} else { # if(is.null(data[[set]]$metrics) || par$force.calc.metrics)

notice("(metrics taken from cache)")

}

} # if(par\$calculate.width)

} # for(set in names(data))

notice("CM.processCenterline() has ended successfully!", TRUE)

# return
return(list(
data = data,
par  = par
))

}

AntoniusGolly/cmgo documentation built on April 16, 2021, 2:52 p.m.