//.........这里部分代码省略.........
    if (method == METHOD_N)
	rtype = CELL_TYPE;


    G_get_window(&region);
    rows = (int)(region.rows * (percent / 100.0));
    cols = region.cols;

    G_debug(2, "region.n=%f  region.s=%f  region.ns_res=%f", region.north,
	    region.south, region.ns_res);
    G_debug(2, "region.rows=%d  [box_rows=%d]  region.cols=%d", region.rows,
	    rows, region.cols);

    npasses = (int)ceil(1.0 * region.rows / rows);

    if (!scan_flag->answer) {
	/* allocate memory (test for enough before we start) */
	if (bin_n)
	    n_array = G_calloc(rows * (cols + 1), G_raster_size(CELL_TYPE));
	if (bin_min)
	    min_array = G_calloc(rows * (cols + 1), G_raster_size(rtype));
	if (bin_max)
	    max_array = G_calloc(rows * (cols + 1), G_raster_size(rtype));
	if (bin_sum)
	    sum_array = G_calloc(rows * (cols + 1), G_raster_size(rtype));
	if (bin_sumsq)
	    sumsq_array = G_calloc(rows * (cols + 1), G_raster_size(rtype));
	if (bin_index)
	    index_array =
		G_calloc(rows * (cols + 1), G_raster_size(CELL_TYPE));

	/* and then free it again */
	if (bin_n)
	    G_free(n_array);
	if (bin_min)
	    G_free(min_array);
	if (bin_max)
	    G_free(max_array);
	if (bin_sum)
	    G_free(sum_array);
	if (bin_sumsq)
	    G_free(sumsq_array);
	if (bin_index)
	    G_free(index_array);

	/** end memory test **/
    }


    /* open input file */
    if (strcmp("-", infile) == 0) {
	from_stdin = TRUE;
	in_fp = stdin;
	infile = G_store("stdin");	/* filename for history metadata */
    }
    else {
	if ((in_fp = fopen(infile, "r")) == NULL)
	    G_fatal_error(_("Unable to open input file <%s>"), infile);
    }

    can_seek = fseek(in_fp, 0, SEEK_SET) == 0;

    /* can't rewind() non-files */
    if (!can_seek && npasses != 1) {
	G_warning(_("If input is not from a file it is only possible to perform a single pass."));
	npasses = 1;

//.........这里部分代码省略.........
		    if (columns[j] == NULL) {
			trans_params[j] = trans_params_def[j];
			continue;
		    }
		    ctype = db_column_Ctype(driver, fi->table, columns[j]);
		    switch (ctype) {
		    case DB_C_TYPE_INT:
		    case DB_C_TYPE_DOUBLE:
		    case DB_C_TYPE_STRING:
			break;
		    case -1:
			G_fatal_error(_("Column <%s> not found in table <%s>"),
				      columns[j], fi->table);
		    default:
			G_fatal_error(_("Unsupported column type of <%s>"),
				      columns[j]);
		    }
		    if (db_select_value
			(driver, fi->table, fi->key, cat, columns[j], &val) != 1
			|| db_test_value_isnull(&val)) {
			trans_params[j] = trans_params_def[j];

			G_warning(_("Unable to select value for category %d from table <%s>, column <%s>. "
				   "For category %d using default transformation parameter %.3f."),
				  cat, fi->table, columns[j], cat,
				  trans_params[j]);
		    }
		    else {
			trans_params[j] = db_get_value_as_double(&val, ctype);
		    }
		}
	    }
	    else {
		G_warning(_("No category number defined. Using default transformation parameters."));

		for (j = 0; j <= IDX_ZROT; j++) {
		    trans_params[j] = trans_params_def[j];
		}
	    }
	    ang = PI * trans_params[IDX_ZROT] / 180;
	}

	/* transform points */
	for (i = 0; i < Points->n_points; i++) {
	    G_debug(3, "idx=%d, cat=%d, xshift=%g, yshift=%g, zshift=%g, "
		    "xscale=%g, yscale=%g, zscale=%g, zrot=%g",
		    i, cat, trans_params[IDX_XSHIFT],
		    trans_params[IDX_YSHIFT], trans_params[IDX_ZSHIFT],
		    trans_params[IDX_XSCALE], trans_params[IDX_YSCALE],
		    trans_params[IDX_ZSCALE], trans_params[IDX_ZROT]);

	    /* transform point */
	    x = trans_params[IDX_XSHIFT] +
		trans_params[IDX_XSCALE] * Points->x[i] * cos(ang)
		- trans_params[IDX_YSCALE] * Points->y[i] * sin(ang);
	    y = trans_params[IDX_YSHIFT] +
		trans_params[IDX_XSCALE] * Points->x[i] * sin(ang)
		+ trans_params[IDX_YSCALE] * Points->y[i] * cos(ang);
	    Points->x[i] = x;
	    Points->y[i] = y;

	    /* ztozero shifts oldmap z to zero, zshift shifts rescaled object
	     * to target elevation: */
	    Points->z[i] =
		((Points->z[i] + ztozero) * trans_params[IDX_ZSCALE]) +
		trans_params[IDX_ZSHIFT];

            if (swap_after) {
                if (swap_xy) {
                    tmp = Points->x[i];
                    Points->x[i] = Points->y[i];
                    Points->y[i] = tmp;
                }
                if (swap_xz) {
                    tmp = Points->z[i];
                    Points->z[i] = Points->x[i];
                    Points->x[i] = tmp;
                }
                if (swap_yz) {
                    tmp = Points->z[i];
                    Points->z[i] = Points->y[i];
                    Points->y[i] = tmp;
                }
            }
	}

	Vect_write_line(New, type, Points, Cats);
	line++;
    }

    if (verbose && format != G_INFO_FORMAT_PLAIN)
	fprintf(stderr, "\r");
    
    if (field > 0) {
	db_close_database_shutdown_driver(driver);
	G_free((void *)trans_params);
    }
    
    return ret;
}

