//.........这里部分代码省略.........
	    if (left == 0 || rite == 0) {
		/* invalid boundary */
		ret = 1;
		break;
	    }
	}
    }
    if (ret) {
	Vect_remove_dangles(Tmp, GV_BOUNDARY, -1, NULL);
	Vect_remove_bridges(Tmp, NULL, NULL, NULL);
    }

    G_set_verbose(0);
    Vect_build_partial(Tmp, GV_BUILD_NONE);
    Vect_build_partial(Tmp, GV_BUILD_BASE);
    G_set_verbose(verbose);
    G_message(_("Merging lines..."));
    Vect_merge_lines(Tmp, GV_BOUNDARY, NULL, NULL);

    /* Attach islands */
    G_message(_("Attaching islands..."));
    /* can take some time, show messages */
    Vect_build_partial(Tmp, GV_BUILD_ATTACH_ISLES);

    /* Calculate new centroids for all areas */
    nareas = Vect_get_num_areas(Tmp);

    Centr = (CENTR *) G_malloc((nareas + 1) * sizeof(CENTR));	/* index from 1 ! */
    for (area = 1; area <= nareas; area++) {
	ret =
	    Vect_get_point_in_area(Tmp, area, &(Centr[area].x),
				   &(Centr[area].y));
	if (ret < 0) {
	    G_warning(_("Cannot calculate area centroid"));
	    Centr[area].valid = 0;
	}
	else {
	    Centr[area].valid = 1;
	}
    }

    /* Query input maps */
    for (input = 0; input < 2; input++) {
	G_message(_("Querying vector map <%s>..."),
		  Vect_get_full_name(&(In[input])));

	for (area = 1; area <= nareas; area++) {
	    Centr[area].cat[input] = Vect_new_cats_struct();

	    G_percent(area, nareas, 1);

	    in_area =
		Vect_find_area(&(In[input]), Centr[area].x, Centr[area].y);
	    if (in_area > 0) {
		in_centr = Vect_get_area_centroid(&(In[input]), in_area);
		if (in_centr > 0) {
		    int i;

		    Vect_read_line(&(In[input]), NULL, Cats, in_centr);
		    /* Add all cats with original field number */
		    for (i = 0; i < Cats->n_cats; i++) {
			if (Cats->field[i] == field[input]) {
			    ATTR *at;

			    Vect_cat_set(Centr[area].cat[input], ofield[input + 1],
					 Cats->cat[i]);

//.........这里部分代码省略.........
    source = parm.source->answer;
    output = parm.output->answer;

    flip = 0;
    if (flag.h->answer)
	flip |= FLIP_H;
    if (flag.v->answer)
	flip |= FLIP_V;

    if (parm.title->answer) {
	title = G_store(parm.title->answer);
	G_strip(title);
    }
    else
	title = NULL;

    if (!input && !source)
	G_fatal_error(_("One of options <%s> or <%s> must be given"),
		      parm.input->key, parm.source->key);

    if (input && source)
	G_fatal_error(_("Option <%s> and <%s> are mutually exclusive"),
		      parm.input->key, parm.source->key);
    
    if (input && !G_is_absolute_path(input)) {
	char path[GPATH_MAX];
	getcwd(path, sizeof(path));
	strcat(path, "/");
	strcat(path, input);
	input = G_store(path);
    }

    if (!input)
	input = source;

    hDS = GDALOpen(input, GA_ReadOnly);
    if (hDS == NULL)
	return 1;

    setup_window(&cellhd, hDS, &flip);

    check_projection(&cellhd, hDS, flag.o->answer);

    Rast_set_window(&cellhd);

    if (parm.band->answer)
	min_band = max_band = atoi(parm.band->answer);
    else
	min_band = 1, max_band = GDALGetRasterCount(hDS);

    G_verbose_message(_("Proceeding with import..."));

    if (max_band > min_band) {
	if (I_find_group(output) == 1)
	    G_warning(_("Imagery group <%s> already exists and will be overwritten."), output);
	I_init_group_ref(&reference);
    }

    for (band = min_band; band <= max_band; band++) {
	char *output2, *title2 = NULL;

	G_message(_("Reading band %d of %d..."),
		  band, GDALGetRasterCount( hDS ));

	hBand = GDALGetRasterBand(hDS, band);
	if (!hBand)
	    G_fatal_error(_("Selected band (%d) does not exist"), band);

	if (max_band > min_band) {
	    G_asprintf(&output2, "%s.%d", output, band);
	    if (title)
		G_asprintf(&title2, "%s (band %d)", title, band);
	    G_debug(1, "Adding raster map <%s> to group <%s>", output2, output);
	    I_add_file_to_group_ref(output2, G_mapset(), &reference);
	}
	else {
	    output2 = G_store(output);
	    if (title)
		title2 = G_store(title);
	}

	query_band(hBand, output2, flag.r->answer, &cellhd, &info);
	create_map(input, band, output2, &cellhd, &info, title, flip);

	G_free(output2);
	G_free(title2);
    }

    if (flag.e->answer)
	update_default_window(&cellhd);

    /* Create the imagery group if multiple bands are imported */
    if (max_band > min_band) {
    	I_put_group_ref(output, &reference);
	I_put_group(output);
	G_message(_("Imagery group <%s> created"), output);
    }

    exit(EXIT_SUCCESS);
}

//.........这里部分代码省略.........
	_("Print path to support files of currently selected monitor and exit");

    select_flag = G_define_flag();
    select_flag->key = 's';
    select_flag->description = _("Do not automatically select when starting");
    select_flag->guisection = _("Manage");

    release_flag = G_define_flag();
    release_flag->key = 'r';
    release_flag->description = _("Release and stop currently selected monitor and exit");
    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) {

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

static int compute_constants(
				/* invert matrix and compute Sig->SubSig[i].cnst          */
				/* Returns singular=1 if a singluar subcluster was found. */
				/* Returns singular=2 if all subclusters were singular.   */
				/* When singular=2 then nsubclasses=0.                    */
				struct ClassSig *Sig, int nbands)
{
    int i, j;
    int b1, b2;
    int singular;
    double det;
    double pi_sum;

    static int first = 1;
    static int *indx;
    static double **y;
    static double *col;


    /* allocate memory first time subroutine is called */
    if (first) {
	indx = G_alloc_ivector(nbands);
	y = G_alloc_matrix(nbands, nbands);
	col = G_alloc_vector(nbands);
	first = 0;
    }

    G_debug(2, "compute_constants()");
    /* invert matrix and compute constant for each subclass */
    i = 0;
    singular = 0;
    do {
	for (b1 = 0; b1 < nbands; b1++)
	    for (b2 = 0; b2 < nbands; b2++)
		Sig->SubSig[i].Rinv[b1][b2] = Sig->SubSig[i].R[b1][b2];

	invert(Sig->SubSig[i].Rinv, nbands, &det, indx, y, col);
	if (det <= ZERO) {
	    if (Sig->nsubclasses == 1) {
		Sig->nsubclasses--;
		singular = 2;
		G_warning(_("Unreliable clustering. "
			    "Try a smaller initial number of clusters"));
	    }
	    else {
		for (j = i; j < Sig->nsubclasses - 1; j++)
		    copy_SubSig(&(Sig->SubSig[j + 1]), &(Sig->SubSig[j]),
				nbands);
		Sig->nsubclasses--;
		singular = 1;
		G_warning(_("Removed a singular subsignature number %d (%d remain)"),
			  i + 1, Sig->nsubclasses);
		if (Sig->nsubclasses < 0)	/* MN added 12/2001: to avoid endless loop */
		    Sig->nsubclasses = 1;
	    }
	}
	else {
	    Sig->SubSig[i].cnst =
		(-nbands / 2.0) * log(2 * M_PI) - 0.5 * log(det);
	    i++;
	}
    } while (i < Sig->nsubclasses);

    /* renormalize pi */
    pi_sum = 0;
    for (i = 0; i < Sig->nsubclasses; i++)
	pi_sum += Sig->SubSig[i].pi;
    for (i = 0; i < Sig->nsubclasses; i++)
	Sig->SubSig[i].pi /= pi_sum;

    return (singular);
}

//.........这里部分代码省略.........
		    c_nbr = c_nbr + asp_c[(int)af.asp];
		    
		    cseg_get(&stream, &stream_nbr, r_nbr, c_nbr);
		    if (stream_nbr <= 0)
			G_fatal_error(_("Stream id not set while tracing"));

		    Vect_append_point(Points, west_offset + c_nbr * ew_res,
				      north_offset - r_nbr * ns_res, 0);
		    if (stream_nbr != stream_id) {
			/* first point of parent stream */
			break;
		    }
		    seg_get(&aspflag, (char *)&af, r_nbr, c_nbr);
		}

		Vect_write_line(&Out, GV_LINE, Points, Cats);

		top--;
	    }
	}
    }
    G_percent(n_outlets, n_outlets, 1);	/* finish it */

    G_message(_("Writing attribute data..."));

    /* Prepeare strings for use in db_* calls */
    db_init_string(&dbsql);
    db_init_string(&valstr);
    db_init_string(&table_name);
    db_init_handle(&handle);

    /* Preparing database for use */
    /* Create database for new vector map */
    Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
    driver = db_start_driver_open_database(Fi->driver,
					   Vect_subst_var(Fi->database,
							          &Out));
    if (driver == NULL) {
	G_fatal_error(_("Unable to start driver <%s>"), Fi->driver);
    }
    db_set_error_handler_driver(driver);

    G_debug(1, "table: %s", Fi->table);
    G_debug(1, "driver: %s", Fi->driver);
    G_debug(1, "database: %s", Fi->database);

    sprintf(buf,
	    "create table %s (%s integer, stream_type varchar(20), type_code integer)",
	    Fi->table, cat_col_name);
    db_set_string(&dbsql, buf);

    if (db_execute_immediate(driver, &dbsql) != DB_OK) {
	db_close_database(driver);
	db_shutdown_driver(driver);
	G_fatal_error(_("Unable to create table: '%s'"), db_get_string(&dbsql));
    }

    if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
	G_warning(_("Unable to create index on table <%s>"), Fi->table);

    if (db_grant_on_table(driver, Fi->table, DB_PRIV_SELECT,
			  DB_GROUP | DB_PUBLIC) != DB_OK)
	G_fatal_error(_("Unable to grant privileges on table <%s>"), Fi->table);

    db_begin_transaction(driver);

    /* stream nodes */
    for (i = 1; i <= n_stream_nodes; i++) {

	sprintf(buf, "insert into %s values ( %lld, \'%s\', %d )",
		Fi->table, i,
		(stream_node[i].n_trib > 0 ? "intermediate" : "start"),
		(stream_node[i].n_trib > 0));

	db_set_string(&dbsql, buf);

	if (db_execute_immediate(driver, &dbsql) != DB_OK) {
	    db_close_database(driver);
	    db_shutdown_driver(driver);
	    G_fatal_error(_("Unable to insert new row: '%s'"),
			  db_get_string(&dbsql));
	}
    }

    db_commit_transaction(driver);
    db_close_database_shutdown_driver(driver);

    Vect_map_add_dblink(&Out, 1, NULL, Fi->table,
			cat_col_name, Fi->database, Fi->driver);

    G_debug(1, "close vector");

    Vect_hist_command(&Out);
    Vect_build(&Out);
    Vect_close(&Out);

    G_free(nodestack);

    return 1;
}

