[GRASS-SVN] r44859 - grass/trunk/lib/vector/Vlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 4 05:11:44 EST 2011
Author: mmetz
Date: 2011-01-04 02:11:44 -0800 (Tue, 04 Jan 2011)
New Revision: 44859
Modified:
grass/trunk/lib/vector/Vlib/break_polygons.c
Log:
add file-based version of Vect_break_polygons()
Modified: grass/trunk/lib/vector/Vlib/break_polygons.c
===================================================================
--- grass/trunk/lib/vector/Vlib/break_polygons.c 2011-01-04 10:10:47 UTC (rev 44858)
+++ grass/trunk/lib/vector/Vlib/break_polygons.c 2011-01-04 10:11:44 UTC (rev 44859)
@@ -16,8 +16,10 @@
\author Update for GRASS 7 Markus Metz
*/
-#include <grass/config.h>
#include <stdlib.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
#include <math.h>
#include <grass/gis.h>
#include <grass/vector.h>
@@ -55,10 +57,29 @@
double x, y; /* coords */
double a1, a2; /* angles */
char cross; /* 0 - do not break, 1 - break */
+ char used; /* 0 - was not used to break line, 1 - was used to break line
+ * this is stored because points are automaticaly marked as cross, even if not used
+ * later to break lines */
} XPNT;
+typedef struct
+{
+ double a1, a2; /* angles */
+ char cross; /* 0 - do not break, 1 - break */
+ char used; /* 0 - was not used to break line, 1 - was used to break line
+ * this is stored because points are automaticaly marked as cross, even if not used
+ * later to break lines */
+} XPNT2;
+
+static int fpoint;
+
+/* Function called from RTreeSearch for point found */
+void srch(int id, int *arg)
+{
+ fpoint = id;
+}
+
/* function used by binary tree to compare items */
-
static int compare_xpnts(const void *Xpnta, const void *Xpntb)
{
XPNT *a, *b;
@@ -83,28 +104,301 @@
return 1;
}
-/*!
- \brief Break polygons in vector map.
+/* break polygons using a file-based search index */
+void
+Vect_break_polygons_file(struct Map_info *Map, int type, struct Map_info *Err)
+{
+ struct line_pnts *BPoints, *Points;
+ struct line_cats *Cats, *ErrCats;
+ int i, j, k, ret, ltype, broken, last, nlines;
+ int nbreaks;
+ struct RTree *RTree;
+ int apoints, npoints, nallpoints, nmarks;
+ XPNT2 XPnt;
+ struct Rect rect;
+ double dx, dy, a1 = 0, a2 = 0;
+ int closed, last_point;
+ char cross;
+ int fd, xpntfd;
+ char *filename;
+
+ G_debug(0, "File-based version of Vect_break_polygons()");
- Breaks lines specified by type in vector map. Points at
- intersections may be optionally written to error map. Input map
- must be opened on level 2 for update at least on GV_BUILD_BASE.
+ filename = G_tempfile();
+ fd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
+ RTree = RTreeNewIndex(fd, 0, 2);
+ remove(filename);
- Function is optimized for closed polygons rigs (e.g. imported from
- OGR) but with clean geometry - adjacent polygons mostly have
- identical boundary. Function creates database of ALL points in the
- map, and then is looking for those where polygons should be broken.
- Lines may be broken only at points existing in input map!
+ filename = G_tempfile();
+ xpntfd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
+ remove(filename);
- \param Map input map where polygons will be broken
- \param type type of line to be broken
- \param Err vector map where points at intersections will be written or NULL
+ BPoints = Vect_new_line_struct();
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+ ErrCats = Vect_new_cats_struct();
- \return
- */
+ nlines = Vect_get_num_lines(Map);
+ G_debug(3, "nlines = %d", nlines);
+ /* Go through all lines in vector, and add each point to structure of points,
+ * if such point already exists check angles of segments and if differ mark for break */
+
+ apoints = 0;
+ nmarks = 0;
+ npoints = 1; /* index starts from 1 ! */
+ nallpoints = 0;
+ XPnt.used = 0;
+
+ G_verbose_message(_("Break polygons Pass 1: select break points"));
+
+ for (i = 1; i <= nlines; i++) {
+ G_percent(i, nlines, 1);
+ G_debug(3, "i = %d", i);
+ if (!Vect_line_alive(Map, i))
+ continue;
+
+ ltype = Vect_read_line(Map, Points, Cats, i);
+ if (!(ltype & type))
+ continue;
+
+ /* This would be confused by duplicate coordinates (angle cannot be calculated) ->
+ * prune line first */
+ Vect_line_prune(Points);
+
+ /* If first and last point are identical it is close polygon, we dont need to register last point
+ * and we can calculate angle for first.
+ * If first and last point are not identical we have to mark for break both */
+ last_point = Points->n_points - 1;
+ if (Points->x[0] == Points->x[last_point] &&
+ Points->y[0] == Points->y[last_point])
+ closed = 1;
+ else
+ closed = 0;
+
+ for (j = 0; j < Points->n_points; j++) {
+ G_debug(3, "j = %d", j);
+ nallpoints++;
+
+ if (j == last_point && closed)
+ continue; /* do not register last of close polygon */
+
+ /* Box */
+ rect.boundary[0] = Points->x[j];
+ rect.boundary[3] = Points->x[j];
+ rect.boundary[1] = Points->y[j];
+ rect.boundary[4] = Points->y[j];
+ rect.boundary[2] = 0;
+ rect.boundary[5] = 0;
+
+ /* Already in DB? */
+ fpoint = -1;
+ RTreeSearch(RTree, &rect, (void *)srch, NULL);
+ G_debug(3, "fpoint = %d", fpoint);
+
+ if (Points->n_points <= 2 ||
+ (!closed && (j == 0 || j == last_point))) {
+ cross = 1; /* mark for cross in any case */
+ }
+ else { /* calculate angles */
+ cross = 0;
+ if (j == 0 && closed) { /* closed polygon */
+ dx = Points->x[last_point] - Points->x[0];
+ dy = Points->y[last_point] - Points->y[0];
+ a1 = atan2(dy, dx);
+ dx = Points->x[1] - Points->x[0];
+ dy = Points->y[1] - Points->y[0];
+ a2 = atan2(dy, dx);
+ }
+ else {
+ dx = Points->x[j - 1] - Points->x[j];
+ dy = Points->y[j - 1] - Points->y[j];
+ a1 = atan2(dy, dx);
+ dx = Points->x[j + 1] - Points->x[j];
+ dy = Points->y[j + 1] - Points->y[j];
+ a2 = atan2(dy, dx);
+ }
+ }
+
+ if (fpoint > 0) { /* Found */
+ /* read point */
+ lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
+ read(xpntfd, &XPnt, sizeof(XPNT2));
+ if (XPnt.cross == 1)
+ continue; /* already marked */
+
+ /* Check angles */
+ if (cross) {
+ XPnt.cross = 1;
+ nmarks++;
+ /* write point */
+ lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
+ write(xpntfd, &XPnt, sizeof(XPNT2));
+ }
+ else {
+ G_debug(3, "a1 = %f xa1 = %f a2 = %f xa2 = %f", a1,
+ XPnt.a1, a2, XPnt.a2);
+ if ((a1 == XPnt.a1 && a2 == XPnt.a2) ||
+ (a1 == XPnt.a2 && a2 == XPnt.a1)) { /* identical */
+
+ }
+ else {
+ XPnt.cross = 1;
+ nmarks++;
+ /* write point */
+ lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
+ write(xpntfd, &XPnt, sizeof(XPNT2));
+ }
+ }
+ }
+ else {
+ /* Add to tree and to structure */
+ RTreeInsertRect(&rect, npoints, RTree);
+ if (j == 0 || j == (Points->n_points - 1) ||
+ Points->n_points < 3) {
+ XPnt.a1 = 0;
+ XPnt.a2 = 0;
+ XPnt.cross = 1;
+ nmarks++;
+ }
+ else {
+ XPnt.a1 = a1;
+ XPnt.a2 = a2;
+ XPnt.cross = 0;
+ }
+ /* write point */
+ lseek(xpntfd, (off_t) (npoints - 1) * sizeof(XPNT2), SEEK_SET);
+ write(xpntfd, &XPnt, sizeof(XPNT2));
+
+ npoints++;
+ }
+ }
+ }
+
+ nbreaks = 0;
+
+ /* Second loop through lines (existing when loop is started, no need to process lines written again)
+ * and break at points marked for break */
+
+ G_verbose_message(_("Break polygons Pass 2: break at selected points"));
+
+ for (i = 1; i <= nlines; i++) {
+ int n_orig_points;
+
+ G_percent(i, nlines, 1);
+ G_debug(3, "i = %d", i);
+ if (!Vect_line_alive(Map, i))
+ continue;
+
+ ltype = Vect_read_line(Map, Points, Cats, i);
+ if (!(ltype & type))
+ continue;
+ if (!(ltype & GV_LINES))
+ continue; /* Nonsense to break points */
+
+ /* Duplicates would result in zero length lines -> prune line first */
+ n_orig_points = Points->n_points;
+ Vect_line_prune(Points);
+
+ broken = 0;
+ last = 0;
+ G_debug(3, "n_points = %d", Points->n_points);
+ for (j = 1; j < Points->n_points; j++) {
+ G_debug(3, "j = %d", j);
+ nallpoints++;
+
+ /* Box */
+ rect.boundary[0] = Points->x[j];
+ rect.boundary[3] = Points->x[j];
+ rect.boundary[1] = Points->y[j];
+ rect.boundary[4] = Points->y[j];
+ rect.boundary[2] = 0;
+ rect.boundary[5] = 0;
+
+ if (Points->n_points <= 1 ||
+ (j == (Points->n_points - 1) && !broken))
+ break;
+ /* One point only or
+ * last point and line is not broken, do nothing */
+
+ RTreeSearch(RTree, &rect, (void *)srch, 0);
+ G_debug(3, "fpoint = %d", fpoint);
+
+ /* read point */
+ lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
+ read(xpntfd, &XPnt, sizeof(XPNT2));
+
+ /* break or write last segment of broken line */
+ if ((j == (Points->n_points - 1) && broken) ||
+ XPnt.cross) {
+ Vect_reset_line(BPoints);
+ for (k = last; k <= j; k++) {
+ Vect_append_point(BPoints, Points->x[k], Points->y[k],
+ Points->z[k]);
+ }
+
+ /* Result may collapse to one point */
+ Vect_line_prune(BPoints);
+ if (BPoints->n_points > 1) {
+ ret = Vect_write_line(Map, ltype, BPoints, Cats);
+ G_debug(3,
+ "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d",
+ ret, j, Points->n_points, BPoints->n_points);
+ }
+
+ if (!broken)
+ Vect_delete_line(Map, i); /* not yet deleted */
+
+ /* Write points on breaks */
+ if (Err) {
+ if (XPnt.cross && !XPnt.used) {
+ Vect_reset_line(BPoints);
+ Vect_append_point(BPoints, Points->x[j], Points->y[j], 0);
+ Vect_write_line(Err, GV_POINT, BPoints, ErrCats);
+ }
+ if (!XPnt.used) {
+ XPnt.used = 1;
+ /* write point */
+ lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
+ write(xpntfd, &XPnt, sizeof(XPNT2));
+ }
+ }
+
+ last = j;
+ broken = 1;
+ nbreaks++;
+ }
+ }
+ if (!broken && n_orig_points > Points->n_points) { /* was pruned before -> rewrite */
+ if (Points->n_points > 1) {
+ Vect_rewrite_line(Map, i, ltype, Points, Cats);
+ G_debug(3, "Line %d pruned, npoints = %d", i,
+ Points->n_points);
+ }
+ else {
+ Vect_delete_line(Map, i);
+ G_debug(3, "Line %d was deleted", i);
+ }
+ }
+ else {
+ G_debug(3, "Line %d was not changed", i);
+ }
+ }
+
+ close(RTree->fd);
+ RTreeFreeIndex(RTree);
+ close(xpntfd);
+ Vect_destroy_line_struct(Points);
+ Vect_destroy_line_struct(BPoints);
+ Vect_destroy_cats_struct(Cats);
+ Vect_destroy_cats_struct(ErrCats);
+ G_verbose_message(_("Breaks: %d"), nbreaks);
+}
+
+
+/* break polygons using a memory-based search index */
void
-Vect_break_polygons(struct Map_info *Map, int type, struct Map_info *Err)
+Vect_break_polygons_mem(struct Map_info *Map, int type, struct Map_info *Err)
{
struct line_pnts *BPoints, *Points;
struct line_cats *Cats, *ErrCats;
@@ -116,6 +410,8 @@
double dx, dy, a1 = 0, a2 = 0;
int closed, last_point, cross;
+ G_debug(0, "Memory-based version of Vect_break_polygons()");
+
RBTree = rbtree_create(compare_xpnts, sizeof(XPNT));
BPoints = Vect_new_line_struct();
@@ -133,6 +429,7 @@
nmarks = 0;
npoints = 0;
nallpoints = 0;
+ XPnt_search.used = 0;
G_verbose_message(_("Break polygons Pass 1: select break points"));
@@ -242,7 +539,7 @@
nbreaks = 0;
nallpoints = 0;
- G_debug(2, "Break polygons: unique vertices: %d", RBTree->count);
+ G_debug(2, "Break polygons: unique vertices: %ld", (long int)RBTree->count);
/* uncomment to check if search tree is healthy */
/* if (rbtree_debug(RBTree, RBTree->root) == 0)
@@ -316,11 +613,12 @@
/* Write points on breaks */
if (Err) {
- if (j < (Points->n_points - 1)) {
+ if (XPnt_found->cross && !XPnt_found->used) {
Vect_reset_line(BPoints);
Vect_append_point(BPoints, Points->x[j], Points->y[j], 0);
Vect_write_line(Err, GV_POINT, BPoints, ErrCats);
}
+ XPnt_found->used = 1;
}
last = j;
@@ -348,6 +646,34 @@
Vect_destroy_line_struct(Points);
Vect_destroy_line_struct(BPoints);
Vect_destroy_cats_struct(Cats);
- Vect_destroy_cats_struct(ErrCats);
G_verbose_message(_("Breaks: %d"), nbreaks);
}
+
+/*!
+ \brief Break polygons in vector map.
+
+ Breaks lines specified by type in vector map. Points at
+ intersections may be optionally written to error map. Input map
+ must be opened on level 2 for update at least on GV_BUILD_BASE.
+
+ Function is optimized for closed polygons rigs (e.g. imported from
+ OGR) but with clean geometry - adjacent polygons mostly have
+ identical boundary. Function creates database of ALL points in the
+ map, and then is looking for those where polygons should be broken.
+ Lines may be broken only at points existing in input map!
+
+ \param Map input map where polygons will be broken
+ \param type type of line to be broken
+ \param Err vector map where points at intersections will be written or NULL
+
+ \return
+ */
+
+void
+Vect_break_polygons(struct Map_info *Map, int type, struct Map_info *Err)
+{
+ if (getenv("GRASS_VECTOR_LOWMEM"))
+ return Vect_break_polygons_file(Map, type, Err);
+ else
+ return Vect_break_polygons_mem(Map, type, Err);
+}
More information about the grass-commit
mailing list