/*************************************************
           concatenation of the adjacent SRTM DEMs  
	         Yuri Fialko 3/24/01
      compile: gcc -O -o SRTMpaste SRTMpaste.c
*****************************************************/
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define null 0
#define nmax 20000
#define maxchar 80
#define pad 1
#define init 1
#define background 30000

FILE *fpin1;
FILE *fpin2;
FILE *fpin1r;
FILE *fpin2r;
FILE *fpout;
FILE *fpoutr;

short intpl=0;
int       ddx,ddy,dxf,dyf;
short int       idem1[nmax][nmax],idem2[nmax][nmax],d;
long            row,col,width,length,widmax,lenmax;
char *lines[maxchar];
char *words[maxchar];
char *units[maxchar];
int xym[5];
int k,i,j;
double xyfirst[5],dx[3],dy[3];
/*char *fdname;*/
char line[maxchar]="";
char name1[maxchar]="";
char name2[maxchar]="";
char name1r[maxchar]="";
char name2r[maxchar]="";
char out_file[maxchar]="";
char out_filer[maxchar]="";
char flag[maxchar]="";

int main(int,char*[]);
int split(char []);
short aver(long,long,long,long);
void readin(char[], char[]);

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

/*     Print error if arguments were not included on command line        */

printf ("\n               DEM Paste Utility");
printf ("\n    Yuri Fialko(yfialko@ucsd.edu), Mar 2001.\n");

/*      Open line input file         */
if (argc < 4)
  {
  printf ("\nUsage: dempaste dem_in1 dem_in2 dem_out [dx] [dy] [-i]\n ");
  printf ("    dem_in1,2:  input DEM files  without .dem extension \n");
  printf ("    dem_out:  output file without .dem extension\n");
  printf ("    dx,dy:  forced horizontal and vertical offset between dem_in1 \n");
  printf ("            and dem_in2, in pixels (default is 0 0)\n");
  printf ("    -i:  interpolation; recommended only if there \n");
  printf ("        are gaps in the SDTS data \n");
  }

if (argc < 2) /*      Prompt for first input file name      */
  {
  printf ("\n\nEnter first input DEM file name:\n ");
  scanf ("%s",name1);
  }
else 
  strcpy (name1, argv[1]);
  strcat(name1,".dem") ;
  strcpy (name1r, name1);
  strcat(name1r,".rsc") ;

/*    
  Prompt for output report file name      */
  
if (argc < 3)  /*      prompt for second input file name      */
  {
    printf ("\nEnter second input DEM file name:\n ");
    scanf ("%s",name2); 
  }
  else
    strcpy (name2, argv[2]);
  strcat(name2,".dem") ;
  strcpy (name2r, name2);
  strcat(name2r,".rsc") ;
  

  if (argc < 4)  /*      prompt for  output dem file name      */
  {
    printf 
("\nEnter output file name:\n ");
    scanf ("%s",out_file); 
  }
  else
    strcpy (out_file, argv[3]);
  strcat(out_file,".dem") ;
  strcpy(out_filer,out_file) ;
  strcat(out_filer,".rsc") ;

  if (argc > 4)  /*      prompt for 2 digit layer number      */
  {
    strcpy (flag, argv[4]);
    dxf=atoi(flag);
  }
  else
   dxf=0; 

  if (argc > 5)  /*      prompt for 2 digit layer number      */
  {
    strcpy (flag, argv[5]);
    dyf=atoi(flag);
  }
  else
   dyf=0; 

  if (argc > 6)  /*      prompt for 2 digit layer number      */
  {
    strcpy (flag, argv[6]);
    if (strcmp(flag,"-i") == 0)
      intpl=1;
    else
    {
      printf ("Correct flag for interpolation is <-i> (no argument for no interpolation)\n");
      exit (1);
    }
  } 
  intpl=0;  /* no interpolation !!! */

/*       Open output .dem file           */
  fpin1 = fopen (name1,"rb");
  fpin2 = fopen (name2,"rb");
  fpin1r = fopen (name1r,"r");
  fpin2r = fopen (name2r,"r");
  if (fpin1 == null)
         {
         printf ("\nERROR OPENING .DEM FILE %s\n",name1);
         printf ("Are you sure you're in the right directory?\n");
         exit (1);
         }

  if (fpin2 == null)
         {
         printf ("\nERROR OPENING .DEM FILE %s\n",name2);
         printf ("Are you sure you're in the right directory?\n");
         exit (1);
         }

  if (fpin1r == null)
         {
         printf ("\nERROR OPENING .rsc FILE %s\n",name1r);
         exit (1);
         }

  if (fpin2r == null)
         {
         printf ("\nERROR OPENING .rsc FILE %s\n",name2r);
         exit (1);
         }


