/* 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,outscale=1.0,offset=40.0,outoffset=40.0;
    float float1, float2, float3;
    int header=512;
    int normalize=0;
    float sum,sumd,sumsq,sumdsq,avg,diff,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], "-outscale") && (argc >= (i+1)))
            {
                outscale = atof(argv[i+1]);
            }
            if(strstr(argv[i], "-normalize") && (argc >= (i+1)))
            {
                normalize = 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];
	avg = ( scale1*(float1-offset) + scale2*(float2-offset) ) / 2.0 ;
	diff = scale1*(float1-offset) - scale2*(float2-offset) ;
        if( normalize )
        {
            float3 = outscale * (diff / avg) + outoffset;
        }
        else
        {
            float3 = outscale * (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;
}