int report(long rectify, int ok)
{
    int minutes, hours;
    long seconds;
    long ncells;

    G_message("%s", ok ? _("complete") : _("failed"));

    if (!ok)
	return 1;

    seconds = rectify;
    minutes = seconds / 60;
    hours = minutes / 60;
    minutes -= hours * 60;
    ncells = target_window.rows * target_window.cols;
    G_verbose_message(_("%d rows, %d cols (%ld cells) completed in"),
			target_window.rows, target_window.cols, ncells);
    if (hours)
	G_verbose_message(_("%d:%02d:%02ld hours"), hours, minutes, seconds % 60);
    else
	G_verbose_message(_("%d:%02ld minutes"), minutes, seconds % 60);
    if (seconds)
	G_verbose_message(_("%.1f cells per minute"),
			  (60.0 * ncells) / ((double)seconds));
		      
    G_message("-----------------------------------------------");

    return 1;
}

static int manage_memory(int srows, int scols, struct globals *globals)
{
    double reg_size_mb, segs_mb;
    int reg_size_count, nseg, nseg_total;

    /* minimum region size to store in search tree */
    reg_size_mb = 2 * globals->datasize +     /* mean, sum */
                  2 * sizeof(int) +           /* id, count */
		  sizeof(unsigned char) + 
		  2 * sizeof(struct REG_NODE *);
    reg_size_mb /= (1024. * 1024.);

    /* put aside some memory for segment structures */
    segs_mb = globals->mb * 0.1;
    if (segs_mb > 10)
	segs_mb = 10;

    /* calculate number of region stats that can be kept in memory */
    reg_size_count = (globals->mb - segs_mb) / reg_size_mb;
    globals->min_reg_size = 3;
    if (reg_size_count < (double) globals->notnullcells / globals->min_reg_size) {
	globals->min_reg_size = (double) globals->notnullcells / reg_size_count;
    }
    else {
	reg_size_count = (double) globals->notnullcells / globals->min_reg_size;
	/* recalculate segs_mb */
	segs_mb = globals->mb - reg_size_count * reg_size_mb;
    }

    G_verbose_message(_("Regions with at least %d cells are stored in memory"),
                      globals->min_reg_size);

    /* calculate number of segments in memory */
    if (globals->bounds_map != NULL) {
	/* input bands, segment ids, bounds map */
	nseg = (1024. * 1024. * segs_mb) /
	       (sizeof(DCELL) * globals->nbands * srows * scols + 
		sizeof(CELL) * 4 * srows * scols);
    }
    else {
	/* input bands, segment ids */
	nseg = (1024. * 1024. * segs_mb) /
	       (sizeof(DCELL) * globals->nbands * srows * scols + 
		sizeof(CELL) * 2 * srows * scols);
    }
    nseg_total = (globals->nrows / srows + (globals->nrows % srows > 0)) *
                 (globals->ncols / scols + (globals->ncols % scols > 0));

    if (nseg > nseg_total)
	nseg = nseg_total;
    
    G_debug(1, "current region:  %d rows, %d cols", globals->nrows, globals->ncols);
    G_debug(1, "segmented to tiles with size:  %d rows, %d cols", srows,
	    scols);
    G_verbose_message(_("Number of segments in memory: %d of %d total"),
                      nseg, nseg_total);
    
    return nseg;
}

/* check compatibility of map header and region header */
void check_header(char* cellname) {

  const char *mapset;
  mapset = G_find_raster(cellname, "");
  if (mapset == NULL) {
    G_fatal_error(_("Raster map <%s> not found"), cellname);
  }
  /* read cell header */
  struct Cell_head cell_hd;
  Rast_get_cellhd (cellname, mapset, &cell_hd);
  
  /* check compatibility with module region */
  if (!((region->ew_res == cell_hd.ew_res)
		&& (region->ns_res == cell_hd.ns_res))) {
    G_fatal_error(_("cell file %s resolution differs from current region"),
				  cellname);
  } else {
    if (opt->verbose) { 
      G_message(_("cell %s header compatible with region header"),
	      cellname);
      fflush(stderr);
    }
  }


  /* check type of input elevation raster and check if precision is lost */
    RASTER_MAP_TYPE data_type;
	data_type = Rast_map_type(opt->elev_grid, mapset);
#ifdef ELEV_SHORT
	G_verbose_message(_("Elevation stored as SHORT (%dB)"),
		sizeof(elevation_type));
	if (data_type == FCELL_TYPE) {
	  G_warning(_("raster %s is of type FCELL_TYPE "
			"--precision may be lost."), opt->elev_grid); 
	}
	if (data_type == DCELL_TYPE) {
	  G_warning(_("raster %s is of type DCELL_TYPE "
			"--precision may be lost."),  opt->elev_grid);
	}
#endif 
#ifdef ELEV_FLOAT
	G_verbose_message( _("Elevation stored as FLOAT (%dB)"), 
			sizeof(elevation_type));
	if (data_type == CELL_TYPE) {
	  G_warning(_("raster %s is of type CELL_TYPE "
		"--you should use r.terraflow.short"), opt->elev_grid); 
	}
	if (data_type == DCELL_TYPE) {
	  G_warning(_("raster %s is of type DCELL_TYPE "
		"--precision may be lost."),  opt->elev_grid);
	}
#endif
	



}