/*  get info from .rsc files */

  i=0;
    while(fgets(line,sizeof(line),fpin1r)!=NULL){
    if ((lines[i] = (char *)malloc(strlen(line)+1)) == NULL)
    {
      printf ("\n Failed to malloc %s\n",line);
      exit (1);
    }
    strcpy(lines[i],line);
    k=split(lines[i]);
    if (strcmp(words[0],"WIDTH") == 0){
      xym[1]=atoi(words[1]);}
    if (strcmp(words[0],"FILE_LENGTH") == 0){
      xym[2]=atoi(words[1]);}
    if (strcmp(words[0],"X_FIRST") == 0){
      xyfirst[1]=atof(words[1]);}
    if (strcmp(words[0],"Y_FIRST") == 0){
      xyfirst[2]=atof(words[1]);}
    if (strcmp(words[0],"X_STEP") == 0){
      dx[1]=(double)atof(words[1]);}
    if (strcmp(words[0],"Y_STEP") == 0){
      dy[1]=atof(words[1]);}
    if (strcmp(words[0],"Y_UNIT") == 0){
     if ((units[1] = (char *)malloc(strlen(words[1])+1)) == NULL)
     {
      printf ("\n Failed to malloc %s\n",words[1]);
      exit (1);
     }
     strcpy(units[1],words[1]);}
    i++;
    }
  i=0;
    while(fgets(line,sizeof(line),fpin2r)!=NULL){
    if ((lines[i] = (char *)malloc(strlen(line)+1)) == NULL)
    {
      printf ("\n Failed to malloc %s\n",line);
      exit (1);
    }
    strcpy(lines[i],line);
    k=split(lines[i]);
    if (strcmp(words[0],"WIDTH") == 0){
      xym[3]=atoi(words[1]);}
    if (strcmp(words[0],"FILE_LENGTH") == 0){
      xym[4]=atoi(words[1]);}
    if (strcmp(words[0],"X_FIRST") == 0){
      xyfirst[3]=atof(words[1]);}
    if (strcmp(words[0],"Y_FIRST") == 0){
      xyfirst[4]=atof(words[1]);}
    if (strcmp(words[0],"X_STEP") == 0){
      dx[2]=atof(words[1]);}
    if (strcmp(words[0],"Y_STEP") == 0){
      dy[2]=atof(words[1]);}
    if (strcmp(words[0],"Y_UNIT") == 0){
     if ((units[2] = (char *)malloc(strlen(words[1])+1)) == NULL)
     {
      printf ("\n Failed to malloc %s\n",words[1]);
      exit (1);
     }
     strcpy(units[2],words[1]);}
    i++;
  }


    readin(units[1],units[2]);
    fclose (fpin1);
    fclose (fpin1r);
    fclose (fpin2);
    fclose (fpin2r);

  fpout = fopen (out_file,"wb");
  fpoutr = fopen (out_filer,"w");
  if (fpout == null)
         {
         printf ("\nERROR OPENING OUTPUT .DEM FILE %s\n",out_file);
         printf ("Is disk full?\n");
         exit (1);
         }

  if (fpoutr == null)
         {
         printf ("\nERROR OPENING .rsc FILE %s\n",out_filer);
         exit (1);
         }


  printf ("\n x,y_start pairs: %f %f %f %f \n",xyfirst[1],xyfirst[2],xyfirst[3],xyfirst[4]);
     
  if ((xym[1]+2*(pad-1) >= nmax) || (xym[2]+2*(pad-1) >= nmax))
         {
         printf ("\n DEM FILE %s is larger than nmax\n",name1);
         printf ("Requested size is %ld x %ld\n",xym[1]+2*(pad-1),xym[2]+2*(pad-1));
         printf ("\n Contact me and I'll try to increase memory size\n");
         exit (1);
         }

  if ((xym[3]+2*(pad-1) >= nmax) || (xym[4]+2*(pad-1) >= nmax))
         {
         printf ("\n DEM FILE %s is larger than nmax\n",name2);
         printf ("Requested size is %ld x %ld\n",xym[3]+2*(pad-1),xym[4]+2*(pad-1));
         printf ("\n Contact me and I'll try to increase memory size\n");
         exit (1);
         }

    ddx=(int)(xyfirst[1]-xyfirst[3])/dx[1];
    ddy=(int)(xyfirst[2]-xyfirst[4])/dy[1];
  /* initialize arrays */
