Closest Point on a Polyline

 

Paul Crisp (Paul.CRISP@syntegra.com) posted ro mapinfo-l on July 3, 2002 some VB code to find the closest point from a given point to a polyline (it could be the "foot" of the shortest perpendicular line from the point to every segment of the polyline, or an end point of the polyline if no perpendicular can be traced that would intersects a segment between its 2 nodes). As the translation to pure MapBasic code is not obvious, I have decided, with Paul's entire support,  to offer my solution for inspection and download.

The essence of Paul's code is contained in 3 functions: returnPointOfShortestPerpendicular(), FindClosestPointOnLineSegment() and PointIsBetween(). These functions need also a sub, ErrorHandler, that is not despite its name a standard MB procedure. It is useful to locate the source of an error in one of the functions.

I have added a small interface that shows a way to use these functions and that allows for some experimenting but is not a "production" tool. It is mainly made of "traps" to make sure that valid data is entered. It does not exclude working with "degrees" tables but the results are probably incorrect then.

The way I have used the SelChangedHandler is most primitive but it suffices for that kind of application. I would be most interested to learn of another way to code for successive selections of different objects.

You can download the following code as a MB listing in the ClstPntOnPline.zip file

'=================================================================================

'Demo program showing a way to find the closest point on a polyline to a given point

'

'using the VB code offered by Paul Crisp (Paul.CRISP@syntegra.com) July 3, 2002

'adapted to MapBasic and presented with a simple interface by

'Jacques Paris (jacques@paris-pc-gis.com) July 3, 2002

'

'=================================================================================

Include "mapbasic.def"

Type tPoint

X as float

Y as float

end type

Declare Function ReturnPointOfShortestPerpendicular(clsPoint As tPoint,

clsNodes() As tPoint, closestPoint As tPoint) As Logical

Declare Function FindClosestPointOnLineSegment(ptNearbyPoint As tPoint,

ptLineSegmentEnd1 As tPoint, ptLineSegmentEnd2 As tPoint,

ptClosestPoint As tPoint) As Logical

Declare Function PointIsBetween(pt1 As tPoint, pt2 As tPoint,

ptInBetween As tPoint) As Logical

Declare sub errorhandler(byval a as string)

Declare sub SelChangedHandler

Declare sub main

dim opline, opoint as object

dim clsPoint as tpoint

dim clsNodes() as tpoint

dim closestPoint as tpoint

dim my_table as string

dim tcall as smallint

 

'=================================================================================

sub main

'=================================================================================

if numwindows()=0 then

note"A window must be open"

exit sub

end if

if int(windowinfo(frontwindow(),3))<>1 then

note "Fontwindow must be a mapper"

exit sub

end if

note "Warning: this demo program will not return correct results if the "+chr$(10)+

"base table uses an 'earth, not projected' (in degrees) coordsys"

tcall=1

run menu command 304

tcall=2

note "Select a PolyLine"

end sub

 

'=================================================================================

sub SelChangedHandler

'=================================================================================

dim nnodes, i as integer

do case tcall

case 1

exit sub

case 2

'A polyline must be selected

opline=selection.obj

if int(objectinfo(opline,1))<>4 then

note "Selected object is not a Polyline"

exit sub

end if

if int(objectinfo(opline,21))>1 then

Note "Works only on single section polylines"

exit sub

end if

my_table=selectioninfo(1)

set coordsys table my_table

nnodes=objectinfo(opline,20)

redim clsNodes(nnodes)

for i=1 to nnodes

clsNodes(i).X=objectnodeX(opline,1,i)

clsNodes(i).Y=objectnodeY(opline,1,i)

next

tcall=3

run menu command 304

exit sub

case 3

tcall=4

note "Select a Point"

exit sub

case 4

'A point must be selected

opoint=selection.obj

if int(objectinfo(opoint,1))<>5 then

note "Selected object is not a Point"

exit sub

end if

if selectioninfo(1)<>my_table then

note "PolyLine and Point must be in the same table"

exit sub

end if

print selectioninfo(1)

clsPoint.X=centroidX(opoint)

clsPoint.Y=centroidY(opoint)

'Calling the minimum distance function

if not ReturnPointOfShortestPerpendicular(clsPoint, clsNodes(), closestPoint) then

note "Function failed to find closest point"

exit sub

end if