static void write_bil_hdr(
    const char *outfile, const struct Cell_head *region,
    int bytes, int order, int header, double null_val)
{
    char out_tmp[GPATH_MAX];
    FILE *fp;

    sprintf(out_tmp, "%s.hdr", outfile);
    G_verbose_message(_("Header File = %s"), out_tmp);

    /* Open Header File */
    fp = fopen(out_tmp, "w");
    if (!fp)
	G_fatal_error(_("Unable to create file <%s>"), out_tmp);

    fprintf(fp, "nrows %d\n", region->rows);
    fprintf(fp, "ncols %d\n", region->cols);
    fprintf(fp, "nbands 1\n");
    fprintf(fp, "nbits %d\n", bytes * 8);
    fprintf(fp, "byteorder %s\n", order == 0 ? "M" : "I");
    fprintf(fp, "layout bil\n");
    fprintf(fp, "skipbytes %d\n", header ? 892 : 0);
    fprintf(fp, "nodata %g\n", null_val);

    fclose(fp);
}

int get_conz_points(void)
{
    char msg[200];
    /* struct Ortho_Control_Points cpz; */

    if (!I_get_con_points(group.name, &group.control_points))
	exit(0);

    sprintf(msg, _("Control Z Point file for group [%s] in [%s] \n \n"),
	    group.name, G_mapset());

    G_verbose_message(_("Computing equations..."));

    Compute_ortho_equation();

    switch (group.con_equation_stat) {
    case -1:
	strcat(msg, _("Poorly placed Control Points!\n"));
	strcat(msg, _("Can not generate the transformation equation.\n"));
	strcat(msg, _("Run OPTION 7 of i.ortho.photo again!\n"));
	break;
    case 0:
	strcat(msg, _("No active Control Points!\n"));
	strcat(msg, _("Can not generate the transformation equation.\n"));
	strcat(msg, _("Run OPTION 7 of i.ortho.photo!\n"));
	break;
    default:
	return 1;
    }
    G_fatal_error("%s", msg);
}

/* Init all costs to/from given node */
int init_node_costs(struct Map_info *Map, int from)
{
    int to, ret, row, col;
    double cost;

    G_verbose_message(_("Init costs from node %d"), from);

    for (to = 1; to <= nnodes; to++) {
	if (from == to)
	    continue;
	ret = Vect_net_shortest_path(Map, from, to, NULL, &cost);
	if (ret == -1) {
	    G_debug(1, "Destination node %d is unreachable from node %d\n", to, from);
	    cost = -2;
	}

	if (from < to) {
	    row = from - 1;
	    col = to - from - 1;
	}
	else {
	    row = to - 1;
	    col = from - to - 1;
	}

	G_debug(3, "init costs %d - > %d = %f\n", from, to, cost);
	nodes_costs[row][col] = cost;
    }

    return 1;
}

static int profile(int coords, const char *map, const char *nulls, char **line)
{
    double e1, n1, e2, n2;
    char buf[1024], profile[1024];
    const char *argv[7];
    int argc = 0;
    int n;
    int projection;

    projection = G_projection();

    argv[argc++] = "r.profile";

    if (coords)
	argv[argc++] = "-g";

    sprintf(buf, "input=%s", map);
    argv[argc++] = G_store(buf);

    argv[argc++] = "output=-";

    sprintf(buf, "null_value=%s", nulls);
    argv[argc++] = G_store(buf);

    strcpy(profile, "coordinates=");
    for (n = 0; line[n]; n += 4) {
	int err = parse_line("line", &line[n], &e1, &n1, &e2, &n2, projection);

	if (err) {
	    G_usage();
	    exit(EXIT_FAILURE);
	}

	if (n > 0)
	    strcat(profile, ",");
	G_format_easting(e1, buf, projection);
	strcat(profile, buf);

	G_format_northing(n1, buf, projection);
	strcat(profile, ",");
	strcat(profile, buf);

	G_format_easting(e2, buf, projection);
	strcat(profile, ",");
	strcat(profile, buf);

	G_format_northing(n2, buf, projection);
	strcat(profile, ",");
	strcat(profile, buf);
    }

    argv[argc++] = profile;

    argv[argc++] = NULL;

    G_verbose_message(_("End coordinate: %.15g, %.15g"), e2, n2);

    return G_vspawn_ex(argv[0], argv);
}

