/*-
 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993
 * University of Illinois
 * US Army Construction Engineering Research Lab
 * Copyright 1993, H. Mitasova (University of Illinois),
 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)
 *
 * modified by McCauley in August 1995
 * modified by Mitasova in August 1995
 *
 */


/*
 * The interpolation library and interpolation programs, both binary and
 * source is copyrighted, but available without fee for education,
 * research and non-commercial purposes. Users may distribute the binary
 * and source code to third parties provided that the copyright notice and
 * this statement appears on all copies and that no charge is made for
 * such copies.  Any entity wishing to integrate all or part of the source
 * code into a product for  commercial use or resale, should contact the
 * U.S.Army CERL and authors of the software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY. THE
 * U.S.Army CERL or authors SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
 * BY THE USER OF THIS SOFTWARE.
 * 
 * By copying this program, you, the user, agree to abide by the copyright
 * conditions and understandings with respect to any software which is
 * marked with a copyright notice.
 */


#include <stdio.h>
#include <math.h>
extern char *calloc ();
extern int free ();
extern char *malloc ();
extern char *realloc ();

#include "gis.h"
#include "linkm.h"
#include "bitmap.h"
#include "surf.h"
#include "tree.h"
#include "dataquad.h"
#include "interpf.h"

/*
 * flexible, normalized segmented processing surface analysis program with
 * tension and smoothing
 * 
 */

double /* pargr */ ns_res, ew_res, inp_ew_res, inp_ns_res;
int inp_rows, inp_cols;
double x_orig, y_orig;
double inp_x_orig, inp_y_orig;
double dmin, ertre, deltx, delty;
int nsizr, nsizc;
int KMAX2 /* , KMIN, KMAX */ ;

double /* datgr */ *az, *adx, *ady, *adxx, *adyy, *adxy;
double /* error */ ertot, ertre, zminac, zmaxac, zmult;

int total = 0;
int NPOINT = 0;
int OUTRANGE = 0;
int NPT = 0;
int deriv, overlap, cursegm;
double fi;

double DETERM;
int NERROR, cond1, cond2;
char fncdsm[32];
char filnam[10];
char msg[1024];
double fstar2, tfsta2, xmin, xmax, ymin, ymax, zmin, zmax, gmin, gmax, c1min,
 c1max, c2min, c2max;
double dnorm;
double smc;

FCELL *zero_array_cell;
struct interp_params params;

FILE *fdredinp, *fdzout, *fddxout, *fddyout, *fdxxout, *fdyyout, *fd4,
*fxyout;
int fdinp, fdsmooth = -1;

/*
 * x,y,z - input data npoint - number of input data fi - tension parameter
 * b - coef. of int. function a - matrix of system of linear equations az-
 * interpolated values z for output grid adx,ady, ... - estimation of
 * derivatives for output grid nsizr,nsizc - number of rows and columns
 * for output grid xmin ... - coordinates of corners of output grid
 * 
 * subroutines input_data - input of data x,y,z (test function or measured
 * data) iterpolate - interpolation of z-values and derivatives to grid
 * secpar_loop- computation of secondary(morphometric) parameters output -
 * output of gridded data and derivatives/sec.parameters check_at_points -
 * interpolation of z-values to given point x,y
 */

char *input;
char *smooth = NULL;
char *mapset;
char *elev = NULL;
char *slope = NULL;
char *aspect = NULL;
char *pcurv = NULL;
char *tcurv = NULL;
char *mcurv = NULL;
char *maskmap = NULL;
char *redinp = NULL;
int sdisk, disk;
FILE *Tmp_fd_z = NULL;
char *Tmp_file_z = NULL;
FILE *Tmp_fd_dx = NULL;
char *Tmp_file_dx = NULL;
FILE *Tmp_fd_dy = NULL;
char *Tmp_file_dy = NULL;
FILE *Tmp_fd_xx = NULL;
char *Tmp_file_xx = NULL;
FILE *Tmp_fd_yy = NULL;
char *Tmp_file_yy = NULL;
FILE *Tmp_fd_xy = NULL;
char *Tmp_file_xy = NULL;

struct BM *bitmask;
struct Cell_head winhd;
struct Cell_head inphd;
struct Cell_head outhd;
struct Cell_head smhd;

void create_temp_files ();
void clean_fatal_error ();


