#! /bin/tcsh -f
#
#	Jiffy script for making a linear combination of two SMV images
#
#	written by James Holton   8-26-12
#

# zero-photon pixel value
set offset = 40
set outoffset = 40
# image header size
set headersize = 512
# swap bytes?
set swab = "N"

set outprefix = mix_

set image1 = ""
set scale1 = ""
set image2 = ""
set scale2 = ""

set tempfile = /tmp/${USER}/temp_mix$$
touch ${tempfile}
if(! -e ${tempfile}) set tempfile = tempfil
touch ${tempfile}images


goto compile_int_add
# jump to bottom to deploy the script
compiled_int_add:

# scan the command line
set i = 0
while ( $i < $#argv )
    @ i = ( $i + 1 )
    @ j = ( $i + 1 )
    set arg = "$argv[$i]"
    set num = `echo $arg | awk '$1+0>-1e-32 && ! /[a-z_A-Z_]/{print ($1+0)}'`
    if("$num" != "") then
	# recognized a number
	if("$scale1" == "") then
	    set scale1 = "$num"
	else
	    if("$scale2" != "") echo "WARNING: more than two scales given! "
	    set scale2 = "$num"
	endif
	continue
    endif
    if("$arg" == "-swap" || "$arg" == "-swab") then
	set swab = "Y"
	continue
    endif
    if("$arg" == "-offset" && $j <= $#argv) then
	set num = `echo $argv[$j] | awk '$1+0!=0 && ! /[a-z_A-Z_]/{print ($1+0)}'`
	if("$num" == "") then
	    echo "WARNING: $argv[$j] does not look like an offset"
	else
	    set offset = $num
	    set outoffset = $num
	    @ i = ( $i + 1 )
	endif
	continue
    endif
    if("$arg" == "-outoffset" && $j <= $#argv) then
	set num = `echo $argv[$j] | awk '$1+0!=0 && ! /[a-z_A-Z_]/{print ($1+0)}'`
	if("$num" == "") then
	    echo "WARNING: $argv[$j] does not look like an offset"
	else
	    set outoffset = $num
	    @ i = ( $i + 1 )
	endif
	continue
    endif
    if("$arg" == "-outprefix" && $j <= $#argv) then
	set outprefix = "$argv[$j]"
	@ i = ( $i + 1 )
	continue
    endif
    if("$arg" =~ *.img) then
	# recognized potential image filename
	ls -1 `echo $arg | awk '{gsub("#","[0-9]");print}'` >> ${tempfile}images
	continue
    endif
end

# defaults for uninitialized variables
if("$scale1" == "" && "$scale2" == "") then
    set scale1 = 0.5
    set scale2 = 0.5
endif
if("$scale2" == "") then
    # must be a value for scale1
    set scale2 = `echo $scale1 | awk '{print 1-$1}'`
endif
if("$scale2" == "") then
    # must be a value for scale2 ? 
    set scale1 = `echo $scale2 | awk '{print 1-$1}'`
endif

echo "mixing ratio: $scale1 $scale2"
echo "zero photon offset: $offset ($outoffset in output)"
echo "swap bytes: $swab"
echo "output prefix: $outprefix"
# fail if not enough info?

set images = `cat ${tempfile}images | wc -l`
if($images < 2) then
    set BAD = "must specify at least two images."
    goto exit
endif

# figure out image prefixes
cat ${tempfile}images |\
awk '{gsub(".img$","");base=$1;\
      while(gsub("[0-9#]$",""));\
      num=substr(base,1+length($0));\
 print $0,num}' |\
cat >! ${tempfile}parsed
awk '{print $1}' ${tempfile}parsed |\
  awk '! seen[$1]{print; ++seen[$1]}' >! ${tempfile}prefixes
awk '{print $2}' ${tempfile}parsed |\
  awk '! seen[$1]{print; ++seen[$1]}' >! ${tempfile}nums
set prefixes = `cat ${tempfile}prefixes | wc -l`
set prefix1 = `head -n 1 ${tempfile}prefixes`
set prefix2 = `tail -n 1 ${tempfile}prefixes`
set nums = `cat ${tempfile}nums`
set nums = `cat ${tempfile}nums`
if($images == 2) then
    set image1 = `head -n 1 ${tempfile}images`
    set image2 = `tail -n 1 ${tempfile}images`
    set prefixes = 2
    set nums = twofiles
endif
if($prefixes < 2) then
    set BAD = "not enough image patterns: $prefix1 $prefix2 $nums"
    goto exit
endif
if($prefixes > 2) then
    set BAD = "too many image patterns, unable to identify pairs."
    goto exit
endif
rm -f ${tempfile}parsed ${tempfile}nums ${tempfile}prefixes
rm -f ${tempfile}images

foreach num ( $nums )

if("$num" != "twofiles") then
    set image1 = ${prefix1}${num}.img
    set image2 = ${prefix2}${num}.img
endif

set outfile = ${outprefix}${num}.img

if(! -e "$image1") then
    set BAD = "unable to find $image1"
    goto exit
endif
if(! -e "$image2") then
    set BAD = "unable to find $image2 for $image1"
    goto exit
endif
if(-s "$outfile") then
    set BAD = "$outfile already exists.  I will not delete it."
    goto exit
endif

echo "mixing:"
echo "  $scale1 x $image1 "
echo "+ $scale2 x $image2 "

set test = `dd bs=512 count=1 if=$image1 |& awk -F "=" '/HEADER/{print $2+0;exit}'`
if("$test" != "") then
    set headersize = $test
endif

int_add -header $headersize -scale1 $scale1 -scale2 $scale2 $image1 $image2 $outfile -outoffset $outoffset


set size1   = `ls -l $image1 | awk '{print $5}'`
set sizeout = `ls -l $outfile | awk '{print $5}'`
if("$size1" != "$sizeout") then
    set BAD = "size mismatch.  something is worng!!! "
    goto exit
endif
echo "= $outfile"

end



exit:

if($?BAD) then
    echo "ERROR: $BAD"
    set me = `basename $0`
    if("$me" == "") set me = img_mix.com
    cat << EOF
usage: 
$me spots_###.img 0.5 background_###.img 0.5
EOF
    exit 9
endif

echo "all done."
exit


compile_int_add:

set path = ( . `dirname $0` $path )

rehash
set test = `int_add |& awk '/usage: int_add/{print 1}'`
if("$test" == "1") goto compiled_int_add

cat << EOF >! int_add.c

/* add two binary "unsigned short int" files together                                           -James Holton  2-14-12

example:

gcc -O -O -o int_add int_add.c -lm
./int_add file1.img file2.img output.img

 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>


char *infile1name = "";
char *infile2name = "";
FILE *infile1 = NULL;
FILE *infile2 = NULL;
char *outfilename = "output.bin\\0";
FILE *outfile = NULL;

int main(int argc, char** argv)
{

    int n,i,j,k,pixels;
    unsigned short int *outimage;
    unsigned short int *inimage1;
    unsigned short int *inimage2;
    float scale1=1.0,scale2=1.0,offset=40.0,outoffset=40.0;
    float float1, float2, float3;
    int header=512;
    float sum,sumd,sumsq,sumdsq,avg,rms,rmsd,max,min;


    /* check argument list */
    for(i=1; i<argc; ++i)
    {
        if(strlen(argv[i]) > 4)
        {
            if(strstr(argv[i]+strlen(argv[i])-4,".bin") || strstr(argv[i]+strlen(argv[i])-4,".img"))
            {
                printf("filename: %s\\n",argv[i]);
                if(infile1 == NULL){
                    infile1 = fopen(argv[i],"r");
                    if(infile1 != NULL) infile1name = argv[i];
                }
                else
                {
                    if(infile2 == NULL){
                        infile2 = fopen(argv[i],"r");
                        if(infile2 != NULL) infile2name = argv[i];
                    }
                    else
                    {
                        outfilename = argv[i];
                    }
                }
            }
        }

        if(argv[i][0] == '-')
        {
            /* option specified */
            if(strstr(argv[i], "-scale1") && (argc >= (i+1)))
            {
                scale1 = atof(argv[i+1]);
            }
            if(strstr(argv[i], "-scale2") && (argc >= (i+1)))
            {
                scale2 = atof(argv[i+1]);
            }
            if(strstr(argv[i], "-offset") && (argc >= (i+1)))
            {
                offset = atof(argv[i+1]);
                outoffset = atof(argv[i+1]);
            }
            if(strstr(argv[i], "-outoffset") && (argc >= (i+1)))
            {
                outoffset = atof(argv[i+1]);
            }
            if(strstr(argv[i], "-header") && (argc >= (i+1)))
            {
                header = atof(argv[i+1]);
            }
        }
    }

    if(infile1 == NULL || infile2 == NULL){
        printf("usage: int_add file1.bin file2.bin [outfile.bin] -scale1 1.0 -scale2 1.0 -offset 40 -header 512\\n");
        printf("options:\\n");\\
        printf("\\t-scale1 \\tscale factor for first file\\n");
        printf("\\t-scale2 \\tscale factor for second file\\n");
        printf("\\t-offset \\tinteger offset to be subtracted from each file\\n");
        printf("\\t-outoffset\\tinteger offset to be added to the output file\\n");
        printf("\\t-header \\tnumber of bytes to ignore in header of each file\\n");
        exit(9);
    }


    /* load first unsigned short int-image */
    fseek(infile1,0,SEEK_END);
    n = ftell(infile1);
    fseek(infile1,0,SEEK_SET);
    fseek(infile2,0,SEEK_SET);
    inimage1 = calloc(n,1);
    inimage2 = calloc(n,1);
    fread(inimage1,n,1,infile1);
    fclose(infile1);
    fread(inimage2,n,1,infile2);
    fclose(infile2);

    pixels = (n-header)/sizeof(unsigned short int);
    outfile = fopen(outfilename,"w");
    if(outfile == NULL)
    {
        printf("ERROR: unable to open %s\\n", outfilename);
        exit(9);
    }

    printf("header = %d bytes pixel zero offset = %g (%g in output)\\n",header,offset,outoffset);

    outimage = calloc(n,1);
    for(i=0;i<=header/sizeof(unsigned short int);++i)
    {
        outimage[i] = inimage1[i];
    }
    sum = sumsq = sumd = sumdsq = 0.0;
    min = 65535.0;
    max = 0.0;
    for(i=header/sizeof(unsigned short int);i<pixels;++i)
    {
        float1 = (float) inimage1[i];
        float2 = (float) inimage2[i];
        float3 = scale1*(float1-offset) + scale2*(float2-offset) + outoffset;
        outimage[i] = (unsigned short int) float3;
        if(outimage[i]>max) max=outimage[i];
        if(outimage[i]<min) min=outimage[i];
        sum += float3;
        sumsq += float3*float3;
    }
    avg = sum/pixels;
    rms = sqrt(sumsq/pixels);
    for(i=0;i<pixels;++i)
    {
        sumd   += outimage[i] - avg;
        sumdsq += (outimage[i] - avg) * (outimage[i] - avg);
    }
    rmsd = sqrt(sumdsq/pixels);
    printf("max = %g min = %g\\n",max,min);
    printf("mean = %g rms = %g rmsd = %g\\n",avg,rms,rmsd);

    printf("writing %s as a %d-byte header with %d %ld-byte unsigned short ints\\n",outfilename,header,pixels,sizeof(unsigned short int));
    outfile = fopen(outfilename,"w");
    fwrite(outimage,n,1,outfile);
    fclose(outfile);


    return;
}

EOF
gcc -o int_add int_add.c -lm

rehash
set test = `int_add | awk '/usage: int_add/{print 1}'`
if("$test" == "1") goto compiled_int_add

set uname = `uname`
curl -o int_add http://bl831.als.lbl.gov/~jamesh/bin_stuff/${uname}/int_add
chmod a+x int_add

rehash
set test = `int_add | awk '/usage: int_add/{print 1}'`
if("$test" == "1") goto compiled_int_add


set BAD = "unable to find/compile or otherwise obtain integer addition program: int_add"
goto exit

goto compiled_int_add