/*****************************************************************************/
/***                                                                       ***/
/***                             process()                                 ***/
/***          Reads in a raster file row by row for processing.		   ***/
/***           Jo Wood, Project ASSIST, V1.0 7th February 1993             ***/
/***                                                                       ***/
/*****************************************************************************/

#include <stdlib.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include "param.h"
#include "nrutil.h"


void process(void)
{

    /*--------------------------------------------------------------------------*/
    /*                              INITIALISE					*/
    /*--------------------------------------------------------------------------*/ 


    DCELL		*row_in,	/* Buffer large enough to hold `wsize' 	*/
			*row_out=NULL,	/* raster rows. When GRASS reads in a	*/
					/* raster row, each element is of type	*/
					/* DCELL	*/
			*window_ptr,	/* Stores local terrain window.		*/
			centre;		/* Elevation of central cell in window.	*/

    CELL                *featrow_out=NULL;   /* store features in CELL */

    struct Cell_head	region;		/* Structure to hold region information	*/

    int			nrows,		/* Will store the current number of 	*/
			ncols,		/* rows and columns in the raster.	*/

			row,col,	/* Counts through each row and column 	*/
					/* of the input raster.			*/

			wind_row,	/* Counts through each row and column	*/
			wind_col,	/* of the local neighbourhood window.	*/

			*index_ptr;	/* Row permutation vector for LU decomp.*/

    double		**normal_ptr,	/* Cross-products matrix.		*/
			*obs_ptr,	/* Observed vector.			*/
			temp;		/* Unused */

    double		*weight_ptr;	/* Weighting matrix for observed values.*/


    /*--------------------------------------------------------------------------*/
    /*               	   GET RASTER AND WINDOW DETAILS			*/
    /*--------------------------------------------------------------------------*/ 

    G_get_window(&region);		/* Fill out the region structure (the	*/
					/* geographical limits etc.)		*/
					
    nrows = G_window_rows();		/* Find out the number of rows and 	*/
    ncols = G_window_cols();		/* columns of the raster.		*/


    if ((region.ew_res/region.ns_res >= 1.01) || /* If EW and NS resolns are	*/
	(region.ns_res/region.ew_res >= 1.01))   /* >1% different, warn user.	*/
    {
	G_warning(_("E-W and N-S grid resolutions are different. Taking average."));
	resoln = (region.ns_res + region.ew_res)/2;
    }
    else
	resoln = region.ns_res;


    /*--------------------------------------------------------------------------*/
    /*              RESERVE MEMORY TO HOLD Z VALUES AND MATRICES		*/
    /*--------------------------------------------------------------------------*/ 

    row_in = (DCELL *) G_malloc(ncols*sizeof(DCELL)*wsize);
					/* Reserve `wsize' rows of memory.	*/

    if(mparam != FEATURE)
      row_out = G_allocate_raster_buf(DCELL_TYPE); /* Initialise output row buffer. 	*/
    else
      featrow_out = G_allocate_raster_buf(CELL_TYPE); /* Initialise output row buffer. 	*/

    window_ptr = (DCELL *) G_malloc(SQR(wsize) * sizeof(DCELL));	
					/* Reserve enough memory for local wind.*/

    weight_ptr = (double *) G_malloc(SQR(wsize) * sizeof(double));	
					/* Reserve enough memory weights matrix.*/

    normal_ptr = dmatrix(0,5,0,5);	/* Allocate memory for 6*6 matrix	*/
    index_ptr  = ivector(0,5);		/* and for 1D vector holding indices	*/
    obs_ptr    = dvector(0,5);		/* and for 1D vector holding observed z */


    /* ---------------------------------------------------------------- */
    /* -            CALCULATE LEAST SQUARES COEFFICIENTS              - */
    /* ---------------------------------------------------------------- */
    
    /*--- Calculate weighting matrix. ---*/

    find_weight(weight_ptr);

    /* Initial coefficients need only be found once since they are 
       constant for any given window size. The only element that 
       changes is the observed vector (RHS of normal equations). */

    /*--- Find normal equations in matrix form. ---*/

    find_normal(normal_ptr,weight_ptr);


    /*--- Apply LU decomposition to normal equations. ---*/

    if (constrained){
	G_ludcmp(normal_ptr,5,index_ptr,&temp);
					   /* To constrain the quadtratic 
					      through the central cell, ignore 
					      the calculations involving the
					      coefficient f. Since these are 
					      all in the last row and column of
					      the matrix, simply redimension.	*/
	/* disp_matrix(normal_ptr,obs_ptr,obs_ptr,5);
	*/
    }
	
    else{
    	G_ludcmp(normal_ptr,6,index_ptr,&temp);
	/* disp_matrix(normal_ptr,obs_ptr,obs_ptr,6);
	*/
    }


    /*--------------------------------------------------------------------------*/
    /*          PROCESS INPUT RASTER AND WRITE OUT RASTER LINE BY LINE		*/
    /*--------------------------------------------------------------------------*/ 

    if(mparam != FEATURE)
      for (wind_row=0; wind_row<EDGE; wind_row++)
	G_put_raster_row(fd_out,row_out,DCELL_TYPE);	/* Write out the edge cells as NULL.	*/
    else
      for (wind_row=0; wind_row<EDGE; wind_row++)
	G_put_raster_row(fd_out,featrow_out,CELL_TYPE);	/* Write out the edge cells as NULL.	*/
    
    for (wind_row=0; wind_row<wsize-1; wind_row++)
    	G_get_raster_row(fd_in,row_in+(wind_row*ncols),wind_row,DCELL_TYPE);
					/* Read in enough of the first rows to	*/
					/* allow window to be examined.		*/

    for (row=EDGE; row<(nrows-EDGE); row++)
    { 
	G_percent(row+1,nrows-EDGE,2);

	G_get_raster_row(fd_in,row_in+((wsize-1)*ncols),row+EDGE,DCELL_TYPE);

	for (col=EDGE; col<(ncols-EDGE); col++)
	{
					/* Find central z value */
	    centre = *(row_in + EDGE*ncols + col);

	    for (wind_row=0; wind_row<wsize; wind_row++)
		for (wind_col=0; wind_col<wsize; wind_col++)

					/* Express all window values relative	*/
					/* to the central elevation.		*/
		    *(window_ptr+(wind_row*wsize)+wind_col) = 
			*(row_in+(wind_row*ncols)+col+wind_col-EDGE) - centre;


    	    /*--- Use LU back substitution to solve normal equations. ---*/
	    find_obs(window_ptr,obs_ptr,weight_ptr);

    	    			/*	disp_wind(window_ptr); 
					disp_matrix(normal_ptr,obs_ptr,obs_ptr,6);
				*/
	
	    if (constrained) 
	    {
	    	G_lubksb(normal_ptr,5,index_ptr,obs_ptr);
		/*
	 	   disp_matrix(normal_ptr,obs_ptr,obs_ptr,5);
		*/
	    }

	    else
	    {
	    	G_lubksb(normal_ptr,6,index_ptr,obs_ptr);
	/*	
	 	  disp_matrix(normal_ptr,obs_ptr,obs_ptr,6);
	*/
		
	    }			

	    /*--- Calculate terrain parameter based on quad. coefficients. ---*/
	    if (mparam == FEATURE)
		*(featrow_out+col) = (CELL)feature(obs_ptr);
	    else
	        *(row_out+col) = param(mparam,obs_ptr);

	    if (mparam == ELEV)
		*(row_out+col) += centre;	/* Add central elevation back */
	}

	if (mparam != FEATURE)
	   G_put_raster_row(fd_out,row_out,DCELL_TYPE); /* Write the row buffer to the output	*/
					/* raster.				*/
	else /* write FEATURE to CELL */
	   G_put_raster_row(fd_out, featrow_out,CELL_TYPE); /* Write the row buffer to the output	*/
					/* raster.				*/         

					/* 'Shuffle' rows down one, and read in	*/
					/*  one new row.			*/
	for (wind_row=0;wind_row<wsize-1;wind_row++)
	   for (col=0; col<ncols; col++)
            *(row_in+(wind_row*ncols)+col) = *(row_in+((wind_row+1)*ncols)+col);
    }

    for (wind_row=0; wind_row<EDGE; wind_row++)
    {
       if (mparam != FEATURE)
 	G_put_raster_row(fd_out,row_out,DCELL_TYPE); /* Write out the edge cells as NULL.*/
       else
        G_put_raster_row(fd_out, featrow_out,CELL_TYPE); /* Write out the edge cells as NULL.*/
    }

    /*--------------------------------------------------------------------------*/
    /*     FREE MEMORY USED TO STORE RASTER ROWS, LOCAL WINDOW AND MATRICES	*/
    /*--------------------------------------------------------------------------*/ 

    G_free(row_in);
    if(mparam != FEATURE)
      G_free(row_out);
    else
      G_free(featrow_out);

    G_free(window_ptr);
    free_dmatrix(normal_ptr,0,5,0,5);
    free_dvector(obs_ptr,0,5);
    free_ivector(index_ptr,0,5);
}