//.........这里部分代码省略.........

	Nviz_get_exag_height(&vp_height, &min, &max);
	G_verbose_message(_("Viewpoint height not given, using calculated "
			    "value %.0f"), vp_height);
    }
    Nviz_set_viewpoint_height(vp_height);

    Nviz_set_viewpoint_position(atof(params->pos->answers[0]),
				atof(params->pos->answers[1]));
    Nviz_set_viewpoint_twist(atoi(params->twist->answer));
    Nviz_set_viewpoint_persp(atoi(params->persp->answer));

    if (params->focus->answer) {
	Nviz_set_focus(&data, atof(params->focus->answers[0]),
		       atof(params->focus->answers[1]),
		       atof(params->focus->answers[2]));
    }

    /* set lights */
    Nviz_set_light_position(&data, 1,
			    atof(params->light_pos->answers[0]),
			    atof(params->light_pos->answers[1]),
			    atof(params->light_pos->answers[2]), 0.0);
    Nviz_set_light_bright(&data, 1,
			  atoi(params->light_bright->answer) / 100.0);
    if (G_str_to_color(params->light_color->answer, &red, &grn, &blu) != 1) {
	red = grn = blu = 255;
    }
    Nviz_set_light_color(&data, 1, red, grn, blu);
    Nviz_set_light_ambient(&data, 1,
			   atof(params->light_ambient->answer) / 100.0);

    /* define fringes */
    if (params->fringe->answer) {
	int nw, ne, sw, se;

	i = 0;
	nw = ne = sw = se = 0;
	while (params->fringe->answers[i]) {
	    const char *edge = params->fringe->answers[i++];

	    if (strcmp(edge, "nw") == 0)
		nw = 1;
	    else if (strcmp(edge, "ne") == 0)
		ne = 1;
	    else if (strcmp(edge, "sw") == 0)
		sw = 1;
	    else if (strcmp(edge, "se") == 0)
		se = 1;
	}
	Nviz_new_fringe(&data, -1,
			Nviz_color_from_str(params->fringe_color->answer),
			atof(params->fringe_elev->answer), nw, ne, sw, se);
    }

    /* draw north arrow */
    if (params->north_arrow->answer) {

	if (!params->north_arrow_size->answer)
	    size = Nviz_get_longdim(&data) / 8.;
	else
	    size = atof(params->north_arrow_size->answer);

	Nviz_set_arrow(&data, atoi(params->north_arrow->answers[0]),
		       atoi(params->north_arrow->answers[1]),
		       size,
		       Nviz_color_from_str(params->north_arrow_color->
					   answer));
	Nviz_draw_arrow(&data);
    }

    GS_clear(data.bgcolor);

    /* cutting planes */
    if (params->cplane->answer)
	draw_cplane(params, &data);

    /* draw */
    Nviz_draw_all(&data);

    /* write to image */
    ret = 0;
    if (strcmp(params->format->answer, "ppm") == 0)
	ret = write_img(output_name, FORMAT_PPM);
    if (strcmp(params->format->answer, "tif") == 0)
	ret = write_img(output_name, FORMAT_TIF);

    if (!ret)
	G_fatal_error(_("Unsupported output format"));

    G_done_msg(_("File <%s> created."), output_name);

    Nviz_destroy_data(&data);
    Nviz_destroy_render_window(offscreen);

    G_free((void *)output_name);
    G_free((void *)params);

    exit(EXIT_SUCCESS);
}

/*!
   \brief Free allocated memory

   \param fv pointer to geovect struct
 */
void gv_free_vectmem(geovect * fv)
{
    geoline *gln, *tmpln;
    
    G_free((void *)fv->filename);
    fv->filename = NULL;
    if (fv->style)
	G_free(fv->style);
    if (fv->hstyle)
	G_free(fv->hstyle);

    if (fv->lines) {
	for (gln = fv->lines; gln;) {
	    if (gln->dims == 2) {
		sub_Vectmem(gln->npts * sizeof(Point2));
		G_free(gln->p2);
	    }

	    if (gln->dims == 3) {
		G_free(gln->p3);
	    }

	    G_free(gln->cats);
	    
	    tmpln = gln;
	    gln = gln->next;
	    sub_Vectmem(sizeof(geoline));
	    G_free(tmpln);
	}

	fv->n_lines = 0;
	fv->lines = NULL;
    }

    if (fv->tstyle) {
	G_free(fv->tstyle->color_column);
	G_free(fv->tstyle->symbol_column);
	G_free(fv->tstyle->size_column);
	G_free(fv->tstyle->width_column);
    }

    return;
}

