v.generalize: add loop support for boyle, chaiken, distance_weighting, hermite, snakes

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@56264 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/v.generalize/displacement.c

@@ -90,7 +90,7 @@ int snakes_displacement(struct Map_info *In, struct Map_info *Out,
 	    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 @@ int snakes_displacement(struct Map_info *In, struct Map_info *Out,
 
 	    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 @@ int snakes_displacement(struct Map_info *In, struct Map_info *Out,
 
     /*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 @@ int snakes_displacement(struct Map_info *In, struct Map_info *Out,
 	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 @@ int snakes_displacement(struct Map_info *In, struct Map_info *Out,
 	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 @@ int snakes_displacement(struct Map_info *In, struct Map_info *Out,
 		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 @@ int snakes_displacement(struct Map_info *In, struct Map_info *Out,
     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;
 }

vector/v.generalize/main.c

@@ -228,7 +228,8 @@ int main(int argc, char *argv[])
 
     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 @@ int main(int argc, char *argv[])
 		    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);
 

vector/v.generalize/matrix.c

@@ -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 @@ int matrix_init(int rows, int cols, MATRIX * res)
     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 matrix_mult(MATRIX a, MATRIX b, MATRIX * res)
 
     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 @@ int matrix_add_identity(double s, MATRIX * m)
 /* 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 @@ void matrix_swap_rows(int x, int y, MATRIX * m)
 
 /* 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 @@ void matrix_row_scalar(int row, double s, MATRIX * m)
  * 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 @@ void matrix_row_add_multiple(int ra, int rb, double s, MATRIX * m)
 }
 
 /* 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 @@ int matrix_inverse(MATRIX a, MATRIX * res, int percents)
      * is row equivalent to the identity matrix.
      */
 
-    int n = a.rows;
+    int n = a->rows;
 
     if (percents)
 	G_percent_reset();
@@ -170,27 +168,27 @@ int matrix_inverse(MATRIX a, MATRIX * res, int percents)
 	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 @@ int matrix_inverse(MATRIX a, MATRIX * res, int percents)
     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 @@ void matrix_mult_scalar(double s, MATRIX * m)
 	    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);
     }

vector/v.generalize/matrix.h

@@ -29,32 +29,32 @@ typedef struct
 } 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

vector/v.generalize/operators.h

@@ -10,15 +10,15 @@ int douglas_peucker_reduction(struct line_pnts *Points, double thresh,
 			      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, 

vector/v.generalize/point.c

@@ -50,23 +50,12 @@ inline double point_dist2(POINT a)
 }
 
 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) {

vector/v.generalize/point.h

@@ -49,13 +49,12 @@ extern inline double point_dist2(POINT a);
  * 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);
 

vector/v.generalize/simplification.c

@@ -239,15 +239,15 @@ int reumann_witkam(struct line_pnts *Points, double thresh, int with_z)
 
     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 @@ int reumann_witkam(struct line_pnts *Points, double thresh, int with_z)
 	 * 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);
 

vector/v.generalize/smoothing.c

@@ -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);
-    return Points->n_points;
+    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 @@ int sliding_averaging(struct line_pnts *Points, double slide, int look_ahead,
     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 @@ int sliding_averaging(struct line_pnts *Points, double slide, int look_ahead,
 	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,39 +148,45 @@ int sliding_averaging(struct line_pnts *Points, double slide, int look_ahead,
 
     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 @@ int sliding_averaging(struct line_pnts *Points, double slide, int look_ahead,
 /* 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]);
-
-    half = look_ahead / 2;
+    point_assign(Points, 0, with_z, &res[0], 0);
 
-    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 @@ int distance_weighting(struct line_pnts *Points, double slide, int look_ahead,
 	    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,46 +306,88 @@ int chaiken(struct line_pnts *Points, double thresh, int with_z)
     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;
-
-    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);
-
-	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 */
+    point_assign(Points, 0, with_z, &p0, 0);
+
+    if (!is_loop) {
+	point_list_add(cur, p0); /* always keep first point */
+	cur = cur->next;
+
+	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);
+
+	    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;
+	    while (cur->next != NULL)
+		cur = cur->next;
 
-	p0 = cur->p;
+	    p0 = cur->p;
+	}
+
+	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);
+
+	    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 */
+		}
+	    }
 
-    point_assign(Points, n - 1, with_z, &p0);
-    point_list_add(cur, p0); /* always keep last point */
+	    while (cur->next != NULL)
+		cur = cur->next;
+
+	    p0 = cur->p;
+	}
+    }
 
     if (point_list_copy_to_line_pnts(head, Points) == -1) {
 	G_fatal_error(_("Out of memory"));
@@ -299,7 +416,7 @@ void refine_tangent(POINT * p)
  * 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 @@ int hermite(struct line_pnts *Points, double step, double angle_thresh,
     double h1, h2, h3, h4;
     int n, i;
     int ni;
+    int is_loop;
 
     n = Points->n_points;
 
@@ -316,44 +434,81 @@ int hermite(struct line_pnts *Points, double step, double angle_thresh,
 	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 @@ int hermite(struct line_pnts *Points, double step, double angle_thresh,
 	}
     }
 
-
-    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 @@ int hermite(struct line_pnts *Points, double step, double angle_thresh,
  * 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 @@ int snakes(struct line_pnts *Points, double alpha, double beta, int with_z)
 	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,27 +636,42 @@ int snakes(struct line_pnts *Points, double alpha, double beta, int with_z)
 	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;
+	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 */
 
+	/* repeat last and first points at the beginning and end
+	 * of each vector respectively */
 
-    /* calculate the matrix */
-    double a = 2.0 * alpha + 6.0 * beta;
-    double b = -alpha - 4.0 * beta;
-    double c = beta;
+	/* 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;
+	}
+    }
 
-    double val[5] = { c, b, a, b, c };
+    /* calculate the matrix */
 
     for (i = 0; i < n + 2 * plus; i++)
 	for (j = 0; j < n + 2 * plus; j++) {
@@ -500,34 +686,51 @@ int snakes(struct line_pnts *Points, double alpha, double beta, int with_z)
     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;
-    }
-
-    matrix_free(g);
-    matrix_free(ginv);
-    matrix_free(xcoord);
-    matrix_free(ycoord);
-    matrix_free(zcoord);
-    matrix_free(xout);
-    matrix_free(yout);
-    matrix_free(zout);
+    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;
+	}
+
+	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;
 }