Updates for correct ROI center, roi getting dtype via plio (17a5cb0f) · Commits · aflab / astrogeology / Autocnet

autocnet/matcher/naive_template.py

+9 −5

Original line number	Diff line number	Diff line
		@@ -107,8 +107,12 @@ def pattern_match(template, image, upsampling=16, metric=cv2.TM_CCOEFF_NORMED, e
		The x offset
		y : float
		The y offset
		strength : float
		max_corr : float
		The strength of the correlation in the range [-1, 1].
		result : ndarray
		(m,n) correlation matrix showing the correlation for all tested coordinates. The
		maximum correlation is reported when the upper left hand corner of the template
		maximally correlates with the image.
		"""
		if upsampling < 1:
		raise ValueError
		@@ -130,12 +134,12 @@ def pattern_match(template, image, upsampling=16, metric=cv2.TM_CCOEFF_NORMED, e
		x, y = (max_loc[0], max_loc[1])

		# Compute the idealized shift (image center)
		ideal_y = u_image.shape[0] / 2
		ideal_x = u_image.shape[1] / 2
		ideal_y = u_image.shape[0] / 2.
		ideal_x = u_image.shape[1] / 2.

		# Compute the shift from template upper left to template center
		y += (u_template.shape[0] / 2)
		x += (u_template.shape[1] / 2)
		y += (u_template.shape[0] / 2.)
		x += (u_template.shape[1] / 2.)

		x = (x - ideal_x) / upsampling
		y = (y - ideal_y) / upsampling

autocnet/matcher/subpixel.py

+35 −131

Original line number	Diff line number	Diff line
		@@ -145,63 +145,6 @@ def check_image_size(imagesize):
		y = floor(y/2)
		return x,y

		def clip_roi(img, center_x, center_y, size_x=200, size_y=200, dtype="uint64"):
		"""
		Given an input image, clip a square region of interest
		centered on some pixel at some size.

		Parameters
		----------
		img : ndarray or object
		The input image to be clipped or an object
		with a read_array method that takes a pixels
		argument in the form [xstart, ystart, xstop, ystop]

		center_x : Numeric
		The x coordinate to the center of the roi

		center_y : Numeric
		The y coordinate to the center of the roi

		img_size : int
		1/2 of the total image size. This value is the
		number of pixels grabbed from each side of the center

		Returns
		-------
		clipped_img : ndarray
		The clipped image
		"""

		try:
		raster_size = img.raster_size
		except:
		# x,y form
		raster_size = img.shape[::-1]
		axr, ax = modf(center_x)
		ayr, ay = modf(center_y)

		if ax + size_x > raster_size[0]:
		size_x = floor(raster_size[0] - center_x)
		if ax - size_x < 0:
		size_x = int(ax)
		if ay + size_y > raster_size[1]:
		size_y =floor(raster_size[1] - center_y)
		if ay - size_y < 0:
		size_y = int(ay)

		# Read from the upper left origin
		pixels = [ax-size_x, ay-size_y, size_x2, size_y2]
		pixels = list(map(int, pixels)) #
		if isinstance(img, np.ndarray):
		subarray = img[pixels[1]:pixels[1] + pixels[3] + 1, pixels[0]:pixels[0] + pixels[2] + 1]
		else:
		try:
		subarray = img.read_array(pixels=pixels, dtype=dtype)
		except:
		return None, 0, 0
		return subarray, axr, ayr

		def subpixel_phase(sx, sy, dx, dy,
		s_img, d_img,
		image_size=(51, 51),
		@@ -237,8 +180,8 @@ def subpixel_phase(sx, sy, dx, dy,
		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()
		s_image = s_roi.array
		d_template = d_roi.array

		if s_image.shape != d_template.shape:

		@@ -263,8 +206,8 @@ def subpixel_phase(sx, sy, dx, dy,
		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()
		s_image = s_roi.array
		d_template = d_roi.array

		if (s_image is None) or (d_template is None):
		return None, None, None
		@@ -374,8 +317,8 @@ def subpixel_transformed_template(sx, sy, dx, dy,
		except:
		d_template_dtype = None

		s_image = bytescale(s_roi.clip(dtype=s_image_dtype))
		d_template = bytescale(d_roi.clip(dtype=d_template_dtype))
		s_image = bytescale(s_roi.array)
		d_template = bytescale(d_roi.array)

		if verbose:
		fig, axs = plt.subplots(1, 5, figsize=(20,10))
		@@ -510,8 +453,8 @@ def subpixel_template_classic(sx, sy, dx, dy,
		print(f'd mm: {d_roi.clip().min()} {d_roi.clip().max()}')"""
		#print(f'{len(isis.get_isis_special_pixels(d_roi.clip()))} chip sps : ', isis.get_isis_special_pixels(d_roi.clip()))

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

		"""print('s shape', s_image.shape)
		print('s mean: ', s_image.mean())
		@@ -526,6 +469,7 @@ def subpixel_template_classic(sx, sy, dx, dy,
		return None, None, None, None

		shift_x, shift_y, metrics, corrmap = func(img_as_float32(d_template), img_as_float32(s_image), **kwargs)
		print(shift_x, shift_y)
		if shift_x is None:
		return None, None, None, None
		# Apply the shift and return
		@@ -622,8 +566,8 @@ def subpixel_template(sx, sy, dx, dy,
		except:
		d_template_dtype = None

		s_image = bytescale(s_roi.clip(dtype=s_image_dtype))
		d_template = bytescale(d_roi.clip(dtype=d_template_dtype))
		s_image = s_roi.array
		d_template = d_roi.array

		if (s_image is None) or (d_template is None):
		return None, None, None, None
		@@ -690,8 +634,8 @@ def subpixel_ciratefi(sx, sy, dx, dy, s_img, d_img, search_size=251, template_si
		size_x=template_size, size_y=template_size)
		s_roi = roi.Roi(s_img, sx, sy,
		size_x=search_size, size_y=search_size)
		template = t_roi.clip()
		search = s_roi.clip()
		template = t_roi.array
		search = s_roi.array

		if template is None or search is None:
		return None, None, None
		@@ -795,8 +739,8 @@ def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=(51, 51), reduction=11, c

		return tf.estimate_transform('affine', destination_coordinates, source_coordinates)
		'''
		def geom_match_simple(base_cube,
		input_cube,
		def geom_match_simple(reference_image,
		moving_image,
		bcenter_x,
		bcenter_y,
		size_x=60,
		@@ -859,63 +803,23 @@ def geom_match_simple(base_cube,
		autocnet.matcher.subpixel.subpixel_phase: for list of kwargs that can be passed to the matcher
		"""
		t1 = time.time()
		if not isinstance(input_cube, GeoDataset):
		raise Exception("input cube must be a geodataset obj")
		if not isinstance(base_cube, GeoDataset):
		raise Exception("match cube must be a geodataset obj")
		if not isinstance(reference_image, GeoDataset):
		raise TypeError("reference_image must be a GeoDataset obj")
		if not isinstance(moving_image, GeoDataset):
		raise TypeError("moving_image must be a GeoDataset obj")

		# Parse the match_func into a function if it comes in as a string
		if not callable(match_func):
		match_func = check_match_func(match_func)

		base_startx = int(bcenter_x - size_x)
		base_starty = int(bcenter_y - size_y)
		base_stopx = int(bcenter_x + size_x)
		base_stopy = int(bcenter_y + size_y)

		image_size = input_cube.raster_size
		match_size = base_cube.raster_size

		# for now, require the entire window resides inside both cubes.
		if base_stopx > match_size[0]:
		raise Exception(f"Window: {base_stopx} > {match_size[0]}, center: {bcenter_x},{bcenter_y}")
		if base_startx < 0:
		raise Exception(f"Window: {base_startx} < 0, center: {bcenter_x},{bcenter_y}")
		if base_stopy > match_size[1]:
		raise Exception(f"Window: {base_stopy} > {match_size[1]}, center: {bcenter_x},{bcenter_y} ")
		if base_starty < 0:
		raise Exception(f"Window: {base_starty} < 0, center: {bcenter_x},{bcenter_y}")

		# specifically not putting this in a try/except, this should never fail
		center_x, center_y = bcenter_x, bcenter_y

		base_corners = [(base_startx,base_starty),
		(base_startx,base_stopy),
		(base_stopx,base_stopy),
		(base_stopx,base_starty),
		(bcenter_x, bcenter_y)]

		dst_corners = []
		for x,y in base_corners:
		try:
		lon, lat = spatial.isis.image_to_ground(base_cube.file_name, x, y)
		dst_corners.append(
		spatial.isis.ground_to_image(input_cube.file_name, lon, lat)
		)
		except: pass

		if len(dst_corners) < 3:
		raise ValueError('Unable to find enough points to compute an affine transformation.')

		base_gcps = np.array([*base_corners])

		dst_gcps = np.array([*dst_corners])

		affine = tf.estimate_transform('affine', np.array([base_gcps]), np.array([dst_gcps]))
		# Estimate the transformation between the two images
		affine = estimate_affine_transformation(reference_image, moving_image, bcenter_x, bcenter_y)
		t2 = time.time()
		print(f'Estimation of the transformation took {t2-t1} seconds.')
		# read_array not getting correct type by default



		# Read the arrays with the correct dtype - needs to be in ROI object for dtype checking.
		base_type = isis.isis2np_types[pvl.load(base_cube.file_name)["IsisCube"]["Core"]["Pixels"]["Type"]]
		base_arr = base_cube.read_array(dtype=base_type)

		@@ -931,9 +835,9 @@ def geom_match_simple(base_cube,
		if verbose:
		fig, axs = plt.subplots(1, 2)
		axs[0].set_title("Base")
		axs[0].imshow(roi.Roi(bytescale(base_arr, cmin=0), bcenter_x, bcenter_y, 25, 25).clip(), cmap="Greys_r")
		axs[0].imshow(roi.Roi(bytescale(base_arr, cmin=0), bcenter_x, bcenter_y, 25, 25).array, cmap="Greys_r")
		axs[1].set_title("Projected Image")
		axs[1].imshow(roi.Roi(bytescale(dst_arr, cmin=0), bcenter_x, bcenter_y, 25, 25).clip(), cmap="Greys_r")
		axs[1].imshow(roi.Roi(bytescale(dst_arr, cmin=0), bcenter_x, bcenter_y, 25, 25).array, cmap="Greys_r")
		plt.show()
		# Run through one step of template matching then one step of phase matching
		# These parameters seem to work best, should pass as kwargs later
		@@ -966,19 +870,19 @@ def geom_match_simple(base_cube,
		fig, axs = plt.subplots(2, 3)
		fig.set_size_inches((30,30))

		oarr = roi.Roi(input_cube.read_array(), sample, line, 150, 150).clip()
		oarr = roi.Roi(input_cube.read_array(), sample, line, 150, 150).array
		axs[0][2].imshow(bytescale(oarr, cmin=0), cmap="Greys_r")
		axs[0][2].axhline(y=oarr.shape[1]/2, color="red", linestyle="-", alpha=1)
		axs[0][2].axvline(x=oarr.shape[1]/2, color="red", linestyle="-", alpha=1)
		axs[0][2].set_title("Original Registered Image")

		barr = roi.Roi(base_arr, bcenter_x, bcenter_y, size_x, size_y).clip()
		barr = roi.Roi(base_arr, bcenter_x, bcenter_y, size_x, size_y).array
		axs[0][0].imshow(bytescale(barr, cmin=0), cmap="Greys_r")
		axs[0][0].axhline(y=barr.shape[1]/2, color="red", linestyle="-", alpha=1)
		axs[0][0].axvline(x=barr.shape[1]/2, color="red", linestyle="-", alpha=1)
		axs[0][0].set_title("Base")

		darr = roi.Roi(dst_arr, x, y, size_x, size_y).clip()
		darr = roi.Roi(dst_arr, x, y, size_x, size_y).array
		axs[0][1].imshow(bytescale(darr, cmin=0), cmap="Greys_r")
		axs[0][1].axhline(y=darr.shape[1]/2, color="red", linestyle="-", alpha=1)
		axs[0][1].axvline(x=darr.shape[1]/2, color="red", linestyle="-", alpha=1)
		@@ -1149,7 +1053,7 @@ def geom_match_classic(base_cube,
		if verbose:
		fig, axs = plt.subplots(1, 3)
		fig.set_size_inches((30,30))
		darr = roi.Roi(input_cube.read_array(dtype=dst_type), sample, line, 100, 100).clip()
		darr = roi.Roi(input_cube.read_array(dtype=dst_type), sample, line, 100, 100).array
		axs[1].imshow(darr, cmap="Greys_r")
		axs[1].scatter(x=[darr.shape[1]/2], y=[darr.shape[0]/2], s=10, c="red")
		axs[1].set_title("Original Registered Image")
		@@ -2145,7 +2049,6 @@ def subpixel_register_point_smart(pointid,
		print('geom_match image:', destination_node['image_path'])
		print('geom_func', geom_func)

		# Apply the transformation
		try:
		affine = estimate_affine_transformation(source_node.geodata,
		destination_node.geodata,
		@@ -2161,6 +2064,7 @@ def subpixel_register_point_smart(pointid,
		updated_measures.append([None, None, m])
		continue

		# Here is where I need to get the two ROIs extracted. Then I need to get the destination ROI affine transformed to the source ROI
		base_arr, dst_arr = affine_warp_image(source_node.geodata,
		destination_node.geodata,
		affine)
		@@ -2175,8 +2079,8 @@ def subpixel_register_point_smart(pointid,
		if match_kwarg['template_size'][1] < size_y:
		size_y = match_kwarg['template_size'][1]

		base_roi = roi.Roi(base_arr, source.apriorisample, source.aprioriline, size_x=size_x, size_y=size_y).clip()
		dst_roi = roi.Roi(dst_arr, source.apriorisample, source.aprioriline, size_x=size_x, size_y=size_y).clip()
		base_roi = roi.Roi(base_arr, source.apriorisample, source.aprioriline, size_x=size_x, size_y=size_y).array
		dst_roi = roi.Roi(dst_arr, source.apriorisample, source.aprioriline, size_x=size_x, size_y=size_y).array

		if np.isnan(base_roi).any() or np.isnan(dst_roi).any():
		print('Unable to process due to NaN values in the input data.')
		@@ -2221,8 +2125,8 @@ def subpixel_register_point_smart(pointid,
		print('Unable to match with this parameter set.')
		continue

		base_roi = roi.Roi(base_arr, source.apriorisample, source.aprioriline, size_x=size_x, size_y=size_y).clip()
		dst_roi = roi.Roi(dst_arr, x, y, size_x=size_x, size_y=size_y).clip()
		base_roi = roi.Roi(base_arr, source.apriorisample, source.aprioriline, size_x=size_x, size_y=size_y).array
		dst_roi = roi.Roi(dst_arr, x, y, size_x=size_x, size_y=size_y).array

		#TODO: When refactored, all this type conversion should happen in the ROI object.
		base_roi = img_as_float32(base_roi)

autocnet/transformation/affine.py

+16 −19

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

		log = logging.getLogger(__name__)

		def estimate_affine_transformation(base_cube,
		input_cube,
		def estimate_affine_transformation(reference_image,
		moving_image,
		bcenter_x,
		bcenter_y,
		size_x=60,
		size_y=60):
		"""
		Using the a priori sensor model, project corner and center points from the base_cube into
		the input_cube and use these points to estimate an affine transformation.
		Using the a priori sensor model, project corner and center points from the reference_image into
		the moving_image and use these points to estimate an affine transformation.

		Parameters
		----------
		base_cube: plio.io.io_gdal.GeoDataset
		reference_image: plio.io.io_gdal.GeoDataset
		source image
		input_cube: plio.io.io_gdal.GeoDataset
		moving_image: plio.io.io_gdal.GeoDataset
		destination image; gets matched to the source image
		bcenter_x: int
		sample location of source measure in base_cube
		sample location of source measure in reference_image
		bcenter_y: int
		line location of source measure in base_cube
		line location of source measure in reference_image
		size_x: int
		half-height of the subimage used in the affine transformation
		size_y: int
		@@ -40,17 +40,17 @@ def estimate_affine_transformation(base_cube,

		"""
		t1 = time.time()
		if not isinstance(input_cube, GeoDataset):
		raise Exception(f"Input cube must be a geodataset obj, but is type {type(input_cube)}.")
		if not isinstance(base_cube, GeoDataset):
		raise Exception(f"Match cube must be a geodataset obj, but is type {type(base_cube)}.")
		if not isinstance(moving_image, GeoDataset):
		raise Exception(f"Input cube must be a geodataset obj, but is type {type(moving_image)}.")
		if not isinstance(reference_image, GeoDataset):
		raise Exception(f"Match cube must be a geodataset obj, but is type {type(reference_image)}.")

		base_startx = int(bcenter_x - size_x)
		base_starty = int(bcenter_y - size_y)
		base_stopx = int(bcenter_x + size_x)
		base_stopy = int(bcenter_y + size_y)

		match_size = base_cube.raster_size
		match_size = reference_image.raster_size

		# for now, require the entire window resides inside both cubes.
		if base_stopx > match_size[0]:
		@@ -66,20 +66,17 @@ def estimate_affine_transformation(base_cube,
		y_coords = [base_starty, base_stopy, base_stopy, base_starty, bcenter_y]

		# Dispatch to the sensor to get the a priori pixel location in the input image
		lons, lats = isis.image_to_ground(base_cube.file_name, x_coords, y_coords, allowoutside=True)
		xs, ys = isis.ground_to_image(input_cube.file_name, lons, lats, allowoutside=True)
		lons, lats = isis.image_to_ground(reference_image.file_name, x_coords, y_coords, allowoutside=True)
		xs, ys = isis.ground_to_image(moving_image.file_name, lons, lats, allowoutside=True)

		# Check for any coords that do not project between images
		base_gcps = []
		dst_gcps = []
		for i, (base_x, base_y) in enumerate(zip(xs, ys)):
		for i, (base_x, base_y) in enumerate(zip(x_coords, y_coords)):
		if xs[i] is not None and ys[i] is not None:
		dst_gcps.append((xs[i], ys[i]))
		base_gcps.append((base_x, base_y))

		base_gcps = np.asarray(base_gcps)
		dst_gcps = np.asarray(dst_gcps)

		log.debug(f'base_gcps: {base_gcps}')
		log.debug(f'dst_gcps: {dst_gcps}')

autocnet/transformation/roi.py

+6 −26

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

		from autocnet.transformation import dtypes
		import pvl

		class Roi():
		"""
		@@ -38,11 +40,10 @@ class Roi():
		bottom_y : int
		The bottom image coordinate in imge space
		"""
		def __init__(self, data, x, y, size_x=200, size_y=200, dtype=None, ndv=None, ndv_threshold=0.5):
		def __init__(self, data, x, y, size_x=200, size_y=200, ndv=None, ndv_threshold=0.5):
		self.data = data
		self.x = x
		self.y = y
		self.dtype = dtype
		self.size_x = size_x
		self.size_y = size_y
		self.ndv = ndv
		@@ -140,7 +141,7 @@ class Roi():
		@property
		def center(self):
		ie = self.image_extent
		return (ie[1] - ie[0])/2, (ie[3]-ie[2])/2
		return (ie[1] - ie[0])/2.+0.5, (ie[3]-ie[2])/2.+0.5

		@property
		def is_valid(self):
		@@ -165,28 +166,7 @@ class Roi():
		data = self.data[pixels[2]:pixels[3]+1,pixels[0]:pixels[1]+1]
		else:
		# Have to reformat to [xstart, ystart, xnumberpixels, ynumberpixels]
		# TODO: I think this will result in an incorrect obj.center when the passed data is a GeoDataset
		pixels = [pixels[0], pixels[2], pixels[1]-pixels[0]+1, pixels[3]-pixels[2]+1]
		data = self.data.read_array(pixels=pixels, dtype=self.dtype)
		data = self.data.read_array(pixels=pixels)
		return data

		def clip(self, dtype=None):
		"""
		Compatibility function that makes a call to the array property.

		Warning: The dtype passed in via this function resets the dtype attribute of this
		instance.

		Parameters
		----------
		dtype : str
		The datatype to be used when reading the ROI information if the read
		occurs through the data object using the read_array method. When using
		this object when the data are a numpy array the dtype has not effect.

		Returns
		-------
		: ndarray
		The array attribute of this object.
		"""
		self.dtype = dtype
		return self.array