int douglas_peucker(struct line_pnts *Points, double thresh, int with_z)
{
    int *stack = G_malloc(sizeof(int) * Points->n_points * 2);

    if (!stack) {
	G_fatal_error(_("Out of memory"));
	return Points->n_points;
    }

    int *index = G_malloc(sizeof(int) * Points->n_points);	/* Indices of points in output line */

    if (!index) {
	G_fatal_error(_("Out of memory"));
	G_free(stack);
	return Points->n_points;
    }

    int top = 2;		/* first free slot in the stack */
    int icount = 1;		/* number of indices stored */
    int i;

    thresh *= thresh;

    index[0] = 0;		/* first point is always in output line */

    /* stack contains pairs of elements: (beginning, end) */

    stack[0] = 0;
    stack[1] = Points->n_points - 1;

    while (top > 0) {		/*there are still segments to consider */
	/*Pop indices of the segment from the stack */
	int last = stack[--top];
	int first = stack[--top];

	double x1 = Points->x[first];
	double y1 = Points->y[first];
	double z1 = Points->z[first];

	double x2 = Points->x[last];
	double y2 = Points->y[last];
	double z2 = Points->z[last];

	int maxindex = -1;
	double maxdist = -1;
	int i;

	for (i = first + 1; i <= last - 1; i++) {	/* Find the furthermost point between first, last */
	    double px, py, pz, pdist;
	    int status;
	    double dist =
		dig_distance2_point_to_line(Points->x[i], Points->y[i],
					    Points->z[i],
					    x1, y1, z1, x2, y2, z2, with_z,
					    &px, &py, &pz, &pdist, &status);

	    if (maxindex == -1 || dist > maxdist) {	/* update the furthermost point so far seen */
		maxindex = i;
		maxdist = dist;
	    }
	}


	if (maxindex == -1 || maxdist <= thresh) {	/* no points between or all point are inside the threshold */
	    index[icount++] = last;
	}
	else {
	    /* break line into two parts, the order of pushing is crucial! It gurantees, that we are going to the left */
	    stack[top++] = maxindex;
	    stack[top++] = last;
	    stack[top++] = first;
	    stack[top++] = maxindex;
	}


    }


    Points->n_points = icount;

    /* finally, select only points marked in the algorithm */
    for (i = 0; i < icount; i++) {
	Points->x[i] = Points->x[index[i]];
	Points->y[i] = Points->y[index[i]];
	Points->z[i] = Points->z[index[i]];
    }

    G_free(stack);
    G_free(index);
    return (Points->n_points);
}

//.........这里部分代码省略.........

	G_debug(1, "DCELL nodata val: %f", dnullval);

	for (row = 0; row < rows; row++) {

	    Rast_get_row(fd, bufer, row, maptype);
	    for (col = 0; col < cols; col++) {
		if (Rast_is_d_null_value(&((DCELL *) bufer)[col])) {
		    ((DCELL *) bufer)[col] = dnullval;
		    n_nulls++;
		}
		else {
		    if (((DCELL *) bufer)[col] == dnullval) {
			nodatavalmatch = 1;
		    }
		    if (dfCellMin > ((DCELL *) bufer)[col])
			dfCellMin = ((DCELL *) bufer)[col];
		    if (dfCellMax < ((DCELL *) bufer)[col])
			dfCellMax = ((DCELL *) bufer)[col];
		}
	    }
	    G_percent(row + 1, rows, 2);
	}
    }
    else {

	CELL inullval = (CELL) nodataval;

	G_debug(1, "CELL nodata val: %d", inullval);

	for (row = 0; row < rows; row++) {

	    Rast_get_row(fd, bufer, row, maptype);
	    for (col = 0; col < cols; col++) {
		if (Rast_is_c_null_value(&((CELL *) bufer)[col])) {
		    ((CELL *) bufer)[col] = inullval;
		    n_nulls++;
		}
		else {
		    if (((CELL *) bufer)[col] == inullval) {
			nodatavalmatch = 1;
		    }
		    if (dfCellMin > ((CELL *) bufer)[col])
			dfCellMin = ((CELL *) bufer)[col];
		    if (dfCellMax < ((CELL *) bufer)[col])
			dfCellMax = ((CELL *) bufer)[col];
		}
	    }
	    G_percent(row + 1, rows, 2);
	}
    }
    G_debug(1, "min %g max %g", dfCellMin, dfCellMax);

    /* can the GDAL datatype hold the data range to be exported ? */
    /* f-flag does not override */
    if (exact_range_check(dfCellMin, dfCellMax, export_datatype, name)) {
	G_warning("Raster export results in data loss.");
	ret = -2;
    }
    G_message(_("Using GDAL data type <%s>"), GDALGetDataTypeName(export_datatype));

    /* a default nodata value was used and NULL cells were present */
    if (n_nulls && default_nodataval) {
	if (maptype == CELL_TYPE)
	    G_important_message(_("Input raster map contains cells with NULL-value (no-data). "
				 "The value %d will be used to represent no-data values in the input map. "
				 "You can specify a nodata value with the %s option."),
				(int)nodataval, nodatakey);
	else
	    G_important_message(_("Input raster map contains cells with NULL-value (no-data). "
				 "The value %g will be used to represent no-data values in the input map. "
				 "You can specify a nodata value with the %s option."),
				nodataval, nodatakey);
    }

    /* the nodata value was present in the exported data */
    if (nodatavalmatch && n_nulls) {
	/* default nodataval didn't work */
	if (default_nodataval) {
	    G_warning(_("The default nodata value is present in raster"
			"band <%s> and would lead to data loss. Please specify a "
			"custom nodata value with the %s parameter."),
		      name, nodatakey);
	}
	/* user-specified nodataval didn't work */
	else {
	    G_warning(_("The user given nodata value %g is present in raster"
			"band <%s> and would lead to data loss. Please specify a "
			"different nodata value with the %s parameter."),
		      nodataval, name, nodatakey);
	}
	ret = -1;
    }

    Rast_close(fd);

    G_free(bufer);

    return ret;
}

/*!
 *  \brief Extract a cell value from raster map (bilinear interpolation).
 *
 *  Extract a cell value from raster map at given northing and easting
 *  with a sampled 3x3 window using a bilinear interpolation.
 *
 *  \param fd file descriptor
 *  \param window region settings
 *  \param cats categories
 *  \param north northing position
 *  \param east easting position
 *  \param usedesc flag to scan category label
 *
 *  \return cell value at given position
 */