int print_time(long *start)
{
    int hours, minutes, seconds;
    long done;

    time(&done);

    seconds = done - *start;
    *start = done;

    hours = seconds / 3600;
    minutes = (seconds - hours * 3600) / 60;
    seconds = seconds % 60;

    if (hours)
	G_verbose_message("%2d:%02d:%02d", hours, minutes, seconds);
    else if (minutes)
	G_verbose_message("%d:%02d", minutes, seconds);
    else
	G_verbose_message("%d seconds", seconds);

    return 0;
}

int rmdac(struct Map_info *Out, struct Map_info *Err)
{
    int i, type, area, ndupl, nlines;

    struct line_pnts *Points;
    struct line_cats *Cats;

    nlines = Vect_get_num_lines(Out);

    Points = Vect_new_line_struct();
    Cats = Vect_new_cats_struct();

    G_debug(1, "nlines =  %d", nlines);

    ndupl = 0;

    for (i = 1; i <= nlines; i++) {
	G_percent(i, nlines, 2);
	if (!Vect_line_alive(Out, i))
	    continue;

	type = Vect_read_line(Out, Points, Cats, i);
	if (!(type & GV_CENTROID))
	    continue;

	area = Vect_get_centroid_area(Out, i);
	G_debug(3, "  area = %d", area);

	if (area < 0) {
	    Vect_delete_line(Out, i);
	    ndupl++;

	    if (Err) {
		Vect_write_line(Err, type, Points, Cats);
	    }
	}
    }

    G_verbose_message(_("Duplicate area centroids: %d"), ndupl);

    Vect_destroy_line_struct(Points);
    Vect_destroy_cats_struct(Cats);

    return ndupl;
}

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

    /* set the window to the cell header */
    if (G_get_cellhd(name, mapset, &cellhd) < 0)
	G_fatal_error(_("Cannot read header of raster map <%s> in <%s>"),
		      name, mapset);

    G_set_window(&cellhd);

    /* open the raster map */
    fd = G_open_cell_old(name, mapset);
    if (fd < 0)
	G_fatal_error(_("Cannot open cell file of raster map <%s> in <%s>"),
		      name, mapset);
    nrows = G_window_rows();
    ncols = G_window_cols();
    cell = G_allocate_cell_buf();
    G_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);
	if (G_get_c_raster_row_nomask(fd, cell, row) < 0)
	    exit(EXIT_SUCCESS);
	G_update_cell_stats(cell, ncols, &statf);
    }
    /* done */
    if (G_verbose() > G_verbose_std())
	G_percent(row, nrows, 2);
    G_close_cell(fd);
    G_free(cell);
    G_rewind_cell_stats(&statf);

    return 0;
}

static void write_bil_wld(const char *outfile, const struct Cell_head *region)
{
    char out_tmp[GPATH_MAX];
    FILE *fp;

    sprintf(out_tmp, "%s.wld", outfile);
    G_verbose_message(_("World File = %s"), out_tmp);

    /* Open World File */
    fp = fopen(out_tmp, "w");
    if (!fp)
	G_fatal_error(_("Unable to create file <%s>"), out_tmp);

    fprintf(fp, "%f\n", region->ew_res);
    fprintf(fp, "0.0\n");
    fprintf(fp, "0.0\n");
    fprintf(fp, "-%f\n", region->ns_res);
    fprintf(fp, "%f\n", region->west + (region->ew_res / 2));
    fprintf(fp, "%f\n", region->north - (region->ns_res / 2));

    fclose(fp);
}

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;
}

