#! /bin/csh -f
# patchcorr3d - to get local 3d cross-correlations between two volumes at
# a grid of positions.
if ($#argv != 14 && $#argv != 16) then
cat <<EOCAT
Usage: patchcorr3d  XSIZE YSIZE ZSIZE  NX NY NZ  XLO XHI YLO YHI ZLO ZHI\
       FILEA FILEB [TRANSFORM_FILE  ORIGINAL_FILEB]
Finds the displacement of the volume in FILEB relative to the volume in FILEA,
 in a series of local volumes of size XSIZE by YSIZE by ZSIZE,
 at grid positions determined by:
   NX, NY, NZ: number of positions in X, Y and Z
   XLO, XHI, YLO, YHI, ZLO, ZHI: minimum and maximum coordinates of the block
       from which local volumes will be extracted 
EOCAT
exit 1
endif

setenv LC_NUMERIC C

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

@ nxl = $1; shift
@ nyl = $1; shift
@ nzl = $1; shift
@ nx = $1; shift
@ ny = $1; shift
@ nz = $1; shift
@ xbl = $1 - 1; shift
@ xhi = $1; shift
@ ybl = $1 - 1; shift
@ yhi = $1; shift
@ zbl = $1 - 1; shift
@ zhi = $1; shift
#
set nxyz =  `header -si $argv[1] | sed '/[[:cntrl:]]/s///g'`
@ nxt = $nxyz[1]
@ nyt = $nxyz[2]
@ nzt = $nxyz[3]

@ xbh = $nxt - $xhi
@ ybh = $nyt - $yhi
@ zbh = $nzt - $zhi

# If a transform file and original volume are entered, use that to modify the
# borders appropriately

if ($#argv == 4) then
    set inxyz =  `header -si $argv[4] | sed '/[[:cntrl:]]/s///g'`

    # Get the transformed coordinates of the corners of the original volume

    set sol = `cat $argv[3]`
    set a11 = $sol[1]
    set a12 = $sol[3]
    set odx = $sol[4]
    set a21 = $sol[9]
    set a22 = $sol[11]
    set ody = $sol[12]
    @ inx = $inxyz[1] / 2
    @ iny = $inxyz[3] / 2
    @ onx = $nxt / 2
    @ ony = $nzt / 2

    @ xll = `echo $a11 $a12 $odx $inx $iny $onx | awk '{print int(-$1 * $4 - $2 * $5 + $3 + $6)}'`
    @ xlr = `echo $a11 $a12 $odx $inx $iny $onx | awk '{print int($1 * $4 - $2 * $5 + $3 + $6)}'`
    @ xul = `echo $a11 $a12 $odx $inx $iny $onx | awk '{print int(-$1 * $4 + $2 * $5 + $3 + $6)}'`
    @ xur = `echo $a11 $a12 $odx $inx $iny $onx | awk '{print int($1 * $4 + $2 * $5 + $3 + $6)}'`
    @ yll = `echo $a21 $a22 $ody $inx $iny $ony | awk '{print int(-$1 * $4 - $2 * $5 + $3 + $6)}'`
    @ ylr = `echo $a21 $a22 $ody $inx $iny $ony | awk '{print int($1 * $4 - $2 * $5 + $3 + $6)}'`
    @ yul = `echo $a21 $a22 $ody $inx $iny $ony | awk '{print int(-$1 * $4 + $2 * $5 + $3 + $6)}'`
    @ yur = `echo $a21 $a22 $ody $inx $iny $ony | awk '{print int($1 * $4 + $2 * $5 + $3 + $6)}'`

    # Sort the values in order to find second and third ones

    @ x1 = $xll
    @ x2 = $xlr
    @ x3 = $xul
    @ x4 = $xur

    if ($x1 > $x2) then
	@ tmp = $x1
	@ x1 = $x2
	@ x2 = $tmp
    endif
    if ($x1 > $x3) then
	@ tmp = $x1
	@ x1 = $x3
	@ x3 = $tmp
    endif
    if ($x1 > $x4) then
	@ tmp = $x1
	@ x1 = $x4
	@ x4 = $tmp
    endif
    if ($x2 > $x3) then
	@ tmp = $x2
	@ x2 = $x3
	@ x3 = $tmp
    endif
    if ($x2 > $x4) then
	@ tmp = $x2
	@ x2 = $x4
	@ x4 = $tmp
    endif
    if ($x3 > $x4) then
	@ tmp = $x3
	@ x3 = $x4
	@ x4 = $tmp
    endif

    @ y1 = $yll
    @ y2 = $ylr
    @ y3 = $yul
    @ y4 = $yur

    if ($y1 > $y2) then
	@ tmp = $y1
	@ y1 = $y2
	@ y2 = $tmp
    endif
    if ($y1 > $y3) then
	@ tmp = $y1
	@ y1 = $y3
	@ y3 = $tmp
    endif
    if ($y1 > $y4) then
	@ tmp = $y1
	@ y1 = $y4
	@ y4 = $tmp
    endif
    if ($y2 > $y3) then
	@ tmp = $y2
	@ y2 = $y3
	@ y3 = $tmp
    endif
    if ($y2 > $y4) then
	@ tmp = $y2
	@ y2 = $y4
	@ y4 = $tmp
    endif
    if ($y3 > $y4) then
	@ tmp = $y3
	@ y3 = $y4
	@ y4 = $tmp
    endif

    # Use the middle values to modify the borders

    if ($x2 > 0) @ xbl += $x2
    if ($x3 < $nxt) @ xbh += ($nxt - $x3)
    if ($y2 > 0) @ zbl += $y2
    if ($y3 < $nzt) @ zbh += ($nzt - $y3)

endif
#
if ($nx == 1) then
  @ dx = 0
  @ xs = $nxt / 2
else
  @ dx = (($nxt - $nxl) - ($xbl + $xbh))/($nx - 1)
  @ xs = $xbl + ($nxl / 2)
endif
#
if ($ny == 1) then
  @ dy = 0
  @ ys = $nyt / 2
else
  @ dy = (($nyt - $nyl) - ($ybl + $ybh))/($ny - 1)
  @ ys = $ybl + ($nyl / 2)
endif
#
if ($nz == 1) then
  @ dz = 0
  @ zs = $nzt / 2
else
  @ dz = (($nzt - $nzl) - ($zbl + $zbh))/($nz - 1)
  @ zs = $zbl + ($nzl / 2)
endif
#
@ nxtap = ($nxl + 9) / 10
@ nytap = ($nyl + 9) / 10
@ nztap = ($nzl + 9) / 10
@ nxpad = ($nxl + 4) / 5
@ nypad = ($nyl + 4) / 5
@ nzpad = ($nzl + 4) / 5
#
set pad1 = "$tmpdir/pad1.$$"
set pad2 = "$tmpdir/pad2.$$"
set cortmp = "$tmpdir/corr3d.$$"
onintr clean
#
@ numtot = $nx * $ny * $nz
echo $numtot "    positions"
while ($nz > 0)
@ iy = $ny
@ yp = $ys
while ($iy > 0)
  @ ix = $nx
  @ xp = $xs
  while ($ix > 0)
tapervoledge >/dev/null <<HERESTRING
$argv[1]
$pad1
$nxl,$nyl,$nzl
$xp,$yp,$zs
$nxpad,$nypad,$nzpad
$nxtap,$nytap,$nztap
HERESTRING
tapervoledge >/dev/null <<HERESTRING
$argv[2]
$pad2
$nxl,$nyl,$nzl
$xp,$yp,$zs
$nxpad,$nypad,$nzpad
$nxtap,$nytap,$nztap
HERESTRING
#    set clipcom = "clip corr -3d -cx $xp -cy $yp -cz $zs -ix $nxl -iy $nyl -iz $nzl $argv corr3d.tmp"
    set displ = `clip corr -3d -n 0 "$pad1" "$pad2" "$cortmp" | grep '^(' | sed '/.*(\(.*\)).*/s//\1/'`
    echo $xp $yp $zs $displ
    @ ix--
    @ xp = $xp + $dx
  end
  @ iy--
  @ yp = $yp + $dy
end
@ nz--
@ zs = $zs + $dz
end

clean:
\rm -f "$cortmp"* "$pad1"* "$pad2"*