DCELL Rast_get_sample_bilinear(int fd,
			       const struct Cell_head * window,
			       struct Categories * cats,
			       double north, double east, int usedesc)
{
    int row, col;
    double grid[2][2];
    DCELL *arow = Rast_allocate_d_buf();
    DCELL *brow = Rast_allocate_d_buf();
    double frow, fcol, trow, tcol;
    DCELL result;

    frow = Rast_northing_to_row(north, window);
    fcol = Rast_easting_to_col(east, window);

    /* convert northing and easting to row and col, resp */
    row = (int)floor(frow - 0.5);
    col = (int)floor(fcol - 0.5);

    trow = frow - row - 0.5;
    tcol = fcol - col - 0.5;

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

    Rast_get_d_row(fd, arow, row);
    Rast_get_d_row(fd, brow, row + 1);

    if (Rast_is_d_null_value(&arow[col]) ||
	Rast_is_d_null_value(&arow[col + 1]) ||
	Rast_is_d_null_value(&brow[col]) ||
	Rast_is_d_null_value(&brow[col + 1])) {
	Rast_set_d_null_value(&result, 1);
	goto done;
    }

    /*-
     * now were ready to do bilinear interpolation over
     * arow[col], arow[col+1],
     * brow[col], brow[col+1]
     */

    if (usedesc) {
	char *buf;

	G_squeeze(buf = Rast_get_c_cat((int *)&(arow[col]), cats));
	grid[0][0] = scancatlabel(buf);
	G_squeeze(buf = Rast_get_c_cat((CELL *) & (arow[col + 1]), cats));
	grid[0][1] = scancatlabel(buf);
	G_squeeze(buf = Rast_get_c_cat((CELL *) & (brow[col]), cats));
	grid[1][0] = scancatlabel(buf);
	G_squeeze(buf = Rast_get_c_cat((CELL *) & (brow[col + 1]), cats));
	grid[1][1] = scancatlabel(buf);
    }
    else {
	grid[0][0] = arow[col];
	grid[0][1] = arow[col + 1];
	grid[1][0] = brow[col];
	grid[1][1] = brow[col + 1];
    }

    result = Rast_interp_bilinear(tcol, trow,
				  grid[0][0], grid[0][1], grid[1][0],
				  grid[1][1]);

  done:
    G_free(arow);
    G_free(brow);

    return result;
}

//.........这里部分代码省略.........
    for (i = 0; i <= max; i++) {
	sum[i] = 0;
	count[i] = 0;
    }

    data_buf = Rast_allocate_c_buf();
    clump_buf = Rast_allocate_c_buf();

    /* get window size */
    rows = window.rows;
    cols = window.cols;

    /* now get the data -- first pass */
    G_message("Complete ...");
    for (row = 0; row < rows; row++) {
	G_percent(row, rows, 2);
	Rast_get_c_row(fd_data, data_buf, row);
	Rast_get_c_row(fd_clump, clump_buf, row);
	for (col = 0; col < cols; col++) {
	    i = clump_buf[col];
	    if (i > max)
		G_fatal_error(
		    "Row=%d Col=%d Cat=%d in clump map [%s]; max=%d.\n"
		    "Cat value > max returned by Rast_get_max_c_cat.",
		    row, col, i, clumpmap, max);
	    if (i < 1)
		continue;	/* ignore zeros and negs */
	    count[i]++;
	    sum[i] += data_buf[col];
	}
    }
    G_percent(row, rows, 2);
    /* free some buffer space */
    G_free(data_buf);
    G_free(clump_buf);

    /* data lists for centroids of clumps */
    e = (int *)G_malloc((max + 1) * sizeof(int));
    n = (int *)G_malloc((max + 1) * sizeof(int));

    i = centroids(fd_clump, e, n, 1, max);

    /* got everything, now do output */
    if (*site_list) {
	char desc[GNAME_MAX * 2 + 40];

	site_info.form = NULL;
	site_info.time = NULL;
	site_info.stime = NULL;
	sprintf(desc, "from %s on map %s using clumps from %s",
		argv[0], datamap, clumpmap);
	site_info.desc = G_store(desc);
	site_info.name = G_store(site_list);
	site_info.labels =
	    G_store("centroid east|centroid north|#cat vol avg t n");
	G_site_put_head(fd_sites, &site_info);
    }
    if (out_mode) {
	fprintf(stdout, "Volume report on data from %s", datamap);
	fprintf(stdout, " using clumps on %s map\n\n", clumpmap);
	fprintf(stdout,
		" Cat    Average   Data   # Cells        Centroid             Total\n");
	fprintf(stdout,
		"Number  in clump  Total  in clump   Easting   Northing       Volume\n\n");
    }
    total_vol = 0.0;

/*!
 * \brief Create window mapping.
 *
 * Creates mapping from cell header into window. The boundaries and 
 * resolution of the two spaces do not have to be the same or aligned in 
 * any way.
 *
 * \param fd file descriptor
 */
