[GRASS-SVN] r56264 - grass/trunk/vector/v.generalize
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 15 07:25:49 PDT 2013
Author: mmetz
Date: 2013-05-15 07:25:44 -0700 (Wed, 15 May 2013)
New Revision: 56264
Modified:
grass/trunk/vector/v.generalize/displacement.c
grass/trunk/vector/v.generalize/main.c
grass/trunk/vector/v.generalize/matrix.c
grass/trunk/vector/v.generalize/matrix.h
grass/trunk/vector/v.generalize/operators.h
grass/trunk/vector/v.generalize/point.c
grass/trunk/vector/v.generalize/point.h
grass/trunk/vector/v.generalize/simplification.c
grass/trunk/vector/v.generalize/smoothing.c
Log:
v.generalize: add loop support for boyle, chaiken, distance_weighting, hermite, snakes
Modified: grass/trunk/vector/v.generalize/displacement.c
===================================================================
--- grass/trunk/vector/v.generalize/displacement.c 2013-05-15 13:35:12 UTC (rev 56263)
+++ grass/trunk/vector/v.generalize/displacement.c 2013-05-15 14:25:44 UTC (rev 56264)
@@ -90,7 +90,7 @@
int q, findex;
POINT cur;
- point_assign(Points, j, with_z, &cur);
+ point_assign(Points, j, with_z, &cur, 0);
/* check whether we alerady have point with the same
* coordinates */
findex = pindex;
@@ -102,12 +102,12 @@
point_index[index] = findex;
if (findex == pindex) {
- point_assign(Points, j, with_z, &pset[pindex]);
+ point_assign(Points, j, with_z, &pset[pindex], 0);
pindex++;
}
first[index] = (j == 0);
line_index[index] = i;
- point_assign(Points, j, with_z, &parray[index]);
+ point_assign(Points, j, with_z, &parray[index], 0);
index++;
}
}
@@ -197,7 +197,7 @@
/*calculate the inverse */
G_message(_("Inverting matrix..."));
- if (!matrix_inverse(k, &kinv, 1))
+ if (!matrix_inverse(&k, &kinv, 1))
G_fatal_error(_("Unable to calculate the inverse matrix"));
G_percent_reset();
@@ -213,8 +213,8 @@
matrix_mult_scalar(0.0, &dx_old);
matrix_mult_scalar(0.0, &dy_old);
- matrix_add(dx_old, dx, &dx_old);
- matrix_add(dy_old, dy, &dy_old);
+ matrix_add(&dx_old, &dx, &dx_old);
+ matrix_add(&dy_old, &dy, &dy_old);
/* calculate force vectors */
for (i = 0; i < index; i++) {
@@ -271,11 +271,11 @@
matrix_mult_scalar(gama, &dx);
matrix_mult_scalar(gama, &dy);
- matrix_add(dx, fx, &fx);
- matrix_add(dy, fy, &fy);
+ matrix_add(&dx, &fx, &fx);
+ matrix_add(&dy, &fy, &fy);
- matrix_mult(kinv, fx, &dx);
- matrix_mult(kinv, fy, &dy);
+ matrix_mult(&kinv, &fx, &dx);
+ matrix_mult(&kinv, &fy, &dy);
for (i = 0; i < index; i++) {
if (point_index[i] == -1)
@@ -286,7 +286,6 @@
dy.a[point_index[i]][0] - dy_old.a[point_index[i]][0];
}
-
}
index = 0;
for (i = 1; i <= n_lines; i++) {
@@ -313,14 +312,14 @@
G_free(need);
G_free(sel);
G_free(tmp_index);
- matrix_free(k);
- matrix_free(kinv);
- matrix_free(dx);
- matrix_free(dy);
- matrix_free(fx);
- matrix_free(fy);
- matrix_free(dx_old);
- matrix_free(dy_old);
+ matrix_free(&k);
+ matrix_free(&kinv);
+ matrix_free(&dx);
+ matrix_free(&dy);
+ matrix_free(&fx);
+ matrix_free(&fy);
+ matrix_free(&dx_old);
+ matrix_free(&dy_old);
return 0;
}
Modified: grass/trunk/vector/v.generalize/main.c
===================================================================
--- grass/trunk/vector/v.generalize/main.c 2013-05-15 13:35:12 UTC (rev 56263)
+++ grass/trunk/vector/v.generalize/main.c 2013-05-15 14:25:44 UTC (rev 56264)
@@ -228,7 +228,8 @@
loop_support_flag = G_define_flag();
loop_support_flag->key = 'l';
- loop_support_flag->description = _("Loop support");
+ loop_support_flag->label = _("Loop support");
+ loop_support_flag->description = _("Modify end points of lines forming a closed loop");
notab_flag = G_define_standard_flag(G_FLG_V_TABLE);
notab_flag->description = _("Do not copy attributes");
@@ -486,39 +487,45 @@
reumann_witkam(Points, thresh, with_z);
break;
case BOYLE:
- boyle(Points, look_ahead, with_z);
+ boyle(Points, look_ahead, loop_support, with_z);
break;
case SLIDING_AVERAGING:
sliding_averaging(Points, slide, look_ahead, loop_support, with_z);
break;
case DISTANCE_WEIGHTING:
- distance_weighting(Points, slide, look_ahead, with_z);
+ distance_weighting(Points, slide, look_ahead, loop_support, with_z);
break;
case CHAIKEN:
- chaiken(Points, thresh, with_z);
+ chaiken(Points, thresh, loop_support, with_z);
break;
case HERMITE:
- hermite(Points, thresh, angle_thresh, with_z);
+ hermite(Points, thresh, angle_thresh, loop_support, with_z);
break;
case SNAKES:
- snakes(Points, alpha, beta, with_z);
+ snakes(Points, alpha, beta, loop_support, with_z);
break;
}
}
-
- if (method != SLIDING_AVERAGING || loop_support == 0){
- /* safety check, BUG in method if not passed */
- if (APoints->x[0] != Points->x[0] ||
- APoints->y[0] != Points->y[0] ||
- APoints->z[0] != Points->z[0])
- G_fatal_error(_("Method '%s' did not preserve first point"), method_opt->answer);
-
- if (APoints->x[APoints->n_points - 1] != Points->x[Points->n_points - 1] ||
- APoints->y[APoints->n_points - 1] != Points->y[Points->n_points - 1] ||
- APoints->z[APoints->n_points - 1] != Points->z[Points->n_points - 1])
- G_fatal_error(_("Method '%s' did not preserve last point"), method_opt->answer);
- }
+ if (loop_support == 0) {
+ /* safety check, BUG in method if not passed */
+ if (APoints->x[0] != Points->x[0] ||
+ APoints->y[0] != Points->y[0] ||
+ APoints->z[0] != Points->z[0])
+ G_fatal_error(_("Method '%s' did not preserve first point"), method_opt->answer);
+
+ if (APoints->x[APoints->n_points - 1] != Points->x[Points->n_points - 1] ||
+ APoints->y[APoints->n_points - 1] != Points->y[Points->n_points - 1] ||
+ APoints->z[APoints->n_points - 1] != Points->z[Points->n_points - 1])
+ G_fatal_error(_("Method '%s' did not preserve last point"), method_opt->answer);
+ }
+ else {
+ /* safety check, BUG in method if not passed */
+ if (Points->x[0] != Points->x[Points->n_points - 1] ||
+ Points->y[0] != Points->y[Points->n_points - 1] ||
+ Points->z[0] != Points->z[Points->n_points - 1])
+ G_fatal_error(_("Method '%s' did not preserve loop"), method_opt->answer);
+ }
Vect_line_prune(Points);
Modified: grass/trunk/vector/v.generalize/matrix.c
===================================================================
--- grass/trunk/vector/v.generalize/matrix.c 2013-05-15 13:35:12 UTC (rev 56263)
+++ grass/trunk/vector/v.generalize/matrix.c 2013-05-15 14:25:44 UTC (rev 56264)
@@ -22,7 +22,7 @@
#include <grass/glocale.h>
#include "matrix.h"
-int matrix_init(int rows, int cols, MATRIX * res)
+int matrix_init(int rows, int cols, MATRIX *res)
{
int i, j;
@@ -46,19 +46,19 @@
return 1;
}
-void matrix_free(MATRIX m)
+void matrix_free(MATRIX *m)
{
int i;
- for (i = 0; i < m.rows; i++)
- G_free(m.a[i]);
- G_free(m.a);
+ for (i = 0; i < m->rows; i++)
+ G_free(m->a[i]);
+ G_free(m->a);
return;
}
-int matrix_mult(MATRIX a, MATRIX b, MATRIX * res)
+int matrix_mult(MATRIX *a, MATRIX *b, MATRIX *res)
{
- if (a.cols != b.rows)
+ if (a->cols != b->rows)
return 0;
/*if (!matrix_init(a.rows, b.cols, res))
@@ -67,17 +67,17 @@
int i, j, k;
- for (i = 0; i < a.rows; i++)
- for (j = 0; j < b.cols; j++) {
+ for (i = 0; i < a->rows; i++)
+ for (j = 0; j < b->cols; j++) {
res->a[i][j] = 0;
- for (k = 0; k < a.cols; k++)
- res->a[i][j] += a.a[i][k] * b.a[k][j];
+ for (k = 0; k < a->cols; k++)
+ res->a[i][j] += a->a[i][k] * b->a[k][j];
}
return 1;
}
-int matrix_add_identity(double s, MATRIX * m)
+int matrix_add_identity(double s, MATRIX *m)
{
if (m->rows != m->cols)
return 0;
@@ -92,7 +92,7 @@
/* three following functions implements elementary row operations on matrices */
/* auxialiry function for matrix_inverse, swaps two rows of given matrix */
-void matrix_swap_rows(int x, int y, MATRIX * m)
+void matrix_swap_rows(int x, int y, MATRIX *m)
{
int i;
@@ -108,7 +108,7 @@
/* auxiliary function for matrix_inverse, multiplies row of a matrix by
* a scalar */
-void matrix_row_scalar(int row, double s, MATRIX * m)
+void matrix_row_scalar(int row, double s, MATRIX *m)
{
int i;
@@ -121,7 +121,7 @@
* one row to another.
* i.e row[ra] = row[ra] + row[rb] * s;
*/
-void matrix_row_add_multiple(int ra, int rb, double s, MATRIX * m)
+void matrix_row_add_multiple(int ra, int rb, double s, MATRIX *m)
{
int i;
@@ -131,26 +131,24 @@
}
/* TODO: don't test directly equality to zero */
-int matrix_inverse(MATRIX a, MATRIX * res, int percents)
-{;
+int matrix_inverse(MATRIX *a, MATRIX *res, int percents)
+{
+ int i, j;
/* not a square matrix */
- if (a.rows != a.cols)
+ if (a->rows != a->cols)
return 0;
- int i, j;
-
/* initialize output matrix to the identity matrix */
- if (!matrix_init(a.rows, a.rows, res)) {
+ if (!matrix_init(a->rows, a->rows, res)) {
G_fatal_error(_("Out of memory"));
return 0;
}
- for (i = 0; i < a.rows; i++) {
- memset(res->a[i], 0, sizeof(double) * a.cols);
+ for (i = 0; i < a->rows; i++) {
+ memset(res->a[i], 0, sizeof(double) * a->cols);
res->a[i][i] = 1;
}
-
/* in order to obtain the inverse of a matrix, we run
* gauss elimination on the matrix and each time we apply
* elementary row operation on this matrix, we apply the
@@ -159,7 +157,7 @@
* is row equivalent to the identity matrix.
*/
- int n = a.rows;
+ int n = a->rows;
if (percents)
G_percent_reset();
@@ -170,27 +168,27 @@
if (percents)
G_percent(i, n, 1);
for (j = i; j < n; j++) {
- if (a.a[j][i] != 0) { /* need to change this row to something */
+ if (a->a[j][i] != 0) { /* need to change this row to something */
found = 1; /* more sensible */
- matrix_swap_rows(i, j, &a);
+ matrix_swap_rows(i, j, a);
matrix_swap_rows(i, j, res);
break;
}
}
if (!found)
return 0;
- double c = (double)1.0 / a.a[i][i];
+ double c = (double)1.0 / a->a[i][i];
- matrix_row_scalar(i, c, &a);
+ matrix_row_scalar(i, c, a);
matrix_row_scalar(i, c, res);
for (j = 0; j < n; j++) {
if (i == j)
continue;
- double c = -a.a[j][i];
+ double c = -a->a[j][i];
if (c == 0.0)
continue;
- matrix_row_add_multiple(j, i, c, &a);
+ matrix_row_add_multiple(j, i, c, a);
matrix_row_add_multiple(j, i, c, res);
}
}
@@ -198,7 +196,7 @@
return 1;
}
-void matrix_mult_scalar(double s, MATRIX * m)
+void matrix_mult_scalar(double s, MATRIX *m)
{
int i, j;
@@ -207,25 +205,25 @@
m->a[i][j] *= s;
}
-void matrix_add(MATRIX a, MATRIX b, MATRIX * res)
+void matrix_add(MATRIX *a, MATRIX *b, MATRIX *res)
{
int i, j;
for (i = 0; i < res->rows; i++)
for (j = 0; j < res->cols; j++)
- res->a[i][j] = a.a[i][j] + b.a[i][j];
+ res->a[i][j] = a->a[i][j] + b->a[i][j];
}
-void matrix_print(MATRIX a)
+void matrix_print(MATRIX *a)
{
int i, j;
- for (i = 0; i < a.rows; i++) {
+ for (i = 0; i < a->rows; i++) {
double s = 0;
- for (j = 0; j < a.cols; j++) {
- printf("%.3lf ", a.a[i][j]);
- s += a.a[i][j];
+ for (j = 0; j < a->cols; j++) {
+ printf("%.3lf ", a->a[i][j]);
+ s += a->a[i][j];
}
printf("|%.5lf\n", s);
}
Modified: grass/trunk/vector/v.generalize/matrix.h
===================================================================
--- grass/trunk/vector/v.generalize/matrix.h 2013-05-15 13:35:12 UTC (rev 56263)
+++ grass/trunk/vector/v.generalize/matrix.h 2013-05-15 14:25:44 UTC (rev 56264)
@@ -29,32 +29,32 @@
} MATRIX;
/* return 1 on success, 0 on error (out of memory) */
-extern int matrix_init(int rows, int cols, MATRIX * res);
+extern int matrix_init(int rows, int cols, MATRIX *res);
/* free the memory occupied by the values of m */
-extern void matrix_free(MATRIX m);
+extern void matrix_free(MATRIX *m);
/* multiply two matrices, Return 1 on success, 0 on failure.
* return value 0 means - bad dimensions */
-extern int matrix_mult(MATRIX a, MATRIX b, MATRIX * res);
+extern int matrix_mult(MATRIX *a, MATRIX *b, MATRIX *res);
/* adds a multiple of the identity matrix to the given matrix
* M = M + s * Id. Returns 1 on success, 0 otherwise */
-extern int matrix_add_identity(double s, MATRIX * m);
+extern int matrix_add_identity(double s, MATRIX *m);
/* calculate the inverse of given (square) matrix. Returns 0 if
* the matrix is not invertible or if an error occurs.
* percents indicates whether we want to show the progress of
* computation
* Otherwise it returns 1 */
-extern int matrix_inverse(MATRIX a, MATRIX * res, int percents);
+extern int matrix_inverse(MATRIX *a, MATRIX *res, int percents);
/* multiplies matrix by a scalar */
-extern void matrix_mult_scalar(double s, MATRIX * m);
+extern void matrix_mult_scalar(double s, MATRIX *m);
/* res = a + b. Does not cheack the dimensions */
-extern void matrix_add(MATRIX a, MATRIX b, MATRIX * res);
+extern void matrix_add(MATRIX *a, MATRIX *b, MATRIX *res);
/* debug function */
-extern void matrix_print(MATRIX a);
+extern void matrix_print(MATRIX *a);
#endif
Modified: grass/trunk/vector/v.generalize/operators.h
===================================================================
--- grass/trunk/vector/v.generalize/operators.h 2013-05-15 13:35:12 UTC (rev 56263)
+++ grass/trunk/vector/v.generalize/operators.h 2013-05-15 14:25:44 UTC (rev 56264)
@@ -10,15 +10,15 @@
double reduction, int with_z);
/* smoothing.c */
-int boyle(struct line_pnts *Points, int look_ahead, int with_z);
+int boyle(struct line_pnts *Points, int look_ahead, int loop_support, int with_z);
int sliding_averaging(struct line_pnts *Points, double slide, int look_ahead,
int loop_support, int with_z);
int distance_weighting(struct line_pnts *Points, double slide, int look_ahead,
- int with_z);
-int chaiken(struct line_pnts *Points, double thresh, int with_z);
+ int loop_support, int with_z);
+int chaiken(struct line_pnts *Points, double thresh, int loop_support, int with_z);
int hermite(struct line_pnts *Points, double step, double angle_thresh,
- int with_z);
-int snakes(struct line_pnts *Points, double alpha, double beta, int with_z);
+ int loop_support, int with_z);
+int snakes(struct line_pnts *Points, double alpha, double beta, int loop_support, int with_z);
/* network.c */
int graph_generalization(struct Map_info *In, struct Map_info *Out,
Modified: grass/trunk/vector/v.generalize/point.c
===================================================================
--- grass/trunk/vector/v.generalize/point.c 2013-05-15 13:35:12 UTC (rev 56263)
+++ grass/trunk/vector/v.generalize/point.c 2013-05-15 14:25:44 UTC (rev 56264)
@@ -50,23 +50,12 @@
}
inline void point_assign(struct line_pnts *Points, int index, int with_z,
- POINT * res)
+ POINT * res, int is_loop)
{
- res->x = Points->x[index];
- res->y = Points->y[index];
- if (with_z) {
- res->z = Points->z[index];
+ if (is_loop) {
+ while (index >= Points->n_points - 1)
+ index -= Points->n_points - 1;
}
- else {
- res->z = 0;
- }
- return;
-}
-
-inline void point_assign_loop(struct line_pnts *Points, int index, int with_z,
- POINT * res)
-{
- index = index % Points->n_points;
res->x = Points->x[index];
res->y = Points->y[index];
if (with_z) {
Modified: grass/trunk/vector/v.generalize/point.h
===================================================================
--- grass/trunk/vector/v.generalize/point.h 2013-05-15 13:35:12 UTC (rev 56263)
+++ grass/trunk/vector/v.generalize/point.h 2013-05-15 14:25:44 UTC (rev 56264)
@@ -49,13 +49,12 @@
* if with z = 0 then res.z = 0
*/
extern inline void point_assign(struct line_pnts *Points, int index,
- int with_z, POINT * res);
+ int with_z, POINT * res, int is_loop);
/* assign point Points[index] to the res
* if with z = 0 then res.z = 0
* loop to infinite
*/
-extern inline void point_assign_loop(struct line_pnts *Points, int index,
- int with_z, POINT * res);
+
/* res = k * a */
extern inline void point_scalar(POINT a, double k, POINT * res);
Modified: grass/trunk/vector/v.generalize/simplification.c
===================================================================
--- grass/trunk/vector/v.generalize/simplification.c 2013-05-15 13:35:12 UTC (rev 56263)
+++ grass/trunk/vector/v.generalize/simplification.c 2013-05-15 14:25:44 UTC (rev 56264)
@@ -239,15 +239,15 @@
count = 1;
- point_assign(Points, 0, with_z, &x1);
- point_assign(Points, 1, with_z, &x2);
+ point_assign(Points, 0, with_z, &x1, 0);
+ point_assign(Points, 1, with_z, &x2, 0);
point_subtract(x2, x1, &sub);
subd = point_dist2(sub);
for (i = 2; i < n; i++) {
- point_assign(Points, i, with_z, &x0);
+ point_assign(Points, i, with_z, &x0, 0);
point_subtract(x1, x0, &diff);
diffd = point_dist2(diff);
sp = point_dot(diff, sub);
@@ -256,8 +256,8 @@
* all variables which do not change for each line-point calculation */
if (dist > thresh) {
- point_assign(Points, i - 1, with_z, &x1);
- point_assign(Points, i, with_z, &x2);
+ point_assign(Points, i - 1, with_z, &x1, 0);
+ point_assign(Points, i, with_z, &x2, 0);
point_subtract(x2, x1, &sub);
subd = point_dist2(sub);
Modified: grass/trunk/vector/v.generalize/smoothing.c
===================================================================
--- grass/trunk/vector/v.generalize/smoothing.c 2013-05-15 13:35:12 UTC (rev 56263)
+++ grass/trunk/vector/v.generalize/smoothing.c 2013-05-15 14:25:44 UTC (rev 56264)
@@ -30,43 +30,78 @@
/* boyle's forward looking algorithm
* return the number of points in the result = Points->n_points
*/
-int boyle(struct line_pnts *Points, int look_ahead, int with_z)
+int boyle(struct line_pnts *Points, int look_ahead, int loop_support, int with_z)
{
- POINT last, npoint, ppoint;
+ POINT last, ppoint;
+ POINT *res;
int next, n, i, p;
double c1, c2;
+ int is_loop, count;
n = Points->n_points;
/* if look_ahead is too small or line too short, there's nothing
* to smooth */
- if (look_ahead < 2 || look_ahead > n) {
+ if (look_ahead < 2 || look_ahead >= n) {
return n;
}
- point_assign(Points, 0, with_z, &last);
+ count = n - 2;
+ is_loop = 0;
+
+ /* is it loop ? */
+ if (Points->x[0] == Points->x[n - 1] &&
+ Points->y[0] == Points->y[n - 1] &&
+ (Points->z[0] == Points->z[n - 1] || with_z == 0) && loop_support){
+ is_loop = 1;
+ count = n + look_ahead;
+ }
+
+ res = G_malloc(sizeof(POINT) * n);
+
+ point_assign(Points, 0, with_z, &last, 0);
+ res[0] = last;
c1 = (double)1 / (double)(look_ahead - 1);
c2 = (double)1 - c1;
next = 1;
- for (i = 0; i < n - 2; i++) {
+ for (i = 0; i < count; i++) {
p = i + look_ahead;
- if (p >= n)
+ if (!is_loop && p >= n)
p = n - 1;
- point_assign(Points, p, with_z, &ppoint);
+ point_assign(Points, p, with_z, &ppoint, is_loop);
point_scalar(ppoint, c1, &ppoint);
point_scalar(last, c2, &last);
- point_add(last, ppoint, &npoint);
- Points->x[next] = npoint.x;
- Points->y[next] = npoint.y;
- Points->z[next] = npoint.z;
+ point_add(last, ppoint, &res[next]);
+ /* original: use smoothed point as new last point */
+ /* last = res[next]; */
next++;
- last = npoint;
+ if (is_loop) {
+ while (next >= n - 1)
+ next -= n - 1;
+ }
+ /* new: smooth with original points */
+ point_assign(Points, p, with_z, &last, is_loop);
}
- points_copy_last(Points, next);
+ for (i = 1; i < n - 1; i++) {
+ Points->x[i] = res[i].x;
+ Points->y[i] = res[i].y;
+ Points->z[i] = res[i].z;
+ }
+ if (is_loop) {
+ Points->x[0] = res[0].x;
+ Points->y[0] = res[0].y;
+ Points->z[0] = res[0].z;
+
+ Points->x[n - 1] = res[0].x;
+ Points->y[n - 1] = res[0].y;
+ Points->z[n - 1] = res[0].z;
+ }
+
+ G_free(res);
+
return Points->n_points;
-
}
/* mcmaster's sliding averaging algorithm. Return the number of points
@@ -83,16 +118,16 @@
POINT *res;
int is_loop, count;
- is_loop=0;
+ is_loop = 0;
n = Points->n_points;
half = look_ahead / 2;
count = n - half;
- /* is it loop ?*/
- if (Points->x[0] == Points->x[n-1]
- && Points->y[0] == Points->y[n-1]
- && (Points->z[0] == Points->z[n-1] || with_z == 0) && loop_support){
+ /* is it loop ? */
+ if (Points->x[0] == Points->x[n - 1] &&
+ Points->y[0] == Points->y[n - 1] &&
+ (Points->z[0] == Points->z[n - 1] || with_z == 0) && loop_support){
is_loop = 1;
count = n + half;
}
@@ -102,7 +137,7 @@
return n;
}
- if (look_ahead >= n || look_ahead == 1)
+ if (look_ahead >= n || look_ahead < 2)
return n;
res = G_malloc(sizeof(POINT) * (n + half));
@@ -113,40 +148,46 @@
sc = (double)1.0 / (double)look_ahead;
- point_assign(Points, 0, with_z, &p);
+ point_assign(Points, 0, with_z, &p, 0);
for (i = 1; i < look_ahead; i++) {
- point_assign(Points, i, with_z, &tmp);
+ point_assign(Points, i, with_z, &tmp, 0);
point_add(p, tmp, &p);
}
/* and calculate the average of remaining points */
for (i = half; i < count; i++) {
- point_assign_loop(Points, i, with_z, &s);
+ point_assign(Points, i, with_z, &s, is_loop);
point_scalar(s, 1.0 - slide, &s);
point_scalar(p, sc * slide, &tmp);
point_add(tmp, s, &res[i]);
- if ((i + half + 1 < n) || is_loop ) {
- point_assign_loop(Points, i - half, with_z, &tmp);
+ if ((i + half + 1 < n) || is_loop) {
+ point_assign(Points, i - half, with_z, &tmp, is_loop);
point_subtract(p, tmp, &p);
- point_assign_loop(Points, i + half + 1, with_z, &tmp);
+ point_assign(Points, i + half + 1, with_z, &tmp, is_loop);
point_add(p, tmp, &p);
}
}
- for (i = half; i < count; i++) {
- Points->x[i] = res[i].x;
- Points->y[i] = res[i].y;
- Points->z[i] = res[i].z;
- }
-
if (is_loop) {
for (i = 0; i < half; i++) {
- Points->x[i] = res[n + i - 1].x;
- Points->y[i] = res[n + i - 1].y;
- Points->z[i] = res[n + i - 1].z;
+ Points->x[i] = res[n + i - 1].x;
+ Points->y[i] = res[n + i - 1].y;
+ Points->z[i] = res[n + i - 1].z;
}
+ for (i = half; i < n; i++) {
+ Points->x[i] = res[i].x;
+ Points->y[i] = res[i].y;
+ Points->z[i] = res[i].z;
+ }
}
+ else {
+ for (i = half; i < n - half; i++) {
+ Points->x[i] = res[i].x;
+ Points->y[i] = res[i].y;
+ Points->z[i] = res[i].z;
+ }
+ }
G_free(res);
@@ -156,40 +197,53 @@
/* mcmaster's distance weighting algorithm. Return the number
* of points in the output line which equals to Points->n_points */
int distance_weighting(struct line_pnts *Points, double slide, int look_ahead,
- int with_z)
+ int loop_support, int with_z)
{
POINT p, c, s, tmp;
int n, i, half, j;
double dists, d;
POINT *res;
+ int is_loop, count;
n = Points->n_points;
+ if (look_ahead >= n || look_ahead < 1)
+ return n;
+
+ half = look_ahead / 2;
+ count = n - half;
+
+ /* is it loop ? */
+ is_loop = 0;
+ if (Points->x[0] == Points->x[n - 1] &&
+ Points->y[0] == Points->y[n - 1] &&
+ (Points->z[0] == Points->z[n - 1] || with_z == 0) && loop_support){
+ is_loop = 1;
+ count = n + half;
+ }
+
if (look_ahead % 2 == 0) {
G_fatal_error(_("Look ahead parameter must be odd"));
return n;
}
- res = (POINT *) G_malloc(sizeof(POINT) * n);
+ res = (POINT *) G_malloc(sizeof(POINT) * (n + half));
if (!res) {
G_fatal_error(_("Out of memory"));
return n;
}
- point_assign(Points, 0, with_z, &res[0]);
+ point_assign(Points, 0, with_z, &res[0], 0);
- half = look_ahead / 2;
-
- for (i = half; i + half < n; i++) {
- point_assign(Points, i, with_z, &c);
+ for (i = half; i < count; i++) {
+ point_assign(Points, i, with_z, &c, is_loop);
s.x = s.y = s.z = 0;
dists = 0;
-
for (j = i - half; j <= i + half; j++) {
if (j == i)
continue;
- point_assign(Points, j, with_z, &p);
+ point_assign(Points, j, with_z, &p, is_loop);
d = point_dist(p, c);
if (d < GRASS_EPSILON)
continue;
@@ -200,30 +254,51 @@
s.y += tmp.y;
s.z += tmp.z;
}
- point_scalar(s, slide / dists, &tmp);
- point_scalar(c, (double)1.0 - slide, &s);
- point_add(s, tmp, &res[i]);
+ if (dists < GRASS_EPSILON) {
+ point_add(c, s, &res[i]);
+ }
+ else {
+ point_scalar(s, slide / dists, &tmp);
+ point_scalar(c, (double)1.0 - slide, &s);
+ point_add(s, tmp, &res[i]);
+ }
}
- for (i = half; i + half < n; i++) {
- Points->x[i] = res[i].x;
- Points->y[i] = res[i].y;
- Points->z[i] = res[i].z;
+ if (is_loop) {
+ for (i = 0; i < half; i++) {
+ Points->x[i] = res[n + i - 1].x;
+ Points->y[i] = res[n + i - 1].y;
+ Points->z[i] = res[n + i - 1].z;
+ }
+ for (i = half; i < n; i++) {
+ Points->x[i] = res[i].x;
+ Points->y[i] = res[i].y;
+ Points->z[i] = res[i].z;
+ }
}
+ else {
+ for (i = half; i < n - half; i++) {
+ Points->x[i] = res[i].x;
+ Points->y[i] = res[i].y;
+ Points->z[i] = res[i].z;
+ }
+ }
G_free(res);
+
return Points->n_points;
}
/* Chaiken's algorithm. Return the number of points in smoothed line
*/
-int chaiken(struct line_pnts *Points, double thresh, int with_z)
+int chaiken(struct line_pnts *Points, double thresh, int loop_support, int with_z)
{
int n, i;
POINT_LIST head, *cur;
POINT p0, p1, p2, m1, tmp;
+ int is_loop;
n = Points->n_points;
@@ -231,47 +306,89 @@
if (n < 3)
return n;
+ is_loop = 0;
+ /* is it loop ? */
+ if (Points->x[0] == Points->x[n - 1] &&
+ Points->y[0] == Points->y[n - 1] &&
+ (Points->z[0] == Points->z[n - 1] || with_z == 0) && loop_support){
+ is_loop = 1;
+ }
+
thresh *= thresh;
head.next = NULL;
cur = &head;
- point_assign(Points, 0, with_z, &p0);
- point_list_add(cur, p0); /* always keep first point */
- cur = cur->next;
+ point_assign(Points, 0, with_z, &p0, 0);
- for (i = 2; i <= n; i++) {
- if (i == n)
- point_assign(Points, i - 1, with_z, &p2);
- else
- point_assign(Points, i, with_z, &p2);
- point_assign(Points, i - 1, with_z, &p1);
+ if (!is_loop) {
+ point_list_add(cur, p0); /* always keep first point */
+ cur = cur->next;
- while (1) {
- point_add(p1, p2, &tmp);
- point_scalar(tmp, 0.5, &m1);
+ for (i = 2; i <= n; i++) {
+ if (i == n)
+ point_assign(Points, i - 1, with_z, &p2, 0);
+ else
+ point_assign(Points, i, with_z, &p2, 0);
+ point_assign(Points, i - 1, with_z, &p1, 0);
- point_list_add(cur, m1);
+ while (1) {
+ point_add(p1, p2, &tmp);
+ point_scalar(tmp, 0.5, &m1);
- if (point_dist_square(p0, m1) > thresh) {
- point_add(p1, m1, &tmp); /* need to refine the partition */
- point_scalar(tmp, 0.5, &p2);
- point_add(p1, p0, &tmp);
- point_scalar(tmp, 0.5, &p1);
+ point_list_add(cur, m1);
+
+ if (point_dist_square(p0, m1) > thresh) {
+ point_add(p1, m1, &tmp); /* need to refine the partition */
+ point_scalar(tmp, 0.5, &p2);
+ point_add(p1, p0, &tmp);
+ point_scalar(tmp, 0.5, &p1);
+ }
+ else {
+ break; /* good approximation */
+ }
}
- else {
- break; /* good approximation */
- }
+
+ while (cur->next != NULL)
+ cur = cur->next;
+
+ p0 = cur->p;
}
- while (cur->next != NULL)
- cur = cur->next;
+ point_assign(Points, n - 1, with_z, &p0, 0);
+ point_list_add(cur, p0); /* always keep last point */
+ }
+ else {
+ for (i = 2; i <= n; i++) {
+ if (i == n)
+ point_assign(Points, 1, with_z, &p2, 0);
+ else
+ point_assign(Points, i, with_z, &p2, 0);
+ point_assign(Points, i - 1, with_z, &p1, 0);
- p0 = cur->p;
+ while (1) {
+ point_add(p1, p2, &tmp);
+ point_scalar(tmp, 0.5, &m1);
+
+ point_list_add(cur, m1);
+
+ if (point_dist_square(p0, m1) > thresh) {
+ point_add(p1, m1, &tmp); /* need to refine the partition */
+ point_scalar(tmp, 0.5, &p2);
+ point_add(p1, p0, &tmp);
+ point_scalar(tmp, 0.5, &p1);
+ }
+ else {
+ break; /* good approximation */
+ }
+ }
+
+ while (cur->next != NULL)
+ cur = cur->next;
+
+ p0 = cur->p;
+ }
}
- point_assign(Points, n - 1, with_z, &p0);
- point_list_add(cur, p0); /* always keep last point */
-
if (point_list_copy_to_line_pnts(head, Points) == -1) {
G_fatal_error(_("Out of memory"));
}
@@ -299,7 +416,7 @@
* returns the number of point in result
*/
int hermite(struct line_pnts *Points, double step, double angle_thresh,
- int with_z)
+ int loop_support, int with_z)
{
POINT_LIST head, *last, *point;
POINT p0, p1, t0, t1, tmp, res;
@@ -308,6 +425,7 @@
double h1, h2, h3, h4;
int n, i;
int ni;
+ int is_loop;
n = Points->n_points;
@@ -316,44 +434,81 @@
return n;
}
+ is_loop = 0;
+
+ /* is it loop ? */
+ if (Points->x[0] == Points->x[n - 1] &&
+ Points->y[0] == Points->y[n - 1] &&
+ (Points->z[0] == Points->z[n - 1] || with_z == 0) && loop_support){
+ is_loop = 1;
+ }
+
/* convert degrees=>radians */
angle_thresh *= M_PI / 180.0;
head.next = NULL;
point = last = &head;
- /* length of p[0]..p[i+1] */
- i = 0;
- point_assign(Points, 0, with_z, &p0);
- point_assign(Points, 1, with_z, &p1);
- /* length of line 0..i */
- length_begin = 0;
- /* length of line from point 0 to i+1 */
- length = point_dist(p0, p1);
- next = 0;
- /* tangent at p0, p1 */
- point_subtract(p1, p0, &t0);
- refine_tangent(&t0);
- point_assign(Points, 2, with_z, &tmp);
- point_subtract(tmp, p0, &t1);
- refine_tangent(&t1);
+ if (!is_loop) {
+ /* length of p[0]..p[i+1] */
+ i = 0;
+ point_assign(Points, 0, with_z, &p0, 0);
+ point_assign(Points, 1, with_z, &p1, 0);
+ /* length of line 0..i */
+ length_begin = 0;
+ /* length of line from point 0 to i+1 */
+ length = point_dist(p0, p1);
+ next = 0;
+ /* tangent at p0, p1 */
+ point_subtract(p1, p0, &t0);
+ refine_tangent(&t0);
+ point_assign(Points, 2, with_z, &tmp, 0);
+ point_subtract(tmp, p0, &t1);
+ refine_tangent(&t1);
+ }
+ else {
+ /* length of p[0]..p[i+1] */
+ i = 0;
+ point_assign(Points, n - 2, with_z, &p0, 0);
+ point_assign(Points, 0, with_z, &p1, 0);
+ /* length of line 0..i */
+ length_begin = 0;
+ /* length of line from point 0 to i+1 */
+ length = point_dist(p0, p1);
+ next = 0;
+ /* tangent at pn-2, p1 */
+ point_assign(Points, 1, with_z, &tmp, 0);
+ point_subtract(tmp, p0, &t0);
+ refine_tangent(&t0);
-
+ point_assign(Points, 0, with_z, &p0, 0);
+ point_assign(Points, 1, with_z, &p1, 0);
+ /* length of line 0..i */
+ length_begin = 0;
+ /* length of line from point 0 to i+1 */
+ length = point_dist(p0, p1);
+ next = 0;
+ /* tangent at p0, p2 */
+ point_assign(Points, 2, with_z, &tmp, 0);
+ point_subtract(tmp, p0, &t1);
+ refine_tangent(&t1);
+ }
+
/* we always operate on the segment point[i]..point[i+1] */
while (i < n - 1) {
if (next > length || (length - length_begin < GRASS_EPSILON)) { /* segmet i..i+1 is finished or too short */
i++;
if (i >= n - 1)
- break; /* we are already out out of line */
- point_assign(Points, i, with_z, &p0);
- point_assign(Points, i + 1, with_z, &p1);
+ break; /* we are already out of line */
+ point_assign(Points, i, with_z, &p0, is_loop);
+ point_assign(Points, i + 1, with_z, &p1, is_loop);
length_begin = length;
length += point_dist(p0, p1);
ni = i + 2;
- if (ni > n - 1)
+ if (!is_loop && ni > n - 1)
ni = n - 1; /* ensure that we are in the line */
t0 = t1;
- point_assign(Points, ni, with_z, &tmp);
+ point_assign(Points, ni, with_z, &tmp, is_loop);
point_subtract(tmp, p0, &t1);
refine_tangent(&t1);
}
@@ -392,14 +547,14 @@
}
}
-
- point_assign(Points, n - 1, with_z, &p0);
+ point_assign(Points, n - 1, with_z, &p0, 0);
point_list_add(last, p0);
if (point_list_copy_to_line_pnts(head, Points) == -1)
G_fatal_error(_("Out of memory"));
point_list_free(head);
+
return Points->n_points;
}
@@ -413,17 +568,33 @@
* instead of O(N^3 * itearations). Probably not needed, for many iterations,
* the result is almost straight line
*/
-int snakes(struct line_pnts *Points, double alpha, double beta, int with_z)
+int snakes(struct line_pnts *Points, double alpha, double beta,
+ int loop_support, int with_z)
{
MATRIX g, ginv, xcoord, ycoord, zcoord, xout, yout, zout;
int n = Points->n_points;
int i, j;
+ double x0, y0, z0;
+ double a = 2.0 * alpha + 6.0 * beta;
+ double b = -alpha - 4.0 * beta;
+ double c = beta;
+ double val[5] = { c, b, a, b, c };
+ int plus = 4;
+ int is_loop = 0;
- if (n < 3)
+ if (n < plus)
return n;
- int plus = 4;
+ /* is it loop ? */
+ if (Points->x[0] == Points->x[n - 1] &&
+ Points->y[0] == Points->y[n - 1] &&
+ (Points->z[0] == Points->z[n - 1] || with_z == 0) && loop_support){
+ is_loop = 1;
+
+ if (n < plus + 2)
+ return n;
+ }
if (!matrix_init(n + 2 * plus, n + 2 * plus, &g)) {
G_fatal_error(_("Out of memory"));
@@ -454,9 +625,9 @@
return n;
}
- double x0 = Points->x[0];
- double y0 = Points->y[0];
- double z0 = Points->z[0];
+ x0 = Points->x[0];
+ y0 = Points->y[0];
+ z0 = Points->z[0];
/* store the coordinates in the column vectors */
for (i = 0; i < n; i++) {
@@ -465,28 +636,43 @@
zcoord.a[i + plus][0] = Points->z[i] - z0;
}
- /* repeat first and last point at the beginning and end
- * of each vector respectively */
- for (i = 0; i < plus; i++) {
- xcoord.a[i][0] = 0;
- ycoord.a[i][0] = 0;
- zcoord.a[i][0] = 0;
+ if (!is_loop) {
+ /* repeat first and last point at the beginning and end
+ * of each vector respectively */
+ for (i = 0; i < plus; i++) {
+ xcoord.a[i][0] = 0;
+ ycoord.a[i][0] = 0;
+ zcoord.a[i][0] = 0;
+ }
+
+ for (i = n + plus; i < n + 2 * plus; i++) {
+ xcoord.a[i][0] = Points->x[n - 1] - x0;
+ ycoord.a[i][0] = Points->y[n - 1] - y0;
+ zcoord.a[i][0] = Points->z[n - 1] - z0;
+ }
}
+ else {
+ /* loop: point 0 and point n - 1 are identical */
- for (i = n + plus; i < n + 2 * plus; i++) {
- xcoord.a[i][0] = Points->x[n - 1] - x0;
- ycoord.a[i][0] = Points->y[n - 1] - y0;
- zcoord.a[i][0] = Points->z[n - 1] - z0;
+ /* repeat last and first points at the beginning and end
+ * of each vector respectively */
+
+ /* add points from n - plus - 1 to n - 2 */
+ for (i = 0, j = n - plus - 1; i < plus; i++, j++) {
+ xcoord.a[i][0] = Points->x[j] - x0;;
+ ycoord.a[i][0] = Points->y[j] - y0;;
+ zcoord.a[i][0] = Points->z[j] - z0;;
+ }
+ /* add points from 1 to plus + 1 */
+ for (i = n + plus, j = 1; i < n + 2 * plus; i++, j++) {
+ xcoord.a[i][0] = Points->x[j] - x0;
+ ycoord.a[i][0] = Points->y[j] - y0;
+ zcoord.a[i][0] = Points->z[j] - z0;
+ }
}
-
/* calculate the matrix */
- double a = 2.0 * alpha + 6.0 * beta;
- double b = -alpha - 4.0 * beta;
- double c = beta;
- double val[5] = { c, b, a, b, c };
-
for (i = 0; i < n + 2 * plus; i++)
for (j = 0; j < n + 2 * plus; j++) {
int index = j - i + 2;
@@ -500,34 +686,51 @@
matrix_add_identity((double)1.0, &g);
/* find its inverse */
- if (!matrix_inverse(g, &ginv, 0)) {
+ if (!matrix_inverse(&g, &ginv, 0)) {
G_fatal_error(_("Unable to find the inverse matrix"));
return n;
}
- if (!matrix_mult(ginv, xcoord, &xout)
- || !matrix_mult(ginv, ycoord, &yout)
- || !matrix_mult(ginv, zcoord, &zout)) {
+ if (!matrix_mult(&ginv, &xcoord, &xout)
+ || !matrix_mult(&ginv, &ycoord, &yout)
+ || !matrix_mult(&ginv, &zcoord, &zout)) {
G_fatal_error(_("Unable to calculate the output vectors"));
return n;
}
- /* copy the new values of coordinates, but
- * never move the last and first point */
- for (i = 1; i < n - 1; i++) {
- Points->x[i] = xout.a[i + plus][0] + x0;
- Points->y[i] = yout.a[i + plus][0] + y0;
- if (with_z)
- Points->z[i] = zout.a[i + plus][0] + z0;
+ if (!is_loop) {
+ /* copy the new values of coordinates, but
+ * never move the last and first point */
+ for (i = 1; i < n - 1; i++) {
+ Points->x[i] = xout.a[i + plus][0] + x0;
+ Points->y[i] = yout.a[i + plus][0] + y0;
+ if (with_z)
+ Points->z[i] = zout.a[i + plus][0] + z0;
+ }
}
+ else {
+ /* copy the new values of coordinates,
+ * also move the last and first point */
+ for (i = 0; i < n; i++) {
+ Points->x[i] = xout.a[i + plus][0] + x0;
+ Points->y[i] = yout.a[i + plus][0] + y0;
+ if (with_z)
+ Points->z[i] = zout.a[i + plus][0] + z0;
+ }
- matrix_free(g);
- matrix_free(ginv);
- matrix_free(xcoord);
- matrix_free(ycoord);
- matrix_free(zcoord);
- matrix_free(xout);
- matrix_free(yout);
- matrix_free(zout);
+ Points->x[n - 1] = Points->x[0];
+ Points->y[n - 1] = Points->y[0];
+ Points->z[n - 1] = Points->z[0];
+ }
+
+ matrix_free(&g);
+ matrix_free(&ginv);
+ matrix_free(&xcoord);
+ matrix_free(&ycoord);
+ matrix_free(&zcoord);
+ matrix_free(&xout);
+ matrix_free(&yout);
+ matrix_free(&zout);
+
return Points->n_points;
}
More information about the grass-commit
mailing list