Mid line between two arbitrary lines defined by points (ie - midline of a racetrack)

Question

Suppose I have two closed loops representing either boundary of a racetrack in two dimensional x,y space. The start of each loop if connected by a line would represent the start/finish line, the midpoint of this line would represent the centerline of the track.

Moving outward from the start line each loop is comprised of line segments which have random lengths with relatively small variations. The number of line segments is high at approximately 20-40k for each loop.

With this in mind the problem I am trying to solve is finding a closed loop representing the midline of the two other loops (ie - the centerline of the track).

At first I thought given that each small line segment is relatively close in length that I could take a step one segment in length which is the smaller of the two segments of either loop, this shorter length is used to inform a point along the longer segment and then a midpoint is found between the two steps of equal length. The new position along each line is updated and the process repeats effectively ping-ponging back and forth between the active loop having the shorter step forward.

I am doing this in python where I have a csv file containing the x and y position (x = *.PositionX and y = *.PositionZ) for the two loops. For reference below is what I was trying to use

import traceback
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt    

fDict = {"Bernese Alps":{"fName":"2021-09-12202910_FM7_DASH","left":4,"right":6},
         "Maple Valley":{"fName":"2021-08-30171416_FM7_DASH","left":4,"right":1},}

track = "Maple Valley"

df = pd.read_csv(fDict[track]["fName"]+".csv")

#The laps which correspond to track limits right and left
lapRight = fDict[track]["right"]
lapLeft  = fDict[track]["left"]

#Position from history for right and left
rightPos = df.loc[df.LapNumber==lapRight,['PositionX','PositionZ']]
leftPos  = df.loc[df.LapNumber==lapLeft,['PositionX','PositionZ']]

#Vector difference from one position to the next both right and left
rightDiff = rightPos.diff()
rightDiff.iloc[0] = rightPos.iloc[0]-rightPos.iloc[-1]
leftDiff = leftPos.diff()
leftDiff.iloc[0] = leftPos.iloc[0]-leftPos.iloc[-1]

#Vector difference magnitudes one position to the next both right and left
rightDiffMag = (rightDiff.PositionX.copy()*0+np.linalg.norm(rightDiff,axis=1)).rename({'PositionX':'PositionMag'})
leftDiffMag  = (leftDiff.PositionX.copy()*0+np.linalg.norm(leftDiff,axis=1)).rename({'PositionX':'PositionMag'})

#Normalized vectors to project onto
rightDiffNormed = rightDiff.div(rightDiffMag,'index')
leftDiffNormed  = leftDiff.div(leftDiffMag,'index')

#Methods to find center point
#(originally avg of x and y was used but change to confirm I was not doing that wrong)
def right_to_left(rightX,rightZ,leftX,leftZ):
    rtl_diff = np.array([leftX-rightX,leftZ-rightZ])
    rtl_mag = np.sqrt(sum(rtl_diff**2))
    rtl_dir = rtl_diff/rtl_mag
    return rtl_dir,rtl_mag

def midPoint(rightX,rightZ,leftX,leftZ):
    rtl_dir,rtl_mag = right_to_left(rightX,rightZ,leftX,leftZ)
    mid_vec = rtl_dir*rtl_mag*.5
    return rightX+mid_vec[0],rightZ+mid_vec[1]

#Create new right,left, and mid line containers
rights = [[rightPos.iloc[0].values[0]],[rightPos.iloc[0].values[1]]]
lefts  = [[leftPos.iloc[0].values[0]],[leftPos.iloc[0].values[1]]]
midP = midPoint(rights[0][0],rights[1][0],lefts[0][0],lefts[1][0])
mids = [[midP[0]],[midP[1]]]

#Initialize which size is longer

#Initialize partial remainders from prior step
partials = {'left':0.0,'right':0.0}
#Initialize indicie tracking for each side
indicies = {'left':1,'right':1}