void Rast__create_window_mapping(int fd)
{
    struct fileinfo *fcb = &R__.fileinfo[fd];
    COLUMN_MAPPING *col;
    int i;
    int x;
    double C1, C2;
    double west;

    if (fcb->open_mode >= 0 && fcb->open_mode != OPEN_OLD)	/* open for write? */
	return;
    if (fcb->open_mode == OPEN_OLD)	/* already open ? */
	G_free(fcb->col_map);

    col = fcb->col_map = alloc_index(R__.rd_window.cols);

    /*
     * for each column in the window, go to center of the cell,
     * compute nearest column in the data file
     * if column is not in data file, set column to 0
     *
     * for lat/lon move window so that west is bigger than
     * cellhd west.
     */
    west = R__.rd_window.west;
    if (R__.rd_window.proj == PROJECTION_LL) {
	while (west > fcb->cellhd.west + 360.0)
	    west -= 360.0;
	while (west < fcb->cellhd.west)
	    west += 360.0;
    }

    C1 = R__.rd_window.ew_res / fcb->cellhd.ew_res;
    C2 = (west - fcb->cellhd.west +
	  R__.rd_window.ew_res / 2.0) / fcb->cellhd.ew_res;
    for (i = 0; i < R__.rd_window.cols; i++) {
	x = C2;
	if (C2 < x)		/* adjust for rounding of negatives */
	    x--;
	if (x < 0 || x >= fcb->cellhd.cols)	/* not in data file */
	    x = -1;
	*col++ = x + 1;
	C2 += C1;
    }

    /* do wrap around for lat/lon */
    if (R__.rd_window.proj == PROJECTION_LL) {
	col = fcb->col_map;
	C2 = (west - 360.0 - fcb->cellhd.west +
	      R__.rd_window.ew_res / 2.0) / fcb->cellhd.ew_res;
	for (i = 0; i < R__.rd_window.cols; i++) {
	    x = C2;
	    if (C2 < x)		/* adjust for rounding of negatives */
		x--;
	    if (x < 0 || x >= fcb->cellhd.cols)	/* not in data file */
		x = -1;
	    if (*col == 0)	/* only change those not already set */
		*col = x + 1;
	    col++;
	    C2 += C1;
	}
    }

    G_debug(3, "create window mapping (%d columns)", R__.rd_window.cols);
    /*  for (i = 0; i < R__.rd_window.cols; i++)
       fprintf(stderr, "%s%ld", i % 15 ? " " : "\n", (long)fcb->col_map[i]);
       fprintf(stderr, "\n");
     */

    /* compute C1,C2 for row window mapping */
    fcb->C1 = R__.rd_window.ns_res / fcb->cellhd.ns_res;
    fcb->C2 =
	(fcb->cellhd.north - R__.rd_window.north +
	 R__.rd_window.ns_res / 2.0) / fcb->cellhd.ns_res;
}

static int load_seeds(struct globals *globals, int srows, int scols, int nseg)
{
    int row, col;
    SEGMENT seeds_seg;
    CELL *seeds_buf, seeds_val;
    int seeds_fd;
    int spos, sneg, have_seeds;
    struct rc Ri;

    G_debug(1, "load_seeds()");
    
    G_message(_("Loading seeds from raster map <%s>..."), globals->seeds);

    if (Segment_open
	(&seeds_seg, G_tempfile(), globals->nrows, globals->ncols,
	 srows, scols, sizeof(CELL), nseg) != TRUE)
	G_fatal_error("Unable to create bounds temporary files");

    seeds_fd = Rast_open_old(globals->seeds, "");
    seeds_buf = Rast_allocate_c_buf();
    
    have_seeds = 0;

    /* load seeds map to segment structure */
    for (row = 0; row < globals->nrows; row++) {
	Rast_get_c_row(seeds_fd, seeds_buf, row);
	for (col = 0; col < globals->ncols; col++) {
	    if (FLAG_GET(globals->null_flag, row, col)) {
		Rast_set_c_null_value(&seeds_val, 1);
	    }
	    else {
		seeds_val = seeds_buf[col];
		if (!Rast_is_c_null_value(&seeds_val))
		    have_seeds = 1;
	    }
	    if (Segment_put(&seeds_seg, &seeds_val, row, col) != 1)
		G_fatal_error(_("Unable to write to temporary file"));
	}
    }

    if (!have_seeds) {
	G_warning(_("No seeds found in '%s'!"), globals->seeds);
	G_free(seeds_buf);
	Rast_close(seeds_fd);
	Segment_close(&seeds_seg);
	return 0;
    }

    spos = 1;
    sneg = -1;

    /* convert seeds to regions */
    G_debug(1, "convert seeds to regions");
    Rast_set_c_null_value(&seeds_val, 1);
    for (row = 0; row < globals->nrows; row++) {
	Rast_get_c_row(seeds_fd, seeds_buf, row);
	for (col = 0; col < globals->ncols; col++) {
	    if (!(FLAG_GET(globals->null_flag, row, col)) && 
	        !(FLAG_GET(globals->candidate_flag, row, col))) {

		if (Rast_is_c_null_value(&(seeds_buf[col]))) {
		    if (Segment_put(&globals->rid_seg, &sneg, row, col) != 1)
			G_fatal_error(_("Unable to write to temporary file"));
		    sneg--;
		    globals->n_regions--;
		}
		else {
		    Ri.row = row;
		    Ri.col = col;
		    read_seed(globals, &seeds_seg, &Ri, spos);
		    spos++;
		}
	    }
	}
    }

    G_free(seeds_buf);
    Rast_close(seeds_fd);
    Segment_close(&seeds_seg);

    globals->n_regions = spos - 1;
    
    flag_clear_all(globals->candidate_flag);
    
    return 1;
}

