# https://dahtah.github.io/imager/imager.html
library(imager)
library(config)
library(PlaneGeometry)

# check ffmpeg present
has.ffmpeg <- function()
{
	Sys.which("ffmpeg") %>% { nchar(.) > 0 }
}

# dump video frames in tempdir
load.video <- function(fname,maxSize=1,skip.to=0,frames=NULL,fps=NULL,extra.args="",verbose=FALSE)
{
	if (!has.ffmpeg()) stop("Can't find ffmpeg. Please install.")
	dd <- paste0(tempdir(),"/vid")
	unlink(dd,recursive=TRUE)
	if (!is.null(frames)) extra.args <- sprintf("%s -vframes %i ",extra.args,frames)
	if (!is.null(fps)) extra.args <- sprintf("%s -vf fps=%.4d ",extra.args,fps)
	
	arg <- sprintf("-i %s %s -ss %s ",fname,extra.args,as.character(skip.to)) %>% paste0(dd,"/image-%d.jpg")
	tryCatch({
	dir.create(dd)
	system2("ffmpeg",arg,stdout=verbose,stderr=verbose)
	fls <- dir(dd,full.names=TRUE)
	if (length(fls)==0) stop("No output was generated")
	ordr <- stringr::str_extract(fls,"(\\d)+\\.jpg") %>% stringr::str_split(stringr::fixed(".")) %>% purrr::map_int(~ as.integer(.[[1]])) %>% order
	fls <- fls[ordr]
	fls
	},
	finally={})
}

#plotImages <- FALSE
plotImages <- TRUE

# camera FOV
fov <- 79

# tape length (cm)
length <- 6.0

# infiltrometer diameter (cm)
diameter <- 10.0

# wall thickness (cm)
wall <- 0.5 

# refractive index
rindex <- 1.49 # methacrylate

# load video frames
#vidf <- 'trama_agr_3.mp4' # centered
vidf <- '15_6_20_PF_NT.mp4' # not-centered
images <- load.video(vidf, extra.args="-q:v 1")

#############################

# read config file or calibrate
ymlf <- gsub('mp4$', 'yml', vidf)
if( file.exists(ymlf) == TRUE )
{
	config <- config::get(file = ymlf)
	start <- config$start
	x1  <- config$x1
	y1  <- config$y1
	x2  <- config$x2
	y2  <- config$y2
	x3  <- config$x3
	y3  <- config$y3
	fov <- config$fov
	wall <- config$wall
	length <- config$length
} else {
	# compute diff of first frames
	motion <- c()
	curr <- load.image(images[1])
	for (image in images[2:50])
	{
		nxt <- load.image(image)
		diff <- sum(abs(nxt - curr))
		curr <- nxt
		motion <- c(motion, diff)
	}
	motion <- c(motion, diff)

	# plot to locate the start frame
	plot(seq(1,50), motion)
	loc <- graphics::locator(1, "n")
	start <- round(loc$x)

	# remove margins in plot
	op <- par(mar = rep(0, 4))

	# click two points over the horizontal line separated 6cm, i.e. from 46 to 40
	image <- load.image(images[start])
	plot(image)
	loc <- graphics::locator(2, "n")

	# zoom and click again two points
	size <- 150
	plot(imsub(image,x>loc$x[1]-size & x<loc$x[1]+size, y>loc$y[1]-size & y<loc$y[1]+size))
	points(size, size, col = rgb(red = 0, green = 1, blue = 0, alpha = 0.2), pch = 16, cex = 4)
	loc2 <- graphics::locator(1, "n")
	plot(imsub(image,x>loc$x[2]-size & x<loc$x[2]+size, y>loc$y[2]-size & y<loc$y[2]+size))
	points(size, size, col = rgb(red = 0, green = 1, blue = 0, alpha = 0.2), pch = 16, cex = 4)
	loc3 <- graphics::locator(1, "n")

	# click one point along the interface
	image <- load.image(images[start])
	plot(image)
	loc4 <- graphics::locator(1, "n")

	# zoom and click again
	size <- 150
	plot(imsub(image,x>loc4$x[1]-size & x<loc4$x[1]+size, y>loc4$y[1]-size & y<loc4$y[1]+size))
	points(size, size, col = rgb(red = 0, green = 1, blue = 0, alpha = 0.2), pch = 16, cex = 4)
	loc5 <- graphics::locator(1, "n")

	x1 <- loc$x[1] + loc2$x[1] - size
	y1 <- loc$y[1] + loc2$y[1] - size
	x2 <- loc$x[2] + loc3$x[1] - size
	y2 <- loc$y[2] + loc3$y[1] - size
	x3 <- loc4$x[1] + loc5$x[1] - size
	y3 <- loc4$y[1] + loc5$y[1] - size

	# write config file
	cat("default:", file=ymlf, sep="\n")
	cat(paste0("  start: ", start), file=ymlf, sep="\n", append=TRUE)
	cat(paste0("  x1: ", x1), file=ymlf, sep="\n", append=TRUE)
	cat(paste0("  y1: ", y1), file=ymlf, sep="\n", append=TRUE)
	cat(paste0("  x2: ", x2), file=ymlf, sep="\n", append=TRUE)
	cat(paste0("  y2: ", y2), file=ymlf, sep="\n", append=TRUE)
	cat(paste0("  x3: ", x3), file=ymlf, sep="\n", append=TRUE)
	cat(paste0("  y3: ", y3), file=ymlf, sep="\n", append=TRUE)
	cat(paste0("  fov: ", fov), file=ymlf, sep="\n", append=TRUE)
	cat(paste0("  wall: ", wall), file=ymlf, sep="\n", append=TRUE)
	cat(paste0("  length: ", length), file=ymlf, sep="\n", append=TRUE)
}

