#! /bin/csh -f
#
#	locate the backstop shadow on an ADSC/SMV image
#
#
#
set tempfile = tempfile

set image = "$1"
set scale = 1
set scale = auto

# everything is in mm
set backstopwidth = 0.3
set backstopdist  = 100.0
set backstopdist  = 60.0



if(! -e "$image" && $?XFORMSTATUSFILE) then
    set image = `awk '{print $2}' $XFORMSTATUSFILE`
endif
if(! -e "$image") then
    echo "usage: find_backstop.com test_0_001.img"
    exit 9
endif

# extract the SMV header
awk 'BEGIN{RS="\f"} {print; exit}' $image >! $tempfile

set wave        = `awk -F "=" '/^WAVE/{print $2+0}' $tempfile | tail -1`
set dist        = `awk -F "=" '/^DIST/{print $2+0}' $tempfile | tail -1`
set beam_center = `awk -F "=" '/^BEAM_CENTER_/{print $2+0}' $tempfile | tail -2`
set pixel       = `awk -F "=" '/^PIXEL_SIZE/{print $2+0}' $tempfile | tail -1`
set width       = `awk -F "=" '/^PIXEL_SIZE/{pixelsize = $2+0} /^SIZE1/{pixels=$2+0} END{print pixelsize*pixels}' $tempfile`

set pix_center  = `echo $beam_center $pixel $width | awk '{print ($4-$2)/$3, $1/$3}'`

# make sure we cut out an area bigger than the backstop
set radius = `echo $backstopwidth $backstopdist $pixel $dist | awk '{r=($1/2)*$4/$2/$3; printf "%d", r}'`
set pix_size = `echo $backstopwidth $backstopdist $pixel $dist | awk '{r=($1/2)*$4/$2/$3; printf "%d", 3.14*r*r}'`
set boxedge = `echo $backstopwidth $backstopdist $pixel $dist | awk '{printf "%d", 4*$1*$4/$2/$3}'`
set size = "${boxedge}x${boxedge}"
set box = `echo $pix_center $boxedge | awk '{printf "%d %d %d %d\n", $1-$NF/2, $2-$NF/2, $1+$NF/2, $2+$NF/2}'`

# create an image of what we EXPECT the beamstop to look like
convert -size $size -depth 8 -fill white -draw "color 0,0 reset" nutn.NULL needle.jpg
set circle = `echo $radius $boxedge | awk '{r=$1; mid=$2/2; printf "%d,%d %d,%d", mid,mid,mid-r,mid}'`
mogrify -pen black -fill black -draw "circle $circle" -blur 10x3 needle.jpg

cutout:
# do the actual pixel cutout
adsc2pgm $image pixels.pgm -scale $scale -box $box
#adsc2pgm $image pixels.pgm -box $box
#adsc2pgm $image pixels.pgm -scale 1 -box 1024 1024 2048 2048 
mogrify -flop -rotate 90 pixels.pgm
convert -filter Point -geometry 400% -negate pixels.pgm pixels.png

# blot out the direct beam
set hdiv = 2.0
set ellipse = `echo $dist $pixel $hdiv $boxedge | awk '{hr=$1*$3/1500/$2+1; vr=$1*0.3/1500/$2+1; mid=$NF/2; printf "%d,%d %d,%d", mid,mid,hr,vr}'`
convert -pen black -fill "rgb(10000,10000,10000)" -draw "ellipse $ellipse 0,360" pixels.pgm nospot.pgm

#make a jpeg of it
convert nospot.pgm haystack.jpg

# use the CC to find best fit
matchup.com needle.jpg haystack.jpg edges | tee matchup.log

set delta  = `awk '/with cc/{print $6,$7;exit}' matchup.log`
rm -f matchup.log

set center = `echo $delta $boxedge | awk '{print $1+$3/2,$2+$3/2}'`

goto report

exit

# pre-defined threshold levels?
convert -segment $pix_size -normalize -negate nospot.pgm segmented.pgm
#convert -threshold 63000 segmented.pgm pane_mask.pgm
#convert -threshold 20 segmented.pgm shadow_mask.pgm
convert -threshold 20 -negate nospot.pgm pane_mask.pgm
convert -threshold 14000 -negate nospot.pgm shadow_mask.pgm
composite -compose Subtract shadow_mask.pgm pane_mask.pgm backstop_mask.pgm
#display pane_mask.pgm
#display shadow_mask.pgm
#display backstop_mask.pgm

convert backstop_mask.pgm backstop_mask.GRAY
od -t u1 -v backstop_mask.GRAY |\
awk -v edge=$boxedge '\
    {for(i=2;i<=NF;++i) {++n; x=n%edge;y=int(n/edge);\
	if($i!=0) print x,y}}' |\
cat >! xy_mask.txt
rm -f backstop_mask.GRAY