print(len(leftPos),len(rightPos))
key = input("PAUSED - 'Y' TO CONTINUE 'N' TO ABORT: ")
if key.lower() ==  'y':
    Done = False
    try:
        while not Done:
            rightIdx = rightPos.index[indicies['right']]
            leftIdx  = leftPos.index[indicies['left']]

            #Project forward one point
            rightIdx_fwd = rightPos.index[0] if indicies['right']>=len(rightPos)-1 else rightPos.index[indicies['right']+1]
            leftIdx_fwd = rightPos.index[0] if indicies['left']>=len(leftPos)-1 else leftPos.index[indicies['left']+1]



            #If remaining left side is longer
            if leftDiffMag.loc[leftIdx]-partials['left'] >= rightDiffMag.loc[rightIdx]-partials['right']:
                active_side = 'left'
                #Step magnitude is full right step
                mag = rightDiffMag.loc[rightIdx]-partials['right']
                if not mag+partials['left'] >= leftDiffMag.loc[leftIdx_fwd]:
                    partials['left']+= mag
                else:
                    partials['left']=0.0
                    indicies['left']+=1

                partials['right']=0.0
                indicies['right']+=1
            else:
                #Step magnitude is full right step
                active_side = 'right'
                mag = leftDiffMag.loc[leftIdx]-partials['left']

                if not mag+partials['right'] >= rightDiffMag.loc[rightIdx_fwd]:
                    partials['right']+= mag
                else:
                    partials['right']=0.0
                    indicies['right']+=1

                partials['left']=0.0
                indicies['left']+=1

            rightUpdate = rightPos.loc[rightIdx]+mag*rightDiffNormed.loc[rightIdx_fwd]
            leftUpdate = leftPos.loc[leftIdx]+mag*leftDiffNormed.loc[leftIdx_fwd]
            midP = midPoint(*rightUpdate.values,*leftUpdate.values)

            rights[0].append(rightUpdate.PositionX)
            rights[1].append(rightUpdate.PositionZ)
            lefts[0].append(leftUpdate.PositionX)
            lefts[1].append(leftUpdate.PositionZ)
            mids[0].append(midP[0])
            mids[1].append(midP[1])

            #print(active_side,indicies['left'],leftIdx,indicies['right'],rightIdx,mag)

            if indicies['right']>=len(rightPos)-1 and indicies['left']>=len(leftPos)-1:
                Done=True

    except Exception as err:
        traceback.print_exception(None, err, err.__traceback__)

        print(len(rights[0]),len(lefts[0]),len(mids[0]))


        axi = plt.subplot()
        axi.plot(*rights,ls='-',marker='.',c='r')
        axi.plot(*lefts,ls='-',marker='.',c='b')
        axi.plot(*mids,ls='-',marker='.',c='k')
        axi.set_aspect('equal')

        fig,axes = plt.subplots(2)
        for i in range(2):
            axes[i].plot(rights[i],c='r')
            axes[i].plot(lefts[i],c='b')
            axes[i].plot(mids[i],c='k')

        #fig,axi = plt.subplots()
        #for rX,rZ,lX,lZ,mX,mZ in zip(*rights,*lefts,*mids):
            #axi.plot([rX,mX,lX],[rZ,mZ,lZ],c='k')
    finally:
        plt.show()

else:
    print("Processing was ABORTED!")

This worked better than I expected but:

• It is slow, and
• Depending on rate of curvature/difference in curvature (I am assuming that is what is driving the error) between the two loops over a given turn results in pulling the midline to one side or the other.
• Also, I am not sure if it is an implementation issue or simplify an artifact of the crude method but the update to the centerline is noisy/jagged (see below comment I believe this is an error in the implementation)

I am looking to find any method which will find the centerline. Bonus points if it lends itself well to obtaining a vector normal to the centerline (and its opposite -- rotated 180 degrees) and corresponding intersection of this normal with either of the bounding loops (ie - a line representing midline distance from start at any point along the track). Further bonus points if centerline step size is normalized.

I was thinking one solution, given that I know the direction of travel and center point at the start of both loops, would be to define a step size for the centerline, then use this to represent a half circle with radius equal the step size, then solve for the shortest line tangent to this circle which is bounded by intersection with the two outer loops. I have crudely illustrated this below. I am not sure there is an efficient/elegant way to solve using this method. ​

Thank you for any feedback/direction you can provide. I think there is probably a known solution to this type of problem and I am simply using the wrong search terms trying to find examples.

PS - Looking more closely and it appears there are aberrations in the update to the outer loops which likely helps explain the jagged update of the midline. Regardless I think the skew of the centerline to one side or the other remains a shortcoming of this approach and not an implementation problem.

