Added mosaic match code (#277) (e93af564) · Commits · aflab / astrogeology / Autocnet

.gitignore

+2 −0

Original line number	Diff line number	Diff line
		@@ -81,10 +81,12 @@ target/
		*.cnet
		*.lis
		*.list
		*.aux.xml

		# Pytest
		lastfailed
		nodeids
		.pytest_cache

		# VSCode
		launch.json

autocnet/init.py

+39 −10

Original line number	Diff line number	Diff line
		import os
		import warnings
		import yaml
		import socket

		from pkg_resources import get_distribution, DistributionNotFound
		from sqlalchemy import create_engine, pool, orm
		from sqlalchemy.event import listen

		import autocnet
		import autocnet.examples
		import autocnet.camera
		import autocnet.cg
		import autocnet.control
		import autocnet.graph
		import autocnet.matcher
		import autocnet.transformation
		import autocnet.utils
		from pkg_resources import get_distribution, DistributionNotFound

		try:
		_dist = get_distribution('autocnet')
		@@ -26,6 +21,40 @@ except DistributionNotFound:
		else:
		__version__ = _dist.version

		#Load the config file and setup a global DB session factory
		try:
		with open(os.environ['autocnet_config'], 'r') as f:
		config = yaml.safe_load(f)
		except:
		warnings.warn('No autocnet_config environment variable set. Defaulting to an en empty configuration.')
		config = {}

		try:
		db_uri = '{}://{}:{}@{}:{}/{}'.format(config['database']['type'],
		config['database']['username'],
		config['database']['password'],
		config['database']['host'],
		config['database']['pgbouncer_port'],
		config['database']['name'])
		hostname = socket.gethostname()
		engine = create_engine(db_uri, poolclass=pool.NullPool,
		connect_args={"application_name":"AutoCNet_{}".format(hostname)},
		isolation_level="AUTOCOMMIT")
		Session = orm.session.sessionmaker(bind=engine)
		except:
		Session = None
		engine = None

		import autocnet.examples
		import autocnet.camera
		import autocnet.cg
		import autocnet.control
		import autocnet.graph
		import autocnet.matcher
		import autocnet.transformation
		import autocnet.utils
		import autocnet.spatial

		# Patch the candidate graph into the root namespace
		from autocnet.graph.network import CandidateGraph

autocnet/cg/cg.py

+187 −2

Original line number	Diff line number	Diff line
		@@ -11,8 +11,6 @@ from scipy.spatial import Voronoi
		import shapely.geometry
		from shapely.geometry import Polygon, Point
		from shapely.affinity import scale
		# from shapely.ops import unary_union
		import cv2

		from autocnet.utils import utils

		@@ -200,3 +198,190 @@ def compute_voronoi(keypoints, intersection=None, geometry=False, s=30): # ADDED
		i += 1

		return voronoi_df

		def single_centroid(geom):
		"""
		For a geom, return the centroid

		Parameters
		----------
		geom : shapely.geom object

		Returns
		-------

		valid : list
		in the form [(x,y)]
		"""
		x, y = geom.centroid.xy
		valid = [(x[0],y[0])]
		return valid

		def nearest(pt, search):
		"""
		Fine the index of nearest (Euclidean) point in a list
		of points.

		Parameters
		----------
		pt : ndarray
		(2,1) array
		search : ndarray
		(n,2) array of points to search within. The
		returned index is the closet point in this set
		to the search

		Returns
		-------
		: int
		The index to the nearest point.
		"""
		return np.argmin(np.sum((search - pt)**2, axis=1))

		def distribute_points(geom, nspts, ewpts):
		"""
		This is a decision tree that attempts to perform a
		very simplistic approximation of the shape
		of the geometry and then place some number of
		north/south and east/west points into the geometry.

		Parameters
		----------
		geom : shapely.geom
		A shapely geometry object

		nspts : int
		The number of points to attempt to place
		in the N/S (up/down) direction

		ewpts : int
		The number of points to attempt to place
		in the E/W (right/left) direction

		Returns
		-------
		valid : list
		of point coordinates in the form [(x1,y1), (x2,y2), ..., (xn, yn)]
		"""
		geom_coords = np.column_stack(geom.exterior.xy)

		coords = np.array(list(zip(*geom.envelope.exterior.xy))[:-1])

		ll = coords[0]
		lr = coords[1]
		ur = coords[2]
		ul = coords[3]

		# Find the points nearest the ul and ur
		ul_actual = geom_coords[nearest(ul, geom_coords)]
		ur_actual = geom_coords[nearest(ur, geom_coords)]
		dist = np.sqrt((ul_actual[1] - ur_actual[1])2 + (ul_actual[0] - ur_actual[0])2)
		m = (ul_actual[1]-ur_actual[1])/(ul_actual[0]-ur_actual[0])
		b = (ul_actual[1] - ul_actual[0] * m)
		newtop = []
		xnodes = np.linspace(ul_actual[0], ur_actual[0], num=ewpts+2)
		for x in xnodes[1:-1]:
		newtop.append((x, m*x+b))

		# Find the points nearest the ll and lr

		ll_actual = geom_coords[nearest(ll, geom_coords)]
		lr_actual = geom_coords[nearest(lr, geom_coords)]
		dist = np.sqrt((ll_actual[1] - lr_actual[1])2 + (ll_actual[0] - lr_actual[0])2)
		m = (ll_actual[1]- lr_actual[1])/(ll_actual[0]-lr_actual[0])
		b = (ll_actual[1] - ll_actual[0] * m)
		newbot = []
		xnodes = np.linspace(ll_actual[0], lr_actual[0], num=ewpts+2)
		for x in xnodes[1:-1]:
		newbot.append((x, m*x+b))
		newpts = []
		for i in range(len(newtop)):
		top = newtop[i]
		bot = newbot[i]
		# Compute the line between top and bottom
		m = (top[1] - bot[1]) / (top[0] - bot[0])
		b = (top[1] - top[0] * m)
		xnodes = np.linspace(bot[0], top[0], nspts+2)
		for x in xnodes[1:-1]:
		newpts.append((x, m*x+b))
		valid = []
		# Perform a spatial intersection check to eject points that are not valid
		for p in newpts:
		pt = Point(p[0], p[1])
		if geom.contains(pt):
		valid.append(p)
		return valid

		def distribute_points_in_geom(geom):
		"""
		Given a geometry, attempt a basic classification of the shape.
		RIght now, this simply attempts to determine if the bounding box
		is generally N/S or generally E/W trending. Once the determination
		is made, the algorithm places points in the geometry and returns
		a list of valid (intersecting) points.

		Parameters
		----------
		geom : shapely.geom object
		The geometry object

		Returns
		-------
		valid : list
		of valid points in the form (x,y) or (lon,lat)

		"""
		coords = list(zip(*geom.envelope.exterior.xy))
		short = np.inf
		long = -np.inf
		shortid = 0
		longid = 0
		for i, p in enumerate(coords[:-1]):
		d = np.sqrt((coords[i+1][0] - p[0])2+(coords[i+1][1]-p[1])2)
		if d < short:
		short = d
		shortid = i
		if d > long:
		long = d
		longid = i
		ratio = short/long
		ns = False
		ew = False
		valid = []

		# The polygons should be encoded with a lower left origin in counter-clockwise direction.
		# Therefore, if the 'bottom' is the short edge it should be id 0 and modulo 2 == 0.
		if shortid % 2 == 0:
		ns = True
		elif longid % 2 == 0:
		ew = True

		# Decision Tree
		if ratio < 0.16 and geom.area < 0.01:
		# Class: Slivers - ignore.
		return
		elif geom.area <= 0.004 and ratio >= 0.25:
		# Single point at the centroid
		valid = single_centroid(geom)
		elif ns==True:
		# Class, north/south poly, multi-point
		nspts = int(round(long, 1) * 10)
		if nspts >= 5:
		nspts = int(nspts/(nspts/3))
		ewpts = max(int(round(short, 1) * 5), 1)
		if nspts == 1 and ewpts == 1:
		valid = single_centroid(geom)
		else:
		valid = distribute_points(geom, nspts, ewpts)
		elif ew == True:
		# Since this is an LS, we should place these diagonally from the 'lower left' to the 'upper right'
		nspts = max(int(round(short, 1) * 5), 1)
		if nspts >= 5:
		nspts = int(nspts/(npts/3))
		ewpts = int(round(long, 1) * 10)
		if nspts == 1 and ewpts == 1:
		valid = single_centroid(geom)
		else:
		valid = distribute_points(geom, nspts, ewpts)

		return valid

autocnet/graph/edge.py

+376 −7

Original line number	Diff line number	Diff line
		from functools import wraps, singledispatch
		import warnings
		from collections import MutableMapping
		from collections import defaultdict, MutableMapping, Counter

		from geoalchemy2.elements import WKBElement
		import numpy as np
		import pandas as pd
		import networkx as nx

		import pyproj
		from scipy.spatial.distance import cdist
		from shapely.geometry import Point
		import sqlalchemy

		import autocnet
		from autocnet import Session, engine
		from autocnet.graph.node import Node
		from autocnet.utils import utils
		from autocnet.matcher import cpu_outlier_detector as od
		@@ -20,6 +23,10 @@ from autocnet.transformation import homography as hm
		from autocnet.transformation import spatial
		from autocnet.vis.graph_view import plot_edge, plot_node, plot_edge_decomposition, plot_matches
		from autocnet.cg import cg
		from autocnet.io.db.model import Images, Keypoints, Matches,\
		Cameras, Base, Overlay, Edges,\
		Costs, Measures, Points, Measures
		from autocnet.io.db.wrappers import DbDataFrame

		from plio.io.io_gdal import GeoDataset
		from csmapi import csmapi
		@@ -340,11 +347,9 @@ class Edge(dict, MutableMapping):
		s_keypoints, d_keypoints = self.get_match_coordinates(clean_keys=clean_keys)
		self.fundamental_matrix, fmask = fm.compute_fundamental_matrix(s_keypoints, d_keypoints, **kwargs)


		if isinstance(self.fundamental_matrix, np.ndarray):
		# Convert the truncated RANSAC mask back into a full length mask
		mask[mask] = fmask

		# Set the initial state of the fundamental mask in the masks
		self.masks[maskname] = mask

		@@ -659,7 +664,7 @@ class Edge(dict, MutableMapping):
		"""
		if not isinstance(self.matches, pd.DataFrame):
		raise AttributeError('Matches have not been computed for this edge')
		voronoi = cg.vor(self, clean_keys, **kwargs)
		voronoi = cg.compute_voronoi(self, clean_keys, **kwargs)
		self.matches = pd.concat([self.matches, voronoi[1]['vor_weights']], axis=1)

		def compute_overlap(self, buffer_dist=0, **kwargs):
		@@ -704,8 +709,372 @@ class Edge(dict, MutableMapping):
		def get_matches(self, clean_keys=[]): # pragma: no cover
		if self.matches.empty:
		return pd.DataFrame()
		self.add_coordinates_to_matches()
		#self.add_coordinates_to_matches()
		matches, _ = self.clean(clean_keys=clean_keys)
		skps = matches[['source_x', 'source_y']]
		dkps = matches[['destination_x', 'destination_y']]
		return matches

		class NetworkEdge(Edge):

		default_msg = {'sidx':None,
		'didx':None,
		'task':None,
		'param_step':0,
		'success':False}

		def __init__(self, args, *kwargs):
		super(NetworkEdge, self).__init__(args, *kwargs)
		self.job_status = defaultdict(dict)

		def _from_db(self, table_obj):
		"""
		Generic database query to pull the row associated with this node
		from an arbitrary table. We assume that the row id matches the node_id.

		Parameters
		----------
		table_obj : object
		The declared table class (from db.model)
		"""
		session = Session()
		res = session.query(table_obj).\
		filter(table_obj.source == self.source['node_id']).\
		filter(table_obj.destination == self.destination['node_id'])
		session.close()
		return res

		@property
		def masks(self):
		res = Session().query(Edges.masks).\
		filter(Edges.source == self.source['node_id']).\
		filter(Edges.destination == self.destination['node_id']).\
		first()

		try:
		df = pd.DataFrame.from_records(res[0])
		df.index = df.index.map(int)
		except:
		ids = list(map(int, self.matches.index.values))
		df = pd.DataFrame(index=ids)
		df.index.name = 'match_id'
		return DbDataFrame(df, parent=self, name='masks')

		@masks.setter
		def masks(self, v):

		def dict_check(input):
		for k, v in input.items():
		if isinstance(v, dict):
		dict_check(v)
		elif v is None:
		continue
		elif np.isnan(v):
		input[k] = None


		df = pd.DataFrame(v)
		session = Session()
		res = session.query(Edges).\
		filter(Edges.source == self.source['node_id']).\
		filter(Edges.destination == self.destination['node_id']).first()
		if res:
		as_dict = df.to_dict()
		dict_check(as_dict)
		# Update the masks
		res.masks = as_dict
		session.add(res)
		session.commit()

		@property
		def costs(self):
		# these are np.float coming out, sqlalchemy needs ints
		ids = list(map(int, self.matches.index.values))
		res = Session().query(Costs).filter(Costs.match_id.in_(ids)).all()
		#qf = q.filter(Costs.match_id.in_(ids))

		if res:
		# Parse the JSON dicts in the cost field into a full dimension dataframe
		costs = {r.match_id:r._cost for r in res}
		df = pd.DataFrame.from_records(costs).T # From records is important because from_dict drops rows with empty dicts
		else:
		df = pd.DataFrame(index=ids)

		df.index.name = 'match_id'
		return DbDataFrame(df, parent=self, name='costs')


		@costs.setter
		def costs(self, v):
		to_db_add = []
		# Get the query obj
		session = Session()
		q = session.query(Costs)
		# Need the new instance here to avoid __setattr__ issues
		df = pd.DataFrame(v)
		for idx, row in df.iterrows():
		# Now invert the expanded dict back into a single JSONB column for storage
		res = q.filter(Costs.match_id == idx).first()
		if res:
		#update the JSON blob
		costs_new_or_updated = row.to_dict()
		for k, v in costs_new_or_updated.items():
		if v is None:
		continue
		elif np.isnan(v):
		v = None
		res._cost[k] = v
		sqlalchemy.orm.attributes.flag_modified(res, '_cost')
		session.add(res)
		session.commit()
		else:
		row = row.to_dict()
		costs = row.pop('_costs', {})
		for k, v in row.items():
		if np.isnan(v):
		v = None
		costs[k] = v
		cost = Costs(match_id=idx, _cost=costs)
		to_db_add.append(cost)
		if to_db_add:
		session.bulk_save_objects(to_db_add)
		session.commit()

		@property
		def matches(self):
		session = Session()
		q = session.query(Matches)
		qf = q.filter(Matches.source == self.source['node_id'],
		Matches.destination == self.destination['node_id'])
		odf = pd.read_sql(qf.statement, q.session.bind).set_index('id')
		df = pd.DataFrame(odf.values, index=odf.index.values, columns=odf.columns.values)
		df.index.name = 'id'
		# Explicit close to get the session cleaned up
		session.close()
		return DbDataFrame(df,
		parent=self,
		name='matches')

		@matches.setter
		def matches(self, v):
		to_db_add = []
		to_db_update = []
		df = pd.DataFrame(v)
		df.index.name = v.index.name
		# Get the query obj
		session = Session()
		q = session.query(Matches)
		for idx, row in df.iterrows():
		# Determine if this is an update or the addition of a new row
		if hasattr(row, 'id'):
		res = q.filter(Matches.id == row.id).first()
		match_id = row.id
		elif v.index.name == 'id':
		res = q.filter(Matches.id == row.name).first()
		match_id = row.name
		else:
		res = None
		if res:
		# update
		mapping = {}
		mapping['id'] = match_id
		for index in row.index:
		row_val = row[index]
		if isinstance(row_val, (np.int,)):
		row_val = int(row_val)
		elif isinstance(row_val, (np.float,)):
		row_val = float(row_val)
		elif isinstance(row_val, WKBElement):
		continue
		mapping[index] = row_val
		to_db_update.append(mapping)
		else:
		match = Matches()
		# Dynamically iterate over the columns and if the match has an
		# attribute with the column name, set it.
		for c in df.columns:
		if hasattr(match, c):
		setattr(match, c, row[c])
		to_db_add.append(match)
		if to_db_add:
		session.bulk_save_objects(to_db_add)
		if to_db_update:
		session.bulk_update_mappings(Matches, to_db_update)
		session.commit()

		@matches.deleter
		def matches(self):
		session = Session()
		session.query(Matches).filter(Matches.source == self.source['node_id'], Matches.destination == self.destination['node_id']).delete()
		session.commit()
		session.close()
		return

		@property
		def ring(self):
		res = self._from_db(Edges).first()
		if res:
		return res.ring
		return

		@ring.setter
		def ring(self, ring):
		# Setters need a single session and so should not make use of the
		# syntax sugar _from_db
		session = Session()
		res = session.query(Edges).\
		filter(Edges.source == self.source['node_id']).\
		filter(Edges.destination == self.destination['node_id']).first()
		if res:
		res.ring = ring
		else:
		edge = Edges(source=self.source['node_id'],
		destination=self.destination['node_id'],
		ring=ring)
		session.add(edge)
		session.commit()
		return

		@property
		def intersection(self):
		if not hasattr(self, '_intersection'):
		s_fp = self.source.footprint
		d_fp = self.destination.footprint
		self._intersection = s_fp.intersection(d_fp)
		return self._intersection

		@property
		def fundamental_matrix(self):
		res = self._from_db(Edges).first()
		if res:
		return np.asarray(res.fundamental)

		@fundamental_matrix.setter
		def fundamental_matrix(self, v):
		session = Session()
		res = session.query(Edges).\
		filter(Edges.source == self.source['node_id']).\
		filter(Edges.destination == self.destination['node_id']).first()
		if res:
		res.fundamental = v
		else:
		edge = Edges(source=self.source['node_id'],
		destination=self.destination['node_id'],
		fundamental = v)
		session.add(edge)
		session.commit()

		def get_overlapping_indices(self, kps):
		ecef = pyproj.Proj(proj='geocent',
		a=self.parent.config['spatial']['semimajor_rad'],
		b=self.parent.config['spatial']['semiminor_rad'])
		lla = pyproj.Proj(proj='longlat',
		a=self.parent.config['spatial']['semiminor_rad'],
		b=self.parent.config['spatial']['semimajor_rad'])
		lons, lats, alts = pyproj.transform(ecef, lla, kps.xm.values, kps.ym.values, kps.zm.values)
		points = [Point(lons[i], lats[i]) for i in range(len(lons))]
		mask = [i for i in range(len(points)) if self.intersection.contains(points[i])]
		return mask

		@property
		def measures(self):
		return Session().query(Measures).filter(sqlalchemy.or_(Measures.imageid == self.source['node_id'], Measures.imageid == self.destination['node_id'])).all()

		def network_to_matches(self, active_point=True, active_measure=True, rejected_jigsaw=False):
		"""
		For the edge, take any points/measures that are in the database and
		convert them into matches on the associated edge.

		Parameters
		----------
		active_point : bool
		If True (default) only select the points that are
		currently set to active.

		active_measure : bool
		If True (default) only add the measures that are
		currently active

		rejected_jigsaw : bool
		If False (default) add any points that are not
		set to jigsaw rejected.

		"""
		source = self.source['node_id']
		destin = self.destination['node_id']

		if source > destin:
		source, destin = destin, source

		q = Session().query(Points.id,
		Points.pointtype,
		Measures.id.label('mid'),
		Measures.sample,
		Measures.line,
		Measures.measuretype,
		Measures.imageid).\
		filter(Points.active==active_point,
		Measures.active==active_measure,
		Measures.jigreject==rejected_jigsaw,
		sqlalchemy.or_(Measures.imageid==source,
		Measures.imageid==destin)).join(Measures)

		df = pd.read_sql(q.statement, engine)
		matches = []
		columns = ['point_id', 'source_measure_id', 'destin_measure_id', 'source', 'source_idx', 'destination', 'destination_idx',
		'lat', 'lon', 'geom', 'source_x', 'source_y', 'destination_x',
		'destination_y', 'shift_x', 'shift_y', 'original_destination_x',
		'original_destination_y']

		def net2matches(grp, matches, source, destin):
		# Grab the image ids and then get the cartesian product of the ids to know which
		# edges to put the matches onto
		if len(grp) != 2:
		return

		imagea = grp[grp['imageid'] == source].iloc[0]
		imageb = grp[grp['imageid'] == destin].iloc[0]
		match = [int(imagea.id), int(imagea.mid), int(imageb.mid),
		source, 0, destin, 0,
		None, None, None,
		imagea['sample'], imagea['line'], imageb['sample'], imageb['line'],
		None, None, None, None]
		matches.append(match)

		df.groupby('id').apply(net2matches, matches, source, destin)
		self.matches = pd.DataFrame(matches, columns=columns)

		def mask_to_counter(self, mask):
		"""
		Take a mask on an edge and convert the mask into a counter where
		the key is the match id and value is 1 (the match is flagged false).

		TODO: Allow the mask to be an iterable (list). The caller of this should
		then worry about normalization as n-mask strings can come in and we
		cannot anticipate how the user might want to normalize the return.

		Parameters
		----------
		mask : str
		The name of the mask

		Returns
		-------
		: collections.Counter
		With keys equal to the indices of the False matches
		and values equal to one

		"""
		mask = self.masks[mask]
		matches_to_disable = mask[mask == False].index
		session = Session()

		bad = {}
		for o in session.query(Matches).filter(Matches.id.in_(matches_to_disable)).all():
		# This can't just set both to False, we loose a ton of good points - a bad point in 1 image is not
		# necessarily bad in all of the other images. Doing it this way assumes that it is...
		bad[o.source_measure_id] = 1
		bad[o.destin_measure_id] = 1

		return Counter(bad)

autocnet/graph/network.py

+312 −11

File changed.

Preview size limit exceeded, changes collapsed.