/* fd.c */
/* computes the flow direction and flow accumulation of a grid */

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

#include "fd.h"
#include "grid.h"

int flowdir(Grid *altG, Grid *fdG) {
  
  int rows = altG->nrows,
    cols = altG->ncols,
    nodata = altG->nodataval,
    **alt = altG->tbl,
    **fd = (int **)malloc(sizeof(int *) * rows);

    
  
  if (fd == NULL)
    return -1;

  int i, j, maxdiff, temp, crtalt, crtfd;

  for (i=0; i<rows; i++) {
    
    if((fd[i] = (int *)malloc(sizeof(int) * cols)) == NULL)
      return -1;


    for (j=0; j < cols; j++){
      if ((crtalt = alt[i][j]) != nodata) {
	maxdiff = 0;
	crtfd = 0;
			
	/* Check cells from 1 to 8 */

	/* 1 */
	if (i != 0 && j != 0 & alt[i-1][j-1] != nodata) {
	  if ((temp = crtalt - alt[i-1][j-1]) > maxdiff) {
	    maxdiff = temp;
	    crtfd = 1;
	  }
	}

	/* 2 */
	if (i != 0 && alt[i-1][j] != nodata) {
	  if ((temp = crtalt - alt[i-1][j]) > maxdiff) {
	    maxdiff = temp;
	    crtfd = 2;
	  }
	}
	
	/* 3 */
	if (i != 0 && j != cols - 1 && alt[i-1][j+1] != nodata) { 
	  if ((temp = crtalt - alt[i-1][j+1]) > maxdiff) {
	    maxdiff = temp;
	    crtfd = 3;
	  }
	}

	/* 4 */
	if (j != 0 && alt[i][j-1] != nodata) {
	  if ((temp = crtalt - alt[i][j-1]) > maxdiff) {
	    maxdiff = temp;
	    crtfd = 4;
	  }
	}

	/* 5 */
	if (j != cols - 1 && alt[i][j+1] != nodata) { 
	  if ((temp = crtalt - alt[i][j+1]) > maxdiff) {
	    maxdiff = temp;
	    crtfd = 5;
	  }
	}

	/* 6 */
	if (i != rows - 1 && j != 0 && alt[i+1][j-1] != nodata) { 
	  if ((temp = crtalt - alt[i+1][j-1]) > maxdiff) {
	    maxdiff = temp;
	    crtfd = 6;
	  }
	}

	/* 7 */
	if (i != rows - 1 && alt[i+1][j] != nodata) { 
	  if ((temp = crtalt - alt[i+1][j]) > maxdiff) {
	    maxdiff = temp;
	    crtfd = 7;
	  }
	}


	/* 8 */
	if (i != rows - 1 && j != cols - 1 && alt[i+1][j+1] != nodata) { 
	  if ((temp = crtalt - alt[i+1][j+1]) > maxdiff) {
	    maxdiff = temp;
	    crtfd = 8;
	  }
	}

	fd[i][j] = crtfd;
	
      }
      else
	fd[i][j] = nodata;
    }
  }

  fdG->ncols = cols;
  fdG->nrows = rows;
  fdG->xllcorner = altG->xllcorner;
  fdG->yllcorner = altG->yllcorner;
  fdG->cellsize = altG->cellsize;
  fdG->nodataval = nodata;
  fdG->minval = 0;
  fdG->maxval = 8;
  
  fdG->tbl = fd;

  return 0;

}


int flowacc(Grid *fdG, Grid *faccG) {

 int rows = fdG->nrows,
    cols = fdG->ncols,
    nodata = fdG->nodataval,
    **fdir = fdG->tbl,
    **fa = (int **)malloc(sizeof(int *) * rows);
 
 if (fa == NULL)
   return -1;
 
 int i, j;
 
 for (i = 0; i<rows; i++) {
   
   if ((fa[i] = (int *)malloc(sizeof(int) * cols)) == NULL)
     return -1;
   
   for (j = 0; j<cols; j++) {
     if (fdir[i][j] == nodata)
       fa[i][j] = nodata;
     else
       fa[i][j] = 0;
   }
   
 }


 int tempdir, temp, r, c, max = 1;


 for (i = 0; i<rows; i++) {
    
   for (j = 0; j<cols; j++) {

     if (fa[i][j] == 0) 
       fac(fdir, fa, i, j, rows, cols, &max);
    
     
   }
 }


 
 faccG->ncols = cols;
 faccG->nrows = rows;
 faccG->xllcorner = fdG->xllcorner;
 faccG->yllcorner = fdG->yllcorner;
 faccG->cellsize = fdG->cellsize;
 faccG->nodataval = nodata;
 faccG->minval = 1;
 faccG->maxval = max;

 faccG->tbl = fa;

 return 0;
 

}

int fac(int **fdir, int **acc, int i, int j, int rows, int cols, int *max) {

  /* Check if this node is the first of the chain */

  int sum = 1;

       /* 1 */
       if ( !(i == 0 || j == 0) && fdir[i-1][j-1] == 8) 
	 sum += fac(fdir, acc, i-1, j-1, rows, cols, max);	        
       /* 2 */
       if ( !(i == 0) && fdir[i-1][j] == 7)
	 sum += fac(fdir, acc, i-1, j, rows, cols, max);
       /* 3 */
       if ( !(i == 0 || j == cols - 1) && fdir[i-1][j+1] == 6)
	 sum += fac(fdir, acc, i-1, j+1, rows, cols, max);
       /* 4 */
       if ( !(j == 0) && fdir[i][j-1] == 5)
	 sum += fac(fdir, acc, i, j-1, rows, cols, max);
       /* 5 */
       if ( !(j == cols - 1) && fdir[i][j+1] == 4)
	 sum += fac(fdir, acc, i, j+1, rows, cols, max);
       /* 6 */
       if ( !(i == rows - 1 || j == 0) && fdir[i+1][j-1] == 3)
	 sum += fac(fdir, acc, i+1, j-1, rows, cols, max);
       /* 7 */
       if ( !(i == rows - 1) && fdir[i+1][j] == 2)
	 sum += fac(fdir, acc, i+1, j, rows, cols, max);
       /* 8 */
       if ( !(i == rows - 1 || j == cols -1) && fdir[i+1][j+1] == 1)
	 sum += fac(fdir, acc, i+1, j+1, rows, cols, max);

       
       if (sum > *max)
	 *max = sum;

       acc[i][j] = sum;

       return sum;
  
}