int read_range(void)
{
    struct FPRange drange;
    struct Range range;
    CELL tmp_min, tmp_max;
    DCELL tmp_dmin, tmp_dmax;
    char buff[1024];
    int i;

    /* read the fpranges and ranges of all input maps */
    for (i = 0; i < noi; i++) {
        if (Rast_read_fp_range(name[i], G_mapset(), &drange) <= 0) {
            sprintf(buff, "Can't read f_range for map %s", name[i]);
            G_fatal_error("%s", buff);
        }
        Rast_get_fp_range_min_max(&drange, &tmp_dmin, &tmp_dmax);

        if (Rast_read_range(name[i], G_mapset(), &range) <= 0) {
            sprintf(buff, "Can't read range for map %s", name[i]);
            G_fatal_error("%s", buff);
        }
        Rast_get_range_min_max(&range, &tmp_min, &tmp_max);
        if (!i || tmp_max > old_max || Rast_is_c_null_value(&old_max))
            old_max = tmp_max;
        if (!i || tmp_min < old_min || Rast_is_c_null_value(&old_min))
            old_min = tmp_min;
        if (!i || tmp_dmax > old_dmax || Rast_is_d_null_value(&old_dmax))
            old_dmax = tmp_dmax;
        if (!i || tmp_dmin < old_dmin || Rast_is_d_null_value(&old_dmin))
            old_dmin = tmp_dmin;
    }				/* for loop */

    return 0;
}

int
Rast3d_read_colors(const char *name, const char *mapset, struct Colors *colors)
 /* adapted from Rast_read_colors */
{
    const char *err;
    struct FPRange drange;
    DCELL dmin, dmax;

    Rast_init_colors(colors);

    Rast_mark_colors_as_fp(colors);

    switch (read_colors(name, mapset, colors)) {
    case -2:
	if (Rast3d_read_range(name, mapset, &drange) >= 0) {
	    Rast_get_fp_range_min_max(&drange, &dmin, &dmax);
	    if (!Rast_is_d_null_value(&dmin) && !Rast_is_d_null_value(&dmax))
		Rast_make_rainbow_fp_colors(colors, dmin, dmax);
	    return 0;
	}
	err = "missing";
	break;
    case -1:
	err = "invalid";
	break;
    default:
	return 1;
    }

    G_warning("color support for [%s] in mapset [%s] %s", name, mapset, err);
    return -1;
}

static int get_cell(DCELL *result, int fd, double x, double y)
{
    static DCELL *row1, *row2;
    static int cur_row = -1;
    static int row, col;
    DCELL *tmp;

    if (!row1) {
	row1 = Rast_allocate_d_buf();
	row2 = Rast_allocate_d_buf();
    }

    col = (int)floor(x - 0.5);
    row = (int)floor(y - 0.5);
    x -= col + 0.5;
    y -= row + 0.5;

    if (row < 0 || row + 1 >= Rast_window_rows() ||
	col < 0 || col + 1 >= Rast_window_cols()) {
	Rast_set_d_null_value(result, 1);
	return 0;
    }

    if (cur_row != row) {
	if (cur_row == row + 1) {
	    tmp = row1; row1 = row2; row2 = tmp;
	    Rast_get_d_row(fd, row1, row);
	}
	else if (cur_row == row - 1) {
	    tmp = row1; row1 = row2; row2 = tmp;
	    Rast_get_d_row(fd, row2, row + 1);
	}
	else {
	    Rast_get_d_row(fd, row1, row);
	    Rast_get_d_row(fd, row2, row + 1);
	}
	cur_row = row;
    }

    if (Rast_is_d_null_value(&row1[col]) || Rast_is_d_null_value(&row1[col+1]) ||
	Rast_is_d_null_value(&row2[col]) || Rast_is_d_null_value(&row2[col+1])) {
	Rast_set_d_null_value(result, 1);
	return 0;
    }

    *result = Rast_interp_bilinear(x, y,
				row1[col], row1[col+1],
				row2[col], row2[col+1]);

    return 1;
}

/************************************************************************
getpoint-- finds crossing point using linear interpolation, 
	  converts from row-column to x-y space, and adds point to current  line.
************************************************************************/
static void getpoint(struct cell *curr, double level,
		     struct Cell_head Cell, struct line_pnts *Points)
{
    double x, y;
    double ratio;
    int p1, p2;

    p1 = curr->edge;
    p2 = (curr->edge + 1) % 4;
    if (Rast_raster_cmp(&curr->z[p1], &curr->z[p2], DCELL_TYPE) == 0)
	ratio = 1;
    else if (Rast_is_d_null_value(&curr->z[p1]))
	ratio = 1 / 2;
    else if (Rast_is_d_null_value(&curr->z[p2]))
	ratio = 1 / 2;
    else
	ratio = (level - curr->z[p1]) / (curr->z[p2] - curr->z[p1]);

