missing pairs in pair selection

[originlake]

I noticed an issue when selecting pairs during feature matching, find_best_altitude returns very large value and hence causing missing pairs.
Here attached the dumped data from my dataset (I'm not allowed to shared the source data). One is origin, the other is directions, which is the exact input for find_best_altitude
https://drive.google.com/file/d/1dYMLx7NsPB8-hYLdarJkrwikSzOej2il/view?usp=sharing

All the images are nadir images, they are pointing down, some of the directions are exact (0,0,-1), while some of them are not but very closed to (0,0,-1), probably due to precision at compute time. The samples_x and samples_y plots like below, as you can see, in order to fit a 2nd degree equation and find the minimum/maximum point, you can get a very large or small extrema, negative extrema will be ignored, so we only consider very large value here. In the above example, I get 8136309.134196892. When multiply this to the directions that are not exact (0,0,-1), for example, (-1.56556042e-05, -6.45686431e-05, -1.00000000e+00), which will get (-1.27378835e+02, -5.25350441e+02, -8.13630913e+06). You can see there will introduce several hundreds of horizontal offsets, this will make it not able to be paired with images that it suppose to be. I got 2 sub reconstructions that are not aligned with each other.

OpenSfM/opensfm/pairs_selection.py

Lines 74 to 102 in c798da1

	def find_best_altitude(
	origin: Dict[str, np.ndarray], directions: Dict[str, np.ndarray]
	) -> float:
	"""Find the altitude that minimize X/Y bounding box. Domain is [0, MAXIMUM_Z].
	'origin' contains per-image positions in worl coordinates
	'directions' viewing directions expected to be homogenized with z=1.
	We sample it every SAMPLE_Z and regress a second polynomial from which we take its extrema.
	"""
	directions_base = np.array([p for p in directions.values()])
	origin_base = np.array([p for p in origin.values()])
	samples_x, samples_y = [], []
	for current_z in range(1, MAXIMUM_Z, SAMPLE_Z):
	scaled = origin_base + directions_base / DEFAULT_Z * current_z
	current_size = (np.max(scaled[:, 0]) - np.min(scaled[:, 0])) ** 2 + (
	np.max(scaled[:, 1]) - np.min(scaled[:, 1])
	) ** 2
	samples_x.append(current_z)
	samples_y.append(current_size)
	coeffs = np.polyfit(samples_x, samples_y, 2)
	extrema = -coeffs[1] / (2 * coeffs[0])
	if extrema < 0:
	logger.info(
	f"Altitude is negative ({extrema}) : viewing directions are probably divergent. Using default altitude of {DEFAULT_Z}"
	)
	extrema = DEFAULT_Z
	return extrema