GRASS Programmer's Manual
6.4.2(2012)
|
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 }