File indexing completed on 2024-05-12 15:27:08
0001 /*************************************************************************** 0002 File : nsl_geom_linesim.c 0003 Project : LabPlot 0004 Description : NSL geometry line simplification functions 0005 -------------------------------------------------------------------- 0006 Copyright : (C) 2016-2019 by Stefan Gerlach (stefan.gerlach@uni.kn) 0007 0008 ***************************************************************************/ 0009 0010 /*************************************************************************** 0011 * * 0012 * This program is free software; you can redistribute it and/or modify * 0013 * it under the terms of the GNU General Public License as published by * 0014 * the Free Software Foundation; either version 2 of the License, or * 0015 * (at your option) any later version. * 0016 * * 0017 * This program is distributed in the hope that it will be useful, * 0018 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 0019 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 0020 * GNU General Public License for more details. * 0021 * * 0022 * You should have received a copy of the GNU General Public License * 0023 * along with this program; if not, write to the Free Software * 0024 * Foundation, Inc., 51 Franklin Street, Fifth Floor, * 0025 * Boston, MA 02110-1301 USA * 0026 * * 0027 ***************************************************************************/ 0028 0029 #include "nsl_geom_linesim.h" 0030 #include "nsl_geom.h" 0031 #include "nsl_common.h" 0032 #include "nsl_sort.h" 0033 #include "nsl_stats.h" 0034 0035 const char* nsl_geom_linesim_type_name[] = {i18n("Douglas-Peucker (number)"), i18n("Douglas-Peucker (tolerance)"), i18n("Visvalingam-Whyatt"), i18n("Reumann-Witkam"), i18n("perpendicular distance"), i18n("n-th point"), 0036 i18n("radial distance"), i18n("Interpolation"), i18n("Opheim"), i18n("Lang")}; 0037 0038 /*********** error calculation functions *********/ 0039 0040 double nsl_geom_linesim_positional_error(const double xdata[], const double ydata[], const size_t n, const size_t index[]) { 0041 double dist = 0; 0042 size_t i = 0, j; /* i: index of index[] */ 0043 do { 0044 /*for every point not in index[] calculate distance to line*/ 0045 /*printf("i=%d (index[i]-index[i+1]=%d-%d)\n", i , index[i], index[i+1]);*/ 0046 for (j = 1; j < index[i+1]-index[i]; j++) { 0047 /*printf("i=%d: j=%d\n", i, j);*/ 0048 dist += nsl_geom_point_line_dist(xdata[index[i]], ydata[index[i]], xdata[index[i+1]], ydata[index[i+1]], xdata[index[i]+j], ydata[index[i]+j]); 0049 /*printf("dist = %g\n", dist);*/ 0050 } 0051 i++; 0052 } while(index[i] != n-1); 0053 0054 return dist/(double)n; 0055 } 0056 0057 double nsl_geom_linesim_positional_squared_error(const double xdata[], const double ydata[], const size_t n, const size_t index[]) { 0058 double dist = 0; 0059 size_t i = 0, j; /* i: index of index[] */ 0060 do { 0061 /*for every point not in index[] calculate distance to line*/ 0062 /*printf("i=%d (index[i]-index[i+1]=%d-%d)\n", i , index[i], index[i+1]);*/ 0063 for (j = 1; j < index[i+1]-index[i]; j++) { 0064 /*printf("i=%d: j=%d\n", i, j);*/ 0065 dist += pow(nsl_geom_point_line_dist(xdata[index[i]], ydata[index[i]], xdata[index[i+1]], ydata[index[i+1]], xdata[index[i]+j], ydata[index[i]+j]), 2); 0066 /*printf("dist = %g\n", dist);*/ 0067 } 0068 i++; 0069 } while(index[i] != n-1); 0070 0071 return dist/(double)n; 0072 } 0073 0074 double nsl_geom_linesim_area_error(const double xdata[], const double ydata[], const size_t n, const size_t index[]) { 0075 double area = 0; 0076 0077 size_t i = 0, j; /* i: index of index[] */ 0078 do { 0079 /* for every point not in index[] calculate area */ 0080 /*printf("i=%d (index[i]-index[i+1]=%d-%d)\n", i , index[i], index[i+1]);*/ 0081 for (j = 1; j < index[i+1]-index[i]; j++) { 0082 /*printf("j=%d: area between %d %d %d\n", j, index[i]+j-1, index[i]+j, index[i+1]);*/ 0083 area += nsl_geom_three_point_area(xdata[index[i]+j-1], ydata[index[i]+j-1], xdata[index[i]+j], ydata[index[i]+j], xdata[index[i+1]], ydata[index[i+1]]); 0084 /*printf("area = %g\n", area);*/ 0085 } 0086 i++; 0087 } while(index[i] != n-1); 0088 0089 0090 return area/(double)n; 0091 } 0092 0093 double nsl_geom_linesim_clip_diag_perpoint(const double xdata[], const double ydata[], const size_t n) { 0094 double dx = nsl_stats_maximum(xdata, n, NULL) - nsl_stats_minimum(xdata, n, NULL); 0095 double dy = nsl_stats_maximum(ydata, n, NULL) - nsl_stats_minimum(ydata, n, NULL); 0096 double d = sqrt(dx*dx+dy*dy); 0097 0098 return d/(double)n; /* per point */ 0099 } 0100 0101 double nsl_geom_linesim_clip_area_perpoint(const double xdata[], const double ydata[], const size_t n) { 0102 double dx = nsl_stats_maximum(xdata, n, NULL) - nsl_stats_minimum(xdata, n, NULL); 0103 double dy = nsl_stats_maximum(ydata, n, NULL) - nsl_stats_minimum(ydata, n, NULL); 0104 double A = dx*dy; 0105 0106 return A/(double)n; /* per point */ 0107 } 0108 0109 double nsl_geom_linesim_avg_dist_perpoint(const double xdata[], const double ydata[], const size_t n) { 0110 double dist = 0; 0111 size_t i; 0112 for (i = 0; i < n-1; i++) { 0113 double dx = xdata[i+1] - xdata[i]; 0114 double dy = ydata[i+1] - ydata[i]; 0115 dist += sqrt(dx*dx + dy*dy); 0116 } 0117 dist /= (double)n; 0118 0119 return dist; 0120 } 0121 0122 /*********** simplification algorithms *********/ 0123 0124 void nsl_geom_linesim_douglas_peucker_step(const double xdata[], const double ydata[], const size_t start, const size_t end, size_t *nout, const double tol, size_t index[]) { 0125 /*printf("DP: %d - %d\n", start, end);*/ 0126 0127 size_t i, nkey = start; 0128 double maxdist = 0; 0129 /* search for key (biggest perp. distance) */ 0130 for (i = start+1; i < end; i++) { 0131 double dist = nsl_geom_point_line_dist(xdata[start], ydata[start], xdata[end], ydata[end], xdata[i], ydata[i]); 0132 if (dist > maxdist) { 0133 maxdist = dist; 0134 nkey = i; 0135 } 0136 } 0137 /*printf("maxdist = %g @ i = %zu\n", maxdist, nkey);*/ 0138 0139 if (maxdist > tol) { 0140 /*printf("take %d\n", nkey);*/ 0141 index[(*nout)++] = nkey; 0142 if (nkey-start > 1) 0143 nsl_geom_linesim_douglas_peucker_step(xdata, ydata, start, nkey, nout, tol, index); 0144 if (end-nkey > 1) 0145 nsl_geom_linesim_douglas_peucker_step(xdata, ydata, nkey, end, nout, tol, index); 0146 } 0147 0148 /*printf("nout=%d\n", *nout);*/ 0149 } 0150 0151 size_t nsl_geom_linesim_douglas_peucker(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { 0152 size_t nout = 0; 0153 0154 /*first point*/ 0155 index[nout++] = 0; 0156 0157 nsl_geom_linesim_douglas_peucker_step(xdata, ydata, 0, n-1, &nout, tol, index); 0158 0159 /* last point */ 0160 if (index[nout-1] != n-1) 0161 index[nout++] = n-1; 0162 0163 /* sort array index */ 0164 nsl_sort_size_t(index, nout); 0165 0166 return nout; 0167 } 0168 size_t nsl_geom_linesim_douglas_peucker_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { 0169 double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); 0170 return nsl_geom_linesim_douglas_peucker(xdata, ydata, n, tol, index); 0171 } 0172 0173 /* 0174 * Douglas-Peucker variant: 0175 * The key of all egdes of the current simplified line is calculated and only the 0176 * largest is added. This is repeated until nout is reached. 0177 * */ 0178 double nsl_geom_linesim_douglas_peucker_variant(const double xdata[], const double ydata[], const size_t n, const size_t nout, size_t index[]) { 0179 size_t i; 0180 if (nout >= n) { /* use all points */ 0181 for (i = 0; i < n; i++) 0182 index[i] = i; 0183 return 0; 0184 } 0185 0186 /* set first and last point in index (other indizes not initialized) */ 0187 size_t ncount = 0; 0188 index[ncount++] = 0; 0189 index[ncount++] = n-1; 0190 0191 if (nout <= 2) /* use only first and last point (perp. dist is zero) */ 0192 return 0.0; 0193 0194 double *dist = (double *)malloc(n * sizeof(double)); 0195 if (dist == NULL) { 0196 /* printf("nsl_geom_linesim_douglas_peucker_variant(): ERROR allocating memory for 'dist'!\n"); */ 0197 return DBL_MAX; 0198 } 0199 0200 double *maxdist = (double *)malloc(nout * sizeof(double)); /* max dist per edge */ 0201 if (maxdist == NULL) { 0202 /* printf("nsl_geom_linesim_douglas_peucker_variant(): ERROR allocating memory for 'maxdist'!\n"); */ 0203 free(dist); 0204 return DBL_MAX; 0205 } 0206 0207 for (i = 0; i < n; i++) { /* initialize dist */ 0208 dist[i] = nsl_geom_point_line_dist(xdata[0], ydata[0], xdata[n-1], ydata[n-1], xdata[i], ydata[i]); 0209 /* printf("nsl_geom_linesim_douglas_peucker_variant(): %zu dist = %g\n", i, dist[i]); */ 0210 } 0211 for (i = 0; i < nout; i++) 0212 maxdist[i] = 0; 0213 0214 double newmaxdist = 0; 0215 do { 0216 /* printf("nsl_geom_linesim_douglas_peucker_variant(): NEXT ITERATION\n"); */ 0217 size_t key = 0; 0218 0219 /* find edge of maximum */ 0220 size_t maxindex; 0221 nsl_stats_maximum(maxdist, ncount, &maxindex); 0222 /* printf("nsl_geom_linesim_douglas_peucker_variant(): found edge with max dist at index %zu\n", maxindex); */ 0223 /*newmaxdist = nsl_stats_maximum(dist, n, &key);*/ 0224 newmaxdist = 0; 0225 for (i = index[maxindex]+1; i < index[maxindex+1]; i++) { 0226 /* printf("nsl_geom_linesim_douglas_peucker_variant(): iterate i=%zu\n", i); */ 0227 if (dist[i] > newmaxdist) { 0228 newmaxdist = dist[i]; 0229 key = i; 0230 } 0231 } 0232 0233 /* printf("nsl_geom_linesim_douglas_peucker_variant(): found key %zu (dist = %g)\n", key, newmaxdist); */ 0234 ncount++; 0235 dist[key] = 0; 0236 0237 /* find index of previous key */ 0238 size_t previndex = 0; 0239 while (index[previndex+1] < key) 0240 previndex++; 0241 /* printf("nsl_geom_linesim_douglas_peucker_variant(): previndex = %zu (update keys %zu - %zu)\n", previndex, index[previndex], index[previndex+1]); */ 0242 0243 size_t v; 0244 /* printf("nsl_geom_linesim_douglas_peucker_variant(): ncount = %zu, previndex = %zu\n", ncount, previndex); */ 0245 /* update dist[]. no update on last key */ 0246 if (ncount < nout) { 0247 /* shift maxdist */ 0248 for (v = ncount; v > previndex; v--) 0249 maxdist[v] = maxdist[v-1]; 0250 0251 double tmpmax = 0; 0252 for (v = index[previndex]+1; v < key; v++) { 0253 /* printf("nsl_geom_linesim_douglas_peucker_variant(): %zu in %zu - %zu", v, index[previndex], key); */ 0254 dist[v] = nsl_geom_point_line_dist(xdata[index[previndex]], ydata[index[previndex]], xdata[key], ydata[key], 0255 xdata[v], ydata[v]); 0256 if (dist[v] > tmpmax) 0257 tmpmax = dist[v]; 0258 0259 /* printf(" dist = %g\n", dist[v]); */ 0260 } 0261 maxdist[previndex] = tmpmax; 0262 0263 tmpmax = 0; 0264 for (v = key+1; v < index[previndex+1]; v++) { 0265 /* printf("nsl_geom_linesim_douglas_peucker_variant(): %zu in %zu - %zu", v, key, index[previndex+1]); */ 0266 dist[v] = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[index[previndex+1]], ydata[index[previndex+1]], 0267 xdata[v], ydata[v]); 0268 if (dist[v] > tmpmax) 0269 tmpmax = dist[v]; 0270 /* printf(" dist = %g\n", dist[v]); */ 0271 } 0272 maxdist[previndex+1] = tmpmax; 0273 } 0274 0275 /* put into index array */ 0276 for (v = ncount; v > previndex+1; v--) 0277 index[v] = index[v-1]; 0278 index[previndex+1] = key; 0279 } while (ncount < nout); 0280 0281 free(maxdist); 0282 free(dist); 0283 0284 return newmaxdist; 0285 } 0286 0287 size_t nsl_geom_linesim_nthpoint(const size_t n, const int step, size_t index[]) { 0288 if (step < 1) { 0289 printf("step size must be > 0 (given: %d)\n", step); 0290 return 0; 0291 } 0292 0293 size_t i, nout = 0; 0294 0295 /*first point*/ 0296 index[nout++] = 0; 0297 0298 for (i = 1; i < n-1; i++) 0299 if (i%step == 0) 0300 index[nout++] = i; 0301 0302 /* last point */ 0303 index[nout++] = n-1; 0304 0305 return nout; 0306 } 0307 0308 size_t nsl_geom_linesim_raddist(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { 0309 size_t i, nout = 0, key = 0; 0310 0311 /*first point*/ 0312 index[nout++] = 0; 0313 0314 for (i = 1; i < n-1; i++) { 0315 /* distance to key point */ 0316 double dist = nsl_geom_point_point_dist(xdata[i], ydata[i], xdata[key], ydata[key]); 0317 /* distance to last point */ 0318 double lastdist = nsl_geom_point_point_dist(xdata[i], ydata[i], xdata[n-1], ydata[n-1]); 0319 /*printf("%d: %g %g\n", i, dist, lastdist);*/ 0320 0321 if (dist > tol && lastdist > tol) { 0322 index[nout++] = i; 0323 key = i; 0324 } 0325 } 0326 0327 /* last point */ 0328 index[nout++] = n-1; 0329 0330 return nout; 0331 } 0332 size_t nsl_geom_linesim_raddist_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { 0333 double tol = 10.*nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); 0334 return nsl_geom_linesim_raddist(xdata, ydata, n, tol, index); 0335 } 0336 0337 size_t nsl_geom_linesim_perpdist(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { 0338 size_t nout = 0, key = 0, i; 0339 0340 /*first point*/ 0341 index[nout++] = 0; 0342 0343 for (i = 1; i < n-1; i++) { 0344 /* distance of point i to line key -- i+1 */ 0345 double dist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[i+1], ydata[i+1], xdata[i], ydata[i]); 0346 /*printf("%d: %g\n", i, dist);*/ 0347 0348 if (dist > tol) { /* take it */ 0349 index[nout++] = i; 0350 key = i; 0351 /*printf("%d: take it (key = %d)\n", i, key);*/ 0352 } else { /* ignore it */ 0353 if (i+1 < n-1) /* last point is taken anyway */ 0354 index[nout++] = i+1; /* take next point in any case */ 0355 /*printf("%d: ignore it (key = %d)\n", i, i+1);*/ 0356 key = ++i; 0357 } 0358 } 0359 0360 /* last point */ 0361 index[nout++] = n-1; 0362 0363 return nout; 0364 } 0365 size_t nsl_geom_linesim_perpdist_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { 0366 double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); 0367 return nsl_geom_linesim_perpdist(xdata, ydata, n, tol, index); 0368 } 0369 0370 size_t nsl_geom_linesim_perpdist_repeat(const double xdata[], const double ydata[], const size_t n, const double tol, const size_t repeat, size_t index[]) { 0371 size_t i, j, nout; 0372 double *xtmp = (double *) malloc(n*sizeof(double)); 0373 if (xtmp == NULL) { 0374 printf("nsl_geom_linesim_perpdist_repeat(): ERROR allocating memory for 'xtmp'!\n"); 0375 return 0; 0376 } 0377 0378 double *ytmp = (double *) malloc(n*sizeof(double)); 0379 if (ytmp == NULL) { 0380 printf("nsl_geom_linesim_perpdist_repeat(): ERROR allocating memory for 'ytmp'!\n"); 0381 free(xtmp); 0382 return 0; 0383 } 0384 0385 size_t *tmpindex = (size_t *) malloc(n*sizeof(size_t)); 0386 if (tmpindex == NULL) { 0387 printf("nsl_geom_linesim_perpdist_repeat(): ERROR allocating memory for 'tmpindex'!\n"); 0388 free(xtmp); 0389 free(ytmp); 0390 return 0; 0391 } 0392 0393 /* first round */ 0394 nout = nsl_geom_linesim_perpdist(xdata, ydata, n, tol, index); 0395 /* repeats */ 0396 for (i = 0; i < repeat - 1; i++) { 0397 for (j = 0; j < nout; j++) { 0398 xtmp[j] = xdata[index[j]]; 0399 ytmp[j] = ydata[index[j]]; 0400 tmpindex[j]= index[j]; 0401 /*printf("%g %g\n", xtmp[j], ytmp[j]);*/ 0402 } 0403 size_t tmpnout = nsl_geom_linesim_perpdist(xtmp, ytmp, nout, tol, tmpindex); 0404 for (j = 0; j < tmpnout; j++) { 0405 index[j] = index[tmpindex[j]]; 0406 /*printf("tmpindex[%d]: %d\n", j, tmpindex[j]);*/ 0407 } 0408 0409 if (tmpnout == nout) /* return if nout does not change anymore */ 0410 break; 0411 else 0412 nout = tmpnout; 0413 } 0414 0415 free(tmpindex); 0416 free(xtmp); 0417 free(ytmp); 0418 0419 return nout; 0420 } 0421 0422 size_t nsl_geom_linesim_interp(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { 0423 size_t i, nout = 0; 0424 0425 /*first point*/ 0426 index[nout++] = 0; 0427 0428 size_t key = 0; 0429 for (i = 1; i < n-1; i++) { 0430 /*printf("%d: %d-%d\n", i, key, i+1);*/ 0431 double dist = nsl_geom_point_line_dist_y(xdata[key], ydata[key], xdata[i+1], ydata[i+1], xdata[i], ydata[i]); 0432 /*printf("%d: dist = %g\n", i, dist);*/ 0433 if (dist > tol) { 0434 /*printf("take it %d\n", i);*/ 0435 index[nout++] = key = i; 0436 } 0437 } 0438 0439 /* last point */ 0440 index[nout++] = n-1; 0441 0442 return nout; 0443 } 0444 size_t nsl_geom_linesim_interp_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { 0445 double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); 0446 return nsl_geom_linesim_interp(xdata, ydata, n, tol, index); 0447 } 0448 0449 size_t nsl_geom_linesim_visvalingam_whyatt(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { 0450 if (n < 3) /* we need at least three points */ 0451 return 0; 0452 0453 size_t i, nout = n; 0454 double *area = (double *) malloc((n-2)*sizeof(double)); /* area associated with every point */ 0455 if (area == NULL) { 0456 printf("nsl_geom_linesim_visvalingam_whyatt(): ERROR allocating memory!\n"); 0457 return 0; 0458 } 0459 for (i = 0; i < n; i++) 0460 index[i] = i; 0461 for (i = 1; i < n-1; i++) 0462 area[i-1] = nsl_geom_three_point_area(xdata[i-1], ydata[i-1], xdata[i], ydata[i], xdata[i+1], ydata[i+1]); 0463 0464 size_t minindex; 0465 while ( nsl_stats_minimum(area, n-2, &minindex) < tol && nout > 2) { 0466 /* double minarea; 0467 while ( (minarea = nsl_stats_minimum(area, n-2, &minindex)) < tol && nout > 2) { */ 0468 /*for (i=0; i < n-2; i++) 0469 if (area[i]<DBL_MAX) 0470 printf("area[%zu] = %g\n", i, area[i]); 0471 */ 0472 /* remove point minindex */ 0473 /*printf("removing point %zu (minindex = %zu, minarea = %g) nout=%zu\n", minindex+1, minindex, minarea, nout-1);*/ 0474 index[minindex+1] = 0; 0475 area[minindex] = DBL_MAX; 0476 double tmparea; 0477 /* update area of neigbor points */ 0478 size_t before = minindex, after = minindex+2; /* find index before and after */ 0479 while (index[before] == 0 && before > 0) 0480 before--; 0481 while (after < n+1 && index[after] == 0) 0482 after++; 0483 if (minindex > 0) { /*before */ 0484 if (before > 0) { 0485 size_t beforebefore=before-1; 0486 while (index[beforebefore] == 0 && beforebefore > 0) 0487 beforebefore--; 0488 /*printf("recalculate area[%zu] from %zu %zu %zu\n", before-1, beforebefore, before, after);*/ 0489 tmparea = nsl_geom_three_point_area(xdata[beforebefore], ydata[beforebefore], xdata[before], ydata[before], xdata[after], ydata[after]); 0490 if (tmparea > area[before-1]) /* take largest value of new and old area */ 0491 area[before-1] = tmparea; 0492 } 0493 } 0494 if (minindex < n-3) { /* after */ 0495 if (after < n-1) { 0496 size_t afterafter = after+1; 0497 while (afterafter < n-1 && index[afterafter] == 0) 0498 afterafter++; 0499 /*printf("recalculate area[%zu] from %zu %zu %zu\n",after-1, before, after, afterafter);*/ 0500 tmparea = nsl_geom_three_point_area(xdata[before], ydata[before], xdata[after], ydata[after], xdata[afterafter], ydata[afterafter]); 0501 if (tmparea > area[after-1]) /* take largest value of new and old area */ 0502 area[after-1] = tmparea; 0503 } 0504 } 0505 nout--; 0506 }; 0507 0508 /*for(i=0;i<n;i++) 0509 printf("INDEX = %d\n", index[i]); 0510 */ 0511 0512 /* condens index */ 0513 i = 1; 0514 size_t newi = 1; 0515 while (newi < n-1) { 0516 while (index[newi] == 0) 0517 newi++; 0518 index[i++] = index[newi++]; 0519 }; 0520 0521 /*for(i=0;i<nout;i++) 0522 printf("INDEX2 = %d\n", index[i]); 0523 */ 0524 0525 free(area); 0526 return nout; 0527 } 0528 size_t nsl_geom_linesim_visvalingam_whyatt_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { 0529 double tol = nsl_geom_linesim_clip_area_perpoint(xdata, ydata, n); 0530 0531 return nsl_geom_linesim_visvalingam_whyatt(xdata, ydata, n, tol, index); 0532 } 0533 0534 size_t nsl_geom_linesim_reumann_witkam(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { 0535 size_t i, nout = 0, key = 0, key2 = 1; 0536 0537 /*first point*/ 0538 index[nout++] = 0; 0539 0540 for (i = 2; i < n-1; i++) { 0541 /* distance to line key -- key2 */ 0542 double dist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[key2], ydata[key2], xdata[i], ydata[i]); 0543 /*printf("%d: %g\n", i, dist);*/ 0544 0545 if (dist > tol) { /* take it */ 0546 /*printf("%d: take it\n", i);*/ 0547 key = i-1; 0548 key2 = i; 0549 index[nout++] = i-1; 0550 } 0551 } 0552 0553 /* last point */ 0554 index[nout++] = n-1; 0555 0556 return nout; 0557 } 0558 size_t nsl_geom_linesim_reumann_witkam_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { 0559 double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); 0560 return nsl_geom_linesim_reumann_witkam(xdata, ydata, n, tol, index); 0561 } 0562 0563 size_t nsl_geom_linesim_opheim(const double xdata[], const double ydata[], const size_t n, const double mintol, const double maxtol, size_t index[]) { 0564 size_t i, nout = 0, key = 0, key2; 0565 0566 /*first point*/ 0567 index[nout++] = 0; 0568 0569 for (i = 1; i < n-1; i++) { 0570 double dist, perpdist; 0571 do { /* find key2 */ 0572 dist = nsl_geom_point_point_dist(xdata[key], ydata[key], xdata[i], ydata[i]); 0573 /*printf("dist = %g: %d-%d\n", dist, key, i);*/ 0574 i++; 0575 } while (dist < mintol); 0576 i--; 0577 if (key == i-1) /*i+1 outside mintol */ 0578 key2 = i; 0579 else 0580 key2 = i-1; /* last point inside */ 0581 /*printf("found key2 @%d\n", key2);*/ 0582 0583 do { /* find next key */ 0584 dist = nsl_geom_point_point_dist(xdata[key], ydata[key], xdata[i], ydata[i]); 0585 perpdist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[key2], ydata[key2], xdata[i], ydata[i]); 0586 /*printf("dist = %g, perpdist=%g: %d\n", dist, perpdist, i);*/ 0587 i++; 0588 } while (dist < maxtol && perpdist < mintol); 0589 i--; 0590 if (key == i-1) /*i+1 outside */ 0591 key = i; 0592 else { 0593 key = i-1; 0594 i--; 0595 } 0596 /*printf("new key: %d\n", key);*/ 0597 index[nout++] = key; 0598 } 0599 0600 /* last point */ 0601 if (index[nout-1] != n-1) 0602 index[nout++] = n-1; 0603 0604 return nout; 0605 } 0606 size_t nsl_geom_linesim_opheim_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { 0607 double mintol = 10.*nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); 0608 /* TODO: calculate max tolerance ? */ 0609 double maxtol = 5.*mintol; 0610 0611 return nsl_geom_linesim_opheim(xdata, ydata, n, mintol, maxtol, index); 0612 } 0613 0614 size_t nsl_geom_linesim_lang(const double xdata[], const double ydata[], const size_t n, const double tol, const size_t region, size_t index[]) { 0615 size_t i, j, nout = 0, key = 0; 0616 0617 /*first point*/ 0618 index[nout++] = 0; 0619 0620 double dist, maxdist; 0621 for (i = 1; i < n-1; i++) { 0622 size_t tmpregion=region; 0623 if (key+tmpregion > n-1) /* end of data set */ 0624 tmpregion = n-1-key; 0625 0626 do { 0627 maxdist = 0; 0628 for (j = 1; j < tmpregion; j++) { 0629 dist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[key+tmpregion], ydata[key+tmpregion], xdata[key+j], ydata[key+j]); 0630 /*printf("%d: dist (%d to %d-%d) = %g\n", j, key+j, key, key+tmpregion, dist);*/ 0631 if (dist > maxdist) 0632 maxdist = dist; 0633 } 0634 /*printf("tol = %g maxdist = %g\n", tol, maxdist);*/ 0635 tmpregion--; 0636 /*printf("region = %d\n", tmpregion);*/ 0637 } while (maxdist>tol && tmpregion>0); 0638 i += tmpregion; 0639 index[nout++] = key = i; 0640 /*printf("take it (%d) key=%d\n", i, key);*/ 0641 } 0642 0643 /* last point */ 0644 if (index[nout-1] != n-1) 0645 index[nout++] = n-1; 0646 0647 return nout; 0648 } 0649 size_t nsl_geom_linesim_lang_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { 0650 double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); 0651 /* TODO: calculate search region */ 0652 size_t region = 0; 0653 printf("nsl_geom_linesim_lang_auto(): Not implemented yet\n"); 0654 0655 return nsl_geom_linesim_lang(xdata, ydata, n, tol, region, index); 0656 } 0657