#!/usr/bin/Rscript --silent
# This script is a splitted second step from frame_proc - the columns with ratios
rem <- function(...) invisible(T)
rem( '
Rscript.exe "%~F0" "%1" "%2" "%3" "%4" "%~F5" "%~F6"
EXIT /B
rem ')
### above=BAT, below=R
library(optparse)

option_list <- list(
	make_option(c("-v", "--verbose"), action="store_true", default=FALSE,
	help="Print extra output [default=false]"),
	make_option(c("-r", "--realsize"), type="double", default="1.0",
	help="real size factor [default %default]",
	metavar="number"),
	make_option(c("-m", "--time"), type="character", default="",
	help="minutes from start [default %default]",
	metavar="number")
)
opt_parser <- OptionParser(usage = "usage: %prog [options]", option_list = option_list, description = "Post-proc object table", epilogue = "Send feedback to mackoel@gmail.com")
flaggs <- parse_args(opt_parser, args = commandArgs(trailingOnly = TRUE), print_help_and_exit = TRUE, positional_arguments = TRUE)
opts <- flaggs$options
args <- flaggs$args
# Input/output
infile <- args[1]
outfile <- args[2]
# Read data
qdata <- read.csv(infile)
# Fix column names
colnames(qdata) <- sub("eea1", "green", colnames(qdata))
colnames(qdata) <- sub("egfr", "red", colnames(qdata))
colnames(qdata) <- sub("cort", "red", colnames(qdata))
# Fix number of cells
if ("cell" %in% colnames(qdata)) {
	cat (infile, "cell column present\n")
	OPTIMALK <- max(qdata$cell)
} else {
	cat (infile, "cell column absent\n")
	OPTIMALK <- 1
	qdata <- cbind(qdata, cell=rep(1, times=dim(qdata)[1]))
}
# Calculate ratios
coloc_red <- qdata$overlap_mask_notempty/qdata$red_mask_notempty
coloc_green <- qdata$overlap_mask_notempty/qdata$green_mask_notempty
coloc_vol <- qdata$overlap_mask_notempty/qdata$m000
coloc_cvol <- qdata$red_mask_notempty/qdata$m000
coloc_evol <- qdata$green_mask_notempty/qdata$m000
# Calculate real sizes
real_area_obj <- qdata$m000 * opts$realsize * opts$realsize
real_radius_obj <- sqrt(qdata$m000/pi) * opts$realsize
real_area_red <- qdata$red_mask_notempty * opts$realsize * opts$realsize
real_radius_red <- sqrt(qdata$red_mask_notempty/pi) * opts$realsize
real_area_green <- qdata$green_mask_notempty * opts$realsize * opts$realsize
real_radius_green <- sqrt(qdata$green_mask_notempty/pi) * opts$realsize
real_area_overlap <- qdata$overlap_mask_notempty * opts$realsize * opts$realsize
real_radius_overlap <- sqrt(qdata$overlap_mask_notempty/pi) * opts$realsize
# Re-Calculate means over channel masks
green_mean <- qdata$green_mask_muc / qdata$green_mask_notempty
red_mean <- qdata$red_mask_muc / qdata$red_mask_notempty
qdata[qdata$green_mask_notempty > 0,]$green_mask_mean <- green_mean[qdata$green_mask_notempty > 0]
qdata[qdata$red_mask_notempty > 0,]$red_mask_mean <- red_mean[qdata$red_mask_notempty > 0]
# Bind new columns
qdata <- cbind(qdata, coloc_red, coloc_green, coloc_vol, coloc_cvol, coloc_evol, real_area_obj, real_radius_obj, real_area_red, real_radius_red, real_area_green, real_radius_green, real_area_overlap, real_radius_overlap)
write.csv(file=outfile, qdata, row.names=F)