Update: It would appear this is frequently a problem in the GIS domain where polylines are used to represent the boundaries of for instance a river and the rivers midline is desired.

Apparently one solution is to generate a Voronoi diagram. Doing this will naturally produce poly vertices and edges representing the midline.

Example: https://gis.stackexchange.com/questions/269238/create-a-center-line-between-2-polyline-qgis

I did a quick test of this following the documentation for scipy.spatial.voronoi_plot_2d (https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.voronoi_plot_2d.html) and it worked well even with the 47k+ points.

Now I just need to detect those Voronoi polygon verticies which are inside the two outer boundaries.

Additional word of caution; As you might expect this method has a hard time with the self intersecting track and so I will probably have to detect this condition or manually prevent it by cutting the track into multiple parts then sewing them back together.

Also some artifacts were produced at the start/finish line.

Answer 1

Ultimately a combination of scipy Voronoi and matplotlib path.contains_points is able to find the midline. I have not addressed the issue of self intersection, and there remains some minor issues with non-midline points at the start-finish line but those issues are fairly minor.

Below is full code to working solution:

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.path as mpltPath
from scipy.spatial import Voronoi
import time

fDict = {"Bernese Alps":{"fName":"2021-09-12202910_FM7_DASH","left":4,"right":6,"correctStart":-1},
         "Maple Valley":{"fName":"2021-08-30171416_FM7_DASH","left":4,"right":1,"correctStart":0},}

track = "Bernese Alps"

df = pd.read_csv(fDict[track]["fName"]+".csv")

# The laps which correspond to track limits right and left
lapRight = fDict[track]["right"]
lapLeft  = fDict[track]["left"]

# Find the lap first point distance offset from actual start
# (below Will break if either lap is 0)
offsetRight = (df.loc[df.LapNumber==lapRight,['DistanceTraveled']]-df.loc[df.LapNumber==lapRight-1,['DistanceTraveled']].max()).iloc[0].values[0]
offsetLeft = (df.loc[df.LapNumber==lapLeft,['DistanceTraveled']]-df.loc[df.LapNumber==lapLeft-1,['DistanceTraveled']].max()).iloc[0].values[0]

# Position history for right and left edges of track
rightPos = df.loc[df.LapNumber==lapRight,['PositionX','PositionZ']].values
leftPos  = df.loc[df.LapNumber==lapLeft,['PositionX','PositionZ']].values

# Finish line crossing vector to determine precise start finish line from distance traveled
finalVectRight = rightPos[-1]-rightPos[0]
finalVectLeft = leftPos[-1]-leftPos[0]

# Find actual start points assuming start and end are inline
# any alignment issues between start and end of lap will cause a skew in the startline
rightStart = rightPos[0]+finalVectRight/np.sqrt(sum(finalVectRight**2))*offsetRight
leftStart = leftPos[0]+finalVectLeft/np.sqrt(sum(finalVectLeft**2))*offsetLeft

if fDict[track] in (0,1):
    # Correct startline skey by assuming start line is normal either x or z axis
    # (this will not work where start line actually is non-normal the principal axes)
    nonStartDirIdx = fDict[track]["correctStart"]
    rightStart[nonStartDirIdx] =  rightPos[0][nonStartDirIdx]
    leftStart[nonStartDirIdx] =  leftPos[0][nonStartDirIdx]
    #average the starline error (appears to be on the order of ~ 2-3 inches ) in direction of travel
    startDirIdx = 1 if not nonStartDirIdx else 0
    avgStartDirPos = (rightStart[startDirIdx]+leftStart[startDirIdx])*.5
    rightStart[startDirIdx] =  avgStartDirPos
    leftStart[startDirIdx] =  avgStartDirPos

# Plot to check
fig,axi = plt.subplots()
plt.plot(*rightPos.T,'r-')
plt.plot(*leftPos.T,'b-')
plt.plot(*np.array([rightStart,leftStart]).T,'k-',marker='.')
axi.set_aspect('equal')
plt.show()

# Generate poly points to compute Voronoi and find midline
# concatenate the left and right side points to create open polygon
polygon = np.concatenate([rightPos,leftPos])
# turn this into matplotlib path object
path = mpltPath.Path(polygon)
# then compute Voronoi diagram for the polygon to obtain the Vornoi verticies
vor = Voronoi(polygon)

