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                     if (Points->n_points > 1) {
00166                         new_short_line =
00167                             Vect_rewrite_line(Map, abs(short_line),
00168                                               short_type, Points, SCats);
00169                     }
00170                     else {
00171                         Vect_delete_line(Map, abs(short_line));
00172                     }
00173 
00174                     /* It may happen that it is one line, in that case we have to take the new
00175                      * short line as long line, orientation is not changed */
00176                     if (abs(line1) == abs(line2)) {
00177                         if (long_line > 0)
00178                             long_line = new_short_line;
00179                         else
00180                             long_line = -new_short_line;
00181                     }
00182 
00183 
00184                     /* Add new line (must be before rewrite of long_line otherwise node could be deleted) */
00185                     long_type =
00186                         Vect_read_line(Map, NULL, LCats, abs(long_line));
00187 
00188                     Vect_reset_cats(OCats);
00189                     for (j = 0; j < SCats->n_cats; j++) {
00190                         Vect_cat_set(OCats, SCats->field[j], SCats->cat[j]);
00191                     }
00192                     for (j = 0; j < LCats->n_cats; j++) {
00193                         Vect_cat_set(OCats, LCats->field[j], LCats->cat[j]);
00194                     }
00195 
00196                     if (long_type == GV_BOUNDARY || short_type == GV_BOUNDARY) {
00197                         type = GV_BOUNDARY;
00198                     }
00199                     else {
00200                         type = GV_LINE;
00201                     }
00202 
00203                     Vect_get_node_coor(Map, node, &nx, &ny, &nz);
00204                     Vect_reset_line(Points);
00205                     Vect_append_point(Points, nx, ny, nz);
00206                     Vect_append_point(Points, x, y, z);
00207                     Vect_write_line(Map, type, Points, OCats);
00208 
00209                     if (Err) {
00210                         Vect_write_line(Err, type, Points, OCats);
00211                     }
00212 
00213                     /* Snap long_line to the new short_line end */
00214                     long_type =
00215                         Vect_read_line(Map, Points, LCats, abs(long_line));
00216                     if (long_line > 0) {
00217                         Points->x[0] = x;
00218                         Points->y[0] = y;
00219                         Points->z[0] = z;
00220                     }
00221                     else {
00222                         Points->x[Points->n_points - 1] = x;
00223                         Points->y[Points->n_points - 1] = y;
00224                         Points->z[Points->n_points - 1] = z;
00225                     }
00226                     Vect_line_prune(Points);
00227                     if (Points->n_points > 1) {
00228                         Vect_rewrite_line(Map, abs(long_line), long_type,
00229                                           Points, LCats);
00230                     }
00231                     else {
00232                         Vect_delete_line(Map, abs(long_line));
00233                     }
00234 
00235                     nmodif += 3;
00236                     clean = 0;
00237 
00238                     break;
00239                 }
00240 
00241                 line1 = line2;
00242                 angle1 = angle2;
00243             }
00244 
00245             if (clean)
00246                 break;
00247         }
00248     }
00249 
00250     return (nmodif);
00251 }
Generated on Tue Apr 6 13:28:10 2010 for GRASS Programmer's Manual by  doxygen 1.6.3