R/cpplot3d.bottom.TSD.R

In arnejohannesholmin/cpplot3d: 3D-plotting and segmentation of sonar data

#*********************************************
#*********************************************
#' Adding bottom surface from segmented bottom voxels.
#'
#' @param data  is the list of required input data, specifically, the bottom segmentation masks (sgbt); the lengths of the beams (lenb); the number of beams (optional, numb); the voxel positions in the x and y direction (psxx and psyx, respectively); and the bottom points in the x, y, and z direction, generated by pplot3d.TSD() (psxr, psyr, and pszr, respectively).
#' @param bottom.res  is a two element vector of bottom.res distances in the x and y direction respectively. Grid points are positioned so that the end points of the range of the data in each direction is kept, and the distance between grid points is approximately equal to the values given in 'bottom.res'. If 'bottom.res' is given as a single value, it is recycled to both directions.
#' @param bottom.smooth  is the width of the Gaussian kernel by which to smooth the bottom surface.
#' @param col  is either a list of parameters used in colscale(), such as x, fun=c("rainbow","grey","heat.colors","terrain.colors","topo.colors","cm.colors","combined"), start=0, end=0.8, h=1/6, s=1, v=1, sc=c(1,0.3), vc=c(0.7,1), flip=FALSE, bias=1, space=c("rgb","Lab"), interpolate=c("linear","spline","fadewhite","fadeblack"), breakpoint=NULL, alpha=1.
#' @param step.size  is the number of points used to generate the bottom at each step, set low to avoid memory limitations slowing down the generation of bottom points.
#'
#' @return
#'
#' @examples
#' \dontrun{}
#'
#' @importFrom akima interp.old
#' @importFrom rgl rgl.surface
#' @importFrom sonR ksmooth.SG
#' @importFrom TSD zeros
#' @importFrom stats quantile
#'
#' @export
#' @rdname cpplot3d.bottom.TSD
#'
cpplot3d.bottom.TSD<-function(data, bottom.res=2, bottom.smooth=NULL, col="topo.col", step.size=1e4, insteps=TRUE, ...){
	
	############ AUTHOR(S): ############
	# Arne Johannes Holmin
	############ LANGUAGE: #############
	# English
	############### LOG: ###############
	# Start: 2013-09-18 - Clean version.
	# Last: 2013-10-07 - Changed name.
	########### DESCRIPTION: ###########
	# Adding compensation for absorption and attenuation for a set of sound beams (Time Varied Gain).
	########## DEPENDENCIES: ###########
	#
	############ VARIABLES: ############
	# ---data--- is the list of required input data, specifically, the bottom segmentation masks (sgbt); the lengths of the beams (lenb); the number of beams (optional, numb); the voxel positions in the x and y direction (psxx and psyx, respectively); and the bottom points in the x, y, and z direction, generated by pplot3d.TSD() (psxr, psyr, and pszr, respectively).
	# ---bottom.res--- is a two element vector of bottom.res distances in the x and y direction respectively. Grid points are positioned so that the end points of the range of the data in each direction is kept, and the distance between grid points is approximately equal to the values given in 'bottom.res'. If 'bottom.res' is given as a single value, it is recycled to both directions.
	# ---bottom.smooth--- is the width of the Gaussian kernel by which to smooth the bottom surface.
	# ---col--- is either a list of parameters used in colscale(), such as x, fun=c("rainbow","grey","heat.colors","terrain.colors","topo.colors","cm.colors","combined"), start=0, end=0.8, h=1/6, s=1, v=1, sc=c(1,0.3), vc=c(0.7,1), flip=FALSE, bias=1, space=c("rgb","Lab"), interpolate=c("linear","spline","fadewhite","fadeblack"), breakpoint=NULL, alpha=1.
	# ---step.size--- is the number of points used to generate the bottom at each step, set low to avoid memory limitations slowing down the generation of bottom points.
	
	
	##################################################
	##################################################
	##### Preparation #####
	# Repeat 'bottom.res':
	if(length(bottom.res)==1){
		bottom.res = c(bottom.res,bottom.res)
		}
	
	# Get the bottom grid:
	if(length(data$numb)==0){
		data$numb = length(data$lenb)
		}
	
	# Get the ranges of the grids:
	if(length(data$psxr)>0){
		minxgrid = min(unlist(lapply(data$psxr,min)))
		maxxgrid = max(unlist(lapply(data$psxr,max)))
		minygrid = min(unlist(lapply(data$psyr,min)))
		maxygrid = max(unlist(lapply(data$psyr,max)))
		# Calculate the rounded number of grid cells for x and y, and create sequences with those number of cells:
		nx = round((maxxgrid-minxgrid)/bottom.res[1])+1
		xgrid = seq(minxgrid, maxxgrid, length=nx)	
		ny = round((maxygrid-minygrid)/bottom.res[2])+1
		ygrid = seq(minygrid, maxygrid, length=ny)	
		
		##### Execution and output #####
		if(insteps){
			# Interpolate the bottom in steps, in order to reduce memory in each step:
			l = sum(sapply(data$psxr,length))
			steps = max(1,round(l/step.size))
			d = sqrt((maxygrid-minygrid) * (maxxgrid-minxgrid) / steps)
			stepsx = ceiling((maxygrid-minygrid)/d)
			stepsy = ceiling((maxxgrid-minxgrid)/d)
			x = unlist(data$psxr)
			y = unlist(data$psyr)
			z = unlist(data$pszr)
			# Put the step separators on the edges between grid cells in the x-direction:
			sepx = unique(round(quantile(x,seq(0,1,length.out=stepsx+1)),bottom.res))
			sepy = unique(round(quantile(y,seq(0,1,length.out=stepsy+1)),bottom.res))
			# Get the indices for each step:
			indicesx = findInterval(x,sepx,all.inside=TRUE)
			indicesx = split(seq_len(l),indicesx)
			indicesy = findInterval(y,sepy,all.inside=TRUE)
			indicesy = split(seq_len(l),indicesy)
			
			# Reserve space for the bottom:
			bottom = list(x=xgrid,y=ygrid,z=vector("list",steps))
			
			# Define parameters used when plotting the time bar for the generation of bottom points:
			infostring = "Generating bottom points:"
			cat(infostring,"\n",sep="")
			totalsteps = stepsx*stepsy
			stepfact = nchar(infostring)/totalsteps
			oldvalue = 0
			
			# Store the valid steps, which are those with more than 'minpoints' points
			minpoints = 5
			validsteps = 0
			npointsinsteps = zeros(totalsteps)
			
			# Run through the steps in both directions and generate bottom points:
			for(i in seq_len(stepsy)){
				for(j in seq_len(stepsx)){
					thisij = (i-1)*stepsx + j
					# Print a dot if the floor of the new value exceeds the old value in:
					thisvalue = floor(thisij*stepfact)
					if(thisvalue > oldvalue){
						cat(rep(".",thisvalue-oldvalue),if(thisij==totalsteps) "\n", sep="")
						oldvalue = thisvalue
						}
					# The the points present in the current field, and store the length:
					thisindices = intersect(indicesx[[j]], indicesy[[i]])
					npointsinsteps[thisij] = length(thisindices)
					# Interpolate if the field contains enough points:
					if(length(thisindices)>minpoints){
						validsteps = validsteps+1
						thisbottom = interp.old(x=x[thisindices], y=y[thisindices], z=z[thisindices], xo=xgrid, yo=ygrid, duplicate="strip")
						bottom$z[[thisij]] = thisbottom$z
						}
					}
				}
			# Combine the bottom points:
			bottom$z = unlist(bottom$z)
			dim(bottom$z) = c(nx,ny,validsteps)
			bottom$z = apply(bottom$z,1:2,mean,na.rm=TRUE)
			}
		else{
			bottom = interp.old(x=unlist(data$psxr), y=unlist(data$psyr), z=unlist(data$pszr), xo=xgrid, yo=ygrid, duplicate="strip")
			}
			
		# Get the number of NAs in the 3 x 3 neighborhood around every bottom surface cell (this is done in an unelegant way...):
		ManyNAs = zeros(dim(bottom$z),9)
		ManyNAs[-1,-1,1] = bottom$z[-nx,-ny]
		ManyNAs[,-1,2] = bottom$z[,-ny]
		ManyNAs[-nx,-1,3] = bottom$z[-1,-ny]
		ManyNAs[-1,,4] = bottom$z[-nx,]
		ManyNAs[,,5] = bottom$z[,]
		ManyNAs[-nx,,6] = bottom$z[-1,]
		ManyNAs[-1,1,7] = bottom$z[-nx,ny]
		ManyNAs[,1,8] = bottom$z[,ny]
		ManyNAs[-nx,1,9] = bottom$z[-1,ny]
		# Count NAs, and convert to TRUE if the number of NAs exceeds 2 in the neighborhood around and including a bottom surface cell:
		dim(ManyNAs) = c(prod(dim(ManyNAs)[1:2]),9)
		ManyNAs = rowSums(is.na(ManyNAs))>2
		
		# Add Gaussian kernel smoother of the bottom:
		if(length(bottom.smooth)>0 && bottom.smooth>0){
			bottom$z = ksmooth.SG(as.matrix(expand.grid(bottom$x,bottom$y)),bottom$z,h=bottom.smooth)$x
			}
			
		# Create the colors:
		if(is.character(col)){
			col = list(fun=col)
			}
		else if(is.function(col)){
			col = list(fun=deparse(substitute(col)))
			}
		# If 'start' is given as a negative value setting the heigth of the span of the colors in the bottom surface, or as a vector of two negative elements, where the first is this heigth and the second is the depth at which the lowest color occurs, the bottom colors are plotted accordingly:
		ll = list(...)
		if(length(ll$start)>0 && ll$start[1]<0){
			# Set the midpoint to 0.5 if not given:
			if(length(ll$start)==1){
				ll$start[2] = min(bottom$z,na.rm=TRUE)
				}
			
			bottom$z_plot = bottom$z
			# Subtract the midpoint:
			bottom$z_plot = bottom$z_plot - ll$start[2]
			# Scale by the depth value:
			bottom$z_plot = bottom$z_plot/(-ll$start[1])
			# If any values are below 0 or above 1, clamp these to [0,1]:
			if(min(bottom$z_plot,na.rm=TRUE)<0){
				bottom$z_plot[bottom$z_plot<0] = 0
				}
			if(max(bottom$z_plot,na.rm=TRUE)>1){
				bottom$z_plot[bottom$z_plot>1] = 1
				}
			ll$start = max(bottom$z_plot,na.rm=TRUE)
			ll$end = min(bottom$z_plot,na.rm=TRUE)
			}
		else{
			bottom$z_plot = bottom$z
			}
		col = do.call("colscale", c(list(x=bottom$z_plot), col, ll))
		# Insert NAs for the edge cells, which were identified by many NAs above:
		col[ManyNAs] = NA
	
		# Plot the bottom:
		rgl.surface(bottom$x, bottom$y, bottom$z, coords=c(1,3,2), back="lines", color=col)
		invisible(bottom)
		}
	##################################################
	##################################################
	}