/*  Program to flag grid cells that meet channel source criterion. 

  David G Tarboton
  Utah State University  
  Started 2/12/98
    */
  
#include "lsm.h"
 
short **src;
int **arr;

/* function to compute area recursively */
void srcarea(i,j)
  int i,j;
  {                                        
    int in,jn,k;
    if(arr[j][i]==0)
    {
      if(i!=0 && i!=ny-1 && j!=0 && j!=nx-1 && dir[j][i]!= -1)
                 /* not on boundary  */
      {
        arr[j][i]= src[j][i];                               
        for(k=1; k<=8; k++)
        {  in=i+d1[k];
           jn=j+d2[k];
/* test if neighbor drains towards cell excluding boundaryies */
           if(dir[jn][in]>0 && (dir[jn][in]-k==4||dir[jn][in]-k==-4))
             {
                srcarea(in,jn);
                arr[j][i]=arr[j][i]+arr[jn][in];
             }
        }
      }
    }
  }

int source(char *areafile,char *slopefile,char *plenfile,char *dirfile, 
		   char *srcfile, char *elvfile, int ipar,float p[3], double x, double y)
{
  float ndva,ndvs,ndvp,ndvd,emax,ndve,wsum,w[9];
  float **aa, **plen, **selev, **elev;
  double tempf;
  int row, col, filetype,err,i,j,iomax,jomax,bound,ik,jk,k;

/* define directions */
   d1[1]=0; d1[2]= -1; d1[3]= -1; d1[4]= -1; d1[5]=0; d1[6]=1; d1[7]=1; d1[8]=1;
   d2[1]=1; d2[2]=1; d2[3]=0; d2[4]= -1; d2[5]= -1; d2[6]= -1; d2[7]=0; d2[8]=1;

/* read grid files */
   if(ipar <= 3)
   {
	   err=gridread(areafile,(void ***)&aa,RPFLTDTYPE,&nx,&ny,&dx,&dy,
             bndbox,&csize,&ndva,&filetype);
   if(err != 0)goto ERROR;
   }

if(ipar == 2)
{
  err=gridread(slopefile,(void ***)&slope,RPFLTDTYPE,&nx,&ny,&dx,&dy,
             bndbox,&csize,&ndvs,&filetype);
  if(err != 0)goto ERROR; 
}
if(ipar == 3)
{
  err=gridread(plenfile,(void ***)&plen,RPFLTDTYPE,&nx,&ny,&dx,&dy,
             bndbox,&csize,&ndvp,&filetype);
  if(err != 0)goto ERROR; 
}
if(ipar == 4)
{
  err=gridread(elvfile,(void ***)&elev,RPFLTDTYPE,&nx,&ny,&dx,&dy,
             bndbox,&csize,&ndve,&filetype);
  if(err != 0)goto ERROR;  
}

if((src = (short **)matalloc(nx,ny, RPSHRDTYPE)) == NULL)
{
  err=7;
  printf("Could not allocate memory for src\n");
  goto ERROR;
}  
    printf ("src allocated \n");
    

/*  Flag sources  */
for(i=0; i < ny; i++)
  for(j=0; j< nx; j++)
  {
	  src[j][i] = 0;
	  if(ipar == 1)   /*  Area threshold   */
	  {
		  src[j][i]= (aa[j][i] >= p[0]) ? 1 : 0;
	  }
	  else if(ipar == 2)   /*  Slope and area combination   */
	  {
		  if(slope[j][i] > 0.)
			  src[j][i] = (aa[j][i]* pow((double)slope[j][i],(double)p[1]) >= p[0])
			  ? 1: 0;
	  }else if(ipar == 3)  /*  Slope and Length combination   */
	  {
		  if(plen[j][i] > 0.)
			  src[j][i] = (aa[j][i] >= p[0]* pow((double)plen[j][i],(double)p[1])) 
			  ? 1: 0;
	  }
  }
  if(ipar == 4)  /* Peuker and Douglas algorithm  */
  {
/*  Initialize internal cells to 1 for Peuker and Douglas algorithm and smooth  */
if((selev = (float **)matalloc(nx,ny, RPFLTDTYPE)) == NULL)
{
  err=7;
  printf("Could not allocate memory for selev\n");
  goto ERROR;
}

  for(i=0; i <ny; i++)
	 for(j=0; j<nx; j++)
	 {
		  if(i == 0 || i == (ny-1) || j == 0 || j == (nx-1) || elev[j][i] <= ndve)
		  {
			  selev[j][i]=elev[j][i];
		  }
		  else
		  {
		    src[j][i] = 1;
/*		    if(i == 162 && j == nx-2)
		    {
			  printf("here");
			  printf("%g\n",elev[j][i]);
			  printf("%g\n",w[0]);
		    }

		    if(i == 163 && j == 1)
		    {
			  printf("here");
			  printf("%g\n",elev[j][i]);   
			  printf("%g\n",w[0]);
		    }
		    tempf = w[0] * elev[j][i];   */
/*		    printf("%f %f %f\n",w[0],elev[j][i],tempf);  */
		    selev[j][i]=p[0] * elev[j][i];
		    wsum=p[0];
		    if(p[1] > 0.)
              for(k=1; k<=7; k=k+2)
			  {
			    if(elev[j+d2[k]][i+d1[k]] > ndve)
				{
			  	  selev[j][i] += elev[j+d2[k]][i+d1[k]]*p[1];
				  wsum += p[1];
				}				  
			  }
		    if(p[2] > 0.)
              for(k=2; k<=8; k=k+2)
			  {
			    if(elev[j+d2[k]][i+d1[k]] > ndve)
				{
			  	  selev[j][i] += elev[j+d2[k]][i+d1[k]]*p[2];
				  wsum += p[2];
				}				  
			  }
		  }
	  }
    free(elev[0]); free(elev);
	elev = selev;   
	for(i=0; i <ny-1; i++)
	  for(j=0; j<nx-1; j++)
	  {
		  emax= elev[j][i];
          iomax=0;
          jomax=0;           
          bound= 0;  /*  .false.  */
/*  --FIRST PASS FLAG MAX ELEVATION IN GROUP OF FOUR  */
          for(ik=0; ik<2; ik++)
			  for(jk=1-ik; jk < 2; jk++)
			  {
                  if(elev[j+jk][i+ik] > emax)
				  {
                      emax=elev[j+jk][i+ik];
                      iomax=ik;
                      jomax=jk;
				  }
                  if(elev[j+jk][i+ik] <= ndve)
                     bound= 1;  /*  .true.  */
			  }
/*  c---Unflag max pixel */
             src[j+jomax][i+iomax] = 0;
/*  c---Unflag pixels where the group of 4 touches a boundary  */
          if(bound == 1)
		  {
            for(ik=0; ik < 2; ik++)
              for(jk=0; jk< 2; jk++)
			  {
				  src[j+jk][i+ik]=0;
			  }
		  }
/* 		  i.e. unflag flats.  */
            for(ik=0; ik < 2; ik++)
              for(jk=0; jk< 2; jk++)
			  {
                 if(elev[j+jk][i+ik] == emax)src[j+jk][i+ik] = 0;
			  }
	  }
  }
  if(ipar <= 3){free(aa[0]);  free(aa);  } 
  if(ipar == 2){free(slope[0]); free(slope);}
  if(ipar == 3){free(plen[0]); free(plen);}
  if(ipar == 4){free(elev[0]); free(elev);}
/*  Debug write */
/*    err = gridwrite("peukdgrid.asc",(void **)src,RPSHRDTYPE,nx,ny,dx,dy,
                bndbox,csize,-999.,0);  */
/*   Now get directions and compute area's  */

  err=gridread(dirfile,(void ***)&dir,RPSHRDTYPE,&nx,&ny,&dx,&dy,
             bndbox,&csize,&ndvd,&filetype);
  if(err != 0)goto ERROR; 
  if((arr = (int **)matalloc(nx,ny, RPINTDTYPE)) == NULL)
  {
    err=7;
	goto ERROR;
  }  

for(i=0; i < ny; i++)
  for(j=0; j<nx; j++)
  {
	  arr[j][i] = 0;
  }
 	        printf("%lf %lf %lf\n", bndbox[0],bndbox[3],csize);
   if(x > 0. && y > 0.)   /*  Only compute area's for designated location  */
   {
		col= (int)floor((x-bndbox[0])/csize);
		row= (int)floor((bndbox[3]-y)/csize);

		if(row <0 || row > ny || col < 0 || col > nx)
		{
			printf("Given coords out of area\n");
			row=0; col=0;
		}
	}
	if(row > 0 && col > 0)
   {
/* call drainage area subroutine for pixel to zero on  */
      srcarea(row,col);
    }
   err = gridwrite(srcfile,(void **)arr,RPINTDTYPE,nx,ny,dx,dy,
                bndbox,csize,-999.,filetype);
   free(src[0]); free(src);
   free(arr[0]); free(arr);
   free(dir[0]); free(dir);
ERROR:
   return(err);
}