int main(int argc, char *argv[])
{
    struct GModule *module;
    struct Option *map_opt, *type_opt, *field_opt, *col_opt, *where_opt,
	*percentile;
    struct Flag *shell_flag, *extended;
    struct Map_info Map;
    struct field_info *Fi;
    dbDriver *Driver;
    dbCatValArray Cvarr;
    struct line_pnts *Points;
    struct line_cats *Cats;
    int otype, ofield;
    int compatible = 1;		/* types are compatible: point+centroid or line+boundary or area */
    int nrec, ctype, nlines, line, nareas, area;
    int nmissing = 0;		/* number of missing atttributes */
    int nnull = 0;		/* number of null values */
    int first = 1;

    /* Statistics */
    int count = 0;		/* number of features with non-null attribute */
    double sum = 0.0;
    double sumsq = 0.0;
    double sumcb = 0.0;
    double sumqt = 0.0;
    double sum_abs = 0.0;
    double min = 0.0 / 0.0;	/* init as nan */
    double max = 0.0 / 0.0;
    double mean, mean_abs, pop_variance, sample_variance, pop_stdev,
	sample_stdev, pop_coeff_variation, kurtosis, skewness;
    double total_size = 0.0;	/* total size: length/area */

    /* Extended statistics */
    int perc;

    module = G_define_module();
    G_add_keyword(_("vector"));
    G_add_keyword(_("statistics"));
    module->label =
	_("Calculates univariate statistics for attribute.");
    module->description = _("Variance and standard "
			    "deviation is calculated only for points if specified.");

    map_opt = G_define_standard_option(G_OPT_V_MAP);

    field_opt = G_define_standard_option(G_OPT_V_FIELD);

    type_opt = G_define_standard_option(G_OPT_V_TYPE);
    type_opt->options = "point,line,boundary,centroid,area";

    col_opt = G_define_standard_option(G_OPT_DB_COLUMN);
    col_opt->required = YES;

    where_opt = G_define_standard_option(G_OPT_DB_WHERE);

    percentile = G_define_option();
    percentile->key = "percentile";
    percentile->type = TYPE_INTEGER;
    percentile->required = NO;
    percentile->options = "0-100";
    percentile->answer = "90";
    percentile->description =
	_("Percentile to calculate (requires extended statistics flag)");

    shell_flag = G_define_flag();
    shell_flag->key = 'g';
    shell_flag->description = _("Print the stats in shell script style");

    extended = G_define_flag();
    extended->key = 'e';
    extended->description = _("Calculate extended statistics");

    G_gisinit(argv[0]);

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

    otype = Vect_option_to_types(type_opt);
    perc = atoi(percentile->answer);

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

    /* open input vector */
    Vect_set_open_level(2);
    Vect_open_old2(&Map, map_opt->answer, "", field_opt->answer);
    ofield = Vect_get_field_number(&Map, field_opt->answer);

    /* Check if types are compatible */
    if ((otype & GV_POINTS) && ((otype & GV_LINES) || (otype & GV_AREA)))
	compatible = 0;
    if ((otype & GV_LINES) && (otype & GV_AREA))
	compatible = 0;

    if (!compatible) {
	G_warning(_("Incompatible vector type(s) specified, only number of features, minimum, maximum and range "
		   "can be calculated"));
    }

    if (extended->answer && !(otype & GV_POINTS)) {
//.........这里部分代码省略.........

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

    add_f = G_define_flag();
    add_f->key = 'a';
    add_f->description = _("Add points on nodes");

    /* options and flags parser */
    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);
    /* TODO: make an option for this */
    mask_type = GV_LINE | GV_BOUNDARY;

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

    Vect_check_input_output_name(map_in->answer, map_out->answer,
				 GV_FATAL_EXIT);

    Vect_set_open_level(2);

    if (1 > Vect_open_old(&In, map_in->answer, ""))
	G_fatal_error(_("Unable to open vector map <%s>"), map_in->answer);

    with_z = Vect_is_3d(&In);

    if (0 > Vect_open_new(&Out, map_out->answer, with_z)) {
	Vect_close(&In);
	G_fatal_error(_("Unable to create vector map <%s>"), map_out->answer);
    }


    if (geo_f->answer) {
	geo = 1;
	if (G_projection() != PROJECTION_LL)
	    G_warning(_("The current projection is not longitude-latitude"));
    }
    else
	geo = 0;

    /* parse filter option and select appropriate lines */
    layer = atoi(field_opt->answer);
    chcat =
	(NetA_initialise_varray
	 (&In, layer, mask_type, where_opt->answer, cat_opt->answer,
	  &varray) == 1);

    /* Create table */
    Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
    Vect_map_add_dblink(&Out, 1, NULL, Fi->table, "cat", Fi->database,
			Fi->driver);
    db_init_string(&sql);
    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);

    db_init_string(&tmp);
    if (deg_opt->answer)
	append_string(&tmp, deg_opt->answer);
    if (close_opt->answer)
	append_string(&tmp, close_opt->answer);
    if (betw_opt->answer)
	append_string(&tmp, betw_opt->answer);
    if (eigen_opt->answer)
	append_string(&tmp, eigen_opt->answer);
    sprintf(buf,
	    "create table %s(cat integer%s)", Fi->table, db_get_string(&tmp));

