GRASS Programmer's Manual  6.4.2(2012)
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, nnodes;
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     nnodes = Vect_get_num_nodes(Map);
00053     for (node = 1; node <= nnodes; node++) {
00054         int i, nlines;
00055 
00056         G_percent(node, nnodes, 1);
00057         G_debug(3, "node = %d", node);
00058         if (!Vect_node_alive(Map, node))
00059             continue;
00060 
00061         while (1) {
00062             float angle1 = -100;
00063             int line1 = -999;   /* value not important, just for debug */
00064             int clean = 1;
00065 
00066             nlines = Vect_get_node_n_lines(Map, node);
00067             G_debug(3, "nlines = %d", nlines);
00068 
00069             for (i = 0; i < nlines; i++) {
00070                 P_LINE *Line;
00071                 int line2;
00072                 float angle2;
00073 
00074                 line2 = Vect_get_node_line(Map, node, i);
00075                 Line = Map->plus.Line[abs(line2)];
00076                 if (!Line)
00077                     continue;
00078                 G_debug(4, "  type = %d", Line->type);
00079                 if (!(Line->type & (otype & GV_LINES)))
00080                     continue;
00081 
00082                 angle2 = Vect_get_node_line_angle(Map, node, i);
00083                 if (angle2 == -9.0)
00084                     continue;   /* Degenerated line */
00085 
00086                 G_debug(4, "  line1 = %d angle1 = %e line2 = %d angle2 = %e",
00087                         line1, angle1, line2, angle2);
00088 
00089                 if (angle2 == angle1) {
00090                     int j;
00091                     double length1, length2;
00092                     int short_line;     /* line with shorter end segment */
00093                     int long_line;      /* line with longer end segment */
00094                     int new_short_line = 0;     /* line number of short line after rewrite */
00095                     int short_type, long_type, type;
00096                     double x, y, z, nx, ny, nz;
00097 
00098                     G_debug(4, "  identical angles -> clean");
00099 
00100                     /* Length of end segments for both lines */
00101                     Vect_read_line(Map, Points, NULL, abs(line1));
00102                     if (line1 > 0) {
00103                         length1 =
00104                             Vect_points_distance(Points->x[0], Points->y[0],
00105                                                  0.0, Points->x[1],
00106                                                  Points->y[1], 0.0, 0);
00107                     }
00108                     else {
00109                         int np;
00110 
00111                         np = Points->n_points;
00112                         length1 =
00113                             Vect_points_distance(Points->x[np - 1],
00114                                                  Points->y[np - 1], 0.0,
00115                                                  Points->x[np - 2],
00116                                                  Points->y[np - 2], 0.0, 0);
00117                     }
00118 
00119                     Vect_read_line(Map, Points, NULL, abs(line2));
00120                     if (line2 > 0) {
00121                         length2 =
00122                             Vect_points_distance(Points->x[0], Points->y[0],
00123                                                  0.0, Points->x[1],
00124                                                  Points->y[1], 0.0, 0);
00125                     }
00126                     else {
00127                         int np;
00128 
00129                         np = Points->n_points;
00130                         length2 =
00131                             Vect_points_distance(Points->x[np - 1],
00132                                                  Points->y[np - 1], 0.0,
00133                                                  Points->x[np - 2],
00134                                                  Points->y[np - 2], 0.0, 0);
00135                     }
00136 
00137                     G_debug(4, "  length1 = %f length2 = %f", length1,
00138                             length2);
00139 
00140                     if (length1 < length2) {
00141                         short_line = line1;
00142                         long_line = line2;
00143                     }
00144                     else {
00145                         short_line = line2;
00146                         long_line = line1;
00147                     }
00148 
00149                     /* Remove end segment from short_line */
00150                     short_type =
00151                         Vect_read_line(Map, Points, SCats, abs(short_line));
00152 
00153                     if (short_line > 0) {
00154                         x = Points->x[1];
00155                         y = Points->y[1];
00156                         z = Points->z[1];
00157                         Vect_line_delete_point(Points, 0);      /* first */
00158                     }
00159                     else {
00160                         x = Points->x[Points->n_points - 2];
00161                         y = Points->y[Points->n_points - 2];
00162                         z = Points->z[Points->n_points - 2];
00163                         Vect_line_delete_point(Points, Points->n_points - 1);   /* last */
00164                     }
00165 
00166                     /* It may happen that it is one line: node could be deleted,
00167                      * in that case we have to read the node coords first */
00168                     Vect_get_node_coor(Map, node, &nx, &ny, &nz);
00169 
00170                     if (Points->n_points > 1) {
00171                         new_short_line =
00172                             Vect_rewrite_line(Map, abs(short_line),
00173                                               short_type, Points, SCats);
00174                     }
00175                     else {
00176                         Vect_delete_line(Map, abs(short_line));
00177                     }
00178 
00179                     /* It may happen that it is one line, in that case we have to take the new
00180                      * short line as long line, orientation is not changed */
00181                     if (abs(line1) == abs(line2)) {
00182                         if (long_line > 0)
00183                             long_line = new_short_line;
00184                         else
00185                             long_line = -new_short_line;
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         nnodes = Vect_get_num_nodes(Map);
00252     }
00253 
00254     return (nmodif);
00255 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines