Improving PyQGIS distance matrix efficiency

Question

I have written a custom processing tool, which I use within a large graphical model in QGIS. I am looking for a way to speed it up.

First step of this processing tool is to connect every point of one input layer to every point of another input layer and measure its line length. This takes quite a while, because usually both layers have more than 10000 points. But I am only interested in point pairs within a given distance to each other. So currently I am going to delete the lines afterwards, which are too long. I think this is quite ineffective, and therefore wanted to ask if there is an easy way to speed up this step, by using a spatial index or something like that. I could create one by using index = QgsSpatialIndex(layer.getFeatures()), but how would I proceed to only connect points I am interested in?

This is a minimal part of my current code:

Memorylayer_VL = QgsVectorLayer(str('LineString?crs=25832')), "tmp_lines", "memory")
Memorylayer_PR = Memorylayer_VL.dataProvider()
maxdistance = 150

# iterate over stop features and create lines
for stop_feat in stop_layer.getFeatures():
    point1 = QgsPoint(stop_feat.geometry().asPoint())
    for source_feat in possibility_layer.getFeatures():
        point2 = QgsPoint(source_feat.geometry().asPoint())
        temp_feat = QgsFeature(fields)
        temp_feat.setGeometry(QgsGeometry.fromPolyline([point1, point2]))
        temp_feat["line_length"] = temp_feat.geometry().length()    
        # Only add lines shorter than maxdistance and longer than 1m (Already simplified like this using @BERAs suggestion in comments) 
        # However, couldn't using a spatial index even have more effect?
        if (temp_feat.geometry().length() < maxdistance) and (tmp_line_feat.geometry().length() > 1):       
            Memorylayer_PR.addFeature(temp_feat)
            Memorylayer_VL.updateFields()
Memorylayer_VL.commitChanges()

# Doing more processing stuff on the just created lines layer....

Answer 1

Below the code for select all solution lines between 0 (excluded) and maxdistance (excluded too). All POSSIBILITY features are inserted inside a spatial index (QgsSpatialIndex). A filter, thanks to the spatial index, with a maxdistance square around the feature of the layer STOP is applied for speed the process and select (intersects) only POSSIBILITY features in the square.

If the layer STOP and POSSIBILITY are the same, the code below makes duplicates and isn't not very efficient because of the duplicate.

qgis_project = QgsProject.instance()
layers = {lyr.name(): lyr for lyr in qgis_project.mapLayers().values()}
lyr_from = layers["STOP"]
lyr_to = layers["POSSIBILITY"]
maxdistance = 150

spatial_idx_to = QgsSpatialIndex()
pr_from = lyr_from.dataProvider()
pr_to = lyr_to.dataProvider()
feats_from = pr_from.getFeatures()
feats_to = pr_to.getFeatures()

# insert destination features to the spatial index
spatial_idx_to.addFeatures(feats_to)

result_lines = []
# iterate over origin features
for ft_from in feats_from:
    # find all the nearest feature thanks to the spatial index
    geom_from = ft_from.geometry().asPoint()
    x_from = geom_from.x()
    y_from = geom_from.y()
    # create a rectangle at maxdistance of the from point
    rect = QgsRectangle(
        QgsPointXY(x_from - maxdistance, y_from - maxdistance),
        QgsPointXY(x_from + maxdistance, y_from + maxdistance)
    )
    # get all points in that intersects the rectangle
    distance_ids = spatial_idx_to.intersects(rect)
    for fid in distance_ids:
        geom_to = lyr_to.getFeature(fid).geometry().asPoint()
        # if the distance between origin and destination < maxdistance
        if 0 < geom_from.distance(geom_to) < maxdistance:
            line_feat = QgsFeature()
            # create a PolylineXY from origin and destination PointXY
            line_feat.setGeometry(QgsGeometry.fromPolylineXY([geom_from, geom_to]))
            # append the line feature to the result list
            result_lines.append(line_feat)

# create the memory layer
memlyr = QgsVectorLayer("LineString?crs=epsg:25832", "tmp_lines", "memory")
memlyr_pr = memlyr.dataProvider()
# add all resulting lines features
memlyr_pr.addFeatures(result_lines)
# finally add the line memory layer to the project
qgis_project.addMapLayer(memlyr)