[GRASS5] [mshapiro@ncsa.uiuc.edu: bug in libes/gis/line_dist.c]

GRASS GRASS developer
1

Hi developers,

find attached a bug report from Michael Shapiro.
Can anyone confirm this bug?

Markus

2

On Thursday 21 February 2002 13:45, Markus Neteler wrote:

Hi developers,

find attached a bug report from Michael Shapiro.
Can anyone confirm this bug?

Markus

Hi Markus,

I was looking through the grass src code and found a bug in the
G_distance2_point_to_line() function in libes/gis/line_dist.c

lines 47,48 should be changed from

  dx = dx * t + x1;
  dy = dy * t + y1;
to
  dx = dx * t;
  dy = dy * t;

This bug has been there ever since I wrote this piece of code,
but only now just discovered it. You should have someone double
check me on this, though, just to be sure.

Michael Shapiro

I think that (dx * t + x1, dy * t + y1) is point x,y projected to line
x1,y1->x2,y2 (comment in code: /* go t from x1,y1 towards x2,y2 */)
and we want distance between this point and and x,y. That is:
dx = dx * t + x1 - x;
dy = dy * t + y1 - y;

Example: line: 1,0 -> 3,0, point 2,3 (distance^2 = 9)
dx = x2 - x1 = 2
dy = y2 - y1 = 0
dx' = x - x1 = 1
dy' = y - y1 = 3
t = (2*1 + 0*3) / (2*2 + 0*0) = 0.5
original:
dx = 2*0.5 + 1 = 2
dy = 0*0.5 + 0 = 0
dx*dx + dy*dy = 2*2 + 0*0 = 4
Shapiro's fix:
dx = 2*0.5 = 1
dy = 0*0.5 = 0
dx*dx + dy*dy = 1*1 + 0*0 = 1
my proposal:
dx = 2*0.5 + 1 - 2 = 0
dy = 0*0.5 + 0 - 3 = -3
dx*dx + dy*dy = (-3)*(-3) + 0*0 = 9

BTW, version from dig_distance2_point_to_line():
  if (t < 0.0) /* go to x1,y1 */
    t = 0.0;
  else if (t > 1.0) /* go to x2,y2 */
    t = 1.0;

  dx = dx * t + x1 - x;
  dy = dy * t + y1 - y;

You should have someone double check me on this, though, just to be sure :slight_smile:

Radim

--
Radim Blazek
Praga, Bohemia

3

On Feb 21, Radim Blazek wrote:
> On Thursday 21 February 2002 13:45, Markus Neteler wrote:
> > Hi developers,
> > I was looking through the grass src code and found a bug in the
> > G_distance2_point_to_line() function in libes/gis/line_dist.c
> >
> > lines 47,48 should be changed from
> >
> > dx = dx * t + x1;
> > dy = dy * t + y1;
> > to
>
> I think that (dx * t + x1, dy * t + y1) is point x,y projected to line
> x1,y1->x2,y2 (comment in code: /* go t from x1,y1 towards x2,y2 */)
> and we want distance between this point and and x,y. That is:
> dx = dx * t + x1 - x;
> dy = dy * t + y1 - y;
>
> You should have someone double check me on this, though, just to be sure :slight_smile:

I've double-checked this. Radim's fix above is correct.

You might even write it in the following (equivalent) way to make a
more clear parallel with the 3 cases above it in the code:

dx = x - (dx * t + x1);
dy = y - (dy * t + y1);

(the change of sign is not significant to the final result since both
dx and dy are squared at the end).

Cheers,

-Carl

--
Carl Worth
USC Information Sciences Institute cworth@east.isi.edu
3811 N. Fairfax Dr. #200, Arlington VA 22203 703-812-3725