int get_range(const char *name, long *min, long *max)
{
    struct Range range;
    int nrows, ncols, row, col;
    CELL *cell;
    int fd;
    CELL cmin, cmax;
    struct Cell_head cellhd;

    if (Rast_read_range(name, "", &range) < 0) {
	Rast_init_range(&range);	/* read the file to get the range */
	Rast_get_cellhd(name, "", &cellhd);
	Rast_set_window(&cellhd);
	cell = Rast_allocate_c_buf();
	fd = Rast_open_old(name, "");
	nrows = Rast_window_rows();
	ncols = Rast_window_cols();
	G_message(_("Reading %s ..."), name);
	for (row = 0; row < nrows; row++) {
	    G_percent(row, nrows, 2);
	    Rast_get_c_row_nomask(fd, cell, row);
	    for (col = 0; col < ncols; col++)
		Rast_update_range(cell[col], &range);
	}
	G_percent(row, nrows, 2);
	Rast_close(fd);
	G_free(cell);
    }

    Rast_get_range_min_max(&range, &cmin, &cmax);
    *min = cmin;
    *max = cmax;

    return 0;
}

int open_map(MAPS* rast) {
	
	int row, col;
	int fd;
	char* mapset;
	struct Cell_head cellhd;
	int bufsize;
	void* tmp_buf;
	
    mapset = (char*)G_find_raster2(rast->elevname, "");	
	
	    if (mapset == NULL)
		G_fatal_error(_("Raster map <%s> not found"), rast->elevname);
	
		rast->fd = Rast_open_old(rast->elevname, mapset);
		Rast_get_cellhd(rast->elevname, mapset, &cellhd);
		rast->raster_type = Rast_map_type(rast->elevname, mapset);

    if (window.ew_res < cellhd.ew_res || window.ns_res < cellhd.ns_res)
	G_warning(_("Region resolution shoudn't be lesser than map %s resolution. Run g.region rast=%s to set proper resolution"),
		      rast->elevname, rast->elevname);

		tmp_buf=Rast_allocate_buf(rast->raster_type);
		rast->elev = (FCELL**) G_malloc((row_buffer_size+1) * sizeof(FCELL*));
	
	for (row = 0; row < row_buffer_size+1; ++row) {
		rast->elev[row] = Rast_allocate_buf(FCELL_TYPE);
		Rast_get_row(rast->fd, tmp_buf,row, rast->raster_type);
				for (col=0;col<ncols;++col)
			get_cell(col, rast->elev[row], tmp_buf, rast->raster_type);
  } /* end elev */

G_free(tmp_buf);
return 0;
}

int get_stats(const char *name, struct Cell_stats *statf)
{
    int fd;
    CELL *cell;
    int row, nrows, ncols;

    fd = Rast_open_old(name, "");
    nrows = Rast_window_rows();
    ncols = Rast_window_cols();
    cell = Rast_allocate_c_buf();

    Rast_init_cell_stats(statf);
    G_message(_("Reading %s ..."), name);
    for (row = 0; row < nrows; row++) {
	G_percent(row, nrows, 2);
	Rast_get_c_row(fd, cell, row);
	Rast_update_cell_stats(cell, ncols, statf);
    }
    if (row < nrows)
	exit(1);
    Rast_close(fd);
    G_free(cell);
    G_percent(row, nrows, 2);

    return 0;
}

int cseg_read_cell(CSEG * cseg, char *map_name, char *mapset)
{
    GW_LARGE_INT row, nrows;
    int map_fd;
    CELL *buffer;

    cseg->name = NULL;
    cseg->mapset = NULL;

    map_fd = Rast_open_old(map_name, mapset);
    nrows = Rast_window_rows();
    buffer = Rast_allocate_c_buf();
    for (row = 0; row < nrows; row++) {
	Rast_get_c_row(map_fd, buffer, row);
	if (Segment_put_row(&(cseg->seg), buffer, row) < 0) {
	    G_free(buffer);
	    Rast_close(map_fd);
	    G_warning("%s(): unable to segment put row for [%s] in [%s]",
		    me, map_name, mapset);
	    return (-1);
	}
    }

    Rast_close(map_fd);
    G_free(buffer);

    cseg->name = G_store(map_name);
    cseg->mapset = G_store(mapset);

    return 0;
}