/*     
  for (i=0; i<=nmax-1; i++)
  for (j=0; j<=nmax-1; j++)
  {
    idem1[i][j]=init;
    idem2[i][j]=init;
  }
*/
  if((xyfirst[2]==xyfirst[4])&&(xym[2]==xym[4])&&(xyfirst[1]<xyfirst[3]))
      /* 1st file is to the west */
    {
        width=xym[1]+xym[3]-1;
        length=xym[2];
        printf ("doing rsc file...\n");
        fprintf(fpoutr,"WIDTH                 %d \n",width);
        fprintf(fpoutr,"FILE_LENGTH           %d \n",length);
        printf("WIDTH         %d\n",width);
        printf("FILE_LENGTH   %d\n",length);
        fprintf(fpoutr,"XMIN                  0 \n");
        fprintf(fpoutr,"XMAX                  %d \n",width-1);
        fprintf(fpoutr,"YMIN                  0 \n");
        fprintf(fpoutr,"YMAX                  %d \n",length-1);
        fprintf(fpoutr,"X_FIRST               %.9f \n",xyfirst[1]);
        fprintf(fpoutr,"Y_FIRST               %.9f \n",xyfirst[2]);
        
        printf ("doing binary dump...\n");
	/*        printf("ddy %d dy %d dyf %d\n",ddy,dy,dyf);
		  printf("dyu %d dyd %d \n",dyu,dyd);*/
        for (row=1; row <= length; row++)
	{
          for (col=1; col <= xym[1]-1; col++)
	  {
            if ((idem1[row][col] == background) && (intpl == 1))
	    {  
	      /*   d=aver(row,col,row-dy,col-xym[1]+dx);
		   fwrite(&d,sizeof(short),1,fpout);*/
            }
            else
              fwrite(&idem1[row][col],sizeof(short),1,fpout);
          }
          for (col=1; col <= xym[3]; col++)
	  {
            if ((idem2[row][col] == background) && (intpl == 1))
	    {  
	      /*   d=aver(row,col,row-dy,col-xym[1]+dx);
		   fwrite(&d,sizeof(short),1,fpout);*/
            }
            else
              fwrite(&idem2[row][col],sizeof(short),1,fpout);
          }
        }
    }
  if((xyfirst[2]==xyfirst[4])&&(xym[2]==xym[4])&&(xyfirst[1]>xyfirst[3]))
      /* 1st file is to the east */
    {
        width=xym[1]+xym[3]-1;
        length=xym[2];
        printf ("doing rsc file...\n");
        fprintf(fpoutr,"WIDTH                 %d \n",width);
        fprintf(fpoutr,"FILE_LENGTH           %d \n",length);
        printf("WIDTH         %d\n",width);
        printf("FILE_LENGTH   %d\n",length);
        fprintf(fpoutr,"XMIN                  0 \n");
        fprintf(fpoutr,"XMAX                  %d \n",width-1);
        fprintf(fpoutr,"YMIN                  0 \n");
        fprintf(fpoutr,"YMAX                  %d \n",length-1);
        fprintf(fpoutr,"X_FIRST               %.9f \n",xyfirst[3]);
        fprintf(fpoutr,"Y_FIRST               %.9f \n",xyfirst[4]);

        printf ("doing binary dump...\n");
	/*        printf("ddy %d dy %d dyf %d\n",ddy,dy,dyf);
		  printf("dyu %d dyd %d \n",dyu,dyd);*/
        for (row=1; row <= length; row++)
	{
          for (col=1; col <= xym[3]-1; col++)
	  {
            if ((idem2[row][col] == background) && (intpl == 1))
	    {  
	      /*   d=aver(row,col,row-dy,col-xym[1]+dx);
		   fwrite(&d,sizeof(short),1,fpout);*/
            }
            else
              fwrite(&idem2[row][col],sizeof(short),1,fpout);
          }
          for (col=1; col <= xym[1]; col++)
	  {
            if ((idem1[row][col] == background) && (intpl == 1))
	    {  
	      /*   d=aver(row,col,row-dy,col-xym[1]+dx);
		   fwrite(&d,sizeof(short),1,fpout);*/
            }
            else
              fwrite(&idem1[row][col],sizeof(short),1,fpout);
          }
        }
    }
  if((xyfirst[1]==xyfirst[3])&&(xym[1]==xym[3])&&(xyfirst[2]>xyfirst[4]))
      /* 1st file is to the north */
    {
        width=xym[1];
        length=xym[2]+xym[4]-1;
        printf ("doing rsc file...\n");
        fprintf(fpoutr,"WIDTH                 %d \n",width);
        fprintf(fpoutr,"FILE_LENGTH           %d \n",length);
        printf("WIDTH         %d\n",width);
        printf("FILE_LENGTH   %d\n",length);
        fprintf(fpoutr,"XMIN                  0 \n");
        fprintf(fpoutr,"XMAX                  %d \n",width-1);
        fprintf(fpoutr,"YMIN                  0 \n");
        fprintf(fpoutr,"YMAX                  %d \n",length-1);
        fprintf(fpoutr,"X_FIRST               %.9f \n",xyfirst[1]);
        fprintf(fpoutr,"Y_FIRST               %.9f \n",xyfirst[2]);

        printf ("doing binary dump...\n");
	/*        printf("ddy %d dy %d dyf %d\n",ddy,dy,dyf);
		  printf("dyu %d dyd %d \n",dyu,dyd);*/
        for (row=1; row <= xym[2]-1; row++)
	{
          for (col=1; col <= width; col++)
	  {
            if ((idem1[row][col] == background) && (intpl == 1))
	    {  
	      /*   d=aver(row,col,row-dy,col-xym[1]+dx);
		   fwrite(&d,sizeof(short),1,fpout);*/
            }
            else
              fwrite(&idem1[row][col],sizeof(short),1,fpout);
          }
        }
        for (row=1; row <= xym[4]; row++)
	{
          for (col=1; col <= width; col++)
	  {
            if ((idem2[row][col] == background) && (intpl == 1))
	    {  
	      /*   d=aver(row,col,row-dy,col-xym[1]+dx);
		   fwrite(&d,sizeof(short),1,fpout);*/
            }
            else
              fwrite(&idem2[row][col],sizeof(short),1,fpout);
          }
        }
    }
  if((xyfirst[1]==xyfirst[3])&&(xym[1]==xym[3])&&(xyfirst[2]<xyfirst[4]))
      /* 1st file is to the north */
    {
        width=xym[1];
        length=xym[2]+xym[4]-1;
        printf ("doing rsc file...\n");
        fprintf(fpoutr,"WIDTH                 %d \n",width);
        fprintf(fpoutr,"FILE_LENGTH           %d \n",length);
        printf("WIDTH         %d\n",width);
        printf("FILE_LENGTH   %d\n",length);
        fprintf(fpoutr,"XMIN                  0 \n");
        fprintf(fpoutr,"XMAX                  %d \n",width-1);
        fprintf(fpoutr,"YMIN                  0 \n");
        fprintf(fpoutr,"YMAX                  %d \n",length-1);
        fprintf(fpoutr,"X_FIRST               %.9f \n",xyfirst[3]);
        fprintf(fpoutr,"Y_FIRST               %.9f \n",xyfirst[4]);

        printf ("doing binary dump...\n");
	/*        printf("ddy %d dy %d dyf %d\n",ddy,dy,dyf);
		  printf("dyu %d dyd %d \n",dyu,dyd);*/
        for (row=1; row <= xym[4]-1; row++)
	{
          for (col=1; col <= width; col++)
	  {
            if ((idem2[row][col] == background) && (intpl == 1))
	    {  
	      /*   d=aver(row,col,row-dy,col-xym[1]+dx);
		   fwrite(&d,sizeof(short),1,fpout);*/
            }
            else
              fwrite(&idem2[row][col],sizeof(short),1,fpout);
          }
        }
        for (row=1; row <= xym[2]; row++)
	{
          for (col=1; col <= width; col++)
	  {
            if ((idem1[row][col] == background) && (intpl == 1))
	    {  
	      /*   d=aver(row,col,row-dy,col-xym[1]+dx);
		   fwrite(&d,sizeof(short),1,fpout);*/
            }
            else
              fwrite(&idem1[row][col],sizeof(short),1,fpout);
          }
        }
    }