main (argc, argv)
  int argc;
  char *argv[];
{
  void lntab ();
  void s_table ();
  int segmen ();
  int interpolate ();
  int output ();
  double amax1 ();
  double amin1 ();
  int max1 ();
  int min1 ();
  int per, npmin;
  int ii, i, m1, ret_val;
  char dminchar[200];
  struct FPRange range;
  DCELL cellmin, cellmax;
  FCELL *cellrow;

  struct
  {
    struct Option *input, *elev, *slope, *aspect, *pcurv, *tcurv, *mcurv, *smooth,
    *maskmap, *zmult, *fi, *segmax, *npmin, *res_ew, *res_ns, *overlap;
  } parm;
  struct
  {
    struct Flag *deriv, *iselev, *cprght;
  } flag;


  G_gisinit (argv[0]);

  if (G_get_set_window (&winhd) == -1)
    G_fatal_error ("G_get_set_window failed");

  inp_ew_res = winhd.ew_res;
  inp_ns_res = winhd.ns_res;
  inp_cols = winhd.cols;
  inp_rows = winhd.rows;
  inp_x_orig = winhd.west;
  inp_y_orig = winhd.south;

  parm.input = G_define_option ();
  parm.input->key = "input";
  parm.input->type = TYPE_STRING;
  parm.input->required = YES;
  parm.input->gisprompt = "old,cell,raster";
  parm.input->description = "Name of the input raster file";

  flag.iselev = G_define_flag ();
  flag.iselev->key = 'r';
  flag.iselev->description = "Do zeroes in input map represent elevation?";

  parm.res_ew = G_define_option ();
  parm.res_ew->key = "ew_res";
  parm.res_ew->type = TYPE_DOUBLE;
  parm.res_ew->required = YES;
  parm.res_ew->description = "Desired East-West resolution";

  parm.res_ns = G_define_option ();
  parm.res_ns->key = "ns_res";
  parm.res_ns->type = TYPE_DOUBLE;
  parm.res_ns->required = YES;
  parm.res_ns->description = "Desired North-South resolution";

  parm.elev = G_define_option ();
  parm.elev->key = "elev";
  parm.elev->type = TYPE_STRING;
  parm.elev->required = NO;
  parm.elev->gisprompt = "new,cell,raster";
  parm.elev->description = "Output z-file (elevation)";

  parm.slope = G_define_option ();
  parm.slope->key = "slope";
  parm.slope->type = TYPE_STRING;
  parm.slope->required = NO;
  parm.slope->gisprompt = "new,cell,raster";
  parm.slope->description = "Slope";

  parm.aspect = G_define_option ();
  parm.aspect->key = "aspect";
  parm.aspect->type = TYPE_STRING;
  parm.aspect->required = NO;
  parm.aspect->gisprompt = "new,cell,raster";
  parm.aspect->description = "Aspect";

  parm.pcurv = G_define_option ();
  parm.pcurv->key = "pcurv";
  parm.pcurv->type = TYPE_STRING;
  parm.pcurv->required = NO;
  parm.pcurv->gisprompt = "new,cell,raster";
  parm.pcurv->description = "Profile curvature";

  parm.tcurv = G_define_option ();
  parm.tcurv->key = "tcurv";
  parm.tcurv->type = TYPE_STRING;
  parm.tcurv->required = NO;
  parm.tcurv->gisprompt = "new,cell,raster";
  parm.tcurv->description = "Tangential curvature";

  parm.mcurv = G_define_option ();
  parm.mcurv->key = "mcurv";
  parm.mcurv->type = TYPE_STRING;
  parm.mcurv->required = NO;
  parm.mcurv->gisprompt = "new,cell,raster";
  parm.mcurv->description = "Mean curvature";

  flag.deriv = G_define_flag ();
  flag.deriv->key = 'd';
  flag.deriv->description = "Output partial derivatives instead";

  parm.smooth = G_define_option ();
  parm.smooth->key = "smooth";
  parm.smooth->type = TYPE_STRING;
  parm.smooth->required = NO;
  parm.smooth->gisprompt = "old,cell,raster";
  parm.smooth->description = "Name of the cell file containing smoothing";
  parm.smooth->description = "Name of the cell file containing smoothing";

  parm.maskmap = G_define_option ();
  parm.maskmap->key = "maskmap";
  parm.maskmap->type = TYPE_STRING;
  parm.maskmap->required = NO;
  parm.maskmap->gisprompt = "old,cell,raster";
  parm.maskmap->description = "Name of the raster file used as mask";

  parm.overlap = G_define_option ();
  parm.overlap->key = "overlap";
  parm.overlap->type = TYPE_INTEGER;
  parm.overlap->required = NO;
  parm.overlap->answer = OVERLAP;
  parm.overlap->description = "Rows/columns overlap for segmentation ";

  parm.zmult = G_define_option ();
  parm.zmult->key = "zmult";
  parm.zmult->type = TYPE_DOUBLE;
  parm.zmult->answer = ZMULT;
  parm.zmult->required = NO;
  parm.zmult->description = "Conversion factor for z-values";

  parm.fi = G_define_option ();
  parm.fi->key = "tension";
  parm.fi->type = TYPE_DOUBLE;
  parm.fi->answer = TENSION;
  parm.fi->required = NO;
  parm.fi->description = "Tension";

  flag.cprght = G_define_flag ();
  flag.cprght->key = 'h';
  flag.cprght->answer = 0;
  flag.cprght->description = "Display reference information";


  if (G_parser (argc, argv))
    exit (1);

  if (flag.cprght->answer)
  {
    fprintf (stderr, "\n");
    fprintf (stderr, "\n");
    fprintf (stderr, "Version: GRASS5.0 alpha,  update: August 1995\n");
    fprintf (stderr, "\n");
    fprintf (stderr, "Authors: original version -  H.Mitasova, L.Mitas\n");
    fprintf (stderr, "         GRASS implementation and segmentation: I.Kosinovsky, D.P. Gerdes\n");
    fprintf (stderr, "\n");
    fprintf (stderr, "Methods used in this program are described in the following papers:\n");
    fprintf (stderr, "Mitasova, H., and  Mitas, L., 1993,\n");
    fprintf (stderr, "Interpolation by Regularized Spline with Tension:\n");
    fprintf (stderr, "I. Theory  and  implementation.  Mathematical Geology, 25, 641-655.\n");
    fprintf (stderr, "\n");
    fprintf (stderr, "Mitasova, H., and Hofierka, L., 1993\n");
    fprintf (stderr, "Interpolation by Regularized Spline with Tension:\n");
    fprintf (stderr, "II. Application to terrain modeling and surface   geometry  analysis.\n");
    fprintf (stderr, "Mathematical Geology, 25, 657-669.\n");
    fprintf (stderr, "\n");
    fprintf (stderr, "Mitasova, H., Mitas, L., Brown, W.M., Gerdes, D.P., Kosinovsky, I.,\n");
    fprintf (stderr, "Baker, T., 1995, Modeling spatially and temporally\n");
    fprintf (stderr, "distributed phenomena: New methods and tools for GRASS GIS.\n");
    fprintf (stderr, "International Journal of Geographic Information Systems,V(9), No(4).\n");
    fprintf (stderr, "(special issue on Integration of GIS and Environmental Modeling)\n");
    fprintf (stderr, "\n");
    fprintf (stderr, "The postscript versions of these papers are available via Internet at\n");
    fprintf (stderr, "http://www.cecer.army.mil/grass/viz/rsst.html\n");
    fprintf (stderr, "\n");
    fprintf (stderr, "Please cite these references in publications where the results of this\n");
    fprintf (stderr, "program were used.\n");
    fprintf (stderr, "\n");
    fprintf (stderr, "\n");


  }

  per = 1;			/* flag.per->answer; */
  input = parm.input->answer;
  smooth = parm.smooth->answer;
  maskmap = parm.maskmap->answer;

  elev = parm.elev->answer;
  slope = parm.slope->answer;
  aspect = parm.aspect->answer;
  pcurv = parm.pcurv->answer;
  tcurv = parm.tcurv->answer;
  mcurv = parm.mcurv->answer;

  cond2 = ((pcurv != NULL) || (tcurv != NULL) || (mcurv != NULL));
  cond1 = ((slope != NULL) || (aspect != NULL) || cond2);
  deriv = flag.deriv->answer;

  ertre = 0.1;
  sscanf (parm.res_ew->answer, "%lf", &ew_res);
  sscanf (parm.res_ns->answer, "%lf", &ns_res);
  sscanf (parm.fi->answer, "%lf", &fi);
  sscanf (parm.zmult->answer, "%lf", &zmult);
  sscanf (parm.overlap->answer, "%d", &overlap);
  /*
   * G_set_embedded_null_value_mode(1);
   */
  outhd.ew_res = ew_res;
  outhd.ns_res = ns_res;
  outhd.east = winhd.east;
  outhd.west = winhd.west;
  outhd.north = winhd.north;
  outhd.south = winhd.south;
  outhd.proj = winhd.proj;
  outhd.zone = winhd.zone;
  G_adjust_Cell_head (&outhd, 0, 0);
  ew_res = outhd.ew_res;
  ns_res = outhd.ns_res;
  nsizc = outhd.cols;
  nsizr = outhd.rows;
  disk = nsizc * nsizr * sizeof (int);

  az = G_alloc_vector (nsizc + 1);
  if (!az)
  {
    G_fatal_error ("Not enough memory for az");
  }
  if (cond1)
  {
    adx = G_alloc_vector (nsizc + 1);
    if (!adx)
    {
      G_fatal_error ("Not enough memory for adx");
    }
    ady = G_alloc_vector (nsizc + 1);
    if (!ady)
    {
      G_fatal_error ("Not enough memory for ady");
    }
    if (cond2)
    {
      adxx = G_alloc_vector (nsizc + 1);
      if (!adxx)
      {
	G_fatal_error ("Not enough memory for adxx");
      }
      adyy = G_alloc_vector (nsizc + 1);
      if (!adyy)
      {
	G_fatal_error ("Not enough memory for adyy");
      }
      adxy = G_alloc_vector (nsizc + 1);
      if (!adxy)
      {
	G_fatal_error ("Not enough memory for adxy");
      }
    }
  }
  mapset = NULL;
  if (smooth != NULL)
  {
    mapset = G_find_file ("cell", smooth, "");
    fprintf(stderr,"mapset for smooth [%s], \n",mapset);
    if (mapset == NULL)
    {
      sprintf (msg, "file [%s] not found", smooth);
      G_fatal_error (msg);
    }
    if ((fdsmooth = G_open_cell_old (smooth, mapset)) < 0)
    {
      sprintf (msg, "Cannot open %s", smooth);
      G_fatal_error (msg);
    }
    if (G_get_cellhd(smooth,mapset,&smhd)<0) 
    { 
      sprintf (msg, "[%s]: Cannot read header", smooth); 
      G_fatal_error (msg); 
    } 
    if ((winhd.ew_res!=smhd.ew_res)||(winhd.ns_res!=smhd.ns_res)) 
    { 
      sprintf (msg, "[%s]: Wrong resolution", smooth); 
      G_fatal_error (msg); 
    } 
    if (G_read_fp_range(smooth,mapset,&range) >=0) 
      G_get_fp_range_min_max(&range,&cellmin,&cellmax); 
    if (G_is_f_null_value(&cellmin) || cellmin < 0)  
    { 
      sprintf(msg, "Incorrect value of smoothing - negative or is NULL"); 
      G_fatal_error(msg); 
    } 
  }
  mapset = NULL;
  mapset = G_find_file ("cell", input, "");
  fprintf(stderr, "mapset for input %s\n",mapset);
  if (mapset == NULL)
  {
    sprintf (msg, "file [%s] not found", input);
    G_fatal_error (msg);
  }
  if (G_get_cellhd (input, mapset, &inphd) < 0)
  {
    sprintf (msg, "[%s]: Cannot read header", input);
    G_fatal_error (msg);
  }
  if ((winhd.ew_res != inphd.ew_res) || (winhd.ns_res != inphd.ns_res))
    fprintf (stderr, "Warning: input map resolution differs from current region resolution!\n");
  if ((fdinp = G_open_cell_old (input, mapset)) < 0)
  {
    sprintf (msg, "Cannot open %s", input);
    G_fatal_error (msg);
  }

  sdisk = 0;
  if (elev != NULL)
    sdisk += disk;
  if (slope != NULL)
    sdisk += disk;
  if (aspect != NULL)
    sdisk += disk;
  if (pcurv != NULL)
    sdisk += disk;
  if (tcurv != NULL)
    sdisk += disk;
  if (mcurv != NULL)
    sdisk += disk;
  fprintf (stderr, "\n");
  fprintf (stderr, "Processing all selected output files will require\n");
  fprintf (stderr, "%d bytes of disk space for temp files\n", sdisk);
  fprintf (stderr, "\n");


  fstar2 = fi * fi / 4.;
  tfsta2 = fstar2 + fstar2;
  deltx = winhd.east - winhd.west;
  delty = winhd.north - winhd.south;
  xmin = winhd.west;
  xmax = winhd.east;
  ymin = winhd.south;
  ymax = winhd.north;
  if ( smooth != NULL)
    smc = -9999;
  else
    smc=0.01;
 

  if (G_read_fp_range (input, mapset, &range) >= 0)
  {
    G_get_fp_range_min_max (&range, &cellmin, &cellmax);
  }

  else
  {
    cellrow = G_allocate_f_raster_buf ();
    for (m1 = 0; m1 < inp_rows; m1++)
    {
      ret_val = G_get_f_raster_row (fdinp, cellrow, m1);
      if (ret_val < 0)
      {
	sprintf (msg, "Cannot get row %d (return value = %d)", m1, ret_val);
	G_fatal_error (msg);
      }
      G_row_update_fp_range (cellrow, m1, &range, FCELL_TYPE);
    }
    G_get_fp_range_min_max (&range, &cellmin, &cellmax);
  }
  if (G_is_f_null_value (&cellmin))
  {
    sprintf (msg, "Max value of a cell file is NULL");
    G_fatal_error (msg);
  }

  zmin = (double) cellmin *zmult;
  zmax = (double) cellmax *zmult;

  fprintf (stderr, "zmin=%lf,zmax=%lf\n", zmin, zmax);

  if (fd4 != NULL)
    fprintf (fd4, "deltx,delty %f %f \n", deltx, delty);
  create_temp_files ();

  IL_init_params_2d (&params, NULL, zmult, KMIN, KMAX, maskmap, outhd.rows,
		     outhd.cols, az, adx, ady, adxx, adyy, adxy, fi, 
                     MAXPOINTS, SCIK1, SCIK2, SCIK3,
		     smc, elev, slope, aspect, pcurv, tcurv, mcurv, dmin, 
                     inp_x_orig, inp_y_orig, deriv,
                     Tmp_fd_z, Tmp_fd_dx, Tmp_fd_dy, Tmp_fd_xx, 
                     Tmp_fd_yy, Tmp_fd_xy, NULL);

  IL_init_func_2d (&params, IL_grid_calc_2d, IL_matrix_create, 
                   IL_check_at_points_2d,
		   IL_secpar_loop_2d, IL_crst, IL_crstg, IL_write_temp_2d);

  if (maskmap != NULL)
  {
    fprintf (stderr, "Temporarily changing the region to desired resolution...\n");
    if (G_set_window (&outhd) < 0)
      G_fatal_error ("Cannot set region to output region!");

    bitmask = BM_create (outhd.cols, outhd.rows);

    if (IL_create_bitmask (&params, bitmask) < 0)
      clean_fatal_error ("Cannot create bitmask");
    /* change region to initial region */
    fprintf (stderr, "Changing the region back to initial...\n");
    if (G_set_window (&winhd) < 0)
      G_fatal_error ("Cannot set region to back to initial region!");
  }

  ertot = 0.;
  if (per)
    fprintf (stderr, "Percent complete: ");
  cursegm = 0;

  NPOINT = IL_resample_interp_segments_2d (&params, bitmask, zmin, zmax, &zminac,
    &zmaxac, &gmin, &gmax, &c1min, &c1max, &c2min, &c2max, &ertot, nsizc,
   &dnorm, flag.iselev->answer, overlap, inp_rows, inp_cols, fdsmooth, fdinp,
   ns_res, ew_res, inp_ns_res, inp_ew_res);
  fprintf (stderr, "dnorm in mainc after grid before out1= %lf \n", dnorm);

  if (NPOINT < 0)
    clean_fatal_error ("split_and_interpolate() failed");

  if (fd4 != NULL)
    fprintf (fd4, "max. error found = %f \n", ertot);
  G_free_vector (az);
  if (cond1)
  {
    G_free_vector (adx);
    G_free_vector (ady);
    if (cond2)
    {
      G_free_vector (adxx);
      G_free_vector (adyy);
      G_free_vector (adxy);
    }
  }
  fprintf (stderr, "dnorm in mainc after grid before out2= %lf \n", dnorm);
  if (IL_resample_output_2d (&params, zmin, zmax, zminac, zmaxac, c1min,
			     c1max, c2min, c2max, gmin, gmax, ertot, input, 
                             &dnorm, &outhd, &winhd, smooth, NPOINT) < 0)
    clean_fatal_error ("Cannot write cell files--try increasing cell size");
  free (zero_array_cell);
  if (elev != NULL)
    fclose (Tmp_fd_z);
  if (slope != NULL)
    fclose (Tmp_fd_dx);
  if (aspect != NULL)
    fclose (Tmp_fd_dy);
  if (pcurv != NULL)
    fclose (Tmp_fd_xx);
  if (tcurv != NULL)
    fclose (Tmp_fd_yy);
  if (mcurv != NULL)
    fclose (Tmp_fd_xy);

  if (elev != NULL)
    unlink (Tmp_file_z);
  if (slope != NULL)
    unlink (Tmp_file_dx);
  if (aspect != NULL)
    unlink (Tmp_file_dy);
  if (pcurv != NULL)
    unlink (Tmp_file_xx);
  if (tcurv != NULL)
    unlink (Tmp_file_yy);
  if (mcurv != NULL)
    unlink (Tmp_file_xy);
  if (fd4)
    fclose (fd4);
  G_close_cell (fdinp);
  if (smooth != NULL)
    G_close_cell (fdsmooth);
}


