Unverified Commit 18778ceb authored by acpaquette's avatar acpaquette Committed by GitHub
Browse files

Phase Roi Update (#629)



* Initial change to get iterative_phase into affine transformation returns

* Final changes to get the phase matcher updated to use Rois

* Adds clipping and addresses the test in subpixel_phase and iterative_phase

* Ignored signs on new translation computation in iterative_phase

Co-authored-by: default avatarjlaura <jlaura@usgs.gov>
parent 7e830432
Loading
Loading
Loading
Loading
+157 −169
Original line number Diff line number Diff line
@@ -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):
    """
+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]

+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