[ create a new paste ] login | about

Link: http://codepad.org/mUVSnJTr    [ raw code | output | fork | 4 comments ]

mohit_at_codepad - C, pasted on Jul 16:
/* Find minimum distance between a point and a line _segment_       */
/* 1. Vector calculation (drop a perpendicular) method              */
/* 2. Ternary search (Enough binary search :) let's move to ternary */
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <assert.h>

#define const

struct point {
  double x, y;
};

struct line {
  struct point start, end;
};

double dist_sq(double x0, double y0, double x1, double y1) {
  const double xDiff = x1 - x0, yDiff = y1 - y0;
  return xDiff * xDiff + yDiff * yDiff;
}

/* Finding perpendicular from pt dropping on the ln */
double find_min_dist1(struct line ln, struct point pt) {
  const double x0 = ln.start.x, y0 = ln.start.y;
  const double x1 = ln.end.x, y1 = ln.start.y;
  const double x2 = pt.x, y2 = pt.y;
  const double xDiff = x1 - x0, yDiff = y1 - y0;
  /* Dividing line into two parts of length in ratio
   * r : (1-r)
   * Joining this point with pt gives the minimum distance
   */
  double r = ( (x2-x0) * xDiff + (y2-y0) * yDiff)
                       / (xDiff * xDiff + yDiff * yDiff);
  /* r = std::min( 1.0, std::max(0.0, r) ); */
  if(r < 0.0) r = 0.0;  // If perpendicular drops off the line,
  if(r > 1.0) r = 1.0;  // pick the nearest point
  {
      const double joinxdif = (x0 + xDiff * r - x2);
      const double joinydif = (y0 + yDiff * r - y2);
      /* printf("DEBUG:: @(%g, %g)\n", x0 + xDiff * r, y0 + yDiff * r); */
      return sqrt( joinxdif * joinxdif + joinydif * joinydif );
  }
}

/* Using ternary search */
double find_min_dist2(struct line ln, struct point pt) {
  double x0 = ln.start.x, y0 = ln.start.y;
  double x1 = ln.end.x, y1 = ln.start.y;
  const double x2 = pt.x, y2 = pt.y;
  const double xDiff = x1 - x0, yDiff = y1 - y0;
  double s = 0.0, e = 1.0;
  while(e - s > DBL_EPSILON) {
    /* Find two ternary points */
    double lft = (s * 2.0 + e) / 3.0;
    double rgt = (e * 2.0 + s) / 3.0;
    if(s == lft) lft = s + DBL_EPSILON;
    if(e == rgt) rgt = e - DBL_EPSILON;
    if(lft >= rgt) break;
    {
        const double xs[] = {ln.start.x + lft * xDiff
                           , ln.start.x + rgt * xDiff};
        const double ys[] = {ln.start.y + lft * yDiff
                           , ln.start.y + rgt * yDiff};
        const double ds[] = { dist_sq(xs[0], ys[0], x2, y2)
                            , dist_sq(xs[1], ys[1], x2, y2) };
        if(ds[0] < ds[1]) e = rgt;
        else s = lft;
    }
  }
  /* printf("DEBUG:: @(%g, %g)\n", ln.start.x + e * xDiff, ln.start.y + e * yDiff); */
  return sqrt( dist_sq(ln.start.x + e * xDiff
                     , ln.start.y + e * yDiff
                     , x2
                     , y2)
              );
}

int main() {
  struct line line1 = { {1.0, 2.0}, {5.0, 6.0} };
  struct point pt1 = {3.0, 4.0};
  struct line line2 = { {1.0, 2.0}, {3.0, 4.0} };
  struct point pt2 = {5.0, 6.0};
  struct line line3 = { {1.0, 2.0}, {11.0, 2.0} };
  struct point pt3 = {5.0, 2.0};
  struct line line4 = { {21.0, 62.0}, {11.0, 5.0} };
  struct point pt4 = {-3.0, -3.0};

  printf( "Distance by method 1 = %.10lf\n", find_min_dist1(line1, pt1) );
  printf( "Distance by method 2 = %.10lf\n", find_min_dist2(line1, pt1) );
  printf("\n");

  printf( "Distance by method 1 = %.10lf\n", find_min_dist1(line2, pt2) );
  printf( "Distance by method 2 = %.10lf\n", find_min_dist2(line2, pt2) );
  printf("\n");

  printf( "Distance by method 1 = %.10lf\n", find_min_dist1(line3, pt3) );
  printf( "Distance by method 2 = %.10lf\n", find_min_dist2(line3, pt3) );
  printf("\n");

  printf( "Distance by method 1 = %.10lf\n", find_min_dist1(line4, pt4) );
  printf( "Distance by method 2 = %.10lf\n", find_min_dist2(line4, pt4) );
  printf("\n");

  return 0;
}

/* End of program */


Output:
1
2
3
4
5
6
7
8
9
10
11
12
Distance by method 1 = 2.0000000000
Distance by method 2 = 2.0000000000

Distance by method 1 = 4.4721359550
Distance by method 2 = 4.4721359550

Distance by method 1 = 0.0000000000
Distance by method 2 = 0.0000000000

Distance by method 1 = 66.4906008395
Distance by method 2 = 66.4906008395



Create a new paste based on this one


Comments:
posted by mohit_at_codepad on Jul 16
In ternary search if you divide the range in golden ration instead of uniformly, you can reuse old (distance) calculation. (I guess) This will reduce some operations.
reply
posted by mohit_at_codepad on Jul 18
If line number 54 behaves as an inifinite loop, change it to something like:
for(i = 0; i < 1000; ++i) ...

This would make sure that the loop terminates.
And 1000 is very very large number (pow(1.5, 1000) can give accuracy of 177 decimal digits) to give a good accuracy.

reply
posted by PhoeniX888 on Jul 19
For line1 and point1 ,
min distance shud b 0 as (3,4) lies on line
{(1,2),(5,6)}...
is my understanding correct ?
reply
posted by mohit_at_codepad on Jul 20
@PhoeniX888
Yes Abhishek you pointed it right. Line number 27 and 50 should be corrected.
Please find the correct answer @
http://codepad.org/gn97VQWx
reply