/* ************************************************************************* */
int open_input_map(const char *name, const char *mapset)
{
    G_debug(3, "Open raster file %s in mapset %s", name, mapset);

    /* open raster map */
    return Rast_open_old(name, mapset);
}

int distdrop ( cell_map *elev,
               cell_map *dist, cell_map *dir,
               cell_map *up, cell_map *dw,
               seg_map *segment_info, move *movements,
               queue **redo_segments)
{
    int all_done = 1;

    int **neighbours = new_int_map ( ( int ) sizeof ( movements ),
                                     ( int ) sizeof ( movements[0] ), NULL );

    elev->fd = Rast_open_old ( elev->name, "" );
    elev->type = Rast_get_map_type( elev->fd );

    init_seg_map(elev, segment_info);
    copy_segment(elev, 0);

    //print_map_with_seg ( elev, segment_info );

    while ( all_done ){ // if all_done != 0? continue: break
        all_done = queue_pixel ( redo_segments,  elev, dist,
                                 dir, up, dw, segment_info,
                                 movements, neighbours);
    }
    return 0;
}

void openfiles(char *r_name, char *g_name, char *b_name,
	       char *h_name, char *i_name, char *s_name,
	       int fd_input[3], int fd_output[3], CELL * rowbuf[3])
{
    fd_input[0] = Rast_open_old(r_name, "");
    fd_input[1] = Rast_open_old(g_name, "");
    fd_input[2] = Rast_open_old(b_name, "");

    /* open output files */
    fd_output[0] = Rast_open_c_new(h_name);
    fd_output[1] = Rast_open_c_new(i_name);
    fd_output[2] = Rast_open_c_new(s_name);

    /* allocate the cell row buffer */
    rowbuf[0] = Rast_allocate_c_buf();
    rowbuf[1] = Rast_allocate_c_buf();
    rowbuf[2] = Rast_allocate_c_buf();
}

/*!
   \brief Load raster map as integer map

   Calling function must have already allocated space in buff for
   struct BM of wind->rows & wind->cols.

   This routine simply loads the map into the bitmap by repetitve calls
   to get_map_row.  Any value other than 0 in the map will set the bitmap.
   (may want to change later to allow specific value to set)

   Changed to use null.

   \param wind current window
   \param map_name raster map name
   \param[out] buff data buffer

   \returns 1 on success
   \return -1 on failure
 */
int Gs_loadmap_as_bitmap(struct Cell_head *wind, const char *map_name,
			 struct BM *buff)
{
    FILEDESC cellfile;
    const char *map_set;
    int *tmp_buf;
    int row, col;

    G_debug(3, "Gs_loadmap_as_bitmap");

    map_set = G_find_raster2(map_name, "");
    if (!map_set) {
	G_warning(_("Raster map <%s> not found"), map_name);
	return -1;
    }

    cellfile = Rast_open_old(map_name, map_set);

    tmp_buf = (int *)G_malloc(wind->cols * sizeof(int));	/* G_fatal_error */
    if (!tmp_buf) {
	return -1;
    }

    G_message(_("Loading raster map <%s>..."),
	      G_fully_qualified_name(map_name, map_set));

    for (row = 0; row < wind->rows; row++) {
	Rast_get_c_row(cellfile, tmp_buf, row);

	for (col = 0; col < wind->cols; col++) {
	    if (Rast_is_c_null_value(&tmp_buf[col])) {
		/* no data */
		BM_set(buff, col, row, 1);
	    }
	    else {
		BM_set(buff, col, row, 0);
	    }
	}
    }

    Rast_close(cellfile);

    G_free(tmp_buf);

    return (1);
}

/*!
   \brief Load raster map as floating point map

   Calling function must have already allocated space in buff for
   wind->rows * wind->cols floats.

   This routine simply loads the map into a 2d array by repetitve calls
   to get_f_raster_row.

   \param wind current window
   \param map_name raster map name
   \param[out] buff data buffer
   \param[out] nullmap null map buffer
   \param[out] has_null indicates if raster map contains null-data

   \return 1 on success
   \return 0 on failure
 */
