#! /bin/csh -f
#  $Id: matchshifts,v bde41caef8de 2009/01/22 06:07:01 mast $
#  Log at end of file
#

nohup
set SIZELIMIT = 1600
set MAXERROR = 3.0

setenv LC_NUMERIC C

if ($#argv == 0)then
  echo 'Usage: matchshifts root_name1 root_name2 NX NY NZ [solve_xf_in] [solve_xf_out]'
  echo ' root_name1 = root name of tomogram being matched TO (omit .rec)'
  echo ' root_name2 = root name of tomogram being transformed (omit .rec)'
  echo ' NX NY NZ = size of volume to correlate'
  echo ' solve_xf_in = file with 3x3 transform only (default solvezero.xf)'
  echo ' solve_xf_out = output file combining transform and shifts (default solve.xf)'
  exit 0
endif

set tmpdir = /usr/tmp
if ($?IMOD_DIR) then
    setenv PATH "$IMOD_DIR/bin:$PATH"
    if (-e "$IMOD_DIR/bin/settmpdir") source "$IMOD_DIR/bin/settmpdir"
endif

set exitstatus = 1
onintr clean

if ($#argv < 5 || $#argv >7) then
  echo 'wrong # of arguments'
  exit 1
endif

set infile = solvezero.xf
set outfile = solve.xf
set recfile1 = $argv[1].rec
set recfile2 = $argv[2].rec
set matfile = "$tmpdir/mat2.$$"
set pad1file = "$tmpdir/pad1.$$"
set pad2file = "$tmpdir/pad2.$$"
set cortmp = "$tmpdir/corr3d.$$"
set proj1file = "$tmpdir/proj1.$$"
set proj2file = "$tmpdir/proj2.$$"
set projst = "$tmpdir/projst.$$"
set rotxf = "$tmpdir/rotxf.$$"
set txcxf = "$tmpdir/txcxf.$$"
set midfile = "$tmpdir/midxf.$$"
set logfile = "$tmpdir/logtmp.$$"
set checkfile1 = matchcheck.rec
set checkfile2 = matchcheck.mat

if($#argv >5 ) then
  set infile = $argv[6]
endif
if($#argv >6 ) then
  set outfile = $argv[7]
endif

# Detect non-zero shifts in input file and just copy the data
#
set deltas = `awk '{print 1000. * $4}' $infile`
@ deltax = $deltas[1]
@ deltay = $deltas[1]
@ deltaz = $deltas[1]
if ($deltax != 0 || $deltay != 0 || $deltaz != 0) then
    if (-e $outfile) \mv -f $outfile $outfile~
    \cp $infile $outfile
    chmod u+rw $outfile
    echo "MATCHSHIFTS: Found non-zero shifts in $infile, copied it to $outfile"
    exit 0
endif

set nx = $argv[3]
set ny = $argv[4]
@ nz = $argv[5]
set nxyz =  `header -si $recfile1 | sed '/[[:cntrl:]]/s///g'`

@ nxuse = $nxyz[1]
if ($nxuse > $SIZELIMIT) @ nxuse = $SIZELIMIT
@ ix0 = ($nxyz[1] - $nxuse) / 2
@ ix1 = ($nxyz[1] + $nxuse) / 2 - 1
@ nzuse = $nxyz[3]
if ($nzuse > $SIZELIMIT) @ nzuse = $SIZELIMIT
@ iz0 = ($nxyz[3] - $nzuse) / 2
@ iz1 = ($nxyz[3] + $nzuse) / 2 - 1

echo " "
echo "MATCHSHIFTS: PROJECTING $recfile1"

xyzproj <<EOF >! "$logfile"
$recfile1
$proj1file
$ix0,$ix1,,,$iz0,$iz1
Z
0,0,0
/
/
/
/
EOF

if ($status) goto error

echo "MATCHSHIFTS: PROJECTING $recfile2"

xyzproj <<EOF  >! "$logfile"
$recfile2
$proj2file
$iz0,$iz1,,,$ix0,$ix1
Z
0,0,0
/
/
/
/
EOF

if ($status) goto error

set sol = `cat $infile`
\echo "1 0 0 1 0 0" > "$rotxf"
\echo "$sol[1] $sol[3] $sol[9] $sol[11] 0 0" >> "$rotxf"

echo "MATCHSHIFTS: STACKING PROJECTIONS WITH ROTATION"

newstack <<EOF >! "$logfile"
2
$proj1file
0
$proj2file
0
1
$projst
/
/
0
1
$rotxf
0,1
0
EOF

if ($status) goto error

@ nxin = $nxyz[1]
@ nyin = $nxyz[2]
@ xtrim = $nxuse / 10
@ ytrim = $nzuse / 10

echo "MATCHSHIFTS: CORRELATING WITH TILTXCORR"
echo " "

tiltxcorr <<EOF
$projst

$txcxf
1
0,0
0
/
0
$xtrim $ytrim
/
/
/
EOF

if ($status) goto error

set txcshift = `tail -n 1 "$txcxf"`

\echo "$sol[1]   $sol[2]   $sol[3]   $txcshift[5]"  > "$midfile"
\echo "$sol[5]   $sol[6]   $sol[7]   0"  >> "$midfile"
\echo "$sol[9]   $sol[10]   $sol[11]   $txcshift[6]"  >> "$midfile"


@ nzt = $nxyz[3]
@ izlo = ($nzt / 2) - ($nz / 2)
@ izhi = $izlo + $nz - 1

echo " "
echo "MATCHSHIFTS: EXTRACTING SMALL MATCHING VOLUME FROM $recfile2"

matchvol -StandardInput <<HERESTRING >! "$logfile"
InputFile $recfile2
OutputFile $matfile
TemporaryDirectory $tmpdir
OutputSizeXYZ $nx $ny $nz
TransformFile $midfile
InterpolationOrder 1
HERESTRING

if ($status) goto error

@ nxtap = ($nx + 9) / 10
@ nytap = ($ny + 9) / 10
@ nztap = ($nz + 9) / 10
@ nxpad = ($nx + 9) / 10
@ nypad = ($ny + 4) / 5
@ nzpad = ($nz + 9) / 10
@ nxcen = $nx / 2
@ nycen = $ny / 2
@ nzcen = $nz / 2

echo "MATCHSHIFTS: GETTING TAPERED, PADDED VOLUMES FOR CORRELATING"

tapervoledge <<HERESTRING >! "$logfile"
$matfile
$pad2file
$nx,$ny,$nz
$nxcen,$nycen,$nzcen
$nxpad,$nypad,$nzpad
$nxtap,$nytap,$nztap
HERESTRING

if ($status) goto error

@ nxcen = $nxyz[1] / 2
@ nycen = $nxyz[2] / 2
@ nzcen = $nxyz[3] / 2

tapervoledge <<HERESTRING >! "$logfile"
$recfile1
$pad1file
$nx,$ny,$nz
$nxcen,$nycen,$nzcen
$nxpad,$nypad,$nzpad
$nxtap,$nytap,$nztap
HERESTRING

if ($status) goto error


#clip resize -ox $nx -oy $ny -oz $nz $recfile1 $clipfile
#newst -se $izlo-$izhi -si $nx,$ny $recfile1 $clipfile

echo "MATCHSHIFTS: CORRELATING VOLUMES TO GET 3-D SHIFTS"

#invert order of files to get amount to shift mat file to match clip file

set displ = `clip corr -3d -n 0 "$pad2file" "$pad1file" "$cortmp" | grep '^(' | sed '/.*(\(.*\)).*/s//\1/'`

set xshift = `echo "$txcshift[5] $displ[1]" | awk '{print $1 + $2}'`
set yshift = `echo "$txcshift[6] $displ[3]" | awk '{print $1 + $2}'`
if (-e $outfile) \mv -f $outfile $outfile~
\echo "$sol[1]   $sol[2]   $sol[3]   $xshift"  > $outfile
\echo "$sol[5]   $sol[6]   $sol[7]   $displ[2]"  >> $outfile
\echo "$sol[9]   $sol[10]   $sol[11]   $yshift"  >> $outfile
chmod u+rw $outfile

echo "MATCHSHIFTS: GETTING VOLUMES TO CHECK THE MATCH"

@ nxextract = $nx * 3

matchvol -StandardInput <<HERESTRING >! "$logfile"
InputFile $recfile2
OutputFile $matfile
TemporaryDirectory $tmpdir
OutputSizeXYZ $nxextract $ny $nz
TransformFile $outfile
InterpolationOrder 1
HERESTRING

if ($status) goto error

@ nxcen = $nxcen - $nx
tapervoledge <<HERESTRING >! "$logfile"
$recfile1
$checkfile1
$nx,$ny,$nz
$nxcen,$nycen,$nzcen
0,0,0
$nxtap,$nytap,$nztap
HERESTRING

if ($status) goto error

@ nxcen = $nx / 2
@ nycen = $ny / 2
@ nzcen = $nz / 2

tapervoledge <<HERESTRING >! "$logfile"
$matfile
$checkfile2
$nx,$ny,$nz
$nxcen,$nycen,$nzcen
0,0,0
$nxtap,$nytap,$nztap
HERESTRING

if ($status) goto error

echo "MATCHSHIFTS: CORRELATING MATCH CHECK VOLUMES"

set displ = `clip corr -3d -n 0 "$checkfile2" "$checkfile1" "$cortmp" | grep '^(' | sed '/.*(\(.*\)).*/s//\1/'`

set checkerror = `echo $displ | awk '{printf "%.2f", sqrt($1 * $1 + $2 * $2 + $3 * $3)}'`

echo "Displacement between check volumes is $displ in X Y Z, total $checkerror"

@ errdiff = `echo "$checkerror $MAXERROR" | awk '{print int(100. * ($1 - $2))}'`
if ($errdiff > 0) then
    echo "ERROR: MATCHSHIFTS - measured displacement between check volumes is $checkerror"
    echo "To see if this is a problem, run:  3dmod -Y $checkfile1 $checkfile2"
    echo In eTomo, press \"View Match Check Volumes\"
    goto clean
endif

set exitstatus = 0
goto clean

error:
if (-e "$logfile") grep ERROR "$logfile"

clean:
\rm -f "$cortmp" "$cortmp~" "$proj1file" "$proj2file" "$projst" "$rotxf" "$txcxf" "$midfile"
\rm -f "$pad1file"* "$pad2file" "$matfile" "$pad2file~" "$matfile~" "$logfile"
exit $exitstatus

#
#  $Log$
#  Revision 3.17  2008/01/04 07:03:38  mast
#  Fixed for space in IMOD_DIR, possibly for space in tmpdir
#
#  Revision 3.16  2006/08/17 18:36:55  mast
#  Added locale statement to keep output of awk from having instead of .
#
#  Revision 3.15  2006/02/16 06:47:09  mast
#  Stripped control chars from output of sed/header etc for Windows
#
#  Revision 3.14  2005/11/19 04:31:26  mast
#  Quote path setting to preserve spaces
#
#  Revision 3.13  2004/08/27 05:46:48  mast
#  Switched to using head -n and tail -n
#
#  Revision 3.12  2004/07/08 22:54:18  mast
#  Switched some echos to \echo to dodge Cygwin tcsh bug
#
#  Revision 3.11  2004/06/29 03:48:54  mast
#  Added nohup
#
#  Revision 3.10  2004/06/15 23:18:05  mast
#  Fix the message that it gives when copying
#
#  Revision 3.9  2004/06/10 05:40:24  mast
#  Made it copy solvezero.xf to solve.xf when it already has shifts, so it can
#  be called in all cases after running solvematch
#
#  Revision 3.8  2003/10/24 02:14:00  mast
#  source settmpdir to get tmpdir, better for Windows
#
#  Revision 3.7  2003/10/16 20:43:34  mast
#  Changed error message upon check volume mismatch to be more helpful
#
#  Revision 3.6  2003/10/11 00:04:21  mast
#  Switched to new input method for matchvol and used linear interpolation
#
#  Revision 3.5  2003/06/20 20:04:59  mast
#  Let tmpdir be defined and fall back to /tmp if /usr/tmp does not exist
#
#  Revision 3.4  2003/03/21 22:25:52  mast
#  Fix error test and output when displacement is too large
#
#  Revision 3.3  2003/03/14 02:11:10  mast
#  Made it produce two checking volumes that the user can inspect, and stopped
#  when displacement between those volumes was too high.  Made it report last
#  error statement and exit with error when a program has an error
#
#  Revision 3.2  2002/07/21 19:34:05  mast
#  Suppressed output of all programs except Tiltxcorr
#
#  Revision 3.1  2002/07/21 05:52:35  mast
#  Decreased padding in X and Z to compensate for increased sizes in Y