/* copy the rest of the 1.rsc contents into out.rsc file */
  fprintf(fpoutr,"X_STEP               %.9f \n",dx[1]);
  fprintf(fpoutr,"Y_STEP               %.9f \n",dy[1]);
  fprintf(fpoutr,"X_UNIT               %s \n",units[1]);
  fprintf(fpoutr,"Y_UNIT               %s \n",units[1]);

  /*  for (i=9; i<=15; i++)
  fputs(lines[i],fpoutr);  outputs the rest info  into out.rsc file */

  printf("\nEnd of Program\n"); 
/* Close files and end      */
  fclose (fpout);
  fclose (fpoutr);

  exit(0);
}

/***********************/

/* readin: reads in binary input and adds padding */
void readin(char u1[], char u2[])
{
        for (row=pad; row <= xym[2]+pad-1; row++)
	{  
          for (col=pad; col <= xym[1]+pad-1; col++)
          {
            fread(&idem1[row][col],sizeof(short),1,fpin1);
            if (abs(idem1[row][col]) >= background)
              idem1[row][col]=background;
	    else
	      {
		/* see if the input data is in feet if (u1[0] == 'f')*/
                if (strcmp(u1,"feet")==0)
		idem1[row][col]= (short)idem1[row][col]*0.3048; 
		}
	  }
        }
        for (row=pad; row <= xym[4]+pad-1; row++)
	{  
          for (col=pad; col <= xym[3]+pad-1; col++)
          {
            fread(&idem2[row][col],sizeof(short),1,fpin2);
            if (abs(idem2[row][col]) >= background)
              idem2[row][col]=background;
	    else
	      {
                if (strcmp(u2,"feet")==0)
		idem2[row][col]= (short)idem2[row][col]*0.3048; 
              }
	  }
        }

/*     add padding     */

        for (row=1; row <= pad-1; row++)
	{
          for (col=1; col <= xym[1]+2*(pad-1); col++)
            idem1[row][col]= background;
          for (col=1; col <= xym[3]+2*(pad-1); col++)
            idem2[row][col]= background;
        }
        for (row=pad; row <= xym[2]+pad-1; row++)
	{
          for (col=1; col <= pad-1; col++)
            idem1[row][col]= background;
          for (col=xym[1]+pad; col <= xym[1]+2*(pad-1); col++)
            idem1[row][col]= background;
        }
        for (row=pad; row <= xym[4]+pad-1; row++)
	{
          for (col=1; col <= pad-1; col++)
            idem2[row][col]= background;
          for (col=xym[3]+pad; col <= xym[3]+2*(pad-1); col++)
            idem2[row][col]= background;
        }
        for (row=xym[2]+pad; row <= xym[2]+2*(pad-1); row++)
	{
          for (col=1; col <= xym[1]+2*(pad-1); col++)
            idem1[row][col]=background;
        }
        for (row=xym[4]+pad; row <= xym[4]+2*(pad-1); row++)
	{
          for (col=1; col <= xym[3]+2*(pad-1); col++)
            idem2[row][col]= background;
         }
}
    
