#!/bin/bash
#
#
# Plot the results from ctffind using gnuplot
#
# Alexis Rohou, June 2014
#
# Copyright 2014 Howard Hughes Medical Institute. All rights reserved.
# Use is subject to Janelia Farm Research Campus Software Copyright 1.1
# license terms ( http://license.janelia.org/license/jfrc_copyright_1_1.html )
#
# Last modification: Tapu Shaik, April 2018
# Alexis Rohou, Dec 2018 (adding score to plot)



# Parse arguments
if [ $# -ne 1 ]
then
	echo "Usage: `basename $0` output_from_ctffind_avrot.txt"
	exit 65
fi

input_fn=$1
output_fn=${1%%.???}.pdf
#input_summary_fn=$(ls ${1%%_avrot.???}.${1##*.})
input_summary_fn="${1%%_avrot.???}.${1##*.}"

# Check whether gnuplot is available
if ! hash gnuplot 2>/dev/null; then
	echo "Gnuplot was not found on your system. Please make sure it is installed."
	exit -1
fi

# Check whether gnuplot >= 4.6 is available
gnuplot_version=$(  gnuplot --version | awk '{print $2}'  )
compare_result=$( echo "$gnuplot_version >= 4.6" | bc )
if [ $compare_result -le 0 ]; then
	echo "This script requires gnuplot version >= 4.6, but you have version $gnuplot_version"
	exit -1
fi

# Define a neat function
function transpose_ignore_comments {
gawk '
{
	if ($1 ~ /^#/) {
		print $0
	} else {
		for (i=1; i<=NF; i++)  {
       		a[NR,i] = $i
    	}
	}
}
NF>p { p = NF }
END {
    for(j=1; j<=p; j++) {
        str=a[1,j]
        for(i=2; i<=NR; i++){
            str=str" "a[i,j];
        }
        print str
    }
}' $1
}

# CTFFind outputs data in lines, but gnuplot wants things in columns
transpose_ignore_comments $input_fn > /tmp/tmp.txt

# Let's grab useful values
pixel_size=$(gawk 'match($0,/Pixel size: ([0-9.]*)/,a) {print a[1]}' $input_fn)
number_of_micrographs=$(gawk 'match($0,/Number of micrographs: ([0-9]*)/,a) {print a[1]}' $input_fn)
mic_name=$(grep "Input file:" $input_fn | awk 'BEGIN {FS="Input file: "} {print $2}' | cut -d " " -f1  | sed 's/_/\\\_/g')
lines_per_micrograph=6

if [ -f $input_summary_fn ]; then
# Let's grab values from the summary file
i=0
while read -a myArray
do
	if [[ ${myArray[0]} != \#* ]]; then
		df_one[++i]=${myArray[1]}
		df_two[i]=${myArray[2]}
		angast[i]=${myArray[3]}
		pshift[i]=${myArray[4]}
		score[i]=${myArray[5]}
		maxres[i]=${myArray[6]}
	fi
done < $input_summary_fn
else
df_one[0]="0"
df_two[0]="0"
angast[0]="0"
pshift[0]="0"
score[0]="0"
maxres[0]="0"
fi

# Run Gnuplot
gnuplot > /dev/null 2>&1  <<EOF
#cat <<EOF > temp.txt
set border linewidth 1.5
if (strstrt(GPVAL_TERMINALS, 'pdfcairo') > 0) {
	set terminal pdfcairo size 10.0cm,4.0cm enhanced font 'Arial,9'
	set output '$output_fn'
} else {
	set terminal postscript size 10,4 enhanced font 'Arial,14'
	set output '${output_fn%%.pdf}.ps'
}

# color definitions
set style line 1 lc rgb '#0060ad' lt 1 lw 2 pt 7 # --- blue
set style line 2 lc rgb 'red'     lt 1 lw 1 pt 7 fill transparent solid 0.5 # --- red
set style line 3 lc rgb 'orange'  lt 3 lw 1 pt 7 # --- orange
set style line 4 lc rgb 'light-blue' lt 1 lw 2 pt 7 # --- light blue
set style line 5 lc rgb 'green'   lt 1 lw 2 pt 7
set style line 6 lc rgb 'gray'    lt 1 lw 2 pt 7

set xlabel 'Spatial frequency(1/Å)'
set ylabel 'Amplitude (or cross-correlation)'
set yrange [-0.1:1.1]
set key outside

defocus_1_values="${df_one[*]}"
defocus_2_values="${df_two[*]}"
angast_values="${angast[*]}"
pshift_values="${pshift[*]}"
score_values="${score[*]}"
maxres_values="${maxres[*]}"

do for [current_micrograph=1:$number_of_micrographs] {
	def_1=sprintf('%.0f Å',word(defocus_1_values,current_micrograph)+0)
	def_2=sprintf('%.0f Å',word(defocus_2_values,current_micrograph)+0)
	angast=sprintf('%.1f °',word(angast_values,current_micrograph)+0)
	pshift=sprintf('%.2f rad',word(pshift_values,current_micrograph)+0)
	score=sprintf('%.3f',word(score_values,current_micrograph)+0)
	maxres=sprintf('%.2f Å',word(maxres_values,current_micrograph)+0)
#	set title 'Micrograph '.current_micrograph." of $number_of_micrographs\n{/*0.5 Defocus 1: ".def_1.' | Defocus 2: '.def_2.' | Azimuth: '.angast.' | Phase shift: '.pshift.' | Score: '.score.'}'
	set title '$mic_name, '.current_micrograph." of $number_of_micrographs micrographs\n{/*0.5 Defocus 1: ".def_1.' | Defocus 2: '.def_2.' | Azimuth: '.angast.' | Phase shift: '.pshift.' | Score: '.score.' | MaxRes: '.maxres.'}'
	plot '/tmp/tmp.txt' using (\$1):(column(4+(current_micrograph-1)*$lines_per_micrograph)) w lines ls 3 title 'CTF fit', \
		 ''             using (\$1):(column(5+(current_micrograph-1)*$lines_per_micrograph)) w lines ls 1 title 'Quality of fit', \
	 	 ''             using (\$1):(column(3+(current_micrograph-1)*$lines_per_micrograph))  w lines ls 5 title 'Amplitude spectrum'
}
EOF

# Convert the postscript file to a pdf file if needed
if [ -f ${output_fn%%.pdf}.ps ]; then
	ps2pdf ${output_fn%%.pdf}.ps $output_fn
fi

# Fix the title of the PDF file
if hash pdftk 2>/dev/null; then
	echo "InfoKey: Title" > /tmp/pdftemp.txt
	echo "InfoValue: $1" >> /tmp/pdftemp.txt
	mv ${output_fn} /tmp/tmp.pdf
	pdftk /tmp/tmp.pdf update_info /tmp/pdftemp.txt output $output_fn
elif hash exiftool 2>/dev/null; then
	exiftool -Title="$1" ${output_fn}
else
	echo "Neither pdftk nor exiftool were found on your system. Please install one of them to get improved metadata in your output PDF file."
fi

# Make a PNG version
#convert -density 300 -flatten ${output_fn} ${output_fn%%.pdf}.png

# Let the user know where to find results
echo " "
echo "Generated plots of 1D profiles and fits: ${output_fn}"
echo " "
