/* ************************************************************************* */
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.
 */
void
remove_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");