void create_temp_files ()
{
  int i;

  zero_array_cell = (FCELL *) malloc (sizeof (FCELL) * nsizc);
  if (!zero_array_cell)
    G_fatal_error ("Not enough memory for zero_array_cell");

  for (i = 0; i < nsizc; i++)
  {
    zero_array_cell[i] = (FCELL) 0;
  }

  if (elev != NULL)
  {
    Tmp_file_z = G_tempfile ();
    if (NULL == (Tmp_fd_z = fopen (Tmp_file_z, "w+")))
    {
      sprintf (msg, "Can't open temp file [%s] ", Tmp_file_z);
      G_fatal_error (msg);
    }
    for (i = 0; i < nsizr; i++)
    {
      if (!(fwrite (zero_array_cell, sizeof (FCELL), nsizc, Tmp_fd_z)))
	clean_fatal_error ("Not enough disk space -- cannot write files");
    }
  }
  if (slope != NULL)
  {
    Tmp_file_dx = G_tempfile ();
    if (NULL == (Tmp_fd_dx = fopen (Tmp_file_dx, "w+")))
    {
      sprintf (msg, "Can't open temp file [%s] ", Tmp_file_dx);
      clean_fatal_error (msg);
    }
    for (i = 0; i < nsizr; i++)
    {
      if (!(fwrite (zero_array_cell, sizeof (FCELL), nsizc, Tmp_fd_dx)))
	clean_fatal_error ("Not enough disk space -- cannot write files");
    }
  }
  if (aspect != NULL)
  {
    Tmp_file_dy = G_tempfile ();
    if (NULL == (Tmp_fd_dy = fopen (Tmp_file_dy, "w+")))
    {
      sprintf (msg, "Can't open temp file [%s] ", Tmp_file_dy);
      clean_fatal_error (msg);
    }
    for (i = 0; i < nsizr; i++)
    {
      if (!(fwrite (zero_array_cell, sizeof (FCELL), nsizc, Tmp_fd_dy)))
	clean_fatal_error ("Not enough disk space -- cannot write files");
    }
  }

  if (pcurv != NULL)
  {
    Tmp_file_xx = G_tempfile ();
    if (NULL == (Tmp_fd_xx = fopen (Tmp_file_xx, "w+")))
    {
      sprintf (msg, "Can't open temp file [%s] ", Tmp_file_xx);
      clean_fatal_error (msg);
    }
    for (i = 0; i < nsizr; i++)
    {
      if (!(fwrite (zero_array_cell, sizeof (FCELL), nsizc, Tmp_fd_xx)))
	clean_fatal_error ("Not enough disk space -- cannot write files");
    }
  }
  if (tcurv != NULL)
  {
    Tmp_file_yy = G_tempfile ();
    if (NULL == (Tmp_fd_yy = fopen (Tmp_file_yy, "w+")))
    {
      sprintf (msg, "Can't open temp file [%s] ", Tmp_file_yy);
      clean_fatal_error (msg);
    }
    for (i = 0; i < nsizr; i++)
    {
      if (!(fwrite (zero_array_cell, sizeof (FCELL), nsizc, Tmp_fd_yy)))
	clean_fatal_error ("Not enough disk space -- cannot write files");
    }
  }
  if (mcurv != NULL)
  {
    Tmp_file_xy = G_tempfile ();
    if (NULL == (Tmp_fd_xy = fopen (Tmp_file_xy, "w+")))
    {
      sprintf (msg, "Can't open temp file [%s] ", Tmp_file_xy);
      clean_fatal_error (msg);
    }
    for (i = 0; i < nsizr; i++)
    {
      if (!(fwrite (zero_array_cell, sizeof (FCELL), nsizc, Tmp_fd_xy)))
	clean_fatal_error ("Not enough disk space -- cannot write files");
    }
  }
}



