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

void
asciiToG3d(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 /