/* Program to Map subwatersheds as by junction area threshold  */

/* Created by David G Tarboton  */

/*    Input is read from file wsheds.in in the following order  */
/*    Record 1: File name of pointer file input.                          */
/*    Record 2: File name of area file input (not used).                  */
/*    Record 3: File name of isolated area file output.                   */
/*    Record 4: Threshold for defining watersheds.                        */

/* This program uses Input/Output routines from fileutil.c  */
/*  Uses  nrutil.c from numerical recipes  */
/*  Uses  indexxi.c modified from numerical recipes to sort integers */

#include "lsm.h" 
#define JMAX 22000      
    int **aa, **imatrix(), indexxi();
	void larea();
	short **aw;
   int **jlist, **jun; 
main()
 {                
   char wfile[MAXLN],pfile[MAXLN],afile[MAXLN]; 
   FILE  *fin;
   int i,j,irz,icz,tresh,jk,k,l,n=0,indl[JMAX],err; 
   double bndbox[4],csize;    
/* 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 in data  */
   fin=fopen("wsheds.in","r");
   fscanf(fin, "%s", pfile);     
   fscanf(fin, "%s", afile); 
   fscanf(fin, "%s", wfile); 
   fscanf(fin, "%d", &tresh);   
/*  read pointers */
   err=gridread(pfile,(void ***)&dir,RPSHRDTYPE,&nx,&ny,&dx,&dy,
   bndbox,&csize,&ndv,&filetype);   
/* allocate memory according to numerical recipes conventions */
   jlist = imatrix(0,2,1,JMAX);
   jun=imatrix(0,2,0,7);   
/*   read areas  */
   err=gridread(afile,(void ***)&aa,RPINTDTYPE,&nx,&ny,&dx,&dy,
   bndbox,&csize,&ndv,&filetype);
/*  allocate watershed matrix file  */
   aw= (short **)  matalloc(nx,ny,RPSHRDTYPE);
/* find all junctions.  Store them in an array jlist */
   for(i=0; i<ny; i++)
   {
      for(j=0; j<nx; j++)
       {
/* initialize area  */
        aw[j][i] = 0;
        jk=junctions(i,j,tresh);
/*  Identify junctions where two streams greater than threshold merge
    or one exits the area  */
        if(jk > 1 || (jk > 0 && aa[j][i] == -1))
          for(k=0; k< jk; k++)
          {
            if(n++ <= JMAX){
            for(l=0; l<3; l++)jlist[l][n]=jun[l][k];
/*            printf("%d %d %d %d %d %d %d\n"
                ,jlist[0][n],jlist[1][n],jlist[2][n],jk,i,j);  */
             }
          }
       }
    }
   printf("%d basins located\n",n);
   printf("#  Row  Column   Area (%5.1f x %5.1f m pixels)\n",dx,dy);
   if(n > JMAX){
      printf("WARNING! only %d could be labeled\n",JMAX);
      printf("Limited by parameter JMAX\n");
      }
/*   sort basins, smallest to largest  */
       indexxi(n,jlist[2],indl-1);
/*   proceed through in reverse order numbering them  */
       for(j=1; j<=n; j++)
        {
         i=indl[n-j];
 /* call recursive subroutine to number all pixels in watershed  */
         irz=jlist[0][i];   
         icz=jlist[1][i];
         printf("%d  %d   %d   %d\n",j,irz+1,icz+1,jlist[2][i]);
/*         printf("%d %d\n",aa[icz][irz],dir[icz][irz]);  */
         larea(irz,icz,j);
        }
	err=gridwrite(wfile, (void **)aw,RPSHRDTYPE,nx,ny,dx,dy,
       bndbox,csize, -1., filetype);
   free(aw[0]);
   free(aw);
   free(dir[0]);
   free(dir);
   free(aa[0]);
   free(aa);
}
 
/* function to identify upstream neighbors that exceed area threshold  */
int junctions(i,j,tresh)
 int i,j,tresh;
 {
  int k,jn,in,jf=0;
  for(k=1; k<=8; k++)
  {  
    in=i+d1[k]; 
    jn=j+d2[k];
/* test if neighbor is in domain and drains towards cell and exceeds threshold */
    if(in >= 0 && in < ny && jn >= 0 && jn < nx)
    if(aa[jn][in]>=tresh && (dir[jn][in]-k==4||dir[jn][in]-k==-4))
    {
      jun[0][jf]=in;
      jun[1][jf]=jn;
      jun[2][jf++]=aa[jn][in];
    }
  }
  return jf;
}
     