void clean_fatal_error (str)
  char *str;
{
  if (Tmp_fd_z)
  {
    fclose (Tmp_fd_z);
    unlink (Tmp_file_z);
  }
  if (Tmp_fd_dx)
  {
    fclose (Tmp_fd_dx);
    unlink (Tmp_file_dx);
  }
  if (Tmp_fd_dy)
  {
    fclose (Tmp_fd_dy);
    unlink (Tmp_file_dy);
  }
  if (Tmp_fd_xx)
  {
    fclose (Tmp_fd_xx);
    unlink (Tmp_file_xx);
  }
  if (Tmp_fd_yy)
  {
    fclose (Tmp_fd_yy);
    unlink (Tmp_file_yy);
  }
  if (Tmp_fd_xy)
  {
    fclose (Tmp_fd_xy);
    unlink (Tmp_file_xy);
  }
  G_fatal_error (str);
}
/*-
 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993
 * University of Illinois
 * US Army Construction Engineering Research Lab  
 * Copyright 1993, H. Mitasova (University of Illinois),
 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)   
 *
 * modified by McCauley in August 1995
 * modified by Mitasova in August 1995  
 *
 */

#define MULT 100000

#include <stdio.h>
#include <math.h>
#include "gis.h"
#include "bitmap.h"
#include "linkm.h"

#include "interpf.h"


/* output cell maps for elevation, aspect, slope and curvatures */

int IL_resample_output_2d (params, zmin, zmax, zminac, zmaxac,
			    c1min, c1max, c2min, c2max,
			    gmin, gmax, ertot, input, dnorm,
			    outhd, winhd, smooth, n_points)

  struct interp_params *params;
  double zmin, zmax;		/* min,max input z-values */
  double zminac, zmaxac;	/* min,max interpolated values */
  double c1min, c1max, c2min, c2max, gmin, gmax;
  double ertot;			/* total interplating func. error */
  char *input;			/* input file name */
  double *dnorm;
  struct Cell_head *outhd;	/* Region with desired resolution */
  struct Cell_head *winhd;	/* Current region */
  char *smooth;
  int n_points;