//.........这里部分代码省略.........
               globals->ncols, globals->col_min, globals->col_max);
    
    globals->row_max++;
    globals->col_max++;
    globals->ncells = (long)(globals->row_max - globals->row_min) *
			    (globals->col_max - globals->col_min);

    /* bounds/constraints */

    Rast_set_c_null_value(&globals->upper_bound, 1);
    Rast_set_c_null_value(&globals->lower_bound, 1);

    if (globals->bounds_map != NULL) {
	if (Segment_open
	    (&globals->bounds_seg, G_tempfile(), globals->nrows, globals->ncols,
	     srows, scols, sizeof(CELL), nseg) != TRUE)
	    G_fatal_error("Unable to create bounds temporary files");

	if (Rast_read_range(globals->bounds_map, globals->bounds_mapset, &range) != 1)
	    G_fatal_error(_("No min/max found in raster map <%s>"),
			  globals->bounds_map);
	Rast_get_range_min_max(&range, &globals->upper_bound,
				       &globals->lower_bound);

	if (Rast_is_c_null_value(&globals->upper_bound) ||
	    Rast_is_c_null_value(&globals->lower_bound)) {
	    
	    G_fatal_error(_("No min/max found in raster map <%s>"),
	                  globals->bounds_map);
	}

	bounds_fd = Rast_open_old(globals->bounds_map, globals->bounds_mapset);
	boundsbuf = Rast_allocate_c_buf();

	for (row = 0; row < globals->nrows; row++) {
	    Rast_get_c_row(bounds_fd, boundsbuf, row);
	    for (col = 0; col < globals->ncols; col++) {
		bounds_val = boundsbuf[col];
		if (FLAG_GET(globals->null_flag, row, col)) {
		    Rast_set_c_null_value(&bounds_val, 1);
		}
		else {
		    if (!Rast_is_c_null_value(&bounds_val)) {
			have_bounds = 1;
			if (globals->lower_bound > bounds_val)
			    globals->lower_bound = bounds_val;
			if (globals->upper_bound < bounds_val)
			    globals->upper_bound = bounds_val;
		    }
		}
		if (Segment_put(&globals->bounds_seg, &bounds_val, row, col) != 1)
		    G_fatal_error(_("Unable to write to temporary file"));
	    }
	}
	Rast_close(bounds_fd);
	G_free(boundsbuf);

	if (!have_bounds) {
	    G_warning(_("There are no boundary constraints in '%s'"), globals->bounds_map);
	    Rast_set_c_null_value(&globals->upper_bound, 1);
	    Rast_set_c_null_value(&globals->lower_bound, 1);
	    Segment_close(&globals->bounds_seg);
	    globals->bounds_map = NULL;
	    globals->bounds_mapset = NULL;
	}
    }
    else {
	G_debug(1, "no boundary constraint supplied.");
    }

    /* other info */
    globals->candidate_count = 0;	/* counter for remaining candidate pixels */

    /* Free memory */

    for (n = 0; n < Ref.nfiles; n++) {
	G_free(inbuf[n]);
	Rast_close(in_fd[n]);
    }

    globals->rs.sum = G_malloc(globals->datasize);
    globals->rs.mean = G_malloc(globals->datasize);

    globals->reg_tree = rgtree_create(globals->nbands, globals->datasize);
    globals->n_regions = s - 1;

    if (globals->seeds) {
	load_seeds(globals, srows, scols, nseg);
    }

    G_debug(1, "Number of initial regions: %d", globals->n_regions);

    G_free(inbuf);
    G_free(in_fd);
    G_free(fp_range);
    G_free(min);
    G_free(max);

    return TRUE;
}

//.........这里部分代码省略.........
    clump_buf = Rast_allocate_c_buf();
    
    /* get window size */
    G_get_window(&window);
    rows = window.rows;
    cols = window.cols;

    /* now get the data -- first pass */
    for (row = 0; row < rows; row++) {
	G_percent(row, rows, 2);
	Rast_get_d_row(fd_data, data_buf, row);
	Rast_get_c_row(fd_clump, clump_buf, row);
	for (col = 0; col < cols; col++) {
	    i = clump_buf[col];
	    if (i > max)
		G_fatal_error(_("Invalid category value %d (max=%d): row=%d col=%d"),
                              i, max, row, col);
	    if (i < 1) {
                G_debug(3, "row=%d col=%d: zero or negs ignored", row, col);
		continue;	/* ignore zeros and negs */
            }
	    if (Rast_is_d_null_value(&data_buf[col])) {
                G_debug(3, "row=%d col=%d: NULL ignored", row, col);
		continue;
            }
            
	    sum[i] += data_buf[col];
	    count[i]++;
	}
    }
    G_percent(1, 1, 1);
    
    /* free some buffer space */
    G_free(data_buf);
    G_free(clump_buf);

    /* data lists for centroids of clumps */
    e = (int *)G_malloc((max + 1) * sizeof(int));
    n = (int *)G_malloc((max + 1) * sizeof(int));

    i = centroids(fd_clump, e, n, 1, max);

    /* close raster maps */
    Rast_close(fd_data);
    Rast_close(fd_clump);
    
    /* got everything, now do output */
    if (centroidsmap) {
        G_message(_("Creating vector point map <%s>..."), centroidsmap);
        /* set comment */
	sprintf(buf, _("From '%s' on raster map <%s> using clumps from <%s>"),
                argv[0], datamap, clumpmap);
        Vect_set_comment(fd_centroids, buf);

        /* create attribute table */        
        Fi = Vect_default_field_info(fd_centroids, 1, NULL, GV_1TABLE);
        
        driver = db_start_driver_open_database(Fi->driver,
                                               Vect_subst_var(Fi->database, fd_centroids));
	if (driver == NULL) {
	    G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
			  Vect_subst_var(Fi->database, fd_centroids), Fi->driver);
	}
        db_set_error_handler_driver(driver);
        
	db_begin_transaction(driver);

