Ramer-Douglas-Peucker line simplification

rdp.c

#include <assert.h>
#include <math.h>
#include <stdio.h>

/* Alternative design: Iterate over path and at each point i check distance to line with points
   i + 1 and i + 2 and skip i + 1 when distance is small than epsilon. */
static float
line_dist(const float *p, const float *p1, const float *p2) {
  float dx = p2[0] - p1[0];
  float dy = p2[1] - p1[1];
  float d = sqrtf(dx * dx + dy * dy);
  return fabs(p[0] * dy - p[1] * dx + p2[0] * p1[1] - p2[1] * p1[0])/d;
}
static void
max_pnt_dist(float *max_dist, int *idx, const float *p, int n) {
  for (int i = 1; i + 1 < n; ++i) {
      float d = line_dist(p+(i*2), p, p+(n-1)*2);
      if (d > *max_dist) {
          *max_dist = d;
          *idx = i;
      }
  }
}
static int
reduct_path(float *dst, const float *p, int n, float e) {
  assert(e >= 0);
  int idx = 0;
  float max_dist = 0;
  max_pnt_dist(&max_dist, &idx, p, n);
  if (max_dist > e) {
    int n1 = reduct_path(dst, p, idx+1, e);
    dst += (n1 - 1) * 2;
    int n2 = reduct_path(dst, p + idx*2, n - idx, e);
    return n1 + n2 - 1;
  }
  dst[0*2+0] = p[0*2+0];
  dst[0*2+1] = p[0*2+1];
  dst[1*2+0] = p[(n-1)*2+0];
  dst[1*2+1] = p[(n-1)*2+1];
  return 2;
}
static void
print_points(const float *p, int n) {
  printf("count: %d\n", n);
  for (int i = 0; i < n; ++i) {
    if (i > 0) printf(" ");
    printf("(%g, %g)", p[i*2+0], p[i*2+1]);
  }
  printf("\n");
}
extern int
main(void) {
  static const float points[20] = {
    0,0, 1,0.1, 2,-0.1, 3,5, 4,6,
    5,7, 6,8.1, 7,9, 8,9, 9,9
  };
  float out[20];
  int n = reduct_path(out, points, 10, 1.0);
  print_points(out, n);
  return 0;
}