@Filename:              focalcurve.c
 
@Author:                Logicon Geodynamics; Carol Gebhardt, 
                        based on an algorithm developed by:
                        Jacek Blaszczynski, Physical Scientist,
                        BLM NARSC and described in:
                        Blaszczynski, J., 1997.  Landform
                        Characterization with GIS, PE&RS February 
                        1997.

@Non Portable Feature:  None
 
@Assumptions:           None

@Inputs:                Ingrid File Name
                        outgrid File Name
                        Average slope grid File Name
                        Negative average slope grid File Name
                        Positive average slope grid File Name
                        Negative count grid File Name
                        Positive count grid File Name
                        Kernel file
                        Neighborhood shape
                        Neighborhood cell width 
                        Neighborhood cell height
                        Neighborhood outer radius
                        Neighborhood inner radius
                        Neighborhood quadrent direction

@Outputs:               Focalcurve grid
                        Average slope grid
                        Negative average slope grid
                        Positive average slope grid
                        Negative count grid
                        Positive count grid

@Description:           Each cell of the input grid is processed 
                        using the specified neighborhood. For
each 
                        neighborhood cell calculate the rise
(value 
                        of center cell minus value of neighbor
cell) 
                        and run (horizontal distance between the 
                        center cell and the neighbor cell). The
slope 
                        is rise over run. Create the output grids 
                        from the accumulated slope values.

-----------------------------------------------------------------
------------ */
 
/*
-----------------------------------------------------------------
------------
        INCLUDE FILES:
  
-----------------------------------------------------------------
--------- */
 
#include <stdio.h>
#include <fcntl.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <ctype.h>
#include <math.h>
 
 
/*
-----------------------------------------------------------------
------------
        EXTERNAL DECLARATIONS:
  
-----------------------------------------------------------------
--------- */
 
/*
-----------------------------------------------------------------
------------
        DEFINES AND CONSTANTS:
  
-----------------------------------------------------------------
--------- */
#define STATUS_OK    0
#define STATUS_ERROR 1

FILE  *in_file;               /* Input grid */
FILE  *out_file;              /* Output focal curve grid */
FILE  *ave_file;              /* Output average slope grid */
FILE  *neg_ave_file;          /* Output negative average ratio
grid */
FILE  *pos_ave_file;          /* Output positive average ratio
grid */
FILE  *neg_count_file;        /* Output negative count ratio grid
*/
FILE  *pos_count_file;        /* Output positive count ratio grid
*/
FILE  *kernel_file;           /* Input kernel file */
 
/*
-----------------------------------------------------------------
------------
        STRUCTURES:
  
-----------------------------------------------------------------
--------- */
 
/*
 * Command Line Parameters structure
 */
typedef struct _sParms {
    char   ingrid [256];               /* Ingrid File Name */
    char   focalcurve_grid [256];      /* outgrid File Name */
    char   ave_slope_grid [256];       /* average slope grid File
Name */
    char   neg_ave_grid [256];         /* negative average gird
File Name */
    char   pos_ave_grid [256];         /* positive average gird
File Name */
    char   neg_count_grid [256];       /* negative count gird
File Name */
    char   pos_count_grid [256];       /* positive count gird
File Name */
    char   kernel_file [256];          /* kernel file */
    int    shape;                      /* neighborhood shape */
    long   width;                      /* neighborhood cell width
*/
    long   height;                     /* neighborhood cell
height */
    float  outer_radius;               /* neighborhood outer
radius */
    float  inner_radius;               /* neighborhood inner
radius */
    int    direction;                  /* neighborhood quadrent
direction */
    int    min_x;                      /* minimum x neighborhood
offset */
    int    max_x;                      /* maximum x neighborhood
offset */
    int    min_y;                      /* minimum y neighborhood
offset */
    int    max_y;                      /* maximum y neighborhood
offset */
    int    ncols;                      /* number of columns in
input grid */
    int    nrows;                      /* number of rows in input
grid */
    float  xllcorner;                  /* x of lower left corner
*/
    float  yllcorner;                  /* y of lower left corner
*/
    float  cellsize;                   /* grid cell size */
    float  NODATA_value;               /* value for no data */
} Parms;

/*
** Neighborhood structure
*/
typedef struct _sNeigh {
    int    row_num;                    /* neighborhood row number
*/
    float  *row_value;                 /* pointer to neighborhood
row */
} Neigh;

/*
** Kernel structure
*/
typedef struct _sKernel {
    float  *kernel_value;               /* kernel value */
} Kernel;
 
/*
-----------------------------------------------------------------
------------
        INTERNAL DECLARATIONS:
  
-----------------------------------------------------------------
--------- */
 
/*
-----------------------------------------------------------------
------------
        PRIVATE FUNCTION PROTOTYPES:
  
-----------------------------------------------------------------
--------- */

