自学教程：C++ G_free函数代码示例

//.........这里部分代码省略.........    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). * *  Extract a cell value from raster map at given northing and easting *  with a sampled 3x3 window using a cubic 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_cubic(int fd,			    const struct Cell_head * window,			    struct Categories * cats,			    double north, double east, int usedesc){    int i, j, row, col;    double grid[4][4];    DCELL *rows[4];    double frow, fcol, trow, tcol;    DCELL result;    for (i = 0; i < 4; i++)	rows[i] = Rast_allocate_d_buf();    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 - 1.5);    col = (int)floor(fcol - 1.5);    trow = frow - row - 1.5;    tcol = fcol - col - 1.5;    if (row < 0 || row + 3 >= Rast_window_rows() ||	col < 0 || col + 3 >= Rast_window_cols()) {	Rast_set_d_null_value(&result, 1);	goto done;    }    for (i = 0; i < 4; i++)	Rast_get_d_row(fd, rows[i], row + i);    for (i = 0; i < 4; i++)	for (j = 0; j < 4; j++)	    if (Rast_is_d_null_value(&rows[i][col + j])) {		Rast_set_d_null_value(&result, 1);		goto done;	    }    /*     * now were ready to do cubic interpolation over     * arow[col], arow[col+1], arow[col+2], arow[col+3],     * brow[col], brow[col+1], brow[col+2], brow[col+3],     * crow[col], crow[col+1], crow[col+2], crow[col+3],     * drow[col], drow[col+1], drow[col+2], drow[col+3],     */    if (usedesc) {	char *buf;	for (i = 0; i < 4; i++) {	    for (j = 0; j < 4; j++) {		G_squeeze(buf =			  Rast_get_c_cat((CELL *) & (rows[i][col + j]),					 cats));		grid[i][j] = scancatlabel(buf);	    }	}    }    else {	for (i = 0; i < 4; i++)	    for (j = 0; j < 4; j++)		grid[i][j] = rows[i][col + j];    }    result = Rast_interp_bicubic(tcol, trow,				 grid[0][0], grid[0][1], grid[0][2],				 grid[0][3], grid[1][0], grid[1][1],				 grid[1][2], grid[1][3], grid[2][0],				 grid[2][1], grid[2][2], grid[2][3],				 grid[3][0], grid[3][1], grid[3][2],				 grid[3][3]);  done:    for (i = 0; i < 4; i++)	G_free(rows[i]);    return result;}

int avl_add(avl_tree * root, const generic_cell k, const long n){    avl_node *p = NULL;    avl_node *node_temp = NULL;    avl_node *critical = NULL;    int d = 0;    int pos1 = 0, pos2 = 0;    int rotation = 0;    if ((root == NULL) || (*root == NULL)) {	G_fatal_error("/navl.c: avl_add: param NULL");	return AVL_ERR;    }    /* search position where insert the new node */    node_temp = avl_individua(*root, k, &p, &d);    if (node_temp != NULL) {	node_temp->counter = node_temp->counter + n;	return AVL_PRES;    }    node_temp = avl_make(k, n);    if (node_temp == NULL) {	G_fatal_error("/navl.c:  avl_add: create node error");	return AVL_ERR;    }    /* link the new node */    node_temp->father = p;    if (d == -1) {	p->left_child = node_temp;    }    else {	if (d == 1) {	    p->right_child = node_temp;	}	else {	    G_free(node_temp);	    G_fatal_error("avl.c: avl_add: new node position unknown");	    return AVL_ERR;	}    }    /* if it's necessary balance the tree */    critical = critical_node(node_temp, &pos1, &pos2, NULL);    if (critical == NULL)	return AVL_ADD;    rotation = (pos1 * 10) + pos2;    switch (rotation) {    case AVL_SS:	avl_rotation_ll(critical);	break;    case AVL_SD:	avl_rotation_lr(critical);	break;    case AVL_DS:	avl_rotation_rl(critical);	break;    case AVL_DD:	avl_rotation_rr(critical);	break;    default:	G_fatal_error("avl, avl_add: balancing error/n");	return AVL_ERR;    }    /* if after rotation the root is changed modufy the pointer to the root */    while ((*root)->father != NULL)	*root = (*root)->father;    return AVL_ADD;}

