Fixed merge conflicts. (bbb37b3c) · Commits · aflab / astrogeology / Autocnet

autocnet/camera/camera.py

+41 −31

Original line number	Diff line number	Diff line
		import numpy as np
		from autocnet.camera.utils import crossform
		try:
		import cv2

		except:
		cv2 = None

		def compute_epipoles(f):
		"""
		@@ -21,9 +24,7 @@ def compute_epipoles(f):
		"""
		u, _, _ = np.linalg.svd(f)
		e = u[:, -1]
		e1 = np.array([[0, -e[2], e[1]],
		[e[2], 0, -e[0]],
		[-e[1], e[0], 0]])
		e1 = crossform(e)

		return e, e1

		@@ -102,24 +103,37 @@ def triangulate(pt, pt1, p, p1):
		pt = pt.T
		if pt1.shape[0] != 3:
		pt1 = pt1.T

		#if cv2:
		X = cv2.triangulatePoints(p, p1, pt[:2], pt1[:2])

		# Homogenize
		X /= X[3]

		X /= X[3] # Homogenize
		return X


		"""
		# Stubbed in for a ticket addressing making OpenCV an optional dependency
		else:
		npts = len(pt)
		a = np.zeros((4, 4))
		coords = np.empty((npts, 4))
		coords[:] = 1
		for i in range(npts):
		# Compute AX = 0
		a[0] = pt[i][0] * p[2] - p[0]
		a[1] = pt[i][1] * p[2] - p[1]
		a[2] = pt1[i][0] * p1[2] - p1[0]
		a[3] = pt1[i][1] * p1[2] - p1[1]
		# v.T is a least squares solution that minimizes the error residual
		u, s, vh = np.linalg.svd(a)
		v = vh.T
		coords[i] = v[:,3] / (v[:,3][-1])
		return coords.T
		"""
		def projection_error(p1, p, pt, pt1):
		"""
		Based on Hartley and Zisserman p.285 this function triangulates
		image correspondences and computes the reprojection error
		by back-projecting the points into the image.

		References
		----------
		.. [Hartley2003]
		This is the classic cost function (minimization problem) into
		the gold standard method for fundamental matrix estimation.

		Parameters
		-----------
		@@ -137,29 +151,25 @@ def projection_error(p1, p, pt, pt1):

		Returns
		-------
		residuals : ndarray
		(n, 1) residuals for each correspondence

		cumulative_error : float
		sum of the residuals
		reproj_error : ndarray
		(n, 1) vector of reprojection errors


		"""
		# SciPy least squares solver needs a vector, so reshape back to a 3x4 c
		# camera matrix at each iteration

		if p1.shape != (3,4):
		p1 = p1.reshape(3,4)

		# Triangulate the correspondences
		xw_est = triangulate(pt, pt1, p, p1)

		# Back project and homogenize
		xhat = np.dot(p, xw_est)
		xhat /= xhat[2]
		x2hat = np.dot(p1, xw_est)
		x2hat /= x2hat[2]
		xhat = triangulate(pt, pt1, p, p1)
		xhat1 = xhat[:3] / xhat[2]
		xhat2 = p1.dot(xhat)
		xhat2 /= xhat2[2]

		# Compute residuals
		dist = (pt.T - xhat)2 + (pt1.T - x2hat)2
		residuals = np.sum(dist, axis=0)
		reproj_error = np.sum(dist)
		# Compute error
		cost = (pt - xhat1)2 + (pt1 - xhat2)2
		cost = np.sqrt(np.sum(cost, axis=0))

		return residuals, reproj_error
		return cost

autocnet/camera/tests/test_camera.py

+2 −3

Original line number	Diff line number	Diff line
		@@ -60,7 +60,6 @@ class TestCamera(unittest.TestCase):
		c = camera.triangulate(coords1, coords2, p, p1)
		np.testing.assert_array_almost_equal(c, truth)

		truth = np.array([ 3.09866357e-02, 2.60295132e-01,
		8.12871690e-02, 5.57281224e-01, 4.72226586e-04])
		residuals, reproj_error = camera.projection_error(p1, p, coords1, coords2)
		truth = np.array([0.17603 , 0.510191, 0.285109, 0.746513, 0.021731])
		residuals = camera.projection_error(p1, p, coords1.T, coords2.T)
		np.testing.assert_array_almost_equal(residuals, truth)

autocnet/camera/utils.py

+8 −0

