This is too long for a comment, but it is not the answer to your issues.
Unfortunately, there is no avoiding that bottleneck.
You could theoretically set limits to the calculation ranges by point numbers (Line & Station, or Track & Bin) if you know the spacing, and make general assumptions with regards to the distance from a receiver point to a source point based upon station interval, and line spacing
You obviously know the following, but I will expound on it anyway.
You are going to have two loops, or do while statements.
One will fetch each receiver point (R File), storing the coordinate values, and other pertinent information.
The second loop will fetch each source point (S File), calculate the inverse or ray distance, and also take the x, and y coordinates of both points, sum the x values and the y values then divide each sum by 2 to derive the mid point coordinates.
I would also take the inverse or ray distance, and divide this by two, this will give you the depth of the record, and I would calculate the azimuth of the ray.
Then write all of these values to another file. (The relations or X file).
So in the third file I would write the Receiver Point Number, Source Point Number, MidpointX, MidpointY, Radial_distance, Half_Distance, Azimuth.
You could add other info if you needed to such as Receiver Coordinates, and Source Coordinates, Preplot or Postplot Designators and so on.
Then after each source point has been looped through, you return to the receiver point file, fetch the next one (first loop), then repeat the process with the source points (second loop) until all of the points in both files are looped through.
It is very calculation intensive, but the math is simple for the most part.
Optimization of the loops is one of the most important things here.
On a seismic grid with 8000 receiver points, and 8000 source points, you will have 64,000,000 mid point locations, and in all honesty, that is not a very intense grid compared to some I have worked on.
The rest of this is not meant to discourage you, but I see the road you are taking, and it will progress to the following.
There is a lot more fun to be had when you start to analyze azimuthal distribution, and create spider graphs for each bin. Then querying specific depth ranges for fold.
Other fun will be the blending of surveyed locations with unsurveyed locations for design purposes and to avoid the striping of the data.
It can then turn into something even more mind boggling when you start turning off specific points within the grid, and have to regenerate the counts for the mid point locations within specific bins based upon the loss of these points.
Have as much fun with this as you can.
Take small steps, and work on optimizing the loops first.
So in a follow up to this, out of morbid curiosity, insanity, and stubbornness, I decided to try this using only the core tools, and processing tools contained within QGIS 3.4.3.
I will not go into full step by step details, but I started by using the regular points generator, and built a grid of points on a 220 ft. by 220 ft. grid for both receiver points, and source points. I then rotated both grids using the processing toolbox rotate function.
After that I built a polygon grid of 110 ft by 110 ft, and then translated and rotated the grid so the points were centered within the polygons.
I then queried out the desired points from the receiver and source point files.
Then I ran the Distance Matrix utility from the Vector Menu (it is also part of the Processing Toolbox)
This created a record for each point on the Receiver Point Grid related to each point on the Source Point Grid.
In all there are 3,128,160 combinations.
The distance matrix file gives you an input and target point id, and the distance between the two points, which in many cases was found to be erroneous unless I am misunderstanding what is happening there.
I then joined the receiver point file with the distance matrix file exporting the coordinates of the receiver points associated with the file.
I repeated the above process using, this time joining the source point file with exported file.
From there I added a mid point x and y field, as well as a mid point distance field, and updated those fields using the SQL functions.
This was a very time consuming process with frequent crashes along the way.
I attempted this using both shape files, and spatialite files. The spatialite files seemed to be more stable at least initially. (After each crash, I rebooted the computer.)
I took the file with the calculated midpoint data, and exported a CSV file of containing the Receiver Point ID, Source Point ID, Mid Point X, Mid Point Y, and the Mid Point Distance.
I then opened Mid Point CSV file, and exported it as a Spatialite file.
The point count was correct, however there seems to be an error in the coordinates as some of the points fall outside of the boundary of the Bin Grid.
I have not tracked this error down as of yet by comparing the joined receiver and source grid coordinates with the original coordinates, that is next on my list.
Out of curiosity I decided to calculate a full fold map with what I came up with. I knew it would be wrong, but I wanted to see how long it would take counting the points in polygons. (Using the first run shape files)
It took a little over eight hours using the Points in Polygons routine.
I am not going to run it again using the spatialite files until I solve the rotational issues which is apparent in both the shape files of the first run, and the spatialite files of the second run.
I am still having problems wrapping my head around using Python and looping through both the receiver and source point files staying with a single point in the receiver file while calculating mid points and distances to each source point from that receiver point, then moving to the next receiver point and repeating this process until all points are calculated.
The rest of the math is easy, but the looping thing has given me fits.
All I want to write out to a CSV file is the Receiver Point ID, Source Point ID, Mid Point X & Y coordinates, and the Distance or Depth. That is the easy part for me so far.
Given what I have experienced with the above experiment, I hope you have a truly monstrous computer, huge hard drives, and a ridiculous amount of RAM, along with an inordinate amount of patience.