/*
================================================================
OpenFiles */
/*
-----------------------------------------------------------------
------------
@Function:      OpenFiles
 
@Description:  Open the input and output files and create output
headers 
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long OpenFiles (
                                 /*RETURN - Completion Status  */
    Parms     * parms);          /*IN/OUT - input parameters */
 

/*
================================================================
GetKernel */
/*
-----------------------------------------------------------------
------------
@Function:      GetKernel
 
@Description:   Get the kernel file and read into an array 
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long GetKernel (
                                 /*RETURN - Completion Status  */
    Parms     * parms,           /*IN  - input parameters */
    Kernel    * kernel)          /*OUT - kernel array */
;

/* ==============================================================
ProcessGrid */
/*
-----------------------------------------------------------------
------------
@Function:      ProcessGrid
 
@Description:  Process each of the grid cells to create new grids 
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long ProcessGrid (
                                 /*RETURN - Completion Status  */
    Parms     * parms,           /*IN  - input parameters */
    Kernel    * kernel)          /*IN  - kernel array */
;

/* ==========================================================
SetNeighborhood */
/*
-----------------------------------------------------------------
------------
@Function:     SetNeighborhood
 
@Description:  Set the neighborhood boundaries 
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long SetNeighborhood (
                                 /*RETURN - Completion Status  */
    Parms     * parms)           /*IN/OUT - input parameters */
;
 
/* ======================================================
GetCommandLineParms */
/*
-----------------------------------------------------------------
------------
@Function:      GetCommandLineParms
 
@Description:   Get and Validate the Command Line Parameters 
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long GetCommandLineParms (
    int         argc,           /*IN   -Number of Arguments on
command line */
    char     ** argv,           /*IN   -Command line arguments */
    Parms     * parms)          /*  OUT-Ptr to parameter
structure */
;
 
/* ======================================================== 
focalcurve-main()*/
/*
-----------------------------------------------------------------
------------
@Function:      focalcurve-main()
 
@Description:
 
 
@Notes:
 
@Parameters:    see below
 
-----------------------------------------------------------------
------------ */
int main (
                        /* RETURN-Completion Status: 0=Success;
1=Error */
    int argc,           /* N/A -Number of command line arguments
(computed
                                by system, not input by user */
    char ** argv )      /*IN   -Command line arguments ( i.e.
usage )
                           1    <Executable Name>
                           2    <Input File name > */
{
 
  Parms s_parms;                /* Parameter structure */
  Parms *parms=&s_parms;        /* Ptr to Parameter struct */
  Kernel s_kernel;              /* Kernel structure */
  Kernel *kernel = &s_kernel;   /* Ptr to kernel array */
  long  status;                 /* Status variable */

  /*
  ** Verify Command Line
  */
  status = GetCommandLineParms( argc, argv, parms );
  if ( status != STATUS_OK )
  {
     return( status );
  }
 
  /*
  ** Open the input and output files 
  */
  status = OpenFiles( parms );
  if ( status != STATUS_OK )
  {
     return( status );
  }
 
  /*
  ** Set neighborhood boundaries 
  */
  status = SetNeighborhood( parms );
  if ( status != STATUS_OK )
  {
     return( status );
  }

  /*
  ** Get kernel file
  */
  status = GetKernel( parms, kernel );
  if ( status != STATUS_OK )
  {
     return( status );
  }
 
  /*
  ** Process the input cells to create output grids 
  */
  status = ProcessGrid( parms, kernel );
  if ( status != STATUS_OK )
  {
     return( status );
  }

  return( STATUS_OK );
}
 
/*
-----------------------------------------------------------------
------------
        PRIVATE FUNCTION BODIES:
  
-----------------------------------------------------------------
--------- */

