#!/bin/tcsh -f

# tomrec
# A script to do a tomographic reconstruction:
#	1.	bmgft: Generate Fourier transforms from all micrographs
#	2.	btomrec: Calculate slabs in Fourier space and back-transform slices
#	3.	bzfft: Assemble slabs and back-transform z-lines
# Usage:
#	tomrec -reconstruction <file> -resolution <angstrom> -size <x,y,z>
#			-thickness <pixels> -scale <number> -rescale <avg,std>
#			-remove -output <file> <file> [<file> ...]
#	Requirement: At least the size must be specified.
#	Note: The disk space requirements are very large - typically 40 times the size of the tomogram!
#	Distributed processing: The lines with "psubmit" can be used with the Peach distributed system.
# Examples:
#	tomrec -rec map.pif -resol 45 -size 1024,1024,120 -thick 20 -scale 1 -rescale 127,25 -remove -out out.star input.star
# Author: Bernard Heymann
# 20060106 - 20071114

# Defaults
set REC = "rec.pif"
set OUTSTAR = "out.star"
set RES = 30
set X = 256
set Y = $X
set Z = 200
set TH = 20
set SC = 1
set RESCALE = 127,20
set REMOVE = " "

# Recreate command line
echo -n tomrec
set a = 1
while ( $a <= $#argv )
	echo -n " $argv[$a]"
	@ a++
end
echo " "

if ( $#argv < 1 ) then
	echo "No input!"
	exit
endif

# Interpret options
set a = 1
while ( "x$argv[$a]" =~ x-* )
#	echo $a $argv[$a]
	if ( "x$argv[$a]" =~ x-rec* ) then
		@ a++
		set REC = $argv[$a]
	endif
	if ( "x$argv[$a]" =~ x-out* ) then
		@ a++
		set OUTSTAR = $argv[$a]
	endif
	if ( "x$argv[$a]" =~ x-reso* ) then
		@ a++
		set RES = $argv[$a]
	endif
	if ( "x$argv[$a]" =~ x-siz* ) then
		@ a++
		set SIZE = $argv[$a]
		set X = `echo $SIZE | cut -f1 -d","`
		set Y = `echo $SIZE | cut -f2 -d","`
		set Z = `echo $SIZE | cut -f3 -d","`
	endif
	if ( "x$argv[$a]" =~ x-thi* ) then
		@ a++
		set TH = $argv[$a]
	endif
	if ( "x$argv[$a]" =~ x-sca* ) then
		@ a++
		set SC = $argv[$a]
	endif
	if ( "x$argv[$a]" =~ x-resc* ) then
		@ a++
		set RESCALE = $argv[$a]
	endif
	if ( "x$argv[$a]" =~ x-rem* ) then
		set REMOVE = $argv[$a]
	endif
	@ a++
end

set INSTAR = $argv[$a]

set SF = "ft.star"
set BASE = "slab"
set DOT = "-trans slices"

rm -rf ${BASE}????.sup

set NMG = `btomo -v 7 $INSTAR | awk '$1=="Micrographs:" { print $2 }'`
set NSEL = `btomo -v 7 $INSTAR | awk '$1=="Micrographs:" { print $4 }'`

echo "Reconstruction parameters:"
echo "-------------------------"
date
echo "Micrographs    = $NMG"
echo "Selected       = $NSEL"
echo "Base           = $BASE"
echo "Size           = ${X},${Y},${Z}"
echo "Scale          = $SC"
echo "Slab           = $TH"
echo "Resolution     = $RES"
echo "Rescale        = $RESCALE"
echo "Output STAR    = $OUTSTAR"
echo "Reconstruction = $REC"
echo ""

# Fourier transform all the micrographs
set N = 0
while ( $N < $NMG )
	set NN = `printf "%03d" $N`
	set JOB = "ft$NN"
	set SFNAME = `echo $SF | cut -f1 -d"."`"_$NN.star"
	if ( ! -e $SFNAME ) then
		echo bmgft -v 1 -select $N -size ${X},${Y},${Z} -scale $SC $REMOVE -out $SFNAME $INSTAR
#		psubmit bmgft -v 1 -select $N -size $X -scale $SC $REMOVE -out $SFNAME $INSTAR -jn $JOB
		bmgft -v 1 -select $N -size ${X},${Y},${Z} -scale $SC $REMOVE -out $SFNAME $INSTAR
	endif
	@ N++
end

set SFNAME = `echo $SF | cut -f1 -d"."`

# Check if all the micrographs have been transformed
set NTR = 0
while ( $NTR < $NSEL )
	sleep 10
	set NTR = `ls -1 ${SFNAME}_???.star | wc -l`
	echo "$NTR micrographs transformed"
end

# Re-integrate all STAR files
echo bmg -v 7 -in $INSTAR -out $SF ${SFNAME}_???.star
bmg -v 7 -in $INSTAR -out $SF ${SFNAME}_???.star

#exit

# Slab reconstructions
set SL = 0
set NSL = 0
while ( $SL < $Z )
	set JOB = ${BASE}`printf %04d ${SL}`
	set FILE = ${JOB}".sup"
	set SFILE = ${JOB}".star"
	@ SLE = $SL + $TH - 1
	if ( ! -e $FILE ) then
		echo btomrec -v 1 -slab ${SL},${SLE} $DOT -scale $SC -rec $FILE -resol $RES -size ${X},${Y},${Z} -out $SFILE $SF
#		psubmit btomrec -v 1 -slab ${SL},${SLE} $DOT -scale $SC -rec $FILE -resol $RES -size ${X},${Y},${Z} -out $SFILE $SF -jn $JOB
		btomrec -v 1 -slab ${SL},${SLE} $DOT -scale $SC -rec $FILE -resol $RES -size ${X},${Y},${Z} -out $SFILE $SF
	endif
	@ SL = $SL + $TH
	@ NSL++
end

# Check if all the slabs have been generated
set NSLG = 0
while ( $NSLG < $NSL )
	sleep 10
	set NSLG = `ls -1 ${BASE}????.sup | wc -l`
	echo "$NSLG slabs generated"
end

# Re-assemble and do z-line back-transforms
bzfft -v 7 -dat byte -rescale $RESCALE -rec $REC -out $OUTSTAR ${BASE}????.star