    switch (curr->edge) {

    case 0:
	y = curr->r;
	x = curr->c + ratio;
	break;
    case 1:
	y = curr->r + ratio;
	x = curr->c + 1;
	break;
    case 2:
	y = curr->r + 1;
	x = curr->c + 1 - ratio;
	break;
    case 3:
	y = curr->r + 1 - ratio;
	x = curr->c;
	break;
    default:
	G_fatal_error(_("Edge number out of range"));
    }
    /* convert r/c values to x/y values */

    y = Cell.north - (y + .5) * Cell.ns_res;
    x = Cell.west + (x + .5) * Cell.ew_res;

    Vect_append_point(Points, x, y, level);

}

/* calculate the running area of null data cells */
void add_null_area(DCELL * rast, struct Cell_head *region, double *area)
{
    int col;

    for (col = 0; col < region->cols; col++) {
	if (Rast_is_d_null_value(&(rast[col]))) {
	    *area += region->ew_res * region->ns_res;
	}
    }
}

/*!
 * \brief Get minumum and maximum value from fp range
 *
 * Extract the min/max from the range structure <i>range</i>.  If the
 * range structure has no defined min/max (first!=0) there will not be
 * a valid range. In this case the min and max returned must be the
 * NULL-value.
 *
 * \param range pointer to FPRange which holds fp range info
 * \param[out] min minimum value
 * \param[out] max maximum value
 */
void Rast_get_fp_range_min_max(const struct FPRange *range,
			       DCELL * min, DCELL * max)
{
    if (range->first_time) {
	Rast_set_d_null_value(min, 1);
	Rast_set_d_null_value(max, 1);
    }
    else {
	if (Rast_is_d_null_value(&(range->min)))
	    Rast_set_d_null_value(min, 1);
	else
	    *min = range->min;

	if (Rast_is_d_null_value(&(range->max)))
	    Rast_set_d_null_value(max, 1);
	else
	    *max = range->max;
    }
}

/**************************************************************
 * apply_filter: apply the filter to a single neighborhood
 *
 *  filter:    filter to be applied
 *  input:     input buffers
 **************************************************************/
DCELL apply_filter(FILTER * filter, DCELL ** input)
{
    int size = filter->size;
    double **matrix = filter->matrix;
    double divisor = filter->divisor;
    int r, c;
    DCELL v;

    v = 0;

    if (divisor == 0) {
	int have_result = 0;

	for (r = 0; r < size; r++)
	    for (c = 0; c < size; c++) {
		if (Rast_is_d_null_value(&input[r][c]))
		    continue;
		v += input[r][c] * matrix[r][c];
		divisor += filter->dmatrix[r][c];
		have_result = 1;
	    }

	if (have_result)
	    v /= divisor;
	else
	    Rast_set_d_null_value(&v, 1);
    }
    else {
	for (r = 0; r < size; r++)
	    for (c = 0; c < size; c++) {
		if (Rast_is_d_null_value(&input[r][c])) {
		    Rast_set_d_null_value(&v, 1);
		    return v;
		}
		v += input[r][c] * matrix[r][c];
	    }

	v /= divisor;
    }

    return v;
}

static double loglike(DCELL * x, struct SubSig *SubSig, int nbands)
{
    int b1, b2;
    double diff1, diff2;
    double sum;

    sum = 0;
    for (b1 = 0; b1 < nbands; b1++)
	for (b2 = 0; b2 < nbands; b2++) {
	    if (Rast_is_d_null_value(&x[b1])
		|| Rast_is_d_null_value(&x[b2]))
		continue;
	    diff1 = x[b1] - SubSig->means[b1];
	    diff2 = x[b2] - SubSig->means[b2];
	    sum += diff1 * diff2 * SubSig->Rinv[b1][b2];
	}

    sum = -0.5 * sum + SubSig->cnst;
    return (sum);
}

void w_kurt(DCELL * result, DCELL(*values)[2], int n, const void *closure)
{
    DCELL sum, ave, sumsq, sumqt, var;
    int count;
    int i;

    sum = 0.0;
    count = 0;

    for (i = 0; i < n; i++) {
	if (Rast_is_d_null_value(&values[i][0]))
	    continue;

	sum += values[i][0] * values[i][1];
	count += values[i][1];
    }

    if (count == 0) {
	Rast_set_d_null_value(result, 1);
	return;
    }

    ave = sum / count;

    sumsq = 0;

    for (i = 0; i < n; i++) {
	DCELL d;

	if (Rast_is_d_null_value(&values[i][0]))
	    continue;

	d = values[i][0] - ave;
	sumsq += d * d * values[i][1];
	sumqt += d * d * d * values[i][1];
    }

    var = sumsq / count;

    *result = sumqt / (count * var * var) - 3;
}