/* Returns 0 - ok , 1 - error */
int
plot(int ctype, struct Map_info *Map, int type, int field,
     char *columns, int ncols, char *sizecol, int size, double scale,
     COLOR * ocolor, COLOR * colors, int y_center, double *max_reference,
     int do3d)
{
    int ltype, nlines, line, col, more, coltype, nselcols;
    double x, y, csize, len;
    struct line_pnts *Points;
    struct line_cats *Cats;
    int cat;
    double *val;
    char buf[2000];
    struct field_info *Fi;
    dbDriver *driver;
    dbValue *value;
    dbString sql;
    dbCursor cursor;
    dbTable *table;
    dbColumn *column;

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

    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_warning(_("Unable to open database <%s> by driver <%s>"),
		  Fi->database,
		  Fi->driver);
	return 1;
    }
    db_set_error_handler_driver(driver);

    val = (double *)G_malloc((ncols + 1) * sizeof(double));	/* + 1 for sizecol */

    Vect_rewind(Map);

    nlines = Vect_get_num_lines(Map);

    /* loop through each vector feature */
    for (line = 1; line <= nlines; line++) {
	G_debug(3, "line = %d", line);
	ltype = Vect_read_line(Map, Points, Cats, line);

	if (!(ltype & type))
	    continue;

	Vect_cat_get(Cats, field, &cat);
	if (cat < 0)
	    continue;

	/* Select values from DB */
	if (ctype == CTYPE_PIE && sizecol != NULL) {
	    sprintf(buf, "select %s, %s from %s where %s = %d", columns,
		    sizecol, Fi->table, Fi->key, cat);
	    nselcols = ncols + 1;
	}
	else {
	    sprintf(buf, "select %s from %s where %s = %d", columns,
		    Fi->table, Fi->key, cat);
	    nselcols = ncols;
	}

	db_set_string(&sql, buf);
	G_debug(3, "SQL: %s", buf);

	if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) !=
	    DB_OK) {
	    G_warning(_("Unable to open select cursor: '%s'"),
		      buf);
	    return 1;
	}

	table = db_get_cursor_table(&cursor);
	if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK || !more)
	    continue;

	for (col = 0; col < nselcols; col++) {
	    column = db_get_table_column(table, col);
	    value = db_get_column_value(column);
	    coltype = db_sqltype_to_Ctype(db_get_column_sqltype(column));
	    switch (coltype) {
	    case DB_C_TYPE_INT:
		val[col] = (double)db_get_value_int(value);
		break;
	    case DB_C_TYPE_DOUBLE:
		val[col] = db_get_value_double(value);
		break;
	    default:
		G_warning("Column type not supported (must be INT or FLOAT)");
		return 1;
	    }
//.........这里部分代码省略.........

int display_label(struct Map_info *Map, int type,
		  struct cat_list *Clist, LATTR *lattr, int chcat)
{
    int ltype;
    struct line_pnts *Points;
    struct line_cats *Cats;
    int ogr_centroids;

    const struct Format_info *finfo;
    
    Points = Vect_new_line_struct();
    Cats = Vect_new_cats_struct();

    Vect_rewind(Map);

    ogr_centroids = FALSE;
    finfo = Vect_get_finfo(Map);
    if (Vect_maptype(Map) == GV_FORMAT_OGR ||
	(Vect_maptype(Map) == GV_FORMAT_POSTGIS &&
         finfo->pg.toposchema_name == NULL)) {
	if (Vect_level(Map) < 2)
	    G_warning(_("Topology level required for drawing centroids "
			"for OGR layers"));
	else if (Vect_get_num_primitives(Map, GV_CENTROID) > 0 &&
		 type & GV_CENTROID)
	    /* label centroids from topo, don't label boundaries */
	    ogr_centroids = TRUE;
    }
    
    while (TRUE) {
	ltype = Vect_read_next_line(Map, Points, Cats);
	if (ltype == -1)
	    G_fatal_error(_("Unable to read vector map"));
	else if (ltype == -2) /* EOF */
	    break;
	
        if (!(type & ltype) && !((type & GV_AREA) && (ltype & GV_CENTROID)))
	    continue;		/* used for both lines and labels */
	
	if (ogr_centroids && ltype == GV_BOUNDARY)
	    /* do not label boundaries */
	    continue;

	process_line(ltype, Points, Cats, lattr, chcat, Clist);
    }

    if (ogr_centroids) {
	/* show label for centroids stored in topo (for OGR layers
	   only) */
	int line, nlines;
	struct bound_box box;
	struct boxlist *list;
	
	list = Vect_new_boxlist(FALSE); /* bboxes not needed */
	Vect_get_constraint_box(Map, &box);
	nlines = Vect_select_lines_by_box(Map, &box, GV_CENTROID, list);
	G_debug(3, "ncentroids (ogr) = %d", nlines);
	
	for (line = 0; line < nlines; line++) {
	    ltype = Vect_read_line(Map, Points, Cats, list->id[line]);
	    process_line(ltype, Points, Cats, lattr, chcat, Clist);
	}
	Vect_destroy_boxlist(list);
    }

    Vect_destroy_line_struct(Points);
    Vect_destroy_cats_struct(Cats);

    return 0;
}

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

    if (strcmp("none", opt4->answer) == 0)
	grid_color = -1;
    else
	grid_color = D_translate_color(opt4->answer);


    if (strcmp("grass", opt2->answer) == 0)
	map_type = 1;
    else if (strcmp("agnps", opt2->answer) == 0)
	map_type = 2;
    else if (strcmp("answers", opt2->answer) == 0)
	map_type = 3;
    else if (strcmp("compass", opt2->answer) == 0)
	map_type = 4;


    scale = atof(opt8->answer);
    if (scale <= 0.0)
	G_fatal_error(_("Illegal value for scale factor"));

    skip = atoi(opt9->answer);
    if (skip <= 0)
	G_fatal_error(_("Illegal value for skip factor"));


    if (opt7->answer) {
	if (map_type != 1 && map_type != 4)
	    G_fatal_error(_("Magnitude is only supported for GRASS and compass aspect maps."));

	mag_map = opt7->answer;
    }
    else if (scale != 1.0)
	G_warning(_("Scale option requires magnitude_map"));


    /* Setup driver and check important information */
    if (D_open_driver() != 0)
      	G_fatal_error(_("No graphics device selected. "
			"Use d.mon to select graphics device."));
    
    D_setup(0);

    /* Read in the map window associated with window */
    G_get_window(&window);

    if (align->answer) {
	struct Cell_head wind;

	Rast_get_cellhd(layer_name, "", &wind);

	/* expand window extent by one wind resolution */
	wind.west += wind.ew_res * ((int)((window.west - wind.west) / wind.ew_res) - (window.west < wind.west));
	wind.east += wind.ew_res * ((int)((window.east - wind.east) / wind.ew_res) + (window.east > wind.east));
	wind.south += wind.ns_res * ((int)((window.south - wind.south) / wind.ns_res) - (window.south < wind.south));
	wind.north += wind.ns_res * ((int)((window.north - wind.north) / wind.ns_res) + (window.north > wind.north));

	wind.rows = (wind.north - wind.south) / wind.ns_res;
	wind.cols = (wind.east - wind.west) / wind.ew_res;

	Rast_set_window(&wind);

	nrows = wind.rows;
	ncols = wind.cols;

	t = (wind.north - window.north) * nrows / (wind.north - wind.south);

//.........这里部分代码省略.........
    Vect_rewind(Map);
    while (1) {
	ltype = Vect_read_next_line(Map, Points, Cats);
	if (ltype == -1)
	    G_fatal_error(_("Unable to read vector map"));
	else if (ltype == -2)		/* EOF */
	    break;

	if (!(type & ltype) && !((type & GV_AREA) && (ltype & GV_CENTROID)))
	    continue;		/* used for both lines and labels */

	D_RGB_color(lattr->color.R, lattr->color.G, lattr->color.B);
	D_text_size(lattr->size, lattr->size);
	if (lattr->font)
	    D_font(lattr->font);
	if (lattr->enc)
	    D_encoding(lattr->enc);

	if (chcat) {
	    int found = 0;

	    for (i = 0; i < Cats->n_cats; i++) {
		if (Cats->field[i] == Clist->field &&
		    Vect_cat_in_cat_list(Cats->cat[i], Clist)) {
		    found = 1;
		    break;
		}
	    }
	    if (!found)
		continue;
	}
	else if (Clist->field > 0) {
	    int found = 0;

	    for (i = 0; i < Cats->n_cats; i++) {
		if (Cats->field[i] == Clist->field) {
		    found = 1;
		    break;
		}
	    }
	    /* lines with no category will be displayed */
	    if (Cats->n_cats > 0 && !found)
		continue;
	}

	if (Vect_cat_get(Cats, lattr->field, &cat)) {
	    int ncats = 0;

	    /* Read attribute from db */
	    db_free_string(&text);
	    for (i = 0; i < Cats->n_cats; i++) {
		int nrows;

		if (Cats->field[i] != lattr->field)
		    continue;
		db_init_string(&stmt);
		sprintf(buf, "select %s from %s where %s = %d", attrcol,
			fi->table, fi->key, Cats->cat[i]);
		G_debug(2, "SQL: %s", buf);
		db_append_string(&stmt, buf);

		if (db_open_select_cursor
		    (driver, &stmt, &cursor, DB_SEQUENTIAL) != DB_OK)
		    G_fatal_error(_("Unable to open select cursor: '%s'"),
				  db_get_string(&stmt));

		nrows = db_get_num_rows(&cursor);

		if (ncats > 0)
		    db_append_string(&text, "/");

		if (nrows > 0) {
		    table = db_get_cursor_table(&cursor);
		    column = db_get_table_column(table, 0);	/* first column */

		    if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK)
			continue;
		    db_convert_column_value_to_string(column, &valstr);

		    db_append_string(&text, db_get_string(&valstr));
		}
		else {
		    G_warning(_("No attribute found for cat %d: %s"), cat,
			      db_get_string(&stmt));
		}

		db_close_cursor(&cursor);
		ncats++;
	    }

	    show_label_line(Points, ltype, lattr, db_get_string(&text));
	}
    }

    db_close_database_shutdown_driver(driver);
    Vect_destroy_line_struct(Points);
    Vect_destroy_cats_struct(Cats);

    return 0;
}

