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;
//.........这里部分代码省略.........