/* split: gets a line and returns array of constituitive words; 
   spaces are ignored */
int split(char s[])
{
   int i=0; int j; int n=0; char word[maxchar]="";
     while (s[i] != '\n') {
       while ((s[i] == ' ') && (s[i] != '\n')) {
         i++;
       }
       j=i;
       while ((s[i] != ' ') && (s[i] != '\n')) {
         word[i-j]=s[i];
         i++;
       }
       word[i-j]='\0';
       if ((words[n] = (char *)malloc(strlen(word)+1)) == NULL)
         return -1;
       strcpy(words[n++],word);
     } 
    return(n-1);
}    

/* aver: interpolates gaps on the bondaries of SDTS quadrangles */
short aver(long x1, long y1, long x2, long y2)
{
   int i; int n; long value; short a[16];
   
   /*     printf("x1,y1,x2,y2: %ld %ld %ld %ld \n",x1,y1,x2,y2);*/
     a[0]=idem1[x1-1][y1]; 
     a[1]=idem1[x1-1][y1+1]; 
     a[2]=idem1[x1-1][y1-1]; 
     a[3]=idem1[x1][y1]; 
     a[4]=idem1[x1][y1+1]; 
     a[5]=idem1[x1][y1-1]; 
     a[6]=idem1[x1+1][y1]; 
     a[7]=idem1[x1+1][y1+1]; 
     a[8]=idem2[x2-1][y2]; 
     a[9]=idem2[x2-1][y2+1]; 
     a[10]=idem2[x2-1][y2-1]; 
     a[11]=idem2[x2][y2]; 
     a[12]=idem2[x2][y2+1]; 
     a[13]=idem2[x2][y2-1]; 
     a[14]=idem2[x2+1][y2]; 
     a[15]=idem2[x2+1][y2+1]; 
     value=0;
     n=0;
     for(i=0; i<= 15; i++)
     { 
       if (a[i] != background)
       {
         n++;
         value=value+a[i];
	 /*        printf("i,a,value: %d    %d    %d\n",i,a[i],value);*/
       }
     }
     if (n < 5)
       value = background;
     else
     {
       value = value/n;
       printf("coord,value:  %ld %ld %ld %ld   %d \n",x1-pad+1,y1-pad+1,x2-pad+1,y2-pad+1,(short)value);
     }
     return((short)value);
}