void c_skew(DCELL * result, DCELL * values, int n, const void *closure)
{
    DCELL sum, ave, sumsq, sumcb, sdev;
    int count;
    int i;

    sum = 0.0;
    count = 0;

    for (i = 0; i < n; i++) {
	if (Rast_is_d_null_value(&values[i]))
	    continue;

	sum += values[i];
	count++;
    }

    if (count == 0) {
	Rast_set_d_null_value(result, 1);
	return;
    }

    ave = sum / count;

    sumsq = 0;

    for (i = 0; i < n; i++) {
	DCELL d;

	if (Rast_is_d_null_value(&values[i]))
	    continue;

	d = values[i] - ave;
	sumsq += d * d;
	sumcb += d * d * d;
    }

    sdev = sqrt(sumsq / count);

    *result = sumcb / (count * sdev * sdev * sdev);
}

int report_range(void)
{
    char buff[300], buff2[300];

    if (Rast_is_d_null_value(&old_dmin) || Rast_is_d_null_value(&old_dmax))
        G_message(_("Old data range is empty"));
    else {
        sprintf(buff, "%.15g", old_dmin);
        sprintf(buff2, "%.15g", old_dmax);
        G_trim_decimal(buff);
        G_trim_decimal(buff2);
        G_message(_("Old data range is %s to %s"), buff, buff2);
    }
    if (Rast_is_c_null_value(&old_min) || Rast_is_c_null_value(&old_max))
        G_message(_("Old integer data range is empty"));
    else
        G_message(_("Old integer data range is %d to %d"),
                  (int)old_min, (int)old_max);

    return 0;
}

void p_cubic_f(struct cache *ibuffer,	/* input buffer                  */
		void *obufptr,	/* ptr in output buffer          */
		int cell_type,	/* raster map type of obufptr    */
		double *row_idx,	/* row index                     */
		double *col_idx,	/* column index          */
	    struct Cell_head *cellhd	/* cell header of input layer    */
    )
{
    /* start nearest neighbor to do some basic tests */
    int row, col;		/* row/col of nearest neighbor   */
    DCELL *cellp, cell;

    /* cut indices to integer */
    row = (int)floor(*row_idx);
    col = (int)floor(*col_idx);

    /* check for out of bounds - if out of bounds set NULL value     */
    if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
        Rast_set_null_value(obufptr, 1, cell_type);
        return;
    }

    cellp = CPTR(ibuffer, row, col);
    /* if nearest is null, all the other interps will be null */
    if (Rast_is_d_null_value(cellp)) {
        Rast_set_null_value(obufptr, 1, cell_type);
        return;
    }
    cell = *cellp;
    
    p_cubic(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
    /* fallback to bilinear if cubic is null */
    if (Rast_is_d_null_value(obufptr)) {
        p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
        /* fallback to nearest if bilinear is null */
	if (Rast_is_d_null_value(obufptr))
	    Rast_set_d_value(obufptr, cell, cell_type);
    }
}

void c_var(DCELL * result, DCELL * values, int n, const void *closure)
{
    DCELL sum, ave, sumsq;
    int count;
    int i;

    sum = 0.0;
    count = 0;

    for (i = 0; i < n; i++) {
	if (Rast_is_d_null_value(&values[i]))
	    continue;

	sum += values[i];
	count++;
    }

    if (count == 0) {
	Rast_set_d_null_value(result, 1);
	return;
    }

    ave = sum / count;

    sumsq = 0;

    for (i = 0; i < n; i++) {
	DCELL d;

	if (Rast_is_d_null_value(&values[i]))
	    continue;

	d = values[i] - ave;
	sumsq += d * d;
    }

    *result = sumsq / count;
}

void c_maxx(DCELL * result, DCELL * values, int n, const void *closure)
{
    DCELL max, maxx;
    int i;

    Rast_set_d_null_value(&max, 1);
    Rast_set_d_null_value(&maxx, 1);

    for (i = 0; i < n; i++) {
	if (Rast_is_d_null_value(&values[i]))
	    continue;

	if (Rast_is_d_null_value(&max) || max < values[i]) {
	    max = values[i];
	    maxx = i;
	}
    }

    if (Rast_is_d_null_value(&maxx))
	Rast_set_d_null_value(result, 1);
    else
	*result = maxx;
}