/*
 * Creates output files as well as history files  and color tables for
 * them.
 */
{
  FCELL *cell1;			/* cell buffer */
  int cf1, cf2, cf3, cf4, cf5, cf6;	/* cell file descriptors */
  int nrows, ncols;		/* current region rows and columns */
  char msg[1024];
  int i, ii;			/* loop counter */
  int zstep;
  char *mapset;
  CELL data1, data2;
  float dat1, dat2;
  struct Colors colors, colors2;
  double value1, value2;
  struct History hist, hist1, hist2, hist3, hist4, hist5;
  struct _Color_Rule_ *rule;
  char *maps, *type;
  int cond1, cond2;

  fprintf (stderr, "dnorm in out at the begin= %lf \n", *dnorm);
  cond2 = ((params->pcurv != NULL) ||
	   (params->tcurv != NULL) ||
	   (params->mcurv != NULL));
  cond1 = ((params->slope != NULL) ||
	   (params->aspect != NULL) || cond2);

  /* change region to output cell file region */
  fprintf (stderr, "Temporarily changing the region to desired resolution...\n");
  if (G_set_window (outhd) < 0)
  {
    fprintf (stderr, "Cannot set region to output region!\n");
    return -1;
  }
  mapset = G_mapset ();

  cell1 = G_allocate_f_raster_buf ();

  if (params->elev != NULL)
  {
    cf1 = G_open_fp_cell_new_uncompressed (params->elev);
    if (cf1 < 0)
    {
      fprintf (stderr, "unable to create raster map %s\n", params->elev);
      return -1;
    }
  }
  if (params->slope != NULL)
  {
    cf2 = G_open_fp_cell_new (params->slope);
    if (cf2 < 0)
    {
      fprintf (stderr, "unable to create raster map %s\n", params->slope);
      return -1;
    }
  }
  if (params->aspect != NULL)
  {
    cf3 = G_open_fp_cell_new (params->aspect);
    if (cf3 < 0)
    {
      fprintf (stderr, "unable to create raster map %s\n", params->aspect);
      return -1;
    }
  }
  if (params->pcurv != NULL)
  {
    cf4 = G_open_fp_cell_new (params->pcurv);
    if (cf4 < 0)
    {
      fprintf (stderr, "unable to create raster map %s\n", params->pcurv);
      return -1;
    }
  }
  if (params->tcurv != NULL)
  {
    cf5 = G_open_fp_cell_new (params->tcurv);
    if (cf5 < 0)
    {
      fprintf (stderr, "unable to create raster map %s\n", params->tcurv);
      return -1;
    }
  }
  if (params->mcurv != NULL)
  {
    cf6 = G_open_fp_cell_new (params->mcurv);
    if (cf6 < 0)
    {
      fprintf (stderr, "unable to create raster map %s\n", params->mcurv);
      return -1;
    }
  }
  nrows = outhd->rows;
  if (nrows != params->nsizr)
  {
    fprintf (stderr, "first change your rows number(%d) to %d!\n",
	     nrows, params->nsizr);
    return -1;
  }
  ncols = outhd->cols;
  if (ncols != params->nsizc)
  {
    fprintf (stderr, "first change your rows number(%d) to %d!\n",
	     ncols, params->nsizc);
    return -1;
  }
  if (params->elev != NULL)
  {
    fseek (params->Tmp_fd_z, 0L, 0);	/* seek to the beginning */
    for (i = 0; i < params->nsizr; i++)
    {
      /* seek to the right row */
      if (fseek (params->Tmp_fd_z, (long)
		 ((params->nsizr - 1 - i) * params->nsizc * sizeof (FCELL)), 0) == -1)
      {
	fprintf (stderr, "cannot fseek to the right spot\n");
	return -1;
      }
      fread (cell1, sizeof (FCELL), params->nsizc, params->Tmp_fd_z);
      if (G_put_f_raster_row (cf1, cell1) < 0)
      {
	fprintf (stderr, "cannot write file\n");
	return -1;
      }
    }
  }
  if (params->slope != NULL)
  {
    fseek (params->Tmp_fd_dx, 0L, 0);	/* seek to the beginning */
    for (i = 0; i < params->nsizr; i++)
    {
      /* seek to the right row */
      if (fseek (params->Tmp_fd_dx, (long)
		 ((params->nsizr - 1 - i) * params->nsizc * sizeof (FCELL)), 0) == -1)
      {
	fprintf (stderr, "cannot fseek to the right spot\n");
	return -1;
      }
      fread (cell1, sizeof (FCELL), params->nsizc, params->Tmp_fd_dx);
      /*
       * for (ii==0;ii<params->nsizc;ii++) { fprintf(stderr,"ii=%d ",ii);
       * fprintf(stderr,"%lf ",cell1[ii]); }
       * fprintf(stderr,"params->nsizc=%d \n",params->nsizc);
       */
      if (G_put_f_raster_row (cf2, cell1) < 0)
      {
	fprintf (stderr, "cannot write file\n");
	return -1;
      }
    }
  }
  if (params->aspect != NULL)
  {
    fseek (params->Tmp_fd_dy, 0L, 0);	/* seek to the beginning */
    for (i = 0; i < params->nsizr; i++)
    {
      /* seek to the right row */
      if (fseek (params->Tmp_fd_dy, (long)
		 ((params->nsizr - 1 - i) * params->nsizc * sizeof (FCELL)), 0) == -1)
      {
	fprintf (stderr, "cannot fseek to the right spot\n");
	return -1;
      }
      fread (cell1, sizeof (FCELL), params->nsizc, params->Tmp_fd_dy);
      if (G_put_f_raster_row (cf3, cell1) < 0)
      {
	fprintf (stderr, "cannot write file\n");
	return -1;
      }
    }
  }
  if (params->pcurv != NULL)
  {
    fseek (params->Tmp_fd_xx, 0L, 0);	/* seek to the beginning */
    for (i = 0; i < params->nsizr; i++)
    {
      /* seek to the right row */
      if (fseek (params->Tmp_fd_xx, (long)
		 ((params->nsizr - 1 - i) * params->nsizc * sizeof (FCELL)), 0) == -1)
      {
	fprintf (stderr, "cannot fseek to the right spot\n");
	return -1;
      }
      fread (cell1, sizeof (FCELL), params->nsizc, params->Tmp_fd_xx);
      if (G_put_f_raster_row (cf4, cell1) < 0)
      {
	fprintf (stderr, "cannot write file\n");
	return -1;
      }
    }
  }
  if (params->tcurv != NULL)
  {
    fseek (params->Tmp_fd_yy, 0L, 0);	/* seek to the beginning */
    for (i = 0; i < params->nsizr; i++)
    {
      /* seek to the right row */
      if (fseek (params->Tmp_fd_yy, (long)
		 ((params->nsizr - 1 - i) * params->nsizc * sizeof (FCELL)), 0) == -1)
      {
	fprintf (stderr, "cannot fseek to the right spot\n");
	return -1;
      }
      fread (cell1, sizeof (FCELL), params->nsizc, params->Tmp_fd_yy);
      if (G_put_f_raster_row (cf5, cell1) < 0)
      {
	fprintf (stderr, "cannot write file\n");
	return -1;
      }
    }
  }
  if (params->mcurv != NULL)
  {
    fseek (params->Tmp_fd_xy, 0L, 0);	/* seek to the beginning */
    for (i = 0; i < params->nsizr; i++)
    {
      /* seek to the right row */
      if (fseek (params->Tmp_fd_xy, (long)
		 ((params->nsizr - 1 - i) * params->nsizc * sizeof (FCELL)), 0) == -1)
      {
	fprintf (stderr, "cannot fseek to the right spot\n");
	return -1;
      }
      fread (cell1, sizeof (FCELL), params->nsizc, params->Tmp_fd_xy);
      if (G_put_f_raster_row (cf6, cell1) < 0)
      {
	fprintf (stderr, "cannot write file\n");
	return -1;
      }
    }
  }

  if (cf1)
    G_close_cell (cf1);
  if (cf2)
    G_close_cell (cf2);
  if (cf3)
    G_close_cell (cf3);
  if (cf4)
    G_close_cell (cf4);
  if (cf5)
    G_close_cell (cf5);
  if (cf6)
    G_close_cell (cf6);


  /* write colormaps and history for output cell files */

  /* colortable for elevations */

  maps = G_find_file ("cell", input, "");

  if (params->elev != NULL)
  {
    if (maps == NULL)
    {
      fprintf (stderr, "file [%s] not found\n", input);
      return -1;
    }
    G_init_colors (&colors2);
    /*
     * G_mark_colors_as_fp(&colors2);
     */

    if (G_read_colors (input, maps, &colors) >= 0)
    {

      if (colors.modular.rules)
      {
	rule = colors.modular.rules;

	while (rule->next)
	  rule = rule->next;

	for (; rule; rule = rule->prev)
	{
	  value1 = rule->low.value * params->zmult;
	  value2 = rule->high.value * params->zmult;
	  G_add_modular_d_raster_color_rule (&value1, rule->low.red,
					     rule->low.grn, rule->low.blu, &value2, rule->high.red, rule->high.grn, rule->high.blu, &colors2);
	}
      }
      if (colors.fixed.rules)
      {
	rule = colors.fixed.rules;

	while (rule->next)
	  rule = rule->next;

	for (; rule; rule = rule->prev)
	{
	  value1 = rule->low.value * params->zmult;
	  value2 = rule->high.value * params->zmult;
	  G_add_d_raster_color_rule (&value1, rule->low.red, rule->low.grn,
				     rule->low.blu, &value2, rule->high.red, rule->high.grn, rule->high.blu, &colors2);
	}
      }
      maps = NULL;
      maps = G_find_file ("cell", params->elev, "");
      if (maps == NULL)
      {
	fprintf (stderr, "file [%s] not found\n", params->elev);
	return -1;
      }
      if (G_write_colors (params->elev, maps, &colors2) < 0)
      {
	fprintf (stderr, "Cannot write color table\n");
	return -1;
      }
      G_quantize_fp_map_range (params->elev, mapset,
			       zminac - 0.5, zmaxac + 0.5,
			   (CELL) (zminac - 0.5), (CELL) (zmaxac + 0.5));

    }
    else
      fprintf (stderr, "No color table for input file -- will not create color table\n");
  }

  /* colortable for slopes */
  if (cond1 & (!params->deriv))
  {

    G_init_colors (&colors);
    G_add_color_rule (0, 255, 255, 255, 2, 255, 255, 0, &colors);
    G_add_color_rule (2, 255, 255, 0, 5, 0, 255, 0, &colors);
    G_add_color_rule (5, 0, 255, 0, 10, 0, 255, 255, &colors);
    G_add_color_rule (10, 0, 255, 255, 15, 0, 0, 255, &colors);
    G_add_color_rule (15, 0, 0, 255, 30, 255, 0, 255, &colors);
    G_add_color_rule (30, 255, 0, 255, 50, 255, 0, 0, &colors);
    G_add_color_rule (50, 255, 0, 0, 90, 0, 0, 0, &colors);


    if (params->slope != NULL)
    {
      maps = NULL;
      maps = G_find_file ("cell", params->slope, "");
      if (maps == NULL)
      {
	fprintf (stderr, "file [%s] not found\n", params->slope);
	return -1;
      }
      G_write_colors (params->slope, maps, &colors);
      G_quantize_fp_map_range (params->slope, mapset, 0., 90., 0, 90);

      type = "raster";
      G_short_history (params->slope, type, &hist1);
      if (params->elev != NULL)
	sprintf (hist1.edhist[0], "The elevation map is %s", params->elev);
      sprintf (hist1.datsrc_1, "raster file %s", input);
      hist1.edlinecnt = 1;

      G_write_history (params->slope, &hist1);

    }


    /* colortable for aspect */

    G_init_colors (&colors);
    G_add_color_rule (0, 255, 255, 255, 0, 255, 255, 255, &colors);
    G_add_color_rule (1, 255, 255, 0, 90, 0, 255, 0, &colors);
    G_add_color_rule (90, 0, 255, 0, 180, 0, 255, 255, &colors);
    G_add_color_rule (180, 0, 255, 255, 270, 255, 0, 0, &colors);
    G_add_color_rule (270, 255, 0, 0, 360, 255, 255, 0, &colors);
    if (params->aspect != NULL)
    {
      maps = NULL;
      maps = G_find_file ("cell", params->aspect, "");
      if (maps == NULL)
      {
	fprintf (stderr, "file [%s] not found\n", params->aspect);
	return -1;
      }
      G_write_colors (params->aspect, maps, &colors);
      G_quantize_fp_map_range (params->aspect, mapset, 0., 360., 0, 360);

      type = "raster";
      G_short_history (params->aspect, type, &hist2);
      if (params->elev != NULL)
	sprintf (hist2.edhist[0], "The elevation map is %s", params->elev);
      sprintf (hist2.datsrc_1, "raster file %s", input);
      hist2.edlinecnt = 1;

      G_write_history (params->aspect, &hist2);

    }


    /* colortable for curvatures */
    if (cond2)
    {
      G_init_colors (&colors);
      dat1 = amin1 (c1min, c2min);
      dat2 = amax1 (c1max, c2max);

      G_add_f_raster_color_rule (dat1, 127, 0, 255,
				 -0.01, 0, 0, 255, &colors);
      G_add_f_raster_color_rule (-0.01, 0, 0, 255,
				 -0.001, 0, 127, 255, &colors);
      G_add_f_raster_color_rule (-0.001, 0, 127, 255,
				 -0.00001, 0, 255, 255, &colors);
      G_add_f_raster_color_rule (-0.00001, 0, 255, 255,
				 0., 200, 255, 200, &colors);
      G_add_f_raster_color_rule (0., 200, 255, 200,
				 0.00001, 255, 255, 0, &colors);
      G_add_f_raster_color_rule (0.00001, 255, 255, 0,
				 0.001, 255, 127, 0, &colors);
      G_add_f_raster_color_rule (0.001, 255, 127, 0,
				 0.01, 255, 0, 0, &colors);
      G_add_f_raster_color_rule (0.01, 255, 0, 0,
				 dat2, 255, 0, 200, &colors);
      maps = NULL;
      if (params->pcurv != NULL)
      {
	maps = G_find_file ("cell", params->pcurv, "");
	if (maps == NULL)
	{
	  fprintf (stderr, "file [%s] not found\n", params->pcurv);
	  return -1;
	}
	G_write_colors (params->pcurv, maps, &colors);
	G_quantize_fp_map_range (params->pcurv, mapset,
				 dat1, dat2,
			     (CELL) (dat1 * MULT), (CELL) (dat2 * MULT));
	type = "raster";
	G_short_history (params->pcurv, type, &hist3);
	if (params->elev != NULL)
	  sprintf (hist3.edhist[0], "The elevation map is %s", params->elev);
	sprintf (hist3.datsrc_1, "raster file %s", input);
	hist3.edlinecnt = 1;

	G_write_history (params->pcurv, &hist3);

      }

      if (params->tcurv != NULL)
      {
	maps = NULL;
	maps = G_find_file ("cell", params->tcurv, "");
	if (maps == NULL)
	{
	  fprintf (stderr, "file [%s] not found\n", params->tcurv);
	  return -1;
	}
	G_write_colors (params->tcurv, maps, &colors);
	G_quantize_fp_map_range (params->tcurv, mapset,
		 dat1, dat2, (CELL) (dat1 * MULT), (CELL) (dat2 * MULT));

	type = "raster";
	G_short_history (params->tcurv, type, &hist4);
	if (params->elev != NULL)
	  sprintf (hist4.edhist[0], "The elevation map is %s", params->elev);
	sprintf (hist4.datsrc_1, "raster file %s", input);
	hist4.edlinecnt = 1;

	G_write_history (params->tcurv, &hist4);

      }

      if (params->mcurv != NULL)
      {
	maps = NULL;
	maps = G_find_file ("cell", params->mcurv, "");
	if (maps == NULL)
	{
	  fprintf (stderr, "file [%s] not found\n", params->mcurv);
	  return -1;
	}
	G_write_colors (params->mcurv, maps, &colors);
	G_quantize_fp_map_range (params->mcurv, mapset,
				 dat1, dat2,
			     (CELL) (dat1 * MULT), (CELL) (dat2 * MULT));

	type = "raster";
	G_short_history (params->mcurv, type, &hist5);
	if (params->elev != NULL)
	  sprintf (hist5.edhist[0], "The elevation map is %s", params->elev);
	sprintf (hist5.datsrc_1, "raster file %s", input);
	hist5.edlinecnt = 1;

	G_write_history (params->mcurv, &hist5);

      }

    }
  }

  if (params->elev != NULL)
  {
    maps = G_find_file ("cell", params->elev, "");
    if (maps == NULL)
    {
      fprintf (stderr, "file [%s] not found \n", params->elev);
      return -1;
    }
    G_short_history (params->elev, "raster", &hist);

/*
should be: (please note that before implementing this, the bug which causes
           dnorm change to 0 hen it leaves ressegm2d.c, should be fixed,
           otherwise this will cause division by 0)
*/
    fprintf (stderr, "dnorm in in out before print= %lf \n", *dnorm);
    if (smooth != NULL) 
      sprintf (hist.edhist[0], "tension=%lf, smoothing=%s", 
               params->fi * 1000. / (*dnorm), smooth); 
    else 
      sprintf (hist.edhist[0], "tension=%lf", params->fi * 1000. /( *dnorm)); 
    sprintf (hist.edhist[1], "dnorm=%lf, zmult=%lf", *dnorm, params->zmult);
    sprintf (hist.edhist[2], "KMAX=%d, KMIN=%d, errtotal=%lf", params->kmax, 
             params->kmin, sqrt (ertot / n_points));
    sprintf (hist.edhist[3], "zmin_data=%lf, zmax_data=%lf", zmin, zmax);
    sprintf (hist.edhist[4], "zmin_int=%lf, zmax_int=%lf", zminac, zmaxac);

    sprintf (hist.datsrc_1, "raster file %s", input);

    hist.edlinecnt = 5;

    G_write_history (params->elev, &hist);
  }

  /* change region to initial region */
  fprintf (stderr, "Changing the region back to initial...\n");
  if (G_set_window (winhd) < 0)
  {
    fprintf (stderr, "Cannot set region to back to initial region!\n");
    return -1;
  }
  return 1;
}
/*-
 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993
 * University of Illinois
 * US Army Construction Engineering Research Lab  
 * Copyright 1993, H. Mitasova (University of Illinois),
 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)   
 *
 * modified by McCauley in August 1995
 * modified by Mitasova in August 1995  
 *
 */