/* function to label area recursively */
void larea(i,j,lab)
  int i,j,lab;
  {                                        
    int in,jn,k,con=0;
      /* con is a flag that signifies possible contaminatin of area
         due to edge effects  */
    if(aw[j][i]!= (short)lab)
    {
      if(i!=0 && i!=ny-1 && j!=0 && j!=nx-1 && dir[j][i]!= -1)
                 /* not on boundary  */
      {
        aw[j][i]=(short)lab;                               
        for(k=1; k<=8; k++)
        {  in=i+d1[k];
           jn=j+d2[k];
/* test if neighbor drains towards cell excluding boundaries */
           if(dir[jn][in]>=0 && (dir[jn][in]-k==4||dir[jn][in]-k==-4))
             {
                larea(in,jn,lab);
                if(aw[jn][in] < 0)
                  con = -1;
             }
           if(dir[jn][in] < 0)
              con = -1;
        }
        if(con == -1)
           aw[j][i]= -1;
      }
    }
  }
                          
/*   Numerical Recipes P 247  */
void isort2(n,ra,rb)
int n;
int ra[],rb[];
{
	int l,j,ir,i;
	float rrb,rra;

	l=(n >> 1)+1;
	ir=n;
	for (;;) {
		if (l > 1) {
			rra=ra[--l];
			rrb=rb[l];
		} else {
			rra=ra[ir];
			rrb=rb[ir];
			ra[ir]=ra[1];
			rb[ir]=rb[1];
			if (--ir == 1) {
				ra[1]=rra;
				rb[1]=rrb;
				return;
			}
		}
		i=l;
		j=l << 1;
		while (j <= ir)	{
			if (j < ir && ra[j] < ra[j+1]) ++j;
			if (rra < ra[j]) {
				ra[i]=ra[j];
				rb[i]=rb[j];
				j += (i=j);
			}
			else j=ir+1;
		}
		ra[i]=rra;
		rb[i]=rrb;
	}
}

#include <malloc.h>
#include <stdio.h>

void nrerror(error_text)
char error_text[];
{
	void exit();

	fprintf(stderr,"Numerical Recipes run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}



float *vector(nl,nh)
int nl,nh;
{
	float *v;

	v=(float *)malloc((unsigned) (nh-nl+1)*sizeof(float));
	if (!v) nrerror("allocation failure in vector()");
	return v-nl;
}

int *ivector(nl,nh)
int nl,nh;
{
	int *v;

	v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
	if (!v) nrerror("allocation failure in ivector()");
	return v-nl;
}

double *dvector(nl,nh)
int nl,nh;
{
	double *v;

	v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
	if (!v) nrerror("allocation failure in dvector()");
	return v-nl;
}



float **matrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
	int i;
	float **m;

	m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float*));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
		if (!m[i]) nrerror("allocation failure 2 in matrix()");
		m[i] -= ncl;
	}
	return m;
}

double **dmatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
	int i;
	double **m;

	m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
	if (!m) nrerror("allocation failure 1 in dmatrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
		if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
		m[i] -= ncl;
	}
	return m;
}

int **imatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
	int i,**m;

	m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
	if (!m) nrerror("allocation failure 1 in imatrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
		if (!m[i]) nrerror("allocation failure 2 in imatrix()");
		m[i] -= ncl;
	}
	return m;
}



float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
float **a;
int oldrl,oldrh,oldcl,oldch,newrl,newcl;
{
	int i,j;
	float **m;

	m=(float **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(float*));
	if (!m) nrerror("allocation failure in submatrix()");
	m -= newrl;

	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl;

	return m;
}



void free_vector(v,nl,nh)
float *v;
int nl,nh;
{
	free((char*) (v+nl));
}

void free_ivector(v,nl,nh)
int *v,nl,nh;
{
	free((char*) (v+nl));
}

void free_dvector(v,nl,nh)
double *v;
int nl,nh;
{
	free((char*) (v+nl));
}



void free_matrix(m,nrl,nrh,ncl,nch)
float **m;
int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}

void free_dmatrix(m,nrl,nrh,ncl,nch)
double **m;
int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}

void free_imatrix(m,nrl,nrh,ncl,nch)
int **m;
int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}



void free_submatrix(b,nrl,nrh,ncl,nch)
float **b;
int nrl,nrh,ncl,nch;
{
	free((char*) (b+nrl));
}



float **convert_matrix(a,nrl,nrh,ncl,nch)
float *a;
int nrl,nrh,ncl,nch;
{
	int i,j,nrow,ncol;
	float **m;

	nrow=nrh-nrl+1;
	ncol=nch-ncl+1;
	m = (float **) malloc((unsigned) (nrow)*sizeof(float*));
	if (!m) nrerror("allocation failure in convert_matrix()");
	m -= nrl;
	for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl;
	return m;
}



void free_convert_matrix(b,nrl,nrh,ncl,nch)
float **b;
int nrl,nrh,ncl,nch;
{
	free((char*) (b+nrl));
}