merge conflict (8fcbafa2) · Commits · aflab / astrogeology / Autocnet

autocnet/matcher/subpixel.py

+162 −174

Original line number	Diff line number	Diff line
		@@ -36,7 +36,7 @@ from autocnet.spatial import isis
		from autocnet.io.db.model import Measures, Points, Images, JsonEncoder
		from autocnet.graph.node import NetworkNode
		from autocnet.transformation import roi
		from autocnet.transformation.affine import estimate_affine_transformation
		from autocnet.transformation.affine import estimate_affine_from_sensors
		from autocnet import spatial
		from autocnet.utils.utils import bytescale

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

		def subpixel_phase(sx, sy, dx, dy,
		s_img, d_img,
		image_size=(51, 51),
		**kwargs):
		def subpixel_phase(reference_roi, walking_roi, affine=None, **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
		@@ -210,37 +207,33 @@ def subpixel_phase(sx, sy, dx, dy,

		Parameters
		----------
		template : ndarray
		The template used to search

		search : ndarray
		The search image
		reference_roi : Object
		An Roi object from autocnet, the reference image to use when computing subpixel offsets.
		Contains either an ndarray or a GeoDataset Object
		walking_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
		affine : Object
		Scikit image affine tranformation. This affine transformation is applied to the walking_roi
		before any comparisons are made

		Returns
		-------
		x_offset : float
		Shift in the x-dimension

		y_offset : float
		Shift in the y-dimension

		strength : tuple
		new_affine : Object
		Scikit image affine transformation. An updated affine transformation that is the new
		translation from the walking_roi to the reference_roi
		: tuple
		With the RMSE error and absolute difference in phase
		"""
		image_size = check_image_size(image_size)
		reference_image = reference_roi.clip()
		walking_template = walking_roi.clip(affine)

		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.array
		d_template = d_roi.array

		if s_image.shape != d_template.shape:
		if reference_image.shape != walking_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]))
		reference_size = reference_image.shape
		walking_size = walking_template.shape
		updated_size_x = int(min(walking_size[1], reference_size[1]))
		updated_size_y = int(min(walking_size[0], reference_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
		@@ -254,21 +247,18 @@ def subpixel_phase(sx, sy, dx, dy,
		# 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.array
		d_template = d_roi.array

		if (s_image is None) or (d_template is None):
		return None, None, None
		reference_roi.size_x, reference_roi.size_y = size
		walking_roi.size_x, walking_roi.size_y = size

		(shift_y, shift_x), error, diffphase = registration.phase_cross_correlation(s_image, d_template, **kwargs)
		dx = d_roi.x - shift_x
		dy = d_roi.y - shift_y
		reference_image = reference_roi.clip()
		walking_template = walking_roi.clip(affine)

		return dx, dy, error, None
		(shift_y, shift_x), error, diffphase = registration.phase_cross_correlation(reference_image,
		walking_template,
		**kwargs)
		new_affine = tf.AffineTransform(translation=(shift_x, shift_y))
		return new_affine, error, diffphase

		def subpixel_transformed_template(sx, sy, dx, dy,
		s_img, d_img,
		@@ -678,7 +668,7 @@ def subpixel_ciratefi(sx, sy, dx, dy, s_img, d_img, search_size=251, template_si
		dy += (y_offset + t_roi.ayr)
		return dx, dy, strength

		def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=(51, 51), reduction=11, convergence_threshold=1.0, max_dist=50, **kwargs):
		def iterative_phase(reference_roi, walking_roi, affine=None, reduction=11, convergence_threshold=0.1, 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
		@@ -689,18 +679,12 @@ def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=(51, 51), reduction=11, c

		Parameters
		----------
		sx : numeric
		The x position of the center of the template to be matched to
		sy : numeric
		The y position of the center of the template to be matched to
		dx : numeric
		The x position of the center of the search to be matched from
		dy : numeric
		The y position of the center of the search to be matched to
		s_img : object
		A plio geodata object from which the template is extracted
		d_img : object
		A plio geodata object from which the search is extracted
		reference_roi : Object
		An Roi object from autocnet, the reference image to use when computing subpixel offsets.
		Contains either an ndarray or a GeoDataset Object
		walking_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
		@@ -710,10 +694,6 @@ def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=(51, 51), reduction=11, c

		Returns
		-------
		dx : float
		The new x value for the match in the destination (d) image
		dy : float
		The new y value for the match in the destination (d) image
		metrics : tuple
		A tuple of metrics. In the case of the phase matcher this are difference
		and RMSE in the phase dimension.
		@@ -722,32 +702,40 @@ def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=(51, 51), reduction=11, c
		--------
		subpixel_phase : the function that applies a single iteration of the phase matcher
		"""

		# get initial destination location
		dsample = dx
		dline = dy
		dsample, dline = walking_roi.x, walking_roi.y
		dx, dy = walking_roi.x, walking_roi.y

		while True:
		shifted_dx, shifted_dy, metrics, _ = subpixel_phase(sx, sy, dx, dy, s_img, d_img, image_size=size, **kwargs)
		size = (walking_roi.size_x, walking_roi.size_y)
		original_size = (walking_roi.size_x, walking_roi.size_y)

		subpixel_affine = tf.AffineTransform()

		while True:
		new_subpixel_affine, error, diffphase = subpixel_phase(reference_roi, walking_roi, affine, **kwargs)
		# 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
		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]
		walking_roi.x, walking_roi.y = dx, dy

		# Break if the solution has converged
		size = (size[0] - reduction, size[1] - reduction)
		walking_roi.size_x, walking_roi.size_y = size
		reference_roi.size_x, reference_roi.size_y = size
		dist = np.linalg.norm([dsample-dx, dline-dy])

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

		return dx, dy, metrics
		walking_roi.size_x, walking_roi.size_y = original_size
		reference_roi.size_x, reference_roi.size_y = original_size
		walking_roi.x, walking_roi.y = dsample, dline
		return subpixel_affine, error, diffphase

		'''def estimate_affine_transformation(destination_coordinates, source_coordinates):
		"""
		@@ -846,7 +834,7 @@ def geom_match_simple(reference_image,
		match_func = check_match_func(match_func)

		# Estimate the transformation between the two images
		affine = estimate_affine_transformation(reference_image, moving_image, bcenter_x, bcenter_y)
		affine = estimate_affine_from_sensors(reference_image, moving_image, bcenter_x, bcenter_y)
		t2 = time.time()
		print(f'Estimation of the transformation took {t2-t1} seconds.')

		@@ -1230,7 +1218,7 @@ def geom_match(destination_cube,


		# Estimate the transformation
		affine = estimate_affine_transformation(destination_corners, source_corners)
		affine = estimate_affine_from_sensors(destination_corners, source_corners)

		# Apply the subpixel matcher with an affine transformation
		restemplate = subpixel_transformed_template(bcenter_x, bcenter_y,
		@@ -2083,7 +2071,7 @@ def subpixel_register_point_smart(pointid,
		print('geom_func', geom_func)

		try:
		affine = estimate_affine_transformation(source_node.geodata,
		affine = estimate_affine_from_sensors(source_node.geodata,
		destination_node.geodata,
		source.apriorisample,
		source.aprioriline)
		@@ -2394,7 +2382,7 @@ def validate_candidate_measure(measure_to_register,

		print(f'Validating measure: {measure_to_register_id} on image: {source_image.name}')
		try:
		affine = estimate_affine_transformation(source_node.geodata, destination_node.geodata, sample, line)
		affine = estimate_affine_from_sensors(source_node.geodata, destination_node.geodata, sample, line)
		except:
		print('Unable to transform image to reference space. Likely too close to the edge of the non-reference image. Setting ignore=True')
		return [np.inf] * len(parameters)

autocnet/matcher/tests/test_subpixel.py

+31 −25

Original line number	Diff line number	Diff line
		@@ -16,6 +16,7 @@ from imageio import imread

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

		@pytest.fixture
		def iris_pair():
		@@ -167,21 +168,26 @@ def test_subpixel_transformed_template_at_edge(apollo_subsets, loc, failure):
		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, (50.49, 52.44, -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),
		reference_image = apollo_subsets[0]
		walking_image = apollo_subsets[1]
		image_size = (51, 51)

		ref_x, ref_y = reference_image.shape[0]/2, reference_image.shape[1]/2
		walk_x, walk_y = walking_image.shape[0]/2, walking_image.shape[1]/2

		reference_roi = roi.Roi(reference_image, ref_x, ref_y, size_x=image_size[0], size_y=image_size[1])
		walking_roi = roi.Roi(walking_image, walk_x, walk_y, size_x=image_size[0], size_y=image_size[1])
		affine, error, diffphase = sp.iterative_phase(reference_roi,
		walking_roi,
		convergence_threshold=convergence_threshold,
		upsample_factor=100)
		assert dx == expected[0]
		assert dy == expected[1]
		ref_to_walk = affine.inverse((ref_x, ref_y))[0]
		assert ref_to_walk[0] == expected[0]
		assert ref_to_walk[1] == expected[1]
		if expected[2] is not None:
		# for i in range(len(strength)):
		assert pytest.approx(strength,6) == expected[2]
		assert pytest.approx(error,6) == expected[2]

		@pytest.mark.parametrize("data, expected", [
		((21,21), (10, 10)),
		@@ -232,13 +238,13 @@ def test_subpixel_template_cooked(x, y, x1, y1, image_size, template_size, expec
		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
		(4, 3, 3, 2, (1,1), (3,2)),
		(4, 3, 3, 2, (2,2), (3,2)), # Increase the search image size
		(4, 3, 3, 2, (2,2), (3,2)), # Increase the template size
		(4, 3, 2, 2, (2,2), (3,2)), # Move point in the x-axis
		(4, 3, 4, 3, (2,2), (3,2)), # Move point in the other x-direction
		(4, 3, 3, 1, (2,2), (3,2)), # Move point negative in the y-axis; also tests size reduction
		(4, 3, 3, 3, (2,2), (3,2)) # Move point positive in the y-axis

		])
		def test_subpixel_phase_cooked(x, y, x1, y1, image_size, expected):
		@@ -264,11 +270,11 @@ def test_subpixel_phase_cooked(x, y, x1, y1, image_size, expected):
		(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)
		reference_roi = roi.Roi(test_image, x, y, size_x=image_size[0], size_y=image_size[1])
		walking_roi = roi.Roi(t_shape, x1, y1, size_x=image_size[0], size_y=image_size[1])

		affine, metrics, _ = sp.subpixel_phase(reference_roi, walking_roi)

		dx, dy = affine.inverse((x1, y1))[0]
		assert dx == expected[0]
		assert dy == expected[1]

autocnet/transformation/affine.py

+7 −7

Original line number	Diff line number	Diff line
		@@ -8,7 +8,7 @@ from autocnet.spatial import isis

		log = logging.getLogger(__name__)

		def estimate_affine_transformation(reference_image,
		def estimate_affine_from_sensors(reference_image,
		moving_image,
		bcenter_x,
		bcenter_y,
		@@ -85,5 +85,5 @@ def estimate_affine_transformation(reference_image,

		affine = tf.estimate_transform('affine', np.array([base_gcps]), np.array([dst_gcps]))
		t2 = time.time()
		print(f'Estimation of the transformation took {t2-t1} seconds.')
		log.debug(f'Estimation of the transformation took {t2-t1} seconds.')
		return affine
		No newline at end of file

autocnet/transformation/roi.py

+31 −8

Original line number	Diff line number	Diff line
		from math import modf, floor
		import numpy as np

		import pvl
		from skimage import transform as tf

		class Roi():
		"""
		@@ -126,13 +126,13 @@ class Roi():
		top_y = self._y - self.size_y
		bottom_y = self._y + self.size_y

		if self._x - self.size_x < 0:
		if left_x < 0:
		left_x = 0
		if self._y - self.size_y < 0:
		if top_y < 0:
		top_y = 0
		if self._x + self.size_x > raster_size[0]:
		if right_x > raster_size[0]:
		right_x = raster_size[0]
		if self._y + self.size_y > raster_size[1]:
		if bottom_y > raster_size[1]:
		bottom_y = raster_size[1]

		return list(map(int, [left_x, right_x, top_y, bottom_y]))
		@@ -170,7 +170,7 @@ class Roi():
		data = self.data.read_array(pixels=pixels)
		return data

		def clip(self, dtype=None):
		def clip(self, affine=None, dtype=None, 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
		@@ -186,5 +186,28 @@ class Roi():
		: ndarray
		The array attribute of this object.
		"""
		#self.dtype = dtype
		self.dtype = dtype

		if affine:
		self.size_x *= 2
		self.size_y *= 2

		array_to_warp = self.array

		# if array_to_warp.shape != ((self.size_y * 2) + 1, (self.size_x * 2) + 1):
		# raise ValueError("Unable to enlarge Roi to apply affine transformation." +
		# f" Was only able to extract {array_to_warp.shape}, when " +
		# f"{((self.size_y * 2) + 1, (self.size_x * 2) + 1)} was asked for. Select, " +
		# "a smaller region of interest" )


		#Affine transform the larger, moving array
		transformed_array = tf.warp(array_to_warp, affine, order=3, mode=mode)
		# print(transformed_array.shape)
		self.size_x = int(self.size_x / 2)
		self.size_y = int(self.size_y / 2)

		new_center = np.array(transformed_array.shape)/2
		return Roi(transformed_array, *new_center, self.size_x, self.size_y).clip()

		return self.array

autocnet/transformation/tests/test_affine.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -6,7 +6,7 @@ from autocnet.transformation import affine
		def test_estimate_affine_transformation():
		gd_base = GeoDataset('tests/test_subpixel_match/B08_012650_1780_XN_02S046W.l1.cal.destriped.crop.cub')
		gd_match = GeoDataset('tests/test_subpixel_match/J04_046447_1777_XI_02S046W.l1.cal.destriped.crop.cub')
		affine_transform = affine.estimate_affine_transformation(gd_base,gd_match, 150, 150)
		affine_transform = affine.estimate_affine_from_sensors(gd_base, gd_match, 150, 150)
		assert affine_transform.rotation == pytest.approx(-0.0012609633370663982, 6)
		assert affine_transform.shear == pytest.approx(0.01262411827876344)
		assert affine_transform.scale[0] == pytest.approx(1.0, 6)