static void convert_row(
    unsigned char *out_buf, const DCELL *raster, int ncols,
    int is_fp, int bytes, int swap_flag, double null_val)
{
    unsigned char *ptr = out_buf;
    int i;

    for (i = 0; i < ncols; i++) {
	DCELL x = Rast_is_d_null_value(&raster[i])
	    ? null_val
	    : raster[i];
	convert_cell(ptr, x, is_fp, bytes, swap_flag);
	ptr += bytes;
    }
}

static void convert_and_write_di(int fd, const void *vbuf)
{
    const DCELL *buf = vbuf;
    struct fileinfo *fcb = &R__.fileinfo[fd];
    CELL *p = (CELL *) fcb->data;
    int i;

    for (i = 0; i < fcb->cellhd.cols; i++)
	if (Rast_is_d_null_value(&buf[i]))
	    Rast_set_c_null_value(&p[i], 1);
	else
	    p[i] = (CELL) buf[i];

    Rast_put_c_row(fd, p);
}

/*!
 * \brief Update range structure (floating-point)
 *
 * Compares the <i>cat</i> value with the minimum and maximum values
 * in the <i>range</i> structure, modifying the range if <i>cat</i>
 * extends the range.
 *
 * NULL-values must be detected and ignored.
 *
 * \param val raster value
 * \param range pointer to Range structure which holds range info
 */
void Rast_update_fp_range(DCELL val, struct FPRange *range)
{
    if (!Rast_is_d_null_value(&val)) {
	if (range->first_time) {
	    range->first_time = 0;
	    range->min = val;
	    range->max = val;
	    return;
	}
	if (val < range->min)
	    range->min = val;
	if (val > range->max)
	    range->max = val;
    }
}

int report_range(void)
{
    struct FPRange drange;
    struct Range range;
    char buff[1024], buff2[300];
    RASTER_MAP_TYPE inp_type;

    inp_type = Rast_map_type(name, "");
    if (inp_type != CELL_TYPE) {
	if (Rast_read_fp_range(name, "", &drange) <= 0)
	    G_fatal_error(_("Unable to read f_range for map %s"), name);

	Rast_get_fp_range_min_max(&drange, &old_dmin, &old_dmax);
	if (Rast_is_d_null_value(&old_dmin) || Rast_is_d_null_value(&old_dmax))
	    G_message(_("Data range is empty"));
	else {
	    sprintf(buff, "%.10f", old_dmin);
	    sprintf(buff2, "%.10f", old_dmax);
	    G_trim_decimal(buff);
	    G_trim_decimal(buff2);
	    G_message(_("Data range of %s is %s to %s (entire map)"), name,
		      buff, buff2);
	}
    }
    if (Rast_read_range(name, "", &range) <= 0)
	G_fatal_error(_("Unable to read range for map <%s>"), name);

    Rast_get_range_min_max(&range, &old_min, &old_max);
    if (Rast_is_c_null_value(&old_min) || Rast_is_c_null_value(&old_max))
	G_message(_("Integer data range of %s is empty"), name);
    else
	G_message(_("Integer data range of %s is %d to %d"),
		  name, (int)old_min, (int)old_max);

    return 0;
}

static int gather_values(RASTER3D_Map * map, DCELL * buff, int x,
			 int y, int z)
{
    int i, j, k, l;
    DCELL value;

    int start_z = z - z_dist;
    int start_y = y - y_dist;
    int start_x = x - x_dist;
    int end_z = start_z + z_size;
    int end_y = start_y + y_size;
    int end_x = start_x + x_size;

    if (start_z < 0)
	start_z = 0;

    if (start_y < 0)
	start_y = 0;

    if (start_x < 0)
	start_x = 0;

    if (end_z > nz)
	end_z = nz;

    if (end_y > ny)
	end_y = ny;

    if (end_x > nx)
	end_x = nx;

    l = 0;

    for (i = start_z; i < end_z; i++) {
	for (j = start_y; j < end_y; j++) {
	    for (k = start_x; k < end_x; k++) {
		value = (DCELL) Rast3d_get_double(map, k, j, i);

		if (Rast_is_d_null_value(&value))
		    continue;

		buff[l] = value;
		l++;
	    }
	}
    }
    return l;
}

int is_null_value(struct RASTER_MAP_PTR buf, int col)
{
    switch (buf.type) {
    case CELL_TYPE:
	return Rast_is_c_null_value(&buf.data.c[col]);
	break;
    case FCELL_TYPE:
	return Rast_is_f_null_value(&buf.data.f[col]);
	break;
    case DCELL_TYPE:
	return Rast_is_d_null_value(&buf.data.d[col]);
	break;
    }

    return -1;
}