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

FILE *openAscii(char *asciiFile, RASTER3D_Region * region){    FILE *fp;    double tmp;    char buff[1024];    char line_buff[1024];    G_debug(3, "openAscii: opens the ascii file and reads the header");    fp = fopen(asciiFile, "r");    if (fp == NULL) {        perror(asciiFile);        G_usage();        exit(EXIT_FAILURE);    }    /* Initialize the default order */    rowOrder = ROW_ORDER_NORTH_TO_SOUTH;    depthOrder = DEPTH_ORDER_BOTTOM_TO_TOP;    /* Read the first line and check for grass version */    G_getl2(line_buff, 1024, fp);    /* First check for new ascii format*/    if (sscanf(line_buff, "version: %s", buff) == 1) {        G_message("Found version information: %s/n", buff);        if (G_strcasecmp(buff, "grass7") == 0) {            /* Parse the row and depth order */            G_getl2(line_buff, 1024, fp);            if (sscanf(line_buff, "order: %s", buff) != 1)                fatalError("Unable to parse the row and depth order");            if (G_strcasecmp(buff, "nsbt") == 0) {                rowOrder = ROW_ORDER_NORTH_TO_SOUTH;                depthOrder = DEPTH_ORDER_BOTTOM_TO_TOP;                G_message("Found north -> south, bottom -> top order (nsbt)");            }            if (G_strcasecmp(buff, "snbt") == 0) {                rowOrder = ROW_ORDER_SOUTH_TO_NORTH;                depthOrder = DEPTH_ORDER_BOTTOM_TO_TOP;                G_message("Found south -> north, bottom -> top order (snbt)");            }            if (G_strcasecmp(buff, "nstb") == 0) {                rowOrder = ROW_ORDER_NORTH_TO_SOUTH;                depthOrder = DEPTH_ORDER_TOP_TO_BOTTOM;                G_message("Found north -> south, top -> bottom order (nstb)");            }            if (G_strcasecmp(buff, "sntb") == 0) {                rowOrder = ROW_ORDER_SOUTH_TO_NORTH;                depthOrder = DEPTH_ORDER_TOP_TO_BOTTOM;                G_message("Found south -> north, top -> bottom order (sntb)");            }        } else {            G_fatal_error(_("Unsupported GRASS version %s"), buff);        }    } else {        /* Rewind the stream if no grass version info found */        rewind(fp);    }    Rast3d_get_window(region);    readHeaderString(fp, "north:", &(region->north));    readHeaderString(fp, "south:", &(region->south));    readHeaderString(fp, "east:", &(region->east));    readHeaderString(fp, "west:", &(region->west));    readHeaderString(fp, "top:", &(region->top));    readHeaderString(fp, "bottom:", &(region->bottom));    readHeaderString(fp, "rows:", &tmp);    region->rows = (int) tmp;    readHeaderString(fp, "cols:", &tmp);    region->cols = (int) tmp;    readHeaderString(fp, "levels:", &tmp);    region->depths = (int) tmp;    return fp;}

/* point_in_buf - test if point px,py is in d buffer of Points ** dalpha is in degrees ** returns:  1 in buffer **           0 not in buffer */static int point_in_buf(struct line_pnts *Points, double px, double py, double da,			double db, double dalpha){    int i, np;    double cx, cy;    double delta, delta_k, k;    double vx, vy, wx, wy, mx, my, nx, ny;    double len, tx, ty, d, da2;    G_debug(3, "point_in_buf()");    dalpha *= PI / 180;		/* convert dalpha from degrees to radians */    np = Points->n_points;    da2 = da * da;    for (i = 0; i < np - 1; i++) {	vx = Points->x[i];	vy = Points->y[i];	wx = Points->x[i + 1];	wy = Points->y[i + 1];	if (da != db) {	    mx = wx - vx;	    my = wy - vy;	    len = LENGTH(mx, my);	    elliptic_tangent(mx / len, my / len, da, db, dalpha, &cx, &cy);	    delta = mx * cy - my * cx;	    delta_k = (px - vx) * cy - (py - vy) * cx;	    k = delta_k / delta;	    /*            G_debug(4, "k = %g, k1 = %g", k, (mx * (px - vx) + my * (py - vy)) / (mx * mx + my * my)); */	    if (k <= 0) {		nx = vx;		ny = vy;	    }	    else if (k >= 1) {		nx = wx;		ny = wy;	    }	    else {		nx = vx + k * mx;		ny = vy + k * my;	    }	    /* inverse transform */	    elliptic_transform(px - nx, py - ny, 1 / da, 1 / db, dalpha, &tx,			       &ty);	    d = dig_distance2_point_to_line(nx + tx, ny + ty, 0, vx, vy, 0,					    wx, wy, 0, 0, NULL, NULL, NULL,					    NULL, NULL);	    /*            G_debug(4, "sqrt(d)*da = %g, len' = %g, olen = %g", sqrt(d)*da, da*LENGTH(tx,ty), LENGTH((px-nx),(py-ny))); */	    if (d <= 1) {		/* G_debug(1, "d=%g", d); */		return 1;	    }	}	else {	    d = dig_distance2_point_to_line(px, py, 0, vx, vy, 0, wx, wy, 0,					    0, NULL, NULL, NULL, NULL, NULL);	    /*            G_debug(4, "sqrt(d)     = %g", sqrt(d)); */	    if (d <= da2) {		return 1;	    }	}    }    return 0;}

int db__driver_open_database(dbHandle * handle){    char *name;    dbConnection default_connection;    MYSQL *res;    db_get_connection(&default_connection);    name = G_store(db_get_handle_dbname(handle));    /* if name is empty use default_connection.databaseName */    if (strlen(name) == 0)	name = default_connection.databaseName;    G_debug(3, "db_driver_open_database() mysql: database definition = '%s'",	    name);    /* Embedded version */    {	char *datadir, *database;	char *server_args[4];	char *buf;	if (!replace_variables(name, &datadir, &database)) {	    db_d_append_error(_("Unable parse MySQL embedded database name"));	    db_d_append_error(mysql_error(connection));	    db_d_report_error();	    return DB_FAILED;	}	server_args[0] = "mesql";	/* this string is not used */	G_asprintf(&buf, "--datadir=%s", datadir);	server_args[1] = buf;	/* With InnoDB it is very slow to close the database */	server_args[2] = "--skip-innodb";	/* OK? */	/* Without --bootstrap it complains about missing 	 * mysql.time_zone_leap_second table */	server_args[3] = "--bootstrap";	/* OK? */	if (mysql_server_init(4, server_args, NULL)) {	    db_d_append_error(_("Cannot initialize MySQL embedded server"));	    db_d_append_error(mysql_error(connection));	    db_d_report_error();	    free(datadir);	    free(database);	    return DB_FAILED;	}	connection = mysql_init(NULL);	mysql_options(connection, MYSQL_OPT_USE_EMBEDDED_CONNECTION, NULL);	res =	    mysql_real_connect(connection, NULL, NULL, NULL, database, 0,			       NULL, 0);	free(datadir);	free(database);	if (res == NULL) {	    db_d_append_error(_("Unable to connect to MySQL embedded server: "));	    db_d_append_error(mysql_error(connection));	    db_d_report_error();	    return DB_FAILED;	}    }    return DB_OK;}

//.........这里部分代码省略.........    buf[1] = GV_TOPO_VER_MINOR;    buf[2] = GV_TOPO_EARLIEST_MAJOR;    buf[3] = GV_TOPO_EARLIEST_MINOR;    buf[4] = ptr->port.byte_order;    if (0 >= dig__fwrite_port_C((char *)buf, 5, fp))	return (-1);    /* determine required offset size from coor file size */    if (ptr->coor_size > (off_t)PORT_LONG_MAX) {	/* can only happen when sizeof(off_t) == 8 */	ptr->off_t_size = 8;    }    else	ptr->off_t_size = 4;    /* add a new field with off_t_size after byte_order? */    /* adjust header size for large files */    if (ptr->off_t_size == 8) {	/* 7 offset values and coor file size: add 8 * 4 */	length += 32;    }    /* bytes 6 - 9 : header size */    if (0 >= dig__fwrite_port_L(&length, 1, fp))	return (0);    /* byte 10 : dimension 2D or 3D */    buf[0] = ptr->with_z;    if (0 >= dig__fwrite_port_C((char *)buf, 1, fp))	return (0);    /* bytes 11 - 58 : bound box */    if (0 >= dig__fwrite_port_D(&(ptr->box.N), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_D(&(ptr->box.S), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_D(&(ptr->box.E), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_D(&(ptr->box.W), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_D(&(ptr->box.T), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_D(&(ptr->box.B), 1, fp))	return (-1);    /* bytes 59 - 86 : number of structures */    if (0 >= dig__fwrite_port_P(&(ptr->n_nodes), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_edges), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_lines), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_areas), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_isles), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_volumes), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_holes), 1, fp))	return (-1);    /* bytes 87 - 110 : number of line types */    if (0 >= dig__fwrite_port_P(&(ptr->n_plines), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_llines), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_blines), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_clines), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_flines), 1, fp))	return (-1);    if (0 >= dig__fwrite_port_P(&(ptr->n_klines), 1, fp))	return (-1);    /* bytes 111 - 138 : Offset */    if (0 >= dig__fwrite_port_O(&(ptr->Node_offset), 1, fp, ptr->off_t_size))	return (-1);    if (0 >= dig__fwrite_port_O(&(ptr->Edge_offset), 1, fp, ptr->off_t_size))	return (-1);    if (0 >= dig__fwrite_port_O(&(ptr->Line_offset), 1, fp, ptr->off_t_size))	return (-1);    if (0 >= dig__fwrite_port_O(&(ptr->Area_offset), 1, fp, ptr->off_t_size))	return (-1);    if (0 >= dig__fwrite_port_O(&(ptr->Isle_offset), 1, fp, ptr->off_t_size))	return (-1);    if (0 >= dig__fwrite_port_O(&(ptr->Volume_offset), 1, fp, ptr->off_t_size))	return (-1);    if (0 >= dig__fwrite_port_O(&(ptr->Hole_offset), 1, fp, ptr->off_t_size))	return (-1);    /* bytes 139 - 142 : Coor size and time */    if (0 >= dig__fwrite_port_O(&(ptr->coor_size), 1, fp, ptr->off_t_size))	return (-1);    G_debug(2, "topo body offset %"PRI_OFF_T, dig_ftell(fp));    return (0);}