/* ==============================================================
ProcessGrid */
/*
-----------------------------------------------------------------
------------
@Function:      ProcessGrid
 
@Description:  Process each of the grid cells to create new grids 
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long ProcessGrid (
    Parms     * parms,           /*IN - input parameters */
    Kernel    * kernel)          /*IN - kernel array */
{
 
  long  status;                 /* Status variable */
  long  status_percent;         /* Percent of rows completed */
  long  old_status;             /* Previous percent of rows
completed */
  long  current_col;            /* Current grid column */
  long  current_row;            /* Current grid row */
  long  neigh_col;              /* Current neighborhood column */
  long  neigh_row;              /* Current neighborhood row */
  long  current_row_index;      /* Array index of current grid
row */
  long  neigh_row_index;        /* Array index of neighborhood
grid row */
  float center_cell;            /* Current center cell value */
  float neigh_cell;             /* Neighborhood cell value */
  long  offset_x;               /* Distance of neighbor from
center */
  long  offset_y;               /* Distance of neighbor from
center */
  Neigh grid_row[100];          /* Array of neighborhood rows */ 
  long  row_index;              /* Index for neighborhood rows */
  long  col_index;              /* Index for neighborhood columns
*/

  long  num_neighbor_cells;     /* Number of cells in
neighborhood */
  long  kernel_index;           /* Kernel location that matches
current */
                                /* cells */
  float total_neg_slope;        /* Total negative slope values */
  float total_pos_slope;        /* Total positive slope values */
  long  num_neg;                /* Number of negative cells */
  long  num_pos;                /* Number of positive cells */
  float total_slope;            /* Total slope values */
  float total_abs_slope;        /* Total absolute slope values */
  float rise;                   /* Rise value between two cells
*/
  float focalcurve_cell;        /* Output cell value */
  float neg_ave_cell;           /* Output cell value */
  float pos_ave_cell;           /* Output cell value */
  float neg_count_cell;         /* Output cell value */
  float pos_count_cell;         /* Output cell value */
  float ave_slope_cell;         /* Output cell value */
  float run;                    /* Distance between cells */
  float slope;                  /* Rise over run */

  /*
  ** Allocate space for neighborhood rows
  */
  for (row_index = 0; row_index < parms->height; row_index++)
  {
     grid_row[row_index].row_value = (float *) malloc 
                                     (sizeof(float) *
parms->ncols);
     if(grid_row[row_index].row_value == 0)
     {
        return ( STATUS_ERROR);
     }
  }

  /*
  ** Initialize the neighborhood
  */
  for (row_index = 0; row_index < parms->height; row_index++)
  {
     grid_row[row_index].row_num = row_index;
  
     for (col_index = 0; col_index < parms->ncols; col_index++)
     {
        fscanf(in_file,"%f ",
&grid_row[row_index].row_value[col_index]);
     }
     fscanf(in_file,"\n");
  }

  /*
  ** Start status meter
  */
  status_percent = 0;
  printf("Focalcurvature: Cells Processed %%%d\n",
status_percent);

  /*
  ** Loop through each input cell and caluclate the output cell
  ** from the neighborhood.
  */
  for (current_row = 0; current_row < parms->nrows;
current_row++)
  {
    /*
    ** Find the current row index (the neighborhood array index
    ** where the current row can be found)
    */
    for (row_index = 0; row_index < parms->height; row_index++)
    {
       if (grid_row[row_index].row_num == current_row)
       {
          current_row_index = row_index;
       }
    }

    /*
    ** Loop through the cells in the current row
    */

    for (current_col = 0; current_col < parms->ncols;
current_col++)
    {

      /*
      ** Initialize the variables
      */
      total_neg_slope = 0;
      total_pos_slope = 0;
      num_neg = 0;
      num_pos = 0;
      total_slope = 0;
      total_abs_slope = 0;
      num_neighbor_cells = 0;
      focalcurve_cell = 0;
      ave_slope_cell = 0;
      neg_count_cell = 0;
      pos_count_cell = 0;
      neg_ave_cell = 0;
      pos_ave_cell = 0;
      kernel_index = 0;

      /*
      ** Get the value of the center cell
      */ 
      center_cell =
grid_row[current_row_index].row_value[current_col];

      if(center_cell != parms->NODATA_value)
      {
        /*
        ** Loop through the neighborhood cells
        */
        for (neigh_row = current_row + parms->min_y; 
             neigh_row <= current_row + parms->max_y;
             neigh_row++)
        {
           /*
           ** Make sure neighborhood cell is on grid
           */
           if (neigh_row > -1 && neigh_row < parms->nrows)
           {
              /*
              ** Find the neighborhood row index (the
neighborhood array index
              ** where the neighborhood row can be found)
              */
              for (row_index = 0; row_index < parms->height;
row_index++)
              {
                 if (grid_row[row_index].row_num == neigh_row)
                 {
                    neigh_row_index = row_index;
                 }
              }

              for (neigh_col = current_col + parms->min_x; 
                   neigh_col <= current_col + parms->max_x;
                   neigh_col++)
              {
                 /*
                 ** Make sure neighborhood cell is on grid
                 */
                 if (neigh_col > -1 && neigh_col < parms->ncols)
                 {
                    neigh_cell =
grid_row[neigh_row_index].row_value[neigh_col];
                    /*
                    ** If center cell then skip it,
                    **  or if kernel = NODATA skip it
                    **  or if neighbor = NODATA skip it
                    */
                    if ((neigh_row == current_row &&
                        neigh_col == current_col) ||
                        kernel->kernel_value[kernel_index] ==
parms->NODATA_value ||
                        neigh_cell == parms->NODATA_value)
                    {
                       kernel_index++;
                    }
                    else
                    {
                       neigh_cell = neigh_cell *
kernel->kernel_value[kernel_index];
                       rise = center_cell - neigh_cell;
  
                       offset_x = fabs(current_col - neigh_col);
                       offset_y = fabs(current_row - neigh_row);
  
                       run = sqrt((parms->cellsize * offset_y) * 
                                  (parms->cellsize * offset_y) +
                                  (parms->cellsize * offset_x) *
                                  (parms->cellsize * offset_x));
   
                       num_neighbor_cells ++;
  
                       if(run != 0)
                       {
                          slope = (rise/run);
                          total_slope = total_slope + slope;
                          total_abs_slope = total_abs_slope +
fabs(slope);
                          if (slope < 0)
                          {
                            num_neg ++;
                            total_neg_slope = total_neg_slope +
slope;
                          }
                          else
                          {
                            num_pos++;
                            total_pos_slope = total_pos_slope +
slope;
                          }
                       }
                       else
                       {
                          printf("Divide by Zero error\n");
                       } 
  
                       kernel_index++;
                    }
                 }
                 else
                 {
                    /*
                    ** Neighborhood cell is not on grid.
                    ** Increment kernel
                    */
                    kernel_index++;
                 }
  
              }  /* end loop through neighborhood columns */
           }  /* end if row in neighborhood */
           else
           {
              /*
              ** Increment by whole row if row off grid
              */
              kernel_index = kernel_index + parms->width;
           }
        }  /* end loop through the neighborhood rows */
  
        /*
        ** Create new grid cell values
        */
        if (num_neighbor_cells > 0)
        {
           focalcurve_cell = total_slope / num_neighbor_cells;
           fprintf(out_file, "%f ", focalcurve_cell);
        }
        else
        {
           fprintf(out_file, "%f ", 0);
        }
  
        if(strncmp(parms->ave_slope_grid, "NO_FILE", 7))
        {
          if(num_neighbor_cells != 0)
          {
            ave_slope_cell = total_abs_slope /
num_neighbor_cells;
            fprintf(ave_file, "%f ", ave_slope_cell);
          }
          else
          {
            fprintf(ave_file, "%f ", 0);
          }
        }
  
        if(strncmp(parms->pos_ave_grid, "NO_FILE", 7))
        {
          if(num_pos == 0)
          {
             pos_ave_cell = 0;
          }
          else
          {
             pos_ave_cell = total_pos_slope / num_pos;
          }
          fprintf(pos_ave_file, "%f ", pos_ave_cell);
        }
  
        if(strncmp(parms->neg_ave_grid, "NO_FILE", 7))
        {
          if(num_neg == 0)
          {
             neg_ave_cell = 0;
          }
          else
          {
             neg_ave_cell = total_neg_slope / num_neg;
          }
          fprintf(neg_ave_file, "%f ", neg_ave_cell);
        }
  
        if(strncmp(parms->neg_count_grid, "NO_FILE", 7))
        {
          fprintf(neg_count_file, "%d ", num_neg);
        }
  
        if(strncmp(parms->pos_count_grid, "NO_FILE", 7))
        {
          fprintf(pos_count_file, "%d ", num_pos);
        }
      }
      else 
      {
        /*
        ** Center is NODATA so output is NODATA
        */
        fprintf(out_file, "%f ", (float)parms->NODATA_value);

        if(strncmp(parms->ave_slope_grid, "NO_FILE", 7))
        {
          fprintf(ave_file, "%f ", (float)parms->NODATA_value);
        }

        if(strncmp(parms->neg_count_grid, "NO_FILE", 7))
        {
          fprintf(neg_count_file, "%f ",
(float)parms->NODATA_value);
        }

        if(strncmp(parms->pos_count_grid, "NO_FILE", 7))
        {
          fprintf(pos_count_file, "%f ",
(float)parms->NODATA_value);
        }

        if(strncmp(parms->neg_ave_grid, "NO_FILE", 7))
        {
          fprintf(neg_ave_file, "%f ",
(float)parms->NODATA_value);
        }

        if(strncmp(parms->pos_ave_grid, "NO_FILE", 7))
        {
          fprintf(pos_ave_file, "%f ",
(float)parms->NODATA_value);
        }
      }
    }   /* end loop through each column  */  

    fprintf(out_file, "\n");

    if(strncmp(parms->ave_slope_grid, "NO_FILE", 7))
    {
      fprintf(ave_file, "\n");
    }

    if(strncmp(parms->neg_count_grid, "NO_FILE", 7))
    {
      fprintf(neg_count_file, "\n");
    }

    if(strncmp(parms->pos_count_grid, "NO_FILE", 7))
    {
      fprintf(pos_count_file, "\n");
    }

    if(strncmp(parms->neg_ave_grid, "NO_FILE", 7))
    {
      fprintf(neg_ave_file, "\n");
    }

    if(strncmp(parms->pos_ave_grid, "NO_FILE", 7))
    {
      fprintf(pos_ave_file, "\n");
    }

    /*
    ** Move neighborhood down one row (if not end of file)
    */
    if (current_row + parms->max_y + 1 < parms->nrows)
    {
       for (row_index = 0; row_index < parms->height;
row_index++)
       {
          if (grid_row[row_index].row_num == current_row +
parms->min_y)
          {
             grid_row[row_index].row_num = current_row +
parms->max_y + 1;
 
             for (col_index = 0; col_index < parms->ncols;
col_index++)
             {
                fscanf(in_file,"%f ", 
                      
&grid_row[row_index].row_value[col_index]);
             }
             fscanf(in_file,"\n");
          }
       }
     }  /* if more rows */

     /*
     ** Print status of processing
     */
     old_status = status_percent;
     status_percent = (int)((((float)current_row + 1) /
                              (float)parms->nrows) * 100);
     if(status_percent != old_status &&
        status_percent % 25 == 0)
     {
        printf("Focalcurvature: Cells Processed %%%d\n",
status_percent);
     }
 
  }  /* loop through rows */  

  /*
  ** Close files
  */
  fclose(in_file);
  fclose(out_file);

  if(strncmp(parms->ave_slope_grid, "NO_FILE", 7))
  {
    fclose(ave_file);
  }

  if(strncmp(parms->neg_count_grid, "NO_FILE", 7))
  {
    fclose(neg_count_file);
  }

  if(strncmp(parms->pos_count_grid, "NO_FILE", 7))
  {
    fclose(pos_count_file);
  }

  if(strncmp(parms->neg_ave_grid, "NO_FILE", 7))
  {
    fclose(neg_ave_file);
  }

  if(strncmp(parms->pos_ave_grid, "NO_FILE", 7))
  {
    fclose(pos_ave_file);
  }

  /*
  ** Free space
  for (row_index = 0; row_index < parms->height; row_index++)
  {
     free (grid_row[row_index].row_value); 
  }
  */

  free (kernel->kernel_value); 

  /*
  ** Everything completed OK
  */
 
  return( STATUS_OK );
}
 