//.........这里部分代码省略.........
	double x, y;
	
	/* Cycle through all lines */
	Vect_rewind(&Map);
	while (1) {
	    type = Vect_read_next_line(&Map, Points, Cats);	/* read line */
	    if (type == 0)
		continue;		/* Dead */

	    if (type == -1)
		G_fatal_error(_("Reading input vector map"));
	    if (type == -2)
		break;
		
	    if (first && Points->n_points > 0) {
		first = 0;
		src_box.E = src_box.W = Points->x[0];
		src_box.N = src_box.S = Points->y[0];
		src_box.T = src_box.B = Points->z[0];
	    }
	    for (i = 0; i < Points->n_points; i++) {
		if (src_box.E < Points->x[i])
		    src_box.E = Points->x[i];
		if (src_box.W > Points->x[i])
		    src_box.W = Points->x[i];
		if (src_box.N < Points->y[i])
		    src_box.N = Points->y[i];
		if (src_box.S > Points->y[i])
		    src_box.S = Points->y[i];
	    }
	    counter++;
	}
	if (counter == 0) {
	    G_warning(_("Input vector map <%s> is empty"), omap_name);
	    exit(EXIT_SUCCESS);
	}
	/* NW corner */
	x = src_box.W;
	y = src_box.N;
	if (pj_do_transform(1, &x, &y, NULL,
			    &info_in, &info_out) < 0) {
	    G_fatal_error(_("Error in pj_do_transform"));
	}
	tgt_box.E = x;
	tgt_box.W = x;
	tgt_box.N = y;
	tgt_box.S = y;
	/* SW corner */
	x = src_box.W;
	y = src_box.S;
	if (pj_do_transform(1, &x, &y, NULL,
			    &info_in, &info_out) < 0) {
	    G_fatal_error(_("Error in pj_do_transform"));
	}
	if (tgt_box.W > x)
	    tgt_box.W = x;
	if (tgt_box.E < x)
	    tgt_box.E = x;
	if (tgt_box.N < y)
	    tgt_box.N = y;
	if (tgt_box.S > y)
	    tgt_box.S = y;
	/* NE corner */
	x = src_box.E;
	y = src_box.N;
	if (pj_do_transform(1, &x, &y, NULL,

/* ************************************************************************* */
int test_geom_data(void)
{
    struct Cell_head region2d;
    G3D_Region region3d;
    N_geom_data *geom = NULL;
    int sum = 0, i;
    double area = 0;

    G_get_set_window(&region2d);

    /*Set the defaults */
    G3d_initDefaults();

    /*get the current region */
    G3d_getWindow(&region3d);

    geom = N_alloc_geom_data();
    if (!geom) {
	G_warning("error in N_alloc_geom_data");
	return 1;
    }
    N_free_geom_data(geom);
    geom = NULL;

    /* ************ 2d region *************** */
    geom = N_init_geom_data_2d(&region2d, geom);
    if (!geom) {
	G_warning("error in N_init_geom_data_2d");
	return 2;
    }

    geom = N_init_geom_data_2d(&region2d, geom);
    if (!geom) {
	G_warning("error in N_init_geom_data_2d");
	return 3;
    }

    if (geom->dim != 2)
	sum++;
    if (geom->planimetric == 0 && geom->area == NULL)
	sum++;
    if (geom->planimetric == 1 && geom->area != NULL)
	sum++;

    /*get areas */
    area = 0.0;
    if (geom->planimetric == 0) {
	for (i = 0; i < geom->rows; i++)
	    area += N_get_geom_data_area_of_cell(geom, i);

	if (area == 0) {
	    G_warning("Wrong area calculation in N_init_geom_data_2d");
	    sum++;
	}
    }

    area = 0.0;
    if (geom->planimetric == 1) {
	for (i = 0; i < geom->rows; i++)
	    area += N_get_geom_data_area_of_cell(geom, i);

	if (area == 0) {
	    G_warning
		("Wrong area calculation in N_get_geom_data_area_of_cell");
	    sum++;
	}
    }


    N_free_geom_data(geom);
    geom = NULL;

    /* ************ 3d region *************** */
    geom = N_init_geom_data_3d(&region3d, geom);
    if (!geom) {
	G_warning("error in N_init_geom_data_3d");
	return 2;
    }

    geom = N_init_geom_data_3d(&region3d, geom);
    if (!geom) {
	G_warning("error in N_init_geom_data_3d");
	return 3;
    }

    if (geom->dim != 3)
	sum++;
    if (geom->planimetric == 0 && geom->area == NULL)
	sum++;

    if (geom->planimetric == 1 && geom->area != NULL)
	sum++;

    /*get areas */
    area = 0.0;
    if (geom->planimetric == 0) {
	for (i = 0; i < geom->rows; i++)
	    area += N_get_geom_data_area_of_cell(geom, i);

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

/*!
   \brief Find area outside island

   \param Map vector map
   \param isle isle id
   \param box isle bbox

   \return area id
   \return 0 if not found
 */
int Vect_isle_find_area(struct Map_info *Map, int isle, const struct bound_box *box)
{
    int i, j, line, sel_area, area, poly;
    const struct Plus_head *plus;
    struct P_line *Line;
    struct P_node *Node;
    struct P_isle *Isle;
    struct P_area *Area;
    struct P_topo_b *topo;
    struct bound_box *abox, nbox;
    static struct boxlist *List = NULL;
    static BOX_SIZE *size_list;
    static int alloc_size_list = 0;

    /* see also Vect_find_area() */

    /* Note: We should check all isle points (at least) because if topology is not clean
     * and two areas overlap, isle which is not completely within area may be attached,
     * but it would take long time */

    G_debug(3, "Vect_isle_find_area () island = %d", isle);
    plus = &(Map->plus);

    if (plus->Isle[isle] == NULL) {
	G_warning(_("Request to find area outside nonexistent isle"));
	return 0;
    }

    if (!List) {
	List = Vect_new_boxlist(1);
	alloc_size_list = 10;
	size_list = G_malloc(alloc_size_list * sizeof(BOX_SIZE));
    }

    Isle = plus->Isle[isle];
    line = abs(Isle->lines[0]);
    Line = plus->Line[line];
    topo = (struct P_topo_b *)Line->topo;
    Node = plus->Node[topo->N1];

    /* select areas by box */
    nbox.E = Node->x;
    nbox.W = Node->x;
    nbox.N = Node->y;
    nbox.S = Node->y;
    nbox.T = PORT_DOUBLE_MAX;
    nbox.B = -PORT_DOUBLE_MAX;
    Vect_select_areas_by_box(Map, &nbox, List);
    G_debug(3, "%d areas overlap island boundary point", List->n_values);

    /* sort areas by bbox size
     * get the smallest area that contains the isle
     * using the bbox size is working because if 2 areas both contain
     * the isle, one of these areas must be inside the other area
     * which means that the bbox of the outer area must be larger than
     * the bbox of the inner area, and equal bbox sizes are not possible */

    if (alloc_size_list < List->n_values) {
	alloc_size_list = List->n_values;
	size_list = G_realloc(size_list, alloc_size_list * sizeof(BOX_SIZE));
    }

    j = 0;
    for (i = 0; i < List->n_values; i++) {
	abox = &List->box[i];

	if (box->E > abox->E || box->W < abox->W || box->N > abox->N ||
	    box->S < abox->S) {
	    G_debug(3, "  isle not completely inside area box");
	    continue;
	}
	
	List->id[j] = List->id[i];
	List->box[j] = List->box[i];
	size_list[j].i = List->id[j];
	size_list[j].box = List->box[j];
	size_list[j].size = (abox->N - abox->S) * (abox->E - abox->W);
	j++;
    }
    List->n_values = j;

    if (List->n_values > 1) {
	if (List->n_values == 2) {
	    /* simple swap */
	    if (size_list[1].size < size_list[0].size) {
		size_list[0].i = List->id[1];
		size_list[1].i = List->id[0];
		size_list[0].box = List->box[1];
		size_list[1].box = List->box[0];
	    }
//.........这里部分代码省略.........

/* code taken from Vect_build_nat() */
int level_one_info(struct Map_info *Map)
{
    struct Plus_head *plus;
    int i, type, first = 1;
    off_t offset;
    struct line_pnts *Points;
    struct line_cats *Cats;
    struct bound_box box;

    int n_primitives, n_points, n_lines, n_boundaries, n_centroids;
    int n_faces, n_kernels;

    G_debug(1, "Count vector objects for level 1");

    plus = &(Map->plus);

    n_primitives = n_points = n_lines = n_boundaries = n_centroids = 0;
    n_faces = n_kernels = 0;

    Points = Vect_new_line_struct();
    Cats = Vect_new_cats_struct();
    
    Vect_rewind(Map);
    /* G_message(_("Registering primitives...")); */
    i = 1;
    while (1) {
	/* register line */
	type = Vect_read_next_line(Map, Points, Cats);

	/* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
	if (type == -1) {
	    G_warning(_("Unable to read vector map"));
	    return 0;
	}
	else if (type == -2) {
	    break;
	}

	/* count features */
	n_primitives++;
	
	if (type & GV_POINT)  /* probably most common */
	    n_points++;
	else if (type & GV_LINE)
	    n_lines++;
	else if (type & GV_BOUNDARY)
	    n_boundaries++;
	else if (type & GV_CENTROID)
	    n_centroids++;
	else if (type & GV_KERNEL)
	    n_kernels++;
	else if (type & GV_FACE)
	    n_faces++;

	offset = Map->head.last_offset;

	G_debug(3, "Register line: offset = %lu", (unsigned long)offset);
	dig_line_box(Points, &box);
	if (first == 1) {
	    Vect_box_copy(&(plus->box), &box);
	    first = 0;
	}
	else
	    Vect_box_extend(&(plus->box), &box);

	/* can't print progress, unfortunately */
/*
	if (G_verbose() > G_verbose_min() && i % 1000 == 0) {
	    if (format == G_INFO_FORMAT_PLAIN)
		fprintf(stderr, "%d..", i);
	    else
		fprintf(stderr, "%11d\b\b\b\b\b\b\b\b\b\b\b", i);
	}
	i++;
*/
    }

    /* save result in plus */
    plus->n_lines = n_primitives;
    plus->n_plines = n_points;
    plus->n_llines = n_lines;
    plus->n_blines = n_boundaries;
    plus->n_clines = n_centroids;
    plus->n_klines = n_kernels;
    plus->n_flines = n_faces;

    return 1;
}

/*!
   \brief Extensive tests for correct topology

   - lines or boundaries of zero length
   - intersecting boundaries, ie. overlapping areas
   - areas without centroids that are not isles

   \param Map vector map
   \param[out] Err vector map where errors will be written or NULL

   \return 1 on success
   \return 0 on error
 */
int Vect_t