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;
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;
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;
00092 int long_line;
00093 int new_short_line = 0;
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
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
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);
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);
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
00175
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
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
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 }