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
'=================================================================================