# Determine which of these lie inside the polygon
# First limit the points under consideration to those within the original physical limits of the track edge extents
# not doing this can cause the inside check to fail
vX = vor.vertices.T[0]
vZ = vor.vertices.T[1]
vorMask = (vX>=polygon.T[0].min())&(vX<=polygon.T[0].max())&(vZ>=polygon.T[1].min())&(vZ<=polygon.T[1].max())
verts = vor.vertices[vorMask]

print("Starting check to find Voronoi vertices inside polygon...")
start=time.time()
insideMask = path.contains_points(verts)
end=time.time()
print("\tChecking which of {:,d} points lie inside the polygon of {:,d} points took {:,.2f} seconds to complete.".format(len(verts),len(polygon),end-start))
midLine = verts[insideMask]

# Plot to check
fig,axi = plt.subplots()
plt.plot(*rightPos.T,'r-')
plt.plot(*leftPos.T,'b-')
plt.plot(*np.array([rightStart,leftStart]).T,'k-',marker='.')
plt.scatter(*verts[insideMask].T,color='orange',marker='.',s=.1)
axi.set_aspect('equal')
plt.show()

And here is an example output:

Starting check to find Voronoi vertices inside polygon...

Checking which of 164,650 points lie inside the polygon of 94,514 points took 143.72 seconds to complete.

PS - Use of matplotlib path.contains_points to find the midline was driven primarily by ease of use. As indicated this method can take a minute or two to find the midline out of 150k+ verticies inside a polygon comprised of a further ~100k points. For polygons and points lists of this size there are likely more efficient solutions but this one meets my needs and was easy to use. If you are interested in something faster see: What's the fastest way of checking if a point is inside a polygon in python

Answer 2

Assume for the moment that at a specific position the outside is two fairly long line segments. The center line would consist of points that have equal distance to their orthogonal projection obituary each of the lines. I would strongly suspect that if you have your lines expressed in Hesse normal form, i.e. as ax+by+c=0 with a²+b²=1, then just adding these two equations will yield the equation of the center line as the angular bisector resp. the corresponding generalisation of that for the case of parallel segments. Make sure the normals of the lines have matching orientation, e.g. all going left from the direction of a vehicle driving along the centerline.

Since you have more than one segment per side, you need to take changes between them into account. On either side you could have a convex or a concave vertex, i.e. a (slight) corner where the boundary bends either interests or outwards. If it bends outwards, then there is a point where the centerline has equal distance to both the segments that meet at that corner. At that point you so using one of them to define the centerline, and start using the other instead. A different perspective would be that each of the two consecutive segments defines its own centerline together with the opposite segment, and the point where the two centerlines intersect is the point where you switch from one to the other. You end up with line segments again.

If a corner bends inwards, then for some points on the centerline, the nearest point on that boundary would not be an orthogonal projection to a line segment but instead that corner. The set of all points that have equal distance to one point and one line is a parabola, with the point its focus and the line its directrix. The orthogonal projections would fall outside both segments on that side. The cut-over between line and segment would be where there orthogonal projections coincide with the end of each segment, i.e. on the line perpendicular to that segment through its endpoint. Are you prepared to deal with parabolic pieces of centerline in your output, or would you need to approximate those somehow?

All of this can become a lot more complicated of you have corners on both sides in close proximity. You might end up in situation where a parabolic piece of centerline gets interrupted by a change in the directrix segment, and stuff like that. Getting a complete solution that is guaranteed to do the right thing for all inputs sounds like a major effort. Something that will work reasonably well fort most typical use cases should be doable.

So far this has been a fairly mathematical answer. In other to make our more algorithmic, suppose that you are currently looking at a pair of opposite line segments and the centerline between them. Your code should determine the next corner on each of the sides, and the segment that comes after that corner. With that information it could determine the point where the current straight centerline would end die to that corner. Take the closer of the two changes, stitch together the centerline segments together in a suitable way (possibly using parabola pieces), then proceed with one of the segments exchanged for its successor. Watch out for centerline segments that have their end point before their starting point, that's where things become messy. I don't have a ready solution for dealing with them, but examine them on test data to decide how likely and how bad those are in practice.