#include <stdio.h>
#include <math.h>
#include "gis.h"

#include "points.h"
#include "dataquad.h"
#include "interpf.h"

int input_data ();
int write_zeros ();

int IL_resample_interp_segments_2d (params, bitmask, zmin, zmax, 
                                    zminac, zmaxac, gmin, gmax,
                                    c1min, c1max, c2min, c2max, ertot, 
                                    offset1, dnorm, iselev, overlap,
                                    inp_rows, inp_cols, fdsmooth, fdinp, 
                                    ns_res, ew_res, inp_ns_res, inp_ew_res)

  struct interp_params *params;
  struct BM *bitmask;		/* bitmask */
  double zmin, zmax;		/* min and max input z-values */
  double *zminac, *zmaxac,	/* min and max interp. z-values */
  *gmin, *gmax,			/* min and max inperp. slope val. */
  *c1min, *c1max, *c2min, *c2max;	/* min and max interp. curv. val. */
  double *ertot;		/* total interplating func. error */
  int offset1;			/* offset for temp file writing */
  double *dnorm;
  int iselev;
  int overlap;
  int inp_rows, inp_cols;
  int fdsmooth, fdinp;
  double ns_res, ew_res;
  double inp_ns_res, inp_ew_res;
{

  int i, j, k, l, m, m1, i1;	/* loop coounters */
  int cursegm = 0;
  int new_comp = 0;
  int n_rows, n_cols, inp_r, inp_c;
  double x_or, y_or, xm, ym;
  static int first = 1, new_first = 1;
  char msg[1024];
  double **matrix, **new_matrix, *b = NULL;
  int *indx, *new_indx;
  static struct fcell_triple *in_points = NULL;	/* input points */
  int inp_check_rows, inp_check_cols,	/* total input rows/cols */
   out_check_rows, out_check_cols;	/* total output rows/cols */
  int first_row, last_row;	/* first and last input row of segment */
  int first_col, last_col;	/* first and last input col of segment */
  int num, prev;
  int div;			/* number of divides */
  int rem_inp_row, rem_inp_col,	/* input rows/cols remainders */
   rem_out_row, rem_out_col;	/* output rows/cols remainders */
  int inp_seg_r, inp_seg_c,	/* # of input rows/cols in segment */
   out_seg_r, out_seg_c;	/* # of output rows/cols in segment */
  int ngstc, nszc		/* first and last output col of the
				 * segment */
  ,ngstr, nszr;			/* first and last output row of the
				 * segment */
  int index;			/* index for input data */
  int c, r;
  int overlap1;
  int p_size;
  struct quaddata *data;
  double xmax, xmin, ymax, ymin;
  int totsegm;			/* total number of segments */
  int total_points = 0;


  xmin = params->x_orig;
  ymin = params->y_orig;
  xmax = xmin + ew_res * params->nsizc;
  ymax = ymin + ns_res * params->nsizr;
  prev = inp_rows * inp_cols;
  if (prev <= params->kmax)
    div = 1;			/* no segmentation */

  else
  {				/* find the number of divides */
    for (i = 2;; i++)
    {
      c = inp_cols / i;
      r = inp_rows / i;
      num = c * r;
      if (num < params->kmin)
      {
	if (((params->kmin - num) > (prev + 1 - params->kmax)) && (prev + 1 < params->KMAX2))
	{
	  div = i - 1;
	  break;
	}
	else
	{
	  div = i;
	  break;
	}
      }
      if ((num > params->kmin) && (num + 1 < params->kmax))
      {
	div = i;
	break;
      }
      prev = num;
    }
  }
  out_seg_r = params->nsizr / div;	/* output rows per segment */
  out_seg_c = params->nsizc / div;	/* output cols per segment */
  inp_seg_r = inp_rows / div;	/* input rows per segment */
  inp_seg_c = inp_cols / div;	/* input rows per segment */
  rem_out_col = params->nsizc % div;
  rem_out_row = params->nsizr % div;
  overlap1 = min1 (overlap, inp_seg_c - 1);
  overlap1 = min1 (overlap1, inp_seg_r - 1);
  out_check_rows = 0;
  out_check_cols = 0;
  inp_check_rows = 0;
  inp_check_cols = 0;

  if (div == 1)
  {
    p_size = inp_seg_c * inp_seg_r;
  }
  else
  {
    p_size = (overlap1 * 2 + inp_seg_c) * (overlap1 * 2 + inp_seg_r);
  }
  if (!in_points)
  {
    if (!(in_points = (struct fcell_triple *) G_malloc (sizeof (struct fcell_triple) * p_size * div)))
    {
      fprintf (stderr, "Cannot allocate memory for in_points\n");
      return -1;
    }
  }

  *dnorm = sqrt (((xmax-xmin)*(ymax-ymin)*p_size)/(inp_rows*inp_cols));
  params->fi = params->fi*(*dnorm)/1000.;

  if (div == 1)
  {				/* no segmentation */
    totsegm = 1;
    cursegm = 1;

    input_data (params, 1, inp_rows, in_points, fdsmooth, fdinp, inp_rows,
		inp_cols, zmin, inp_ns_res, inp_ew_res);

    x_or = 0.;
    y_or = 0.;
    xm = params->nsizc * ew_res;
    ym = params->nsizr * ns_res;

    data = (struct quaddata *) quad_data_new (x_or, y_or, xm, ym,
			 params->nsizr, params->nsizc, 0, params->KMAX2);
    m1 = 0;
    for (k = 1; k <= p_size; k++)
    {
      if (!G_is_f_null_value (&(in_points[k - 1].z)))
      {
	if ((iselev) || (in_points[k - 1].z + zmin != 0.))
	{
	  data->points[m1].x = in_points[k - 1].x / (*dnorm);
	  data->points[m1].y = in_points[k - 1].y / (*dnorm);
	  data->points[m1].z = (double) (in_points[k - 1].z) / (*dnorm);
	  data->points[m1].smooth = in_points[k - 1].smooth;
	  m1++;
	}
      }
    }
    data->n_points = m1;
    total_points = m1;
    if (!(indx = G_alloc_ivector (params->KMAX2 + 1)))
    {
      fprintf (stderr, "Cannot allocate memory for indx\n");
      return -1;
    }
    if (!(matrix = G_alloc_matrix (params->KMAX2 + 1, params->KMAX2 + 1)))
    {
      fprintf (stderr, "Cannot allocate memory for matrix\n");
      return -1;
    }
    if (!(b = G_alloc_vector (params->KMAX2 + 2)))
    {
      fprintf (stderr, "Cannot allocate memory for b\n");
      return -1;
    }

    if (params->matrix_create (params, data->points, m1, matrix, indx) < 0)
      return -1;
    for (i = 0; i < m1; i++)
    {
      b[i + 1] = data->points[i].z;
    }
    b[0] = 0.;
    G_lubksb (matrix, m1 + 1, indx, b);

    params->check_points (params, data, b, ertot, zmin, *dnorm);

    if (params->grid_calc (params, data, bitmask,
			   zmin, zmax, zminac, zmaxac, gmin, gmax,
	       c1min, c1max, c2min, c2max, ertot, b, offset1, *dnorm) < 0)
    {
      fprintf (stderr, "interpolation failed\n");
      return -1;
    }
    else
    {
      if (totsegm != 0)
      {
	G_percent (cursegm, totsegm, 1);
      }
      /*
       * if (b) G_free_vector(b); if (matrix) G_free_matrix(matrix); if
       * (indx) G_free_ivector(indx);
       */
      fprintf (stderr, "dnorm in ressegm after grid before out= %lf \n", *dnorm);
      return total_points;
    }
  }

  out_seg_r = params->nsizr / div;	/* output rows per segment */
  out_seg_c = params->nsizc / div;	/* output cols per segment */
  inp_seg_r = inp_rows / div;	/* input rows per segment */
  inp_seg_c = inp_cols / div;	/* input rows per segment */
  rem_out_col = params->nsizc % div;
  rem_out_row = params->nsizr % div;
  overlap1 = min1 (overlap, inp_seg_c - 1);
  overlap1 = min1 (overlap1, inp_seg_r - 1);
  out_check_rows = 0;
  out_check_cols = 0;
  inp_check_rows = 0;
  inp_check_cols = 0;

  totsegm = div * div;

  /* set up a segment */
  for (i = 1; i <= div; i++)
  {				/* input and output rows */
    if (i <= div - rem_out_row)
      n_rows = out_seg_r;
    else
      n_rows = out_seg_r + 1;
    inp_r = inp_seg_r;
    out_check_cols = 0;
    inp_check_cols = 0;
    ngstr = out_check_rows + 1;	/* first output row of the segment */
    nszr = ngstr + n_rows - 1;	/* last output row of the segment */
    y_or = (ngstr - 1) * ns_res;/* y origin of the segment */
    /*
     * Calculating input starting and ending rows and columns of this
     * segment
     */
    first_row = (int) (y_or / inp_ns_res) + 1;
    if (first_row > overlap1)
    {
      first_row -= overlap1;	/* middle */
      last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
      if (last_row > inp_rows)
      {
	first_row -= (last_row - inp_rows);	/* bottom */
	last_row = inp_rows;
      }
    }
    else
    {
      first_row = 1;		/* top */
      last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
    }
    if ((last_row > inp_rows) || (first_row < 1))
    {
      fprintf (stderr, "Row overlap too large!\n");
      return -1;
    }
    input_data (params, first_row, last_row, in_points, fdsmooth, fdinp,
		inp_rows, inp_cols, zmin, inp_ns_res, inp_ew_res);

    for (j = 1; j <= div; j++)
    {				/* input and output cols */
      if (j <= div - rem_out_col)
	n_cols = out_seg_c;
      else
	n_cols = out_seg_c + 1;
      inp_c = inp_seg_c;

      ngstc = out_check_cols + 1;	/* first output col of the segment */
      nszc = ngstc + n_cols - 1;/* last output col of the segment */
      x_or = (ngstc - 1) * ew_res;	/* x origin of the segment */

      first_col = (int) (x_or / inp_ew_res) + 1;
      if (first_col > overlap1)
      {
	first_col -= overlap1;	/* middle */
	last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
	if (last_col > inp_cols)
	{
	  first_col -= (last_col - inp_cols);	/* right */
	  last_col = inp_cols;
	}
      }
      else
      {
	first_col = 1;		/* left */
	last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
      }
      if ((last_col > inp_cols) || (first_col < 1))
      {
	fprintf (stderr, "Column overlap too large!\n");
	return -1;
      }
      m = 0;
      /* Getting points for interpolation (translated) */

      xm = nszc * ew_res;
      ym = nszr * ns_res;
      data = (struct quaddata *) quad_data_new (x_or, y_or, xm, ym,
		   nszr - ngstr + 1, nszc - ngstc + 1, 0, params->KMAX2);
      new_comp = 0;

      for (k = 0; k <= last_row - first_row; k++)
      {
	for (l = first_col - 1; l < last_col; l++)
	{
	  index = k * inp_cols + l;
	  if (!G_is_f_null_value (&(in_points[k - 1].z)))
	  {
	    if ((iselev) || (in_points[index].z + zmin != 0.))
	    {
	      /* if the point is inside the segment (not overlapping) */
	      if ((in_points[index].x - x_or >= 0) &&
		  (in_points[index].y - y_or >= 0) &&
		  ((nszc - 1) * ew_res - in_points[index].x >= 0) &&
		  ((nszr - 1) * ns_res - in_points[index].y >= 0))
		total_points += 1;
	      data->points[m].x = (in_points[index].x - x_or) / (*dnorm);
	      data->points[m].y = (in_points[index].y - y_or) / (*dnorm);
	      data->points[m].z = (double) (in_points[index].z) / (*dnorm);
	      data->points[m].smooth = in_points[index].smooth;
	      m++;
	    }
	    else
	      new_comp = 1;
	  }
	  else
	    new_comp = 1;
	  /*
	   * fprintf(stderr,"%lf,%lf,%lf
	   * zmin=%lf\n",in_points[index].x,in_points[index].y,in_points[in
	   * dex].z,zmin);
	   */
	}
      }
      if (m <= params->KMAX2)
	data->n_points = m;
      else
	data->n_points = params->KMAX2;
      out_check_cols += n_cols;
      inp_check_cols += inp_c;
      cursegm = (i - 1) * div + j - 1;

      /* show before to catch 0% */
      if (totsegm != 0)
      {
	G_percent (cursegm, totsegm, 1);
      }
      if (m == 0)
      {
	/*
	 * fprintf(stderr,"Warning: segment with zero points encountered,
	 * insrease overlap\n");
	 */
	write_zeros (params, data, offset1);
      }
      else
      {

	if (new_comp)
	{
	  if (new_first)
	  {
	    new_first = 0;
	    if (!b)
	    {
	      if (!(b = G_alloc_vector (params->KMAX2 + 2)))
	      {
		fprintf (stderr, "Cannot allocate memory for b\n");
		return -1;
	      }
	    }
	    if (!(new_indx = G_alloc_ivector (params->KMAX2 + 1)))
	    {
	      fprintf (stderr, "Cannot allocate memory for new_indx\n");
	      return -1;
	    }
	    if (!(new_matrix = G_alloc_matrix (params->KMAX2 + 1, params->KMAX2 + 1)))
	    {
	      fprintf (stderr, "Cannot allocate memory for new_matrix\n");
	      return -1;
	    }
	  }
	  if (params->matrix_create (params, data->points, data->n_points,
				     new_matrix, new_indx) < 0)
	    return -1;
	  for (i1 = 0; i1 < m; i1++)
	  {
	    b[i1 + 1] = data->points[i1].z;
	  }
	  b[0] = 0.;
	  G_lubksb (new_matrix, data->n_points + 1, new_indx, b);

	  params->check_points (params, data, b, ertot, zmin, *dnorm);

	  if (params->grid_calc (params, data, bitmask,
				 zmin, zmax, zminac, zmaxac, gmin, gmax,
	       c1min, c1max, c2min, c2max, ertot, b, offset1, *dnorm) < 0)
	  {

	    fprintf (stderr, "interpolate() failed\n");
	    return -1;
	  }
	}
	else
	{
	  if (first)
	  {
	    first = 0;
	    if (!b)
	    {
	      if (!(b = G_alloc_vector (params->KMAX2 + 2)))
	      {
		fprintf (stderr, "Cannot allocate memory for b\n");
		return -1;
	      }
	    }
	    if (!(indx = G_alloc_ivector (params->KMAX2 + 1)))
	    {
	      fprintf (stderr, "Cannot allocate memory for indx\n");
	      return -1;
	    }
	    if (!(matrix = G_alloc_matrix (params->KMAX2 + 1, params->KMAX2 + 1)))
	    {
	      fprintf (stderr, "Cannot allocate memory for matrix\n");
	      return -1;
	    }
	    if (params->matrix_create (params, data->points, data->n_points,
				       matrix, indx) < 0)
	      return -1;
	  }
	  for (i1 = 0; i1 < m; i1++)
	    b[i1 + 1] = data->points[i1].z;
	  b[0] = 0.;
	  G_lubksb (matrix, data->n_points + 1, indx, b);

	  params->check_points (params, data, b, ertot, zmin, *dnorm);

	  if (params->grid_calc (params, data, bitmask,
				 zmin, zmax, zminac, zmaxac, gmin, gmax,
	       c1min, c1max, c2min, c2max, ertot, b, offset1, *dnorm) < 0)
	  {

	    fprintf (stderr, "interpolate() failed\n");
	    return -1;
	  }
	}
      }
      if (data)
      {
	free (data->points);
	free (data);
      }
      /*
       * cursegm++;
       */
      /* show after to catch 100% */
      if (totsegm != 0)
      {
	G_percent (cursegm, totsegm, 1);
      }
    }
    inp_check_rows += inp_r;
    out_check_rows += n_rows;
  }
  /*
   * if (b) G_free_vector(b); if (indx) G_free_ivector(indx); if (matrix)
   * G_free_matrix(matrix);
   */
  fprintf (stderr, "dnorm in ressegm after grid before out2= %lf \n", *dnorm);
  return total_points;
}