//.........这里部分代码省略.........		    else {			if (wat_map_type == CELL_TYPE) {			    wat_value = *((CELL *)watptr);			}			else if (wat_map_type == FCELL_TYPE) {			    wat_value = *((FCELL *)watptr);			}			else if (wat_map_type == DCELL_TYPE) {			    wat_value = *((DCELL *)watptr);			}		    }		}		else {		    wat_value = 1;		}	    }	    wabuf[c].wat = wat_value;	    wabuf[c].ele = alt_value;	    alt_value_buf[c] = alt_value;	    ptr = G_incr_void_ptr(ptr, ele_size);	    if (run_flag) {		watptr = G_incr_void_ptr(watptr, wat_size);	    }	}	seg_put_row(&watalt, (char *) wabuf, r);	seg_put_row(&aspflag, (char *)afbuf, r);		if (er_flag) {	    cseg_put_row(&r_h, alt_value_buf, r);	}    }    G_percent(nrows, nrows, 1);    /* finish it */    Rast_close(ele_fd);    G_free(wabuf);    G_free(afbuf);        if (run_flag) {	Rast_close(wat_fd);	G_free(watbuf);    }    MASK_flag = (do_points < nrows * ncols);        /* do RUSLE */    if (er_flag) {	if (ob_flag) {	    fd = Rast_open_old(ob_name, "");	    buf = Rast_allocate_c_buf();	    for (r = 0; r < nrows; r++) {		G_percent(r, nrows, 1);		Rast_get_c_row(fd, buf, r);		for (c = 0; c < ncols; c++) {		    block_value = buf[c];		    if (!Rast_is_c_null_value(&block_value) && block_value) {			seg_get(&aspflag, (char *)&af, r, c);			FLAG_SET(af.flag, RUSLEBLOCKFLAG);			seg_put(&aspflag, (char *)&af, r, c);		    }		}	    }	    G_percent(nrows, nrows, 1);    /* finish it */	    Rast_close(fd);	    G_free(buf);	}	if (ril_flag) {