/* ======================================================
GetCommandLineParms */
/*
-----------------------------------------------------------------
------------
@Function:      GetCommandLineParms
 
@Description:   Get the arguments from the command line
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long GetCommandLineParms (
                        /*RETURN-Completion Status  */
  int     argc,         /*IN   - Number of Arguments in Command
Line */
  char ** argv,         /*IN   - Command Line Arguments */
  Parms * parms)        /*OUT  - Ptr to parameter struct */
{
 
  long  status;         /* Status returned from functions */
 
  /*
  ** Make sure there are enough parameters
  */ 
  if ( argc != 15 )
  {
    fprintf(stderr, "Insufficient parameters for focalcurve\n");
    return( STATUS_ERROR );
  }

  /*
  ** Process the Command Line
  */
  strcpy(parms->ingrid, argv[1]);
  strcpy(parms->focalcurve_grid, argv[2]);
  strcpy(parms->ave_slope_grid, argv[3]);
  strcpy(parms->neg_ave_grid, argv[4]);
  strcpy(parms->pos_ave_grid, argv[5]);
  strcpy(parms->neg_count_grid, argv[6]);
  strcpy(parms->pos_count_grid, argv[7]);
  parms->shape = atoi(argv[8]);
  parms->width = atoi(argv[9]);
  parms->height = atoi(argv[10]);
  parms->outer_radius = atof(argv[11]);
  parms->inner_radius = atof(argv[12]);
  parms->direction = atoi(argv[13]);
  strcpy(parms->kernel_file, argv[14]);
 
  return( STATUS_OK );
}
 
