Merge pull request #170 from acpaquette/voronoi (b4c1ed4c) · Commits · aflab / astrogeology / Autocnet

autocnet/cg/cg.py

+117 −0

Original line number	Diff line number	Diff line
		import pandas as pd
		import numpy as np

		from scipy.spatial import ConvexHull
		from scipy.spatial import Voronoi
		import cv2

		from autocnet.utils import utils


		def convex_hull_ratio(points, ideal_area):
		@@ -99,3 +104,115 @@ def get_area(poly1, poly2):
		"""
		intersection_area = poly1.Intersection(poly2).GetArea()
		return intersection_area


		def vor(edge, clean_keys=[], s=30):
		"""
		Creates a voronoi diagram for an edge using either the coordinate
		transformation or using the homography between source and destination.

		The coordinate transformation uses the footprint of source and destination to
		calculate an intersection between the two images, then transforms the vertices of
		the intersection back into pixel space.

		If a coordinate transform does not exist, use the homography to project the destination image
		onto the source image, producing an area of intersection.

		The intersection vertices are then scaled by a factor of s (default 30), this accounts for the
		areas of the voronoi that would be missed if the scaled vertices were not included into the
		voronoi calculation.


		Parameters
		----------
		edge : object
		An edge object

		clean_keys : list
		Of strings used to apply masks to omit correspondences

		s : int
		offset for the corners of the image

		Returns
		-------
		vor : object
		Scipy Voronoi object

		voronoi_df : dataframe
		3 column pandas dataframe of x, y, and weights

		"""
		source_corners = edge.source.geodata.xy_corners
		destination_corners = edge.destination.geodata.xy_corners

		matches, _ = edge.clean(clean_keys=clean_keys)

		source_keypoints_pd = edge.source.get_keypoint_coordinates(index=matches['source_idx'],
		homogeneous=True)
		destination_keypoints_pd = edge.destination.get_keypoint_coordinates(index=matches['destination_idx'],
		homogeneous=True)

		if edge.source.geodata.coordinate_transformation.this is not None:
		source_footprint_poly = edge.source.geodata.footprint
		destination_footprint_poly = edge.destination.geodata.footprint

		intersection_poly = destination_footprint_poly.Intersection(source_footprint_poly)
		intersection_geom = intersection_poly.GetGeometryRef(0)
		intersect_points = intersection_geom.GetPoints()

		intersection_points = [edge.source.geodata.latlon_to_pixel(lat, lon) for lat, lon in intersect_points]

		else:
		H, mask = cv2.findHomography(destination_keypoints_pd.values,
		source_keypoints_pd.values,
		cv2.RANSAC,
		2.0)

		proj_corners = []
		for c in destination_corners:
		x, y, h = utils.reproj_corner(H, c)
		x /= h
		y /= h
		h /= h
		proj_corners.append((x, y))

		orig_poly = utils.array_to_poly(source_corners)
		proj_poly = utils.array_to_poly(proj_corners)

		intersection_poly = orig_poly.Intersection(proj_poly)
		intersection_geom = intersection_poly.GetGeometryRef(0)
		intersection_points = intersection_geom.GetPoints()

		centroid = intersection_poly.Centroid().GetPoint()

		voronoi_df = pd.DataFrame(data=source_keypoints_pd, columns=["x", "y", "vor_weights"])

		voronoi_df["x"] = source_keypoints_pd['x']
		voronoi_df["y"] = source_keypoints_pd['y']
		keypoints = np.asarray(source_keypoints_pd)

		inters = np.empty((len(intersection_points), 2))
		for g, (i, j) in enumerate(intersection_points):
		scaledx, scaledy = utils.scale_point((i, j), centroid, s)
		point = np.array([scaledx, scaledy])
		inters[g] = point

		keypoints = np.vstack((keypoints[:, :2], inters))
		vor = Voronoi(keypoints)

		poly_array = []
		for i, region in enumerate(vor.regions):
		region_point = vor.points[np.argwhere(vor.point_region==i)]
		if -1 not in region:
		polygon_points = [vor.vertices[i] for i in region]
		if len(polygon_points) != 0:
		polygon = utils.array_to_poly(polygon_points)
		intersection = polygon.Intersection(intersection_poly)
		poly_array = np.append(poly_array, intersection)
		polygon_area = intersection.GetArea()
		voronoi_df.loc[(voronoi_df["x"] == region_point[0][0][0]) &
		(voronoi_df["y"] == region_point[0][0][1]),
		'vor_weights'] = polygon_area

		return vor, voronoi_df

autocnet/examples/Apollo15/cube_adjacency.json

0 → 100644

+3 −0

{"AS15-M-0297_crop.cub": ["AS15-M-0298_crop.cub", "AS15-M-0299_crop.cub"],

 "AS15-M-0298_crop.cub": ["AS15-M-0299_crop.cub", "AS15-M-0297_crop.cub"],

 "AS15-M-0299_crop.cub": ["AS15-M-0298_crop.cub", "AS15-M-0297_crop.cub"]}

autocnet/graph/edge.py

+19 −0

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
		@@ -769,3 +770,21 @@ class Edge(dict, MutableMapping):
		total_overlap_coverage = (convex_poly.GetArea()/intersection_area)

		return total_overlap_coverage

		def compute_weights(self, clean_keys, **kwargs):
		"""
		Computes a voronoi diagram for the overlap between two images
		then gets the area of each polygon resulting in a voronoi weight.
		These weights are then appended to the matches dataframe.

		Parameters
		----------
		clean_keys : list
		Of strings used to apply masks to omit correspondences

		"""
		if self.matches is None:
		raise AttributeError('Matches have not been computed for this edge')
		voronoi = cg.vor(self, clean_keys, **kwargs)
		self.matches = pd.concat([self.matches, voronoi[1]['vor_weights']], axis=1)

autocnet/graph/tests/test_edge.py

+126 −0

Original line number	Diff line number	Diff line
		@@ -8,6 +8,7 @@ from plio.io import io_gdal

		from autocnet.examples import get_path
		from autocnet.graph.network import CandidateGraph
		from autocnet.utils.utils import array_to_poly

		from .. import edge
		from .. import node
		@@ -129,3 +130,128 @@ class TestEdge(unittest.TestCase):

		self.assertRaises(AttributeError, cg.edge[0][1].coverage)
		self.assertEqual(e.coverage(), 0.3)

		def test_voronoi_transform(self):
		keypoint_df = pd.DataFrame({'x': (15, 18, 18, 12, 12), 'y': (5, 10, 15, 15, 10)})
		keypoint_matches = [[0, 0, 1, 0],
		[0, 1, 1, 1],
		[0, 2, 1, 2],
		[0, 3, 1, 3],
		[0, 4, 1, 4]]

		matches_df = pd.DataFrame(data=keypoint_matches, columns=['source_image', 'source_idx',
		'destination_image', 'destination_idx'])
		e = edge.Edge()

		e.clean = MagicMock(return_value=(matches_df, None))
		e.matches = matches_df

		source_node = MagicMock(spec=node.Node())
		destination_node = MagicMock(spec=node.Node())

		source_node.get_keypoint_coordinates = MagicMock(return_value=keypoint_df)
		destination_node.get_keypoint_coordinates = MagicMock(return_value=keypoint_df)

		e.source = source_node
		e.destination = destination_node

		source_geodata = Mock(spec=io_gdal.GeoDataset)
		destination_geodata = Mock(spec=io_gdal.GeoDataset)

		e.source.geodata = source_geodata
		e.destination.geodata = destination_geodata

		source_corners = [(0, 0),
		(20, 0),
		(20, 20),
		(0, 20)]

		destination_corners = [(10, 5),
		(30, 5),
		(30, 25),
		(10, 25)]

		source_poly = array_to_poly(source_corners)
		destination_poly = array_to_poly(destination_corners)

		def latlon_to_pixel(i, j):
		return vals[(i, j)]

		e.source.geodata.latlon_to_pixel = MagicMock(side_effect=latlon_to_pixel)
		e.destination.geodata.latlon_to_pixel = MagicMock(side_effect=latlon_to_pixel)

		e.source.geodata.footprint = source_poly
		e.source.geodata.xy_corners = source_corners
		e.destination.geodata.footprint = destination_poly
		e.destination.geodata.xy_corners = destination_corners

		vals = {(10, 5): (10, 5), (20, 5): (20, 5), (20, 20): (20, 20), (10, 20): (10, 20)}

		weights = pd.DataFrame({"vor_weights": (19, 28, 37.5, 37.5, 28)})

		e.compute_weights(clean_keys=[])

		k = 0
		for i in e.matches['vor_weights']:
		self.assertAlmostEquals(i, weights['vor_weights'][k])
		k += 1

		def test_voronoi_homography(self):
		source_keypoint_df = pd.DataFrame({'x': (15, 18, 18, 12, 12), 'y': (5, 10, 15, 15, 10)})
		destination_keypoint_df = pd.DataFrame({'x': (5, 8, 8, 2, 2), 'y': (0, 5, 10, 10, 5)})
		keypoint_matches = [[0, 0, 1, 0],
		[0, 1, 1, 1],
		[0, 2, 1, 2],
		[0, 3, 1, 3],
		[0, 4, 1, 4]]

		matches_df = pd.DataFrame(data = keypoint_matches, columns=['source_image', 'source_idx',
		'destination_image', 'destination_idx'])
		e = edge.Edge()

		e.clean = MagicMock(return_value=(matches_df, None))
		e.matches = matches_df

		source_node = MagicMock(spec=node.Node())
		destination_node = MagicMock(spec=node.Node())

		source_node.get_keypoint_coordinates = MagicMock(return_value=source_keypoint_df)
		destination_node.get_keypoint_coordinates = MagicMock(return_value=destination_keypoint_df)

		e.source = source_node
		e.destination = destination_node

		source_geodata = Mock(spec=io_gdal.GeoDataset)
		destination_geodata = Mock(spec=io_gdal.GeoDataset)

		e.source.geodata = source_geodata
		e.destination.geodata = destination_geodata

		source_corners = [(0, 0),
		(20, 0),
		(20, 20),
		(0, 20)]

		destination_corners = [(0, 0),
		(20, 0),
		(20, 20),
		(0, 20)]

		e.source.geodata.coordinate_transformation.this = None
		e.destination.geodata.coordinate_transformation.this = None

		e.source.geodata.xy_corners = source_corners
		e.destination.geodata.xy_corners = destination_corners

		weights = pd.DataFrame({"vor_weights": (19, 28, 37.5, 37.5, 28)})

		e.compute_weights(clean_keys=[])

		k = 0
		for i in e.matches['vor_weights']:
		self.assertAlmostEquals(i, weights['vor_weights'][k])
		k += 1

autocnet/transformation/transformations.py

+4 −0

Original line number	Diff line number	Diff line
		@@ -411,6 +411,10 @@ class Homography(TransformationMatrix):
		def error(self):
		return self.compute_error(self.x1, self.x2)

		@property
		def inverse(self):
		return np.linalg.inv(self)

		def compute_error(self, a, b, mask=None):
		"""
		Give this homography, compute the planar reprojection error