//.........这里部分代码省略.........
    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);

    /* get entered parameters */
    RNET = input_RNET->answer;
    TEMPKA = input_TEMPKA->answer;
    PATM = input_PATM->answer;
    G0 = input_G0->answer;
    d_pt_alpha = atof(input_PT->answer);

    ETa = output->answer;

    /* open pointers to input raster files */
    infd_RNET = Rast_open_old(RNET, "");
    infd_TEMPKA = Rast_open_old(TEMPKA, "");
    infd_PATM = Rast_open_old(PATM, "");
    infd_G0 = Rast_open_old(G0, "");

    /* read headers of raster files */
    Rast_get_cellhd(RNET, "", &cellhd);
    Rast_get_cellhd(TEMPKA, "", &cellhd);
    Rast_get_cellhd(PATM, "", &cellhd);
    Rast_get_cellhd(G0, "", &cellhd);

    /* Allocate input buffer */
    inrast_RNET = Rast_allocate_d_buf();
    inrast_TEMPKA = Rast_allocate_d_buf();
    inrast_PATM = Rast_allocate_d_buf();
    inrast_G0 = Rast_allocate_d_buf();

    /* get rows and columns number of the current region */
    nrows = Rast_window_rows();
    ncols = Rast_window_cols();

    /* allocate output buffer */
    outrast = Rast_allocate_d_buf();

    /* open pointers to output raster files */
    outfd = Rast_open_new(ETa, DCELL_TYPE);

    /* start the loop through cells */
    for (row = 0; row < nrows; row++) {

	G_percent(row, nrows, 2);
	/* read input raster row into line buffer */
	Rast_get_d_row(infd_RNET, inrast_RNET, row);
	Rast_get_d_row(infd_TEMPKA, inrast_TEMPKA, row);
	Rast_get_d_row(infd_PATM, inrast_PATM, row);
	Rast_get_d_row(infd_G0, inrast_G0, row);

	for (col = 0; col < ncols; col++) {
	    /* read current cell from line buffer */
            d_rnet = ((DCELL *) inrast_RNET)[col];
            d_tempka = ((DCELL *) inrast_TEMPKA)[col];
            d_pt_patm = ((DCELL *) inrast_PATM)[col];
            d_g0 = ((DCELL *) inrast_G0)[col];

	    /*Delta_pt and Ghamma_pt */
	    d_pt_delta = pt_delta(d_tempka);
	    d_pt_ghamma = pt_ghamma(d_tempka, d_pt_patm);

	    /*Calculate ET */
	    d_daily_et =
		pt_daily_et(d_pt_alpha, d_pt_delta, d_pt_ghamma, d_rnet, d_g0,
			    d_tempka);
	    if (zero->answer && d_daily_et < 0)
		d_daily_et = 0.0;

	    /* write calculated ETP to output line buffer */
	    outrast[col] = d_daily_et;
	}

	/* write output line buffer to output raster file */
	Rast_put_d_row(outfd, outrast);
    }
    /* free buffers and close input maps */

    G_free(inrast_RNET);
    G_free(inrast_TEMPKA);
    G_free(inrast_PATM);
    G_free(inrast_G0);
    Rast_close(infd_RNET);
    Rast_close(infd_TEMPKA);
    Rast_close(infd_PATM);
    Rast_close(infd_G0);

    /* generate color table between -20 and 20 */
    Rast_make_rainbow_colors(&color, -20, 20);
    Rast_write_colors(ETa, G_mapset(), &color);

    Rast_short_history(ETa, "raster", &history);
    Rast_command_history(&history);
    Rast_write_history(ETa, &history);

    /* free buffers and close output map */
    G_free(outrast);
    Rast_close(outfd);

    return (EXIT_SUCCESS);
}

//.........这里部分代码省略.........
    int n_nulls = 0;

    /* Better use selected GDAL datatype instead of 
     * the best match with GRASS raster map types ? */

    if (maptype == FCELL_TYPE) {

	/* Source datatype understandable by GDAL */
	GDALDataType datatype = GDT_Float32;
	FCELL fnullval = (FCELL) nodataval;

	G_debug(1, "FCELL nodata val: %f", fnullval);

	for (row = 0; row < rows; row++) {

	    Rast_get_row(fd, bufer, row, maptype);
	    for (col = 0; col < cols; col++) {
		if (Rast_is_f_null_value(&((FCELL *) bufer)[col])) {
		    ((FCELL *) bufer)[col] = fnullval;
		    if (n_nulls == 0) {
			GDALSetRasterNoDataValue(hBand, nodataval);
		    }
		    n_nulls++;
		}
	    }

	    if (GDALRasterIO
		(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
		 0, 0) >= CE_Failure) {
		G_warning(_("Unable to write GDAL raster file"));
		return -1;
	    }
	    G_percent(row + 1, rows, 2);
	}
    }
    else if (maptype == DCELL_TYPE) {

	GDALDataType datatype = GDT_Float64;
	DCELL dnullval = (DCELL) nodataval;

	G_debug(1, "DCELL nodata val: %f", dnullval);

	for (row = 0; row < rows; row++) {

	    Rast_get_row(fd, bufer, row, maptype);
	    for (col = 0; col < cols; col++) {
		if (Rast_is_d_null_value(&((DCELL *) bufer)[col])) {
		    ((DCELL *) bufer)[col] = dnullval;
		    if (n_nulls == 0) {
			GDALSetRasterNoDataValue(hBand, nodataval);
		    }
		    n_nulls++;
		}
	    }

	    if (GDALRasterIO
		(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
		 0, 0) >= CE_Failure) {
		G_warning(_("Unable to write GDAL raster file"));
		return -1;
	    }
	    G_percent(row + 1, rows, 2);
	}
    }
    else {

	GDALDataType datatype = GDT_Int32;
	CELL inullval = (CELL) nodataval;

	G_debug(1, "CELL nodata val: %d", inullval);

	for (row = 0; row < rows; row++) {

	    Rast_get_row(fd, bufer, row, maptype);
	    for (col = 0; col < cols; col++) {
		if (Rast_is_c_null_value(&((CELL *) bufer)[col])) {
		    ((CELL *) bufer)[col] = inullval;
		    if (n_nulls == 0) {
			GDALSetRasterNoDataValue(hBand, nodataval);
		    }
		    n_nulls++;
		}
	    }

	    if (GDALRasterIO
		(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
		 0, 0) >= CE_Failure) {
		G_warning(_("Unable to write GDAL raster file"));
		return -1;
	    }
	    G_percent(row + 1, rows, 2);
	}
    }

    Rast_close(fd);

    G_free(bufer);

    return ret;
}

