labplot  2.8.2
About: LabPlot is an application for plotting and analysis of 2D and 3D functions and data. It is a complete rewrite of LabPlot1 and lacks in the first release a lot of features available in the predecessor. On the other hand, the GUI and the usability is more superior.
  Fossies Dox: labplot-2.8.2.tar.gz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

nsl_geom_linesim.c
Go to the documentation of this file.
1 /***************************************************************************
2  File : nsl_geom_linesim.c
3  Project : LabPlot
4  Description : NSL geometry line simplification functions
5  --------------------------------------------------------------------
6  Copyright : (C) 2016-2019 by Stefan Gerlach (stefan.gerlach@uni.kn)
7 
8  ***************************************************************************/
9 
10 /***************************************************************************
11  * *
12  * This program is free software; you can redistribute it and/or modify *
13  * it under the terms of the GNU General Public License as published by *
14  * the Free Software Foundation; either version 2 of the License, or *
15  * (at your option) any later version. *
16  * *
17  * This program is distributed in the hope that it will be useful, *
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
20  * GNU General Public License for more details. *
21  * *
22  * You should have received a copy of the GNU General Public License *
23  * along with this program; if not, write to the Free Software *
24  * Foundation, Inc., 51 Franklin Street, Fifth Floor, *
25  * Boston, MA 02110-1301 USA *
26  * *
27  ***************************************************************************/
28 
29 #include "nsl_geom_linesim.h"
30 #include "nsl_geom.h"
31 #include "nsl_common.h"
32 #include "nsl_sort.h"
33 #include "nsl_stats.h"
34 
35 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"),
36  i18n("radial distance"), i18n("Interpolation"), i18n("Opheim"), i18n("Lang")};
37 
38 /*********** error calculation functions *********/
39 
40 double nsl_geom_linesim_positional_error(const double xdata[], const double ydata[], const size_t n, const size_t index[]) {
41  double dist = 0;
42  size_t i = 0, j; /* i: index of index[] */
43  do {
44  /*for every point not in index[] calculate distance to line*/
45  /*printf("i=%d (index[i]-index[i+1]=%d-%d)\n", i , index[i], index[i+1]);*/
46  for (j = 1; j < index[i+1]-index[i]; j++) {
47  /*printf("i=%d: j=%d\n", i, j);*/
48  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]);
49  /*printf("dist = %g\n", dist);*/
50  }
51  i++;
52  } while(index[i] != n-1);
53 
54  return dist/(double)n;
55 }
56 
57 double nsl_geom_linesim_positional_squared_error(const double xdata[], const double ydata[], const size_t n, const size_t index[]) {
58  double dist = 0;
59  size_t i = 0, j; /* i: index of index[] */
60  do {
61  /*for every point not in index[] calculate distance to line*/
62  /*printf("i=%d (index[i]-index[i+1]=%d-%d)\n", i , index[i], index[i+1]);*/
63  for (j = 1; j < index[i+1]-index[i]; j++) {
64  /*printf("i=%d: j=%d\n", i, j);*/
65  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);
66  /*printf("dist = %g\n", dist);*/
67  }
68  i++;
69  } while(index[i] != n-1);
70 
71  return dist/(double)n;
72 }
73 
74 double nsl_geom_linesim_area_error(const double xdata[], const double ydata[], const size_t n, const size_t index[]) {
75  double area = 0;
76 
77  size_t i = 0, j; /* i: index of index[] */
78  do {
79  /* for every point not in index[] calculate area */
80  /*printf("i=%d (index[i]-index[i+1]=%d-%d)\n", i , index[i], index[i+1]);*/
81  for (j = 1; j < index[i+1]-index[i]; j++) {
82  /*printf("j=%d: area between %d %d %d\n", j, index[i]+j-1, index[i]+j, index[i+1]);*/
83  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]]);
84  /*printf("area = %g\n", area);*/
85  }
86  i++;
87  } while(index[i] != n-1);
88 
89 
90  return area/(double)n;
91 }
92 
93 double nsl_geom_linesim_clip_diag_perpoint(const double xdata[], const double ydata[], const size_t n) {
94  double dx = nsl_stats_maximum(xdata, n, NULL) - nsl_stats_minimum(xdata, n, NULL);
95  double dy = nsl_stats_maximum(ydata, n, NULL) - nsl_stats_minimum(ydata, n, NULL);
96  double d = sqrt(dx*dx+dy*dy);
97 
98  return d/(double)n; /* per point */
99 }
100 
101 double nsl_geom_linesim_clip_area_perpoint(const double xdata[], const double ydata[], const size_t n) {
102  double dx = nsl_stats_maximum(xdata, n, NULL) - nsl_stats_minimum(xdata, n, NULL);
103  double dy = nsl_stats_maximum(ydata, n, NULL) - nsl_stats_minimum(ydata, n, NULL);
104  double A = dx*dy;
105 
106  return A/(double)n; /* per point */
107 }
108 
109 double nsl_geom_linesim_avg_dist_perpoint(const double xdata[], const double ydata[], const size_t n) {
110  double dist = 0;
111  size_t i;
112  for (i = 0; i < n-1; i++) {
113  double dx = xdata[i+1] - xdata[i];
114  double dy = ydata[i+1] - ydata[i];
115  dist += sqrt(dx*dx + dy*dy);
116  }
117  dist /= (double)n;
118 
119  return dist;
120 }
121 
122 /*********** simplification algorithms *********/
123 
124 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[]) {
125  /*printf("DP: %d - %d\n", start, end);*/
126 
127  size_t i, nkey = start;
128  double maxdist = 0;
129  /* search for key (biggest perp. distance) */
130  for (i = start+1; i < end; i++) {
131  double dist = nsl_geom_point_line_dist(xdata[start], ydata[start], xdata[end], ydata[end], xdata[i], ydata[i]);
132  if (dist > maxdist) {
133  maxdist = dist;
134  nkey = i;
135  }
136  }
137  /*printf("maxdist = %g @ i = %zu\n", maxdist, nkey);*/
138 
139  if (maxdist > tol) {
140  /*printf("take %d\n", nkey);*/
141  index[(*nout)++] = nkey;
142  if (nkey-start > 1)
143  nsl_geom_linesim_douglas_peucker_step(xdata, ydata, start, nkey, nout, tol, index);
144  if (end-nkey > 1)
145  nsl_geom_linesim_douglas_peucker_step(xdata, ydata, nkey, end, nout, tol, index);
146  }
147 
148  /*printf("nout=%d\n", *nout);*/
149 }
150 
151 size_t nsl_geom_linesim_douglas_peucker(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) {
152  size_t nout = 0;
153 
154  /*first point*/
155  index[nout++] = 0;
156 
157  nsl_geom_linesim_douglas_peucker_step(xdata, ydata, 0, n-1, &nout, tol, index);
158 
159  /* last point */
160  if (index[nout-1] != n-1)
161  index[nout++] = n-1;
162 
163  /* sort array index */
164  nsl_sort_size_t(index, nout);
165 
166  return nout;
167 }
168 size_t nsl_geom_linesim_douglas_peucker_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) {
169  double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n);
170  return nsl_geom_linesim_douglas_peucker(xdata, ydata, n, tol, index);
171 }
172 
173 /*
174  * Douglas-Peucker variant:
175  * The key of all egdes of the current simplified line is calculated and only the
176  * largest is added. This is repeated until nout is reached.
177  * */
178 double nsl_geom_linesim_douglas_peucker_variant(const double xdata[], const double ydata[], const size_t n, const size_t nout, size_t index[]) {
179  size_t i;
180  if (nout >= n) { /* use all points */
181  for (i = 0; i < n; i++)
182  index[i] = i;
183  return 0;
184  }
185 
186  /* set first and last point in index (other indizes not initialized) */
187  size_t ncount = 0;
188  index[ncount++] = 0;
189  index[ncount++] = n-1;
190 
191  if (nout <= 2) /* use only first and last point (perp. dist is zero) */
192  return 0.0;
193 
194  double *dist = (double *)malloc(n * sizeof(double));
195  if (dist == NULL) {
196  /* printf("nsl_geom_linesim_douglas_peucker_variant(): ERROR allocating memory for 'dist'!\n"); */
197  return DBL_MAX;
198  }
199 
200  double *maxdist = (double *)malloc(nout * sizeof(double)); /* max dist per edge */
201  if (maxdist == NULL) {
202  /* printf("nsl_geom_linesim_douglas_peucker_variant(): ERROR allocating memory for 'maxdist'!\n"); */
203  free(dist);
204  return DBL_MAX;
205  }
206 
207  for (i = 0; i < n; i++) { /* initialize dist */
208  dist[i] = nsl_geom_point_line_dist(xdata[0], ydata[0], xdata[n-1], ydata[n-1], xdata[i], ydata[i]);
209  /* printf("nsl_geom_linesim_douglas_peucker_variant(): %zu dist = %g\n", i, dist[i]); */
210  }
211  for (i = 0; i < nout; i++)
212  maxdist[i] = 0;
213 
214  double newmaxdist = 0;
215  do {
216  /* printf("nsl_geom_linesim_douglas_peucker_variant(): NEXT ITERATION\n"); */
217  size_t key = 0;
218 
219  /* find edge of maximum */
220  size_t maxindex;
221  nsl_stats_maximum(maxdist, ncount, &maxindex);
222  /* printf("nsl_geom_linesim_douglas_peucker_variant(): found edge with max dist at index %zu\n", maxindex); */
223  /*newmaxdist = nsl_stats_maximum(dist, n, &key);*/
224  newmaxdist = 0;
225  for (i = index[maxindex]+1; i < index[maxindex+1]; i++) {
226  /* printf("nsl_geom_linesim_douglas_peucker_variant(): iterate i=%zu\n", i); */
227  if (dist[i] > newmaxdist) {
228  newmaxdist = dist[i];
229  key = i;
230  }
231  }
232 
233  /* printf("nsl_geom_linesim_douglas_peucker_variant(): found key %zu (dist = %g)\n", key, newmaxdist); */
234  ncount++;
235  dist[key] = 0;
236 
237  /* find index of previous key */
238  size_t previndex = 0;
239  while (index[previndex+1] < key)
240  previndex++;
241  /* printf("nsl_geom_linesim_douglas_peucker_variant(): previndex = %zu (update keys %zu - %zu)\n", previndex, index[previndex], index[previndex+1]); */
242 
243  size_t v;
244  /* printf("nsl_geom_linesim_douglas_peucker_variant(): ncount = %zu, previndex = %zu\n", ncount, previndex); */
245  /* update dist[]. no update on last key */
246  if (ncount < nout) {
247  /* shift maxdist */
248  for (v = ncount; v > previndex; v--)
249  maxdist[v] = maxdist[v-1];
250 
251  double tmpmax = 0;
252  for (v = index[previndex]+1; v < key; v++) {
253  /* printf("nsl_geom_linesim_douglas_peucker_variant(): %zu in %zu - %zu", v, index[previndex], key); */
254  dist[v] = nsl_geom_point_line_dist(xdata[index[previndex]], ydata[index[previndex]], xdata[key], ydata[key],
255  xdata[v], ydata[v]);
256  if (dist[v] > tmpmax)
257  tmpmax = dist[v];
258 
259  /* printf(" dist = %g\n", dist[v]); */
260  }
261  maxdist[previndex] = tmpmax;
262 
263  tmpmax = 0;
264  for (v = key+1; v < index[previndex+1]; v++) {
265  /* printf("nsl_geom_linesim_douglas_peucker_variant(): %zu in %zu - %zu", v, key, index[previndex+1]); */
266  dist[v] = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[index[previndex+1]], ydata[index[previndex+1]],
267  xdata[v], ydata[v]);
268  if (dist[v] > tmpmax)
269  tmpmax = dist[v];
270  /* printf(" dist = %g\n", dist[v]); */
271  }
272  maxdist[previndex+1] = tmpmax;
273  }
274 
275  /* put into index array */
276  for (v = ncount; v > previndex+1; v--)
277  index[v] = index[v-1];
278  index[previndex+1] = key;
279  } while (ncount < nout);
280 
281  free(maxdist);
282  free(dist);
283 
284  return newmaxdist;
285 }
286 
287 size_t nsl_geom_linesim_nthpoint(const size_t n, const int step, size_t index[]) {
288  if (step < 1) {
289  printf("step size must be > 0 (given: %d)\n", step);
290  return 0;
291  }
292 
293  size_t i, nout = 0;
294 
295  /*first point*/
296  index[nout++] = 0;
297 
298  for (i = 1; i < n-1; i++)
299  if (i%step == 0)
300  index[nout++] = i;
301 
302  /* last point */
303  index[nout++] = n-1;
304 
305  return nout;
306 }
307 
308 size_t nsl_geom_linesim_raddist(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) {
309  size_t i, nout = 0, key = 0;
310 
311  /*first point*/
312  index[nout++] = 0;
313 
314  for (i = 1; i < n-1; i++) {
315  /* distance to key point */
316  double dist = nsl_geom_point_point_dist(xdata[i], ydata[i], xdata[key], ydata[key]);
317  /* distance to last point */
318  double lastdist = nsl_geom_point_point_dist(xdata[i], ydata[i], xdata[n-1], ydata[n-1]);
319  /*printf("%d: %g %g\n", i, dist, lastdist);*/
320 
321  if (dist > tol && lastdist > tol) {
322  index[nout++] = i;
323  key = i;
324  }
325  }
326 
327  /* last point */
328  index[nout++] = n-1;
329 
330  return nout;
331 }
332 size_t nsl_geom_linesim_raddist_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) {
333  double tol = 10.*nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n);
334  return nsl_geom_linesim_raddist(xdata, ydata, n, tol, index);
335 }
336 
337 size_t nsl_geom_linesim_perpdist(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) {
338  size_t nout = 0, key = 0, i;
339 
340  /*first point*/
341  index[nout++] = 0;
342 
343  for (i = 1; i < n-1; i++) {
344  /* distance of point i to line key -- i+1 */
345  double dist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[i+1], ydata[i+1], xdata[i], ydata[i]);
346  /*printf("%d: %g\n", i, dist);*/
347 
348  if (dist > tol) { /* take it */
349  index[nout++] = i;
350  key = i;
351  /*printf("%d: take it (key = %d)\n", i, key);*/
352  } else { /* ignore it */
353  if (i+1 < n-1) /* last point is taken anyway */
354  index[nout++] = i+1; /* take next point in any case */
355  /*printf("%d: ignore it (key = %d)\n", i, i+1);*/
356  key = ++i;
357  }
358  }
359 
360  /* last point */
361  index[nout++] = n-1;
362 
363  return nout;
364 }
365 size_t nsl_geom_linesim_perpdist_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) {
366  double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n);
367  return nsl_geom_linesim_perpdist(xdata, ydata, n, tol, index);
368 }
369 
370 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[]) {
371  size_t i, j, nout;
372  double *xtmp = (double *) malloc(n*sizeof(double));
373  if (xtmp == NULL) {
374  printf("nsl_geom_linesim_perpdist_repeat(): ERROR allocating memory for 'xtmp'!\n");
375  return 0;
376  }
377 
378  double *ytmp = (double *) malloc(n*sizeof(double));
379  if (ytmp == NULL) {
380  printf("nsl_geom_linesim_perpdist_repeat(): ERROR allocating memory for 'ytmp'!\n");
381  free(xtmp);
382  return 0;
383  }
384 
385  size_t *tmpindex = (size_t *) malloc(n*sizeof(size_t));
386  if (tmpindex == NULL) {
387  printf("nsl_geom_linesim_perpdist_repeat(): ERROR allocating memory for 'tmpindex'!\n");
388  free(xtmp);
389  free(ytmp);
390  return 0;
391  }
392 
393  /* first round */
394  nout = nsl_geom_linesim_perpdist(xdata, ydata, n, tol, index);
395  /* repeats */
396  for (i = 0; i < repeat - 1; i++) {
397  for (j = 0; j < nout; j++) {
398  xtmp[j] = xdata[index[j]];
399  ytmp[j] = ydata[index[j]];
400  tmpindex[j]= index[j];
401  /*printf("%g %g\n", xtmp[j], ytmp[j]);*/
402  }
403  size_t tmpnout = nsl_geom_linesim_perpdist(xtmp, ytmp, nout, tol, tmpindex);
404  for (j = 0; j < tmpnout; j++) {
405  index[j] = index[tmpindex[j]];
406  /*printf("tmpindex[%d]: %d\n", j, tmpindex[j]);*/
407  }
408 
409  if (tmpnout == nout) /* return if nout does not change anymore */
410  break;
411  else
412  nout = tmpnout;
413  }
414 
415  free(tmpindex);
416  free(xtmp);
417  free(ytmp);
418 
419  return nout;
420 }
421 
422 size_t nsl_geom_linesim_interp(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) {
423  size_t i, nout = 0;
424 
425  /*first point*/
426  index[nout++] = 0;
427 
428  size_t key = 0;
429  for (i = 1; i < n-1; i++) {
430  /*printf("%d: %d-%d\n", i, key, i+1);*/
431  double dist = nsl_geom_point_line_dist_y(xdata[key], ydata[key], xdata[i+1], ydata[i+1], xdata[i], ydata[i]);
432  /*printf("%d: dist = %g\n", i, dist);*/
433  if (dist > tol) {
434  /*printf("take it %d\n", i);*/
435  index[nout++] = key = i;
436  }
437  }
438 
439  /* last point */
440  index[nout++] = n-1;
441 
442  return nout;
443 }
444 size_t nsl_geom_linesim_interp_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) {
445  double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n);
446  return nsl_geom_linesim_interp(xdata, ydata, n, tol, index);
447 }
448 
449 size_t nsl_geom_linesim_visvalingam_whyatt(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) {
450  if (n < 3) /* we need at least three points */
451  return 0;
452 
453  size_t i, nout = n;
454  double *area = (double *) malloc((n-2)*sizeof(double)); /* area associated with every point */
455  if (area == NULL) {
456  printf("nsl_geom_linesim_visvalingam_whyatt(): ERROR allocating memory!\n");
457  return 0;
458  }
459  for (i = 0; i < n; i++)
460  index[i] = i;
461  for (i = 1; i < n-1; i++)
462  area[i-1] = nsl_geom_three_point_area(xdata[i-1], ydata[i-1], xdata[i], ydata[i], xdata[i+1], ydata[i+1]);
463 
464  size_t minindex;
465  while ( nsl_stats_minimum(area, n-2, &minindex) < tol && nout > 2) {
466  /* double minarea;
467  while ( (minarea = nsl_stats_minimum(area, n-2, &minindex)) < tol && nout > 2) { */
468  /*for (i=0; i < n-2; i++)
469  if (area[i]<DBL_MAX)
470  printf("area[%zu] = %g\n", i, area[i]);
471  */
472  /* remove point minindex */
473  /*printf("removing point %zu (minindex = %zu, minarea = %g) nout=%zu\n", minindex+1, minindex, minarea, nout-1);*/
474  index[minindex+1] = 0;
475  area[minindex] = DBL_MAX;
476  double tmparea;
477  /* update area of neigbor points */
478  size_t before = minindex, after = minindex+2; /* find index before and after */
479  while (index[before] == 0 && before > 0)
480  before--;
481  while (after < n+1 && index[after] == 0)
482  after++;
483  if (minindex > 0) { /*before */
484  if (before > 0) {
485  size_t beforebefore=before-1;
486  while (index[beforebefore] == 0 && beforebefore > 0)
487  beforebefore--;
488  /*printf("recalculate area[%zu] from %zu %zu %zu\n", before-1, beforebefore, before, after);*/
489  tmparea = nsl_geom_three_point_area(xdata[beforebefore], ydata[beforebefore], xdata[before], ydata[before], xdata[after], ydata[after]);
490  if (tmparea > area[before-1]) /* take largest value of new and old area */
491  area[before-1] = tmparea;
492  }
493  }
494  if (minindex < n-3) { /* after */
495  if (after < n-1) {
496  size_t afterafter = after+1;
497  while (afterafter < n-1 && index[afterafter] == 0)
498  afterafter++;
499  /*printf("recalculate area[%zu] from %zu %zu %zu\n",after-1, before, after, afterafter);*/
500  tmparea = nsl_geom_three_point_area(xdata[before], ydata[before], xdata[after], ydata[after], xdata[afterafter], ydata[afterafter]);
501  if (tmparea > area[after-1]) /* take largest value of new and old area */
502  area[after-1] = tmparea;
503  }
504  }
505  nout--;
506  };
507 
508  /*for(i=0;i<n;i++)
509  printf("INDEX = %d\n", index[i]);
510  */
511 
512  /* condens index */
513  i = 1;
514  size_t newi = 1;
515  while (newi < n-1) {
516  while (index[newi] == 0)
517  newi++;
518  index[i++] = index[newi++];
519  };
520 
521  /*for(i=0;i<nout;i++)
522  printf("INDEX2 = %d\n", index[i]);
523  */
524 
525  free(area);
526  return nout;
527 }
528 size_t nsl_geom_linesim_visvalingam_whyatt_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) {
529  double tol = nsl_geom_linesim_clip_area_perpoint(xdata, ydata, n);
530 
531  return nsl_geom_linesim_visvalingam_whyatt(xdata, ydata, n, tol, index);
532 }
533 
534 size_t nsl_geom_linesim_reumann_witkam(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) {
535  size_t i, nout = 0, key = 0, key2 = 1;
536 
537  /*first point*/
538  index[nout++] = 0;
539 
540  for (i = 2; i < n-1; i++) {
541  /* distance to line key -- key2 */
542  double dist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[key2], ydata[key2], xdata[i], ydata[i]);
543  /*printf("%d: %g\n", i, dist);*/
544 
545  if (dist > tol) { /* take it */
546  /*printf("%d: take it\n", i);*/
547  key = i-1;
548  key2 = i;
549  index[nout++] = i-1;
550  }
551  }
552 
553  /* last point */
554  index[nout++] = n-1;
555 
556  return nout;
557 }
558 size_t nsl_geom_linesim_reumann_witkam_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) {
559  double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n);
560  return nsl_geom_linesim_reumann_witkam(xdata, ydata, n, tol, index);
561 }
562 
563 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[]) {
564  size_t i, nout = 0, key = 0, key2;
565 
566  /*first point*/
567  index[nout++] = 0;
568 
569  for (i = 1; i < n-1; i++) {
570  double dist, perpdist;
571  do { /* find key2 */
572  dist = nsl_geom_point_point_dist(xdata[key], ydata[key], xdata[i], ydata[i]);
573  /*printf("dist = %g: %d-%d\n", dist, key, i);*/
574  i++;
575  } while (dist < mintol);
576  i--;
577  if (key == i-1) /*i+1 outside mintol */
578  key2 = i;
579  else
580  key2 = i-1; /* last point inside */
581  /*printf("found key2 @%d\n", key2);*/
582 
583  do { /* find next key */
584  dist = nsl_geom_point_point_dist(xdata[key], ydata[key], xdata[i], ydata[i]);
585  perpdist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[key2], ydata[key2], xdata[i], ydata[i]);
586  /*printf("dist = %g, perpdist=%g: %d\n", dist, perpdist, i);*/
587  i++;
588  } while (dist < maxtol && perpdist < mintol);
589  i--;
590  if (key == i-1) /*i+1 outside */
591  key = i;
592  else {
593  key = i-1;
594  i--;
595  }
596  /*printf("new key: %d\n", key);*/
597  index[nout++] = key;
598  }
599 
600  /* last point */
601  if (index[nout-1] != n-1)
602  index[nout++] = n-1;
603 
604  return nout;
605 }
606 size_t nsl_geom_linesim_opheim_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) {
607  double mintol = 10.*nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n);
608  /* TODO: calculate max tolerance ? */
609  double maxtol = 5.*mintol;
610 
611  return nsl_geom_linesim_opheim(xdata, ydata, n, mintol, maxtol, index);
612 }
613 
614 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[]) {
615  size_t i, j, nout = 0, key = 0;
616 
617  /*first point*/
618  index[nout++] = 0;
619 
620  double dist, maxdist;
621  for (i = 1; i < n-1; i++) {
622  size_t tmpregion=region;
623  if (key+tmpregion > n-1) /* end of data set */
624  tmpregion = n-1-key;
625 
626  do {
627  maxdist = 0;
628  for (j = 1; j < tmpregion; j++) {
629  dist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[key+tmpregion], ydata[key+tmpregion], xdata[key+j], ydata[key+j]);
630  /*printf("%d: dist (%d to %d-%d) = %g\n", j, key+j, key, key+tmpregion, dist);*/
631  if (dist > maxdist)
632  maxdist = dist;
633  }
634  /*printf("tol = %g maxdist = %g\n", tol, maxdist);*/
635  tmpregion--;
636  /*printf("region = %d\n", tmpregion);*/
637  } while (maxdist>tol && tmpregion>0);
638  i += tmpregion;
639  index[nout++] = key = i;
640  /*printf("take it (%d) key=%d\n", i, key);*/
641  }
642 
643  /* last point */
644  if (index[nout-1] != n-1)
645  index[nout++] = n-1;
646 
647  return nout;
648 }
649 size_t nsl_geom_linesim_lang_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) {
650  double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n);
651  /* TODO: calculate search region */
652  size_t region = 0;
653  printf("nsl_geom_linesim_lang_auto(): Not implemented yet\n");
654 
655  return nsl_geom_linesim_lang(xdata, ydata, n, tol, region, index);
656 }
657 
#define i18n(m)
Definition: nsl_common.h:38
double nsl_geom_point_line_dist(double x1, double y1, double x2, double y2, double xp, double yp)
Definition: nsl_geom.c:36
double nsl_geom_point_point_dist(double x1, double y1, double x2, double y2)
Definition: nsl_geom.c:32
double nsl_geom_three_point_area(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: nsl_geom.c:44
double nsl_geom_point_line_dist_y(double x1, double y1, double x2, double y2, double xp, double yp)
Definition: nsl_geom.c:40
double nsl_geom_linesim_positional_error(const double xdata[], const double ydata[], const size_t n, const size_t index[])
double nsl_geom_linesim_positional_squared_error(const double xdata[], const double ydata[], const size_t n, const size_t index[])
size_t nsl_geom_linesim_raddist_auto(const double xdata[], const double ydata[], const size_t n, size_t index[])
double nsl_geom_linesim_clip_area_perpoint(const double xdata[], const double ydata[], const size_t n)
size_t nsl_geom_linesim_douglas_peucker_auto(const double xdata[], const double ydata[], const size_t n, size_t index[])
size_t nsl_geom_linesim_raddist(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[])
size_t nsl_geom_linesim_perpdist_auto(const double xdata[], const double ydata[], const size_t n, size_t index[])
size_t nsl_geom_linesim_reumann_witkam_auto(const double xdata[], const double ydata[], const size_t n, size_t index[])
double nsl_geom_linesim_avg_dist_perpoint(const double xdata[], const double ydata[], const size_t n)
size_t nsl_geom_linesim_opheim_auto(const double xdata[], const double ydata[], const size_t n, size_t index[])
size_t nsl_geom_linesim_interp(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[])
size_t nsl_geom_linesim_perpdist(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[])
double nsl_geom_linesim_clip_diag_perpoint(const double xdata[], const double ydata[], const size_t n)
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[])
size_t nsl_geom_linesim_douglas_peucker(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[])
size_t nsl_geom_linesim_reumann_witkam(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[])
size_t nsl_geom_linesim_visvalingam_whyatt_auto(const double xdata[], const double ydata[], const size_t n, size_t index[])
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[])
size_t nsl_geom_linesim_lang_auto(const double xdata[], const double ydata[], const size_t n, size_t index[])
size_t nsl_geom_linesim_nthpoint(const size_t n, const int step, size_t index[])
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[])
double nsl_geom_linesim_douglas_peucker_variant(const double xdata[], const double ydata[], const size_t n, const size_t nout, size_t index[])
const char * nsl_geom_linesim_type_name[]
size_t nsl_geom_linesim_interp_auto(const double xdata[], const double ydata[], const size_t n, size_t index[])
size_t nsl_geom_linesim_visvalingam_whyatt(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[])
double nsl_geom_linesim_area_error(const double xdata[], const double ydata[], const size_t n, const size_t index[])
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[])
void nsl_sort_size_t(size_t array[], const size_t n)
Definition: nsl_sort.c:40
double nsl_stats_minimum(const double data[], const size_t n, size_t *index)
Definition: nsl_stats.c:35
double nsl_stats_maximum(const double data[], const size_t n, size_t *index)
Definition: nsl_stats.c:53