insert into my_table (obj) Values ( createpoint(closestPoint.X, closestPoint.Y) )

end case

end program

end sub

 

'=================================================================================

sub errorhandler (byval a as string)

'=================================================================================

note "Error encountered in "+a+chr$(10)+chr$(10)+err()+" "+error$()

end sub

 

'=================================================================================

Function ReturnPointOfShortestPerpendicular(clsPoint As tPoint,

clsNodes() As tPoint, closestPoint As tPoint) As Logical

'=================================================================================

' Purpose: Routine to return the end point of the perpendicular

' on to nearest line segment from a point

' Inputs: Point, points collection from line

'

' Returns: point

'*************************

OnError GoTo labelError

Dim i As Integer

Dim ptTemp As tPoint

Dim intCount As Integer

Dim iClosestPoint As Integer

Dim dblShortestDist As Float

Dim dblTempDist As Float

Dim ptLineSegmentEnd1, ptLineSegmentEnd2 as tPoint

intCount = ubound(clsNodes)

'Treat line feature line as a series of lines and get the perpendicular onto each

'take the shortest as the return value

'Initialise the values from the first line segment then compare each of the rest

'Give dblShortestDist a value just in case

closestPoint.X=clsNodes(1).X

closestPoint.Y=clsNodes(1).Y

dblShortestDist = Sqr((clsPoint.X - clsNodes(1).X) ^ 2 + (clsPoint.Y - clsNodes(1).Y) ^ 2)

iClosestPoint = 1

For i = 1 To intCount - 1

ptLineSegmentEnd1.X = clsNodes(i).X

ptLineSegmentEnd2.X = clsNodes(i + 1).X

ptLineSegmentEnd1.Y = clsNodes(i).Y

ptLineSegmentEnd2.Y = clsNodes(i + 1).Y

If FindClosestPointOnLineSegment(clsPoint, ptLineSegmentEnd1, ptLineSegmentEnd2, ptTemp) Then

dblTempDist = Sqr((clsPoint.X - ptTemp.X) ^ 2 + (clsPoint.Y - ptTemp.Y) ^ 2)

If dblTempDist < dblShortestDist Then

dblShortestDist = dblTempDist

closestPoint.X = ptTemp.X

closestPoint.Y = ptTemp.Y

iClosestPoint = i

End If

End If

Next

ReturnPointOfShortestPerpendicular = True

Exit Function

labelError:

ReturnPointOfShortestPerpendicular = False

Call ErrorHandler("GenericUtils: ReturnPointOfShortestPerpendicular")

End Function

 

'=================================================================================

Function FindClosestPointOnLineSegment(ptNearbyPoint As tPoint,

ptLineSegmentEnd1 As tPoint, ptLineSegmentEnd2 As tPoint,

ptClosestPoint As tPoint) As Logical

'=================================================================================

' Purpose: Generic Routine to take three points, two on a line and

' the other to one side, and return the closest point on the line

' i.e. give a point where the perpendicular from the third point will fall

' on the line

' Inputs: three points

'

' Returns: Logical - closest point

'*************************

Dim dblX As Float

Dim dblX2 As Float

Dim dblY As Float

Dim dblY2 As Float

Dim dblDeltaX As Float

Dim dblDeltaY As Float

Dim dblReturnX As Float

Dim dblReturnY As Float

Dim dblM As Float

Dim dblC As Float

Dim dblD As Float

Dim dblTempDist1 As Float

Dim dblTempDist2 As Float

 

OnError GoTo labelError

'If the point is close enough to either of the segment ends then return them

 

FindClosestPointOnLineSegment = False

If Abs(ptLineSegmentEnd1.X - ptNearbyPoint.X) < 0.1 And Abs(ptLineSegmentEnd1.Y - ptNearbyPoint.Y) < 0.1 Then

'its within 10cm - this will do

ptClosestPoint.X = ptLineSegmentEnd1.X

ptClosestPoint.Y = ptLineSegmentEnd1.Y

FindClosestPointOnLineSegment = True

Exit Function

End If

If Abs(ptLineSegmentEnd2.X - ptNearbyPoint.X) < 0.1 And Abs(ptLineSegmentEnd2.Y - ptNearbyPoint.Y) < 0.1 Then

'its within 10cm - this will do

ptClosestPoint.X = ptLineSegmentEnd2.X