Original line number	Diff line number	Diff line
		import math
		import numpy as np

		def crossform(a):
		"""
		Convert a three element vector into a 3 x 3 skew matrix as per
		Hartley and Zisserman pg. 581
		"""
		return np.array([[0, -a[2], a[1]],
		[a[2], 0, -a[0]],
		[-a[1], a[0], 0]])

		def normalize(a):
		"""

autocnet/graph/edge.py

+295 −9

Original line number	Diff line number	Diff line
		@@ -3,6 +3,8 @@ from collections import MutableMapping

		import numpy as np
		import pandas as pd

		from scipy.spatial.distance import cdist
		from scipy.spatial import Voronoi
		import cv2

		@@ -12,9 +14,9 @@ from autocnet.matcher import outlier_detector as od
		from autocnet.matcher import suppression_funcs as spf
		from autocnet.matcher import subpixel as sp
		from autocnet.matcher.feature import FlannMatcher
		from autocnet.transformation.decompose import coupled_decomposition
		from autocnet.transformation.transformations import FundamentalMatrix, Homography
		from autocnet.vis.graph_view import plot_edge
		from autocnet.vis.graph_view import plot_node
		from autocnet.vis.graph_view import plot_edge, plot_node, plot_edge_decomposition
		from autocnet.cg import cg


		@@ -93,32 +95,298 @@ class Edge(dict, MutableMapping):
		def health(self):
		return self._health.health

		def match(self, k=2):
		def decompose_and_match(self, k=2, maxiteration=3, size=18, buf_dist=3,**kwargs):
		"""
		Similar to match, this method first decomposed the image into
		$4^{maxiteration}$ subimages and applys matching between each sub-image.

		This method is potential slower than the standard match due to the
		overhead in matching, but can be significantly more accurate. The
		increase in accuracy is a function of the total image size. Suggested
		values for maxiteration are provided below.

		Parameters
		----------
		k : int
		The number of neighbors to find

		method : {'coupled', 'whole'}
		whether to utilize coupled decomposition
		or match the whole image

		maxiteration : int
		When using coupled decomposition, the number of recursive
		divisions to apply. The total number of resultant
		sub-images will be 4 ** maxiteration. Approximate values:

		\| Number of megapixels \| maxiteration \|
		\|----------------------\|--------------\|
		\| m < 10 \|1-2\|
		\| 10 < m < 30 \| 3 \|
		\| 30 < m < 100 \| 4 \|
		\| 100 < m < 1000 \| 5 \|
		\| m > 1000 \| 6 \|

		size : int
		When using coupled decomposition, the total number of points
		to check in each sub-image to try and find a match.
		Selection of this number is a balance between seeking a
		representative mid-point and computational cost.

		buf_dist : int
		When using coupled decomposition, the distance from the edge of
		the (sub)image a point must be in order to be used as a
		partioning point. The smaller the distance, the more likely
		percision errors can results in erroneous partitions.
		"""
		def mono_matches(a, b, aidx=None, bidx=None):
		"""
		Apply the FLANN match_features

		Parameters
		----------
		a : object
		A node object

		b : object
		A node object

		aidx : iterable
		An index for the descriptors to subset

		bidx : iterable
		An index for the descriptors to subset
		"""
		# Subset if requested
		if aidx is not None:
		ad = a.descriptors[aidx]
		else:
		ad = a.descriptors

		if bidx is not None:
		bd = b.descriptors[bidx]
		else:
		bd = b.descriptors

		# Load, train, and match
		fl.add(ad, a.node_id, index=aidx)
		fl.train()
		matches = fl.query(bd, b.node_id, k, index=bidx)
		self._add_matches(matches)
		fl.clear()

		def func(group):
		ratio = 0.8
		res = [False] * len(group)
		if len(res) == 1:
		return [single]
		if group.iloc[0] < group.iloc[1] * ratio:
		res[0] = True
		return res

		# Grab the original image arrays
		sdata = self.source.get_array()
		ddata = self.destination.get_array()

		ssize = sdata.shape
		dsize = ddata.shape

		# Grab all the available candidate keypoints
		skp = self.source.get_keypoints()
		dkp = self.destination.get_keypoints()

		# Set up the membership arrays
		self.smembership = np.zeros(sdata.shape, dtype=np.int16)
		self.dmembership = np.zeros(ddata.shape, dtype=np.int16)
		self.smembership[:] = -1
		self.dmembership[:] = -1
		pcounter = 0

		# FLANN Matcher
		fl= FlannMatcher()

		for k in range(maxiteration):
		partitions = np.unique(self.smembership)
		for p in partitions:
		sy_part, sx_part = np.where(self.smembership == p)
		dy_part, dx_part = np.where(self.dmembership == p)

		# Get the source extent
		minsy = np.min(sy_part)
		maxsy = np.max(sy_part) + 1
		minsx = np.min(sx_part)
		maxsx = np.max(sx_part) + 1

		# Get the destination extent
		mindy = np.min(dy_part)
		maxdy = np.max(dy_part) + 1
		mindx = np.min(dx_part)
		maxdx = np.max(dx_part) + 1

		# Clip the sub image from the full images
		asub = sdata[minsy:maxsy, minsx:maxsx]
		bsub = ddata[mindy:maxdy, mindx:maxdx]

		# Utilize the FLANN matcher to find a match to approximate a center
		fl.add(self.destination.descriptors, self.destination.node_id)
		fl.train()

		scounter = 0
		decompose = False
		while True:
		sub_skp = skp.query('x >= {} and x <= {} and y >= {} and y <= {}'.format(minsx, maxsx, minsy, maxsy))
		# Check the size to ensure a valid return
		if len(sub_skp) == 0:
		break # No valid keypoints in this (sub)image
		if size > len(sub_skp):
		size = len(sub_skp)
		candidate_idx = np.random.choice(sub_skp.index, size=size, replace=False)
		candidates = self.source.descriptors[candidate_idx]
		matches = fl.query(candidates, self.source.node_id, k=3, index=candidate_idx)

		# Apply Lowe's ratio test to try to find a 'good' starting point
		mask = matches.groupby('source_idx')['distance'].transform(func).astype('bool')
		candidate_matches = matches[mask]
		match_idx = candidate_matches['source_idx']

		# Extract those matches that pass the ratio check
		sub_skp = skp.iloc[match_idx]

		# Check that valid points remain
		if len(sub_skp) == 0:
		break

		# Locate the candidate closest to the middle of all of the matches
		smx, smy = sub_skp[['x', 'y']].mean()
		mid = np.array([[smx, smy]])
		dists = cdist(mid, sub_skp[['x', 'y']])
		closest = sub_skp.iloc[np.argmin(dists)]
		closest_idx = closest.name
		soriginx, soriginy = closest[['x', 'y']]

		# Grab the corresponding point in the destination
		q = candidate_matches.query('source_idx == {}'.format(closest.name))
		dest_idx = q['destination_idx'].iat[0]
		doriginx = dkp.at[dest_idx, 'x']
		doriginy = dkp.at[dest_idx, 'y']

		if mindy + buf_dist <= doriginy <= maxdy - buf_dist\
		and mindx + 3 <= doriginx <= maxdx - 3:
		# Point is good to split on
		decompose = True
		break
		else:
		scounter += 1
		if scounter >= maxiteration:
		break

		# Clear the Flann matcher for reuse
		fl.clear()

		# Check that the identified match falls within the (sub)image
		# This catches most bad matches that have passed the ratio check
		if not (buf_dist <= doriginx - mindx <= bsub.shape[1] - buf_dist) or not\
		(buf_dist <= doriginy - mindy <= bsub.shape[0] - buf_dist):
		decompose = False

		if decompose:
		# Apply coupled decomposition, shifting the origin to the sub-image
		s_submembership, d_submembership = coupled_decomposition(asub, bsub,
		sorigin=(soriginx - minsx, soriginy - minsy),
		dorigin=(doriginx - mindx, doriginy - mindy),
		**kwargs)

		# Shift the returned membership counters to a set of unique numbers
		s_submembership += pcounter
		d_submembership += pcounter

		# And assign membership
		self.smembership[minsy:maxsy,
		minsx:maxsx] = s_submembership
		self.dmembership[mindy:maxdy,
		mindx:maxdx] = d_submembership
		pcounter += 4

		# Now match the decomposed segments to one another
		for p in np.unique(self.smembership):
		sy_part, sx_part = np.where(self.smembership == p)
		dy_part, dx_part = np.where(self.dmembership == p)

		# Get the source extent
		minsy = np.min(sy_part)
		maxsy = np.max(sy_part) + 1
		minsx = np.min(sx_part)
		maxsx = np.max(sx_part) + 1

		# Get the destination extent
		mindy = np.min(dy_part)
		maxdy = np.max(dy_part) + 1
		mindx = np.min(dx_part)
		maxdx = np.max(dx_part) + 1

		# Get the indices of the candidate keypoints within those regions / variables are pulled before decomp.
		sidx = skp.query('x >= {} and x <= {} and y >= {} and y <= {}'.format(minsx, maxsx, minsy, maxsy)).index
		didx = dkp.query('x >= {} and x <= {} and y >= {} and y <= {}'.format(mindx, maxdx, mindy, maxdy)).index
		# If the candidates < k, OpenCV throws an error
		if len(sidx) >= k and len(didx) >=k:
		mono_matches(self.source, self.destination, sidx, didx)
		mono_matches(self.destination, self.source, didx, sidx)

		def match(self, k=2, **kwargs):
		"""
		Given two sets of descriptors, utilize a FLANN (Approximate Nearest
		Neighbor KDTree) matcher to find the k nearest matches. Nearness is
		the euclidean distance between descriptors.

		The matches are then added as an attribute to the edge object.

		Parameters
		----------
		k : int
		The number of neighbors to find
		"""
		def mono_matches(a, b, aidx=None, bidx=None):
		"""
		Apply the FLANN match_features

		Returns
		-------
		Parameters
		----------
		a : object
		A node object

		b : object
		A node object

		aidx : iterable
		An index for the descriptors to subset

		bidx : iterable
		An index for the descriptors to subset
		"""
		def mono_matches(a, b):
		fl.add(a.descriptors, a.node_id)
		# Subset if requested
		if aidx is not None:
		ad = a.descriptors[aidx]
		else:
		ad = a.descriptors

		if bidx is not None:
		bd = b.descriptors[bidx]
		else:
		bd = b.descriptors

		# Load, train, and match
		fl.add(ad, a.node_id, index=aidx)
		fl.train()
		self._add_matches(fl.query(b.descriptors, b.node_id, k))
		matches = fl.query(bd, b.node_id, k, index=bidx)
		self._add_matches(matches)
		fl.clear()

		fl = FlannMatcher()
		mono_matches(self.source, self.destination)
		mono_matches(self.destination, self.source)



		def _add_matches(self, matches):
		"""
		Given a dataframe of matches, either append to an existing
		@@ -177,7 +445,7 @@ class Edge(dict, MutableMapping):
		See Also
		--------
		autocnet.transformation.transformations.FundamentalMatrix
		:

		"""
		if not hasattr(self, 'matches'):
		raise AttributeError('Matches have not been computed for this edge')
		@@ -208,6 +476,21 @@ class Edge(dict, MutableMapping):
		# Set the initial state of the fundamental mask in the masks
		self.masks = ('fundamental', mask)

		def refine_fundamental_matrix_matches(self, **kwargs): # pragma: no cover
		"""
		Given an estimated fundamental matrix, refine the correspondences based
		on the reprojective error.

		See Also
		--------
		autocnet.transformation.transformations.FundamentalMatrix.refine_matches
		"""
		if not hasattr(self, 'fundamental_matrix'):
		raise AttributeError('No fundamental matrix exists for this edge.')
		return

		self.fundamental_matrix.refine_matches(**kwargs)

		def compute_homography(self, method='ransac', clean_keys=[], pid=None, **kwargs):
		"""
		For each edge in the (sub) graph, compute the homography
		@@ -406,6 +689,9 @@ class Edge(dict, MutableMapping):
		# Else, plot the whole edge
		return plot_edge(self, ax=ax, clean_keys=clean_keys, **kwargs)

		def plot_decomposition(self, args, *kwargs): #pragma: no cover
		return plot_edge_decomposition(self, args, *kwargs)

		def clean(self, clean_keys, pid=None):
		"""
		Given a list of clean keys and a provenance id compute the

autocnet/graph/network.py

+21 −0

Original line number	Diff line number	Diff line
		@@ -296,6 +296,17 @@ class CandidateGraph(nx.Graph):
		"""
		self.apply_func_to_edges('match', args, *kwargs)

		def decompose_and_match_features(self, args, *kwargs):
		"""
		For all edges in the graph, apply coupled decomposition followed by
		feature matching.

		See Also
		--------
		autocnet.graph.edge.Edge.decompose_and_match
		"""
		self.apply_func_to_edges('decompose_and_match', args, *kwargs)

		def compute_clusters(self, func=markov_cluster.mcl, args, *kwargs):
		"""
		Apply some graph clustering algorithm to compute a subset of the global
		@@ -402,6 +413,16 @@ class CandidateGraph(nx.Graph):
		'''
		self.apply_func_to_edges('compute_fundamental_matrix', args, *kwargs)

		def refine_fundamental_matrix_matches(self, args, *kwargs):
		"""
		Refine the fundamental matrix matches using reprojective error

		See Also
		--------
		autocnet.transformation.transformations.FundamentalMatrix.refine_matches
		"""
		self.apply_func_to_edges('refine_fundamental_matrix_matches', args, *kwargs)

		def subpixel_register(self, args, *kwargs):
		'''
		Compute subpixel offsets for all edges using identical parameters