Initial coupled decomp (f76bed03) · Commits · aflab / astrogeology / Autocnet

autocnet/graph/edge.py

+236 −4

Changes for autocnet/graph/edge.py: 236 added lines, 4 removed lines.
Original line number	Diff line number	Diff line
		@@ -3,6 +3,7 @@ from collections import MutableMapping

		import numpy as np
		import pandas as pd
		from scipy.spatial.distance import cdist

		from autocnet.utils import utils
		from autocnet.matcher import health
		@@ -10,6 +11,7 @@ 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.cg import cg
		@@ -90,7 +92,7 @@ class Edge(dict, MutableMapping):
		def health(self):
		return self._health.health

		def match(self, k=2):
		def match(self, k=2, method='coupled', maxiteration=3, size=18, **kwargs):
		"""
		Given two sets of descriptors, utilize a FLANN (Approximate Nearest
		Neighbor KDTree) matcher to find the k nearest matches. Nearness is
		@@ -102,20 +104,228 @@ class Edge(dict, MutableMapping):
		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.

		size : int
		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.

		Returns
		-------

		"""
		def mono_matches(a, b):
		fl.add(a.descriptors, a.node_id)

		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:
		bidx = 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()

		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

		if method == 'whole':
		fl = FlannMatcher()
		mono_matches(self.source, self.destination)
		mono_matches(self.destination, self.source)

		elif method == 'coupled':
		# Grab the matches data frame and identify the source and destination images and keypoints
		e = self

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

		ssize = sdata.shape
		dsize = ddata.shape

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

		smembership = np.zeros(sdata.shape, dtype=np.int16)
		dmembership = np.zeros(ddata.shape, dtype=np.int16)
		smembership[:] = -1
		dmembership[:] = -1
		maxiterations = 3
		pcounter = 0

		fl= FlannMatcher()

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

		"""
		Debug: Why is it that sometimes dy, dx is empty?
		"""

		# 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(e.destination.descriptors, e.destination.node_id)
		fl.train()

		searching = True
		scounter = 0
		while searching:
		sub_skp = skp.query('x >= {} and x <= {} and y >= {} and y <= {}'.format(minsx, maxsx, minsy, maxsy))
		size = 18
		if size > len(sub_skp):
		size = len(sub_skp)
		candidate_idx = np.random.choice(sub_skp.index, size=size, replace=False)
		candidates = e.source.descriptors[candidate_idx]
		matches = fl.query(candidates, e.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]

		### FLANN FINISHED ###

		# 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']])
		try:
		closest = sub_skp.iloc[np.argmin(dists)]
		except:
		continue
		closest_idx = closest.name
		soriginx, soriginy = closest[['x', 'y']]

		# Grab the corresponding point in the destination
		dest_idx = candidate_matches[candidate_matches['source_idx'] == closest.name]['destination_idx']
		doriginx, doriginy = dkp.loc[dest_idx][['x', 'y']].values[0]

		if not mindy + 1 <= doriginy <= maxdy - 1 or not mindx + 1 <= doriginx <= maxdx - 1:
		scounter += 1
		if scounter >= 10:
		searching = False
		else:
		searching = False
		# Clear the Flann matcher for reuse
		fl.clear()

		if scounter >= 10:
		break

		# 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
		smembership[minsy:maxsy,
		minsx:maxsx] = s_submembership
		dmembership[mindy:maxdy,
		mindx:maxdx] = d_submembership
		pcounter += 4

		smembership -= np.min(smembership)
		dmembership -= np.min(dmembership)

		if len(np.unique(smembership)) != len(np.unique(dmembership)):
		return smembership, dmembership

		# Now match the decomposed segments to one another
		for p in np.unique(smembership):
		sy_part, sx_part = np.where(smembership == p)
		dy_part, dx_part = np.where(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(e.source, e.destination, sidx, didx)
		mono_matches(e.destination, e.source, didx, sidx)

		def _add_matches(self, matches):
		"""
		Given a dataframe of matches, either append to an existing
		@@ -462,3 +672,25 @@ class Edge(dict, MutableMapping):
		total_overlap_coverage = (convex_poly.GetArea()/intersection_area)

		return total_overlap_coverage

		def decompose(self, maxiterations=3):
		"""
		Apply coupled decomposition to the images and
		match identified sub-images

		Parameters
		----------
		maxiterations : int
		The number of iterations. Appropriate values:

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


		"""
		pass

autocnet/graph/node.py

+14 −2

Changes for autocnet/graph/node.py: 14 added lines, 2 removed lines.
Original line number	Diff line number	Diff line
		@@ -154,7 +154,7 @@ class Node(dict, MutableMapping):

		return self.coverage_area

		def get_array(self, band=1):
		def get_byte_array(self, band=1):
		"""
		Get a band as a 32-bit numpy array

		@@ -167,6 +167,19 @@ class Node(dict, MutableMapping):
		array = self.geodata.read_array(band=band)
		return bytescale(array)

		def get_array(self, band=1):
		"""
		Get a band as a 32-bit numpy array

		Parameters
		----------
		band : int
		The band to read, default 1
		"""

		array = self.geodata.read_array(band=band)
		return array

		def get_keypoints(self, index=None):
		"""
		Return the keypoints for the node. If index is passed, return
		@@ -460,4 +473,3 @@ class Node(dict, MutableMapping):
		mask = panel[clean_keys].all(axis=1)
		matches = self._keypoints[mask]
		return matches, mask

autocnet/matcher/feature.py

+25 −11

Changes for autocnet/matcher/feature.py: 25 added lines, 11 removed lines.
Original line number	Diff line number	Diff line
		@@ -28,9 +28,10 @@ class FlannMatcher(object):
		def __init__(self, flann_parameters=DEFAULT_FLANN_PARAMETERS):
		self._flann_matcher = cv2.FlannBasedMatcher(flann_parameters, {})
		self.nid_lookup = {}
		self.search_idx = {}
		self.node_counter = 0

		def add(self, descriptor, nid):
		def add(self, descriptor, nid, index=None):
		"""
		Add a set of descriptors to the matcher and add the image
		index key to the image_indices attribute
		@@ -46,6 +47,10 @@ class FlannMatcher(object):
		self._flann_matcher.add([descriptor])
		self.nid_lookup[self.node_counter] = nid
		self.node_counter += 1
		if index is not None:
		self.search_idx = dict((i, j) for i, j in enumerate(index))
		else:
		self.search_idx = dict((i,i) for i in range(len(descriptor)))

		def clear(self):
		"""
		@@ -55,6 +60,7 @@ class FlannMatcher(object):
		self._flann_matcher.clear()
		self.nid_lookup = {}
		self.node_counter = 0
		self.search_idx = {}

		def train(self):
		"""
		@@ -62,7 +68,7 @@ class FlannMatcher(object):
		"""
		self._flann_matcher.train()

		def query(self, descriptor, query_image, k=3):
		def query(self, descriptor, query_image, k=3, index=None):
		"""

		Parameters
		@@ -76,6 +82,10 @@ class FlannMatcher(object):
		k : int
		The number of nearest neighbors to search for

		index : iterable
		An iterable of observation indices to utilize for
		the input descriptors

		Returns
		-------
		matched : dataframe
		@@ -86,22 +96,26 @@ class FlannMatcher(object):

		matches = self._flann_matcher.knnMatch(descriptor, k=k)
		matched = []
		for m in matches:
		for i in m:
		for i, m in enumerate(matches):
		for j in m:
		if index is not None:
		qid = index[i]
		else:
		qid = j.queryIdx
		source = query_image
		destination = self.nid_lookup[i.imgIdx]
		destination = self.nid_lookup[j.imgIdx]
		if source < destination:
		matched.append((query_image,
		i.queryIdx,
		qid,
		destination,
		i.trainIdx,
		i.distance))
		self.search_idx[j.trainIdx],
		j.distance))
		elif source > destination:
		matched.append((destination,
		i.trainIdx,
		self.search_idx[j.trainIdx],
		query_image,
		i.queryIdx,
		i.distance))
		qid,
		j.distance))
		else:
		warnings.warn('Likely self neighbor in query!')
		return pd.DataFrame(matched, columns=['source_image', 'source_idx',

autocnet/transformation/decompose.py

0 → 100644

+119 −0

Changes for autocnet/transformation/decompose.py: 119 added lines, 0 removed lines.
Original line number	Diff line number	Diff line
		import numpy as np
		from scipy.stats import pearsonr

		RADIAL_SIZE = 720
		RADIAL_STEP = 2 * np.pi / RADIAL_SIZE
		THETAS = np.round(np.arange(0, 2 * np.pi, RADIAL_STEP), 5)

		def cart2polar(x, y):
		theta = np.arctan2(y, x)
		return theta

		def index_coords(data, origin=None):
		"""Creates x & y coords for the indicies in a numpy array "data".
		"origin" defaults to the center of the image. Specify origin=(0,0)
		to set the origin to the lower left corner of the image."""
		ny, nx = data.shape[:2]
		if origin is None:
		origin_x, origin_y = nx // 2, ny // 2
		else:
		origin_x, origin_y = origin
		x, y = np.meshgrid(np.arange(nx), np.arange(ny))
		x -= origin_x
		y -= origin_y
		return x, y

		def reproject_image_into_polar(data, origin=None):
		"""Reprojects a 3D numpy array ("data") into a polar coordinate system.
		"origin" is a tuple of (x0, y0) and defaults to the center of the image."""
		ny, nx = data.shape[:2]
		if origin is None:
		origin = (nx//2, ny//2)

		# Determine that the theta coords will be
		x, y = index_coords(data, origin=origin)
		theta = cart2polar(x, y)
		# -180 to 180 conversion to 0 to 360
		theta[theta < 0] += 2 * np.pi
		return theta

		def coupled_decomposition(sdata, ddata, sorigin=(), dorigin=(), M=4, sub_skp=None):
		"""
		Apply coupled decomposition to two 2d images.

		sdata : ndarray
		(n,m) array of values to decompose

		ddata : ndarray
		(j,k) array of values to decompose

		sorigin : tuple
		in the form (x,y)

		dorigin : tuple
		in the form (x,y)

		"""

		soriginx, soriginy = sorigin
		doriginx, doriginy = dorigin

		# Create membership arrays for each input image
		smembership = np.ones(sdata.shape)
		dmembership = np.ones(ddata.shape)

		# Project the image into a polar coordinate system centered on p_{1}
		stheta = reproject_image_into_polar(sdata, origin=(int(soriginx), int(soriginy)))
		dtheta = reproject_image_into_polar(ddata, origin=(int(doriginx), int(doriginy)))

		# Compute the mean profiles for each radial slice
		smean = np.empty(RADIAL_SIZE)
		dmean = np.empty(RADIAL_SIZE)
		for i, t in enumerate(THETAS):
		# The way this method words, it is possible to get nan values in some of the steps as this is discrete
		smean[i] = np.mean(sdata[(t <= stheta) & (stheta <= t + RADIAL_STEP)])
		dmean[i] = np.mean(ddata[(t <= dtheta) & (dtheta <= t + RADIAL_STEP)])

		# Rotate the second image around the origin and compute the correlation coeff. for each 0.5 degree rotation.
		maxp = -1
		maxidx = 0
		for j in range(RADIAL_SIZE):
		dsearch = np.concatenate((dmean[j:], dmean[:j]))
		r, p = pearsonr(smean, dsearch)
		if r >= maxp:
		maxp = r
		maxidx = j

		# Maximum correlation (theta) defines the angle of rotation for the destination image
		theta = THETAS[maxidx]

		if theta <= np.pi:
		lam = theta
		else:
		lam = 2 * np.pi - theta

		# Classify the sub-images based on the decomposition size (M) and theta
		breaks = np.linspace(0, 2 * np.pi, M + 1)
		for i, t in enumerate(breaks[:-1]):
		smembership[(t <= stheta) & ( stheta <= breaks[i+1])] = i

		for i, t in enumerate(breaks[:-1]):
		# Handle the boundary crossers
		start_theta = t + theta
		stop_theta = breaks[i + 1] + theta

		if stop_theta > 2 * np.pi:
		stop_theta -= 2 * np.pi
		if start_theta > 2 * np.pi:
		start_theta -= 2 * np.pi

		if start_theta > stop_theta:
		# Handles the case where theta is a negative rotation
		dmembership[(start_theta <= dtheta) & (dtheta <= 2 * np.pi)] = i
		dmembership[(0 <= dtheta) * dtheta <= stop_theta + lam] = i
		dmembership[(start_theta <= dtheta) & (dtheta <= stop_theta)] = i
		else:
		# Handles the standard case without boundary crossers
		dmembership[(start_theta <= dtheta) & (dtheta <= stop_theta)] = i

		return smembership, dmembership