ptClosestPoint.Y = ptLineSegmentEnd2.Y

FindClosestPointOnLineSegment = True

Exit Function

End If

'Otherwise we need to do some maths

 

'Get the numbers into the right order - mdblMapX2 needs to be bigger

If ptLineSegmentEnd1.X > ptLineSegmentEnd2.X Then

dblX = ptLineSegmentEnd2.X

dblX2 = ptLineSegmentEnd1.X

dblY = ptLineSegmentEnd2.Y

dblY2 = ptLineSegmentEnd1.Y

Else

dblX = ptLineSegmentEnd1.X

dblX2 = ptLineSegmentEnd2.X

dblY = ptLineSegmentEnd1.Y

dblY2 = ptLineSegmentEnd2.Y

End If

dblDeltaX = dblX2 - dblX

dblDeltaY = dblY2 - dblY

 

If dblDeltaX = 0 Then 'Vertical line

dblReturnX = dblX2

'Y value can be same as the existing point

dblReturnY = ptNearbyPoint.Y

 

ElseIf dblDeltaY = 0 Then 'horizontal line

dblReturnY = dblY2

'X value can be the same as the existing point

dblReturnX = ptNearbyPoint.X

 

Else

dblM = dblDeltaY / dblDeltaX

dblC = dblY - (dblM * dblX)

 

dblD = ptNearbyPoint.Y + (ptNearbyPoint.X / dblM)

 

'Now do the new coords

dblReturnX = (dblD - dblC) / (dblM + (1 / dblM))

dblReturnY = (dblM * dblReturnX) + dblC

End If

ptClosestPoint.X = dblReturnX

ptClosestPoint.y = dblReturnY

'Check this closest point is actually between the other two

If Not PointIsBetween(ptLineSegmentEnd1, ptLineSegmentEnd2, ptClosestPoint) Then

'The nearest point on that segment is actually one of the vertices

dblTempDist1 = sqr ((ptClosestPoint.X - ptLineSegmentEnd1.X) ^ 2 +(ptClosestPoint.Y - ptLineSegmentEnd1.Y) ^ 2)

dblTempDist2 = sqr ((ptClosestPoint.X - ptLineSegmentEnd2.X) ^ 2 +(ptClosestPoint.Y - ptLineSegmentEnd2.Y) ^ 2)

 

If dblTempDist1 < dblTempDist2 Then

ptClosestPoint.X = ptLineSegmentEnd1.X

ptClosestPoint.Y = ptLineSegmentEnd1.Y

Else

ptClosestPoint.X = ptLineSegmentEnd2.X

ptClosestPoint.Y = ptLineSegmentEnd2.Y

End If 'Vertex 1 is close

End If 'Returned point is not on line segment

FindClosestPointOnLineSegment = True

Exit Function

labelError:

FindClosestPointOnLineSegment = False

Call ErrorHandler("FindClosestPointOnLineSegment")

 

End Function

 

'=================================================================================

Function PointIsBetween(pt1 As tPoint, pt2 As tPoint,

ptInBetween As tPoint) As Logical

'=================================================================================

' Purpose: Routine to check whether a point falls between two others

' Inputs: Three Points

'

' Returns: Logical

'*************************

OnError GoTo labelError

Dim dblMapX As Float

Dim dblMapX2 As Float

Dim dblMapY As Float

Dim dblMapY2 As Float

 

PointIsBetween = False

'Sort the values so x2, y2 are the larger

If pt1.X > pt2.X Then

dblMapX = pt2.X

dblMapX2 = pt1.X

Else

dblMapX = pt1.X

dblMapX2 = pt2.X

End If

If pt1.Y > pt2.Y Then

dblMapY = pt2.Y

dblMapY2 = pt1.Y

Else

dblMapY = pt1.Y

dblMapY2 = pt2.Y

End If

'Use a 10cm tolerance and see if the point falls between the other two

If (ptInBetween.X + 0.1 >= dblMapX) And (ptInBetween.X - 0.1 <= dblMapX2) And

(ptInBetween.Y + 0.1 >= dblMapY) And (ptInBetween.Y - 0.1 <= dblMapY2) Then

PointIsBetween = True

End If

Exit Function

labelError:

PointIsBetween = False

Call ErrorHandler("GenericUtils: PointIsBetween")

End Function

 

'=================================================================================