int area_area(struct Map_info *In, int *field, struct Map_info *Tmp,
	      struct Map_info *Out, struct field_info *Fi,
	      dbDriver * driver, int operator, int *ofield,
	      ATTRIBUTES * attr, struct ilist *BList, double snap)
{
    int ret, input, line, nlines, area, nareas;
    int in_area, in_centr, out_cat;
    struct line_pnts *Points;
    struct line_cats *Cats;
    CENTR *Centr;
    char buf[1000];
    dbString stmt;
    int nmodif;
    int verbose;

    verbose = G_verbose();

    Points = Vect_new_line_struct();
    Cats = Vect_new_cats_struct();

    /* optional snap */
    if (snap > 0) {
	int i, j, snapped_lines = 0;
	struct bound_box box;
	struct boxlist *boxlist = Vect_new_boxlist(0);
	struct ilist *reflist = Vect_new_list();
	
	G_message(_("Snapping boundaries with %g ..."), snap);

	/* snap boundaries in B to boundaries in A,
	 * not modifying boundaries in A */

	if (BList->n_values > 1)
	    qsort(BList->value, BList->n_values, sizeof(int), cmp_int);

	snapped_lines = 0;
	nlines = BList->n_values;
	for (i = 0; i < nlines; i++) {
	    line = BList->value[i];
	    Vect_read_line(Tmp, Points, Cats, line);
	    /* select lines by box */
	    Vect_get_line_box(Tmp, line, &box);
	    box.E += snap;
	    box.W -= snap;
	    box.N += snap;
	    box.S -= snap;
	    box.T = 0.0;
	    box.B = 0.0;
	    Vect_select_lines_by_box(Tmp, &box, GV_BOUNDARY, boxlist);
	    
	    if (boxlist->n_values > 0) {
		Vect_reset_list(reflist);
		for (j = 0; j < boxlist->n_values; j++) {
		    int aline = boxlist->id[j];

		    if (!bsearch(&aline, BList->value, BList->n_values,
			sizeof(int), cmp_int)) {
			G_ilist_add(reflist, aline);
		    }
		}
		
		/* snap bline to alines */
		if (Vect_snap_line(Tmp, reflist, Points, snap, 0, NULL, NULL)) {
		    /* rewrite bline*/
		    Vect_delete_line(Tmp, line);
		    ret = Vect_write_line(Tmp, GV_BOUNDARY, Points, Cats);
		    G_ilist_add(BList, ret);
		    snapped_lines++;
		    G_debug(3, "line %d snapped", line);
		}
	    }
	}
	Vect_destroy_boxlist(boxlist);
	Vect_destroy_list(reflist);

	G_verbose_message(n_("%d boundary snapped",
                             "%d boundaries snapped",
                             snapped_lines), snapped_lines);
    }

    /* same procedure like for v.in.ogr:
     * Vect_clean_small_angles_at_nodes() can change the geometry so that new intersections
     * are created. We must call Vect_break_lines(), Vect_remove_duplicates()
     * and Vect_clean_small_angles_at_nodes() until no more small dangles are found */
    do {
	G_message(_("Breaking lines..."));
	Vect_break_lines_list(Tmp, NULL, BList, GV_BOUNDARY, NULL);

	/* Probably not necessary for LINE x AREA */
	G_message(_("Removing duplicates..."));
	Vect_remove_duplicates(Tmp, GV_BOUNDARY, NULL);

	G_message(_("Cleaning boundaries at nodes..."));
	nmodif =
	    Vect_clean_small_angles_at_nodes(Tmp, GV_BOUNDARY, NULL);
    } while (nmodif > 0);

    /* ?: May be result of Vect_break_lines() + Vect_remove_duplicates() any dangle or bridge?
     * In that case, calls to Vect_remove_dangles() and Vect_remove_bridges() would be also necessary */

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

//.........这里部分代码省略.........
	    /* check if stats column exists */
	    G_debug(1, "check if stats column exists");
	    db_get_column(Adriver, AFi->table, stats_column_opt->answer,
			  &column);
	    if (column) {
		/* check stats column type */
		if (db_column_Ctype
		    (Adriver, AFi->table,
		     stats_column_opt->answer) != DB_C_TYPE_DOUBLE)
		    G_fatal_error(_("scolumn must be of type double"));

		db_free_column(column);
		column = NULL;
	    }
	    else {
		/* create stats column */
		/* db_add_column() exists but is not implemented,
		 * see lib/db/stubs/add_col.c */
		sprintf(buf, "alter table %s add column %s double",
				AFi->table, stats_column_opt->answer);
		db_set_string(&stmt, buf);
		if (db_execute_immediate(Adriver, &stmt) != DB_OK)
		    G_fatal_error(_("Unable to add column <%s>"),
				  stats_column_opt->answer);
	    }
	}
    }
    else
	AFi = NULL;

    Pdriver = NULL;
    if (method_opt->answer) {

	G_verbose_message(_("collecting attributes from points vector..."));

	PFi = Vect_get_field(&PIn, point_field);
	if (PFi == NULL)
	    G_fatal_error(_("Database connection not defined for layer %d"),
			  point_field);

	Pdriver = db_start_driver_open_database(PFi->driver, PFi->database);
	if (Pdriver == NULL)
	    G_fatal_error(_("Unable to open database <%s> with driver <%s>"),
			  PFi->database, PFi->driver);

	/* check if point column exists */
	db_get_column(Pdriver, PFi->table, point_column_opt->answer, &column);
	if (column) {
	    db_free_column(column);
	    column = NULL;
	}
	else {
	    G_fatal_error(_("Column <%s> not found in table <%s>"),
			  point_column_opt->answer, PFi->table);
	}

	/* Check column type */
	ctype =
	    db_column_Ctype(Pdriver, PFi->table, point_column_opt->answer);

	if (ctype == DB_C_TYPE_INT)
	    half = menu[method].half;
	else if (ctype == DB_C_TYPE_DOUBLE)
	    half = 0;
	else
	    G_fatal_error(_("column for points vector must be numeric"));

//.........这里部分代码省略.........
#endif
    G_get_window(&outcellhd);

    if(gprint_bounds->answer && !print_bounds->answer)
	print_bounds->answer = gprint_bounds->answer;
    curr_proj = G_projection();

    /* Get projection info for output mapset */
    if ((out_proj_info = G_get_projinfo()) == NULL)
	G_fatal_error(_("Unable to get projection info of output raster map"));

    if ((out_unit_info = G_get_projunits()) == NULL)
	G_fatal_error(_("Unable to get projection units of output raster map"));

    if (pj_get_kv(&oproj, out_proj_info, out_unit_info) < 0)
	G_fatal_error(_("Unable to get projection key values of output raster map"));

    /* Change the location           */
    G__create_alt_env();
    G__setenv("GISDBASE", indbase->answer ? indbase->answer : G_gisdbase());
    G__setenv("LOCATION_NAME", inlocation->answer);

    permissions = G__mapset_permissions(setname);
    if (permissions < 0)	/* can't access mapset       */
	G_fatal_error(_("Mapset <%s> in input location <%s> - %s"),
		      setname, inlocation->answer,
		      permissions == 0 ? _("permission denied")
		      : _("not found"));

    /* if requested, list the raster maps in source location - MN 5/2001 */
    if (list->answer) {
	int i;
	char **list;
	G_verbose_message(_("Checking location <%s> mapset <%s>"),
			  inlocation->answer, setname);
	list = G_list(G_ELEMENT_RASTER, G__getenv("GISDBASE"),
		      G__getenv("LOCATION_NAME"), setname);
	for (i = 0; list[i]; i++) {
	    fprintf(stdout, "%s\n", list[i]);
	}
	fflush(stdout);
	exit(EXIT_SUCCESS);	/* leave r.proj after listing */
    }

    if (!inmap->answer)
	G_fatal_error(_("Required parameter <%s> not set"), inmap->key);

    if (!G_find_raster(inmap->answer, setname))
	G_fatal_error(_("Raster map <%s> in location <%s> in mapset <%s> not found"),
		      inmap->answer, inlocation->answer, setname);

    /* Read input map colour table */
    have_colors = Rast_read_colors(inmap->answer, setname, &colr);

    /* Get projection info for input mapset */
    if ((in_proj_info = G_get_projinfo()) == NULL)
	G_fatal_error(_("Unable to get projection info of input map"));

    if ((in_unit_info = G_get_projunits()) == NULL)
	G_fatal_error(_("Unable to get projection units of input map"));

    if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
	G_fatal_error(_("Unable to get projection key values of input map"));

    G_free_key_value(in_proj_info);
    G_free_key_value(in_unit_info);

//.........这里部分代码省略.........
    release_flag->guisection = _("Manage");

    truecolor_flag = G_define_flag();
    truecolor_flag->key = 't';
    truecolor_flag->description = _("Disable true colors");
    truecolor_flag->guisection = _("Settings");

    update_flag = G_define_flag();
    update_flag->key = 'u';
    update_flag->label = _("Open output file in update mode");
    update_flag->description = _("Requires --overwrite flag");
    update_flag->guisection = _("Settings");

    x_flag = G_define_flag();
    x_flag->key = 'x';
    x_flag->label = _("Launch light-weight wx monitor without toolbars and statusbar");
    x_flag->description = _("Requires 'start=wx0-7'");
    x_flag->guisection = _("Settings");

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

    if (x_flag->answer && start_opt->answer && strncmp(start_opt->answer, "wx", 2) != 0)
        G_warning(_("Flag -%c has effect only for wx monitors (%s=wx0-7)"),
                  x_flag->key, start_opt->key);
            
    if (selected_flag->answer || release_flag->answer ||
        cmd_flag->answer || sfile_flag->answer) {
	if (list_flag->answer)
	    G_warning(_("Flag -%c ignored"), list_flag->key);
	mon = G_getenv_nofatal("MONITOR");
	if (mon) {
	    if (selected_flag->answer) {
		G_verbose_message(_("Currently selected monitor:"));
		fprintf(stdout, "%s\n", mon);
	    }
	    else if (cmd_flag->answer) {
		G_message(_("List of commands for monitor <%s>:"), mon);
		list_cmd(mon, stdout);
	    }
            else if (sfile_flag->answer) {
                list_files(mon, stdout);
            }
	    else if (mon) { /* release */
		G_unsetenv("MONITOR");
		G_verbose_message(_("Monitor <%s> released"), mon); 
                ret = stop_mon(mon);
	    }
	}
	else
	    G_important_message(_("No monitor selected"));
	
	exit(EXIT_SUCCESS);
    }

    if (list_flag->answer) {
	print_list(stdout);
	exit(EXIT_SUCCESS);
    }
	
    nopts = 0;
    if (start_opt->answer)
	nopts++;
    if (stop_opt->answer)
	nopts++;
    if (select_opt->answer)

void
subcluster(struct SigSet *S, int Class_Index, int *Max_num, int maxsubclasses)
{
    int nparams_clust;
    int ndata_points;
    int min_i, min_j;
    int nbands;
    double rissanen;
    double min_riss;
    struct ClassSig *Sig;

    static int first = 1;
    static struct SigSet min_S;
    static struct ClassSig *min_Sig;

    /* set class pointer */
    Sig = &(S->ClassSig[Class_Index]);

    /* set number of bands */
    nbands = S->nbands;

    /* allocate scratch class first time subroutine is called */
    if (first) {
	int i;

	I_InitSigSet(&min_S);
	I_SigSetNBands(&min_S, nbands);
	min_Sig = I_NewClassSig(&min_S);

	/* allocate enough subsignatures in scratch space */
	for (i = 0; i < maxsubclasses; i++)
	    I_NewSubSig(&min_S, min_Sig);

	first = 0;
    }

    /* compute number of parameters per cluster */
    nparams_clust = 1 + nbands + 0.5 * (nbands + 1) * nbands;

    /* compute number of data points */
    ndata_points = Sig->ClassData.npixels * nbands - total_nulls;
    if (ndata_points <= 1)
	G_fatal_error("Not enough data points");

    /* Check for too few pixels */
    *Max_num = (ndata_points + 1) / nparams_clust - 1;
    if (maxsubclasses > *Max_num / 2)
	maxsubclasses = *Max_num / 2;
    if (maxsubclasses < 1) {
	G_warning(_("Not enough pixels in class %d"),
		  Class_Index + 1);
	Sig->nsubclasses = 0;
	Sig->used = 0;
	return;
    }

    /* check for too many subclasses */
    if (Sig->nsubclasses > maxsubclasses) {
	Sig->nsubclasses = maxsubclasses;
	G_warning(_("Too many subclasses for class index %d"),
		  Class_Index + 1);
	G_message(_("Number of subclasses set to %d"),
		  Sig->nsubclasses);
    }


    /* initialize clustering */
    seed(Sig, nbands);

    /* EM algorithm */
    min_riss = refine_clusters(Sig, nbands);
    G_debug(1, "Subclasses = %d Rissanen = %f", Sig->nsubclasses,
	      min_riss);
    copy_ClassSig(Sig, min_Sig, nbands);

    G_debug(1, "combine classes");
    while (Sig->nsubclasses > 1) {
	min_i = min_j = 0;
	reduce_order(Sig, nbands, &min_i, &min_j);
	G_verbose_message(_("Combining subclasses (%d,%d)..."), min_i + 1,
			  min_j + 1);
	
	rissanen = refine_clusters(Sig, nbands);
	G_debug(1, "Subclasses = %d; Rissanen = %f", Sig->nsubclasses,
		rissanen);
	if (rissanen < min_riss) {
	    min_riss = rissanen;
	    copy_ClassSig(Sig, min_Sig, nbands);
	}
    }

    copy_ClassSig(min_Sig, Sig, nbands);
}

/*!
   \brief Build partial topology for vector map.

   Should only be used in special cases of vector processing.

   This functions optionally builds only some parts of
   topology. Highest level is specified by build parameter which may
   be:
    - GV_BUILD_NONE - nothing is build
    - GV_BUILD_BASE - basic topology, nodes, lines, spatial index;
    - GV_BUILD_AREAS - build areas and islands, but islands are not attached to areas;
    - GV_BUILD_ATTACH_ISLES - attach islands to areas;
    - GV_BUILD_CENTROIDS - assign centroids to areas, build category index;
    - GV_BUILD_ALL - top level, the same as GV_BUILD_CENTROIDS.
    
   If functions is called with build lower than current value of the
   Map, the level is downgraded to requested value.

   All calls to Vect_write_line(), Vect_rewrite_line(),
   Vect_delete_line() respect the last value of build used in this
   function.

   Note that the functions has effect only if requested level is
   higher than current level, to rebuild part of topology, call first
   downgrade and then upgrade, for example:

   - Vect_build()
   - Vect_build_partial(,GV_BUILD_BASE,)
   - Vect_build_partial(,GV_BUILD_AREAS,) 

   \param Map vector map
   \param build highest level of build

   \return 1 on success
   \return 0 on error
 */
int Vect_build_partial(struct Map_info *Map, int build)
{
    struct Plus_head *plus;
    int ret;

    G_debug(3, "Vect_build(): build = %d", build);

    /* If topology is already build (map on > level 2), set level to 1
     * so that lines will be read by V1_read_ (all lines) */
    Map->level = LEVEL_1; /* may be not needed, because V1_read is used
                             directly by Vect_build_ */

    if (Map->format != GV_FORMAT_OGR_DIRECT &&
        !(Map->format == GV_FORMAT_POSTGIS && Map->fInfo.pg.toposchema_name))
        /* don't write support files for OGR direct and PostGIS Topology */
	Map->support_updated = TRUE;

    if (!Map->plus.Spidx_built) {
	if (Vect_open_sidx(Map, 2) < 0)
	    G_fatal_error(_("Unable to open spatial index file for vector map <%s>"),
			    Vect_get_full_name(Map));
    }

    plus = &(Map->plus);
    if (build > GV_BUILD_NONE && !Map->temporary) {
	G_message(_("Building topology for vector map <%s>..."),
		  Vect_get_full_name(Map));
    }
    plus->with_z = Map->head.with_z;
    plus->spidx_with_z = Map->head.with_z;

    if (build == GV_BUILD_ALL && plus->built < GV_BUILD_ALL) {
	dig_cidx_free(plus);	/* free old (if any) category index */
	dig_cidx_init(plus);
    }

    ret = ((*Build_array[Map->format]) (Map, build));
    if (ret == 0) {
	return 0;
    }

    if (build > GV_BUILD_NONE) {
        Map->level = LEVEL_2;
	G_verbose_message(_("Topology was built"));
    }

    plus->mode = GV_MODE_WRITE;

    if (build == GV_BUILD_ALL) {
	plus->cidx_up_to_date = TRUE;	/* category index was build */
	dig_cidx_sort(plus);
    }

    if (build > GV_BUILD_NONE) {
	G_message(_("Number of nodes: %d"), plus->n_nodes);
	G_message(_("Number of primitives: %d"), plus->n_lines);
	G_message(_("Number of points: %d"), plus->n_plines);
	G_message(_("Number of lines: %d"), plus->n_llines);
	G_message(_("Number of boundaries: %d"), plus->n_blines);
	G_message(_("Number of centroids: %d"), plus->n_clines);

	if (plus->n_flines > 0)
	    G_message(_("Number of faces: %d"), plus->n_flines);

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

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

    /*  
     *  fprintf(stdout, maxdistance->answer);
     *  fprintf(stdout, "Maxd is %f", maxd);
     *  fprintf(stdout, xcoord->answer);
     *  fprintf(stdout, "xval is %f", xval);
     *  fprintf(stdout, ycoord->answer);
     *  fprintf(stdout, "yval is %f", yval);
     */

    if (maxd == 0.0) {
	G_get_window(&window);
	x = window.proj;
	G_format_resolution(window.ew_res, ewres, x);
	G_format_resolution(window.ns_res, nsres, x);
	EW_DIST1 =
	    G_distance(window.east, window.north, window.west, window.north);
	/* EW Dist at South Edge */
	EW_DIST2 =
	    G_distance(window.east, window.south, window.west, window.south);
	/* NS Dist at East edge */
	NS_DIST1 =
	    G_distance(window.east, window.north, window.east, window.south);
	/* NS Dist at West edge */
	NS_DIST2 =
	    G_distance(window.west, window.north, window.west, window.south);
	xres = ((EW_DIST1 + EW_DIST2) / 2) / window.cols;
	yres = ((NS_DIST1 + NS_DIST2) / 2) / window.rows;
	if (xres > yres)
	    maxd = xres;
	else
	    maxd = yres;
    }

    /* Look at maps given on command line */
    if (vect) {
	for (i = 0; vect[i]; i++) ;
	nvects = i;

	Map = (struct Map_info *)G_malloc(nvects * sizeof(struct Map_info));

	width = mwidth = 0;
	for (i = 0; i < nvects; i++) {
	    str = strchr(vect[i], '@');
	    if (str)
		j = str - vect[i];
	    else
		j = strlen(vect[i]);
	    if (j > width)
		width = j;

	    mapset = G_find_vector2(vect[i], "");
	    if (!mapset)
		G_fatal_error(_("Vector map <%s> not found"), vect[i]);

	    j = strlen(mapset);
	    if (j > mwidth)
		mwidth = j;

	    level = Vect_open_old(&Map[i], vect[i], mapset);
	    if (level < 2)
		G_fatal_error(_("You must build topology on vector map <%s>"),
			      vect[i]);

	    G_verbose_message(_("Building spatial index..."));
	    Vect_build_spatial_index(&Map[i]);
	}
    }

    if (!coords_opt->answer) {
	/* if coords are not given on command line, read them from stdin */
	setvbuf(stdin, NULL, _IOLBF, 0);
	setvbuf(stdout, NULL, _IOLBF, 0);
	while (fgets(buf, sizeof(buf), stdin) != NULL) {
	    ret = sscanf(buf, "%lf%c%lf", &xval, &ch, &yval);
	    if (ret == 3 && (ch == ',' || ch == ' ' || ch == '\t')) {
		what(xval, yval, maxd, width, mwidth, topo_flag->answer,
		     printattributes->answer, shell_flag->answer);
	    }
	    else {
		G_warning(_("Unknown input format, skipping: '%s'"), buf);
		continue;
	    }
	}
    }
    else {
	/* use coords given on command line */
	for (i = 0; coords_opt->answers[i] != NULL; i += 2) {
	    xval = atof(coords_opt->answers[i]);
	    yval = atof(coords_opt->answers[i + 1]);
	    what(xval, yval, maxd, width, mwidth, topo_flag->answer,
		 printattributes->answer, shell_flag->answer);
	}
    }

    for (i = 0; i < nvects; i++)
	Vect_close(&Map[i]);

    exit(EXIT_SUCCESS);
}

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

    iloc_name = ilocopt->answer;

    if (ibaseopt->answer)
	gbase = ibaseopt->answer;
    else
	gbase = G_store(G_gisdbase());

    if (!ibaseopt->answer && strcmp(iloc_name, G_location()) == 0)
	G_fatal_error(_("Input and output locations can not be the same"));

    lmax = atof(smax->answer);
    if (lmax < 0)
	lmax = 0;

    Out_proj = G_projection();
    if (Out_proj == PROJECTION_LL && flag.wrap->answer)
	nowrap = 1;
    
    G_begin_distance_calculations();

    /* Change the location here and then come back */

    select_target_env();
    G_setenv_nogisrc("GISDBASE", gbase);
    G_setenv_nogisrc("LOCATION_NAME", iloc_name);
    stat = G_mapset_permissions(iset_name);
    
    if (stat >= 0) {		/* yes, we can access the mapset */
	/* if requested, list the vector maps in source location - MN 5/2001 */
	if (flag.list->answer) {
	    int i;
	    char **list;
	    G_verbose_message(_("Checking location <%s> mapset <%s>"),
			      iloc_name, iset_name);
	    list = G_list(G_ELEMENT_VECTOR, G_getenv_nofatal("GISDBASE"),
			  G_getenv_nofatal("LOCATION_NAME"), iset_name);
	    if (list[0]) {
		for (i = 0; list[i]; i++) {
		    fprintf(stdout, "%s\n", list[i]);
		}
		fflush(stdout);
	    }
	    else {
		G_important_message(_("No vector maps found"));
	    }
	    exit(EXIT_SUCCESS);	/* leave v.proj after listing */
	}

	if (mapopt->answer == NULL) {
	    G_fatal_error(_("Required parameter <%s> not set"), mapopt->key);
	}

	G_setenv_nogisrc("MAPSET", iset_name);
	/* Make sure map is available */
	mapset = G_find_vector2(map_name, iset_name);
	if (mapset == NULL)
	    G_fatal_error(_("Vector map <%s> in location <%s> mapset <%s> not found"),
			  map_name, iloc_name, iset_name);

	 /*** Get projection info for input mapset ***/
	in_proj_keys = G_get_projinfo();
	if (in_proj_keys == NULL)
	    exit(EXIT_FAILURE);

	/* apparently the +over switch must be set in the input projection,