You could potentially convert the DEM object into a Numpy array, and do your tests on that. For example, sorting by the height difference.
Unfortunately, the distance between each pixel to another pixel is something you have to do a lot of times, and as such, you can't really do anything to reduce the number of tests, aside from making a few rules based on difference in the indices between cells.
Having said that, what about this (in pseudo code):
- Get your raster as an array
- For each ROW and COL (this is your source)
- For each other ROW and COL (this is your destination)
- Ignore the destination if it's this cell
- Ignore the destination if the index is more than nMinX away in rows (find this by getting the your 500 m / cell width) and less than nMaxX
- Ignore the destination if the index is more than nMinY away in cols (similar to 5.) and nMaxY
- optional, but work out a diagonal ignore test if you need to?
- Check this cell is within the height constraint, and if it is, check the distance constraint.
- If checks in 8 pass, create your line and move along.
- You may want to store the line as a cell index, and also ignore the returning line from the destination cell (otherwise you'll have a from/to and to/from line for each pair.
I think doing this in Python will be fairly straight forward.