void main(int argc, char **argv)
{
  char infile[MAXLN],areafile[MAXLN],slopefile[MAXLN],plenfile[MAXLN],
	  dirfile[MAXLN],srcfile[MAXLN],elvfile[MAXLN];
  int i,ipar,err;
  float p[3];
  double xr, yr;

  FILE *fp;
  
  if(argc < 2)
  {
     printf("\nUsage:\n"); 
     printf("\n  %s <filename> \n",argv[0]);      
     exit(0);
  }
   sprintf(infile,"%s",argv[1]);
   fp = fopen(infile,"r");
   fscanf(fp, "%s",areafile); eol(fp);
   fscanf(fp, "%s",slopefile); eol(fp);
   fscanf(fp, "%s",plenfile); eol(fp);   
   fscanf(fp, "%s",dirfile); eol(fp);
   fscanf(fp, "%s",srcfile); eol(fp);
   fscanf(fp, "%s",elvfile); eol(fp);
/*  Skip 4 lines  */
   for( i =0; i < 4; i++)eol(fp);
   fscanf(fp, "%d%f%f%f",&ipar,&p[0],&p[1],&p[2]); eol(fp);
   fscanf(fp,"%lf%lf",&xr,&yr);

   if(err=source(areafile,slopefile,plenfile,dirfile,srcfile,elvfile,
	   ipar,p,xr,yr) != 0)
    printf("Source Error %d\n",err);
}