Commit afbcd51d authored by Lauren Adoram-Kershner's avatar Lauren Adoram-Kershner Committed by GitHub
Browse files

Null fit (#508) 

* initial commit

* adding type check bug fix

* addressing comments round 1

* removing iterative kwargs from function calls

* addressing comments round 2

* adding note, fixing typo

* bug fix & adding classic matching methods
parent 46ff71ea

autocnet/matcher/subpixel.py

+333 −1
Original line number Diff line number Diff line
@@ -399,6 +399,70 @@ def subpixel_transformed_template(sx, sy, dx, dy, 


    return dx, dy, metrics, corrmap



def subpixel_template_classic(sx, sy, dx, dy,

                              s_img, d_img,

                              image_size=(251, 251),

                              template_size=(51,51),

                              func=pattern_match,

                              **kwargs):

    """

    Uses a pattern-matcher on subsets of two images determined from the passed-in keypoints and optional sizes to

    compute an x and y offset from the search keypoint to the template keypoint and an associated strength.

    Parameters

    ----------

    sx : Numeric

         Source X coordinate

    sy : Numeric

         Source y coordinate

    dx : Numeric

         The desintation x coordinate

    dy : Numeric

         The destination y coordinate

    s_img : GeoDataset

            The source image GeoDataset

    d_img : GeoDataset

            The destination image GeoDataset

    image_size : tuple

                 (xsize, ysize) of the image that is searched within (this should be larger

                 than the template size)

    template_size : tuple

                    (xsize, ysize) of the template to iterate over the image in order

                    to identify the area(s) of highest correlation.



    Returns

    -------

    x_shift : float

              Shift in the x-dimension

    y_shift : float

              Shift in the y-dimension

    strength : float

               Strength of the correspondence in the range [-1, 1]

    See Also

    --------

    autocnet.matcher.naive_template.pattern_match : for the kwargs that can be passed to the matcher

    autocnet.matcher.naive_template.pattern_match_autoreg : for the jwargs that can be passed to the autoreg style matcher

    """



    image_size = check_image_size(image_size)

    template_size = check_image_size(template_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=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, None



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



    dx = d_roi.x - shift_x

    dy = d_roi.y - shift_y



    return dx, dy, metrics, corrmap





def subpixel_template(sx, sy, dx, dy,

                      s_img, d_img,

                      image_size=(251, 251),
@@ -656,8 +720,267 @@ def estimate_affine_transformation(destination_coordinates, source_coordinates): 
    source_coordinates = np.asarray(source_coordinates)



    return tf.estimate_transform('affine', destination_coordinates, source_coordinates)

def geom_match_classic(base_cube,

                       input_cube,

                       bcenter_x,

                       bcenter_y,

                       size_x=60,

                       size_y=60,

                       template_kwargs={"image_size":(59,59), "template_size":(31,31)},

                       phase_kwargs=None,

                       verbose=True):

    """

    Propagates a source measure into destination images and then perfroms subpixel registration.

    Measure creation is done by projecting the (lon, lat) associated with the source measure into the

    destination image. The created measure is then matched to the source measure using a quick projection

    of the destination image into source image space (using an affine transformation) and a naive

    template match with optional phase template match.

    Parameters

    ----------

    base_cube:  plio.io.io_gdal.GeoDataset

                source image

    input_cube: plio.io.io_gdal.GeoDataset

                destination image; gets matched to the source image

    bcenter_x:  int

                sample location of source measure in base_cube

    bcenter_y:  int

                line location of source measure in base_cube

    size_x:     int

                half-height of the subimage used in the affine transformation

    size_y:     int

                half-width of the subimage used in affine transformation

    template_kwargs: dict

                     contains keywords necessary for autocnet.matcher.subpixel.subpixel_template

    phase_kwargs:    dict

                     contains kwargs for autocnet.matcher.subpixel.subpixel_phase

    verbose:    boolean

                indicates level of print out desired. If True, two subplots are output; the first subplot contains

                the source subimage and projected destination subimage, the second subplot contains the registered

                measure's location in the base subimage and the unprojected destination subimage with the corresponding

                template metric correlation map.

    Returns

    -------

    sample: int

            sample of new measure in destination image space

    line:   int

            line of new measures in destination image space

    dist:   np.float or tuple of np.float

            distance matching algorithm moved measure

            if template matcher only (default): returns dist_template

            if template and phase matcher:      returns (dist_template, dist_phase)

    metric: np.float or tuple of np.float

            matching metric output by the matcher

            if template matcher only (default): returns maxcorr

            if template and phase matcher:      returns (maxcorr, perror, pdiff)

    temp_corrmap: np.ndarray

            correlation map of the naive template matcher

    See Also

    --------

    autocnet.matcher.subpixel.subpixel_template: for list of kwargs that can be passed to the matcher

    autocnet.matcher.subpixel.subpixel_phase: for list of kwargs that can be passed to the matcher

    """



    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")



    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

    mlat, mlon = spatial.isis.image_to_ground(base_cube.file_name, bcenter_x, bcenter_y)

    center_x, center_y = spatial.isis.ground_to_image(input_cube.file_name, mlon, mlat)[::-1]



    base_corners = [(base_startx,base_starty),

                    (base_startx,base_stopy),

                    (base_stopx,base_stopy),

                    (base_stopx,base_starty)]



    dst_corners = []

    for x,y in base_corners:

        try:

            lat, lon = 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)[::-1])

        except ProcessError as e:

            if 'Requested position does not project in camera model' in e.stderr:

                print(f'Skip geom_match; Region of interest corner located at ({lon}, {lat}) does not project to image {input_cube.base_name}')

                return None, None, None, None, None



    base_gcps = np.array([*base_corners])

    base_gcps[:,0] -= base_startx

    base_gcps[:,1] -= base_starty



    dst_gcps = np.array([*dst_corners])

    start_x = dst_gcps[:,0].min()

    start_y = dst_gcps[:,1].min()

    stop_x = dst_gcps[:,0].max()

    stop_y = dst_gcps[:,1].max()

    dst_gcps[:,0] -= start_x

    dst_gcps[:,1] -= start_y



    affine = tf.estimate_transform('affine', np.array([*base_gcps]), np.array([*dst_gcps]))



    # read_array not getting correct type by default

    isis2np_types = {

                    "UnsignedByte" : "uint8",

                    "SignedWord" : "int16",

                    "Real" : "float64"

    }



    base_pixels = list(map(int, [base_corners[0][0], base_corners[0][1], size_x*2, size_y*2]))

    base_type = isis2np_types[pvl.load(base_cube.file_name)["IsisCube"]["Core"]["Pixels"]["Type"]]

    base_arr = base_cube.read_array(pixels=base_pixels, dtype=base_type)



    dst_pixels = list(map(int, [start_x, start_y, stop_x-start_x, stop_y-start_y]))

    dst_type = isis2np_types[pvl.load(input_cube.file_name)["IsisCube"]["Core"]["Pixels"]["Type"]]

    dst_arr = input_cube.read_array(pixels=dst_pixels, dtype=dst_type)



    dst_arr = tf.warp(dst_arr, affine)

    dst_arr = dst_arr[:size_y*2, :size_x*2]



    if verbose:

        fig, axs = plt.subplots(1, 2)

        axs[0].set_title("Base")

        axs[0].imshow(bytescale(base_arr), cmap="Greys_r")

        axs[1].set_title("Projected Image")

        axs[1].imshow(bytescale(dst_arr), 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

    restemplate = subpixel_template_classic(size_x, size_y, size_x, size_y, bytescale(base_arr), bytescale(dst_arr), **template_kwargs)



    x,y,maxcorr,temp_corrmap = restemplate

    if x is None or y is None:

        return None, None, None, None, None

    metric = maxcorr

    sample, line = affine([x, y])[0]

    sample += start_x

    line += start_y

    dist = np.linalg.norm([center_x-sample, center_y-line])



    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()

        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")



        axs[0].imshow(base_arr, cmap="Greys_r")

        axs[0].scatter(x=[base_arr.shape[1]/2], y=[base_arr.shape[0]/2], s=10, c="red")

        axs[0].set_title("Base")



        pcm = axs[2].imshow(temp_corrmap**2, interpolation=None, cmap="coolwarm")

        plt.show()



    return sample, line, dist, metric, temp_corrmap





def geom_match(base_cube,

               input_cube,

               bcenter_x,

               bcenter_y,

               size_x=60,

               size_y=60,

               template_kwargs={"image_size":(59,59), "template_size":(31,31)},

               phase_kwargs=None,

               verbose=True):

    """

    Propagates a source measure into destination images and then perfroms subpixel registration.

    Measure creation is done by projecting the (lon, lat) associated with the source measure into the

    destination image. The created measure is then matched to the source measure using a quick projection

    of the destination image into source image space (using an affine transformation) and a naive

    template match with optional phase template match.

    Parameters

    ----------

    base_cube:  plio.io.io_gdal.GeoDataset

                source image

    input_cube: plio.io.io_gdal.GeoDataset

                destination image; gets matched to the source image

    bcenter_x:  int

                sample location of source measure in base_cube

    bcenter_y:  int

                line location of source measure in base_cube

    size_x:     int

                half-height of the subimage used in the affine transformation

    size_y:     int

                half-width of the subimage used in affine transformation

    template_kwargs: dict

                    contains keywords necessary for autocnet.matcher.subpixel.subpixel_template

    phase_kwargs:   dict

                    contains kwargs for autocnet.matcher.subpixel.subpixel_phase

    verbose: boolean

             indicates level of print out desired. If True, two subplots are output; the first subplot contains

             the source subimage and projected destination subimage, the second subplot contains the registered

             measure's location in the base subimage and the unprojected destination subimage with the corresponding

             template metric correlation map.



    Returns

    -------

    sample: int

            sample of new measure in destination image space

    line:   int

            line of new measures in destination image space

    dist:   np.float or tuple of np.float

            distance matching algorithm moved measure

            if template matcher only (default): returns dist_template

            if template and phase matcher:      returns (dist_template, dist_phase)

    metric: np.float or tuple of np.float

            matching metric output by the matcher

            if template matcher only (default): returns maxcorr

            if template and phase matcher:      returns (maxcorr, perror, pdiff)

    temp_corrmap: np.ndarray

                  correlation map of the naive template matcher

    See Also

    --------

    autocnet.matcher.subpixel.subpixel_template: for list of kwargs that can be passed to the matcher

    autocnet.matcher.subpixel.subpixel_phase: for list of kwargs that can be passed to the matcher

    """



    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")



    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

    mlat, mlon = spatial.isis.image_to_ground(base_cube.file_name, bcenter_x, bcenter_y)

    center_x, center_y = spatial.isis.ground_to_image(input_cube.file_name, mlon, mlat)[::-1]



def geom_match(destination_cube,

               source_cube,

               bcenter_x,
@@ -923,6 +1246,7 @@ def subpixel_register_point(pointid, 
                            cost_func=lambda x,y: 1/x**2 * y,

                            threshold=0.005,

                            ncg=None,

                            version='new',

                            **kwargs):



    """
@@ -954,6 +1278,14 @@ def subpixel_register_point(pointid, 
    if not ncg.Session:

        raise BrokenPipeError('This func requires a database session from a NetworkCandidateGraph.')



    version = version.lower()

    geom_funcs = {"classic": geom_match_classic,

                "new": geom_match

                }

    if version not in geom_funcs.keys():

        raise Exception(f"{version} not a valid geom_match function version.")

    geom_func = geom_funcs[version]



    if isinstance(pointid, Points):

        pointid = pointid.id
@@ -987,7 +1319,7 @@ def subpixel_register_point(pointid, 


            print('geom_match image:', res.path)

            try:

                new_x, new_y, dist, metric,  _ = geom_match(source_node.geodata, destination_node.geodata,

                new_x, new_y, dist, metric,  _ = geom_func(source_node.geodata, destination_node.geodata,

                                                        source.apriorisample, source.aprioriline,

                                                        template_kwargs=subpixel_template_kwargs)

            except Exception as e:

autocnet/transformation/roi.py

+23 −3
Original line number Diff line number Diff line
@@ -78,6 +78,26 @@ class Roi(): 
    def ndv(self, ndv):

        self._ndv = ndv



    @property

    def size_x(self):

        return self._size_x



    @size_x.setter

    def size_x(self, size_x):

        if not isinstance(size_x, int):

            raise TypeError(f'size_x must be type integer, not {type(size_x)}')

        self._size_x = size_x



    @property

    def size_y(self):

        return self._size_y



    @size_y.setter

    def size_y(self, size_y):

        if not isinstance(size_y, int):

            raise TypeError(f'size_y must be type integer, not {type(size_y)}')

        self._size_y = size_y



    @property

    def image_extent(self):

        """
@@ -95,7 +115,7 @@ class Roi(): 
        left_x = self._x - self.size_x

        right_x = self._x + self.size_x

        top_y = self._y - self.size_y

        bottom_y = self.y + self.size_y

        bottom_y = self._y + self.size_y



        if self._x - self.size_x < 0:

            left_x = 0
@@ -125,14 +145,14 @@ class Roi(): 
    @property

    def array(self):

        """

        The clopped array associated with this ROI.

        The clipped array associated with this ROI.

        """

        pixels = self.image_extent

        if isinstance(self.data, np.ndarray):

             return self.data[pixels[2]:pixels[3]+1,pixels[0]:pixels[1]+1]

        else:

            # Have to reformat to [xstart, ystart, xnumberpixels, ynumberpixels]

            pixels = [pixels[0], pixels[2], pixels[1]-pixels[0], pixels[3]-pixels[2]]

            pixels = [pixels[0], pixels[2], pixels[1]-pixels[0]+1, pixels[3]-pixels[2]+1]

            return self.data.read_array(pixels=pixels, dtype=self.dtype)



    def clip(self, dtype=None):