# compute horizontal angle
angle <- atan2(y2-y1, x2-x1)

# rotate reference points
image <- load.image(images[start])
width <- dim(image)[1]
height <- dim(image)[2]
points.rotate <- function(x,y) list(x=width/2+(x-width/2)*cos(-angle)-(y-height/2)*sin(-angle),y=height/2+(x-width/2)*sin(-angle)+(y-height/2)*cos(-angle))
r <- points.rotate(x1, y1)
x1 <- r$x[1]
y1 <- r$y[1]
r <- points.rotate(x2, y2)
x2 <- r$x[1]
y2 <- r$y[1]
r <- points.rotate(x3, y3)
x3 <- r$x[1]
y3 <- r$y[1]
scale <- ( x2 - x1 ) / length

# analysis windows
horiz1 <- y3 - 5 - 20
horiz2 <- y3 + 5
x1 <- x1 - 20
x2 <- x2 + 20

# geometry wall thickness from visible angle
dist <- (width/2/scale) / tan(fov/2*pi/180)
hdist <- abs(y3-height/2)/scale
circ1 <- Circle$new(c(0,dist+diameter/2), diameter/2)
circ2 <- Circle$new(c(0,dist+diameter/2), diameter/2-wall)
line <- Line$new(c(0,0), c(hdist, dist))
point1 <- intersectionCircleLine(circ1, line)[[1]]
point2 <- intersectionCircleLine(circ1, line)[[2]]
if(point2[2] < point1[2]){ point1 <- point2 }
point2 <- intersectionCircleLine(circ2, line)[[1]]
point3 <- intersectionCircleLine(circ2, line)[[2]]
if(point3[2] < point2[2]){ point2 <- point3 }
wall <- sqrt(sum((point1 - point2)^2))
print(paste0('Geometric wall thickness: ', wall, ' cm'))

# correct apparent wall thickness
wall <- wall / rindex
print(paste0('Optic wall thickness: ', wall, ' cm'))

# empty results
resultx <- c()
resulty <- c()
resultz <- c()

#############################