/* input of data for interpolation and smoothing parameters */

input_data (params, first_row, last_row, points, fdsmooth, fdinp,
	    inp_rows, inp_cols, zmin, inp_ns_res, inp_ew_res)
  struct interp_params *params;
  int first_row, last_row;
  struct fcell_triple *points;
  int fdsmooth, fdinp;
  int inp_rows, inp_cols;
  double zmin;
  double inp_ns_res, inp_ew_res;
{
  double x, y, sm;		/* input data and smoothing */
  FCELL z;
  int ddisk = 0, sddisk = 0;	/* amount of disk space needed */
  double deltx, delty;
  char buf[1024];
  int m1, m2;			/* loop counters */
  static int first_time = 1;
  int ret_val, ret_val1;	/* return values of G_get_map_row */
  static FCELL *cellinp = NULL;	/* cell buffer for input data */
  static FCELL *cellsmooth = NULL;	/* cell buffer for smoothing */
  int ind = 0;


  if (!cellinp)
    cellinp = G_allocate_f_raster_buf ();
  if (!cellsmooth)
    cellsmooth = G_allocate_f_raster_buf ();

  for (m1 = 0; m1 <= last_row - first_row; m1++)
  {
    ret_val = G_get_f_raster_row (fdinp, cellinp, inp_rows - m1 - first_row);
    if (ret_val < 0)
    {
      fprintf (stderr, "Cannot get row %d (return value = %d)\n", m1, ret_val);
      return -1;
    }
    if (fdsmooth >= 0)
    {
      ret_val1 = G_get_f_raster_row (fdsmooth, cellsmooth, inp_rows - m1 - first_row);
      if (ret_val1 < 0)
      {
	fprintf (stderr, "Cannot get smoothing row\n");
      }
    }
    y = params->y_orig + (m1 + first_row - 1 + 0.5) * inp_ns_res;
    for (m2 = 0; m2 < inp_cols; m2++)
    {
      x = params->x_orig + (m2 + 0.5) * inp_ew_res;
      /*
       * z = cellinp[m2]*params->zmult;
       */
      if (fdsmooth >= 0)
	sm = (double) cellsmooth[m2];
      else
	sm = 0.01;

      points[m1 * inp_cols + m2].x = x - params->x_orig;
      points[m1 * inp_cols + m2].y = y - params->y_orig;
      if (!G_is_f_null_value (cellinp + m2))
      {
	points[m1 * inp_cols + m2].z = cellinp[m2] * params->zmult - zmin;
      }
      else
      {
	G_set_f_null_value (&(points[m1 * inp_cols + m2].z), 1);
      }
      points[m1 * inp_cols + m2].smooth = sm;
    }
  }
  return 1;
}