/*
================================================================
GetKernel */
/*
-----------------------------------------------------------------
------------
@Function:      GetKernel
 
@Description:   The kernel array has two uses:
                1) If irregular or weighted, the user will supply
                   a kernel file which is read into the kernel
array.
                2) If no kernel file is supplied, the kernel
array is
                   used to mask the shape desired by the user.
NODATA 
                   is put into the array for cells that are not
to be
                   concidered part of the neighborhood. 
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long GetKernel (
                                 /*RETURN-Completion Status  */
    Parms     * parms,           /*IN  - input parameters */
    Kernel     * kernel)          /*OUT - kernel array */
{
 
  long    status;                /* Status returned from
functions */
  int     num_kernels;           /* Number of kernels */
  int     kernel_index;          /* Kernel index */
  int     current_y;             /* Current kernel y location */
  int     current_x;             /* Current kernel x location */
  float   run;                   /* Distance between cells */
  Kernel s_kernel;               /* Kernel structure */
  Kernel * my_kernel = &s_kernel;/* Ptr to kernel structure */

  if (strncmp(parms->kernel_file, "NO_FILE", 7))
  {
     /*
     ** Open the kernel File.
     */
     if( (kernel_file = fopen( parms->kernel_file, "r" )) ==
NULL)
     {
       printf("Cannot open '%s' \n", parms->kernel_file );
       return( STATUS_ERROR );
     }
 
     /*
     ** Read the header 
     */
     fscanf(kernel_file, "%d %d\n", &parms->width,
&parms->height);

     num_kernels = parms->width * parms->height;
   
     /*
     ** Create the space for the kernel array
     */
     my_kernel->kernel_value = (float *) malloc
(sizeof(float)*num_kernels);
     if (my_kernel->kernel_value == 0)
     {
        return ( STATUS_ERROR);
     }

     /*
     ** If shape is WEIGHT then read in kernel as written
     */
     if (parms->shape == 6)
     {
        kernel_index = 0;

        /*
        ** Read file into array
        */
        for (current_y = 0; current_y < parms->height;
current_y++)
        {
           for (current_x = 0; current_x < parms->width;
current_x++)
           {
              fscanf(kernel_file, "%f ",
&my_kernel->kernel_value[kernel_index]);

              kernel_index++;
           }
           fscanf(kernel_file,"\n");
        }
     }
     else
     {
        kernel_index = 0;

        /*
        ** If IRREGULAR and kernel value is zero, change to the
        ** NODATA_value. If kernel value is non-zero, value is 1.
        */
        for (current_y = 0; current_y < parms->height;
current_y++)
        {
           for (current_x = 0; current_x < parms->width;
current_x++)
           {
              fscanf(kernel_file, "%f ",
&my_kernel->kernel_value[kernel_index]);
              if(my_kernel->kernel_value[kernel_index] == 0)
              {
                 my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
              }
              else
              {
                 my_kernel->kernel_value[kernel_index] = 1;
              }
              kernel_index++;
           }
           fscanf(kernel_file,"\n");
        }
     }

     /*
     ** Close file
     */
     fclose(kernel_file);

  }
  else   /* no kernel file */
  {

     num_kernels = parms->width * parms->height;
   
     /*
     ** Create the space for the kernel array
     */
     if (NULL==(my_kernel->kernel_value = (float
*)malloc(num_kernels * sizeof(float))))
     {
        return ( STATUS_ERROR);
     }

    kernel_index = 0;

    /*
    ** If rectangle, all kernel values are 1
    */
    if (parms->shape == 1)
    {
      for (kernel_index = 0; 
           kernel_index < parms->width * parms->height;
kernel_index++)
      {
         my_kernel->kernel_value[kernel_index] = 1.0;
      }
    } 
 
    /*
    ** If circle, kernel values are 1 if within outer radius
    */
    if (parms->shape == 2)
    {
      for (current_y = parms->min_y; current_y <= parms->max_y;
current_y++)
      {
        for (current_x = parms->min_x; current_x <= parms->max_x;
current_x++)
        { 
           run = sqrt(current_x * current_x + current_y *
current_y);
           if (run <= parms->outer_radius)
           {
              my_kernel->kernel_value[kernel_index] = 1.0;
              kernel_index++;
           }
           else
           {
              my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
              kernel_index++;
           }
        }
      }
    } 
 
    /*
    ** If annulus, kernel values are 1 if smaller than outer
radius
    ** and bigger than inner radius
    */
    if (parms->shape == 3)
    {
      for (current_y = parms->min_y; current_y <= parms->max_y;
current_y++)
      {
        for (current_x = parms->min_x; current_x <= parms->max_x;
current_x++)
        { 
            run = sqrt(current_x * current_x + current_y *
current_y);
            if (run <= parms->outer_radius && run >=
parms->inner_radius)
            {
               my_kernel->kernel_value[kernel_index] = 1.0;
               kernel_index++;
            }
            else
            {
               my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
               kernel_index++;
            }
        }
      }
    } 
 
    /*
    ** If wedge, kernel values are 1 if smaller than outer radius
    ** and within the specified quarter circle 
    */
    if (parms->shape == 4)
    {
      for (current_y = parms->min_y; current_y <= parms->max_y;
current_y++)
      {
        for (current_x = parms->min_x; current_x <= parms->max_x;
current_x++)
        { 
            run = sqrt(current_x * current_x + current_y *
current_y);
            if (run <= parms->outer_radius)
            {
              if(parms->direction == 1)  /* east */
              {
                if(current_x > 0)
                {
                   if(fabs(atan(current_y / current_x)) <=
atan(1))
                   {
                      my_kernel->kernel_value[kernel_index] =
1.0;
                      kernel_index++;
                   }
                   else
                   {
                      my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
                      kernel_index++;
                   }
                }
                else
                {
                   my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
                   kernel_index++;
                }
              }
              else if(parms->direction == 2)
              {
                if(current_y >= 0 && current_x != 0)
                {
                   if(fabs(atan(current_y / current_x)) >=
atan(1) &&
                      fabs(atan(current_y / current_x)) <=
atan(-1))
                   {
                      my_kernel->kernel_value[kernel_index] =
1.0;
                      kernel_index++;
                   }
                   else
                   {
                      my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
                      kernel_index++;
                   }
                }
                else
                {
                   my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
                   kernel_index++;
                }
              }
              else if(parms->direction == 3)
              {
                if(current_x < 0)
                {
                   if(fabs(atan(current_y / current_x)) <=
atan(1))
                   {
                      my_kernel->kernel_value[kernel_index] =
1.0;
                      kernel_index++;
                   }
                   else
                   {
                      my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
                      kernel_index++;
                   }
                }
                else
                {
                   my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
                   kernel_index++;
                }
              }
              else if(parms->direction == 4)
              {
                if(current_y <= 0 && current_x != 0)
                {
                  if(fabs(atan(current_y / current_x)) >=
fabs(atan(1)) &&
                     fabs(atan(current_y / current_x)) <=
fabs(atan(-1)))
                  {
                     my_kernel->kernel_value[kernel_index] = 1.0;
                     kernel_index++;
                  }
                  else
                  {
                     my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
                     kernel_index++;
                  }
                }
                else
                {
                   my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
                   kernel_index++;
                }
              }
            }
            else
            {
               my_kernel->kernel_value[kernel_index] =
parms->NODATA_value;
               kernel_index++;
            }
        }
      }
    } 
  }

  kernel->kernel_value = my_kernel->kernel_value;

  /*
  KERNEL DEBUG PRINT
  for (kernel_index = 0; kernel_index < num_kernels;
kernel_index++)
  { 
     printf("kernel %f \n", kernel->kernel_value[kernel_index]); 
  }
  */

  return( STATUS_OK );
}
 