int Gs_loadmap_as_float(struct Cell_head *wind, const char *map_name,
			float *buff, struct BM *nullmap, int *has_null)
{
    FILEDESC cellfile;
    const char *map_set;
    int offset, row, col;

    G_debug(3, "Gs_loadmap_as_float(): name=%s", map_name);

    map_set = G_find_raster2(map_name, "");
    if (!map_set) {
	G_warning(_("Raster map <%s> not found"), map_name);
	return 0;
    }
    *has_null = 0;

    cellfile = Rast_open_old(map_name, map_set);

    G_message(_("Loading raster map <%s>..."),
	      G_fully_qualified_name(map_name, map_set));

    for (row = 0; row < wind->rows; row++) {
	offset = row * wind->cols;
	Rast_get_f_row(cellfile, &(buff[offset]), row);

	G_percent(row, wind->rows, 2);

	for (col = 0; col < wind->cols; col++) {
	    if (Rast_is_f_null_value(buff + offset + col)) {
		*has_null = 1;
		BM_set(nullmap, col, row, 1);
	    }
	    /* set nm */
	}
    }
    G_percent(1, 1, 1);

    G_debug(4, "  has_null=%d", *has_null);

    Rast_close(cellfile);

    return (1);
}

/*!
   \brief Open and allocate space for the group band files.

   \param refer pointer to band files structure
   \param[out] band_buffer buffer to read one row of each band
   \param[out] band_fd band files descriptors
 */
void open_band_files(struct Ref *refer, CELL *** band_buffer, int **band_fd)
{
    int n, nbands;

    char *name, *mapset;

    G_debug(3, "open_band_files()");

    /* allocate row buffers and open raster maps */
    nbands = refer->nfiles;
    *band_buffer = (CELL **) G_malloc(nbands * sizeof(CELL *));
    *band_fd = (int *)G_malloc(nbands * sizeof(int));

    for (n = 0; n < nbands; n++) {
	(*band_buffer)[n] = Rast_allocate_c_buf();
	name = refer->file[n].name;
	mapset = refer->file[n].mapset;
	(*band_fd)[n] = Rast_open_old(name, mapset);
    }
}

void rast_segment_open(SEGMENT * segment, const char *name,
                       RASTER_MAP_TYPE * map_type)
{
    /* TODO: check if not passing the mapset is OK */
    int rowio = Rast_open_old(name, "");

    *map_type = Rast_get_map_type(rowio);
    int segment_rows = 64;

    /* we use long segments because this is how the values a binned */
    int segment_cols = Rast_input_window_cols();
    int segments_in_memory = 4;

    if (Segment_open(segment, G_tempfile(), Rast_input_window_rows(),
                     Rast_input_window_cols(), segment_rows, segment_cols,
                     Rast_cell_size(*map_type), segments_in_memory) != 1)
        G_fatal_error(_("Cannot create temporary file with segments of a raster map"));
    rast_segment_load(segment, rowio, *map_type);
    Rast_close(rowio);          /* we won't need the raster again */
}

int get_cats(const char *name, const char *mapset)
{
    int fd;
    int row, nrows, ncols;
    CELL *cell;
    struct Cell_head cellhd;

    /* set the window to the cell header */
    Rast_get_cellhd(name, mapset, &cellhd);

    Rast_set_window(&cellhd);

    /* open the raster map */
    fd = Rast_open_old(name, mapset);
    nrows = Rast_window_rows();
    ncols = Rast_window_cols();
    cell = Rast_allocate_c_buf();
    Rast_init_cell_stats(&statf);

    /* read the raster map */
    G_verbose_message(_("Reading <%s> in <%s>"), name, mapset);
    for (row = 0; row < nrows; row++) {
	if (G_verbose() > G_verbose_std())
	    G_percent(row, nrows, 2);
	Rast_get_c_row_nomask(fd, cell, row);
	Rast_update_cell_stats(cell, ncols, &statf);
    }
    /* done */
    if (G_verbose() > G_verbose_std())
	G_percent(row, nrows, 2);
    Rast_close(fd);
    G_free(cell);
    Rast_rewind_cell_stats(&statf);

    return 0;
}

/* 
 * do_histogram() - Creates histogram for CELL
 *
 * RETURN: EXIT_SUCCESS / EXIT_FAILURE
 */
