clean_nodes.c

Go to the documentation of this file.
00001 
00020 #include <stdlib.h>
00021 #include <grass/gis.h>
00022 #include <grass/Vect.h>
00023 #include <grass/glocale.h>
00024 
00038 int
00039 Vect_clean_small_angles_at_nodes(struct Map_info *Map, int otype,
00040                                  struct Map_info *Err)
00041 {
00042     int node;
00043     int nmodif = 0;
00044     struct line_pnts *Points;
00045     struct line_cats *SCats, *LCats, *OCats;
00046 
00047     Points = Vect_new_line_struct();
00048     SCats = Vect_new_cats_struct();
00049     LCats = Vect_new_cats_struct();
00050     OCats = Vect_new_cats_struct();
00051 
00052     for (node = 1; node <= Vect_get_num_nodes(Map); node++) {
00053         int i, nlines;
00054 
00055         G_debug(3, "node = %d", node);
00056         if (!Vect_node_alive(Map, node))
00057             continue;
00058 
00059         while (1) {
00060             float angle1 = -100;
00061             int line1 = -999;   /* value not important, just for debug */
00062             int clean = 1;
00063 
00064 
00065             nlines = Vect_get_node_n_lines(Map, node);
00066             G_debug(3, "nlines = %d", nlines);
00067 
00068             for (i = 0; i < nlines; i++) {
00069                 P_LINE *Line;
00070                 int line2;
00071                 float angle2;
00072 
00073                 line2 = Vect_get_node_line(Map, node, i);
00074                 Line = Map->plus.Line[abs(line2)];
00075                 if (!Line)
00076                     continue;
00077                 G_debug(4, "  type = %d", Line->type);
00078                 if (!(Line->type & (otype & GV_LINES)))
00079                     continue;
00080 
00081                 angle2 = Vect_get_node_line_angle(Map, node, i);
00082                 if (angle2 == -9.0)
00083                     continue;   /* Degenerated line */
00084 
00085                 G_debug(4, "  line1 = %d angle1 = %e line2 = %d angle2 = %e",
00086                         line1, angle1, line2, angle2);
00087 
00088                 if (angle2 == angle1) {
00089                     int j;
00090                     double length1, length2;
00091                     int short_line;     /* line with shorter end segment */
00092                     int long_line;      /* line with longer end segment */
00093                     int new_short_line = 0;     /* line number of short line after rewrite */
00094                     int short_type, long_type, type;
00095                     double x, y, z, nx, ny, nz;
00096 
00097                     G_debug(4, "  identical angles -> clean");
00098 
00099                     /* Length of end segments for both lines */
00100                     Vect_read_line(Map, Points, NULL, abs(line1));
00101                     if (line1 > 0) {
00102                         length1 =
00103                             Vect_points_distance(Points->x[0], Points->y[0],
00104                                                  0.0, Points->x[1],
00105                                                  Points->y[1], 0.0, 0);
00106                     }
00107                     else {
00108                         int np;
00109 
00110                         np = Points->n_points;
00111                         length1 =
00112                             Vect_points_distance(Points->x[np - 1],
00113                                                  Points->y[np - 1], 0.0,
00114                                                  Points->x[np - 2],
00115                                                  Points->y[np - 2], 0.0, 0);
00116                     }
00117 
00118                     Vect_read_line(Map, Points, NULL, abs(line2));
00119                     if (line2 > 0) {
00120                         length2 =
00121                             Vect_points_distance(Points->x[0], Points->y[0],
00122                                                  0.0, Points->x[1],
00123                                                  Points->y[1], 0.0, 0);
00124                     }
00125                     else {
00126                         int np;
00127 
00128                         np = Points->n_points;
00129                         length2 =
00130                             Vect_points_distance(Points->x[np - 1],
00131                                                  Points->y[np - 1], 0.0,
00132                                                  Points->x[np - 2],
00133                                                  Points->y[np - 2], 0.0, 0);
00134                     }
00135 
00136                     G_debug(4, "  length1 = %f length2 = %f", length1,
00137                             length2);
00138 
00139                     if (length1 < length2) {
00140                         short_line = line1;
00141                         long_line = line2;
00142                     }
00143                     else {
00144                         short_line = line2;
00145                         long_line = line1;
00146                     }
00147 
00148                     /* Remove end segment from short_line */
00149                     short_type =
00150                         Vect_read_line(Map, Points, SCats, abs(short_line));
00151 
00152                     if (short_line > 0) {
00153                         x = Points->x[1];
00154                         y = Points->y[1];
00155                         z = Points->z[1];
00156                         Vect_line_delete_point(Points, 0);      /* first */
00157                     }
00158                     else {
00159                         x = Points->x[Points->n_points - 2];
00160                         y = Points->y[Points->n_points - 2];
00161                         z = Points->z[Points->n_points - 2];
00162                         Vect_line_delete_point(Points, Points->n_points - 1);   /* last */
00163                     }
00164 
00165             /* It may happen that it is one line: node could be deleted,
00166              * in that case we have to read the node coords first */
00167             Vect_get_node_coor(Map, node, &nx, &ny, &nz);
00168 
00169                     if (Points->n_points > 1) {
00170                         new_short_line =
00171                             Vect_rewrite_line(Map, abs(short_line),
00172                                               short_type, Points, SCats);
00173                     }
00174                     else {
00175                         Vect_delete_line(Map, abs(short_line));
00176                     }
00177 
00178                     /* It may happen that it is one line, in that case we have to take the new
00179                      * short line as long line, orientation is not changed */
00180                     if (abs(line1) == abs(line2)) {
00181                         if (long_line > 0)
00182                             long_line = new_short_line;
00183                         else
00184                             long_line = -new_short_line;
00185                     }
00186 
00187 
00188                     /* Add new line (must be before rewrite of long_line otherwise node could be deleted) */
00189                     long_type =
00190                         Vect_read_line(Map, NULL, LCats, abs(long_line));
00191 
00192                     Vect_reset_cats(OCats);
00193                     for (j = 0; j < SCats->n_cats; j++) {
00194                         Vect_cat_set(OCats, SCats->field[j], SCats->cat[j]);
00195                     }
00196                     for (j = 0; j < LCats->n_cats; j++) {
00197                         Vect_cat_set(OCats, LCats->field[j], LCats->cat[j]);
00198                     }
00199 
00200                     if (long_type == GV_BOUNDARY || short_type == GV_BOUNDARY) {
00201                         type = GV_BOUNDARY;
00202                     }
00203                     else {
00204                         type = GV_LINE;
00205                     }
00206 
00207                     Vect_reset_line(Points);
00208                     Vect_append_point(Points, nx, ny, nz);
00209                     Vect_append_point(Points, x, y, z);
00210                     Vect_write_line(Map, type, Points, OCats);
00211 
00212                     if (Err) {
00213                         Vect_write_line(Err, type, Points, OCats);
00214                     }
00215 
00216                     /* Snap long_line to the new short_line end */
00217                     long_type =
00218                         Vect_read_line(Map, Points, LCats, abs(long_line));
00219                     if (long_line > 0) {
00220                         Points->x[0] = x;
00221                         Points->y[0] = y;
00222                         Points->z[0] = z;
00223                     }
00224                     else {
00225                         Points->x[Points->n_points - 1] = x;
00226                         Points->y[Points->n_points - 1] = y;
00227                         Points->z[Points->n_points - 1] = z;
00228                     }
00229                     Vect_line_prune(Points);
00230                     if (Points->n_points > 1) {
00231                         Vect_rewrite_line(Map, abs(long_line), long_type,
00232                                           Points, LCats);
00233                     }
00234                     else {
00235                         Vect_delete_line(Map, abs(long_line));
00236                     }
00237 
00238                     nmodif += 3;
00239                     clean = 0;
00240 
00241                     break;
00242                 }
00243 
00244                 line1 = line2;
00245                 angle1 = angle2;
00246             }
00247 
00248         if (clean || !Vect_node_alive(Map, node))
00249                 break;
00250         }
00251     }
00252 
00253     return (nmodif);
00254 }