Fixes #426 (#431) (76a91a76) · Commits · aflab / astrogeology / Autocnet

autocnet/matcher/subpixel.py

+86 −56

Original line number	Diff line number	Diff line
		@@ -11,6 +11,7 @@ from autocnet.matcher.naive_template import pattern_match, pattern_match_autoreg
		from autocnet.matcher import ciratefi
		from autocnet.io.db.model import Measures, Points, Images, JsonEncoder
		from autocnet.graph.node import NetworkNode
		from autocnet.transformation import roi

		import geopandas as gpd
		import pandas as pd
		@@ -68,12 +69,15 @@ def check_image_size(imagesize):
		imagesize : tuple
		in the form (size_x, size_y)
		"""
		x = imagesize[0] / 2
		y = imagesize[1] / 2
		x = imagesize[0]
		y = imagesize[1]

		if x % 2 == 0:
		x += 1
		if y % 2 == 0:
		y += 1
		x = floor(x/2)
		y = floor(y/2)
		return x,y

		def clip_roi(img, center_x, center_y, size_x=200, size_y=200):
		@@ -133,7 +137,10 @@ def clip_roi(img, center_x, center_y, size_x=200, size_y=200):
		return None, 0, 0
		return subarray, axr, ayr

		def subpixel_phase(template, search, **kwargs):
		def subpixel_phase(sx, sy, dx, dy,
		s_img, d_img,
		image_size=(251, 251),
		**kwargs):
		"""
		Apply the spectral domain matcher to a search and template image. To
		shift the images, the x_shift and y_shift, need to be subtracted from
		@@ -149,7 +156,7 @@ def subpixel_phase(template, search, **kwargs):
		search : ndarray
		The search image

		Returnsslurm-2235260_89.out.2235260_89.out
		Returns
		-------
		x_offset : float
		Shift in the x-dimension
		@@ -160,10 +167,48 @@ def subpixel_phase(template, search, **kwargs):
		strength : tuple
		With the RMSE error and absolute difference in phase
		"""
		if not template.shape == search.shape:
		raise ValueError('Both the template and search images must be the same shape.')
		(y_shift, x_shift), error, diffphase = register_translation(search, template, **kwargs)
		return x_shift, y_shift, (error, diffphase)
		image_size = check_image_size(image_size)

		s_roi = roi.Roi(s_img, sx, sy, size_x=image_size[0], size_y=image_size[1])
		d_roi = roi.Roi(d_img, dx, dy, size_x=image_size[0], size_y=image_size[1])

		s_image = s_roi.clip()
		d_template = d_roi.clip()

		if s_image.shape != d_template.shape:

		s_size = s_image.shape
		d_size = d_template.shape
		updated_size_x = int(min(s_size[1], d_size[1]))
		updated_size_y = int(min(s_size[0], d_size[0]))

		# Have to subtract 1 from even entries or else the round up that
		# occurs when the size is split over the midpoint causes the
		# size to be too large by 1.
		if updated_size_x % 2 == 0:
		updated_size_x -= 1
		if updated_size_y % 2 == 0:
		updated_size_y -= 1

		# Since the image is smaller than the requested size, set the size to
		# the current maximum image size and reduce from there on potential
		# future iterations.
		size = check_image_size((updated_size_x, updated_size_y))
		s_roi = roi.Roi(s_img, sx, sy,
		size_x=size[0], size_y=size[1])
		d_roi = roi.Roi(d_img, dx, dy,
		size_x=size[0], size_y=size[1])
		s_image = s_roi.clip()
		d_template = d_roi.clip()

		if (s_image is None) or (d_template is None):
		return None, None, None

		(shift_y, shift_x), error, diffphase = register_translation(s_image, d_template, **kwargs)
		dx = d_roi.x - shift_x
		dy = d_roi.y - shift_y

		return dx, dy, (error, diffphase)

		def subpixel_template(sx, sy, dx, dy,
		s_img, d_img,
		@@ -223,16 +268,19 @@ def subpixel_template(sx, sy, dx, dy,
		image_size = check_image_size(image_size)
		template_size = check_image_size(template_size)

		s_image, _, _ = clip_roi(s_img, sx, sy, size_x=image_size[0], size_y=image_size[1])
		d_template, dxr, dyr = clip_roi(d_img, dx, dy, size_x=template_size[0], size_y=template_size[1])
		s_roi = roi.Roi(s_img, sx, sy, size_x=image_size[0], size_y=image_size[1])
		d_roi = roi.Roi(d_img, dx, dy, size_x=template_size[0], size_y=template_size[1])

		s_image = s_roi.clip()
		d_template = d_roi.clip()

		if (s_image is None) or (d_template is None):
		return None, None, None
		return None, None, None, None

		shift_x, shift_y, metrics, corrmap = func(d_template, s_image, **kwargs)

		dx = (dx - shift_x + dxr)
		dy = (dy - shift_y + dyr)
		dx = d_roi.x - shift_x
		dy = d_roi.y - shift_y

		return dx, dy, metrics, corrmap

		@@ -272,19 +320,22 @@ def subpixel_ciratefi(sx, sy, dx, dy, s_img, d_img, search_size=251, template_si
		strength : float
		Strength of the correspondence in the range [-1, 1]
		"""
		template, _, _ = clip_roi(d_img, dx, dy,
		t_roi = roi.Roi(d_img, dx, dy,
		size_x=template_size, size_y=template_size)
		search, dxr, dyr = clip_roi(s_img, sx, sy,
		s_roi = roi.Roi(s_img, sx, sy,
		size_x=search_size, size_y=search_size)
		template = t_roi.clip()
		search = s_roi.clip()

		if template is None or search is None:
		return None, None, None

		x_offset, y_offset, strength = ciratefi.ciratefi(template, search, **kwargs)
		dx += (x_offset + dxr)
		dy += (y_offset + dyr)
		dx += (x_offset + t_roi.axr)
		dy += (y_offset + t_roi.ayr)
		return dx, dy, strength

		def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=251, reduction=11, convergence_threshold=1.0, max_dist=50, **kwargs):
		def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=(251, 251), reduction=11, convergence_threshold=1.0, max_dist=50, **kwargs):
		"""
		Iteratively apply a subpixel phase matcher to source (s_img) and destination (d_img)
		images. The size parameter is used to set the initial search space. The algorithm
		@@ -307,10 +358,8 @@ def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=251, reduction=11, conver
		A plio geodata object from which the template is extracted
		d_img : object
		A plio geodata object from which the search is extracted
		size : int, tuple
		One half of the total size of the template, so a 251 default results in a 502 pixel search space.
		If an int, the template is square. If a tuple, in the form (x,y), is passed an
		irregularly shaped template can be used.
		size : tuple
		Size of the template in the form (x,y)
		reduction : int
		With each recursive call to this func, the size is reduced by this amount
		convergence_threshold : float
		@@ -334,50 +383,31 @@ def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=251, reduction=11, conver
		# get initial destination location
		dsample = dx
		dline = dy
		if isinstance(size, int):
		size = (size, size)
		while True:
		s_template, _, _ = clip_roi(s_img, sx, sy,
		size_x=size[0], size_y=size[1])
		d_search, dxr, dyr = clip_roi(d_img, dx, dy,
		size_x=size[0], size_y=size[1])

		if (s_template is None) or (d_search is None):
		return None, None, None
		if s_template.shape != d_search.shape:
		s_size = s_template.shape
		d_size = d_search.shape
		updated_size_x = int(min(s_size[1], d_size[1])) # Why is this /2?
		updated_size_y = int(min(s_size[0], d_size[0]))
		# Since the image is smaller than the requested size, set the size to
		# the current maximum image size and reduce from there on potential
		# future iterations.
		size = (updated_size_x, updated_size_y)
		s_template, _, _ = clip_roi(s_template, sx, sy,
		size_x=size[0], size_y=size[1])
		d_search, dxr, dyr = clip_roi(d_search, dx, dy,
		size_x=size[0], size_y=size[1])
		if (s_template is None) or (d_search is None):
		return None, None, None

		# Apply the phase matcher
		while True:
		try:
		shift_x, shift_y, metrics = subpixel_phase(s_template, d_search, **kwargs)
		shifted_dx, shifted_dy, metrics = subpixel_phase(sx, sy, dx, dy, s_img, d_img, image_size=size, **kwargs)
		except:
		return None, None, None
		# Apply the shift to d_search and compute the new correspondence location
		dx += shift_x # The implementation already applies the dxr, dyr shifts
		dy += shift_y


		# Compute the amount of move the matcher introduced
		delta_dx = abs(shifted_dx - dx)
		delta_dy = abs(shifted_dy - dy)
		dx = shifted_dx
		dy = shifted_dy

		# Break if the solution has converged
		size = (size[0] - reduction, size[1] - reduction)

		dist = np.linalg.norm([dsample-dx, dline-dy])

		if min(size) < 1:
		return None, None, None
		if abs(shift_x) <= convergence_threshold and\
		abs(shift_y) <= convergence_threshold and\
		if delta_dx <= convergence_threshold and\
		delta_dy<= convergence_threshold and\
		abs(dist) <= max_dist:
		break

		return dx, dy, metrics

		def subpixel_register_measure(measureid, iterative_phase_kwargs={}, subpixel_template_kwargs={},

autocnet/matcher/tests/test_ciratefi.py

+4 −4

Original line number	Diff line number	Diff line
		@@ -7,7 +7,7 @@ from imageio import imread
		from scipy.ndimage.interpolation import rotate

		from autocnet.examples import get_path
		from autocnet.matcher import subpixel as sp
		from autocnet.transformation import roi
		from .. import ciratefi

		import pytest
		@@ -32,7 +32,7 @@ def img_coord():
		@pytest.fixture
		def template(img, img_coord):
		coord_x, coord_y = img_coord
		template, _, _ = sp.clip_roi(img, coord_x, coord_y, 5, 5)
		template= roi.Roi(img, coord_x, coord_y, 5, 5).clip()
		template = rotate(template, 90)
		return template

		@@ -40,7 +40,7 @@ def template(img, img_coord):
		def search():
		coord_x, coord_y = (482.09783936, 652.40679932)
		img = imread(get_path('AS15-M-0298_SML.png'), as_gray=True)
		search, _, _ = sp.clip_roi(img, coord_x, coord_y, 21, 21)
		search = roi.Roi(img, coord_x, coord_y, 21, 21).clip()
		return search

		@pytest.fixture
		@@ -48,7 +48,7 @@ def offset_template(img, img_coord):
		coord_x, coord_y = img_coord
		coord_x += 1
		coord_y += 1
		offset_template, _, _ = sp.clip_roi(img, coord_x, coord_y, 5, 5)
		offset_template = roi.Roi(img, coord_x, coord_y, 5, 5).clip()
		return offset_template

		def test_cifi_radii_too_large(template, search):

autocnet/matcher/tests/test_subpixel.py

+97 −32

Original line number	Diff line number	Diff line
		@@ -36,16 +36,6 @@ def test_clip_roi(center_x, center_y, size, expected):

		assert clip.mean() == expected


		def test_subpixel_phase(apollo_subsets):
		a = apollo_subsets[0]
		b = apollo_subsets[1]

		xoff, yoff, err = sp.subpixel_phase(a, b)
		assert xoff == 0
		assert yoff == 2
		assert len(err) == 2

		def test_subpixel_template(apollo_subsets):
		def clip_side_effect(args, *kwargs):
		if np.array_equal(a, args[0]):
		@@ -66,34 +56,109 @@ def test_subpixel_template(apollo_subsets):
		a, b, upsampling=16)

		assert strength >= 0.99
		assert nx == 50.9375
		assert ny == 53.0625
		assert nx == 50.5
		assert ny == 52.4375

		@pytest.mark.parametrize("convergence_threshold, expected", [(1.0, (None, None, None)),
		(2.0, (50.49, 52.44, (0.039507, -9.5e-20)))])
		@pytest.mark.parametrize("convergence_threshold, expected", [(2.0, (50.49, 52.08, (0.039507, -9.5e-20)))])
		def test_iterative_phase(apollo_subsets, convergence_threshold, expected):
		def clip_side_effect(args, *kwargs):
		if np.array_equal(a, args[0]):
		return a, 0, 0
		else:
		return b, 0, 0
		a = apollo_subsets[0]
		b = apollo_subsets[1]
		with patch('autocnet.matcher.subpixel.clip_roi', side_effect=clip_side_effect):
		nx, ny, strength = sp.iterative_phase(a.shape[1]/2, a.shape[0]/2,
		dx, dy, strength = sp.iterative_phase(a.shape[1]/2, a.shape[0]/2,
		b.shape[1]/2, b.shape[1]/2,
		a, b, convergence_threshold=convergence_threshold,
		a, b,
		size=(51,51),
		convergence_threshold=convergence_threshold,
		upsample_factor=100)
		assert nx == expected[0]
		assert ny == expected[1]
		assert dx == expected[0]
		assert dy == expected[1]
		if expected[2] is not None:
		for i in range(len(strength)):
		assert pytest.approx(strength[i],6) == expected[2][i]

		@pytest.mark.parametrize("data, expected", [
		((21,21), (10.5, 10.5)),
		((20,20), (11,11)),
		((0,0), (1,1))
		((21,21), (10, 10)),
		((20,20), (10,10))
		])
		def test_check_image_size(data, expected):
		assert sp.check_image_size(data) == expected

		@pytest.mark.parametrize("x, y, x1, y1, image_size, template_size, expected",[
		(4, 3, 3, 2, (3,3), (3,3), (3,2)),
		(4, 3, 3, 2, (7,7), (3,3), (3,2)), # Increase the search image size
		(4, 3, 3, 2, (7,7), (5,5), (3,2)), # Increase the template size
		(4, 3, 2, 2, (7,7), (3,3), (3,2)), # Move point in the x-axis
		(4, 3, 4, 3, (7,7), (3,3), (3,2)), # Move point in the other x-direction
		(4, 3, 3, 1, (7,7), (3,3), (3,2)), # Move point negative in the y-axis
		(4, 3, 3, 3, (7,7), (3,3), (3,2)) # Move point positive in the y-axis

		])
		def test_subpixel_template_cooked(x, y, x1, y1, image_size, template_size, expected):
		test_image = np.array(((0, 0, 0, 0, 0, 0, 0, 1, 0),
		(0, 0, 0, 0, 0, 0, 0, 1, 0),
		(0, 0, 0, 1, 1, 1, 0, 1, 0),
		(0, 0, 0, 0, 1, 0, 0, 0, 0),
		(0, 0, 0, 0, 1, 0, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 1, 1, 1),
		(0, 1, 1, 1, 0, 0, 1, 0, 1),
		(0, 1, 0, 1, 0, 0, 1, 0, 1),
		(0, 1, 1, 1, 0, 0, 1, 0, 1),
		(0, 0, 0, 0, 0, 0, 1, 1, 1)), dtype=np.uint8)

		# Should yield (-3, 3) offset from image center
		t_shape = np.array(((0, 0, 0, 0, 0, 0, 0),
		(0, 0, 1, 1, 1, 0, 0),
		(0, 0, 0, 1, 0, 0, 0),
		(0, 0, 0, 1, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 0)), dtype=np.uint8)

		dx, dy, corr, corrmap = sp.subpixel_template(x, y, x1, y1,
		test_image, t_shape,
		image_size=image_size,
		template_size=template_size,
		upsampling=1)
		assert corr >= 1.0 # geq because sometime returning weird float > 1 from OpenCV
		assert dx == expected[0]
		assert dy == expected[1]

		@pytest.mark.parametrize("x, y, x1, y1, image_size, expected",[
		(4, 3, 3, 2, (3,3), (3,2)),
		(4, 3, 3, 2, (5,5), (3,2)), # Increase the search image size
		(4, 3, 3, 2, (5,5), (3,2)), # Increase the template size
		(4, 3, 2, 2, (5,5), (3,2)), # Move point in the x-axis
		(4, 3, 4, 3, (5,5), (3,2)), # Move point in the other x-direction
		(4, 3, 3, 1, (5,5), (3,2)), # Move point negative in the y-axis; also tests size reduction
		(4, 3, 3, 3, (5,5), (3,2)) # Move point positive in the y-axis

		])
		def test_subpixel_phase_cooked(x, y, x1, y1, image_size, expected):
		test_image = np.array(((0, 0, 0, 0, 0, 0, 0, 1, 0),
		(0, 0, 0, 0, 0, 0, 0, 1, 0),
		(0, 0, 0, 1, 1, 1, 0, 1, 0),
		(0, 0, 0, 0, 1, 0, 0, 0, 0),
		(0, 0, 0, 0, 1, 0, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 1, 1, 1),
		(0, 1, 1, 1, 0, 0, 1, 0, 1),
		(0, 1, 0, 1, 0, 0, 1, 0, 1),
		(0, 1, 1, 1, 0, 0, 1, 0, 1),
		(0, 0, 0, 0, 0, 0, 1, 1, 1)), dtype=np.uint8)

		# Should yield (-3, 3) offset from image center
		t_shape = np.array(((0, 0, 0, 0, 0, 0, 0),
		(0, 0, 1, 1, 1, 0, 0),
		(0, 0, 0, 1, 0, 0, 0),
		(0, 0, 0, 1, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 0),
		(0, 0, 0, 0, 0, 0, 0)), dtype=np.uint8)

		dx, dy, metrics = sp.subpixel_phase(x, y, x1, y1,
		test_image, t_shape,
		image_size=image_size)

		assert dx == expected[0]
		assert dy == expected[1]
		No newline at end of file

autocnet/spatial/overlap.py

+11 −7

Original line number	Diff line number	Diff line
		@@ -12,9 +12,9 @@ from autocnet import config, dem, Session
		from autocnet.cg import cg as compgeom
		from autocnet.io.db.model import Images, Measures, Overlay, Points, JsonEncoder
		from autocnet.spatial import isis
		from autocnet.matcher.subpixel import clip_roi
		from autocnet.matcher.cpu_extractor import extract_most_interesting
		from autocnet.transformation.spatial import reproject
		from autocnet.transformation import roi

		from plurmy import Slurm
		import csmapi
		@@ -163,6 +163,7 @@ def place_points_in_overlap(nodes, geom, cam_type="csm",
		points = []
		semi_major = config['spatial']['semimajor_rad']
		semi_minor = config['spatial']['semiminor_rad']

		valid = compgeom.distribute_points_in_geom(geom, **distribute_points_kwargs)
		if not valid:
		warnings.warn('Failed to distribute points in overlap')
		@@ -194,17 +195,20 @@ def place_points_in_overlap(nodes, geom, cam_type="csm",
		sample, line = image_coord.samp, image_coord.line

		# Extract ORB features in a sub-image around the desired point
		image, _, _ = clip_roi(node.geodata, sample, line, size_x=size, size_y=size)
		image_roi = roi.Roi(node.geodata, sample, line, size_x=size, size_y=size)
		image = image_roi.clip()
		try:
		interesting = extract_most_interesting(image)
		except:
		warnings.warn('Could not find an interesting feature around point')
		continue

		# kps are in the image space with upper left origin, so convert to
		# center origin and then convert back into full image space
		newsample = sample + (interesting.x - size)
		newline = line + (interesting.y - size)
		# kps are in the image space with upper left origin and the roi
		# could be the requested size or smaller if near an image boundary.
		# So use the roi upper left_x and top_y for the actual origin.
		left_x, _, top_y, _ = image_roi.image_extent
		newsample = left_x + interesting.x
		newline = top_y + interesting.y

		# Get the updated lat/lon from the feature in the node
		if cam_type == "isis":

autocnet/spatial/tests/test_overlap_points.py

+7 −5

Original line number	Diff line number	Diff line
		@@ -9,36 +9,39 @@ from autocnet.graph.node import Node
		import csmapi

		MockKeypoints = namedtuple('Keypoints', ['x', 'y'])
		mockkeypoints = MockKeypoints(0,0)
		mockkeypoints = MockKeypoints(5,5)

		@patch('autocnet.cg.cg.distribute_points_in_geom', return_value=[(0, 0), (5, 5), (10, 10)])
		@patch('autocnet.spatial.overlap.clip_roi', return_value=np.zeros((3,3)))
		@patch('autocnet.spatial.overlap.extract_most_interesting', return_value=mockkeypoints)
		def test_place_points_in_overlap(point_distributer, clip_roi, extractor):
		def test_place_points_in_overlap(point_distributer, extractor):
		# Mock setup
		first_node = MagicMock()
		first_node.camera = MagicMock()
		first_node.camera.groundToImage.return_value = csmapi.ImageCoord(1.0, 0.0)
		first_node.camera.imageToGround.return_value = csmapi.EcefCoord(1.0,1.0,1.0)
		first_node.isis_serial = '1'
		first_node.geodata = np.zeros((100,100))
		first_node.__getitem__.return_value = 1
		second_node = MagicMock()
		second_node.camera = MagicMock()
		second_node.camera.groundToImage.return_value = csmapi.ImageCoord(1.0, 1.0)
		second_node.camera.imagetoground.return_value = csmapi.EcefCoord(1.0,1.0,1.0)
		second_node.isis_serial = '2'
		second_node.geodata = np.zeros((100,100))
		second_node.__getitem__.return_value = 2
		third_node = MagicMock()
		third_node.camera = MagicMock()
		third_node.camera.groundToImage.return_value = csmapi.ImageCoord(0.0, 1.0)
		third_node.camera.imageToGround.return_value = csmapi.EcefCoord(1.0,1.0,1.0)
		third_node.isis_serial = '3'
		third_node.geodata = np.zeros((100,100))
		third_node.__getitem__.return_value = 3
		fourth_node = MagicMock()
		fourth_node.camera = MagicMock()
		fourth_node.camera.groundToImage.return_value = csmapi.ImageCoord(0.0, 0.0)
		first_node.camera.imageToGround.return_value = csmapi.EcefCoord(1.0,0,1.0)
		fourth_node.isis_serial = '4'
		fourth_node.geodata = np.zeros((100,100))
		fourth_node.__getitem__.return_value = 4

		# Actual function being tested
		@@ -55,7 +58,6 @@ def test_place_points_in_overlap(point_distributer, clip_roi, extractor):
		assert measure_serials == ['1', '2', '3', '4']

		# Check the mocks
		np.testing.assert_array_equal(point_distributer.call_args[0][0], np.array([0.0, 0.0, 0.0]))
		first_node.camera.groundToImage.assert_called()
		second_node.camera.groundToImage.assert_called()
		third_node.camera.groundToImage.assert_called()