How can I align roads to follow hex bisection paths in QGIS?

I’ve created a hex grid and I would like to align and simplify each road such that each road segment has basically 2-3 vertices within each hex:

1. Entry side (centered on the side it enters on)
2. Center of hex
3. Exit side (center on the side it exists on) — Unless the road terminates within that hex.

This answer provided some python I’ve adapted. I had to do some refactoring for version differences, also. It takes a pretty long time to execute on what I consider a well-specced machine, probably due to the size of my map, but that’s fine. My goal is to reduce my time manual doing things, so I don’t mind the execution time so much, though I wouldn’t mind learning tips for parallelization of foreach processing the way you can easily do in newer versions of powershell.

As of this edit of the question, those three goals are met, but what are left are two cases I’d like to "clean up", but I’m having a hard time conceptualizing solutions in ways that don’t unintentionally break parts I want unchanged.

Here are my current results. The orange lines are the original road, white is "snapped", and the circled areas are segments I don’t want to have (ignore the two vertices with tiny red circles, I had the layer in vertex edit mode). The ones in the purple circle can just be deleted, they represent a road that enters and leaves in the same side. The ones in red should be deleted and replaced with two lines directly connecting the hexes that were connected, instead of using the intermediate hex where the road passes through very near and tight to one vertex.

Here is the python I used. I changed it to check each line of each hex and if the road crosses it, create a segment from the centroid of the line to the centroid of the hex.

I’m not sure yet how to begin thinking about how to remove the elements I don’t want (red and purple circles), and how to put in replacement ones for the red circles.

hexgrid = QgsVectorLayer(r"C:UsersBrendanDesktopGenerated LayersCore RegionamericanSouthwestGridLC.shp", "hexgrid", "ogr")
roads = QgsVectorLayer(r"C:UsersBrendanDesktopGenerated LayersCore RegionroadsMultipart.shp", "roads", "ogr")  # Must be multipart!

roadFeat = next(roads.getFeatures()) # We just have 1 geometry
road = roadFeat.geometry() 
indicesHexSides = ((0,1), (1,2), (2,3), (3,4), (4,5), (5,0))

epsilon = 0.01
# Function to compare whether 2 segments are equal (even if inverted)
def isSegmentAlreadySaved(v1, v2):
    for segment in listSegments:        
        p1 = QgsPointXY(segment[0][0], segment[0][1])
        p2 = QgsPointXY(segment[1][0], segment[1][1])
        if v1.compare(p1, epsilon) and v2.compare(p2, epsilon) 
            or v1.compare(p2, epsilon) and v2.compare(p1, epsilon):
            return True
    return False

# Let's find the nearest sides of hexagons where routes cross
listSegments = []
for hexFeat in hexgrid.getFeatures():
    hex = hexFeat.geometry()
    if hex.intersects( road ):
        for side in indicesHexSides:
            line = QgsGeometry.fromPolyline([hex.vertexAt(side[0]), hex.vertexAt(side[1])])
            if line.intersects( road ):
                # Only append new lines, we don't want duplicates!!!
                if not isSegmentAlreadySaved(hex.centroid().asPoint(), line.centroid().asPoint()): 
                    listSegments.append( [[hex.centroid().asPoint().x(), hex.centroid().asPoint().y()], [line.centroid().asPoint().x(),line.centroid().asPoint().y()]] )