/* douglas-peucker algorithm which simplifies a line to a line with * at most reduction% of points. * returns the number of points in the output line. It is approx  * reduction/100 * Points->n_points. */int douglas_peucker_reduction(struct line_pnts *Points, double thresh,			      double reduction, int with_z){    int i;    int n = Points->n_points;    /* the maximum number of points  which may be      * included in the output */    int nexp = n * (reduction / (double)100.0);    /* line too short */    if (n < 3)	return n;    /* indicates which point were selected by the algorithm */    int *sel;    sel = G_calloc(sizeof(int), n);    if (sel == NULL) {	G_fatal_error(_("Out of memory"));	return n;    }    /* array used for storing the indices of line segments+furthest point */    int *index;    index = G_malloc(sizeof(int) * 3 * n);    if (index == NULL) {	G_fatal_error(_("Out of memory"));	G_free(sel);	return n;    }    int indices;    indices = 0;    /* preserve first and last point */    sel[0] = sel[n - 1] = 1;    nexp -= 2;    thresh *= thresh;    double d;    int mid = get_furthest(Points, 0, n - 1, with_z, &d);    int em;    /* priority queue of line segments,     * key is the distance of the furthest point */    binary_heap pq;    if (!binary_heap_init(n, &pq)) {	G_fatal_error(_("Out of memory"));	G_free(sel);	G_free(index);	return n;    }    if (d > thresh) {	index[0] = 0;	index[1] = n - 1;	index[2] = mid;	binary_heap_push(d, 0, &pq);	indices = 3;    }    /* while we can add new points and queue is non-empty */    while (nexp > 0) {	/* empty heap */	if (!binary_heap_extract_max(&pq, &em))	    break;	int left = index[em];	int right = index[em + 1];	int furt = index[em + 2];	/*mark the furthest point */	sel[furt] = 1;	nexp--;	/* consider left and right segment */	mid = get_furthest(Points, left, furt, with_z, &d);	if (d > thresh) {	    binary_heap_push(d, indices, &pq);	    index[indices++] = left;	    index[indices++] = furt;	    index[indices++] = mid;	}	mid = get_furthest(Points, furt, right, with_z, &d);	if (d > thresh) {	    binary_heap_push(d, indices, &pq);	    index[indices++] = furt;	    index[indices++] = right;	    index[indices++] = mid;	}//.........这里部分代码省略.........

//.........这里部分代码省略.........	   row buffer to 0; this row buffer	   will hold one row of the clipped	   raster map.  Then read row i of the	   map and the corresponding null values	   into tmp and nulltmp buffers */	switch (data_type) {	case CELL_TYPE:	    G_zero_raster_buf(tmp, data_type);	    G_get_raster_row(finput, tmp, i, CELL_TYPE);	    break;	case FCELL_TYPE:	    G_zero_raster_buf(ftmp, data_type);	    G_get_raster_row(finput, ftmp, i, FCELL_TYPE);	    break;	case DCELL_TYPE:	    G_zero_raster_buf(dtmp, data_type);	    G_get_raster_row(finput, dtmp, i, DCELL_TYPE);	    break;	}	G_get_null_value_row(finput, nulltmp, i);	/* for all the columns one by one */	for (j = col0; j < col0 + ncols; j++) {	    /* if circles are used for sampling */	    if ((int)radius) {		dist = sqrt(((double)i - center_row) *			    ((double)i - center_row) +			    ((double)j - center_col) *			    ((double)j - center_col));		/* copy the contents of tmp into the		   appropriate cell in buf */		if (dist < radius) {		    switch (data_type) {		    case CELL_TYPE:			*(*(buf + i + 1 - row0) + j + 1 - col0) = *(tmp + j);			break;		    case FCELL_TYPE:			*(*(buf + i + 1 - row0) + j + 1 - col0) = *(ftmp + j);			break;		    case DCELL_TYPE:			*(*(buf + i + 1 - row0) + j + 1 - col0) = *(dtmp + j);			break;		    }		    *(*(null_buf + i + 1 - row0) + j + 1 - col0) =			*(nulltmp + j);		}	    }	    /* if circles are not used and	       if the choice is not "by region" or	       if this column is in region "index" */	    else if (choice->wrum != 'r' || *(tmp1 + j) == index) {		/* copy the contents of the correct tmp		   into the appropriate cell in the buf		   and the corresponding null values into		   the appropriate cell in null_buf */		switch (data_type) {		case CELL_TYPE:		    *(*(buf + i + 1 - row0) + j + 1 - col0) = *(tmp + j);		    break;		case FCELL_TYPE:		    *(*(buf + i + 1 - row0) + j + 1 - col0) = *(ftmp + j);		    break;		case DCELL_TYPE:		    *(*(buf + i + 1 - row0) + j + 1 - col0) = *(dtmp + j);		    break;		}		*(*(null_buf + i + 1 - row0) + j + 1 - col0) = *(nulltmp + j);	    }	}    }    switch (data_type) {    case CELL_TYPE:	G_free(tmp);	break;    case FCELL_TYPE:	G_free(ftmp);	break;    case DCELL_TYPE:	G_free(dtmp);	break;    }    if (choice->wrum == 'r') {	G_free(tmp1);	G_close_cell(fr);    }    G_free(nulltmp);    return;}

//.........这里部分代码省略.........    current =        (dglEdgesetTraverser_s *) G_calloc(nnodes + 1,                                           sizeof(dglEdgesetTraverser_s));    tin = (int *)G_calloc(nnodes + 1, sizeof(int));    min_tin = (int *)G_calloc(nnodes + 1, sizeof(int));    parent = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));    stack = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));    current_edge = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));    mark = (int *)G_calloc(nnodes + 1, sizeof(int));    if (!tin || !min_tin || !parent || !stack || !current || !mark) {        G_fatal_error(_("Out of memory"));        return -1;    }    for (i = 1; i <= nnodes; i++) {        dglEdgeset_T_Initialize(&current[i], graph,                                dglNodeGet_OutEdgeset(graph,                                        dglGetNode(graph, i)));        current_edge[i] = dglEdgeset_T_First(&current[i]);        tin[i] = mark[i] = 0;    }    dglNode_T_Initialize(&nt, graph);    time = 0;    for (current_node = dglNode_T_First(&nt); current_node;            current_node = dglNode_T_Next(&nt)) {        dglInt32_t current_id = dglNodeGet_Id(graph, current_node);        if (tin[current_id] == 0) {            int children = 0;	/*number of subtrees rooted at the root/current_node */            stack[0] = current_node;            stack_size = 1;            parent[current_id] = NULL;            while (stack_size) {                dglInt32_t *node = stack[stack_size - 1];                dglInt32_t node_id = dglNodeGet_Id(graph, node);                if (tin[node_id] == 0)	/*vertex visited for the first time */                    min_tin[node_id] = tin[node_id] = ++time;                else {		/*return from the recursion */                    dglInt32_t to = dglNodeGet_Id(graph,                                                  dglEdgeGet_Tail(graph,                                                          current_edge                                                          [node_id]));                    if (min_tin[to] >= tin[node_id])	/*no path from the subtree above the current node */                        mark[node_id] = 1;	/*so the current node must be an articulation point */                    if (min_tin[to] < min_tin[node_id])                        min_tin[node_id] = min_tin[to];                    current_edge[node_id] = dglEdgeset_T_Next(&current[node_id]);	/*proceed to the next edge */                }                /*try next edges */                for (; current_edge[node_id]; current_edge[node_id] = dglEdgeset_T_Next(&current[node_id])) {                    dglInt32_t *to =                        dglEdgeGet_Tail(graph, current_edge[node_id]);                    if (to == parent[node_id])                        continue;	/*skip parent */                    int to_id = dglNodeGet_Id(graph, to);                    if (tin[to_id]) {	/*back edge, cannot be a bridge/articualtion point */                        if (tin[to_id] < min_tin[node_id])                            min_tin[node_id] = tin[to_id];                    }                    else {	/*forward edge */                        if (node_id == current_id)                            children++;	/*if root, increase number of children */                        parent[to_id] = node;                        stack[stack_size++] = to;                        break;                    }                }                if (!current_edge[node_id])                    stack_size--;	/*current node completely processed */            }            if (children > 1)                mark[current_id] = 1;	/*if the root has more than 1 subtrees rooted at it, then it is an					 * articulation point */        }    }    for (i = 1; i <= nnodes; i++)        if (mark[i]) {            points++;            Vect_list_append(articulation_list, i);        }    dglNode_T_Release(&nt);    for (i = 1; i <= nnodes; i++)        dglEdgeset_T_Release(&current[i]);    G_free(current);    G_free(tin);    G_free(min_tin);    G_free(parent);    G_free(stack);    G_free(current_edge);    return points;}

void cell_clip_drv(int col0, int row0, int ncols, int nrows, double **value,		   int index, int cntwhole, float radius){    register int i, j;    int cnt = 0, p;    double *rich, *richtmp;    char *name, *mapset;    DCELL **buf;    DCELL **null_buf;    RASTER_MAP_TYPE data_type;    /*       col0 = starting column for area to be clipped       row0 = starting row for area to be clipped       ncols = number of columns in area to be clipped       nrows = number of rows in area to be clipped       value =       index = number of the region to be clipped, if there's a region map       buf = pointer to array containing the clipped area, a smaller area       than the original raster map to be read from finput                            printf("h2/n");       pat = pointer to array containing the map of patch numbers       cor = pointer to array containing the map of interior area     */    name = choice->fn;    mapset = G_mapset();    data_type = G_raster_map_type(name, mapset);    /* dynamically allocate storage for the       buffer that will hold the contents of       the window */    buf = (DCELL **) G_calloc(nrows + 3, sizeof(DCELL *));    for (i = 0; i < nrows + 3; i++) {	buf[i] = (DCELL *) G_calloc(ncols + 3, sizeof(DCELL));    }    /* dynamically allocate storage for the       buffer that will hold the null values for       the clipped area */    null_buf = (DCELL **) G_calloc(nrows + 3, sizeof(DCELL *));    for (i = 0; i < nrows + 3; i++)	null_buf[i] = (DCELL *) G_calloc(ncols + 3, sizeof(DCELL));    /* call the cell_clip routine */    cell_clip(buf, null_buf, row0, col0, nrows, ncols, index, radius);    /* dynamically allocate memory for       the richness array */    richtmp = (double *)G_calloc(MAX, sizeof(double));    /* go through the sampling area       pixel by pixel */    for (i = 1; i < nrows + 1; i++) {	for (j = 1; j < ncols + 1; j++) {	    /* if buf[i][j] is not a null value,	       call get_rich to tally up the	       number of different attributes	       in the sampling area and fill	       the richness array with those	       attributes */	    if ((buf[i][j] || buf[i][j] == 0.0) && null_buf[i][j] == 0.0) {		/*printf("buf[%d][%d] = %f/n",i,j,buf[i][j]); */		get_rich(buf[i][j], richtmp, &cnt);	    }	}    }    if (cnt) {	rich = (double *)G_calloc(cnt, sizeof(double));	for (i = 0; i < cnt; i++) {	    rich[i] = richtmp[i];	}	G_free(richtmp);	/* call ANSI C runtime library	   function qsort to sort the 	   richness array into ascending	   order */	qsort(rich, cnt, sizeof(double), compar);	/* moving window */	if (choice->wrum == 'm') {	    if (is_not_empty_buffer(buf, null_buf, nrows + 1, ncols + 1)) {		if (center_is_not_null(buf, null_buf, nrows, ncols))		    mv_texture(nrows, ncols, buf, null_buf, value, index,			       rich, cnt, cntwhole);		else {//.........这里部分代码省略.........

static int write_smooth_bnd(struct COOR *line_begin, struct COOR *line_end,	/* start and end point of line */			    int n	/* number of points to write */    ){    static struct line_pnts *points = NULL;    double x, y;    double dx, dy;    int idx, idy;    struct COOR *p, *last;    int i, total;    if (!points)	points = Vect_new_line_struct();    Vect_reset_line(points);    n++;			/* %% 6.4.88 */    p = line_begin;    /* allocate the arrays and get the first point */    y = cell_head.north - (double)p->row * cell_head.ns_res;    x = cell_head.west + (double)p->col * cell_head.ew_res;    Vect_append_point(points, x, y, 0.0);    /* generate the list of smoothed points, may be duplicate points */    total = 1;    for (i = 1; i < n; i++) {	if (i < 10)	    G_debug(3, " row: %d col: %d/n", p->row, p->col);	last = p;	if ((p = move(p)) == NULPTR) {	/* this should NEVER happen */	    G_debug(3, "write_line:  line terminated unexpectedly/n");	    G_debug(3, "  previous (%d) point %p (%d,%d,%d) %p %p/n",		    direction, last, last->row, last->col, last->node,		    last->fptr, last->bptr);	    exit(EXIT_FAILURE);	}	idy = (p->row - last->row);	idx = (p->col - last->col);	dy = (idy > 0) ? 0.5 : ((idy < 0) ? -0.5 : 0.0);	/* dy = 0.0, 0.5, or -0.5 */	dx = (idx > 0) ? 0.5 : ((idx < 0) ? -0.5 : 0.0);	/* dx = 0.0, 0.5, or -0.5 */	y = cell_head.north - (last->row + dy) * cell_head.ns_res;	x = cell_head.west + (last->col + dx) * cell_head.ew_res;	total++;	Vect_append_point(points, x, y, 0.0);	y = cell_head.north - (p->row - dy) * cell_head.ns_res;	x = cell_head.west + (p->col - dx) * cell_head.ew_res;	total++;	Vect_append_point(points, x, y, 0.0);    }				/* end of for i */    y = cell_head.north - (double)p->row * cell_head.ns_res;    x = cell_head.west + (double)p->col * cell_head.ew_res;    total++;    Vect_append_point(points, x, y, 0.0);    /* strip out the duplicate points from the list */    Vect_line_prune(points);    G_debug(3, "removed duplicates: %d", total - points->n_points);    /* write files */    Vect_write_line(&Map, GV_BOUNDARY, points, Cats);    /* now free all the pointers */    p = line_begin;    for (i = 1; i < n; i++) {	if (i < 10)	    G_debug(3, " row: %d col: %d/n", p->row, p->col);	last = p;	if ((p = move(p)) == NULPTR)	    break;	if (last == p)	    break;	if (last->fptr != NULPTR)	    if (last->fptr->fptr == last)		last->fptr->fptr = NULPTR;	/* now it can already ne NULL */	if (last->fptr != NULPTR)	    if (last->fptr->bptr == last)		last->fptr->bptr = NULPTR;	if (last->bptr != NULPTR)	    if (last->bptr->fptr == last)		last->bptr->fptr = NULPTR;	if (last->bptr != NULPTR)	    if (last->bptr->bptr == last)		last->bptr->bptr = NULPTR;	G_free(last);    }				/* end of for i */    if (p != NULPTR)	G_free(p);//.........这里部分代码省略.........

//.........这里部分代码省略.........			    "  pnt2 = %d dist2 = %f iso2 = %d max iso1 = %d",			    i, pnts2[i].distance, pnts2[i].iso,			    pnts1[npnts1 - 1].iso);		    if (pnts1[npnts1 - 1].iso < pnts2[i].iso) {			G_debug(3, "    -> cut here");			npnts2 = i;			break;		    }		}	    }	    G_debug(3, "  npnts2 cut = %d", npnts2);	    /* Biggest cost shoud be equal if exist (npnts > 0). Cut out overlapping segments,	     *  this can cut only points on line but not first points */	    if (npnts1 > 1 && npnts2 > 1) {		while (npnts1 > 1 && npnts2 > 1) {		    if (pnts1[npnts1 - 1].distance >= pnts2[npnts2 - 1].distance) {	/* overlap */			npnts1--;			npnts2--;		    }		    else {			break;		    }		}	    }	    G_debug(3, "  npnts1 2. cut = %d", npnts1);	    G_debug(3, "  npnts2 2. cut = %d", npnts2);	    /* Now we have points in both directions which may not overlap, npoints in one	     *  direction may be 0 but not both */	    /* Join both arrays, iso of point is for next segment (point is at the beginning) */	    /* In case npnts1 == 0 add point at distance 0 */	    if (npnts1 == 0) {		G_debug(3,			"  npnts1 = 0 -> add first at distance 0, cat = %d",			pnts2[npnts2 - 1].iso);		pnts1[0].iso = pnts2[npnts2 - 1].iso;	/* use last point iso in reverse direction */		pnts1[0].distance = 0;		npnts1++;	    }	    for (i = npnts2 - 1; i >= 0; i--) {		/* Check if identical */		if (pnts1[npnts1 - 1].distance == pnts2[i].distance)		    continue;		if (npnts1 == apnts1) {		    apnts1 += 1;		    pnts1 =			(ISOPOINT *) G_realloc(pnts1,					       apnts1 * sizeof(ISOPOINT));		}		pnts1[npnts1].iso = pnts2[i].iso - 1;	/* last may be -1, but it is not used */		pnts1[npnts1].distance = pnts2[i].distance;		npnts1++;	    }	    /* In case npnts2 == 0 add point at the end */	    if (npnts2 == 0) {		pnts1[npnts1].iso = 0;	/* not used */		pnts1[npnts1].distance = l;		npnts1++;	    }	    /* Create line segments. */	    for (i = 1; i < npnts1; i++) {		cat = pnts1[i - 1].iso + 1;		G_debug(3, "  segment %f - %f cat %d", pnts1[i - 1].distance,			pnts1[i].distance, cat);		ret =		    Vect_line_segment(Points, pnts1[i - 1].distance,				      pnts1[i].distance, SPoints);		if (ret == 0) {		    G_warning(_			      ("Cannot get line segment, segment out of line"));		}		else {		    Vect_reset_cats(Cats);		    Vect_cat_set(Cats, 1, cat);		    Vect_write_line(&Out, ltype, SPoints, Cats);		}	    }	}	else {	    /* arc is not reachable */	    G_debug(3, "  -> arc is not reachable");	    Vect_reset_cats(Cats);	    Vect_write_line(&Out, ltype, Points, Cats);	}    }    Vect_build(&Out);    /* Free, ... */    G_free(Nodes);    G_free(Centers);    Vect_close(&Map);    Vect_close(&Out);    exit(EXIT_SUCCESS);}