set center = `awk '{++n;sumx+=$1;sumy+=$2} END{if(! n) exit; print sumx/n, sumy/n, sqrt(n/3.14)}' xy_mask.txt`
if($#center < 3) then
    echo "cannot see backstop..."
    if ( "$scale" != "auto") then
	set scale = auto
	goto cutout
    endif
endif
set radius = $center[3]
set center = ( $center[1] $center[2] )

goto report

convert -threshold 20 -edge 0 nospot.pgm pane_edge.pgm

set threshold = 13000
foreach threshold ( 20 100 1000 10000 20000 )
convert -threshold $threshold -edge 0 nospot.pgm shadow_edge.pgm
composite -compose Subtract shadow_edge.pgm pane_edge.pgm backstop_edge.pgm
#display backstop_edge.pgm

# now extract all this as a list of numbers
convert backstop_edge.pgm backstop_edge.GRAY
od -t u1 -v backstop_edge.GRAY |\
awk -v edge=$boxedge '\
    {for(i=2;i<=NF;++i) {++n; x=n%edge;y=int(n/edge);\
	if($i!=0) print x,y}}' |\
cat >! xy_edge.txt
rm -f backstop_edge.GRAY

set center = `awk '{++n;sumx+=$1;sumy+=$2} END{if(! n) exit; print sumx/n, sumy/n}' xy_mask.txt`
set radius = `awk -v xc=$center[1] -v yc=$center[2] '{++n;sumr+=sqrt(($1-xc)^2+($2-yc)^2)} END{if(! n) exit; print sumr/n}' xy_edge.txt`

echo $threshold $radius
end

# fit the edge detection to an elipse
cat << EOF | gnuplot -persist
plot '${tempfile}xy' using 1:2

xc1=$center[1]
yc1=$center[2]
r1=$radius
xc1=$center[1]
yc1=$center[2]
r1=33
xc1=66
yc1=56
f(x,y) = ( r1**2 -(x-xc1)**2 -(y-yc1)**2 )**2
fit 1-f(x,y) 'xy_edge.txt' using 1:2:(1) via xc1,yc1
fit 1-f(x,y) 'xy_edge.txt' using 1:2:(1) via xc1,yc1,r1

splot 'xy_edge.txt' using 1:2:(1), f(x,y)
exit

xc2=xc1-2
yc2=yc1-2
f(x,y) = ( 2*(r1**2) -(x-xc1)**2 -(y-yc1)**2 -(x-xc2)**2 -(y-yc2)**2 )**2

fit 1-f(x,y) 'xy_edge.txt' using 1:2:(1) via xc1,yc1
fit 1-f(x,y) 'xy_edge.txt' using 1:2:(1) via xc1,yc1,xc2,yc2
fit 1-f(x,y) 'xy_edge.txt' using 1:2:(1) via xc1,yc1,xc2,yc2,r1
splot 'xy_edge.txt' using 1:2:(1), f(x,y)

r2=r1
f(x,y) = ( r1**2 -(x-xc1)**2 -(y-yc1)**2 )**2 + ( r2**2 -(x-xc2)**2 -(y-yc2)**2 )**2
fit 1-f(x,y) 'xy_edge.txt' using 1:2:(1) via xc1,yc1,r1,xc2,yc2,r2
splot 'xy_edge.txt' using 1:2:(1), f(x,y)

EOF

# get the center and radius from the gnuplot fit.log
cat fit.log |\
awk '/^xc1 / && $2=="="{xc1=$3}\
     /^yc1 / && $2=="="{yc1=$3}\
     /^r1 /  && $2=="="{r1=$3}\
  END{print xc1,yc1,r1}' |\
cat >! backstop.txt
set center = `awk '{print $1,$2}' backstop.txt`
set radius = `awk '{print $3}' backstop.txt`


report:
# calculate vector between the beam and the backstop center
set delta = `echo $center $boxedge | awk '{print $1-$3/2,$2-$3/2}'`

# convert this to real position of the backstop
set offset = `echo $delta $pixel $backstopdist $dist | awk '{printf "%.3f %.3f", -$1*$3*$4/$5,$2*$3*$4/$5}'`
set correction = `echo $offset | awk '{print $1,$2}'`

echo ""
echo "beam     is at $beam_center mm"
echo "backstop is $offset mm off-center"
echo ""
echo "beamstop.com $correction"

# decorate the image with our fitted beamstop
set drawcenter = `echo $center | awk '{printf "%d,%d", $1, $2}'`
set edge1  = `echo $center $radius | awk '{printf "%d,%d", $1, $2-$3}'`
set edge2  = `echo $center $radius | awk '{printf "%d,%d", $1, $2-$3+1}'`
convert -pen black -fill black -draw "rectangle 0,0 $boxedge,$boxedge" \
      -fill white -draw "circle $drawcenter $edge1" -fill black \
      -fill black -draw "circle $drawcenter $edge2" -fill white \
      pixels.pgm circle.pgm
composite -compose Plus pixels.pgm circle.pgm marked.pgm

# decorate the image with the direct beam position
set hline  = `echo $boxedge | awk '{$1/=2; printf "%d,%d %d,%d", $1-5,$1,$1+5,$1}'`
set vline  = `echo $boxedge | awk '{$1/=2; printf "%d,%d %d,%d", $1,$1-5,$1,$1+5}'`
mogrify -pen white -fill white\
   -draw "line $hline[1] $hline[2]" \
   -draw "line $vline[1] $vline[2]" \
   -negate \
   marked.pgm

convert -filter Point -geometry 400% marked.pgm marked.png

set sane = `echo $offset $backstopwidth | awk '{print (sqrt($1^2+$2^2) < $3)}'`
if( ! $sane) then
    echo "ERROR: this is too far off to be for real! "
    exit 9
endif

#display marked.png