/* input line must be looped */static void convolution_line(struct line_pnts *Points, double da, double db,			     double dalpha, int side, int round, int caps,			     double tol, struct line_pnts *nPoints){    int i, j, res, np;    double *x, *y;    double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;    double vx1, vy1, wx1, wy1;    double a0, b0, c0, a1, b1, c1;    double phi1, phi2, delta_phi;    double nsegments, angular_tol, angular_step;    double angle0, angle1;    int inner_corner, turns360;    G_debug(3, "convolution_line() side = %d", side);    np = Points->n_points;    x = Points->x;    y = Points->y;    if ((np == 0) || (np == 1))	return;    if ((x[0] != x[np - 1]) || (y[0] != y[np - 1])) {	G_fatal_error(_("Line is not looped"));	return;    }    Vect_reset_line(nPoints);    if ((da == 0) || (db == 0)) {	Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);	return;    }    side = (side >= 0) ? (1) : (-1);	/* normalize variable */    dalpha *= PI / 180;		/* convert dalpha from degrees to radians */    angular_tol = angular_tolerance(tol, da, db);    i = np - 2;    norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);    elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);    angle1 = atan2(ty, tx);    nx = x[i] + vx;    ny = y[i] + vy;    mx = x[i + 1] + vx;    my = y[i + 1] + vy;    if (!round)	line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);    for (i = 0; i <= np - 2; i++) {	G_debug(4, "point %d, segment %d-%d", i, i, i + 1);	/* save the old values */	if (!round) {	    a0 = a1;	    b0 = b1;	    c0 = c1;	}	wx = vx;	wy = vy;	angle0 = angle1;	norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);	if ((tx == 0) && (ty == 0))	    continue;	elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);	angle1 = atan2(ty, tx);	nx = x[i] + vx;	ny = y[i] + vy;	mx = x[i + 1] + vx;	my = y[i + 1] + vy;	if (!round)	    line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);	delta_phi = angle1 - angle0;	if (delta_phi > PI)	    delta_phi -= 2 * PI;	else if (delta_phi <= -PI)	    delta_phi += 2 * PI;	/* now delta_phi is in [-pi;pi] */	turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15);	inner_corner = (side * delta_phi <= 0) && (!turns360);	/* if <line turns 360> and (<caps> and <not round>) */	if (turns360 && caps && (!round)) {	    norm_vector(0, 0, vx, vy, &tx, &ty);	    elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx, &ty);	    Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0);	    G_debug(4, " append point (c) x=%.16f y=%.16f", x[i] + wx + tx,		    y[i] + wy + ty);	    Vect_append_point(nPoints, nx + tx, ny + ty, 0);	/* nx == x[i] + vx, ny == y[i] + vy */	    G_debug(4, " append point (c) x=%.16f y=%.16f", nx + tx, ny + ty);	}	if ((!turns360) && (!round) && (!inner_corner)) {	    res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);	    if (res == 1) {		Vect_append_point(nPoints, rx, ry, 0);		G_debug(4, " append point (o) x=%.16f y=%.16f", rx, ry);//.........这里部分代码省略.........

void draw_histogram(const char *map_name, int x0, int y0, int width,		    int height, int color, int flip, int horiz,		    int map_type, int is_fp, struct FPRange render_range){    int i, nsteps, ystep;    long cell_count = 0;    double max_width, width_mult, dx;    double dy, y0_adjust;	/* only needed for CELL maps */    struct stat_list dist_stats;    struct stat_node *ptr;    struct Range range;    struct FPRange fprange;    CELL c_map_min, c_map_max;    DCELL d_map_min, d_map_max;    double map_min, map_max, map_range, user_range;    double crop_min_perc = 0.0, crop_max_perc = 1.0, pad_min_perc = 0.0;    if (horiz) {	max_width = height * 1.75;	nsteps = width - 3;    }    else {	max_width = width * 1.75;	nsteps = height - 3;    }    if (render_range.first_time) {	/* user specified range, can be either larger	    or smaller than actual map's range */	if (is_fp) {	    Rast_read_fp_range(map_name, "", &fprange);	    Rast_get_fp_range_min_max(&fprange, &d_map_min, &d_map_max);	    map_min = (double)d_map_min;	    map_max = (double)d_map_max;	}	else {	    Rast_read_range(map_name, "", &range);	    Rast_get_range_min_max(&range, &c_map_min, &c_map_max);	    map_min = (double)c_map_min;	    map_max = (double)c_map_max;	}	map_range = map_max - map_min;	user_range = render_range.max - render_range.min;	if (horiz)	    nsteps = (int)(0.5 + (map_range * (width - 3) / user_range));	else	    nsteps = (int)(0.5 + (map_range * (height - 3) / user_range));	G_debug(1, "number of steps for r.stats = %d, height-3=%d  width-3=%d",		nsteps, height - 3, width - 3);	/* need to know the % of the MAP range where user range starts and stops.	 *   note that MAP range can be fully inside user range, in which case	 *   keep 0-100% aka 0,nsteps, i.e. the step number in the nsteps range */	if (render_range.min > map_min) {	   crop_min_perc = (render_range.min - map_min) / map_range;	   G_debug(3, "min: %.02f vs. %.02f (%.02f) ... %.02f%%",	   	   render_range.min, map_min, map_range, 100 * crop_min_perc);	}	if (render_range.max > map_max) {	    crop_max_perc = 1.0 - ((render_range.max - map_max) / user_range);	    G_debug(3, "max: %.02f vs. %.02f (%.02f) ... %.02f%%",		    map_max, render_range.max, map_range, 100 * crop_max_perc);	}	if (render_range.min < map_min) {	   pad_min_perc = (map_min - render_range.min) / user_range;	   G_debug(3, "Min: %.02f vs. %.02f (%.02f) ... %.02f%%",	   	   map_min, render_range.min, user_range, 100 * pad_min_perc);	}#ifdef amplify_gain	/* proportion of nsteps to width, use as mult factor to boost the 1.75x	    when spread out over more nsteps than we are displaying */	G_debug(0, "max_width was: %.2f  (nsteps=%d)", max_width, nsteps);	if (nsteps > ((horiz ? width : height) - 3.0))	    max_width *= nsteps / ((horiz ? width : height) - 3.0);	G_debug(0, "max_width now: %.2f", max_width);#endif    }    /* TODO */    if (!is_fp && render_range.first_time) {	G_warning(_("Histogram constrained by range not yet implemented for "		  "categorical rasters"));	return;    }    /* get the distribution statistics */    get_stats(map_name, &dist_stats, nsteps, map_type);//.........这里部分代码省略.........

int dig_Wr_P_line(struct Plus_head *Plus, int n, struct gvfile * fp){    int n_edges = 0;    char ch;    struct P_line *ptr;    G_debug(4, "dig_Wr_P_line() line = %d", n);    ptr = Plus->Line[n];    /* if NULL i.e. dead write just 0 instead of type */    if (ptr == NULL) {	G_debug(4, "    line is dead -> write 0 only");	ch = 0;	if (0 >= dig__fwrite_port_C(&ch, 1, fp))	    return (-1);	return 0;    }    /* type */    ch = (char)dig_type_to_store(ptr->type);    G_debug(5, "    line type  %d -> %d", ptr->type, ch);    if (0 >= dig__fwrite_port_C(&ch, 1, fp))	return (-1);    /* offset */    if (0 >= dig__fwrite_port_O(&(ptr->offset), 1, fp, Plus->off_t_size))	return (-1);	    if (!ptr->topo)	return (0);	    /* nothing else for points */    /* centroids */    if (ptr->type & GV_CENTROID) {	struct P_topo_c *topo = (struct P_topo_c *)ptr->topo;		if (0 >= dig__fwrite_port_P(&(topo->area), 1, fp))	    return (-1);    }    /* lines */    else if (ptr->type & GV_LINE) {	struct P_topo_l *topo = (struct P_topo_l *)ptr->topo;	if (0 >= dig__fwrite_port_P(&(topo->N1), 1, fp))	    return (-1);	if (0 >= dig__fwrite_port_P(&(topo->N2), 1, fp))	    return (-1);    }    /* boundaries */    else if (ptr->type & GV_BOUNDARY) {	struct P_topo_b *topo = (struct P_topo_b *)ptr->topo;	if (0 >= dig__fwrite_port_P(&(topo->N1), 1, fp))	    return (-1);	if (0 >= dig__fwrite_port_P(&(topo->N2), 1, fp))	    return (-1);	if (0 >= dig__fwrite_port_P(&(topo->left), 1, fp))	    return (-1);	if (0 >= dig__fwrite_port_P(&(topo->right), 1, fp))	    return (-1);    }    /* faces */    else if ((ptr->type & GV_FACE) && Plus->with_z) {	/* reserved for face */	struct P_topo_f *topo = (struct P_topo_f *)ptr->topo;	if (0 >= dig__fwrite_port_I(&n_edges, 1, fp))	    return (-1);	/* here will be list of edges */	/* left / right volume / hole */	if (0 >= dig__fwrite_port_P(&(topo->left), 1, fp))	    return (-1);	if (0 >= dig__fwrite_port_P(&(topo->right), 1, fp))	    return (-1);    }    /* kernels */    else if ((ptr->type & GV_KERNEL) && Plus->with_z) {	/* reserved for kernel (volume number) */	struct P_topo_k *topo = (struct P_topo_k *)ptr->topo;	/* volume */	if (0 >= dig__fwrite_port_P(&(topo->volume), 1, fp))	    return (-1);    }    return (0);}

int close_array_seg(void){    struct Colors colors;    int incr, max, red, green, blue, rd, gr, bl, flag;    int c, r, map_fd;    CELL *cellrow, value;    CELL *theseg;    RAMSEG thesegseg;    cellrow = Rast_allocate_c_buf();    if (seg_flag || bas_flag || haf_flag) {	if (seg_flag) {	    theseg = bas;	    thesegseg = bas_seg;	}	else if (bas_flag) {	    theseg = bas;	    thesegseg = bas_seg;	}	else {	    theseg = haf;	    thesegseg = haf_seg;	}	max = n_basins;	G_debug(1, "%d basins created", max);	Rast_init_colors(&colors);	if (max > 0)	    Rast_make_random_colors(&colors, 1, max);	else {	    G_warning(_("No basins were created. Verify threshold and region settings."));	    Rast_make_random_colors(&colors, 1, 2);	}	if (max < 1000 && max > 0) {	    Rast_set_c_color((CELL) 0, 0, 0, 0, &colors);	    r = 1;	    incr = 0;	    while (incr >= 0) {		G_percent(r, max, 2);		for (gr = 130 + incr; gr <= 255; gr += 20) {		    for (rd = 90 + incr; rd <= 255; rd += 30) {			for (bl = 90 + incr; bl <= 255; bl += 40) {			    flag = 1;			    while (flag) {				Rast_get_c_color(&r, &red, &green, &blue, &colors);				/* if existing rule is too dark then append a new				   rule to override it */				if ((blue * .11 + red * .30 + green * .59) <				    100) {				    Rast_set_c_color(r, rd, gr, bl, &colors);				    flag = 0;				}				if (++r > max) {				    gr = rd = bl = 300;				    flag = 0;				    incr = -1;				}			    }			}		    }		}		if (incr >= 0) {		    incr += 15;		    if (incr > 120)			incr = 7;		}	    }	    G_percent(r - 1, max, 3);	/* finish it */	}	else	    G_debug(1,		    "Too many subbasins to reasonably check for color brightness");	/* using the existing stack of while/for/for/for/while loops ... */    }    /* stream segments map */    if (seg_flag) {	map_fd = Rast_open_c_new(seg_name);	for (r = 0; r < nrows; r++) {	    Rast_set_c_null_value(cellrow, ncols);	/* reset row to all NULL */	    for (c = 0; c < ncols; c++) {		value = FLAG_GET(swale, r, c);		if (value)		    cellrow[c] = bas[SEG_INDEX(bas_seg, r, c)];	    }	    Rast_put_row(map_fd, cellrow, CELL_TYPE);	}	Rast_close(map_fd);	Rast_write_colors(seg_name, this_mapset, &colors);    }    /* basins map */    if (bas_flag) {	map_fd = Rast_open_c_new(bas_name);	for (r = 0; r < nrows; r++) {	    for (c = 0; c < ncols; c++) {		cellrow[c] = bas[SEG_INDEX(bas_seg, r, c)];		if (cellrow[c] == 0)		    Rast_set_c_null_value(cellrow + c, 1);	    }//.........这里部分代码省略.........

/* clean_parallel - clean parallel line created by parallel_line: ** - looking for loops and if loop doesn't contain any other loop **   and centroid of loop is in buffer removes this loop (repeated) ** - optionally removes all end points in buffer *    parameters: *      Points - parallel line *      origPoints - original line *      d - offset *      rm_end - remove end points in buffer ** note1: on some lines (multiply selfcrossing; lines with end points **        in buffer of line other; some shapes of ends ) may create nosense ** note2: this function is stupid and slow, somebody more clever **        than I am should write paralle_line + clean_parallel **        better;    RB March 2000 */static void clean_parallel(struct line_pnts *Points,			   struct line_pnts *origPoints, double d, int rm_end){    int i, j, np, npn, sa, sb;    int sa_max = 0;    int first = 0, current, last, lcount;    double *x, *y, px, py, ix, iy;    static struct line_pnts *sPoints = NULL;    G_debug(4, "clean_parallel(): npoints = %d, d = %f, rm_end = %d",	    Points->n_points, d, rm_end);    x = Points->x;    y = Points->y;    np = Points->n_points;    if (sPoints == NULL)	sPoints = Vect_new_line_struct();    Vect_reset_line(sPoints);    npn = 1;    /* remove loops */    while (first < np - 2) {	/* find first loop which doesn't contain any other loop */	current = first;	last = Points->n_points - 2;	lcount = 0;	while (find_cross	       (Points, current, last - 1, current + 1, last, &sa,		&sb) != 0) {	    if (lcount == 0) {		first = sa;	    }			/* move first forward */	    current = sa + 1;	    last = sb;	    lcount++;	    G_debug(5, "  current = %d, last = %d, lcount = %d", current,		    last, lcount);	}	if (lcount == 0) {	    break;	}			/* loop not found */	/* ensure sa is monotonically increasing, so npn doesn't reset low */	if (sa > sa_max)	    sa_max = sa;	if (sa < sa_max)	    break;	/* remove loop if in buffer */	if ((sb - sa) == 1) {	/* neighbouring lines overlap */	    j = sb + 1;	    npn = sa + 1;	}	else {	    Vect_reset_line(sPoints);	    dig_find_intersection(x[sa], y[sa], x[sa + 1], y[sa + 1], x[sb],				  y[sb], x[sb + 1], y[sb + 1], &ix, &iy);	    Vect_append_point(sPoints, ix, iy, 0);	    for (i = sa + 1; i < sb + 1; i++) {	/* create loop polygon */		Vect_append_point(sPoints, x[i], y[i], 0);	    }	    Vect_find_poly_centroid(sPoints, &px, &py);	    if (point_in_buf(origPoints, px, py, d)) {	/* is loop in buffer ? */		npn = sa + 1;		x[npn] = ix;		y[npn] = iy;		j = sb + 1;		npn++;		if (lcount == 0) {		    first = sb;		}	    }	    else {		/* loop is not in buffer */		first = sb;		continue;	    }	}	for (i = j; i < Points->n_points; i++) {	/* move points down */	    x[npn] = x[i];	    y[npn] = y[i];//.........这里部分代码省略.........

int thin_streams(void){    int i, j, r, c, done;    CELL stream_id;    int next_node;    struct sstack    {	int stream_id;	int next_trib;    } *nodestack;    int top = 0, stack_step = 1000;    int n_trib_total;    int n_thinned = 0;    G_message(_("Thinning stream segments..."));    nodestack = (struct sstack *)G_malloc(stack_step * sizeof(struct sstack));    for (i = 0; i < n_outlets; i++) {	G_percent(i, n_outlets, 2);	r = outlets[i].r;	c = outlets[i].c;	cseg_get(&stream, &stream_id, r, c);	if (stream_id == 0)	    continue;	/* add root node to stack */	G_debug(2, "add root node");	top = 0;	nodestack[top].stream_id = stream_id;	nodestack[top].next_trib = 0;	/* depth first post order traversal */	G_debug(2, "traverse");	while (top >= 0) {	    done = 1;	    stream_id = nodestack[top].stream_id;	    G_debug(3, "stream_id %d, top %d", stream_id, top);	    if (nodestack[top].next_trib < stream_node[stream_id].n_trib) {		/* add to stack */		G_debug(3, "get next node");		next_node =		    stream_node[stream_id].trib[nodestack[top].next_trib];		G_debug(3, "add to stack: next %d, trib %d, n trib %d",			next_node, nodestack[top].next_trib,			stream_node[stream_id].n_trib);		nodestack[top].next_trib++;		top++;		if (top >= stack_step) {		    /* need more space */		    stack_step += 1000;		    nodestack =			(struct sstack *)G_realloc(nodestack,						   stack_step *						   sizeof(struct sstack));		}		nodestack[top].next_trib = 0;		nodestack[top].stream_id = next_node;		done = 0;		G_debug(3, "go further down");	    }	    if (done) {		/* thin stream segment */		G_debug(3, "thin stream segment %d", stream_id);		if (thin_seg(stream_id) == 0)		    G_debug(3, "segment %d not thinned", stream_id);		else {		    G_debug(3, "segment %d thinned", stream_id);		    n_thinned++;		}		top--;		/* count tributaries */		if (top >= 0) {		    n_trib_total = 0;		    stream_id = nodestack[top].stream_id;		    for (j = 0; j < stream_node[stream_id].n_trib; j++) {			/* intermediate */			if (stream_node[stream_node[stream_id].trib[j]].			    n_trib > 0)			    n_trib_total +=				stream_node[stream_node[stream_id].trib[j]].				n_trib_total;			/* start */			else			    n_trib_total++;		    }		    stream_node[stream_id].n_trib_total = n_trib_total;		}	    }	}    }    G_percent(n_outlets, n_outlets, 1);	/* finish it */    G_free(nodestack);    //.........这里部分代码省略.........

//.........这里部分代码省略.........		temp_clump++;	/* skip col */		if (*temp_clump == OLD)		    *temp_clump = NEW;	    }	    /* modify the OLD index */	    index[OLD] = NEW;	}	/* write initial clump IDs */	/* this works also with writing out cur_clump, but only 	 * prev_clump is complete and will not change any more */	if (row > 0) {	    if (write(cfd, prev_clump + 1, csize) != csize)		G_fatal_error(_("Unable to write to temp file"));	}	/* switch the buffers so that the current buffer becomes the previous */	temp_cell = cur_in;	cur_in = prev_in;	prev_in = temp_cell;	temp_clump = cur_clump;	cur_clump = prev_clump;	prev_clump = temp_clump;    }    /* write last row with initial clump IDs */    if (write(cfd, prev_clump + 1, csize) != csize)	G_fatal_error(_("Unable to write to temp file"));    G_percent(1, 1, 1);    /* generate a renumbering scheme */    G_message(_("Generating renumbering scheme..."));    G_debug(1, "%d initial labels", label);    /* allocate final clump ID */    renumber = (CELL *) G_malloc((label + 1) * sizeof(CELL));    renumber[0] = 0;    cat = 1;    G_percent(0, label, 1);    for (n = 1; n <= label; n++) {	G_percent(n, label, 1);	OLD = n;	NEW = index[n];	if (OLD != NEW) {	    renumber[n] = 0;	    /* find valid clump ID */	    while (OLD != NEW) {		OLD = NEW;		NEW = index[OLD];	    }	    index[n] = NEW;	}	else	    /* set final clump id */	    renumber[n] = cat++;    }        /* rewind temp file */    lseek(cfd, 0, SEEK_SET);    if (print) {	fprintf(stdout, "clumps=%d/n", cat - 1);    }    else {	/****************************************************	 *                      PASS 2                      *

int query(struct Map_info *Map){    int i, j, idx, cat_no, nlines, type;    register int line_num;    struct line_pnts *Points;    struct line_cats *Cats;    struct field_info *Fi;    dbString stmt, value_string;    dbDriver *driver;    /* Initialize the Point struct */    Points = Vect_new_line_struct();    Cats = Vect_new_cats_struct();    G_message(_("Reading features..."));    /* Cycle through all lines and make a list of categories of 'qfield' for each category given by 'field' */    nlines = Vect_get_num_lines(Map);    for (line_num = 1; line_num <= nlines; line_num++) {	type = Vect_read_line(Map, Points, Cats, line_num);	if (!(type & options.type))	    continue;	for (i = 0; i < Cats->n_cats; i++) {	    if (Cats->field[i] == options.field) {		cat_no = Cats->cat[i];		idx = find_cat(cat_no, 1);		for (j = 0; j < Cats->n_cats; j++) {		    if (Cats->field[j] == options.qfield) {	/* Add to list */			if (Values[idx].nqcats == Values[idx].aqcats) {	/* Alloc space */			    Values[idx].aqcats += 2;			    Values[idx].qcat =				(int *)G_realloc(Values[idx].qcat,						 Values[idx].aqcats *						 sizeof(int));			}			Values[idx].qcat[Values[idx].nqcats] = Cats->cat[j];			Values[idx].nqcats++;		    }		}	    }	}	/* If there is no field cat add cat -1, values for cat -1 are reported at the end  */	Vect_cat_get(Cats, options.field, &cat_no);	if (cat_no == -1) {	    idx = find_cat(cat_no, 1);	    for (j = 0; j < Cats->n_cats; j++) {		if (Cats->field[j] == options.qfield) {	/* Add to list */		    if (Values[idx].nqcats == Values[idx].aqcats) {	/* Alloc space */			Values[idx].aqcats += 2;			Values[idx].qcat =			    (int *)G_realloc(Values[idx].qcat,					     Values[idx].aqcats *					     sizeof(int));		    }		    Values[idx].qcat[Values[idx].nqcats] = Cats->cat[j];		    Values[idx].nqcats++;		}	    }	}	G_percent(line_num, nlines, 2);    }    db_init_string(&stmt);    db_init_string(&value_string);    if ((Fi = Vect_get_field(Map, options.qfield)) == NULL)	G_fatal_error(_("Database connection not defined for layer %d. Use v.db.connect first."),		      options.qfield);    /* Open driver */    driver = db_start_driver_open_database(Fi->driver, Fi->database);    if (driver == NULL)	G_fatal_error(_("Unable to open database <%s> by driver <%s>"),		      Fi->database, Fi->driver);    /* Query the database for each category */    G_message(_("Querying database... "));    for (i = 0; i < vstat.rcat; i++) {	int j, ctype, nrows, more;	char buf[2000];	dbCursor cursor;	dbTable *table;	dbColumn *column;	dbValue *value;	G_debug(3, "cat %d", Values[i].cat);	G_percent(i + 1, vstat.rcat, 1);	/* Skip if cat is zero and large number of query categories (many features without category).	 * It would cause problems on server side and take long time. Postgres limit is 10000 */	/* TODO: verify because no category is encoded as cat = -1, not cat = zero */	if (Values[i].cat == 0 && Values[i].nqcats > 1000) {//.........这里部分代码省略.........

/*!   /brief Load vector map to memory   The other alternative may be to load to a tmp file   /param grassname vector map name   /param[out] number of loaded features   /return pointer to geoline struct   /return NULL on failure */geoline *Gv_load_vect(const char *grassname, int *nlines){    struct Map_info map;    struct line_pnts *points;    struct line_cats *Cats = NULL;    geoline *top, *gln, *prev;    int np, i, n, nareas, nl = 0, area, type, is3d;    struct Cell_head wind;    float vect[2][3];    const char *mapset;    mapset = G_find_vector2(grassname, "");    if (!mapset) {	G_warning(_("Vector map <%s> not found"), grassname);	return NULL;    }    Vect_set_open_level(2);    if (Vect_open_old(&map, grassname, "") == -1) {	G_warning(_("Unable to open vector map <%s>"),		  G_fully_qualified_name(grassname, mapset));	return NULL;    }    top = gln = (geoline *) G_malloc(sizeof(geoline));	/* G_fatal_error */    if (!top) {	return NULL;    }    prev = top;#ifdef TRAK_MEM    Tot_mem += sizeof(geoline);#endif    points = Vect_new_line_struct();    Cats = Vect_new_cats_struct();    G_get_set_window(&wind);    Vect_set_constraint_region(&map, wind.north, wind.south, wind.east,			       wind.west, PORT_DOUBLE_MAX, -PORT_DOUBLE_MAX);    is3d = Vect_is_3d(&map);    /* Read areas */    n = Vect_get_num_areas(&map);    nareas = 0;    G_debug(3, "Reading vector areas (nareas = %d)", n);    for (area = 1; area <= n; area++) {	G_debug(3, " area %d", area);	Vect_get_area_points(&map, area, points);	if (points->n_points < 3)	    continue;	/* initialize style */	gln->highlighted = 0;	gln->type = OGSF_POLYGON;	gln->npts = np = points->n_points;	G_debug(3, "  np = %d", np);	if (is3d) {	    gln->dims = 3;	    gln->p3 = (Point3 *) G_calloc(np, sizeof(Point3));	/* G_fatal_error */	    if (!gln->p3) {		return (NULL);	    }#ifdef TRAK_MEM	    Tot_mem += (np * sizeof(Point3));#endif	}	else {	    gln->dims = 2;	    gln->p2 = (Point2 *) G_calloc(np, sizeof(Point2));	/* G_fatal_error */	    if (!gln->p2) {		return (NULL);	    }#ifdef TRAK_MEM	    Tot_mem += (np * sizeof(Point2));#endif	}	for (i = 0; i < np; i++) {	    if (is3d) {		gln->p3[i][X] = points->x[i];		gln->p3[i][Y] = points->y[i];		gln->p3[i][Z] = points->z[i];	    }	    else {//.........这里部分代码省略.........

voidasciiToG3d(FILE * fp, RASTER3D_Region * region, int convertNull, char *nullValue){    int x, y, z;    int col, row, depth;    double value;    char buff[256];    int tileX, tileY, tileZ;    Rast3d_get_tile_dimensions_map(map, &tileX, &tileY, &tileZ);    Rast3d_min_unlocked(map, RASTER3D_USE_CACHE_X);    Rast3d_autolock_on(map);    Rast3d_unlock_all(map);    G_message(_("Loading data ...  (%dx%dx%d)"), region->cols, region->rows,              region->depths);    G_debug(3,            "asciiToG3d: writing the 3D raster map, with rows %i cols %i depths %i",            region->rows, region->cols, region->depths);    for (z = 0; z < region->depths; z++) {        G_percent(z, region->depths, 1);        if ((z % tileZ) == 0)            Rast3d_unlock_all(map);        for (y = 0; y < region->rows; y++) /* go south to north */            for (x = 0; x < region->cols; x++) {                /* From west to east */                col = x;                /* The default is to read rows from north to south */                row = y;                /* From bottom to the top */                depth = z;                /* Read rows as from south to north */                if (rowOrder == ROW_ORDER_SOUTH_TO_NORTH)                    row = region->rows - y - 1;                /* Read XY layer from top to bottom */                if (depthOrder == DEPTH_ORDER_TOP_TO_BOTTOM)                    depth = region->depths - z - 1;                if (fscanf(fp, "%s", buff) != 1) {                    if (feof(fp))                        G_warning(_("End of file reached while still loading data."));                    G_debug(3,                            "missing data at col=%d row=%d depth=%d last_value=[%.4f]",                            x + 1, y + 1, z + 1, value);                    fatalError("asciiToG3d: read failed");                }                /* Check for null value */                if (convertNull && strncmp(buff, nullValue, strlen(nullValue)) == 0) {                    Rast3d_set_null_value(&value, 1, DCELL_TYPE);                } else {                    if (sscanf(buff, "%lf", &value) != 1) {                        G_warning(_("Invalid value detected"));                        G_debug(1, "invalid value at col=%d row=%d depth=%d last_value=[%s]",                                x + 1, y + 1, z + 1, buff);                        fatalError("asciiToG3d: read failed");                    }                }                /* Write the data */                Rast3d_put_double(map, col, row, depth, value);            }    }    if (fscanf(fp, "%lf", &value) == 1) {        G_warning(_("Data exists in input file after fully importing "                    "expected data.  [%.4f ...]"), value);    }    if (!Rast3d_flush_all_tiles(map))        fatalError("asciiTog3d: error flushing tiles");    Rast3d_autolock_off(map);    Rast3d_unlock_all(map);    G_percent(1, 1, 1);}

/*!   /brief Calculate normals   OPTIMIZED for constant dy & dx   The norm array is always the same size, but diff resolutions   force resampled data points to have their normals recalculated,   then only those norms are passed to n3f during drawing.   Norms are converted to a packed unsigned int for storage,   must be converted back at time of use.   /todo fix to correctly calculate norms when mapped to sphere!   Uses the previous and next cells (when available) for normal    calculations to produce smoother normals   /param gs surface (geosurf)   /return 1 on success   /return 0 on failure */int gs_calc_normals(geosurf * gs){    int row, col;    int xcnt, ycnt;    int xmod, ymod;    if (!gs->norm_needupdate || !gs->norms) {	return (0);    }    gs->norm_needupdate = 0;    gs_update_curmask(gs);    xmod = gs->x_mod;    ymod = gs->y_mod;    xcnt = VCOLS(gs);    ycnt = VROWS(gs);    init_vars(gs);    G_debug(5, "gs_calc_normals(): id=%d", gs->gsurf_id);    /* first row - just use single cell */    /* first col - use bottom & right neighbors */    calc_norm(gs, 0, 0, NBR);    for (col = 1; col < xcnt; col++) {	/* turn off top neighbor for first row */	calc_norm(gs, 0, col * xmod, ~NTOP);    }    /* use bottom & left neighbors for last col */    calc_norm(gs, 0, col * xmod, NBL);    /* now use four neighboring points for rows 1 - (n-1) */    for (row = 1; row < ycnt; row++) {	if (!(row % 100))	    G_debug(5, "gs_calc_normals(): row=%d", row);	/* turn off left neighbor for first col */	calc_norm(gs, row * ymod, 0, ~NLFT);	/* use all 4 neighbors until last col */	for (col = 1; col < xcnt; col++) {	    calc_norm(gs, row * ymod, col * xmod, NALL);	}	/* turn off right neighbor for last col */	calc_norm(gs, row * ymod, col * xmod, ~NRGT);    }    /* last row */    /* use top & right neighbors for first col */    calc_norm(gs, row * ymod, 0, NTR);    for (col = 1; col < xcnt; col++) {	/* turn off bottom neighbor for last row */	calc_norm(gs, row * ymod, col * xmod, ~NBOT);    }    /* use top & left neighbors for last column */    calc_norm(gs, row * ymod, col * xmod, NTL);    return (1);}

/* parallel_line - remove duplicate points from input line and *  creates new parallel line in 'd' offset distance; *  'tol' is tolerance between arc and polyline; *  this function doesn't care about created loops; * *  New line is written to existing nPoints structure. */static void parallel_line(struct line_pnts *Points, double d, double tol,			  struct line_pnts *nPoints){    int i, j, np, na, side;    double *x, *y, nx, ny, tx, ty, vx, vy, ux, uy, wx, wy;    double atol, atol2, a, av, aw;    G_debug(4, "parallel_line()");    Vect_reset_line(nPoints);    Vect_line_prune(Points);    np = Points->n_points;    x = Points->x;    y = Points->y;    if (np == 0)	return;    if (np == 1) {	Vect_append_point(nPoints, x[0], y[0], 0);	/* ? OK, should make circle for points ? */	return;    }    if (d == 0) {	Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);	return;    }    side = (int)(d / fabs(d));    atol = 2 * acos(1 - tol / fabs(d));    for (i = 0; i < np - 1; i++) {	vect(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);	vx = ty * d;	vy = -tx * d;	nx = x[i] + vx;	ny = y[i] + vy;	Vect_append_point(nPoints, nx, ny, 0);	nx = x[i + 1] + vx;	ny = y[i + 1] + vy;	Vect_append_point(nPoints, nx, ny, 0);	if (i < np - 2) {	/* use polyline instead of arc between line segments */	    vect(x[i + 1], y[i + 1], x[i + 2], y[i + 2], &ux, &uy);	    wx = uy * d;	    wy = -ux * d;	    av = atan2(vy, vx);	    aw = atan2(wy, wx);	    a = (aw - av) * side;	    if (a < 0)		a += 2 * PI;	    /* TODO: a <= PI can probably fail because of representation error */	    if (a <= PI && a > atol) {		na = (int)(a / atol);		atol2 = a / (na + 1) * side;		for (j = 0; j < na; j++) {		    av += atol2;		    nx = x[i + 1] + fabs(d) * cos(av);		    ny = y[i + 1] + fabs(d) * sin(av);		    Vect_append_point(nPoints, nx, ny, 0);		}	    }	}    }    Vect_line_prune(nPoints);}

/* Same as path() but get start/stop from the command line (for non-interactive use)   Hamish Bowman March 2007 */int coor_path(struct Map_info *Map, const struct color_rgb *hcolor,	      int be_bold, double start_x, double start_y,	      double end_x, double end_y){    int ret;    double nx, ny, fx, fy, tx, ty, msize, maxdist;    struct line_pnts *Points;    int start_node, end_node;    double fdist, tdist, cost;    Points = Vect_new_line_struct();    msize = 10 * (D_d_to_u_col(2.0) - D_d_to_u_col(1.0));	/* do it better */    G_debug(1, "msize = %f/n", msize);    /*       x1 = D_d_to_u_col ((double)(screen_x-WDTH));       y1 = D_d_to_u_row ((double)(screen_y-WDTH));       x2 = D_d_to_u_col ((double)(screen_x+WDTH));       y2 = D_d_to_u_row ((double)(screen_y+WDTH));       x1 = fabs ( x2 - x1 );       y1 = fabs ( y2 - y1 );       if ( x1 > y1 ) maxdist = x1;       else maxdist = y1;     *//**  maxdist = 10 pixels on the display (WDTH*2); ? **   ie related to zoom level ??  just use msize ?? **/    maxdist = msize;    G_debug(1, "Maximum distance in map units = %f/n", maxdist);    /* Vect_find_node(): find number of nearest node, 0 if not found */    start_node = Vect_find_node(Map, start_x, start_y, 0.0, maxdist, 0);    if (start_node > 0) {	Vect_get_node_coor(Map, start_node, &nx, &ny, NULL);	fprintf(stderr, _("Node %d: %f %f/n"), start_node, nx, ny);    }    if (start_node > 0) {	fx = nx;	fy = ny;    }    else {	fx = start_x;	fy = start_y;    }    D_RGB_color(hcolor->r, hcolor->g, hcolor->b);    D_plot_icon(fx, fy, G_ICON_BOX, 0.0, msize);    end_node = Vect_find_node(Map, end_x, end_y, 0.0, maxdist, 0);    if (end_node > 0) {	Vect_get_node_coor(Map, end_node, &nx, &ny, NULL);	fprintf(stderr, _("Node %d: %f %f/n"), end_node, nx, ny);    }    if (end_node > 0) {	tx = nx;	ty = ny;    }    else {	tx = end_x;	ty = end_y;    }    D_RGB_color(hcolor->r, hcolor->g, hcolor->b);    D_plot_icon(tx, ty, G_ICON_CROSS, 0.0, msize);    G_debug(2, "find path %f %f -> %f %f", fx, fy, tx, ty);    ret =	Vect_net_shortest_path_coor(Map, fx, fy, 0.0, tx, ty, 0.0,				    5 * maxdist, 5 * maxdist, &cost, Points,				    NULL, NULL, NULL, &fdist, &tdist);    if (ret == 0) {	fprintf(stdout, _("Destination unreachable/n"));    }    else {	fprintf(stdout, _("Costs on the network = %f/n"), cost);	fprintf(stdout, _("  Distance to the network = %f, "			  "distance from the network = %f/n/n"),		fdist, tdist);	display(Map, Points, hcolor, 1, 1, be_bold);    }    return 0;}

//.........这里部分代码省略.........    rows = Rast_window_rows();    cols = Rast_window_cols();    /* If we use x,y as seed... */    if (sdxy_opt->answer) {	G_get_window(&window);	east = window.east;	north = window.north;	G_scan_easting(sdxy_opt->answers[0], &east, G_projection());	G_scan_northing(sdxy_opt->answers[1], &north, G_projection());	start_col = (int)Rast_easting_to_col(east, &window);	start_row = (int)Rast_northing_to_row(north, &window);	if (start_row < 0 || start_row > rows ||	    start_col < 0 || start_col > cols)	    G_fatal_error(_("Seed point outside the current region"));    }    /* Open terrain map */    in_terran_fd = Rast_open_old(terrainmap, "");    /* Open seed map */    if (smap_opt->answer)	out_fd = Rast_open_old(seedmap, "");    /* Pointers to rows. Row = ptr to 'col' size array. */    in_terran = (FCELL **) G_malloc(rows * sizeof(FCELL *));    out_water = (FCELL **) G_malloc(rows * sizeof(FCELL *));    if (in_terran == NULL || out_water == NULL)	G_fatal_error(_("G_malloc: out of memory"));    G_debug(1, "Loading maps...");    /* foo_rows[row] == array with data (2d array). */    for (row = 0; row < rows; row++) {	in_terran[row] = (FCELL *) G_malloc(cols * sizeof(FCELL));	out_water[row] = (FCELL *) G_calloc(cols, sizeof(FCELL));	/* In newly created space load data from file. */	Rast_get_f_row(in_terran_fd, in_terran[row], row);	if (smap_opt->answer)	    Rast_get_f_row(out_fd, out_water[row], row);	G_percent(row + 1, rows, 5);    }    /* Set seed point */    if (sdxy_opt->answer)	/* Check is water level higher than seed point */	if (in_terran[start_row][start_col] >= water_level)	    G_fatal_error(_("Given water level at seed point is below earth surface. "			   "Increase water level or move seed point."));    out_water[start_row][start_col] = 1;    /* Close seed map for reading. */    if (smap_opt->answer)	Rast_close(out_fd);    /* Open output map for writing. */    if (lakemap)	out_fd = lake_fd;    else	out_fd = Rast_open_new(seedmap, 1);

int dig_Rd_P_line(struct Plus_head *Plus, int n, struct gvfile * fp){    int n_edges;    char tp;    struct P_line *ptr;    G_debug(4, "dig_Rd_P_line()");    if (0 >= dig__fread_port_C(&tp, 1, fp))	return (-1);    if (tp == 0) {		/* dead */	G_debug(4, "    line is dead");	Plus->Line[n] = NULL;	return 0;    }    ptr = dig_alloc_line();    /* type */    ptr->type = dig_type_from_store(tp);    G_debug(5, "    line type  %d -> %d", tp, ptr->type);    /* offset */    if (0 >= dig__fread_port_O(&(ptr->offset), 1, fp, Plus->off_t_size))	return (-1);    if (ptr->type == GV_POINT) {	ptr->topo = NULL;    }    else {	ptr->topo = dig_alloc_topo(ptr->type);    }    /* centroids */    if (ptr->type & GV_CENTROID) {	struct P_topo_c *topo = (struct P_topo_c *)ptr->topo;	if (0 >= dig__fread_port_P(&(topo->area), 1, fp))	    return -1;    }    /* lines */    else if (ptr->type & GV_LINE) {	struct P_topo_l *topo = (struct P_topo_l *)ptr->topo;	if (0 >= dig__fread_port_P(&(topo->N1), 1, fp))	    return -1;	if (0 >= dig__fread_port_P(&(topo->N2), 1, fp))	    return -1;    }    /* boundaries */    else if (ptr->type & GV_BOUNDARY) {	struct P_topo_b *topo = (struct P_topo_b *)ptr->topo;	if (0 >= dig__fread_port_P(&(topo->N1), 1, fp))	    return -1;	if (0 >= dig__fread_port_P(&(topo->N2), 1, fp))	    return -1;	if (0 >= dig__fread_port_P(&(topo->left), 1, fp))	    return -1;	if (0 >= dig__fread_port_P(&(topo->right), 1, fp))	    return -1;    }    /* faces */    else if ((ptr->type & GV_FACE) && Plus->with_z) {	/* reserved for face edges */	struct P_topo_f *topo = (struct P_topo_f *)ptr->topo;	if (0 >= dig__fread_port_I(&n_edges, 1, fp))	    return -1;	/* here will be list of edges */	/* left / right volume */	if (0 >= dig__fread_port_P(&(topo->left), 1, fp))	    return -1;	if (0 >= dig__fread_port_P(&(topo->left), 1, fp))	    return -1;    }    /* kernels */    else if ((ptr->type & GV_KERNEL) && Plus->with_z) {	/* reserved for kernel (volume number) */	struct P_topo_k *topo = (struct P_topo_k *)ptr->topo;	if (0 >= dig__fread_port_P(&(topo->volume), 1, fp))	    return -1;    }    Plus->Line[n] = ptr;    return (0);}

/* ************************************************************************* */void write_vtk_points(input_maps * in, FILE * fp, RASTER3D_Region region, int dp,                      int type, double scale){    int x, y, z, percentage = 0;    int rows, cols, depths;    void *rast_top = NULL;    void *rast_bottom = NULL;    void *ptr_top = NULL;    void *ptr_bottom = NULL;    double topval = 0, bottomval = 0;    double zcoor, ycoor, xcoor;    double zcoor1, ycoor1, xcoor1;    rows = region.rows;    cols = region.cols;    depths = region.depths;    rast_top = Rast_allocate_buf(in->topMapType);    rast_bottom = Rast_allocate_buf(in->bottomMapType);    G_debug(3, _("write_vtk_points: Writing point coordinates"));    for (z = 0; z < depths; z++) {        for (y = 0; y < rows; y++) {            G_percent(percentage, (rows * depths - 1), 10);            percentage++;            Rast_get_row(in->top, rast_top, y, in->topMapType);            Rast_get_row(in->bottom, rast_bottom, y, in->bottomMapType);            for (x = 0, ptr_top = rast_top, ptr_bottom = rast_bottom;                    x < cols;                    x++, ptr_top =                        G_incr_void_ptr(ptr_top, Rast_cell_size(in->topMapType)),                    ptr_bottom =                        G_incr_void_ptr(ptr_bottom,                                        Rast_cell_size(in->bottomMapType))) {                /*Get the values */                topval =                    get_raster_value_as_double(in->topMapType, ptr_top, 0.0);                bottomval =                    get_raster_value_as_double(in->bottomMapType, ptr_bottom,                                               0.0);                if (type == 1) { /*Structured Grid */                    /*Calculate the coordinates */                    xcoor =                        region.west + (region.ew_res / 2 +                                       region.ew_res * (x));                    /* Here the raster3d north->south coordinate system is used */                    ycoor =                        region.north - (region.ns_res / 2 +                                        region.ns_res * (y));                    zcoor =                        (bottomval +                         z * (topval - bottomval) / (depths - 1)) * scale;                    xcoor -= x_extent;                    ycoor -= y_extent;                    fprintf(fp, "%.*f ", dp, xcoor);                    fprintf(fp, "%.*f ", dp, ycoor);                    fprintf(fp, "%.*f/n", dp, zcoor);                } else { /*Unstructured Grid */                    /*Write for every cell the coordinates for a hexahedron -> 8 points */                    /*VTK Hexaeder */                    /* bottom                     * 3 --- 2                     * |     |                     * 0 --- 1                     * top                     * 7 --- 6                     * |     |                     * 4 --- 5                     */                    xcoor = region.west + (region.ew_res * (x)); /*0, 3, 4, 7 */                    /* Here the raster3d north->south coordinate system is used */                    ycoor = region.north - (region.ns_res * (y)); /*2, 3, 6, 7 */                    zcoor = (bottomval + z * (topval - bottomval) / (depths)) * scale; /*0, 1, 2, 3 */                    xcoor1 = region.west + (region.ew_res + region.ew_res * (x)); /*1, 2, 5, 6 */                    /* Here the raster3d north->south coordinate system is used */                    ycoor1 = region.north - (region.ns_res + region.ns_res * (y)); /*0, 1, 4, 5 */                    zcoor1 = (bottomval + z * (topval - bottomval) / (depths) + (topval - bottomval) / (depths)) * scale; /*4, 5, ,6 ,7 */                    xcoor -= x_extent;                    ycoor -= y_extent;                    xcoor1 -= x_extent;                    ycoor1 -= y_extent;                    /*0 */                    fprintf(fp, "%.*f ", dp, xcoor);//.........这里部分代码省略.........

/*!  /brief Read Plus_head from file  /param fp pointer to gvfile structure  /param[in,out] ptr pointer to Plus_head structure  /return -1 error  /return  0 OK */int dig_Rd_Plus_head(struct gvfile * fp, struct Plus_head *ptr){    unsigned char buf[5];    int byte_order;    dig_rewind(fp);    /* bytes 1 - 5 */    if (0 >= dig__fread_port_C((char *)buf, 5, fp))	return (-1);    ptr->version.topo.major = buf[0];    ptr->version.topo.minor = buf[1];    ptr->version.topo.back_major = buf[2];    ptr->version.topo.back_minor = buf[3];    byte_order = buf[4];    G_debug(2,	    "Topo header: file version %d.%d , supported from GRASS version %d.%d",	    ptr->version.topo.major, ptr->version.topo.minor, ptr->version.topo.back_major,	    ptr->version.topo.back_minor);    G_debug(2, "  byte order %d", byte_order);    /* check version numbers */    if (ptr->version.topo.major > GV_TOPO_VER_MAJOR ||	ptr->version.topo.minor > GV_TOPO_VER_MINOR) {	/* The file was created by GRASS library with higher version than this one */	if (ptr->version.topo.back_major > GV_TOPO_VER_MAJOR ||	    ptr->version.topo.back_minor > GV_TOPO_VER_MINOR) {	    /* This version of GRASS lib is lower than the oldest which can read this format */	    G_debug(1, "Topology format version %d.%d",		    ptr->version.topo.major, ptr->version.topo.minor);	    G_fatal_error		(_("This version of GRASS (%d.%d) is too old to read this topology format."		 " Try to rebuild topology or upgrade GRASS to at least version %d."),		 GRASS_VERSION_MAJOR, GRASS_VERSION_MINOR, GRASS_VERSION_MAJOR + 1);	    return (-1);	}	G_warning(_("Your GRASS version does not fully support topology format %d.%d of the vector."		    " Consider to rebuild topology or upgrade GRASS."),		  ptr->version.topo.major, ptr->version.topo.minor);    }    if (ptr->version.topo.major < GV_TOPO_VER_MAJOR ||	(ptr->version.topo.major == GV_TOPO_VER_MAJOR &&	 ptr->version.topo.minor < GV_TOPO_VER_MINOR)) {	/* The file was created by GRASS library with lower version than this one */	/* This version of GRASS lib can not read this old format */	G_warning(_("Old topology format version %d.%d is not supported by this release."		    " Try to rebuild topology."),		  ptr->version.topo.major, ptr->version.topo.minor);	return (-1);    }    /* init Port_info structure and set as default */    dig_init_portable(&(ptr->port), byte_order);    dig_set_cur_port(&(ptr->port));    /* bytes 6 - 9 : header size */    if (0 >= dig__fread_port_L(&(ptr->head_size), 1, fp))	return (-1);    G_debug(2, "  header size %ld", ptr->head_size);    /* determine required offset size from header size */    /* this is not safe in case new fields get added in later versions */    /* better: add a new field with off_t_size after byte_order? */    if (ptr->head_size >= 142 + 32) /* keep in sync with dig_Wr_Plus_head() */	ptr->off_t_size = 8;    else	ptr->off_t_size = 4;    if (sizeof(off_t) < ptr->off_t_size) {	G_warning(_("Vector exceeds supported file size limit"));	return (-1);    }    G_debug(2, "topo off_t size = %d", ptr->off_t_size);    /* byte 10 : dimension 2D or 3D */    if (0 >= dig__fread_port_C((char *)buf, 1, fp))	return (-1);    ptr->with_z = buf[0];    G_debug(2, "  with_z %d", ptr->with_z);    /* bytes 11 - 58 : bound box */    if (0 >= dig__fread_port_D(&(ptr->box.N), 1, fp))	return (-1);    if (0 >= dig__fread_port_D(&(ptr->box.S), 1, fp))	return (-1);//.........这里部分代码省略.........

/* ************************************************************************* */void write_vtk_unstructured_grid_cells(FILE * fp, RASTER3D_Region region){    int x, y, z, percentage;    int rows, cols, depths, count;    rows = region.rows;    cols = region.cols;    depths = region.depths;    G_debug(3, _("write_vtk_unstructured_grid_cells: Writing the cells"));    fprintf(fp, "CELLS %i %i/n", region.cols * region.rows * region.depths,            region.cols * region.rows * region.depths * 9);    count = 0;    percentage = 0;    /*The point - cell links */    for (z = 0; z < depths; z++) {        for (y = 0; y < rows; y++) {            G_percent(percentage, (rows * depths - 1), 10);            percentage++;            for (x = 0; x < cols; x++) {                /*Voxel */                fprintf(fp, "%i %i %i %i %i %i %i %i %i/n", 8,                        count * 8, count * 8 + 1, count * 8 + 3,                        count * 8 + 2, count * 8 + 4, count * 8 + 5,                        count * 8 + 7, count * 8 + 6);                /*Hexaeder                 * fprintf(fp, "%i %i %i %i %i %i %i %i %i/n", 8,                 * count * 8, count * 8 + 1, count * 8 + 2, count * 8 + 3,                 * count * 8 + 4, count * 8 + 5, count * 8 + 6,                 * count * 8 + 7);                 */                count++;            }        }    }    percentage = 0;    fprintf(fp, "CELL_TYPES %i/n", region.cols * region.rows * region.depths);    /*the cell types */    for (z = 0; z < depths; z++) {        for (y = 0; y < rows; y++) {            G_percent(percentage, (rows * depths - 1), 10);            percentage++;            for (x = 0; x < cols; x++) {                /*Voxel */                fprintf(fp, "11/n");                /*Hexaeder                 * fprintf(fp, "12/n");                 */            }        }    }    fprintf(fp, "CELL_DATA %i/n", region.cols * region.rows * region.depths); /*We have celldata  */    return;}

/* * This function generates parallel line (with loops, but not like the old ones). * It is not to be used directly for creating buffers. * + added elliptical buffers/par.lines support * * dalpha - direction of elliptical buffer major axis in degrees * da - distance along major axis * db: distance along minor (perp.) axis * side: side >= 0 - right side, side < 0 - left side * when (da == db) we have plain distances (old case) * round - 1 for round corners, 0 for sharp corners. (tol is used only if round == 1) */static void parallel_line(struct line_pnts *Points, double da, double db,			  double dalpha, int side, int round, int caps, int looped,			  double tol, struct line_pnts *nPoints){    int i, j, res, np;    double *x, *y;    double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;    double vx1, vy1, wx1, wy1;    double a0, b0, c0, a1, b1, c1;    double phi1, phi2, delta_phi;    double nsegments, angular_tol, angular_step;    int inner_corner, turns360;    G_debug(3, "parallel_line()");    if (looped && 0) {	/* start point != end point */	return;    }    Vect_reset_line(nPoints);    if (looped) {	Vect_append_point(Points, Points->x[1], Points->y[1], Points->z[1]);    }    np = Points->n_points;    x = Points->x;    y = Points->y;    if ((np == 0) || (np == 1))	return;    if ((da == 0) || (db == 0)) {	Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);	return;    }    side = (side >= 0) ? (1) : (-1);	/* normalize variable */    dalpha *= PI / 180;		/* convert dalpha from degrees to radians */    angular_tol = angular_tolerance(tol, da, db);    for (i = 0; i < np - 1; i++) {	/* save the old values */	a0 = a1;	b0 = b1;	c0 = c1;	wx = vx;	wy = vy;	norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);	if ((tx == 0) && (ty == 0))	    continue;	elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);	nx = x[i] + vx;	ny = y[i] + vy;	mx = x[i + 1] + vx;	my = y[i + 1] + vy;	line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);	if (i == 0) {	    if (!looped)		Vect_append_point(nPoints, nx, ny, 0);	    continue;	}	delta_phi = atan2(ty, tx) - atan2(y[i] - y[i - 1], x[i] - x[i - 1]);	if (delta_phi > PI)	    delta_phi -= 2 * PI;	else if (delta_phi <= -PI)	    delta_phi += 2 * PI;	/* now delta_phi is in [-pi;pi] */	turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15);	inner_corner = (side * delta_phi <= 0) && (!turns360);	if ((turns360) && (!(caps && round))) {	    if (caps) {		norm_vector(0, 0, vx, vy, &tx, &ty);		elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx,				 &ty);	    }	    else {		tx = 0;		ty = 0;//.........这里部分代码省略.........

/* ************************************************************************* */void write_vtk_data(FILE * fp, void *map, RASTER3D_Region region, char *varname,                    int dp){    double value;    double nullvalue;    int x, y, z, percentage;    int rows, cols, depths, typeIntern;    rows = region.rows;    cols = region.cols;    depths = region.depths;    /*the nullvalue */    if (!sscanf(param.null_val->answer, "%lf", &nullvalue)) {        G_warning("Null value is not valid, using 0 instead.");        nullvalue = 0;    }    G_debug(3,            _("write_vtk_data: Writing Celldata %s with rows %i cols %i depths %i to vtk-ascii file"),            varname, rows, cols, depths);    fprintf(fp, "SCALARS %s float 1/n", varname);    fprintf(fp, "LOOKUP_TABLE default/n");    typeIntern = Rast3d_tile_type_map(map);    percentage = 0;    for (z = 0; z < depths; z++) {        /* In case of structured grid data, the point/cell coordinates           are computed based on the default north->south raster3d coordinate system.           We need to compute south -> north ordering for image data.         */        if (!param.structgrid->answer) {            for (y = rows - 1; y >= 0; y--) {                G_percent(percentage, (rows * depths - 1), 10);                percentage++;                for (x = 0; x < cols; x++) {                    value =                        get_g3d_raster_value_as_double(map, x, y, z,                                                       typeIntern, nullvalue);                    fprintf(fp, "%.*f ", dp, value);                }                fprintf(fp, "/n");            }        } else {            for (y = 0; y < rows; y++) {                G_percent(percentage, (rows * depths - 1), 10);                percentage++;                for (x = 0; x < cols; x++) {                    value =                        get_g3d_raster_value_as_double(map, x, y, z,                                                       typeIntern, nullvalue);                    fprintf(fp, "%.*f ", dp, value);                }                fprintf(fp, "/n");            }        }    }}

/* * side: side >= 0 - extracts contour on right side of edge, side < 0 - extracts contour on left side of edge * if the extracted contour is the outer contour, it is returned in ccw order * else if it is inner contour, it is returned in cw order */static void extract_contour(struct planar_graph *pg, struct pg_edge *first,			    int side, int winding, int stop_at_line_end,			    struct line_pnts *nPoints){    int j;    int v;			/* current vertex number */    int v0;    int eside;			/* side of the current edge */    double eangle;		/* current edge angle with Ox (according to the current direction) */    struct pg_vertex *vert;	/* current vertex */    struct pg_vertex *vert0;	/* last vertex */    struct pg_edge *edge;	/* current edge; must be edge of vert */    /*    int cs; *//* on which side are we turning along the contour */    /* we will always turn right and dont need that one */    double opt_angle, tangle;    int opt_j, opt_side, opt_flag;    G_debug(3, "extract_contour(): v1=%d, v2=%d, side=%d, stop_at_line_end=%d",	    first->v1, first->v2, side, stop_at_line_end);    Vect_reset_line(nPoints);    edge = first;    if (side >= 0) {	eside = 1;	v0 = edge->v1;	v = edge->v2;    }    else {	eside = -1;	v0 = edge->v2;	v = edge->v1;    }    vert0 = &(pg->v[v0]);    vert = &(pg->v[v]);    eangle = atan2(vert->y - vert0->y, vert->x - vert0->x);    while (1) {	Vect_append_point(nPoints, vert0->x, vert0->y, 0);	G_debug(4, "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", v0,		v, eside, edge->v1, edge->v2);	G_debug(4, "ec: append point x=%.18f y=%.18f", vert0->x, vert0->y);	/* mark current edge as visited on the appropriate side */	if (eside == 1) {	    edge->visited_right = 1;	    edge->winding_right = winding;	}	else {	    edge->visited_left = 1;	    edge->winding_left = winding;	}	opt_flag = 1;	for (j = 0; j < vert->ecount; j++) {	    /* exclude current edge */	    if (vert->edges[j] != edge) {		tangle = vert->angles[j] - eangle;		if (tangle < -PI)		    tangle += 2 * PI;		else if (tangle > PI)		    tangle -= 2 * PI;		/* now tangle is in (-PI, PI) */		if (opt_flag || (tangle < opt_angle)) {		    opt_j = j;		    opt_side = (vert->edges[j]->v1 == v) ? (1) : (-1);		    opt_angle = tangle;		    opt_flag = 0;		}	    }	}	/* 	G_debug(4, "ec: opt: side=%d opt_flag=%d opt_angle=%.18f opt_j=%d opt_step=%d",	        side, opt_flag, opt_angle, opt_j, opt_step);	*/	/* if line end is reached (no other edges at curr vertex) */	if (opt_flag) {	    if (stop_at_line_end) {		G_debug(3, "    end has been reached, will stop here");		break;	    }	    else {		opt_j = 0;	/* the only edge of vert is vert->edges[0] */		opt_side = -eside;	/* go to the other side of the current edge */		G_debug(3, "    end has been reached, turning around");	    }	}	/* break condition */	if ((vert->edges[opt_j] == first) && (opt_side == side))	    break;//.........这里部分代码省略.........

/* ************************************************************************* */void write_vtk_rgb_data(void *map_r, void *map_g, void *map_b,                        FILE * fp, const char *varname,                        RASTER3D_Region region, int dp){    double value = 0;    int x, y, z, percentage, k;    int rows, cols, depths;    int typeIntern[3];    void *maprgb = NULL;    G_debug(3, "write_vtk_rgb_data: Writing RGB data");    rows = region.rows;    cols = region.cols;    depths = region.depths;    typeIntern[0] = Rast3d_tile_type_map(map_r);    typeIntern[1] = Rast3d_tile_type_map(map_g);    typeIntern[2] = Rast3d_tile_type_map(map_b);    percentage = 0;    /********************** WRITE RGB VOXEL DATA; CELL OR POINT ****************/    fprintf(fp, "COLOR_SCALARS %s 3/n", varname);    for (z = 0; z < depths; z++) {        for (y = 0; y < rows; y++) {            G_percent(percentage, (rows * depths - 1), 10);            percentage++;            for (x = 0; x < cols; x++) {                for (k = 0; k < 3; k++) {                    if (k == 0)                        maprgb = map_r;                    if (k == 1)                        maprgb = map_g;                    if (k == 2)                        maprgb = map_b;                    /* In case of structured grid data, the point/cell coordinates                       are computed based on the default north->south raster3d coordinate system.                       We need to compute south -> north ordering for image data.                     */                    if (!param.structgrid->answer)                        value =                            get_g3d_raster_value_as_double(maprgb, x, rows - y - 1, z,                                                           typeIntern[k],                                                           0.0);                    else                        value =                            get_g3d_raster_value_as_double(maprgb, x, y, z,                                                           typeIntern[k],                                                           0.0);                    /*Test of value range, the data should be 1 byte gray values */                    if (value > 255 || value < 0) {                        G_warning(_("Wrong 3D raster map values! Values should in between 0 and 255!"));                        fprintf(fp, "0 ");                    } else {                        fprintf(fp, "%.*f ", dp, (value / 255));                    }                }                fprintf(fp, "/n");            }        }    }    return;}

/* area_outer and area_isles[i] must be closed non self-intersecting lines   side: 0 - auto, 1 - right, -1 left */static void buffer_lines(struct line_pnts *area_outer, struct line_pnts **area_isles,			 int isles_count, int side, double da, double db,			 double dalpha, int round, int caps, double tol,			 struct line_pnts **oPoints, struct line_pnts ***iPoints,			 int *inner_count){    struct planar_graph *pg2;    struct line_pnts *sPoints, *cPoints;    struct line_pnts **arrPoints;    int i, count = 0;    int res, winding;    int auto_side;    int more = 8;    int allocated = 0;    double px, py;    G_debug(3, "buffer_lines()");    auto_side = (side == 0);    /* initializations */    sPoints = Vect_new_line_struct();    cPoints = Vect_new_line_struct();    arrPoints = NULL;    /* outer contour */    G_debug(3, "    processing outer contour");    *oPoints = Vect_new_line_struct();    if (auto_side)	side =	    get_polygon_orientation(area_outer->x, area_outer->y,				    area_outer->n_points -				    1) ? LEFT_SIDE : RIGHT_SIDE;    convolution_line(area_outer, da, db, dalpha, side, round, caps, tol,		     sPoints);    pg2 = pg_create(sPoints);    extract_outer_contour(pg2, 0, *oPoints);    res = extract_inner_contour(pg2, &winding, cPoints);    while (res != 0) {	if (winding == 0) {	    int check_poly = 1;	    double area_size;	    dig_find_area_poly(cPoints, &area_size);	    if (area_size == 0) {		G_warning(_("zero area size"));		check_poly = 0;	    }	    if (cPoints->x[0] != cPoints->x[cPoints->n_points - 1] ||		cPoints->y[0] != cPoints->y[cPoints->n_points - 1]) {		G_warning(_("Line was not closed"));		check_poly = 0;	    }	    if (check_poly && !Vect_point_in_poly(cPoints->x[0], cPoints->y[0], area_outer)) {		if (Vect_get_point_in_poly(cPoints, &px, &py) == 0) {		    if (!point_in_buf(area_outer, px, py, da, db, dalpha)) {			add_line_to_array(cPoints, &arrPoints, &count, &allocated,					  more);			cPoints = Vect_new_line_struct();		    }		}		else {		    G_warning(_("Vect_get_point_in_poly() failed"));		}	    }	}	res = extract_inner_contour(pg2, &winding, cPoints);    }    pg_destroy_struct(pg2);    /* inner contours */    G_debug(3, "    processing inner contours");    for (i = 0; i < isles_count; i++) {	if (auto_side)	    side =		get_polygon_orientation(area_isles[i]->x, area_isles[i]->y,					area_isles[i]->n_points -					1) ? RIGHT_SIDE : LEFT_SIDE;	convolution_line(area_isles[i], da, db, dalpha, side, round, caps,			 tol, sPoints);	pg2 = pg_create(sPoints);	extract_outer_contour(pg2, 0, cPoints);	res = extract_inner_contour(pg2, &winding, cPoints);	while (res != 0) {	    if (winding == -1) {		int check_poly = 1;		double area_size;		dig_find_area_poly(cPoints, &area_size);		if (area_size == 0) {		    G_warning(_("zero area size"));		    check_poly = 0;		}		if (cPoints->x[0] != cPoints->x[cPoints->n_points - 1] ||		    cPoints->y[0] != cPoints->y[cPoints->n_points - 1]) {//.........这里部分代码省略.........

/* ************************************************************************* */void write_vtk_vector_data(void *map_x, void *map_y, void *map_z,                           FILE * fp, const char *varname,                           RASTER3D_Region region, int dp){    double value = 0;    int x, y, z, percentage, k;    int rows, cols, depths;    int typeIntern[3];    void *mapvect = NULL;    G_debug(3, "write_vtk_vector_data: Writing vector data");    rows = region.rows;    cols = region.cols;    depths = region.depths;    typeIntern[0] = Rast3d_tile_type_map(map_x);    typeIntern[1] = Rast3d_tile_type_map(map_y);    typeIntern[2] = Rast3d_tile_type_map(map_z);    percentage = 0;    /********************** WRITE VECTOR DATA; CELL OR POINT ****************/    fprintf(fp, "VECTORS %s float/n", varname);    for (z = 0; z < depths; z++) { /*From the bottom to the top */        for (y = 0; y < rows; y++) {            G_percent(percentage, (rows * depths - 1), 10);            percentage++;            for (x = 0; x < cols; x++) {                for (k = 0; k < 3; k++) {                    if (k == 0)                        mapvect = map_x;                    if (k == 1)                        mapvect = map_y;                    if (k == 2)                        mapvect = map_z;                    /* In case of structured grid data, the point/cell coordinates                       are computed based on the default north->south raster3d coordinate system.                       We need to compute south -> north ordering for image data.                     */                    if (!param.structgrid->answer)                        value =                            get_g3d_raster_value_as_double(mapvect, x, rows - y - 1, z,                                                           typeIntern[k],                                                           0.0);                    else                        value =                            get_g3d_raster_value_as_double(mapvect, x, y, z,                                                           typeIntern[k],                                                           0.0);                    fprintf(fp, "%.*f ", dp, value);                }                fprintf(fp, "/n");            }        }    }    return;}

//.........这里部分代码省略.........    if (nsubregion_row < 0)	nsubregion_row = 0;    nsubregions = nsubregion_row * nsubregion_col;    elaboration_reg.south = original_reg.north;    last_row = FALSE;    while (last_row == FALSE) {	/* For each row */	P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,		      GENERAL_ROW);	if (elaboration_reg.north > original_reg.north) {	/* First row */	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,			  FIRST_ROW);	}	if (elaboration_reg.south <= original_reg.south) {	/* Last row */	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,			  LAST_ROW);	    last_row = TRUE;	}	nsply =	    ceil((elaboration_reg.north -		  elaboration_reg.south) / stepN) + 0.5;	/*	if (nsply > NSPLY_MAX)	    nsply = NSPLY_MAX;	*/	G_debug(1, "nsply = %d", nsply);	elaboration_reg.east = original_reg.west;	last_column = FALSE;	while (last_column == FALSE) {	/* For each column */	    subregion++;	    if (nsubregions > 1)		G_message(_("Processing subregion %d of %d..."), subregion, nsubregions);	    else /* v.outlier -e will report mean point distance: */		G_warning(_("No subregions found! Check values for 'ew_step' and 'ns_step' parameters"));	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,			  GENERAL_COLUMN);	    if (elaboration_reg.west < original_reg.west) {	/* First column */		P_set_regions(&elaboration_reg, &general_box, &overlap_box,			      dims, FIRST_COLUMN);	    }	    if (elaboration_reg.east >= original_reg.east) {	/* Last column */		P_set_regions(&elaboration_reg, &general_box, &overlap_box,			      dims, LAST_COLUMN);		last_column = TRUE;	    }	    nsplx =		ceil((elaboration_reg.east -		      elaboration_reg.west) / stepE) + 0.5;	    /*	    if (nsplx > NSPLX_MAX)

//.........这里部分代码省略.........	    }	    if (ivis_fd >= 0) {		if (vis[col] < 0)		    vis[col] = 0; /* negative visibility is invalid, print a WARNING ? */		/* round to nearest visibility bin */		/* rounding result: watch out for fp representation error */		vis[col] = ((int) (vis[col] + 0.5));	    }	    /* check if both maps are active and if whether any value has changed */	    if ((ialt_fd >= 0) && (ivis_fd >= 0) &&		((prev_vis != vis[col]) || (prev_alt != alt[col]))) {		prev_alt = alt[col];	/* update new values */		prev_vis = vis[col];		if (optimize) {		    int in_cache = ticache.search(alt[col], vis[col], &ti);		    if (!in_cache) {			pre_compute_hv(alt[col], vis[col]);	/* re-compute transformation inputs */			ti = compute();	/* ... */			ticache.add(ti, alt[col], vis[col]);		    }		}		else {		    pre_compute_hv(alt[col], vis[col]);	/* re-compute transformation inputs */		    ti = compute();	/* ... */		}	    }	    else {		/* only one of the maps is being used */		if ((ivis_fd >= 0) && (prev_vis != vis[col])) {		    prev_vis = vis[col];	/* keep track of previous visibility */		    if (optimize) {			int in_cache = ticache.search(0, vis[col], &ti);			if (!in_cache) {			    pre_compute_v(vis[col]);	/* re-compute transformation inputs */			    ti = compute();	/* ... */			    ticache.add(ti, 0, vis[col]);			}		    }		    else {			pre_compute_v(vis[col]);	/* re-compute transformation inputs */			ti = compute();	/* ... */		    }		}		if ((ialt_fd >= 0) && (prev_alt != alt[col])) {		    prev_alt = alt[col];	/* keep track of previous altitude */		    if (optimize) {			int in_cache = ticache.search(alt[col], 0, &ti);			if (!in_cache) {			    pre_compute_h(alt[col]);	/* re-compute transformation inputs */			    ti = compute();	/* ... */			    ticache.add(ti, alt[col], 0);			}		    }		    else {			pre_compute_h(alt[col]);	/* re-compute transformation inputs */			ti = compute();	/* ... */		    }		}	    }	    G_debug(3, "Computed r%d (%d), c%d (%d)", row, nrows, col, ncols);	    /* transform from iscale.[min,max] to [0,1] */	    buf[col] =		(buf[col] - iscale.min) / ((float)iscale.max -					   (float)iscale.min);	    buf[col] = transform(ti, imask, buf[col]);	    /* transform from [0,1] to oscale.[min,max] */	    buf[col] =		buf[col] * ((float)oscale.max - (float)oscale.min) +		oscale.min;	    if (oint && (buf[col] > (float)oscale.max))		G_warning(_("The output data will overflow. Reflectance > 100%%"));	}	/* write output */	if (oint)	    write_fp_to_cell(ofd, buf);	else	    Rast_put_row(ofd, buf, FCELL_TYPE);    }    G_percent(1, 1, 1);    /* free allocated memory */    G_free(buf);    if (ialt_fd >= 0)	G_free(alt);    if (ivis_fd >= 0)	G_free(vis);}