int do_histogram(const char *name)
{
    CELL *cell;
    struct Cell_head cellhd;
    struct Cell_stats statf;
    int nrows, ncols;
    int row;
    int fd;

    Rast_get_cellhd(name, "", &cellhd);

    Rast_set_window(&cellhd);
    fd = Rast_open_old(name, "");

    nrows = Rast_window_rows();
    ncols = Rast_window_cols();
    cell = Rast_allocate_c_buf();

    Rast_init_cell_stats(&statf);
    for (row = 0; row < nrows; row++) {
	Rast_get_c_row_nomask(fd, cell, row);
	Rast_update_cell_stats(cell, ncols, &statf);
    }

    if (row == nrows)
	Rast_write_histogram_cs(name, &statf);

    Rast_free_cell_stats(&statf);
    Rast_close(fd);
    G_free(cell);

    if (row < nrows)
	return -1;

    return 0;
}

int camera_angle(char *name)
{
    int row, col, nrows, ncols;
    double XC = group.XC;
    double YC = group.YC;
    double ZC = group.ZC;
    double c_angle, c_angle_min, c_alt, c_az, slope, aspect;
    double radians_to_degrees = 180.0 / M_PI;
    /* double degrees_to_radians = M_PI / 180.0; */
    DCELL e1, e2, e3, e4, e5, e6, e7, e8, e9;
    double factor, V, H, dx, dy, dz, key;
    double north, south, east, west, ns_med;
    FCELL *fbuf0, *fbuf1, *fbuf2, *tmpbuf, *outbuf;
    int elevfd, outfd;
    struct Cell_head cellhd;
    struct Colors colr;
    FCELL clr_min, clr_max;
    struct History hist;
    char *type;

    G_message(_("Calculating camera angle to local surface..."));
    
    select_target_env();
    
    /* align target window to elevation map, otherwise we get artefacts
     * like in r.slope.aspect -a */
     
    Rast_get_cellhd(elev_name, elev_mapset, &cellhd);

    Rast_align_window(&target_window, &cellhd);
    Rast_set_window(&target_window);
    
    elevfd = Rast_open_old(elev_name, elev_mapset);
    if (elevfd < 0) {
	G_fatal_error(_("Could not open elevation raster"));
	return 1;
    }

    nrows = target_window.rows;
    ncols = target_window.cols;
    
    outfd = Rast_open_new(name, FCELL_TYPE);
    fbuf0 = Rast_allocate_buf(FCELL_TYPE);
    fbuf1 = Rast_allocate_buf(FCELL_TYPE);
    fbuf2 = Rast_allocate_buf(FCELL_TYPE);
    outbuf = Rast_allocate_buf(FCELL_TYPE);
    
    /* give warning if location units are different from meters and zfactor=1 */
    factor = G_database_units_to_meters_factor();
    if (factor != 1.0)
	G_warning(_("Converting units to meters, factor=%.6f"), factor);

    G_begin_distance_calculations();
    north = Rast_row_to_northing(0.5, &target_window);
    ns_med = Rast_row_to_northing(1.5, &target_window);
    south = Rast_row_to_northing(2.5, &target_window);
    east = Rast_col_to_easting(2.5, &target_window);
    west = Rast_col_to_easting(0.5, &target_window);
    V = G_distance(east, north, east, south) * 4;
    H = G_distance(east, ns_med, west, ns_med) * 4;
    
    c_angle_min = 90;
    Rast_get_row(elevfd, fbuf1, 0, FCELL_TYPE);
    Rast_get_row(elevfd, fbuf2, 1, FCELL_TYPE);

    for (row = 0; row < nrows; row++) {
	G_percent(row, nrows, 2);
	
	Rast_set_null_value(outbuf, ncols, FCELL_TYPE);

	/* first and last row */
	if (row == 0 || row == nrows - 1) {
	    Rast_put_row(outfd, outbuf, FCELL_TYPE);
	    continue;
	}
	
	tmpbuf = fbuf0;
	fbuf0 = fbuf1;
	fbuf1 = fbuf2;
	fbuf2 = tmpbuf;
	
	Rast_get_row(elevfd, fbuf2, row + 1, FCELL_TYPE);

	north = Rast_row_to_northing(row + 0.5, &target_window);

	for (col = 1; col < ncols - 1; col++) {
	    
	    e1 = fbuf0[col - 1];
	    if (Rast_is_d_null_value(&e1))
		continue;
	    e2 = fbuf0[col];
	    if (Rast_is_d_null_value(&e2))
		continue;
	    e3 = fbuf0[col + 1];
	    if (Rast_is_d_null_value(&e3))
		continue;
	    e4 = fbuf1[col - 1];
	    if (Rast_is_d_null_value(&e4))
		continue;
	    e5 = fbuf1[col];
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    int i, j, nlines, type, field, cat;
    int fd;

    /* struct Categories RCats; *//* TODO */
    struct Cell_head window;
    RASTER_MAP_TYPE out_type;
    CELL *cell;
    DCELL *dcell;
    double drow, dcol;
    char buf[2000];
    struct Option *vect_opt, *rast_opt, *field_opt, *col_opt, *where_opt;
    int Cache_size;
    struct order *cache;
    int cur_row;
    struct GModule *module;

    struct Map_info Map;
    struct line_pnts *Points;
    struct line_cats *Cats;
    int point;
    int point_cnt;		/* number of points in cache */
    int outside_cnt;		/* points outside region */
    int nocat_cnt;		/* points inside region but without category */
    int dupl_cnt;		/* duplicate categories */
    struct bound_box box;

    int *catexst, *cex;
    struct field_info *Fi;
    dbString stmt;
    dbDriver *driver;
    int select, norec_cnt, update_cnt, upderr_cnt, col_type;

    G_gisinit(argv[0]);

    module = G_define_module();
    G_add_keyword(_("vector"));
    G_add_keyword(_("raster"));
    G_add_keyword(_("position"));
    G_add_keyword(_("querying"));
    G_add_keyword(_("attribute table"));
    module->description =
	_("Uploads raster values at positions of vector points to the table.");

    vect_opt = G_define_standard_option(G_OPT_V_INPUT);
    vect_opt->key = "vector";
    vect_opt->description =
	_("Name of input vector points map for which to edit attribute table");

    rast_opt = G_define_standard_option(G_OPT_R_INPUT);
    rast_opt->key = "raster";
    rast_opt->description = _("Name of existing raster map to be queried");

    field_opt = G_define_standard_option(G_OPT_V_FIELD);

    col_opt = G_define_option();
    col_opt->key = "column";
    col_opt->type = TYPE_STRING;
    col_opt->required = YES;
    col_opt->description =
	_("Column name (will be updated by raster values)");

    where_opt = G_define_standard_option(G_OPT_DB_WHERE);

    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);


    field = atoi(field_opt->answer);

    db_init_string(&stmt);
    Points = Vect_new_line_struct();
    Cats = Vect_new_cats_struct();

    G_get_window(&window);
    Vect_region_box(&window, &box);	/* T and B set to +/- PORT_DOUBLE_MAX */

    /* Open vector */
    Vect_set_open_level(2);
    Vect_open_old(&Map, vect_opt->answer, "");

    Fi = Vect_get_field(&Map, field);
    if (Fi == NULL)
	G_fatal_error(_("Database connection not defined for layer %d"),
		      field);

    /* Open driver */
    driver = db_start_driver_open_database(Fi->driver, Fi->database);
    if (driver == NULL) {
	G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
		      Fi->database, Fi->driver);
    }

    /* Open raster */
    fd = Rast_open_old(rast_opt->answer, "");

    out_type = Rast_get_map_type(fd);

    /* TODO: Later possibly category labels */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	b = nrows;
	l = 0;
	r = ncols;
    }

    D_set_src(t, b, l, r);
    D_update_conversions();

    /* figure out arrow scaling if using a magnitude map */
    if (opt7->answer) {
	Rast_init_fp_range(&range);	/* really needed? */
	if (Rast_read_fp_range(mag_map, "", &range) != 1)
	    G_fatal_error(_("Problem reading range file"));
	Rast_get_fp_range_min_max(&range, &mag_min, &mag_max);

	scale *= 1.5 / fabs(mag_max);
	G_debug(3, "scaling=%.2f  rast_max=%.2f", scale, mag_max);
    }

    if (grid_color > 0) {	/* ie not "none" */
	/* Set color */
	D_use_color(grid_color);

	/* Draw vertical grids */
	for (col = 0; col < ncols; col++)
	    D_line_abs(col, 0, col, nrows);

	/* Draw horizontal grids */
	for (row = 0; row < nrows; row++)
	    D_line_abs(0, row, ncols, row);
    }

    /* open the raster map */
    layer_fd = Rast_open_old(layer_name, "");

    raster_type = Rast_get_map_type(layer_fd);

    /* allocate the cell array */
    raster_row = Rast_allocate_buf(raster_type);


    if (opt7->answer) {
	/* open the magnitude raster map */
	mag_fd = Rast_open_old(mag_map, "");

	mag_raster_type = Rast_get_map_type(mag_fd);

	/* allocate the cell array */
	mag_raster_row = Rast_allocate_buf(mag_raster_type);
    }


    /* loop through cells, find value, determine direction (n,s,e,w,ne,se,sw,nw),
       and call appropriate function to draw an arrow on the cell */

    for (row = 0; row < nrows; row++) {
	Rast_get_row(layer_fd, raster_row, row, raster_type);
	ptr = raster_row;

	if (opt7->answer) {
	    Rast_get_row(mag_fd, mag_raster_row, row, mag_raster_type);
	    mag_ptr = mag_raster_row;
	}

	for (col = 0; col < ncols; col++) {

int main(int argc, char *argv[])
{
    struct GModule *module;
    int infile;
    const char *mapset;
    size_t cell_size;
    int ytile, xtile, y, overlap;
    int *outfiles;
    void *inbuf;

    G_gisinit(argv[0]);

    module = G_define_module();
    G_add_keyword(_("raster"));
    G_add_keyword(_("tiling"));
    module->description =
	_("Splits a raster map into tiles.");

    parm.rastin = G_define_standard_option(G_OPT_R_INPUT);

    parm.rastout = G_define_option();
    parm.rastout->key = "output";
    parm.rastout->type = TYPE_STRING;
    parm.rastout->required = YES;
    parm.rastout->multiple = NO;
    parm.rastout->description = _("Output base name");

    parm.width = G_define_option();
    parm.width->key = "width";
    parm.width->type = TYPE_INTEGER;
    parm.width->required = YES;
    parm.width->multiple = NO;
    parm.width->description = _("Width of tiles (columns)");

    parm.height = G_define_option();
    parm.height->key = "height";
    parm.height->type = TYPE_INTEGER;
    parm.height->required = YES;
    parm.height->multiple = NO;
    parm.height->description = _("Height of tiles (rows)");

    parm.overlap = G_define_option();
    parm.overlap->key = "overlap";
    parm.overlap->type = TYPE_INTEGER;
    parm.overlap->required = NO;
    parm.overlap->multiple = NO;
    parm.overlap->description = _("Overlap of tiles");

    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);

    G_get_set_window(&src_w);
    overlap = parm.overlap->answer ? atoi(parm.overlap->answer) : 0;

    mapset = G_find_raster2(parm.rastin->answer, "");
    if (mapset == NULL)
        G_fatal_error(_("Raster map <%s> not found"), parm.rastin->answer);

    /* set window to old map */
    Rast_get_cellhd(parm.rastin->answer, "", &src_w);
    dst_w = src_w;
    dst_w.cols = atoi(parm.width->answer);
    dst_w.rows = atoi(parm.height->answer);
    G_adjust_Cell_head(&dst_w, 1, 1);

    xtiles = (src_w.cols + dst_w.cols - 1) / dst_w.cols;
    ytiles = (src_w.rows + dst_w.rows - 1) / dst_w.rows;

    G_debug(1, "X: %d * %d, Y: %d * %d",
	    xtiles, dst_w.cols, ytiles, dst_w.rows);

    src_w.cols = xtiles * dst_w.cols + 2 * overlap;
    src_w.rows = ytiles * dst_w.rows + 2 * overlap;
    src_w.west = src_w.west - overlap * src_w.ew_res;
    src_w.east = src_w.west + (src_w.cols + 2 * overlap) * src_w.ew_res;
    src_w.north = src_w.north + overlap * src_w.ns_res;
    src_w.south = src_w.north - (src_w.rows + 2 * overlap) * src_w.ns_res;

    Rast_set_input_window(&src_w);

    /* set the output region */
    ovl_w = dst_w;
    ovl_w.cols = ovl_w.cols + 2 * overlap;
    ovl_w.rows = ovl_w.rows + 2 * overlap;

    G_adjust_Cell_head(&ovl_w, 1, 1);
    Rast_set_output_window(&ovl_w);

    infile = Rast_open_old(parm.rastin->answer, "");
    map_type = Rast_get_map_type(infile);
    cell_size = Rast_cell_size(map_type);

    inbuf = Rast_allocate_input_buf(map_type);

    outfiles = G_malloc(xtiles * sizeof(int));

    G_debug(1, "X: %d * %d, Y: %d * %d",
	    xtiles, dst_w.cols, ytiles, dst_w.rows);

    G_message(_("Generating %d x %d = %d tiles..."), xtiles, ytiles, xtiles * ytiles);
//.........这里部分代码省略.........

/* SAMPLING UNIT SETUP DRIVER */
void sample(int t0, int b0, int l0, int r0, char *name, char *name1,
	    char *name2, double *msc)
{
    int btn, d, fmask;
    double tmp;

    /*
       IN
       name = name of raster map to be set up
       name1 = name of overlay vector map
       name2 = name of overlay site map
       msc[0]= cols of region/width of screen
       msc[1]= rows of region/height of screen
       t0 = top row of sampling frame
       b0 = bottom row of sampling frame
       l0 = left col of sampling frame
       r0 = right col of sampling frame
     */

    /* determine whether the user will use
       the keyboard or mouse to setup the 
       sampling units */

  keyboard:
    fprintf(stderr, "\n\n    HOW WILL YOU SPECIFY SAMPLING UNITS?\n");
    fprintf(stderr,
	    "\n       Use keyboard to enter sampling unit dimensions   1");
    fprintf(stderr,
	    "\n       Use mouse to draw sampling units                 2\n");
    fprintf(stderr,
	    "\n                                            Which Number?  ");
    numtrap(1, &tmp);
    d = (int)(tmp);
    if (d < 1 || d > 2) {
	fprintf(stderr, "     You did not enter a 1 or 2, try again\n");
	goto keyboard;
    }

    if (d == 1 || d == 2) {
	/* return a value > 0 to fmask if there is
	   a MASK present */

	fprintf(stderr,
		"\n    If a MASK is not present (see r.mask) a beep may sound\n");
	fprintf(stderr,
		"    and a WARNING may be printed that can be ignored.\n");
	fprintf(stderr,
		"    If a MASK is present there will be no warning.\n");
	fmask = Rast_open_old("MASK", G_mapset());
	fprintf(stderr, "\n");


	/* call the routine to setup sampling
	   units manually */

	if (d == 1)
	    man_unit(t0, b0, l0, r0, name, name1, name2, msc, fmask);

	/* call the routine to setup sampling
	   units graphically */

	else if (d == 2)
	    graph_unit(t0, b0, l0, r0, name, name1, name2, msc, fmask);

	Rast_close(fmask);
    }
    /* if neither, then exit */
    else
	exit(0);

    return;
}

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

	    fprintf(stderr, "\n    Do you want to sample using rectangles");

	    if (!G_yes
		("\n       (including squares) (y) or circles (n)?   ", 1)) {
	      getradius:
		fprintf(stderr,
			"\n    What radius do you want for the circles?  Radius");
		fprintf(stderr,
			"\n       is in pixels; add 0.5 pixels, for the center");
		fprintf(stderr,
			"\n       pixel, to the number of pixels outside the");
		fprintf(stderr,
			"\n       center pixel.  Type a real number with one");
		fprintf(stderr,
			"\n       decimal place ending in .5 (e.g., 4.5):        ");

		numtrap(1, &radius);
		if (radius > 100.0) {
		    fprintf(stderr,
			    "\n    Are you sure that you want such a large");

		    if (!G_yes("\n       radius (> 100 pixels)?   ", 1))
			goto getradius;
		}

		ratio = 1.0;
		u_w = (int)(2 * radius);
		u_l = (int)(2 * radius);

		if (fmask > 0) {
		    count = 0;
		    row_buf = Rast_allocate_buf(CELL_TYPE);
		    fr = Rast_open_old(n1, G_mapset());
		    for (j = t; j < b; j++) {
			Rast_zero_buf(row_buf, CELL_TYPE);
			Rast_get_row(fr, row_buf, j, CELL_TYPE);
			for (k = l; k < r; k++) {
			    if (*(row_buf + k))
				count++;
			}
		    }
		    G_free(row_buf);
		    Rast_close(fr);
		    cnt = (double)(count);
		    if (cnt)
			cnt = sqrt(cnt);
		    else
			cnt = 0;
		}
		else {
		    count = (w_l - (int)(start[0])) * (w_w - (int)(start[1]));
		}
	    }

	    /* if sampling using rectangles/squares */

	    else {

		/* get the width/length ratio */

	      getratio:
		fprintf(stderr,
			"\n    Sampling unit SHAPE (aspect ratio, #cols/#rows) "
			"expressed as real number"
			"\n    (e.g., 10 cols/5 rows = 2.0) for sampling units "

void filter_holes(Gfile * out)
{
    int row, col, nrows, ncols;

    void *arast, *brast, *crast;
    int i, pixel[9], cold, warm, shadow, nulo, lim;

    Gfile tmp;

    nrows = Rast_window_rows();
    ncols = Rast_window_cols();

    if (nrows < 3 || ncols < 3)
        return;

    /* Open to read */
    if ((out->fd = Rast_open_old(out->name, "")) < 0)
        G_fatal_error(_("Unable to open raster map <%s>"), out->name);

    arast = Rast_allocate_buf(CELL_TYPE);
    brast = Rast_allocate_buf(CELL_TYPE);
    crast = Rast_allocate_buf(CELL_TYPE);

    /* Open to write */
    sprintf(tmp.name, "_%d.BBB", getpid());
    tmp.rast = Rast_allocate_buf(CELL_TYPE);
    if ((tmp.fd = Rast_open_new(tmp.name, CELL_TYPE)) < 0)
        G_fatal_error(_("Unable to create raster map <%s>"), tmp.name);

    G_important_message(_("Filling small holes in clouds..."));

    /* Se puede acelerar creandolos nulos y luego arast = brast
       brast = crast y cargando crast solamente
       G_set_f_null_value(cell[2], ncols);
     */

    for (row = 0; row < nrows; row++) {
        G_percent(row, nrows, 2);
        /* Read row values */
        if (row != 0) {
            Rast_get_c_row(out->fd, arast, row - 1);
        }
        Rast_get_c_row(out->fd, brast, row);
        if (row != (nrows - 1)) {
            Rast_get_c_row(out->fd, crast, row + 1);
        }
        /* Analysis of all pixels */
        for (col = 0; col < ncols; col++) {
            pixel[0] = pval(brast, col);
            if (pixel[0] == 0) {
                if (row == 0) {
                    pixel[1] = -1;
                    pixel[2] = -1;
                    pixel[3] = -1;
                    if (col == 0) {
                        pixel[4] = -1;
                        pixel[5] = pval(brast, col + 1);
                        pixel[6] = -1;
                        pixel[7] = pval(crast, col);
                        pixel[8] = pval(crast, col + 1);
                    }
                    else if (col != (ncols - 1)) {
                        pixel[4] = pval(brast, col - 1);
                        pixel[5] = pval(brast, col + 1);
                        pixel[6] = pval(crast, col - 1);
                        pixel[7] = pval(crast, col);
                        pixel[8] = pval(crast, col + 1);
                    }
                    else {
                        pixel[4] = pval(brast, col - 1);
                        pixel[5] = -1;
                        pixel[6] = pval(crast, col - 1);
                        pixel[7] = pval(crast, col);
                        pixel[8] = -1;
                    }
                }
                else if (row != (nrows - 1)) {
                    if (col == 0) {
                        pixel[1] = -1;
                        pixel[2] = pval(arast, col);
                        pixel[3] = pval(arast, col + 1);
                        pixel[4] = -1;
                        pixel[5] = pval(brast, col + 1);
                        pixel[6] = -1;
                        pixel[7] = pval(crast, col);
                        pixel[8] = pval(crast, col + 1);
                    }
                    else if (col != (ncols - 1)) {
                        pixel[1] = pval(arast, col - 1);
                        pixel[2] = pval(arast, col);
                        pixel[3] = pval(arast, col + 1);
                        pixel[4] = pval(brast, col - 1);
                        pixel[5] = pval(brast, col + 1);
                        pixel[6] = pval(crast, col - 1);
                        pixel[7] = pval(crast, col);
                        pixel[8] = pval(crast, col + 1);
                    }
                    else {
                        pixel[1] = pval(arast, col - 1);
                        pixel[2] = pval(arast, col);
//.........这里部分代码省略.........