/*
================================================================
OpenFiles */
/*
-----------------------------------------------------------------
------------
@Function:      OpenFiles
 
@Description:  Open the input and output files and create output
headers
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long OpenFiles (
                                 /*RETURN-Completion Status  */
    Parms     * parms)           /*IN  - input parameters */
{
 
  long    status;                /* Status returned from
functions */
  long    index;                 /* String index */

  /*
  ** Open the ingrid File.
  */
  if( (in_file = fopen( parms->ingrid, "r" )) == NULL)
  {
    /*
    ** If can't open file, try converting to lower case first.
    */
    for(index = 0; index <= strlen(parms->ingrid); index++)
    {
       parms->ingrid[index] = tolower(parms->ingrid[index]);
    }

    if( (in_file = fopen( parms->ingrid, "r" )) == NULL)
    {
      printf("Cannot open '%s' \n", parms->ingrid );
      return( STATUS_ERROR );
    }
  }

  /*
  ** Open the outgrid File.
  */
  if( (out_file = fopen( parms->focalcurve_grid, "w" )) == NULL)
  {
    printf("Cannot open '%s' \n", parms->focalcurve_grid );
    return( STATUS_ERROR );
  }

  /*
  ** Create the output file headers
  */
  fscanf(in_file,"ncols         %d\n", &parms->ncols);
  fprintf(out_file,"ncols         %d\n",parms->ncols);

  fscanf(in_file,"nrows         %d\n", &parms->nrows);
  fprintf(out_file,"nrows         %d\n",parms->nrows);
 
  fscanf(in_file,"xllcorner     %f\n", &parms->xllcorner);
  fprintf(out_file,"xllcorner     %.3f\n",parms->xllcorner);
 
  fscanf(in_file,"yllcorner     %f\n", &parms->yllcorner);
  fprintf(out_file,"yllcorner     %.3f\n",parms->yllcorner);
 
  fscanf(in_file,"cellsize      %f\n", &parms->cellsize);
  fprintf(out_file,"cellsize      %.3f\n",parms->cellsize);
 
  fscanf(in_file,"NODATA_value  %f\n", &parms->NODATA_value);
  fprintf(out_file,"NODATA_value  %.3f\n",parms->NODATA_value);

  /*
  ** Open the negative ratio File.
  */
  if (strncmp(parms->neg_ave_grid, "NO_FILE", 7))
  {
    if( (neg_ave_file = fopen( parms->neg_ave_grid, "w" )) ==
NULL)
    {
      printf("Cannot open '%s' \n", parms->neg_ave_grid );
      return( STATUS_ERROR );
    }

    /*
    ** Create the output file headers
    */
    fprintf(neg_ave_file,"ncols         %d\n",parms->ncols);
 
    fprintf(neg_ave_file,"nrows         %d\n",parms->nrows);
 
    fprintf(neg_ave_file,"xllcorner    
%.3f\n",parms->xllcorner);
 
    fprintf(neg_ave_file,"yllcorner    
%.3f\n",parms->yllcorner);
 
    fprintf(neg_ave_file,"cellsize      %.3f\n",parms->cellsize);
 
    fprintf(neg_ave_file,"NODATA_value 
%.3f\n",parms->NODATA_value);
  }

  /*
  ** Open the positive ratio File.
  */
  if (strncmp(parms->pos_ave_grid, "NO_FILE", 7))
  {
    if( (pos_ave_file = fopen( parms->pos_ave_grid, "w" )) ==
NULL)
    {
      printf("Cannot open '%s' \n", parms->pos_ave_grid );
      return( STATUS_ERROR );
    }

    /*
    ** Create the output file headers
    */
    fprintf(pos_ave_file,"ncols         %d\n",parms->ncols);
 
    fprintf(pos_ave_file,"nrows         %d\n",parms->nrows);
 
    fprintf(pos_ave_file,"xllcorner    
%.3f\n",parms->xllcorner);
 
    fprintf(pos_ave_file,"yllcorner    
%.3f\n",parms->yllcorner);
 
    fprintf(pos_ave_file,"cellsize      %.3f\n",parms->cellsize);
 
    fprintf(pos_ave_file,"NODATA_value 
%.3f\n",parms->NODATA_value);
  }

  /*
  ** Open the negative count File.
  */
  if (strncmp(parms->neg_count_grid, "NO_FILE", 7))
  {
    if( (neg_count_file = fopen( parms->neg_count_grid, "w" )) ==
NULL)
    {
      printf("Cannot open '%s' \n", parms->neg_count_grid );
      return( STATUS_ERROR );
    }

    /*
    ** Create the output file headers
    */
    fprintf(neg_count_file,"ncols         %d\n",parms->ncols);
 
    fprintf(neg_count_file,"nrows         %d\n",parms->nrows);
 
    fprintf(neg_count_file,"xllcorner    
%.3f\n",parms->xllcorner);

    fprintf(neg_count_file,"yllcorner    
%.3f\n",parms->yllcorner);
 
    fprintf(neg_count_file,"cellsize     
%.3f\n",parms->cellsize);
 
    fprintf(neg_count_file,"NODATA_value 
%.3f\n",parms->NODATA_value);
  }

  /*
  ** Open the positive count File.
  */
  if (strncmp(parms->pos_count_grid, "NO_FILE", 7))
  {
    if( (pos_count_file = fopen( parms->pos_count_grid, "w" )) ==
NULL)
    {
      printf("Cannot open '%s' \n", parms->pos_count_grid );
      return( STATUS_ERROR );
    }

    /*
    ** Create the output file headers
    */
    fprintf(pos_count_file,"ncols         %d\n",parms->ncols);
 
    fprintf(pos_count_file,"nrows         %d\n",parms->nrows);
 
    fprintf(pos_count_file,"xllcorner    
%.3f\n",parms->xllcorner);
 
    fprintf(pos_count_file,"yllcorner    
%.3f\n",parms->yllcorner);
 
    fprintf(pos_count_file,"cellsize     
%.3f\n",parms->cellsize);
 
    fprintf(pos_count_file,"NODATA_value 
%.3f\n",parms->NODATA_value);
  }

  /*
  ** Open the average slope File.
  */
  if (strncmp(parms->ave_slope_grid, "NO_FILE", 7))
  {
    if( (ave_file = fopen( parms->ave_slope_grid, "w" )) == NULL)
    {
      printf("Cannot open '%s' \n", parms->ave_slope_grid );
      return( STATUS_ERROR );
    }

    /*
    ** Create the output file headers
    */
    fprintf(ave_file,"ncols         %d\n",parms->ncols);
 
    fprintf(ave_file,"nrows         %d\n",parms->nrows);
 
    fprintf(ave_file,"xllcorner     %.3f\n",parms->xllcorner);
 
    fprintf(ave_file,"yllcorner     %.3f\n",parms->yllcorner);
 
    fprintf(ave_file,"cellsize      %.3f\n",parms->cellsize);
 
    fprintf(ave_file,"NODATA_value  %.3f\n",parms->NODATA_value);
  }

  return( STATUS_OK );
}

/* ==========================================================
SetNeighborhood */
/*
-----------------------------------------------------------------
------------
@Function:     SetNeighborhood
 
@Description:  Set the neighborhood boundaries 
 
@Notes:
 
@Parameters:    See below
 
-----------------------------------------------------------------
------------ */
static long SetNeighborhood (
                                 /*RETURN - Completion Status  */
    Parms     * parms)           /*IN/OUT - input parameters */
{
 
  long    status;                /* Status returned from
functions */
  long    x_offset;              /* X offset from center */
  long    y_offset;              /* Y offset from center */

  /*
  ** Set the neighborhood boundaries 
  */
  x_offset = (int)((parms->width -1) / 2);
  y_offset = (int)((parms->height -1) / 2);

  parms->min_x = x_offset * (-1);
  parms->min_y = y_offset * (-1);

  if((parms->width -1) % 2 > 0)
  {
    parms->max_x = x_offset + 1;
  }
  else
  {
    parms->max_x = x_offset;
  }

  if((parms->height -1) % 2 > 0)
  {
    parms->max_y = y_offset + 1;
  }
  else
  {
    parms->max_y = y_offset;
  }

  return( STATUS_OK );
}