[GRASS-dev] [INFO] Line simplification algorithm
tlaronde at polynum.com
tlaronde at polynum.com
Tue Nov 21 09:59:19 EST 2006
On Tue, Nov 21, 2006 at 03:31:52PM +1300, Hamish wrote:
>
> "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
I'm talking about the main pruning routing dig_prune() to be found in
VECTORLIB.
Since dig_prune dated back 1987, and was unchanged during 10 years it is
unlikely that it was not correctly doing its work.
BTW, since I'm using memories from 3 years ago (I haven't used GRASS
GPL since 5.0) am I confusing v.prune with another module (v.clean)?
But I don't think so.
I think I used v.rm.dangles (or something like that), v.prune (and drop
it since the curves disappeared) and v.spag (that was not doing as much
work as it shall do).
I attach the original version of dig_prune.
> 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
> -------------------------------------------------------------------
That's puzzled me, since the 4.1 version is short, commented and make
sense. Does something happen between 4.1 and 4.2, that is without a
track in a CVS (not here?). Did someone mess with the code previously to
this change?
The "not only for linux version" is strange: there is nothing system
dependant in the original.
>
> You don't have to name the guilty if you don't want to :), but who do
> you consider to be really good?
>
OK :) That's quite simple: you can trust the CERL (Michael Shapiro,
James Westervelt, Dave Gerdes). There has been programs in alpha coming
from outside the CERL. In this case...
BTW, some alpha programs for vector (namely v.spag and v.cutter) were
buggy. But this was not due to the incompetence of the programmer (well:
developer). This is an historic consequence: Dave Gerdes was in a hurry
just before the team was disbanded, and the alpha remained alpha.
FWIW, I have fixed v.spag(1). There were two main problems:
1) To be sure that the program ends, you must ensure that a line is
whether unchanged but never touched again, or split but resulting in
another line (not the same or you will enter an infinite loop). There
were cases where lines were cut in two new lines, a null length one, and
the remaining, that is the same as the previous;
2) Nodes created were not correctly snapped to existing node, leading to
the creation of several nodes with the exact same coordinates =>
visually, lines were connected, but not to the same node so without
snapping nodes areas could not be created.
There were, too, missing handling of colinear cases.
v.cutter(1) is an ambitious complex monster. I have chosen to replace it
with an extended version of v.spag(1) , since v.spag is trying to cut a
map agains itself when you think at it that is it's a particular case of
v.cutter.
Cheers,
--
Thierry Laronde (Alceste) <tlaronde +AT+ polynum +dot+ com>
http://www.kergis.com/
Key fingerprint = 0FF7 E906 FBAF FE95 FD89 250D 52B1 AE95 6006 F40C
-------------- next part --------------
/* $Id: prune.c,v 1.7 2006/10/31 14:23:31 tlaronde Exp $
* Copyright 2004 Thierry LARONDE <tlaronde at polynum.com>
* All rights reserved.
*
* See the COPYRIGHTS at the root of the source directory.
*
*!!!THIS SOFTWARE IS PROVIDED ``AS IS'' WITHOUT ANY WARRANTIES!!!
* USE IT AT YOUR OWN RISK
********************************************************************/
/* @(#)prune.c 2.1 6/26/87 */
/*
* This subroutine resamples a dense string of x,y coordinates to produce a
* set of coordinates that approaches hand digitizing. That is, the density
* of points is very low on straight lines, and highest on tight curves.
*
* xarr and yarr - double precision sets of coordinates. num - the total
* number of points in the set. thresh - the distance that a string must
* wander from a straight line before another point is selected.
*/
#include <float.h>
#include <stdio.h>
#include <math.h>
#include "vector.h"
int
dig_prune(struct line_pnts * points, double thresh)
{
double *ox, *oy, *nx, *ny;
double thresh_sq;
double half_thresh;
double half_thresh_sq;
double cur_x, cur_y;
double dx, dy;
double last_x, last_y;
double a1, b1;
double a2, b2;
double int_x, int_y;
double tst_dist;
int o_num;
int n_num;
int at_num;
/* Automatic variables must be initialize before used */
/*
* But here dx and dy are set before used since the while loop is
* true OR the function returns without using them
*/
dx = dy = 0.0;
if (points->n_points <= 2)
return (points->n_points);
ox = points->x;
oy = points->y;
nx = points->x;
ny = points->y;
o_num = points->n_points;
n_num = 0;
half_thresh = thresh / 2;
half_thresh_sq = half_thresh * half_thresh;
thresh_sq = thresh * thresh;
/* First point is retained */
cur_x = *ox++;
cur_y = *oy++;
n_num++;
nx++;
ny++;
at_num = 1;
/* Search for first point > half_thresh from current point */
while (at_num < o_num) {
dx = *ox - cur_x;
dy = *oy - cur_y;
if (dx * dx + dy * dy > half_thresh_sq)
break;
at_num++;
ox++;
oy++;
}
/* If all points in thresh just output last point */
if (at_num == o_num) {
n_num = 2;
*nx = *(ox - 1);
*ny = *(oy - 1);
return (n_num);
}
/* if all but last point in thresh, just output last point */
if (at_num == o_num - 1) {
n_num = 2;
*nx = *(ox);
*ny = *(oy);
return (n_num);
}
/* calculate line equation coefficients */
if (dx == 0.0)
a1 = DBL_MAX;
else
a1 = dy / dx; /* y = ax + b */
b1 = *oy - a1 * *ox;
last_x = *ox++;
last_y = *oy++;
at_num++;
/*
* Main loop. Look for first point outside of thresh limit around line;
* save last point; look for next point to define new line; repeat
*/
while (at_num < o_num) {
/* Look for first point outside of thresh limit */
if (a1 == 0.0)
a2 = DBL_MAX;
else
a2 = -1 / a1;
while (at_num < o_num) {
/* Calc. coef. of perp. line through current point */
b2 = *oy - a2 * *ox;
/* Solve two equations for common point */
int_x = (b2 - b1) / (a1 - a2);
int_y = a2 * int_x + b2;
/* Calc. dist between intersection and cur. point */
dx = *ox - int_x;
dy = *oy - int_y;
tst_dist = dx * dx + dy * dy;
/* Is this > thresh (within thresh band) ? */
if (tst_dist > thresh_sq)
break;
/*
* If not, this point is next candidate for next saved point
*/
last_x = *ox++;
last_y = *oy++;
at_num++;
}
*nx++ = last_x;
*ny++ = last_y;
n_num++;
if (at_num == o_num)
return (n_num);
cur_x = last_x;
cur_y = last_y;
/* Search for next point > half_thresh from current point */
while (at_num < o_num) {
dx = *ox - cur_x;
dy = *oy - cur_y;
if (dx * dx + dy * dy > half_thresh_sq)
break;
at_num++;
ox++;
oy++;
}
if (at_num == o_num) {
n_num++;
*nx = *(ox - 1);
*ny = *(oy - 1);
return (n_num);
}
if (dx == 0.0)
a1 = DBL_MAX;
else
a1 = dy / dx; /* y = ax + b */
b1 = cur_y - a1 * cur_x;
last_x = *ox;
last_y = *oy;
}
return (n_num);
}
More information about the grass-dev
mailing list