int write_zeros (params, data, offset1)

  struct interp_params *params;
  struct quaddata *data;	/* given segment */
  int offset1;			/* offset for temp file writing */

{

  /*
   * C C       INTERPOLATION BY FUNCTIONAL METHOD : TPS + complete regul.
   * c
   */
  double x_or = data->x_orig;
  double y_or = data->y_orig;
  int n_rows = data->n_rows;
  int n_cols = data->n_cols;
  int cond1, cond2;
  double dx, dy, dxx, dyy, dxy;
  int n1, k1, k2, k, i1, l, m, i;
  int ngstc, nszc, ngstr, nszr;
  double zz;
  int offset, offset2;
  double ns_res, ew_res;

  ns_res = (((struct quaddata *) (data))->ymax -
	    ((struct quaddata *) (data))->y_orig) / data->n_rows;
  ew_res = (((struct quaddata *) (data))->xmax -
	    ((struct quaddata *) (data))->x_orig) / data->n_cols;

  cond2 = ((params->adxx != NULL) || (params->adyy != NULL) || (params->adxy != NULL));
  cond1 = ((params->adx != NULL) || (params->ady != NULL) || cond2);

  ngstc = (int) (x_or / ew_res + 0.5) + 1;
  nszc = ngstc + n_cols - 1;
  ngstr = (int) (y_or / ns_res + 0.5) + 1;
  nszr = ngstr + n_rows - 1;

  for (k = ngstr; k <= nszr; k++)
  {
    offset = offset1 * (k - 1);	/* rows offset */
    for (l = ngstc; l <= nszc; l++)
    {
      /*
       * params->az[l] = 0.;
       */
      G_set_f_null_value (params->az + l, 1);
      if (cond1)
      {
	/*
	 * params->adx[l] = (FCELL)0.; params->ady[l] = (FCELL)0.;
	 */
	G_set_f_null_value (params->adx + l, 1);
	G_set_f_null_value (params->ady + l, 1);
	if (cond2)
	{
	  G_set_f_null_value (params->adxx + l, 1);
	  G_set_f_null_value (params->adyy + l, 1);
	  G_set_f_null_value (params->adxy + l, 1);
	  /*
	   * params->adxx[l] = (FCELL)0.; params->adyy[l] = (FCELL)0.;
	   * params->adxy[l] = (FCELL)0.;
	   */
	}
      }
    }
    offset2 = (offset + ngstc - 1) * sizeof (FCELL);
    if (params->wr_temp (params, ngstc, nszc, offset2) < 0)
      return -1;
  }
  return 1;
}

