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