Test Updates and Iterative Phase fix (32b55a4b) · Commits · aflab / astrogeology / Autocnet

autocnet/matcher/subpixel.py

+37 −40

Original line number	Diff line number	Diff line
		@@ -100,7 +100,7 @@ def subpixel_phase(reference_roi, moving_roi, affine=tf.AffineTransform(), **kwa
		"""
		reference_roi.clip()
		reference_image = reference_roi.clipped_array
		moving_roi.clip(affine=affine)
		moving_roi.clip(affine=affine, warp_mode="constant", coord_mode="constant")
		walking_template = moving_roi.clipped_array

		if reference_image.shape != walking_template.shape:
		@@ -125,8 +125,11 @@ def subpixel_phase(reference_roi, moving_roi, affine=tf.AffineTransform(), **kwa
		reference_roi.size_x, reference_roi.size_y = size
		moving_roi.size_x, moving_roi.size_y = size

		reference_image = reference_roi.clip()
		walking_template = moving_roi.clip(affine)
		reference_roi.clip()
		reference_image = reference_roi.clipped_array
		moving_roi.clip(affine=affine)
		walking_template = moving_roi.clipped_array

		(shift_y, shift_x), error, diffphase = registration.phase_cross_correlation(reference_image,
		walking_template,
		**kwargs)
		@@ -236,8 +239,6 @@ def iterative_phase(reference_roi, moving_roi, affine=tf.AffineTransform(), redu
		moving_roi : Object
		An Roi object from autocnet, the walking image to move around and make comparisons to
		the reference roi. Contains either an ndarray or a GeoDataset Object
		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
		@@ -267,7 +268,7 @@ def iterative_phase(reference_roi, moving_roi, affine=tf.AffineTransform(), redu
		# Compute the amount of move the matcher introduced
		delta_dx, delta_dy = abs(new_subpixel_affine.translation)
		subpixel_affine += new_subpixel_affine
		dx, dy = subpixel_affine.inverse((reference_roi.x, reference_roi.y))[0]
		dx, dy = subpixel_affine.inverse((moving_roi.x, moving_roi.y))[0]
		moving_roi.x, moving_roi.y = dx, dy

		# Break if the solution has converged
		@@ -793,43 +794,39 @@ def fourier_mellen(ref_image, moving_image, affine=tf.AffineTransform(), verbose
		Error returned by the iterative phase matcher
		"""
		# Get the affine transformation for scale + rotation invariance
		affine = estimate_logpolar_transform(ref_image.array, moving_image.array, verbose=verbose)

		# warp the source image to match the destination
		ref_warped = ref_image.clip(affine)
		sx, sy = affine.inverse(np.asarray(ref_image.array.shape)/2)[0]
		affine = estimate_logpolar_transform(ref_image.clipped_array, moving_image.clipped_array, verbose=verbose)

		# get translation with iterative phase
		newx, newy, error = iterative_phase(sx, sy, sx, sy, ref_warped, moving_image, **phase_kwargs)
		subpixel_affine, error, diffphase = iterative_phase(ref_image, moving_image, affine=affine, **phase_kwargs)

		# if verbose:
		# fig, axes = plt.subplots(2, 2, figsize=(8, 8))
		# ax = axes.ravel()
		#
		# ax[0].imshow(ref_warped)
		# ax[0].set_title("Image 1 Transformed")
		# ax[0].axhline(y=sy, color="red", linestyle="-", alpha=1, linewidth=1)
		# ax[0].axvline(x=sx, color="red", linestyle="-", alpha=1, linewidth=1)
		#
		# ax[2].imshow(ref_image)
		# ax[2].set_title("Image 1")
		# ax[2].axhline(y=ref_image.array.shape[0]/2, color="red", linestyle="-", alpha=1, linewidth=1)
		# ax[2].axvline(x=ref_image.array.shape[1]/2, color="red", linestyle="-", alpha=1, linewidth=1)
		#
		# ax[1].imshow(moving_image)
		# ax[3].imshow(moving_image)
		#
		# if not newx or not newy:
		# ax[1].set_title("Image 2 REGISTRATION FAILED")
		# ax[3].set_title("Image 2 REGISTRATION FAILED")
		# else :
		# ax[3].set_title("Image 2 Registered")
		# ax[1].axhline(y=newy, color="red", linestyle="-", alpha=1, linewidth=1)
		# ax[1].axvline(x=newx, color="red", linestyle="-", alpha=1, linewidth=1)
		# ax[3].axhline(y=newy, color="red", linestyle="-", alpha=1, linewidth=1)
		# ax[3].axvline(x=newx, color="red", linestyle="-", alpha=1, linewidth=1)

		if verbose:
		fig, axes = plt.subplots(2, 2, figsize=(8, 8))
		ax = axes.ravel()

		ax[0].imshow(ref_warped)
		ax[0].set_title("Image 1 Transformed")
		ax[0].axhline(y=sy, color="red", linestyle="-", alpha=1, linewidth=1)
		ax[0].axvline(x=sx, color="red", linestyle="-", alpha=1, linewidth=1)

		ax[2].imshow(ref_image)
		ax[2].set_title("Image 1")
		ax[2].axhline(y=ref_image.array.shape[0]/2, color="red", linestyle="-", alpha=1, linewidth=1)
		ax[2].axvline(x=ref_image.array.shape[1]/2, color="red", linestyle="-", alpha=1, linewidth=1)

		ax[1].imshow(moving_image)
		ax[3].imshow(moving_image)

		if not newx or not newy:
		ax[1].set_title("Image 2 REGISTRATION FAILED")
		ax[3].set_title("Image 2 REGISTRATION FAILED")
		else :
		ax[3].set_title("Image 2 Registered")
		ax[1].axhline(y=newy, color="red", linestyle="-", alpha=1, linewidth=1)
		ax[1].axvline(x=newx, color="red", linestyle="-", alpha=1, linewidth=1)
		ax[3].axhline(y=newy, color="red", linestyle="-", alpha=1, linewidth=1)
		ax[3].axvline(x=newx, color="red", linestyle="-", alpha=1, linewidth=1)

		return newx, newy, error
		return subpixel_affine, error, diffphase

		def subpixel_register_point_smart(pointid,
		cost_func=lambda x,y: 1/x*2 y,

autocnet/matcher/tests/test_subpixel.py

+76 −111

Original line number	Diff line number	Diff line
		@@ -2,21 +2,21 @@ import math
		import os
		import sys
		import unittest
		from unittest.mock import patch, Mock
		from unittest.mock import patch, MagicMock, Mock, PropertyMock
		import logging


		from skimage import transform as tf
		from skimage.util import img_as_float
		from skimage import color
		from skimage import data

		import pytest
		import tempfile

		import numpy as np
		from imageio import imread

		from plio.io.io_gdal import GeoDataset
		from plio.io.io_gdal import GeoDataset, array_to_raster

		from autocnet.examples import get_path
		import autocnet.matcher.subpixel as sp
		@@ -35,48 +35,48 @@ def iris_pair():
		rescaled = tf.rescale(rotated, scale)
		sizer, sizec = image.shape
		rts_image = rescaled[:sizer, :sizec]
		return image, rts_image

		@pytest.fixture
		def apollo_subsets():
		# These need to be geodata sets or just use mocks...
		arr1 = imread(get_path('AS15-M-0295_SML(1).png'))[100:201, 123:224]
		arr2 = imread(get_path('AS15-M-0295_SML(2).png'))[235:336, 95:196]
		return arr1, arr2

		roi_raster1 = tempfile.NamedTemporaryFile()
		roi_raster2 = tempfile.NamedTemporaryFile()

		@pytest.mark.parametrize("center_x, center_y, size, expected", [(4, 4, 9, 404),
		(55.4, 63.1, 27, 6355)])
		def test_clip_roi(center_x, center_y, size, expected):
		img = np.arange(10000).reshape(100, 100)
		array_to_raster(image, roi_raster1.name)
		array_to_raster(rts_image, roi_raster2.name)

		clip, axr, ayr = sp.clip_roi(img, center_x, center_y, size)
		roi1 = roi.Roi(GeoDataset(roi_raster1.name), x=705, y=705, size_x=50, size_y=50)
		roi2 = roi.Roi(GeoDataset(roi_raster2.name), x=705, y=705, size_x=50, size_y=50)
		return roi1, roi2

		assert clip.mean() == expected

		def clip_side_effect(args, *kwargs):
		if np.array_equal(a, args[0]):
		return a, 0, 0
		def clip_side_effect(arr, clip=False):
		if not clip:
		return arr
		else:
		center_y = b.shape[0] / 2
		center_x = b.shape[1] / 2
		bxr, bx = math.modf(center_x)
		byr, by = math.modf(center_y)
		bx = int(bx)
		by = int(by)
		return b[by-10:by+11, bx-10:bx+11], bxr, byr
		center_y = arr.shape[0] / 2
		center_x = arr.shape[1] / 2
		xr, x = math.modf(center_x)
		yr, y = math.modf(center_y)
		x = int(x)
		y = int(y)
		return arr[y-10:y+11, x-10:x+11]

		@pytest.fixture
		def apollo_subsets():
		roi1 = roi.Roi(GeoDataset(get_path('AS15-M-0295_SML(1).png')), x=173, y=150, size_x=50, size_y=50)
		roi2 = roi.Roi(GeoDataset(get_path('AS15-M-0295_SML(2).png')), x=145, y=285, size_x=50, size_y=50)
		roi1.clip()
		roi2.clip()
		return roi1, roi2

		def test_subpixel_template(apollo_subsets):
		a = apollo_subsets[0]
		b = apollo_subsets[1]
		with patch('autocnet.matcher.subpixel.clip_roi', side_effect=clip_side_effect):
		nx, ny, strength, _ = sp.subpixel_template(a.shape[1]/2, a.shape[0]/2,
		b.shape[1]/2, b.shape[0]/2,
		a, b, upsampling=16)

		assert strength >= 0.99
		assert nx == 50.5
		assert ny == 52.4375
		b.size_x = 10
		b.size_y = 10
		affine, metrics, corr_map = sp.subpixel_template(a, b, upsampling=16)
		nx, ny = affine.translation
		assert nx == -0.3125
		assert ny == 1.5
		assert np.max(corr_map) >= 0.9367293
		assert metrics >= 0.9367293

		@pytest.mark.parametrize("loc, failure", [((0,4), True),
		((4,0), True),
		@@ -84,44 +84,33 @@ def test_subpixel_template(apollo_subsets):
		def test_subpixel_template_at_edge(apollo_subsets, loc, failure):
		a = apollo_subsets[0]
		b = apollo_subsets[1]
		b.size_x = 10
		b.size_y = 10

		def func(args, *kwargs):
		corr = np.zeros((10,10))
		corr[loc[0], loc[1]] = 10
		return 0, 0, 0, corr

		with patch('autocnet.matcher.subpixel.clip_roi', side_effect=clip_side_effect):
		if failure:
		# with pytest.warns(UserWarning, match=r'Maximum correlation \S+'):
		# match=r'Maximum correlation \S+'
		nx, ny, strength, _ = sp.subpixel_template(a.shape[1]/2, a.shape[0]/2,
		b.shape[1]/2, b.shape[0]/2,
		a, b, upsampling=16,
		affine, metrics, corr_map = sp.subpixel_template(a, b, upsampling=16,
		func=func)
		else:
		nx, ny, strength, _ = sp.subpixel_template(a.shape[1]/2, a.shape[0]/2,
		b.shape[1]/2, b.shape[0]/2,
		a, b, upsampling=16,
		affine, metrics, corr_map = sp.subpixel_template(a, b, upsampling=16,
		func=func)
		assert nx == 50.5

		def test_subpixel_transformed_template(apollo_subsets):
		a = apollo_subsets[0]
		b = apollo_subsets[1]
		transform = tf.AffineTransform(rotation=math.radians(1), scale=(1.1,1.1))
		with patch('autocnet.matcher.subpixel.clip_roi', side_effect=clip_side_effect):
		nx, ny, strength, _ = sp.subpixel_transformed_template(a.shape[1]/2, a.shape[0]/2,
		b.shape[1]/2, b.shape[0]/2,
		a, b, transform, upsampling=16)

		assert strength >= 0.83
		assert nx == pytest.approx(50.576284)
		assert ny == pytest.approx(54.0081)
		nx, ny = affine.translation
		assert nx == 0


		def test_estimate_logpolar_transform(iris_pair):
		img1, img2 = iris_pair
		affine = sp.estimate_logpolar_transform(img1, img2)
		roi1, roi2 = iris_pair
		roi1.size_x = 705
		roi1.size_y = 705
		roi2.size_x = 705
		roi2.size_y = 705
		roi1.clip()
		roi2.clip()
		affine = sp.estimate_logpolar_transform(roi1.clipped_array, roi2.clipped_array)

		assert pytest.approx(affine.scale, 0.1) == 0.71
		assert pytest.approx(affine.rotation, 0.1) == 0.34
		@@ -130,57 +119,33 @@ def test_estimate_logpolar_transform(iris_pair):


		def test_fourier_mellen(iris_pair):
		img1, img2 = iris_pair
		nx, ny, error = sp.fourier_mellen(img1, img2, phase_kwargs = {"reduction" : 11, "size":(401, 401), "convergence_threshold" : 1, "max_dist":100})
		roi1, roi2 = iris_pair
		roi1.size_x = 200
		roi1.size_y = 200
		roi2.size_x = 200
		roi2.size_y = 200
		roi1.clip()
		roi2.clip()
		affine, metrics, corrmap = sp.fourier_mellen(roi1, roi2, phase_kwargs = {"reduction" : 11, "convergence_threshold" : 1, "max_dist":100})

		assert pytest.approx(nx, 0.01) == 996.39
		assert pytest.approx(ny, 0.01) == 984.912
		assert pytest.approx(error, 0.01) == 0.0422


		@pytest.mark.parametrize("loc, failure", [((0,4), True),
		((4,0), True),
		((1,1), False)])
		def test_subpixel_transformed_template_at_edge(apollo_subsets, loc, failure):
		a = apollo_subsets[0]
		b = apollo_subsets[1]

		def func(args, *kwargs):
		corr = np.zeros((5,5))
		corr[loc[0], loc[1]] = 10
		return 0, 0, 0, corr

		transform = tf.AffineTransform(rotation=math.radians(1), scale=(1.1,1.1))
		with patch('autocnet.matcher.subpixel.clip_roi', side_effect=clip_side_effect):
		if failure:
		print(a.shape[1]/2, a.shape[0]/2,b.shape[1]/2, b.shape[0]/2,
		a, b)
		nx, ny, strength, _ = sp.subpixel_transformed_template(a.shape[1]/2, a.shape[0]/2,
		b.shape[1]/2, b.shape[0]/2,
		a, b, transform, upsampling=16,
		func=func)
		else:
		nx, ny, strength, _ = sp.subpixel_transformed_template(a.shape[1]/2, a.shape[0]/2,
		b.shape[1]/2, b.shape[0]/2,
		a, b, transform, upsampling=16,
		func=func)
		assert nx == 50.5

		@pytest.mark.parametrize("convergence_threshold, expected", [(2.0, (50.49, 52.08, -9.5e-20))])
		@pytest.mark.parametrize("convergence_threshold, expected", [(2.0, (-0.32, 1.66, -9.5e-20))])
		def test_iterative_phase(apollo_subsets, convergence_threshold, expected):
		a = apollo_subsets[0]
		b = apollo_subsets[1]
		dx, dy, strength = sp.iterative_phase(a.shape[1]/2, a.shape[0]/2,
		b.shape[1]/2, b.shape[1]/2,
		a, b,
		size=(51,51),
		affine, metrics, corr_map = sp.iterative_phase(a, b,
		convergence_threshold=convergence_threshold,
		upsample_factor=100)
		dx, dy = affine.translation
		assert dx == expected[0]
		assert dy == expected[1]
		if expected[2] is not None:
		# for i in range(len(strength)):
		assert pytest.approx(strength,6) == expected[2]
		assert pytest.approx(metrics,6) == expected[2]

		@pytest.mark.parametrize("data, expected", [
		((21,21), (10, 10)),

autocnet/transformation/roi.py

+5 −15

Original line number	Diff line number	Diff line
		@@ -74,20 +74,10 @@ class Roi():

		@property
		def clip_center(self):
		if not getattr(self, '_clip_center', None):
		if not hasattr(self, '_clip_center'):
		self.clip()
		return self._clip_center

		@property
		def clipped_array(self):
		"""
		The clipped array associated with this ROI.
		"""
		if not hasattr(self, "_clipped_array"):
		self.clip()
		return self._clipped_array


		@property
		def affine(self):
		return self._affine
		@@ -236,7 +226,7 @@ class Roi():

		return x_in_image_space, y_in_image_space

		def clip(self, size_x=None, size_y=None, affine=None, dtype=None, mode="reflect"):
		def clip(self, size_x=None, size_y=None, affine=None, dtype=None, warp_mode="reflect", coord_mode="reflect"):
		"""
		Compatibility function that makes a call to the array property.
		Warning: The dtype passed in via this function resets the dtype attribute of this
		@@ -299,7 +289,7 @@ class Roi():
		warped_data = tf.warp(data,
		self.affine,
		order=3,
		mode='reflect')
		mode=warp_mode)

		self.warped_array_center = self.affine.inverse(data_center)[0]

		@@ -318,7 +308,7 @@ class Roi():
		# the xi, yi are intentionally handed in backward, because the map_coordinates indexes column major
		pixel_locked = ndimage.map_coordinates(warped_data,
		np.meshgrid(yi, xi, indexing='ij'),
		mode=mode,
		mode=coord_mode,
		order=3)

		self._clip_center = (np.array(pixel_locked.shape)[::-1]) / 2.0
		@@ -337,7 +327,7 @@ class Roi():
		# the xi, yi are intentionally handed in backward, because the map_coordinates indexes column major
		pixel_locked = ndimage.map_coordinates(data,
		np.meshgrid(yi, xi, indexing='ij'),
		mode=mode,
		mode=coord_mode,
		order=3)

		if self.buffer != 0: