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

/* ************************************************************************* */void write_vtk_vector_data(void *map_x, void *map_y, void *map_z,                           FILE * fp, const char *varname,                           RASTER3D_Region region, int dp){    double value = 0;    int x, y, z, percentage, k;    int rows, cols, depths;    int typeIntern[3];    void *mapvect = NULL;    G_debug(3, "write_vtk_vector_data: Writing vector data");    rows = region.rows;    cols = region.cols;    depths = region.depths;    typeIntern[0] = Rast3d_tile_type_map(map_x);    typeIntern[1] = Rast3d_tile_type_map(map_y);    typeIntern[2] = Rast3d_tile_type_map(map_z);    percentage = 0;    /********************** WRITE VECTOR DATA; CELL OR POINT ****************/    fprintf(fp, "VECTORS %s float/n", varname);    for (z = 0; z < depths; z++) { /*From the bottom to the top */        for (y = 0; y < rows; y++) {            G_percent(percentage, (rows * depths - 1), 10);            percentage++;            for (x = 0; x < cols; x++) {                for (k = 0; k < 3; k++) {                    if (k == 0)                        mapvect = map_x;                    if (k == 1)                        mapvect = map_y;                    if (k == 2)                        mapvect = map_z;                    /* In case of structured grid data, the point/cell coordinates                       are computed based on the default north->south raster3d coordinate system.                       We need to compute south -> north ordering for image data.                     */                    if (!param.structgrid->answer)                        value =                        get_g3d_raster_value_as_double(mapvect, x, rows - y - 1, z,                                                       typeIntern[k],                                                       0.0);                    else                        value =                        get_g3d_raster_value_as_double(mapvect, x, y, z,                                                       typeIntern[k],                                                       0.0);                    fprintf(fp, "%.*f ", dp, value);                }                fprintf(fp, "/n");            }        }    }    return;}

int do_cum(void){    int r, c, dr, dc;    char asp_val, asp_val_down;    char is_swale, this_flag_value, flag_value;    DCELL value, valued;    POINT point;    int killer, threshold;    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };    WAT_ALT wa, wadown;    G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));    if (bas_thres <= 0)	threshold = 60;    else	threshold = bas_thres;    for (killer = 0; killer < do_points; killer++) {	G_percent(killer, do_points, 1);	seg_get(&astar_pts, (char *)&point, 0, killer);	r = point.r;	c = point.c;	bseg_get(&asp, &asp_val, r, c);	if (asp_val) {	    dr = r + asp_r[ABS(asp_val)];	    dc = c + asp_c[ABS(asp_val)];	}	/* skip user-defined depressions */	else	    dr = dc = -1;	bseg_get(&bitflags, &this_flag_value, r, c);	FLAG_UNSET(this_flag_value, WORKEDFLAG);	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) {	    /* TODO: do not distribute flow along edges, this causes artifacts */	    seg_get(&watalt, (char *)&wa, r, c);	    value = wa.wat;	    is_swale = FLAG_GET(this_flag_value, SWALEFLAG);	    if (fabs(value) >= threshold && !is_swale) {		is_swale = 1;		FLAG_SET(this_flag_value, SWALEFLAG);	    }	    seg_get(&watalt, (char *)&wadown, dr, dc);	    valued = wadown.wat;	    if (value > 0) {		if (valued > 0)		    valued += value;		else		    valued -= value;	    }	    else {		if (valued < 0)		    valued += value;		else		    valued = value - valued;	    }	    wadown.wat = valued;	    seg_put(&watalt, (char *)&wadown, dr, dc);	    /* update asp for depression */	    if (is_swale || fabs(valued) >= threshold) {		bseg_get(&bitflags, &flag_value, dr, dc);		FLAG_SET(flag_value, SWALEFLAG);		bseg_put(&bitflags, &flag_value, dr, dc);		is_swale = 1;	    }	    else {		if (er_flag && !is_swale && !FLAG_GET(this_flag_value, RUSLEBLOCKFLAG))		    slope_length(r, c, dr, dc);	    }	}	bseg_put(&bitflags, &this_flag_value, r, c);    }    G_percent(do_points, do_points, 1);	/* finish it */    seg_close(&astar_pts);    return 0;}

/* write_area - make table of area equivalences and write attribute file */int write_area(struct area_table *a_list,	/* list of areas */	       struct equiv_table *e_list,	/* list of equivalences between areas */	       int n_areas,	/* lengths of e_list, a_list */	       int n_equiv){    struct line_pnts *points = Vect_new_line_struct();    int n, i;    struct area_table *p;    char *temp_buf;    int *equivs;    int cat;    int catNum;    double x, y;    equivs = NULL;    total_areas = 0;    if (n_equiv < n_areas) {	equivs = (int *)G_malloc(n_areas * sizeof(int));	n = n_equiv;    }    else {	equivs = (int *)G_malloc(n_equiv * sizeof(int));	n = n_areas;    }    for (i = 0; i < n; i++) {	if ((e_list + i)->mapped)	    equivs[i] = (e_list + i)->where;	else {	    total_areas++;	    equivs[i] = i;	}    }    if (n < n_areas) {	for (i = n; i < n_areas; i++) {	    total_areas++;	    equivs[i] = i;	}    }    catNum = 1;    G_important_message(_("Writing areas..."));    for (i = 0, p = a_list; i < n_areas; i++, p++) {	G_percent(i, n_areas, 3);	if (equivs[i] == i && p->width > 0 && !Rast_is_d_null_value(&(p->cat))) {	    char buf[1000];	    if (value_flag) {	/* raster value */		cat = (int)p->cat;	    }	    else {		/* sequence */		cat = catNum;		catNum++;	    }	    x = cell_head.west + (p->col +				  (p->width / 2.0)) * cell_head.ew_res;	    y = cell_head.north - (p->row + 0.5) * cell_head.ns_res;	    switch (data_type) {	    case CELL_TYPE:		G_debug(3,			"vector x = %.3f, y = %.3f, cat = %d; raster cat = %d",			x, y, cat, (int)p->cat);		break;	    case FCELL_TYPE:		G_debug(3,			"vector x = %.3f, y = %.3f, cat = %d; raster cat = %f",			x, y, cat, (float)p->cat);		break;	    case DCELL_TYPE:		G_debug(3,			"vector x = %.3f, y = %.3f, cat = %d; raster cat = %lf",			x, y, cat, p->cat);		break;	    }	    Vect_reset_line(points);	    Vect_append_point(points, x, y, 0.0);	    Vect_reset_cats(Cats);	    Vect_cat_set(Cats, 1, cat);	    Vect_write_line(&Map, GV_CENTROID, points, Cats);	    if (driver != NULL && !value_flag) {		sprintf(buf, "insert into %s values (%d, ", Fi->table, cat);		db_set_string(&sql, buf);		switch (data_type) {		case CELL_TYPE:		    sprintf(buf, "%d", (int)p->cat);		    break;		case FCELL_TYPE:		case DCELL_TYPE:		    sprintf(buf, "%f", p->cat);		    break;//.........这里部分代码省略.........

/* * check_stats() - Check and update statistics  * * RETURN: 0 on success / 1 on failure */int check_stats(const char *name){    RASTER_MAP_TYPE data_type;    struct Histogram histogram;    struct Categories cats;    struct Range range;    struct FPRange fprange;    int i, histo_num;    int cats_ok;    int max;    data_type = Rast_map_type(name, "");    G_message(_("Updating statistics for [%s]..."), name);    if (!do_histogram(name))	return 1;    if (Rast_read_histogram(name, "", &histogram) <= 0)	return 1;    /* Init histogram range */    if (data_type == CELL_TYPE)	Rast_init_range(&range);    else	Rast_init_fp_range(&fprange);    G_message(_("Updating histogram range..."));    i = histo_num = Rast_get_histogram_num(&histogram);    while (i >= 0) {	G_percent(i, histo_num, 2);	if (data_type == CELL_TYPE)	    Rast_update_range(Rast_get_histogram_cat(i--, &histogram), &range);	else	    Rast_update_fp_range((DCELL) Rast_get_histogram_cat(i--, &histogram),			      &fprange);    }    /* Write histogram range */    if (data_type == CELL_TYPE)	Rast_write_range(name, &range);    else	Rast_write_fp_range(name, &fprange);    /* Get category status and max */    cats_ok = (Rast_read_cats(name, "", &cats) >= 0);    max = (data_type == CELL_TYPE ? range.max : fprange.max);    /* Further category checks */    if (!cats_ok)	Rast_init_cats("", &cats);    else if (cats.num != max) {	cats.num = max;	cats_ok = 0;    }    /* Update categories if needed */    if (!cats_ok) {	G_message(_("Updating the number of categories for [%s]..."), name);	Rast_write_cats(name, &cats);    }    Rast_free_histogram(&histogram);    Rast_free_cats(&cats);    return 0;}

//.........这里部分代码省略.........    do_stdout = strcmp("-", outfile) == 0;    if (flag.int_out->answer && flag.float_out->answer)	G_fatal_error(_("-i and -f are mutually exclusive"));    fd = Rast_open_old(name, "");    if (flag.int_out->answer)	is_fp = 0;    else if (flag.float_out->answer)	is_fp = 1;    else	is_fp = Rast_get_map_type(fd) != CELL_TYPE;    if (parm.bytes->answer)	bytes = atoi(parm.bytes->answer);    else if (is_fp)	bytes = 4;    else	bytes = 2;    if (is_fp && bytes < 4)	G_fatal_error(_("Floating-point output requires bytes=4 or bytes=8"));#ifndef HAVE_LONG_LONG_INT    if (!is_fp && bytes > 4)	G_fatal_error(_("Integer output doesn't support bytes=8 in this build"));#endif    G_get_window(&region);    /* open bin file for writing */    if (do_stdout)	fp = stdout;    else if (NULL == (fp = fopen(outfile, "w")))	G_fatal_error(_("Unable to create file <%s>"), outfile);    /* Set up Parameters for GMT header */    if (flag.gmt_hd->answer) {	if (!is_fp && bytes > 4)	    G_fatal_error(_("GMT grid doesn't support 64-bit integers"));	make_gmt_header(&header, name, outfile, &region, null_val);    }    /* Write out BIL support files compatible with Arc-View */    if (flag.bil_hd->answer) {	G_message(_("Creating BIL support files..."));	write_bil_hdr(outfile, &region,		      bytes, order, flag.gmt_hd->answer, null_val);	write_bil_wld(outfile, &region);    }    /* Write out GMT Header if required */    if (flag.gmt_hd->answer)	write_gmt_header(&header, swap_flag, fp);    nrows = Rast_window_rows();    ncols = Rast_window_cols();    in_buf = Rast_allocate_d_buf();    out_buf = G_malloc(ncols * bytes);    if (is_fp) {	G_message(_("Exporting raster as floating values (bytes=%d)"), bytes);	if (flag.gmt_hd->answer)	    G_message(_("Writing GMT float format ID=1"));    }    else {	G_message(_("Exporting raster as integer values (bytes=%d)"), bytes);	if (flag.gmt_hd->answer)	    G_message(_("Writing GMT integer format ID=2"));    }    G_verbose_message(_("Using the current region settings..."));    G_verbose_message(_("north=%f"), region.north);    G_verbose_message(_("south=%f"), region.south);    G_verbose_message(_("east=%f"), region.east);    G_verbose_message(_("west=%f"), region.west);    G_verbose_message(_("r=%d"), region.rows);    G_verbose_message(_("c=%d"), region.cols);    for (row = 0; row < nrows; row++) {	G_percent(row, nrows, 2);	Rast_get_d_row(fd, in_buf, row);	convert_row(out_buf, in_buf, ncols, is_fp, bytes, swap_flag, null_val);	if (fwrite(out_buf, bytes, ncols, fp) != ncols)	    G_fatal_error(_("Error writing data"));    }    G_percent(row, nrows, 2);	/* finish it off */    Rast_close(fd);    fclose(fp);    return EXIT_SUCCESS;}

int test_large_file_sparse_random(int depths, int rows, int cols, int tile_size){    int sum = 0;    int x, y, z, i;    DCELL value, random_value;    DCELL *random_value_vector = G_calloc(RAND_VALUE_VECTOR_SIZE, sizeof(DCELL));    G_message("Testing DCELL put function for large files filled with sparse random values");    RASTER3D_Region region;    RASTER3D_Map *map = NULL;    /* We need to set up a specific region for the new raster3d map.     * First we safe the default region. */    Rast3d_get_window(&region);    region.bottom = -365.5;    region.top = 365.5;    region.south = -90;    region.north = 90;    region.west = -180;    region.east = 180;    region.rows = rows;    region.cols = cols;    region.depths = depths;    Rast3d_adjust_region(&region);    G_message("Creating 3D raster map filled with sparse random values");    map = Rast3d_open_new_opt_tile_size("test_put_get_value_dcell_large_sparse_random",                                        RASTER3D_USE_CACHE_XY, &region, DCELL_TYPE, tile_size);    /* The window is the same as the map region ... of course */    Rast3d_set_window_map(map, &region);    srand(1);    /* We fill the random value vector */    for(i = 0; i < RAND_VALUE_VECTOR_SIZE; i++) {        /* Put the counter as cell value */        value = (DCELL)rand();        value /= RAND_MAX;        if(value <= 0.7)            value = 0.0;        else if(value <= 0.8)            value = 1.0;        else if(value <= 0.9)            value = 2.0;        else if(value <= 1.0)            value = 3.0;        else            value = 4.0;        random_value_vector[i] = value;    }    i = 0;    for(z = 0; z < region.depths; z++) {        G_percent(z, region.depths, 1);        for(y = 0; y < region.rows; y++) {            for(x = 0; x < region.cols; x++) {                /* Put the counter as cell value */                value = random_value_vector[i];                Rast3d_put_value(map, x, y, z, &value, DCELL_TYPE);                i++;                if(i == RAND_VALUE_VECTOR_SIZE)                    i = 0;            }        }    }    G_percent(1, 1, 1);    /* Write everything to the disk */    Rast3d_flush_all_tiles(map);    Rast3d_close(map);    G_message("Verifying 3D raster map filled with sparse random values");    map = Rast3d_open_cell_old("test_put_get_value_dcell_large_sparse_random",                               G_mapset(), &region, DCELL_TYPE, RASTER3D_USE_CACHE_XY);    i = 0;    for(z = 0; z < region.depths; z++) {        G_percent(z, region.depths, 1);        for(y = 0; y < region.rows; y++) {            for(x = 0; x < region.cols; x++) {                /* Check the counter as cell value */                Rast3d_get_value(map, x, y, z, &value, DCELL_TYPE);                if(fabs(value - random_value_vector[i]) > EPSILON) {                    G_message("At: z %i y %i x %i -- value %.14lf != %.14lf/n",                              z, y, x, value, random_value);                    sum++;                }                i++;                if(i == RAND_VALUE_VECTOR_SIZE)                    i = 0;            }        }//.........这里部分代码省略.........

//.........这里部分代码省略.........    kexo[5] = KEXO5;    kexo[6] = KEXO6;    kexo[7] = KEXO7;    kexo[8] = KEXO8;    kexo[9] = KEXO9;    /********************/    /********************/    for (; *ptr != NULL; ptr++) {	if (nfiles > MAXFILES)	    G_fatal_error(_("Too many input maps. Only %d allowed."),			  MAXFILES);	name = *ptr;	/* Allocate input buffer */	in_data_type[nfiles-1] = Rast_map_type(name, "");	/* For some strange reason infd[0] cannot be used later */	/* So nfiles is initialized with nfiles = 1 */	infd[nfiles] = Rast_open_old(name, "");	Rast_get_cellhd(name, "", &cellhd);	inrast[nfiles-1] = Rast_allocate_buf(in_data_type[nfiles-1]);	nfiles++;    }    nfiles--;    if (nfiles < MAXFILES) 	G_fatal_error(_("The input band number should be 15"));    /***************************************************/    /* Allocate output buffer, use input map data_type */    nrows = Rast_window_rows();    ncols = Rast_window_cols();    out_data_type = DCELL_TYPE;    for (i = 0; i < MAXFILES; i++)	outrast[i] = Rast_allocate_buf(out_data_type);    outfd[1] = Rast_open_new(result0, 1);    outfd[2] = Rast_open_new(result1, 1);    outfd[3] = Rast_open_new(result2, 1);    outfd[4] = Rast_open_new(result3, 1);    outfd[5] = Rast_open_new(result4, 1);    outfd[6] = Rast_open_new(result5, 1);    outfd[7] = Rast_open_new(result6, 1);    outfd[8] = Rast_open_new(result7, 1);    outfd[9] = Rast_open_new(result8, 1);    outfd[10] = Rast_open_new(result9, 1);    outfd[11] = Rast_open_new(result10, 1);    outfd[12] = Rast_open_new(result11, 1);    outfd[13] = Rast_open_new(result12, 1);    outfd[14] = Rast_open_new(result13, 1);    outfd[15] = Rast_open_new(result14, 1);    /* Process pixels */    DCELL dout[MAXFILES];    DCELL d[MAXFILES];    for (row = 0; row < nrows; row++) {	G_percent(row, nrows, 2);	/* read input map */	for (i = 1; i <= MAXFILES; i++)	    Rast_get_row(infd[i], inrast[i-1], row, in_data_type[i-1]);	/*process the data */	for (col = 0; col < ncols; col++) {	    for (i = 0; i < MAXFILES; i++) {		switch (in_data_type[i]) {		case CELL_TYPE:		    d[i] = (double)((CELL *) inrast[i])[col];		    break;		case FCELL_TYPE:		    d[i] = (double)((FCELL *) inrast[i])[col];		    break;		case DCELL_TYPE:		    d[i] = (double)((DCELL *) inrast[i])[col];		    break;		}		/* if radiance mode or Thermal band */		if (radiance || i >= 10) {		    dout[i] = gain[i] * (d[i] - 1.0);		}		/* if reflectance default mode and Not Thermal Band */		else {		    dout[i] = gain[i] * (d[i] - 1.0);		    dout[i] =			rad2ref_aster(dout[i], doy, sun_elevation, kexo[i]);		}		outrast[i][col] = dout[i];	    }	}	for (i = 1; i <= MAXFILES; i++)	    Rast_put_row(outfd[i], outrast[i-1], out_data_type);    }    for (i = 1; i <= MAXFILES; i++) {	G_free(inrast[i-1]);	Rast_close(infd[i]);	G_free(outrast[i-1]);	Rast_close(outfd[i]);    }    exit(EXIT_SUCCESS);}

//.........这里部分代码省略.........    infile = param.inputfile->answer;    if (NULL == G_find_raster3d(infile, ""))	Rast3d_fatal_error(_("3D raster map <%s> not found"), infile);    map =	Rast3d_open_cell_old(infile, G_find_raster3d(infile, ""), &region,			RASTER3D_TILE_SAME_AS_FILE, RASTER3D_USE_CACHE_DEFAULT);    if (map == NULL)	Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"), infile);    map_type = Rast3d_tile_type_map(map);    i = 0;    while (param.percentile->answers[i])	i++;     n_zones = zone_info.n_zones;    if (n_zones == 0)        n_zones = 1;    stats = create_univar_stat_struct(map_type, i);    for (i = 0; i < n_zones; i++) {	unsigned int j;	for (j = 0; j < stats[i].n_perc; j++) {	    sscanf(param.percentile->answers[j], "%lf", &(stats[i].perc[j]));	}    }    for (z = 0; z < depths; z++) {	/* From the bottom to the top */	if (!(param.shell_style->answer))	    G_percent(z, depths - 1, 10);	for (y = 0; y < rows; y++) {	    for (x = 0; x < cols; x++) {		zone = 0;		if (zone_info.n_zones) {		    if (zmap_type == FCELL_TYPE) {			Rast3d_get_value(zmap, x, y, z, &val_f, FCELL_TYPE);			if (Rast3d_is_null_value_num(&val_f, FCELL_TYPE))			    continue;			if (val_f < 0)			    zone = val_f - 0.5;			else			    zone = val_f + 0.5;		    }		    else if (zmap_type == DCELL_TYPE) {			Rast3d_get_value(zmap, x, y, z, &val_d, DCELL_TYPE);			if (Rast3d_is_null_value_num(&val_d, DCELL_TYPE))			    continue;			if (val_d < 0)			    zone = val_d - 0.5;			else			    zone = val_d + 0.5;		    }                    zone -= zone_info.min;                }		if (map_type == FCELL_TYPE) {		    Rast3d_get_value(map, x, y, z, &val_f, map_type);		    if (!Rast3d_is_null_value_num(&val_f, map_type)) {			if (param.extended->answer) {			    if (stats[zone].n >= stats[zone].n_alloc) {				size_t msize;				stats[zone].n_alloc += 1000;				msize = stats[zone].n_alloc * sizeof(FCELL);

int null_distance(const char *name1, const char *name2, int *zerro_row, int *zerro_col){    RASTER_MAP_TYPE maptype1, maptype2;    const char *mapset;    int mapd1, mapd2;    void *inrast1, *inrast2;    int nrows, ncols, row, col;    void *cell1, *cell2;    /* NOTE: no need to controll, if the map exists. it should be checked in edge.c */    mapset = G_find_raster2(name1, "");    maptype1 = Rast_map_type(name1, mapset);    mapd1 = Rast_open_old(name1, mapset);    inrast1 = Rast_allocate_buf(maptype1);    mapset = G_find_raster2(name2, "");    maptype2 = Rast_map_type(name2, mapset);    mapd2 = Rast_open_old(name2, mapset);    inrast2 = Rast_allocate_buf(maptype2);    G_message(_("Reading maps  <%s,%s> while finding 0 distance ..."), name1,              name2);    ncols = Rast_window_cols();    nrows = Rast_window_rows();    for (row = 0; row < nrows; row++) {        G_percent(row, nrows, 2);        Rast_get_row(mapd1, inrast1, row, maptype1);        Rast_get_row(mapd2, inrast2, row, maptype2);        for (col = 0; col < ncols; col++) {            /* first raster */            switch (maptype1) {            case CELL_TYPE:                cell1 = ((CELL **) inrast1)[col];                break;            case FCELL_TYPE:                cell1 = ((FCELL **) inrast1)[col];                break;            case DCELL_TYPE:                cell1 = ((DCELL **) inrast1)[col];                break;            }            /* second raster */            switch (maptype2) {            case CELL_TYPE:                cell2 = ((CELL **) inrast2)[col];                break;            case FCELL_TYPE:                cell2 = ((FCELL **) inrast2)[col];                break;            case DCELL_TYPE:                cell2 = ((DCELL **) inrast2)[col];                break;            }            if (!Rast_is_null_value(&cell1, maptype1) &&                    !Rast_is_null_value(&cell2, maptype2)) {                *zerro_row = row;                *zerro_col = col;                /* memory cleanup */                G_free(inrast1);                G_free(inrast2);                /* closing raster maps */                Rast_close(mapd1);                Rast_close(mapd2);                return 1;            }        }    }    /* memory cleanup */    G_free(inrast1);    G_free(inrast2);    /* closing raster maps */    Rast_close(mapd1);    Rast_close(mapd2);    return 0;}

//.........这里部分代码省略.........    ymax = cellhd.north;    nrows = Rast_window_rows();    ncols = Rast_window_cols();    /***************************************************/    /* Allocate output buffer */    /***************************************************/    outrast = Rast_allocate_d_buf();    outfd = Rast_open_new(h0, DCELL_TYPE);    /***************************************************/    /* Allocate memory for temporary images            */    double **d_Roh, **d_Rah;    if ((d_Roh = G_alloc_matrix(nrows, ncols)) == NULL)	G_message("Unable to allocate memory for temporary d_Roh image");    if ((d_Rah = G_alloc_matrix(nrows, ncols)) == NULL)	G_message("Unable to allocate memory for temporary d_Rah image");    /***************************************************/    /* MANUAL T0DEM WET/DRY PIXELS */    DCELL d_Rn_dry,d_g0_dry;    DCELL d_t0dem_dry,d_t0dem_wet;    if (flag2->answer) {	/* Process tempk min / max pixels */	/* Internal use only */	DCELL d_Rn_wet,d_g0_wet;	DCELL d_Rn,d_g0,d_h0;	DCELL t0dem_min,t0dem_max;        /*********************/	for (row = 0; row < nrows; row++) {	    DCELL d_t0dem;	    G_percent(row, nrows, 2);	    Rast_get_d_row(infd_t0dem,inrast_t0dem,row);	    Rast_get_d_row(infd_Rn,inrast_Rn,row);	    Rast_get_d_row(infd_g0,inrast_g0,row);	    /*process the data */	    for (col = 0; col < ncols; col++) {		d_t0dem = ((DCELL *) inrast_t0dem)[col];		d_Rn = ((DCELL *) inrast_Rn)[col];		d_g0 = ((DCELL *) inrast_g0)[col];		if (Rast_is_d_null_value(&d_t0dem) ||		    Rast_is_d_null_value(&d_Rn) ||                     Rast_is_d_null_value(&d_g0)) {		    /* do nothing */		}		else {		    if (d_t0dem <= 250.0) {			/* do nothing */		    }		    else {			d_h0 = d_Rn - d_g0;			if (d_t0dem < t0dem_min &&			    d_Rn > 0.0 && d_g0 > 0.0 && d_h0 > 0.0 &&			    d_h0 < 100.0) {			    t0dem_min = d_t0dem;			    d_t0dem_wet = d_t0dem;			    d_Rn_wet = d_Rn;			    d_g0_wet = d_g0;			    m_col_wet = col;			    m_row_wet = row;			}			if (d_t0dem > t0dem_max &&			    d_Rn > 0.0 && d_g0 > 0.0 && d_h0 > 100.0 &&			    d_h0 < 500.0) {

//.........这里部分代码省略.........	ret = get_node_costs(terms[0], terms[i], &cost);	if (ret == 0) {            /* GTC Terminal refers to an Steiner tree endpoint */	    G_fatal_error(_("Terminal at node [%d] cannot be connected "			    "to terminal at node [%d]"), terms[0], terms[i]);	}    }    /* Remove not reachable from list of SP candidates */    j = 0;    for (i = 1; i <= nnodes; i++) {	ret = get_node_costs(terms[0], i, &cost);	if (ret == 0) {	    testnode[i] = 0;	    G_debug(2, "node %d removed from list of Steiner point candidates/n", i );	    j++;	}    }    G_message(_("[%d] (not reachable) nodes removed from list "		"of Steiner point candidates"), j);    /* calc costs for terminals MST */    ret = mst(&Map, terms, nterms, &cost, PORT_DOUBLE_MAX, NULL, NULL, 0, 1);	/* no StP, rebuild */    G_message(_("MST costs = %f"), cost);    /* Go through all nodes and try to use as steiner points -> find that which saves most costs */    nspused = 0;    for (j = 0; j < nsp; j++) {	sp = 0;	G_verbose_message(_("Search for [%d]. Steiner point"), j + 1);	for (i = 1; i <= nnodes; i++) {	    G_percent(i, nnodes, 1);	    if (testnode[i] == 0) {		G_debug(3, "skip test for %d/n", i);		continue;	    }	    ret =		mst(&Map, terms, nterms + j, &tmpcost, cost, NULL, NULL, i,		    0);	    G_debug(2, "cost = %f x %f/n", tmpcost, cost);	    if (tmpcost < cost) {	/* sp candidate */		G_debug(3,			"  steiner candidate node = %d mst = %f (x last = %f)/n",			i, tmpcost, cost);		sp = i;		cost = tmpcost;	    }	}	if (sp > 0) {	    G_message(_("Steiner point at node [%d] was added "			"to terminals (MST costs = %f)"), sp, cost);	    terms[nterms + j] = sp;	    init_node_costs(&Map, sp);	    testnode[sp] = 0;	    nspused++;	    /* rebuild for nex cycle */	    ret =		mst(&Map, terms, nterms + nspused, &tmpcost, PORT_DOUBLE_MAX,		    NULL, NULL, 0, 1);	}	else {			/* no steiner found */	    G_message(_("No Steiner point found -> leaving cycle"));	    break;	}

void CalcSurface(void){    int Count, OutRows, OutCols;    int Row, Row2, Col, Col2, RanRows, RanCols;    int owC, oeC, onR, osR, wC, eC, nR, sR;    double **Randoms;    G_debug(2, "CalcSurface()");    OutRows = BigF.RowPlus;    OutCols = BigF.ColPlus;    RanRows = (2 * OutRows + Rs);    RanCols = (2 * OutCols + Cs);    owC = osR = 0;    wC = OutCols;    sR = OutRows;    oeC = RanCols - 1;    onR = RanRows - 1;    if (OutCols > 0)	eC = RanCols - (OutCols + 1);    else	eC = oeC;    if (OutRows > 0)	nR = RanRows - (OutRows + 1);    else	nR = onR;    Randoms = (double **)G_malloc(sizeof(double *) * RanRows);    for (Row = 0; Row < RanRows; Row++)	Randoms[Row] = (double *)G_malloc(RanCols * sizeof(double));    /* OLD       for( Row = OutRows; Row < RanRows; Row++) {       for( Col = OutCols; Col < OutCols + Cs; Col++) {       Randoms[ Row][ Col] = GasDev();       }       }       for( Row = OutRows - 1; Row >= 0; Row--) {       for( Col = OutCols; Col < OutCols + Cs; Col++) {       Randoms[ Row][ Col] = GasDev();       }       }       for( Row = 0; Row < RanRows; Row++) {       for( Col = 0; Col < OutCols; Col++)       Randoms[ Row][ Col] = GasDev();       for( Col = OutCols + Cs; Col < RanCols; Col++)       Randoms[ Row][ Col] = GasDev();       }       end OLD */    for (Row = sR; Row <= nR; Row++) {	for (Col = wC; Col <= eC; Col++) {	    Randoms[Row][Col] = GasDev();	}    }    Col = wC - 1;    Col2 = eC + 1;    while (Col >= 0) {	for (Row = sR; Row <= nR; Row++) {	    Randoms[Row][Col] = GasDev();	    Randoms[Row][Col2] = GasDev();	}	Col--;	Col2++;    }    Row = sR - 1;    Row2 = nR + 1;    while (Row >= 0) {	for (Col = 0; Col < RanCols; Col++) {	    Randoms[Row][Col] = GasDev();	    Randoms[Row2][Col] = GasDev();	}	Row--;	Row2++;    }    Count = 0;    if (FDM == -1) {	for (Row = 0; Row < Rs; Row++) {	    if (ODD(Row)) {		for (Col = Cs - 1; Col >= 0; Col--) {		    G_percent(Count++, MapCount, 1);		    Surface[Row][Col] =			MakePP(Row, Col, OutRows, OutCols, Randoms, BigF);		}	    }	    else {		for (Col = 0; Col < Cs; Col++) {		    G_percent(Count++, MapCount, 1);		    Surface[Row][Col] =			MakePP(Row, Col, OutRows, OutCols, Randoms, BigF);		}	    }	}//.........这里部分代码省略.........

//.........这里部分代码省略.........	for (i = 0; i < rows; i++) {	    double *data0, *data1;	    data0 = data[0] + i * cols;	    data1 = data[1] + i * cols;	    G_get_map_row(maskfd, maskbuf, i);	    for (j = 0; j < cols; j++, data0++, data1++) {		if (maskbuf[j] == (CELL) 0) {		    *(data0) = 0.0;		    *(data1) = 0.0;		}	    }	}    }    G_message(_("Rotating data..."));    /* rotate the data array for standard display */    for (i = 0; i < rows; i++) {	double temp;	for (j = 0; j < halfcols; j++) {	    temp = *(data[0] + i * cols + j);	    *(data[0] + i * cols + j) = *(data[0] + i * cols + j + halfcols);	    *(data[0] + i * cols + j + halfcols) = temp;	    temp = *(data[1] + i * cols + j);	    *(data[1] + i * cols + j) = *(data[1] + i * cols + j + halfcols);	    *(data[1] + i * cols + j + halfcols) = temp;	}    }    for (i = 0; i < halfrows; i++) {	double temp;	for (j = 0; j < cols; j++) {	    temp = *(data[0] + i * cols + j);	    *(data[0] + i * cols + j) =		*(data[0] + (i + halfrows) * cols + j);	    *(data[0] + (i + halfrows) * cols + j) = temp;	    temp = *(data[1] + i * cols + j);	    *(data[1] + i * cols + j) =		*(data[1] + (i + halfrows) * cols + j);	    *(data[1] + (i + halfrows) * cols + j) = temp;	}    }    /* close input cell maps and release the row buffers */    fclose(realfp);    fclose(imagfp);    if (maskfd >= 0) {	G_close_cell(maskfd);	G_free(maskbuf);    }    /* perform inverse FFT */    G_message(_("Starting Inverse FFT..."));    fft(1, data, totsize, cols, rows);    /* set up a window for the transform cell map */    G_set_window(&orig_wind);    /* open the output cell map and allocate a cell row buffer */    if ((outputfd = G_open_cell_new(Cellmap_orig)) < 0)	G_fatal_error(_("Unable to create raster map <%s>"),		      Cellmap_orig);    cell_row = G_allocate_cell_buf();    /* Write out result to a new cell map */    G_message(_("Writing data..."));    for (i = 0; i < or; i++) {	for (j = 0; j < oc; j++) {	    *(cell_row + j) = (CELL) (*(data[0] + i * cols + j) + 0.5);	}	G_put_raster_row(outputfd, cell_row, CELL_TYPE);	G_percent(i+1, or, 2);    }    G_close_cell(outputfd);    G_free(cell_row);    {	struct Colors colors;	struct Range range;	CELL min, max;	/* make a real component color table */	G_read_range(Cellmap_orig, G_mapset(), &range);	G_get_range_min_max(&range, &min, &max);	G_make_grey_scale_colors(&colors, min, max);	G_write_colors(Cellmap_orig, G_mapset(), &colors);    }    /* Release memory resources */    G_free(data[0]);    G_free(data[1]);    G_done_msg(" ");    exit(EXIT_SUCCESS);}

//.........这里部分代码省略.........	exit(EXIT_FAILURE);        rnet = input1->answer;    g0 = input2->answer;    h0 = input3->answer;    result1 = output1->answer;    result2 = output2->answer;    makin = flag1->answer;        /***************************************************/     infd_rnet = Rast_open_old(rnet, "");    inrast_rnet = Rast_allocate_d_buf();        /***************************************************/     infd_g0 = Rast_open_old(g0, "");    inrast_g0 = Rast_allocate_d_buf();        /***************************************************/     infd_h0 = Rast_open_old(h0, "");    inrast_h0 = Rast_allocate_d_buf();        /***************************************************/     nrows = Rast_window_rows();    ncols = Rast_window_cols();    outrast1 = Rast_allocate_d_buf();    if (makin) 	outrast2 = Rast_allocate_d_buf();        /* Create New raster files */     outfd1 = Rast_open_new(result1, DCELL_TYPE);    if (makin) 	outfd2 = Rast_open_new(result2, DCELL_TYPE);            /* Process pixels */     for (row = 0; row < nrows; row++)    {	DCELL d;	DCELL d_rnet;	DCELL d_g0;	DCELL d_h0;	G_percent(row, nrows, 2);	        /* read input maps */         Rast_get_d_row(infd_rnet, inrast_rnet, row);	Rast_get_d_row(infd_g0, inrast_g0, row);	Rast_get_d_row(infd_h0, inrast_h0, row);	        /*process the data */         for (col = 0; col < ncols; col++)        {            d_rnet = ((DCELL *) inrast_rnet)[col];            d_g0 = ((DCELL *) inrast_g0)[col];            d_h0 = ((DCELL *) inrast_h0)[col];	    if (Rast_is_d_null_value(&d_rnet) ||                 Rast_is_d_null_value(&d_g0) ||		Rast_is_d_null_value(&d_h0)) {		Rast_set_d_null_value(&outrast1[col], 1);		if (makin) 		    Rast_set_d_null_value(&outrast2[col], 1);	    }	    else {                /* calculate evaporative fraction       */                 d = evap_fr(d_rnet, d_g0, d_h0);		outrast1[col] = d;                /* calculate soil moisture              */                 if (makin)                 {                    d = soilmoisture(d);                    outrast2[col] = d;		}	    }        }	Rast_put_d_row(outfd1, outrast1);	if (makin)             Rast_put_d_row(outfd2, outrast2);    }    G_free(inrast_rnet);    G_free(inrast_g0);    G_free(inrast_h0);    Rast_close(infd_rnet);    Rast_close(infd_g0);    Rast_close(infd_h0);    G_free(outrast1);    Rast_close(outfd1);    if (makin) {        G_free(outrast2);	Rast_close(outfd2);    }    Rast_short_history(result1, "raster", &history);    Rast_command_history(&history);    Rast_write_history(result1, &history);    if (makin) {	Rast_short_history(result2, "raster", &history);	Rast_command_history(&history);	Rast_write_history(result2, &history);    }    exit(EXIT_SUCCESS);}

int main(int argc, char *argv[]){    int i, j, centroid, otype, count;    int nlines, nareas;    int field;    struct GModule *module;    struct Option *in_opt, *field_opt, *out_opt, *type_opt;    struct Option *size_opt, *zmod_opt, *objmod_opt;    FILE *fd;    /* Vector */    struct Map_info In;    struct line_pnts *Points;    struct line_cats *Cats;    int type;    G_gisinit(argv[0]);    /* Module options */    module = G_define_module();    G_add_keyword(_("vector"));    G_add_keyword(_("export"));    module->description =	_("Converts GRASS x,y,z points to POV-Ray x,z,y format.");    in_opt = G_define_standard_option(G_OPT_V_INPUT);    field_opt = G_define_standard_option(G_OPT_V_FIELD_ALL);        type_opt = G_define_standard_option(G_OPT_V3_TYPE);    type_opt->answer = "point,line,area,face";    out_opt = G_define_standard_option(G_OPT_F_OUTPUT);    out_opt->required = YES;    out_opt->description = _("Name for output POV file");    size_opt = G_define_option();    size_opt->key = "size";    size_opt->type = TYPE_STRING;    size_opt->required = NO;    size_opt->answer = "10";    size_opt->label = _("Radius of sphere for points and tube for lines");    size_opt->description = _("May be also variable, e.g. grass_r.");    zmod_opt = G_define_option();    zmod_opt->key = "zmod";    zmod_opt->type = TYPE_STRING;    zmod_opt->required = NO;    zmod_opt->description = _("Modifier for z coordinates");    zmod_opt->description = _("This string is appended to each z coordinate. "			      "Examples: '*10', '+1000', '*10+100', '*exaggeration'");    objmod_opt = G_define_option();    objmod_opt->key = "objmod";    objmod_opt->type = TYPE_STRING;    objmod_opt->required = NO;    objmod_opt->label = _("Object modifier (OBJECT_MODIFIER in POV-Ray documentation)");    objmod_opt->description = _("Example: /"pigment { color red 0 green 1 blue 0 }/"");        if (G_parser(argc, argv))	exit(EXIT_FAILURE);    /* Check output type */    otype = Vect_option_to_types(type_opt);    Points = Vect_new_line_struct();    Cats = Vect_new_cats_struct();    /* open input vector */    Vect_set_open_level(2);    if (Vect_open_old2(&In, in_opt->answer, "", field_opt->answer) < 0)	G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);        field = Vect_get_field_number(&In, field_opt->answer);        /* Open output file */    if ((fd = fopen(out_opt->answer, "w")) == NULL) {	Vect_close(&In);	G_fatal_error(_("Unable to create output file <%s>"), out_opt->answer);    }    if (zmod_opt->answer == NULL)	    zmod_opt->answer = G_store("");    if (objmod_opt->answer == NULL)	    objmod_opt->answer = G_store("");    nlines = Vect_get_num_lines(&In);    nareas = Vect_get_num_areas(&In);    count = 0;    /* Lines */    if ((otype &	 (GV_POINTS | GV_LINES | GV_BOUNDARY | GV_CENTROID | GV_FACE |	  GV_KERNEL))) {	for (i = 1; i <= nlines; i++) {	    G_percent(i, nlines, 2);	    type = Vect_read_line(&In, Points, Cats, i);	    G_debug(2, "line = %d type = %d", i, type);	    	    if (field != -1 && Vect_cat_get(Cats, field, NULL) == 0)		continue;//.........这里部分代码省略.........

/*    Called by Vect_remove_bridges() and Vect_chtype_bridges():   chtype = 0 -> works like Vect_remove_bridges()   chtype = 1 -> works like Vect_chtype_bridges()   Algorithm: Go thorough all lines,    if both sides of the line have left and side 0 (candidate) do this check:   follow adjacent lines in one direction (nearest to the right at the end node),   if we reach this line again without dangle in the way, but with this line    traversed from other side it is a bridge.   List of all lines in chain is created during the cycle. */voidremove_bridges(struct Map_info *Map, int chtype, struct Map_info *Err){    int i, type, nlines, line;    int left, right, node1, node2, current_line, next_line;    int bridges_removed = 0;	/* number of removed bridges */    int lines_removed = 0;	/* number of lines removed */    char *lmsg;    struct Plus_head *Plus;    struct line_pnts *Points;    struct line_cats *Cats;    struct ilist *CycleList;    struct ilist *BridgeList;    int dangle, other_side;    if (chtype)	lmsg = "changed lines";    else	lmsg = "removed lines";    Plus = &(Map->plus);    Points = Vect_new_line_struct();    Cats = Vect_new_cats_struct();    CycleList = Vect_new_list();    BridgeList = Vect_new_list();    nlines = Vect_get_num_lines(Map);    G_debug(1, "nlines =  %d", nlines);    for (line = 1; line <= nlines; line++) {	G_percent(line, nlines, 1);	if (!Vect_line_alive(Map, line))	    continue;	type = Vect_read_line(Map, NULL, NULL, line);	if (!(type & GV_BOUNDARY))	    continue;	Vect_get_line_areas(Map, line, &left, &right);	if (left != 0 || right != 0)	    continue;		/* Cannot be bridge */	G_debug(2, "line %d - bridge candidate", line);	Vect_get_line_nodes(Map, line, &node1, &node2);	if (abs(node1) == abs(node2))	    continue;		/* either zero length or loop -> cannot be a bridge */	current_line = -line;	/* we start with negative (go forward, node2 ) */	dangle = 0;	other_side = 0;	Vect_reset_list(CycleList);	Vect_reset_list(BridgeList);	while (1) {	    next_line =		dig_angle_next_line(Plus, current_line, GV_RIGHT,				    GV_BOUNDARY);	    /* Add this line to the list */	    /* TODO: Vect_val_in_list() and Vect_list_append() behave O(n)	     * change to O(log n) */	    if (Vect_val_in_list(CycleList, abs(next_line)))	/* other side -> bridge chain */		Vect_list_append(BridgeList, abs(next_line));	    else		dig_list_add(CycleList, abs(next_line)); /* not in list, can add new line fast */	    if (abs(next_line) == abs(current_line)) {		G_debug(4, "  dangle -> no bridge");		dangle = 1;		break;	    }	    if (abs(next_line) == line) {	/* start line reached */		/* which side */		if (next_line < 0) {	/* other side (connected by node 2) */		    G_debug(5, "  other side reached");		    other_side = 1;		}		else {		/* start side */		    break;		}	    }//.........这里部分代码省略.........

int test_large_file_zeros(int depths, int rows, int cols, int tile_size){    int sum = 0;    int x, y, z;    DCELL value;    G_message("Testing DCELL put function for large files filled with zeros");    RASTER3D_Region region;    RASTER3D_Map *map = NULL;    /* We need to set up a specific region for the new raster3d map.     * First we safe the default region. */    Rast3d_get_window(&region);    region.bottom = -365.5;    region.top = 365.5;    region.south = -90;    region.north = 90;    region.west = -180;    region.east = 180;    region.rows = rows;    region.cols = cols;    region.depths = depths;    Rast3d_adjust_region(&region);    G_message("Creating 3D raster map filled with zeros");    map = Rast3d_open_new_opt_tile_size("test_put_get_value_dcell_large_zeros",                                        RASTER3D_USE_CACHE_XY, &region, DCELL_TYPE, tile_size);    /* The window is the same as the map region ... of course */    Rast3d_set_window_map(map, &region);    for(z = 0; z < region.depths; z++) {        G_percent(z, region.depths, 1);        for(y = 0; y < region.rows; y++) {            for(x = 0; x < region.cols; x++) {                /* Put the counter as cell value */                value = 0.0;                Rast3d_put_value(map, x, y, z, &value, DCELL_TYPE);            }        }    }    G_percent(1, 1, 1);    /* Write everything to the disk */    Rast3d_flush_all_tiles(map);    Rast3d_close(map);    G_message("Verifying 3D raster map filled with zeros");    map = Rast3d_open_cell_old("test_put_get_value_dcell_large_zeros",                               G_mapset(), &region, DCELL_TYPE, RASTER3D_USE_CACHE_XY);    for(z = 0; z < region.depths; z++) {        G_percent(z, region.depths, 1);        for(y = 0; y < region.rows; y++) {            for(x = 0; x < region.cols; x++) {                /* Check the counter as cell value */                Rast3d_get_value(map, x, y, z, &value, DCELL_TYPE);                if(value > EPSILON) {                    G_message("At: z %i y %i x %i -- value %.14lf != 0.0/n",                              z, y, x, value);                    sum++;                }            }        }    }    G_percent(1, 1, 1);    Rast3d_close(map);    G_remove("grid3", "test_put_get_value_dcell_large_zeros");    return sum;}

//.........这里部分代码省略.........	fd = fopen(input, "r");    if (fd == NULL) {	G_fatal_error(_("Unable to read input from <%s>"), input);    }    direction = 1;    sz = 0;    if (flag.s->answer) {	sz = getgrdhead(fd, &cellhd);	/* for Surfer files, the data type is always FCELL_TYPE,	   the multiplier and the null_val_str are never used */	data_type = FCELL_TYPE;	mult = 1.;	null_val_str = "";	/* rows in surfer files are ordered from bottom to top,	   opposite of normal GRASS ordering */	direction = -1;    }    else	sz = gethead(fd, &cellhd, &data_type, &mult, &null_val_str);    if (!sz)	G_fatal_error(_("Can't get cell header"));    nrows = cellhd.rows;    ncols = cellhd.cols;    Rast_set_window(&cellhd);    if (nrows != Rast_window_rows())	G_fatal_error(_("OOPS: rows changed from %d to %d"), nrows,		      Rast_window_rows());    if (ncols != Rast_window_cols())	G_fatal_error(_("OOPS: cols changed from %d to %d"), ncols,		      Rast_window_cols());    rast_ptr = Rast_allocate_buf(data_type);    rast = rast_ptr;    cf = Rast_open_new(output, data_type);    for (row = 0; row < nrows; row++) {	G_percent(row, nrows, 2);	for (col = 0; col < ncols; col++) {	    if (fscanf(fd, "%s", y) != 1) {		Rast_unopen(cf);		G_fatal_error(_("Data conversion failed at row %d, col %d"),			      row + 1, col + 1);	    }	    if (strcmp(y, null_val_str)) {		x = atof(y);		if ((float)x == GS_BLANK) {		    Rast_set_null_value(rast_ptr, 1, data_type);		}		else {		    Rast_set_d_value(rast_ptr,					 (DCELL) (x * mult), data_type);		}	    }	    else {		Rast_set_null_value(rast_ptr, 1, data_type);	    }	    rast_ptr = G_incr_void_ptr(rast_ptr, Rast_cell_size(data_type));	}	fwrite(rast, Rast_cell_size(data_type), ncols, ft);	rast_ptr = rast;    }    G_percent(nrows, nrows, 2);    G_debug(1, "Creating support files for %s", output);    sz = 0;    if (direction < 0) {	sz = -ncols * Rast_cell_size(data_type);	G_fseek(ft, sz, SEEK_END);	sz *= 2;    }    else {	G_fseek(ft, 0L, SEEK_SET);    }    for (row = 0; row < nrows; row += 1) {	fread(rast, Rast_cell_size(data_type), ncols, ft);	Rast_put_row(cf, rast, data_type);	G_fseek(ft, sz, SEEK_CUR);    }    fclose(ft);    unlink(temp);    Rast_close(cf);    if (title)	Rast_put_cell_title(output, title);    Rast_short_history(output, "raster", &history);    Rast_command_history(&history);    Rast_write_history(output, &history);    G_done_msg(" ");    exit(EXIT_SUCCESS);}

int test_large_file(int depths, int rows, int cols, int tile_size){    int sum = 0;    int x, y, z, count;    DCELL value;    G_message("Testing DCELL put function for large files");    RASTER3D_Region region;    RASTER3D_Map *map = NULL;    /* We need to set up a specific region for the new raster3d map.     * First we safe the default region. */    Rast3d_get_window(&region);    region.bottom = -365.5;    region.top = 365.5;    region.south = -90;    region.north = 90;    region.west = -180;    region.east = 180;    region.rows = rows;    region.cols = cols;    region.depths = depths;    Rast3d_adjust_region(&region);    G_message("Creating 3D raster map");    map = Rast3d_open_new_opt_tile_size("test_put_get_value_dcell_large",                                        RASTER3D_USE_CACHE_XY, &region, DCELL_TYPE, tile_size);    /* The window is the same as the map region ... of course */    Rast3d_set_window_map(map, &region);    /*Write -1 first to see if the tile handling works correctly */    for(z = 0; z < region.depths; z++) {        G_percent(z, region.depths, 1);        for(y = 0; y < region.rows; y++) {            for(x = 0; x < region.cols; x++) {                /* Put the counter as cell value */                value = -1;                Rast3d_put_value(map, x, y, z, &value, DCELL_TYPE);            }        }    }    G_percent(1, 1, 1);    G_message("Rewriting the values");    /* Now write the values to be evaluated */    count = 1;    for(z = 0; z < region.depths; z++) {        G_percent(z, region.depths, 1);        for(y = 0; y < region.rows; y++) {            for(x = 0; x < region.cols; x++) {                /* Put the counter as cell value */                value = count;                Rast3d_put_value(map, x, y, z, &value, DCELL_TYPE);                count++;            }        }    }    G_percent(1, 1, 1);    G_message("Verifying 3D raster map");    count = 1;    for(z = 0; z < region.depths; z++) {        G_percent(z, region.depths, 1);        for(y = 0; y < region.rows; y++) {            for(x = 0; x < region.cols; x++) {                /* Check the counter as cell value */                Rast3d_get_value(map, x, y, z, &value, DCELL_TYPE);                if(fabs(value - (double)(count) > EPSILON)) {                    G_message("At: z %i y %i x %i -- value %.14lf != %.14lf/n",                              z, y, x, value, (double)(count));                    sum++;                }                count++;            }        }    }    G_percent(1, 1, 1);    Rast3d_close(map);    G_remove("grid3", "test_put_get_value_dcell_large");    return sum;}

//.........这里部分代码省略.........    /* segment lib segment size */    srows = 64;    scols = 64;    nseg = manage_memory(srows, scols, globals);    /* create segment structures */    if (Segment_open	(&globals->bands_seg, G_tempfile(), globals->nrows, globals->ncols, srows,	 scols, inlen, nseg) != 1)	G_fatal_error("Unable to create input temporary files");    if (Segment_open	(&globals->rid_seg, G_tempfile(), globals->nrows, globals->ncols, srows,	 scols, outlen, nseg * 2) != 1)	G_fatal_error("Unable to create input temporary files");    /* load input bands to segment structure */    if (Ref.nfiles > 1)	G_message(_("Loading input bands..."));    else	G_message(_("Loading input band..."));    globals->bands_val = (double *)G_malloc(inlen);    globals->second_val = (double *)G_malloc(inlen);    /* initial segment ID */    s = 1;    globals->row_min = globals->nrows;    globals->row_max = 0;    globals->col_min = globals->ncols;    globals->col_max = 0;    for (row = 0; row < globals->nrows; row++) {	G_percent(row, globals->nrows, 4);	for (n = 0; n < Ref.nfiles; n++) {	    Rast_get_d_row(in_fd[n], inbuf[n], row);	}	for (col = 0; col < globals->ncols; col++) {	    is_null = 0;	/*Assume there is data */	    for (n = 0; n < Ref.nfiles; n++) {		globals->bands_val[n] = inbuf[n][col];		if (Rast_is_d_null_value(&inbuf[n][col])) {		    is_null = 1;		}		else {		    if (globals->weighted == FALSE)		    	/* scaled version */			globals->bands_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]);		}	    }	    if (Segment_put(&globals->bands_seg,	                    (void *)globals->bands_val, row, col) != 1)		G_fatal_error(_("Unable to write to temporary file"));	    if (!is_null) {		if (!globals->seeds) {		    /* sequentially number all cells with a unique segment ID */		    id = s;		    s++;		}		/* get min/max row/col to narrow the processing window */		if (globals->row_min > row)		    globals->row_min = row;		if (globals->row_max < row)

//.........这里部分代码省略.........        data_type_bluechan = Rast_map_type(bluechan, "");        inrast_bluechan = Rast_allocate_buf(data_type_bluechan);    }    if (chan5chan) {        infd_chan5chan = Rast_open_old(chan5chan, "");        data_type_chan5chan = Rast_map_type(chan5chan, "");        inrast_chan5chan = Rast_allocate_buf(data_type_chan5chan);    }    if (chan7chan) {        infd_chan7chan = Rast_open_old(chan7chan, "");        data_type_chan7chan = Rast_map_type(chan7chan, "");        inrast_chan7chan = Rast_allocate_buf(data_type_chan7chan);    }    nrows = Rast_window_rows();    ncols = Rast_window_cols();    /* Create New raster files */    outfd = Rast_open_new(result, DCELL_TYPE);    outrast = Rast_allocate_d_buf();    /* Process pixels */    for (row = 0; row < nrows; row++)    {	DCELL d_bluechan;	DCELL d_greenchan;	DCELL d_redchan;	DCELL d_nirchan;	DCELL d_chan5chan;	DCELL d_chan7chan;	G_percent(row, nrows, 2);	/* read input maps */	Rast_get_row(infd_redchan,inrast_redchan,row,data_type_redchan);	if (nirchan) {	    Rast_get_row(infd_nirchan,inrast_nirchan,row,data_type_nirchan);	}	if (bluechan) {	    Rast_get_row(infd_bluechan,inrast_bluechan,row,data_type_bluechan);	}	if (greenchan) {	    Rast_get_row(infd_greenchan,inrast_greenchan,row,data_type_greenchan);	}	if (chan5chan) {	    Rast_get_row(infd_chan5chan,inrast_chan5chan,row,data_type_chan5chan);	}	if (chan7chan) {	    Rast_get_row(infd_chan7chan,inrast_chan7chan,row,data_type_chan7chan);	}	/* process the data */	for (col = 0; col < ncols; col++)	{	    switch(data_type_redchan){		    case CELL_TYPE:			d_redchan   = (double) ((CELL *) inrast_redchan)[col];			if(opt.bits->answer)				d_redchan *= 1.0/(pow(2,dnbits)-1);				break;		    case FCELL_TYPE:			d_redchan   = (double) ((FCELL *) inrast_redchan)[col];			break;		    case DCELL_TYPE:			d_redchan   = ((DCELL *) inrast_redchan)[col];

//.........这里部分代码省略.........	Vect_reset_line(Points);	Vect_append_point(Points, coor[i].x, coor[i].y, 0.0);	Vect_append_point(Points, coor[i - 1].x, coor[i - 1].y, 0.0);	Vect_write_line(&Out, Type, Points, Cats);    }    G_free(coor);    /* Copy input points as centroids */    nfields = Vect_cidx_get_num_fields(&In);    cats = (int **)G_malloc(nfields * sizeof(int *));    ncats = (int *)G_malloc(nfields * sizeof(int));    fields = (int *)G_malloc(nfields * sizeof(int));    for (i = 0; i < nfields; i++) {	ncats[i] = 0;	cats[i] =	    (int *)G_malloc(Vect_cidx_get_num_cats_by_index(&In, i) *			    sizeof(int));	fields[i] = Vect_cidx_get_field_number(&In, i);    }    if (line_flag->answer)	ctype = GV_POINT;    else	ctype = GV_CENTROID;    nlines = Vect_get_num_lines(&In);    G_message(_("Writing sites to output..."));    for (line = 1; line <= nlines; line++) {	G_percent(line, nlines, 2);	type = Vect_read_line(&In, Points, Cats, line);	if (!(type & GV_POINTS))	    continue;	if (!Vect_point_in_box(Points->x[0], Points->y[0], 0.0, &Box))	    continue;	Vect_write_line(&Out, ctype, Points, Cats);	for (i = 0; i < Cats->n_cats; i++) {	    int f, j;	    f = -1;	    for (j = 0; j < nfields; j++) {	/* find field */		if (fields[j] == Cats->field[i]) {		    f = j;		    break;		}	    }	    if (f > -1) {		cats[f][ncats[f]] = Cats->cat[i];		ncats[f]++;	    }	}    }    /* Copy tables */    if (!(table_flag->answer)) {	int ttype, ntabs = 0;	struct field_info *IFi, *OFi;

int do_cum_mfd(void){    int r, c, dr, dc;    DCELL value, valued, *wat_nbr;    POINT point;    WAT_ALT wa;    int killer, threshold;    /* MFD */    int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;    double *dist_to_nbr, *weight, sum_weight, max_weight;    int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side, max_side;    double dx, dy;    CELL ele, *ele_nbr;    char asp_val, asp_val_down;    double prop, max_acc;    int workedon, edge, is_swale, flat;    char *flag_nbr, this_flag_value, flag_value;    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };    G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));    G_debug(1, "MFD convergence factor set to %d.", c_fac);    /* distances to neighbours */    dist_to_nbr = (double *)G_malloc(sides * sizeof(double));    weight = (double *)G_malloc(sides * sizeof(double));    for (ct_dir = 0; ct_dir < sides; ct_dir++) {	/* get r, c (r_nbr, c_nbr) for neighbours */	r_nbr = nextdr[ct_dir];	c_nbr = nextdc[ct_dir];	/* account for rare cases when ns_res != ew_res */	dy = ABS(r_nbr) * window.ns_res;	dx = ABS(c_nbr) * window.ew_res;	if (ct_dir < 4)	    dist_to_nbr[ct_dir] = dx + dy;	else	    dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);    }    flag_nbr = (char *)G_malloc(sides * sizeof(char));    wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));    ele_nbr = (CELL *)G_malloc(sides * sizeof(CELL));    workedon = 0;    if (bas_thres <= 0)	threshold = 60;    else	threshold = bas_thres;    for (killer = 0; killer < do_points; killer++) {	G_percent(killer, do_points, 1);	seg_get(&astar_pts, (char *)&point, 0, killer);	r = point.r;	c = point.c;	bseg_get(&asp, &asp_val, r, c);	if (asp_val) {	    dr = r + asp_r[ABS(asp_val)];	    dc = c + asp_c[ABS(asp_val)];	}	/* skip user-defined depressions */	else	    dr = dc = -1;	/* WORKEDFLAG has been set during A* Search	 * reversed meaning here: 0 = done, 1 = not yet done */	bseg_get(&bitflags, &this_flag_value, r, c);	FLAG_UNSET(this_flag_value, WORKEDFLAG);		if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) {	    r_max = dr;	    c_max = dc;	    seg_get(&watalt, (char *)&wa, r, c);	    value = wa.wat;	    /* get weights */	    max_weight = 0;	    sum_weight = 0;	    np_side = -1;	    mfd_cells = 0;	    stream_cells = 0;	    swale_cells = 0;	    astar_not_set = 1;	    ele = wa.ele;	    is_null = 0;	    edge = 0;	    flat = 1;	    /* this loop is needed to get the sum of weights */	    for (ct_dir = 0; ct_dir < sides; ct_dir++) {		/* get r, c (r_nbr, c_nbr) for neighbours */		r_nbr = r + nextdr[ct_dir];		c_nbr = c + nextdc[ct_dir];		weight[ct_dir] = -1;		wat_nbr[ct_dir] = 0;		ele_nbr[ct_dir] = 0;//.........这里部分代码省略.........

int write_grid(struct grid_description *grid_info, struct Map_info *Map, int nbreaks, int out_type){    int i, k, j;    int rows, cols;    int num_v_rows, num_v_cols;    double x, y, x_len, y_len;    double sx, sy;    double width, length;    double next_x, next_y;    double snext_x, snext_y;    double angle, dum;    struct line_pnts *Points;    Points = Vect_new_line_struct();    num_v_rows = grid_info->num_vect_rows;    num_v_cols = grid_info->num_vect_cols;    rows = grid_info->num_rows;    cols = grid_info->num_cols;    width = grid_info->width;    length = grid_info->length;    angle = grid_info->angle;    /*     * For latlon, must draw in shorter sections     * to make sure that each section of the grid     * line is less than half way around the globe     */     x_len = length / (1. * nbreaks + 1);     y_len = width / (1. * nbreaks + 1);    /* write out all the vector lengths (x vectors) of the entire grid  */    G_verbose_message(_("Writing out vector rows..."));    y = grid_info->origin_y;    for (i = 0; i < num_v_rows; ++i) {	double startx;	startx = grid_info->origin_x;	G_percent(i, num_v_rows, 2);	for (k = 0; k < cols; k++) {	    x = startx;            j = 0;	    do {		if (j < nbreaks)		    next_x = x + x_len;		else		    next_x = startx + length;		sx = x;		sy = y;		snext_x = next_x;		dum = y;		rotate(&x, &y, grid_info->origin_x, grid_info->origin_y,		       angle);		rotate(&next_x, &dum, grid_info->origin_x,		       grid_info->origin_y, angle);		write_vect(x, y, next_x, dum, Map, Points, out_type);		y = sy;		x = next_x = snext_x;                j++;	    } while (j <= nbreaks);	    startx += length;	}	y += width;    }    /* write out all the vector widths (y vectors) of the entire grid  */    G_verbose_message(_("Writing out vector columns..."));    x = grid_info->origin_x;    for (i = 0; i < num_v_cols; ++i) {        double starty;	starty = grid_info->origin_y;	G_percent(i, num_v_cols, 2);	for (k = 0; k < rows; k++) {	  y = starty;	  j = 0;	  do {	      if (j < nbreaks)		  next_y = y + y_len;	      else		  next_y = starty + width;	      sx = x;	      sy = y;	      snext_y = next_y;	      dum = x;	      rotate(&x, &y, grid_info->origin_x, grid_info->origin_y, angle);	      rotate(&dum, &next_y, grid_info->origin_x, grid_info->origin_y,		    angle);	      write_vect(x, y, dum, next_y, Map, Points, out_type);	      x = sx;	      y = next_y = snext_y;//.........这里部分代码省略.........

//.........这里部分代码省略.........	    G_fatal_error(_("Too many ETo files. Only %d allowed."),			  MAXFILES);	name2 = *ptr2;	/* Allocate input buffer */	infd2[nfiles2] = Rast_open_old(name2, "");	Rast_get_cellhd(name2, "", &cellhd);	inrast2[nfiles2] = Rast_allocate_d_buf();	nfiles2++;    }    nfiles2--;    if (nfiles2 <= 1)	G_fatal_error(_("The min specified input map is two"));    /* Allocate output buffer, use input map data_type */    nrows = Rast_window_rows();    ncols = Rast_window_cols();    outrast = Rast_allocate_d_buf();    /* Create New raster files */    outfd = Rast_open_new(result, 1);    /*******************/    /* Process pixels */    double doy[MAXFILES];    double sum[MAXFILES];    for (row = 0; row < nrows; row++)     {	DCELL d_out;	DCELL d_ETrF[MAXFILES];	DCELL d[MAXFILES];	DCELL d1[MAXFILES];	DCELL d2[MAXFILES];	G_percent(row, nrows, 2);	/* read input map */	for (i = 1; i <= nfiles; i++) 	    Rast_get_d_row(infd[i], inrast[i], row);		for (i = 1; i <= nfiles1; i++) 	    Rast_get_d_row(infd1[i], inrast1[i], row);	for (i = 1; i <= nfiles2; i++) 	    Rast_get_d_row (infd2[i], inrast2[i], row);	/*process the data */	for (col = 0; col < ncols; col++)         {            int	d1_null=0;            int	d_null=0;	    for (i = 1; i <= nfiles; i++)             {		    if (Rast_is_d_null_value(&((DCELL *) inrast[i])[col]))		    	d_null=1;		    else	                d[i] = ((DCELL *) inrast[i])[col];	    }	    for (i = 1; i <= nfiles1; i++)             {		    if (Rast_is_d_null_value(&((DCELL *) inrast1[i])[col]))			d1_null=1;		    else	                d1[i] = ((DCELL *) inrast1[i])[col];	    }	    for (i = 1; i <= nfiles2; i++) 

/*!   /brief Remove duplicate features from vector map.   Remove duplicate lines of given types from vector map. Duplicate   lines may be optionally written to error map. Input map must be   opened on level 2 for update. Categories are merged.   GV_BUILD_BASE is sufficient.   /param[in,out] Map vector map where duplicate lines will be deleted   /param type type of line to be delete   /param[out] Err vector map where duplicate lines will be written or NULL   /return void */voidVect_remove_duplicates(struct Map_info *Map, int type, struct Map_info *Err){    struct line_pnts *APoints, *BPoints;    struct line_cats *ACats, *BCats, *Cats;    int i, j, c, atype, btype, bline;    int nlines, nbcats_orig;    struct bound_box ABox;    struct boxlist *List;    int ndupl;    APoints = Vect_new_line_struct();    BPoints = Vect_new_line_struct();    ACats = Vect_new_cats_struct();    BCats = Vect_new_cats_struct();    Cats = Vect_new_cats_struct();    List = Vect_new_boxlist(0);    nlines = Vect_get_num_lines(Map);    G_debug(1, "nlines =  %d", nlines);    /* Go through all lines in vector, for each select lines which overlap MBR of     *  this line and check if some of them is identical. If someone is identical     *  remove current line. (In each step just one line is deleted)     */    ndupl = 0;    for (i = 1; i <= nlines; i++) {	G_percent(i, nlines, 1);	if (!Vect_line_alive(Map, i))	    continue;	atype = Vect_read_line(Map, APoints, ACats, i);	if (!(atype & type))	    continue;	Vect_line_box(APoints, &ABox);	Vect_select_lines_by_box(Map, &ABox, type, List);	G_debug(3, "  %d lines selected by box", List->n_values);	for (j = 0; j < List->n_values; j++) {	    bline = List->id[j];	    G_debug(3, "  j = %d bline = %d", j, bline);	    if (i == bline)		continue;	    btype = Vect_read_line(Map, BPoints, BCats, bline);	    /* check for duplicates */	    if (!Vect_line_check_duplicate(APoints, BPoints, Vect_is_3d(Map)))		continue;	    /* Lines area identical -> remove current */	    if (Err) {		Vect_write_line(Err, atype, APoints, ACats);	    }	    Vect_delete_line(Map, i);	    /* Merge categories */	    nbcats_orig = BCats->n_cats;	    for (c = 0; c < ACats->n_cats; c++)		Vect_cat_set(BCats, ACats->field[c], ACats->cat[c]);	    if (BCats->n_cats > nbcats_orig) {		G_debug(4, "cats merged: n_cats %d -> %d", nbcats_orig,			BCats->n_cats);		Vect_rewrite_line(Map, bline, btype, BPoints, BCats);	    }	    ndupl++;	    break;		/* line was deleted -> take the next one */	}	nlines = Vect_get_num_lines(Map);	/* For future when lines with cats will be rewritten */	G_debug(3, "nlines =  %d/n", nlines);    }    G_verbose_message("Removed duplicates: %d", ndupl);}

/* ************************************************************************* */void write_vtk_rgb_data(void *map_r, void *map_g, void *map_b,                        FILE * fp, const char *varname,                        RASTER3D_Region region, int dp){    double value = 0;    int x, y, z, percentage, k;    int rows, cols, depths;    int typeIntern[3];    void *maprgb = NULL;    G_debug(3, "write_vtk_rgb_data: Writing RGB data");    rows = region.rows;    cols = region.cols;    depths = region.depths;    typeIntern[0] = Rast3d_tile_type_map(map_r);    typeIntern[1] = Rast3d_tile_type_map(map_g);    typeIntern[2] = Rast3d_tile_type_map(map_b);    percentage = 0;    /********************** WRITE RGB VOXEL DATA; CELL OR POINT ****************/    fprintf(fp, "COLOR_SCALARS %s 3/n", varname);    for (z = 0; z < depths; z++) {        for (y = 0; y < rows; y++) {            G_percent(percentage, (rows * depths - 1), 10);            percentage++;            for (x = 0; x < cols; x++) {                for (k = 0; k < 3; k++) {                    if (k == 0)                        maprgb = map_r;                    if (k == 1)                        maprgb = map_g;                    if (k == 2)                        maprgb = map_b;                    /* In case of structured grid data, the point/cell coordinates                       are computed based on the default north->south raster3d coordinate system.                       We need to compute south -> north ordering for image data.                     */                    if (!param.structgrid->answer)                        value =                        get_g3d_raster_value_as_double(maprgb, x, rows - y - 1, z,                                                       typeIntern[k],                                                       0.0);                    else                        value =                        get_g3d_raster_value_as_double(maprgb, x, y, z,                                                       typeIntern[k],                                                       0.0);                    /*Test of value range, the data should be 1 byte gray values */                    if (value > 255 || value < 0) {                        G_warning(_("Wrong 3d raster map values! Values should in between 0 and 255!/n"));                        fprintf(fp, "0 ");                    } else {                        fprintf(fp, "%.*f ", dp, (value / 255));                    }                }                fprintf(fp, "/n");            }        }    }    return;}