自学教程：C++ G_warning函数代码示例

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

voidsubcluster(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 */intplot(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_topo_check(struct Map_info *Map, struct Map_info *Err){    int line, nlines;    int nerrors, n_zero_lines, n_zero_boundaries;    struct line_pnts *Points;    struct line_cats *Cats;    /* rebuild topology if needed */    if (Vect_get_built(Map) != GV_BUILD_ALL) {	Vect_build_partial(Map, GV_BUILD_NONE);	Vect_build(Map);    }    G_message(_("Checking for topological errors..."));    Points = Vect_new_line_struct();    Cats = Vect_new_cats_struct();    /* lines or boundaries of zero length */    n_zero_lines = n_zero_boundaries = 0;    nlines = Vect_get_num_lines(Map);    for (line = 1; line <= nlines; line++) {	int type;	if (!Vect_line_alive(Map, line))	    continue;	    	type = Vect_get_line_type(Map, line);	if (type & GV_LINES) {	    double len;	    	    Vect_read_line(Map, Points, Cats, line);	    len = Vect_line_length(Points);	    	    if (len == 0) {		if (type & GV_LINE)		    n_zero_lines++;		else if (type & GV_BOUNDARY)		    n_zero_boundaries++;		    		if (Err)		    Vect_write_line(Err, type, Points, Cats);	    }	}    }        if (n_zero_lines)	G_warning(_("Number of lines of length zero: %d"), n_zero_lines);    if (n_zero_boundaries)	G_warning(_("Number of boundaries of length zero: %d"), n_zero_boundaries);    /* remaining checks are for areas only */    if (Vect_get_num_primitives(Map, GV_BOUNDARY) == 0)	return 1;    /* intersecting boundaries -> overlapping areas */    nerrors = Vect_check_line_breaks(Map, GV_BOUNDARY, Err);    if (nerrors)	G_warning(_("Number of boundary intersections: %d"), nerrors);    /* areas without centroids that are not isles     * only makes sense if all boundaries are correct */    nerrors = 0;    for (line = 1; line <= nlines; line++) {	int type;		if (!Vect_line_alive(Map, line))	    continue;	    	type = Vect_get_line_type(Map, line);	if (type == GV_BOUNDARY) {	    struct P_topo_b *topo = (struct P_topo_b *)Map->plus.Line[line]->topo;	    if (topo->left == 0 || topo->right == 0) {		G_debug(3, "line = %d left = %d right = %d", line, 			topo->left, topo->right);		nerrors++;	    }	}    }    if (nerrors)	G_warning(_("Skipping further checks because of incorrect boundaries"));    else {	int i, area, left, right, neighbour;	int nareas = Vect_get_num_areas(Map);//.........这里部分代码省略.........

//.........这里部分代码省略.........				G_fatal_error("avl_make unknown error");				return RLI_ERRORE;			    }			}		    }		    totCorr = 1;		}		/* endif not equal fcells */		else {		/*equal fcells */		    totCorr++;		}		precCell = corrCell;	    }	}    }    /*last closing */    if (a != 0) {	if (albero == NULL) {	    uc.val.fc = precCell;	    albero = avl_make(uc, totCorr);	    if (albero == NULL) {		G_fatal_error("avl_make error");		return RLI_ERRORE;	    }	    m++;	}	else {	    uc.val.fc = precCell;	    ris = avl_add(&albero, uc, totCorr);	    switch (ris) {	    case AVL_ERR:		{		    G_fatal_error("avl_add error");		    return RLI_ERRORE;		}	    case AVL_ADD:		{		    m++;		    break;		}	    case AVL_PRES:		{		    break;		}	    default:		{		    G_fatal_error("avl_add unknown error");		    return RLI_ERRORE;		}	    }	}    }    array = G_malloc(m * sizeof(AVL_tableRow));    if (array == NULL) {	G_fatal_error("malloc array failed");	return RLI_ERRORE;    }    tot = avl_to_array(albero, zero, array);    if (tot != m) {	G_warning("avl_to_array unaspected value. the result could be wrong");	return RLI_ERRORE;    }    char *sval;    sval = par[0];    double alpha_double;    alpha_double = (double)atof(sval);    /* calculate index summary */    for (i = 0; i < m; i++) {	t = (double)(array[i]->tot);	p = t / area;	G_debug(1, "Valore p: %g, valore pow: %g", p, pow(p, alpha_double));	somma = somma + pow(p, alpha_double);    }    indice = (1 / (1 - alpha_double)) * log(somma);    if (isnan(indice) || isinf(indice)) {	indice = -1;    }    G_debug(1, "Valore somma: %g Valore indice: %g", somma, indice);    *result = indice;    G_free(array);    if (masked)	G_free(mask_buf);    return RLI_OK;}

/*!   /brief Build area on given side of line (GV_LEFT or GV_RIGHT)   /param Map pointer to Map_info structure   /param iline line id   /param side side (GV_LEFT or GV_RIGHT)   /return > 0 area id   /return < 0 isle id   /return 0 not created (may also already exist) */int Vect_build_line_area(struct Map_info *Map, int iline, int side){    int area, isle, n_lines;    struct Plus_head *plus;    struct bound_box box;    static struct line_pnts *APoints = NULL;    plus_t *lines;    double area_size;    plus = &(Map->plus);    G_debug(3, "Vect_build_line_area() line = %d, side = %d", iline, side);    if (!APoints)	APoints = Vect_new_line_struct();        /* get area */    area = dig_line_get_area(plus, iline, side);    if (area != 0) {        /* -> there is already an area on this side of the line, skip */        G_debug(3, "  area/isle = %d -> skip", area);        return 0;    }        /* build an area with this line */    n_lines = dig_build_area_with_line(plus, iline, side, &lines);    G_debug(3, "  n_lines = %d", n_lines);    if (n_lines < 1) {	return 0;    }				/* area was not built */    /* get line points which forms a boundary of an area */    Vect__get_area_points(Map, lines, n_lines, APoints);    dig_line_box(APoints, &box);    Vect_line_prune(APoints);    if (APoints->n_points < 4) {	G_warning(_("Area of size = 0.0 (less than 4 vertices) ignored"));	return 0;    }    /* Area or island ? */    dig_find_area_poly(APoints, &area_size);    /* area_size = dig_find_poly_orientation(APoints); */    /* area_size is not real area size, we are only interested in the sign */    G_debug(3, "  area/isle size = %f", area_size);    if (area_size > 0) {	/* CW: area */	/* add area structure to plus */	area = dig_add_area(plus, n_lines, lines, &box);	if (area == -1) {	/* error */	    G_fatal_error(_("Unable to add area (map closed, topo saved)"));	}	G_debug(3, "  -> area %d", area);	return area;    }    else if (area_size < 0) {	/* CCW: island */	isle = dig_add_isle(plus, n_lines, lines, &box);	if (isle == -1) {	/* error */	    G_fatal_error(_("Unable to add isle (map closed, topo saved)"));	}	G_debug(3, "  -> isle %d", isle);	return -isle;    }    else {	/* TODO: What to do with such areas? Should be areas/isles of size 0 stored,	 *        so that may be found and cleaned by some utility	 *  Note: it would be useful for vertical closed polygons, but such would be added twice	 *        as area */	G_warning(_("Area of size = 0.0 ignored"));    }    return 0;}

static int reestimate(struct ClassSig *Sig, int nbands){    int i;    int s;    int b1, b2;    int singular;    double pi_sum;    double diff1, diff2;    struct ClassData *Data;    G_debug(2, "reestimate()");    /* set data pointer */    Data = &(Sig->ClassData);    pi_sum = 0;    for (i = 0; i < Sig->nsubclasses; i++) {	/* Compute N */	Sig->SubSig[i].N = 0;	for (s = 0; s < Data->npixels; s++)	    Sig->SubSig[i].N += Data->p[s][i];	Sig->SubSig[i].pi = Sig->SubSig[i].N;	/* Compute means and variances for each subcluster, */	/* and remove small clusters.                       */	/* For large subclusters */	if (Sig->SubSig[i].N > SMALLEST_SUBCLUST) {	    /* Compute mean */	    for (b1 = 0; b1 < nbands; b1++) {		Sig->SubSig[i].means[b1] = 0;		for (s = 0; s < Data->npixels; s++)		    if (!Rast_is_d_null_value(&Data->x[s][b1]))			Sig->SubSig[i].means[b1] +=			    Data->p[s][i] * Data->x[s][b1];		Sig->SubSig[i].means[b1] /= (Sig->SubSig[i].N);		/* Compute R */		for (b2 = 0; b2 <= b1; b2++) {		    Sig->SubSig[i].R[b1][b2] = 0;		    for (s = 0; s < Data->npixels; s++) {			if (!Rast_is_d_null_value(&Data->x[s][b1])			    && !Rast_is_d_null_value(&Data->x[s][b2])) {			    diff1 = Data->x[s][b1] - Sig->SubSig[i].means[b1];			    diff2 = Data->x[s][b2] - Sig->SubSig[i].means[b2];			    Sig->SubSig[i].R[b1][b2] +=				Data->p[s][i] * diff1 * diff2;			}		    }		    Sig->SubSig[i].R[b1][b2] /= (Sig->SubSig[i].N);		    Sig->SubSig[i].R[b2][b1] = Sig->SubSig[i].R[b1][b2];		}	    }	}	/* For small subclusters */	else {	    G_warning(_("Subsignature %d only contains %.0f pixels"),		      i, Sig->SubSig[i].N);	    	    Sig->SubSig[i].pi = 0;	    for (b1 = 0; b1 < nbands; b1++) {		Sig->SubSig[i].means[b1] = 0;		for (b2 = 0; b2 < nbands; b2++)		    Sig->SubSig[i].R[b1][b2] = 0;	    }	}	pi_sum += Sig->SubSig[i].pi;    }    /* Normalize probabilities for subclusters */    if (pi_sum > 0) {	for (i = 0; i < Sig->nsubclasses; i++)	    Sig->SubSig[i].pi /= pi_sum;    }    else {	for (i = 0; i < Sig->nsubclasses; i++)	    Sig->SubSig[i].pi = 0;    }    /* Compute constants and reestimate if any singular subclusters occur */    singular = compute_constants(Sig, nbands);    return (singular);}

//.........这里部分代码省略.........	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);	if (plus->n_klines > 0)	    G_message(_("Number of kernels: %d"), plus->n_klines);    }    if (plus->built >= GV_BUILD_AREAS) {	int line, nlines, area, nareas, err_boundaries, err_centr_out,	    err_centr_dupl, err_nocentr;	struct P_line *Line;	struct Plus_head *Plus;	/* Count errors (it does not take much time comparing to build process) */	Plus = &(Map->plus);	nlines = Vect_get_num_lines(Map);	err_boundaries = err_centr_out = err_centr_dupl = 0;	for (line = 1; line <= nlines; line++) {	    Line = Plus->Line[line];	    if (!Line)		continue;	    if (Line->type == GV_BOUNDARY) {		struct P_topo_b *topo = (struct P_topo_b *)Line->topo;		if (topo->left == 0 || topo->right == 0) {		    G_debug(3, "line = %d left = %d right = %d", line, 			    topo->left, topo->right);		    err_boundaries++;		}	    }	    if (Line->type == GV_CENTROID) {		struct P_topo_c *topo = (struct P_topo_c *)Line->topo;		if (topo->area == 0)		    err_centr_out++;		else if (topo->area < 0)		    err_centr_dupl++;	    }	}	err_nocentr = 0;	nareas = Vect_get_num_areas(Map);	for (area = 1; area <= nareas; area++) {	    if (!Vect_area_alive(Map, area))		continue;	    line = Vect_get_area_centroid(Map, area);	    if (line == 0)		err_nocentr++;	}	G_message(_("Number of areas: %d"), plus->n_areas);	G_message(_("Number of isles: %d"), plus->n_isles);#if 0	/* not an error, message disabled to avoid confusion */	if (err_nocentr)	    G_message(_("Number of areas without centroid: %d"),		      err_nocentr);#endif	if (plus->n_clines > plus->n_areas)	    G_warning(_("Number of centroids exceeds number of areas: %d > %d"),		      plus->n_clines, plus->n_areas);	if (err_boundaries)	    G_warning(_("Number of incorrect boundaries: %d"),		      err_boundaries);	if (err_centr_out)	    G_warning(_("Number of centroids outside area: %d"),		      err_centr_out);	if (err_centr_dupl)	    G_warning(_("Number of duplicate centroids: %d"), err_centr_dupl);    }    else if (build > GV_BUILD_NONE) {	G_message(_("Number of areas: -"));	G_message(_("Number of isles: -"));    }    return 1;}

//.........这里部分代码省略.........    if(num_raster_inputs > 0)        wps_print_literal_input_output(WPS_INPUT, 0, 1, "grass_band_number", "Band to select for processing (default is all bands)",            "This parameter defines band number of the input raster files which should be processed. As default all bands are processed and used as single and multiple inputs for raster modules.",            "integer", 0, NULL, 0, NULL, TYPE_OTHER);    /* End of inputs */    wps_print_data_inputs_end();    /* Start of the outputs */    wps_print_process_outputs_begin();    found_output = 0;    /*parse the output. only raster maps, vector maps, space time raster and vector datasets plus stdout are supported */    if (st->n_opts) {	opt = &st->first_option;	while (opt != NULL) {            identifier = NULL;            title = NULL;            abstract = NULL;            value = NULL;            is_output = 0;            data_type = TYPE_OTHER;	    if (opt->gisprompt) {		const char *atts[] = { "age", "element", "prompt", NULL };		top = G_calloc(strlen(opt->gisprompt) + 1, 1);		strcpy(top, opt->gisprompt);		s = strtok(top, ",");		for (i = 0; s != NULL && atts[i] != NULL; i++) {                    char *token = G_store(s);                    /* we print only the output parameter */                    if(strcmp(token, "new") == 0)                        is_output = 1;                    if(strcmp(token, "raster") == 0)                    {                        data_type = TYPE_RASTER;                    }                    if(strcmp(token, "vector") == 0)                    {                        data_type = TYPE_VECTOR;                    }                    if(strcmp(token, "stds") == 0)                    {                        data_type = TYPE_STDS;                    }                    if(strcmp(token, "strds") == 0)                    {                        data_type = TYPE_STRDS;                    }                    if(strcmp(token, "stvds") == 0)                    {                        data_type = TYPE_STVDS;                    }                    if(strcmp(token, "file") == 0)                    {                        data_type = TYPE_PLAIN_TEXT;                    }                    s = strtok(NULL, ",");                    G_free(token);		}		G_free(top);	    }            /* Only single module output is supported!! */            if(is_output == 1)            {		if(opt->multiple == YES)		    G_warning(_("Multiple outputs are not supported by WPS 1.0.0"));                identifier = opt->key;                 if(opt->label) {                    title = opt->label;		}                if (opt->description) {		    if(!opt->label)			title = opt->description;		    else			abstract = opt->description;		}                if(data_type == TYPE_RASTER || data_type == TYPE_VECTOR || 		   data_type == TYPE_STRDS  || data_type == TYPE_STVDS || 		   data_type == TYPE_STDS  || data_type == TYPE_PLAIN_TEXT) {                    wps_print_complex_output(identifier, title, abstract, data_type);                    found_output = 1;                }            }	    opt = opt->next_opt;	}        /* we assume the computatuon output on stdout, if no raster/vector output was found*/        if(found_output == 0)            wps_print_complex_output("stdout", "Module output on stdout", "The output of the module written to stdout", TYPE_PLAIN_TEXT);    }    wps_print_process_outputs_end();    wps_print_process_description_end();    wps_print_process_descriptions_end();}

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

//.........这里部分代码省略.........    bxma = 2093461. * conv;    byma = 731529. * conv;    bresx = 2. * conv;    bresy = 2. * conv;    maxwab = 100000;    mx2o = (int)((bxma - bxmi) / bresx);    my2o = (int)((byma - bymi) / bresy);    /* relative small box coordinates: leave 1 grid layer for overlap */    bxmi = bxmi - mixx + stepx;    bymi = bymi - miyy + stepy;    bxma = bxma - mixx - stepx;    byma = byma - miyy - stepy;    mx2 = mx2o - 2 * ((int)(stepx / bresx));    my2 = my2o - 2 * ((int)(stepy / bresy));#endif    elevin = parm.elevin->answer;    wdepth = parm.wdepth->answer;    dxin = parm.dxin->answer;    dyin = parm.dyin->answer;    detin = parm.detin->answer;    tranin = parm.tranin->answer;    tauin = parm.tauin->answer;    manin = parm.manin->answer;    tc = parm.tc->answer;    et = parm.et->answer;    conc = parm.conc->answer;    flux = parm.flux->answer;    erdep = parm.erdep->answer;    outwalk = parm.outwalk->answer;     /*      sscanf(parm.nwalk->answer, "%d", &maxwa); */    sscanf(parm.niter->answer, "%d", &timesec);    sscanf(parm.outiter->answer, "%d", &iterout);/*    sscanf(parm.density->answer, "%d", &ldemo); */    sscanf(parm.diffc->answer, "%lf", &frac);    sscanf(parm.maninval->answer, "%lf", &manin_val);    /* Recompute timesec from user input in minutes     * to real timesec in seconds */    timesec = timesec * 60.0;    iterout = iterout * 60.0;    if ((timesec / iterout) > 100.0)	G_message(_("More than 100 files are going to be created !!!!!"));    /* compute how big the raster is and set this to appr 2 walkers per cell */    if (parm.nwalk->answer == NULL) {	maxwa = mx * my * 2;	rwalk = (double)(mx * my * 2.);	G_message(_("default nwalk=%d, rwalk=%f"), maxwa, rwalk);    }    else {	sscanf(parm.nwalk->answer, "%d", &maxwa);	rwalk = (double)maxwa;    }    /*rwalk = (double) maxwa; */    if (conv != 1.0)	G_message(_("Using metric conversion factor %f, step=%f"), conv,		  step);    if ((tc == NULL) && (et == NULL) && (conc == NULL) && (flux == NULL) &&	(erdep == NULL))	G_warning(_("You are not outputting any raster or site files"));    ret_val = input_data();    if (ret_val != 1)	G_fatal_error(_("Input failed"));    /* mandatory for si,sigma */    si = G_alloc_matrix(my, mx);    sigma = G_alloc_matrix(my, mx);    /* memory allocation for output grids */    dif = G_alloc_fmatrix(my, mx);    if (erdep != NULL || et != NULL)	er = G_alloc_fmatrix(my, mx);    seeds(rand1, rand2);    grad_check();    if (et != NULL)	erod(si);    /* treba dat output pre topoerdep */    main_loop();    if (tserie == NULL) {	ii = output_data(0, 1.);	if (ii != 1)	    G_fatal_error(_("Cannot write raster maps"));    }    /* Exit with Success */    exit(EXIT_SUCCESS);}

//.........这里部分代码省略.........	}	if (KEY("lpos")) {	    if (sscanf(data, "%d", &itmp) < 1 || itmp < 0) {		itmp = -1;		error(key, data, "illegal lpos");		continue;	    }	    vector.layer[vec].lpos = itmp;	    continue;	}	if (KEY("symbol")) {	    /* TODO: test here if isymbol exists */	    vector.layer[vec].symbol = G_store(data);	    continue;	}	if (KEY("eps")) {	    char *cc;	    G_chop(data);	    vector.layer[vec].epspre = G_store(data);	    /* epstype: 0 - no eps, 1 - common eps, 2 - eps for each category */	    vector.layer[vec].epstype = 1;	    /* find dynamic filename by cat number character */	    /* pre is filename before the $, suf is filename after the $ */	    cc = (char *)strchr(vector.layer[vec].epspre, '$');	    if (cc != NULL) {		*cc = '/0';		vector.layer[vec].epssuf = G_store(cc + sizeof(char));		vector.layer[vec].epstype = 2;		G_debug(2, "epstype=%d, pre=[%s], suf=[%s]",			vector.layer[vec].epstype, vector.layer[vec].epspre,			vector.layer[vec].epssuf);	    }	    else {		G_debug(2, "epstype=%d, eps file=[%s]",			vector.layer[vec].epstype, vector.layer[vec].epspre);	    }	    continue;	}	if (KEY("size")) {	    if (sscanf(data, "%lf", &size) != 1 || size <= 0.0) {		size = 1.0;		error(key, data, "illegal size request (vpoints)");	    }	    vector.layer[vec].size = size;	    continue;	}	/* 	   GRASS 6.3: sizecol renamed to sizecolumn	   remove sizecol test and the warning in GRASS7	 */	if (KEY("sizecol")) {	    G_warning(_("The mapping instruction <%s> will be renamed to <%s> "		       "in future versions of GRASS. Please use <%s> instead."),		      "sizecol", "sizecolumn", "sizecolumn");	}	if (KEY("sizecol") || KEY("sizecolumn")) {	    G_strip(data);	    vector.layer[vec].sizecol = G_store(data);	    continue;	}	if (KEY("scale")) {	    if (sscanf(data, "%lf", &scale) != 1 || scale <= 0.0) {		scale = 1.0;		error(key, data, "illegal scale request (vpoints)");	    }	    vector.layer[vec].scale = scale;	    continue;	}	if (KEY("rotate")) {	    if (sscanf(data, "%lf", &rotate) != 1) {		rotate = 0.0;		error(key, data, "illegal rotation request (vpoints)");	    }	    vector.layer[vec].rotate = rotate;	    continue;	}	if (KEY("rotatecolumn")) {	    G_strip(data);	    vector.layer[vec].rotcol = G_store(data);	    continue;	}	error(key, "", "illegal request (vpoints)");    }    vector.count++;    return 1;}

int main(int argc, char **argv){    char title[80];    char buf[80], *p;    struct Ortho_Image_Group group;    struct GModule *module;    struct Option *group_opt;    /* must run in a term window */    G_putenv("GRASS_UI_TERM", "1");    /* initialize grass */    G_gisinit(argv[0]);    module = G_define_module();    module->keywords = _("imagery, orthorectify");    module->description = _("Menu driver for the photo imagery programs.");    group_opt = G_define_standard_option(G_OPT_I_GROUP);    group_opt->description =	_("Name of imagery group for ortho-rectification");    if (G_parser(argc, argv))	exit(EXIT_FAILURE);    strncpy(group.name, group_opt->answer, 99);    group.name[99] = '/0';    /* strip off mapset if it's there: I_() fns only work with current mapset */    if ((p = strchr(group.name, '@')))	*p = 0;    /* get and check the group reference files */    if (!I_get_group_ref(group.name, &group.group_ref)) {	G_warning(_("Pre-selected group <%s> not found"), group.name);	/* clean the wrong name in GROUPFILE */	I_put_group("");	/* ask for new group name */	if (!I_ask_group_old(	    _("Enter imagery group for ortho-rectification"), group.name))	    exit(EXIT_SUCCESS);	I_get_group_ref(group.name, &group.group_ref);    }    if (group.group_ref.nfiles <= 0)	G_fatal_error(_("Group [%s] contains no files"), group.name);    I_put_group(group.name);    while (1) {	if (!I_get_group(group.name)) {	    exit(EXIT_SUCCESS);	}	/* print the screen full of options */	sprintf(title, "i.ortho.photo -- /tImagery Group = %s ", group.name);	G_clear_screen();	fprintf(stderr, "%s/n/n", title);	fprintf(stderr, "Initialization Options:/n");	fprintf(stderr, "/n");	fprintf(stderr, "   1.     Select/Modify imagery group/n");	fprintf(stderr, "   2.     Select/Modify imagery group target/n");	fprintf(stderr, "   3.     Select/Modify target elevation model/n");	fprintf(stderr, "   4.     Select/Modify imagery group camera/n");	fprintf(stderr, "/n");	fprintf(stderr, "Transformation Parameter Computations:/n");	fprintf(stderr, "/n");	fprintf(stderr, "   5.     Compute image-to-photo transformation/n");	fprintf(stderr, "   6.     Initialize exposure station parameters/n");	fprintf(stderr, "   7.     Compute ortho-rectification parameters/n");	fprintf(stderr, "/n");	fprintf(stderr, "Ortho-rectification Option:/n");	fprintf(stderr, "/n");	fprintf(stderr, "   8.     Ortho-rectify imagery files/n");	fprintf(stderr, "/n");	fprintf(stderr, "RETURN   exit/n");	fprintf(stderr, "/n> ");	/* Get the option */	if (!G_gets(buf))	    continue;	if (*buf == 0)		/* exit */	    exit(EXIT_SUCCESS);	/* run the program chosen */	G_strip(buf);	fprintf(stderr, "<%s>/n", buf);	if (strcmp(buf, "1") == 0)	    run_system("i.group");	if (strcmp(buf, "2") == 0)	    run_etc_imagery("i.photo.target", group.name);	if (strcmp(buf, "3") == 0)	    run_etc_imagery("i.photo.elev", group.name);	if (strcmp(buf, "4") == 0)	    run_etc_imagery("i.photo.camera", group.name);	if (strcmp(buf, "5") == 0)	    run_etc_imagery("i.photo.2image", group.name);	if (strcmp(buf, "6") == 0)//.........这里部分代码省略.........