[GRASS-dev] [INFO] Line simplification algorithm
Hamish
hamish_nospam at yahoo.com
Mon Nov 20 21:31:52 EST 2006
tlaronde polynum.com wrote:
> FWIW, I seem to remember a thread on your list about the line
> simplification algorithm, and the need of implementation of
> Douglas/Peuker kind of algorithm.
>
> ISTR too that when I used GRASS GPL (back in 2003), v.prune(1) was
> doing its job the worst way possible (making local simplification,
> with a sliding pair, i.e. cumulating errors by comparing a vertex
> against the previous one and discarding the previous if it was in the
> threshold, and then comparing with next, leading to the disappearance
> of the curves).
"v.prune" in GRASS 5 is essentially unchanged since moving into CVS in
1999:
http://freegis.org/cgi-bin/viewcvs.cgi/grass/src/mapdev/v.prune/v.prune.c.diff?r1=1.1&r2=1.3
(and the diff between GRASS 4.3 and 5.4 versions attached; not much
there either)
In GRASS 6 there is
vector/v.clean/prune.c by Radim Blazek
which calls
lib/vector/Vlib/line.c
Vect_line_prune() # removes duplicate verticies
Vect_line_prune_thresh() # this calls dig_prune()
which leads to:
lib/vector/diglib/prune.c by Michel Wurtz for GRASS 4.2.1
great header comments, Dave Gerdes is a "probable" in the copyright,
but.."by Michel Wurtz" ??
"This is a complete rewriting of the previous dig_prune subroutine."
there was a thread a couple of months ago about GRASS 6's v.clean pruning
problems and the meaning of thresh:
http://thread.gmane.org/gmane.comp.gis.grass.user/15726
more on GRASS 4.3's version:
grass43/src/mapdev/v.prune/CHANGES:
-----------------------------------------------------------------
So for the benefit of the grass community, there is a new
dig_prune function, commented (in bad english), that should
replace src421/src/mapdev/diglib/prune.c (and not only for
linux version, the code is no linux dependant)
Also edited: src421/src/mapdev/v.prune/v.prune.c, replacing
the comparaison "==" by "<" in line 114 (test for Threshold != 0)
A null value just eliminates duplicate points.
The code is shorter, but the source bigger (many comments).
Michel Wurtz mw at engees.u-strasbg.fr Feb 18, 1998
-------------------------------------------------------------------
and these files-
$ grass43/src$ find | grep prune
./general/g.help/help/Commands.def/vprune.def
./general/g.help/help/17.manual/Help.pages/v.prune
./mapdev/libes/old_dev/prune.c # /* @(#)prune.c 2.1 6/26/87 */
+/* added Feb 26, 1987 -dpg */
+dig_prune() and dlg_prune()
./mapdev/diglib/prune.c # /* @(#)prune.c 3.0 2/19/98 */
/* by Michel Wurtz for GRASS 4.2.1
./mapdev/diglib/prune.c.veryold # /* @(#)prune.c 2.1 6/26/87 */
(without Feb '87 updates)
./mapdev/v.prune
./mapdev/v.prune/CHANGES # cropped section above
./mapdev/v.prune/Gmakefile
./mapdev/v.prune/prune.c # by Dave Gerdes 6/89
./mapdev/v.prune/v.prune.c # see attached diff to grass5
/*
** - Written by Dave Gerdes 6/89
** US Army Construction Engineering Research Lab
** - New parser installed, 12/90 David Stigberg
** - corrected line 115 (and rewritten ../diglib/prune.c) 2/98 Michel Wurtz
** - added v.support check line 233 Markus Neteler 7/98
*/
./mapdev/OLDlib/prune.c # /* @(#)prune.c 2.1 6/26/87 */
# bugfix better than ./mapdev/diglib/prune.c.veryold
./mapdev/v.out.moss/prune_points.c
./xgrass/xclip/v.prune
./xgrass/xhelp/17.manual/Help.pages/v.prune
./xgrass/xhelp/Commands.def/vprune.def
./tcltkgrass/module/v.prune
which version have you used?
Perhaps some insights in there, please let us know!
> Looking at what I have in KerGIS (that is the original program with
> Dave Gerdes implementation) I'm puzzled... since this _is_ a
> Douglas/Peuker kind of algorithm with some slight modifications that
> make sense.
>
> So if you want it, just go back to the origin (adapting to your new
> vector engine).
I'm not sure how far we are from there actually?
> One question remains: why and who discard a working program to replace
> it with a disaster?
maybe something in the above,
> To be honest, working with the original sources, I know who I can
> trust (developer or programmer name)---and I know I have (had) only
> cosmetic work to fix the program---and who I can highly distrust---and
> who has made such a poor job that is eating my time not in proportion
> to his contribution (I will have spent 90% of the time fixing 10%
> of the sources and know I will have to redo later everything from
> scratch). Blunders may occur---been there, done that. But there
> are very suspicious blunders...
You don't have to name the guilty if you don't want to :), but who do
you consider to be really good?
Hamish
-------------- next part --------------
--- grass43/src/mapdev/v.prune/v.prune.c 1999-11-02 05:22:39.000000000 +1300
+++ grass-5.4.0/src/mapdev/v.prune/v.prune.c 2000-12-30 02:56:22.000000000 +1300
@@ -6,17 +6,16 @@
** - added v.support check line 233 Markus Neteler 7/98
*/
#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
#include "gis.h"
#include "Vect.h"
#include "dig_atts.h"
-/*#include "dig_head.h" */
#define MAIN
#define USAGE "v.prune [-i] input=name output=name thresh=value\n"
-long ftell ();
-double atof ();
double threshold;
int *pruned;
int inches;
@@ -29,8 +28,6 @@
static char *dig_name = NULL ;
static char *out_name = NULL ;
-static int load_args() ;
-
#ifdef OLDPARSE
struct Command_keys vars[] = {
{ "dig_in", 1 },
@@ -45,21 +42,29 @@
#endif /*OLDPARSE*/
double dig_unit_conversion ();
+int debugf(char *, ...);
+int export(char *, char *, char *);
+int doit(struct Map_info *, struct Map_info *);
+int get_line_center(double *, double *, struct line_pnts *);
+int get_pruned_area_points(struct Map_info *, int, struct line_pnts *);
+int get_pruned_isle_points(struct Map_info *, int, struct line_pnts *);
static int snapped = 0 ;
-main (argc, argv)
-int argc;
-char **argv;
+int main (int argc, char **argv)
{
- int ret;
char *mapset;
char errmsg[100];
+ struct GModule *module;
struct Option *old, *new, *thresh;
struct Flag *inch_flag;
G_gisinit(argv[0]);
+ module = G_define_module();
+ module->description =
+ "Prunes points from binary GRASS vector data files.";
+
inch_flag = G_define_flag();
inch_flag->key = 'i';
inch_flag->description = "set threshold value in inches";
@@ -125,7 +130,7 @@
/*DEBUG*/
/*
fprintf (stderr, "in = '%s' out = '%s'\n", dig_name, out_name);
- fprintf (stderr, "thresh = %lf inches = %d\n", threshold, inches);
+ fprintf (stderr, "thresh = %f inches = %d\n", threshold, inches);
exit(0);
*/
#ifdef OLDPARSE
@@ -150,7 +155,7 @@
#endif /*OLDPARSE*/
/* Show advertising */
- printf("\n\n Prune: threshold = %lf\n\n", threshold) ;
+ fprintf (stdout,"\n\n Prune: threshold = %f\n\n", threshold) ;
if ((mapset = G_find_vector2 (dig_name, "")) == NULL)
{
@@ -165,10 +170,7 @@
}
#ifdef OLDPARSE
-static
-load_args (position, str)
-int position;
-char *str;
+static int load_args (int position, char *str)
{
switch(position)
{
@@ -194,11 +196,14 @@
#endif /*OLDPARSE*/
#ifdef DEBUG
-debugf (format, a, b, c, d, e, f, g, h, i, j, k, l)
-char *format;
-int a, b, c, d, e, f, g, h, i, j, k, l;
+int debugf (char *format, a)
{
- fprintf (stderr, format, a, b, c, d, e, f, g, h, i, j, k, l);
+ va_list a;
+ va_start(format,a);
+ fprintf (stderr, format, a);
+ va_end(format,a);
+
+ return 0;
}
#endif
@@ -208,13 +213,12 @@
/*ORIGINAL struct head Head; */
struct dig_head d_head;
-export(dig_name, mapset, out_name)
-char *dig_name, *mapset, *out_name;
+int export (char *dig_name, char *mapset, char *out_name)
{
FILE *In;
- char buf[1024];
double x,y;
- int type, cat;
+ int cat;
+ char type;
long offset;
char errmsg[200];
struct Map_info InMap, OutMap;
@@ -285,18 +289,15 @@
return(0) ;
}
-doit (InMap, OutMap)
-struct Map_info *InMap, *OutMap;
+int doit (struct Map_info *InMap, struct Map_info *OutMap)
{
struct line_pnts *Points, *AreaPoints, **IslesPoints;
register int old_n_points, area, line, type;
- int att_ind, i, j;
- double l, delta, att_x, att_y, x1, y1, x2, y2;
- int point_found, binary;
- long offset;
+ int att_ind, i;
+ double att_x, att_y, x1;
+ int point_found;
int diff, max_n_isles=0;
int left;
- extern double sqrt();
/*Points.alloc_points = 0; */
@@ -305,7 +306,7 @@
AreaPoints = Vect_new_line_struct ();
IslesPoints = NULL;
- /*DEBUG*/ fprintf (stderr, "Resultant threshold = %lf\n", threshold);
+ /*DEBUG*/ fprintf (stderr, "Resultant threshold = %f\n", threshold);
left = diff = 0;
pruned = (int *) G_malloc(sizeof( int) * (InMap->n_lines + 1));
/* flags for each line indicating if the line was pruned or not */
@@ -365,7 +366,7 @@
att_ind = InMap->Area[area].att;
if(0>(AreaPoints->n_points=get_pruned_area_points (InMap, area, AreaPoints)))
{
- printf("Warning ! Can't read area %d");
+ fprintf (stdout,"Warning ! Can't read area %d",area);
continue;
}
if((x1=dig_point_in_poly(InMap->Att[att_ind].x, InMap->Att[att_ind].y, AreaPoints)) == 0.0)
@@ -385,7 +386,7 @@
{
IslesPoints[i]->alloc_points = 0;
IslesPoints[i]->n_points = get_pruned_isle_points(InMap, InMap->Area[area].isles[i], IslesPoints[i]);
- if(IslesPoints[i]->n_points==0) printf ("WARNING ");
+ if(IslesPoints[i]->n_points==0) fprintf (stdout,"WARNING ");
} /* isles loop */
if(InMap->Area[area].n_isles)
@@ -401,11 +402,11 @@
{
write_att(Out, FILE_AREA, att_x, att_y, InMap->Att[att_ind].cat);
}
- else printf("WARNING: the area is empty!\n");
+ else fprintf (stdout,"WARNING: the area is empty!\n");
} /* if not inside area */
else
{
- if(x1<0.0) printf("couldn't read the line!\n");
+ if(x1<0.0) fprintf (stdout,"couldn't read the line!\n");
write_att(Out, FILE_AREA,
InMap->Att[att_ind].x, InMap->Att[att_ind].y,
InMap->Att[att_ind].cat);
@@ -416,6 +417,8 @@
Vect_destroy_line_struct(AreaPoints);
for (i=0;i<max_n_isles;i++)
Vect_destroy_line_struct(IslesPoints[i]);
+
+ return 0;
}
/*
@@ -429,9 +432,9 @@
** return 0 on success ,-1 on error
*/
-get_line_center (x, y, Points)
- double *x, *y;
- struct line_pnts *Points;
+int get_line_center (
+ double *x,double *y,
+ struct line_pnts *Points)
{
register int i;
register int n_points;
@@ -440,7 +443,6 @@
double half_dist; /* half total line length */
double len; /* tmp length of current line seg */
double frac; /* overshoot / line length */
- double fabs ();
n_points = Points->n_points;
ux = Points->x;
@@ -507,11 +509,8 @@
static int first_time = 1; /* zero at startup */
static struct line_pnts TPoints;
-int
-get_pruned_area_points (Map, area, BPoints)
- struct Map_info *Map;
- int area;
- struct line_pnts *BPoints;
+int
+get_pruned_area_points (struct Map_info *Map, int area, struct line_pnts *BPoints)
{
register int i, line;
int start, end, to, from, inc;
@@ -569,11 +568,8 @@
return (BPoints->n_points);
}
-int
-get_pruned_isle_points (Map, isle, BPoints)
- struct Map_info *Map;
- int isle;
- struct line_pnts *BPoints;
+int
+get_pruned_isle_points (struct Map_info *Map, int isle, struct line_pnts *BPoints)
{
register int i, line;
int start, end, to, from, inc;
More information about the grass-dev
mailing list