# iterate images
#for (time in c(start:length(images)))
for (time in c(start))
{
	# rotate image
	image <- load.image(images[time])
	map.rotate <- function(x,y) list(x=width/2+(x-width/2)*cos(angle)-(y-height/2)*sin(angle),y=height/2+(x-width/2)*sin(angle)+(y-height/2)*cos(angle))
	im <- imwarp(image,map=map.rotate, interpolation = "nearest", direction="backward", boundary = "neumann")
	if(plotImages == TRUE)
	{
		plot(im, axes=FALSE)
	}

	# draw analysis windows
	if(plotImages == TRUE)
	{
		lines(c(x1,x2,x2,x1,x1), c(horiz1,horiz1,horiz1+20,horiz1+20,horiz1), col = rgb(red = 0, green = 1, blue = 0, alpha = 1.0))
		lines(c(x1,x2,x2,x1,x1), c(horiz2,horiz2,horiz2+20,horiz2+20,horiz2), col = rgb(red = 1, green = 0, blue = 0, alpha = 1.0))
	}

	col1 <- rep(0, x2-x1+1)
	for (val in rep(1:20)) {
		col1 <- col1 + imrow(R(im),horiz1+val-1)[rep(x1:x2)] # red channel
	}
	col1 <- col1 / 20.0

	col2 <- rep(0, x2-x1+1)
	for (val in rep(1:20)) {
		col2 <- col2 + imrow(R(im),horiz2+val-1)[rep(x1:x2)] # red channel
	}
	col2 <- col2 / 20.0

	if(plotImages == TRUE)
	{
		lines(rep(x1:x2), horiz1-20-col1*200, col = rgb(red = 0, green = 1, blue = 0, alpha = 1.0))
		lines(rep(x1:x2), horiz2+20+20+max(col2)*200-col2*200, col = rgb(red = 1, green = 0, blue = 0, alpha = 1.0))
	}

	#############################
	# water curve analysis

	# locate the step
	q20 <- quantile(col1, c(0.20))[[1]]
	q80 <- quantile(col1, c(0.80))[[1]]
	level = 0.5 * q20 + 0.5 * q80

	# Search from right to left
	for (i in rep(length(col1):1)) {
		if(col1[i] < level)
		{
			break
		}
	}

	# refine and search again
	q10 <- quantile(col1[rep((i-100):(i+100))], c(0.10))[[1]]
	q90 <- quantile(col1[rep((i-100):(i+100))], c(0.90))[[1]]
	level =  0.9 * q10 + 0.1 * q90

	# this time from left to right
	for (j in rep((i-100):(i+100)))
	{
		if(col1[j] > level)
		{
			break
		}
	}
	water <- j

	# plot water level
	if(plotImages == TRUE)
	{
		lines(c(x1+j-1,x1+j-1), c(horiz1+1, horiz1+19), col = rgb(red = 1, green = 1, blue = 0, alpha = 1.0))
	}

	#############################
	# tape analysis

	# locate the marks as minimums
	div <- length(col2)/12
	min <- c()
	min[1] <- which.min(col2[rep(                 1:as.integer( 1*div))])
	min[2] <- which.min(col2[rep(as.integer( 1*div):as.integer( 3*div))]) + as.integer( 1*div) - 1
	min[3] <- which.min(col2[rep(as.integer( 3*div):as.integer( 5*div))]) + as.integer( 3*div) - 1
	min[4] <- which.min(col2[rep(as.integer( 5*div):as.integer( 7*div))]) + as.integer( 5*div) - 1
	min[5] <- which.min(col2[rep(as.integer( 7*div):as.integer( 9*div))]) + as.integer( 7*div) - 1
	min[6] <- which.min(col2[rep(as.integer( 9*div):as.integer(11*div))]) + as.integer( 9*div) - 1
	min[7] <- which.min(col2[rep(as.integer(11*div):length(col2))])       + as.integer(11*div) - 1

	for (j in 1:7)
	{
		# center the mark
		for (k in 1:10)
		{
			range <- rep((min[j]-10):(min[j]+10))
			if( sum(col2[range]*range)/sum(col2[range]) > min[j])
			{
				min[j] <- min[j] - 1
			}
			else
			{
				min[j] <- min[j] + 1
			}
		}

		# fit to a quadratic and compute its minimum
		range <- rep((min[j]-10):(min[j]+10))
		d <- data.frame(x=range, y=col2[range])
		model <- lm(y~x+I(x^2), data=d);
		min[j] <- -model$coefficients[[2]]/2/model$coefficients[[3]]
	}

	# plot tape marks
	if(plotImages == TRUE)
	{
		lines(c(x1+min[1]-1,x1+min[1]-1), c(horiz2+1, horiz2+19), col = rgb(red = 1, green = 1, blue = 0, alpha = 1.0))
		lines(c(x1+min[2]-1,x1+min[2]-1), c(horiz2+1, horiz2+19), col = rgb(red = 1, green = 1, blue = 0, alpha = 1.0))
		lines(c(x1+min[3]-1,x1+min[3]-1), c(horiz2+1, horiz2+19), col = rgb(red = 1, green = 1, blue = 0, alpha = 1.0))
		lines(c(x1+min[4]-1,x1+min[4]-1), c(horiz2+1, horiz2+19), col = rgb(red = 1, green = 1, blue = 0, alpha = 1.0))
		lines(c(x1+min[5]-1,x1+min[5]-1), c(horiz2+1, horiz2+19), col = rgb(red = 1, green = 1, blue = 0, alpha = 1.0))
		lines(c(x1+min[6]-1,x1+min[6]-1), c(horiz2+1, horiz2+19), col = rgb(red = 1, green = 1, blue = 0, alpha = 1.0))
		lines(c(x1+min[7]-1,x1+min[7]-1), c(horiz2+1, horiz2+19), col = rgb(red = 1, green = 1, blue = 0, alpha = 1.0))
	}

	# compute linear regression
	d <- data.frame(x=min, y=c(46, 45, 44, 43, 42, 41, 40))
	model <- lm(y~x, data=d)
	ca <- model$coefficients[[1]]
	cb <- model$coefficients[[2]]

	# apply regression to the water level
	watercm <- water * cb + ca

	# correct water level and apply regression
	water2 <- (x1+water-1)-wall*((width-(x1+water-1))/scale-width/scale*0.5)/(width/scale/tan(fov/2*pi/180)/2)*scale
	water2 <- water2 - x1 + 1
	water2cm <- water2 * cb + ca

	# store results
	resultx <- c(resultx, time)
	resulty <- c(resulty, watercm)
	resultz <- c(resultz, water2cm)
}

data.frame(x=resultx, y=resulty, z=resultz)