//.........这里部分代码省略.........

    /* Copy all rows */
    while (1) {
	int select;

	if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK) {
	    G_warning(_("Unable to fetch data from table <%s>"),
		      from_tblname);
	    db_close_cursor(&cursor);
	    db_close_database_shutdown_driver(to_driver);
	    if (from_driver != to_driver) {
		db_close_database_shutdown_driver(from_driver);
	    }
	    return DB_FAILED;
	}
	if (!more)
	    break;

	sprintf(buf, "insert into %s values ( ", to_tblname);
	db_set_string(&sql, buf);
	select = 1;
	for (col = 0; col < ncols; col++) {
	    column = db_get_table_column(table, col);
	    colname = db_get_column_name(column);
	    sqltype = db_get_column_sqltype(column);
	    ctype = db_sqltype_to_Ctype(sqltype);
	    value = db_get_column_value(column);

	    if (selcol && G_strcasecmp(colname, selcol) == 0) {
		if (db_test_value_isnull(value))
		    continue;
		if (!bsearch(&(value->i), ivalues, nvals, sizeof(int), cmp)) {
		    select = 0;
		    break;
		}
	    }
	    if (col > 0)
		db_append_string(&sql, ", ");
	    db_convert_value_to_string(value, sqltype, &value_string);
	    switch (ctype) {
	    case DB_C_TYPE_STRING:
	    case DB_C_TYPE_DATETIME:
		if (db_test_value_isnull(value)) {
		    db_append_string(&sql, "null");
		}
		else {
		    db_double_quote_string(&value_string);
		    db_append_string(&sql, "'");
		    db_append_string(&sql, db_get_string(&value_string));
		    db_append_string(&sql, "'");
		}
		break;
	    case DB_C_TYPE_INT:
	    case DB_C_TYPE_DOUBLE:
		if (db_test_value_isnull(value)) {
		    db_append_string(&sql, "null");
		}
		else {
		    db_append_string(&sql, db_get_string(&value_string));
		}
		break;
	    default:
		G_warning(_("Unknown column type (column <%s>)"),
			  colname);
		db_close_cursor(&cursor);
		db_close_database_shutdown_driver(to_driver);
		if (from_driver != to_driver) {
		    db_close_database_shutdown_driver(from_driver);
		}
		return DB_FAILED;
	    }
	}
	if (!select)
	    continue;
	db_append_string(&sql, ")");
	G_debug(3, db_get_string(&sql));
	if (db_execute_immediate(to_driver, &sql) != DB_OK) {
	    G_warning("Unable to insert new record: '%s'",
		      db_get_string(&sql));
	    db_close_cursor(&cursor);
	    db_close_database_shutdown_driver(to_driver);
	    if (from_driver != to_driver) {
		db_close_database_shutdown_driver(from_driver);
	    }
	    return DB_FAILED;
	}
    }
    if (selcol)
	G_free(ivalues);
    G_debug(3, "Table copy OK");

    db_close_cursor(&cursor);
    db_commit_transaction(to_driver);
    db_close_database_shutdown_driver(to_driver);
    if (from_driver != to_driver) {
	db_close_database_shutdown_driver(from_driver);
    }

    return DB_OK;
}

int read_rgb(char *key, char *data)
{
    char names[3][100];
    char fullname[100];
    int i;

    if (sscanf(data, "%s %s %s", names[0], names[1], names[2]) != 3) {
	error(key, data, "illegal request (rgb)");
	return 0;
    }

    PS.do_raster = 0;
    PS.do_colortable = 0;
    if (PS.cell_fd >= 0) {
	G_close_cell(PS.cell_fd);
	G_free(PS.cell_name);
	G_free(PS.cell_mapset);
	G_free_colors(&PS.colors);
	PS.cell_fd = -1;
    }

    /* initialize group structure (for compatibility with PS_raster_plot()) */
    I_init_group_ref(&grp.ref);

    /*
     * not relevant here
     *
     if (I_get_group_ref(grp.group_name, &grp.ref) == 0)
     G_fatal_error(_("Can't get group information"));
     */

    grp.group_name = "RGB Group";

    /* get file names for R, G, & B */

    for (i = 0; i < 3; i++) {
	char *mapset, *name, *p;

	name = names[i];

	p = strchr(name, '@');
	if (p) {
	    *p = '\0';
	    mapset = p + 1;
	}
	else {
	    mapset = G_find_file2("cell", name, "");
	    if (!mapset) {
		error(name, "", "not found");
		return 0;
	    }
	}

	grp.name[i] = G_store(name);
	grp.mapset[i] = G_store(mapset);

	/* read in colors */
	if (G_read_colors(grp.name[i], grp.mapset[i], &(grp.colors[i])) == -1) {
	    sprintf(fullname, "%s in %s", grp.name[i], grp.mapset[i]);
	    error(fullname, "", "can't read color table");
	    return 0;
	}

	/* open raster maps for reading */
	if ((grp.fd[i] = G_open_cell_old(grp.name[i], grp.mapset[i])) < 0) {
	    sprintf(fullname, "%s in %s", grp.name[i], grp.mapset[i]);
	    error(fullname, "", "can't open raster map");
	    G_free_colors(&(grp.colors[i]));
	    return 0;
	}
    }

    strcpy(PS.celltitle, grp.group_name);
    G_strip(PS.celltitle);
    return 1;
}

/*!
 *  \brief Extract a cell value from raster map (cubic interpolation).
 *
 *