Simplified voronoi to always map source onto destination when using homography... (d0a641ee) · Commits · aflab / astrogeology / Autocnet

autocnet/cg/cg.py

+38 −0

Changes for autocnet/cg/cg.py: 38 added lines, 0 removed lines.
Original line number	Diff line number	Diff line
		import pandas as pd
		import numpy as np

		from scipy.spatial import ConvexHull

		from autocnet.utils import utils

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


		def compute_vor_weight(vor, voronoi_df, intersection_poly, verbose):
		"""

		Parameters
		----------
		vor : Voronoi
		Scipy Voronoi object

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

		intersection_poly : polygon
		Intersection polygon to use for
		clipping the voronoi diagram

		verbose : boolean
		Set to True to display the calculated voronoi diagram
		to the user
		"""
		i = 0
		poly_array = []
		for region in 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]),
		'weights'] = polygon_area
		i += 1

autocnet/graph/edge.py

+102 −0

Changes for autocnet/graph/edge.py: 102 added lines, 0 removed lines.
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 import Voronoi
		import cv2

		from autocnet.utils import utils
		from autocnet.matcher import health
		@@ -484,3 +486,103 @@ class Edge(dict, MutableMapping):
		total_overlap_coverage = (convex_poly.GetArea()/intersection_area)

		return total_overlap_coverage

		def vor(self, clean_keys=[], s=30, verbose=False):
		"""
		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
		----------
		clean_keys : list
		Of strings used to apply masks to omit correspondences

		s : int
		offset for the corners of the image

		verbose : boolean
		Set to True plots the voronoi diagram

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

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

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

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

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

		if self.source.geodata.coordinate_transformation.this is not None:
		print("Image has coordinate transform.")
		source_footprint_poly = self.source.geodata.footprint
		destination_footprint_poly = self.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 = [self.source.geodata.latlon_to_pixel(lat, lon) for lat, lon in intersect_points]

		else:
		print("Other")
		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", "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)
		cg.compute_vor_weight(vor, voronoi_df, intersection_poly, verbose)

		return vor, voronoi_df

autocnet/graph/tests/test_edge.py

+137 −0

Changes for autocnet/graph/tests/test_edge.py: 137 added lines, 0 removed lines.
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,139 @@ 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))

		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)}

		data_frame = pd.DataFrame({"x": (15, 18, 18, 12, 12),
		"y": (5, 10, 15, 15, 10)})
		weights = pd.DataFrame({"weights": (19, 28, 37.5, 37.5, 28)})

		frames = [data_frame, weights]
		weight_pd = pd.concat(frames, axis=1)

		vor = e.vor(clean_keys=[])

		for i in vor[1]:
		k = 0
		for j in vor[1][i]:
		print(i, k, j)
		self.assertAlmostEquals(j, weight_pd[i][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))

		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

		data_frame = pd.DataFrame({"x": (15, 18, 18, 12, 12),
		"y": (5, 10, 15, 15, 10)})
		weights = pd.DataFrame({"weights": (19, 28, 37.5, 37.5, 28)})

		frames = [data_frame, weights]
		weight_pd = pd.concat(frames, axis=1)

		vor = e.vor(clean_keys=[])

		for i in vor[1]:
		k = 0
		for j in vor[1][i]:
		self.assertAlmostEquals(j, weight_pd[i][k])
		k += 1

autocnet/transformation/transformations.py

+4 −0

Changes for autocnet/transformation/transformations.py: 4 added lines, 0 removed lines.
Original line number	Diff line number	Diff line
		@@ -457,6 +457,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

autocnet/utils/utils.py

+46 −0

Changes for autocnet/utils/utils.py: 46 added lines, 0 removed lines.
Original line number	Diff line number	Diff line
		@@ -334,3 +334,49 @@ def array_to_poly(array):
		poly = ogr.CreateGeometryFromJson(json.dumps(geom))
		return poly


		def reproj_corner(H, corner):
		"""
		Reproject a pixel in one image into another image

		Parameters
		----------
		H : object
		(3,3) ndarray or Homography object

		corner : iterable
		A 2 element iterable in the form x, y
		"""
		if len(corner) == 2:
		coords = np.array([corner[0], corner[1], 1])
		elif len(corner) == 3:
		coords = np.asarray(corner)
		coords *= coords[-1]
		return H.dot(coords)


		def scale_point(point, centroid, scale):
		"""
		Given a point, centroid, and a scalar scales the given pointer
		around the given centroid

		Parameters
		----------
		point : tuple
		(x, y) coordinates for a given point

		centroid : tuple
		(x, y, 0) coordinates for the centroid of a polygon

		scale : int
		The multiplier to scale by

		Returns
		-------
		vector : ndarray
		(2, 1) array of the scaled point returned in x, y form
		"""
		point = np.asarray(point)
		centroid = centroid[:2]
		vector = ((point - centroid)*scale) + centroid
		return vector