Unverified Commit fb6e1831 authored by Kelvin Rodriguez's avatar Kelvin Rodriguez Committed by GitHub
Browse files

geom_match_simple (#511) 

* simp

* forced simple method

* updated simple algo

* forced simple method

* removed prints

* removed more prints
parent 6e514c96

autocnet/matcher/cross_instrument_matcher.py

+6 −6
Original line number Diff line number Diff line
@@ -50,7 +50,7 @@ from autocnet.spatial import isis 
from autocnet.transformation.spatial import reproject, oc2og

from autocnet.matcher.cpu_extractor import extract_most_interesting

from autocnet.transformation import roi

from autocnet.matcher.subpixel import geom_match

from autocnet.matcher.subpixel import geom_match_simple

from autocnet.utils.utils import bytescale



import warnings
@@ -277,7 +277,7 @@ def propagate_point(Session, 
                print(f'prop point: base_image: {base_image}')

                print(f'prop point: dest_image: {dest_image}')

                print(f'prop point: (sx, sy): ({sx}, {sy})')

                x,y, dist, metrics, corrmap = geom_match(base_image, dest_image, sx, sy, \

                x,y, dist, metrics, corrmap = geom_match_simple(base_image, dest_image, sx, sy, 16, 16, \

                        template_kwargs=template_kwargs, \

                        verbose=verbose)

            except Exception as e:

autocnet/matcher/subpixel.py

+199 −5
Original line number Diff line number Diff line
@@ -3,6 +3,7 @@ from math import modf, floor 
import numpy as np

import warnings



import sys

from skimage.feature import register_translation

from skimage import transform as tf
@@ -11,6 +12,11 @@ from matplotlib import pyplot as plt 
from plio.io.io_gdal import GeoDataset

from pysis.exceptions import ProcessError



import pvl



import PIL

from PIL import Image



from autocnet.matcher.naive_template import pattern_match, pattern_match_autoreg

from autocnet.matcher import ciratefi

from autocnet.io.db.model import Measures, Points, Images, JsonEncoder
@@ -19,10 +25,7 @@ from autocnet.transformation import roi 
from autocnet import spatial

from autocnet.utils.utils import bytescale



import geopandas as gpd

import pandas as pd



import pvl

PIL.Image.MAX_IMAGE_PIXELS = sys.float_info.max



isis2np_types = {

        "UnsignedByte" : "uint8",
@@ -720,6 +723,197 @@ 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_simple(base_cube,

                       input_cube,

                       bcenter_x,

                       bcenter_y,

                       size_x=60,

                       size_y=60,

                       template_kwargs={"image_size":(101,101), "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.



    This version projects the entirity of the input cube onto the base



    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

    """

    print("in geommatch")

    print("subpixel kwargs", template_kwargs)



    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

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



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



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



    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(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(dtype=dst_type)



    box = (0, 0, max(dst_arr.shape[1], base_arr.shape[1]), max(dst_arr.shape[0], base_arr.shape[0]))

    dst_arr = np.array(Image.fromarray(dst_arr).crop(box))



    dst_arr = tf.warp(dst_arr, affine)



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

        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(bcenter_x, bcenter_y, bcenter_x, bcenter_y, bytescale(base_arr, cmin=0), bytescale(dst_arr, cmin=0), **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]

    dist = np.linalg.norm([bcenter_x-x, bcenter_y-y])



    if verbose:

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

        fig.set_size_inches((30,30))



        oarr = roi.Roi(input_cube.read_array(), sample, line, 150, 150).clip()

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



        darr = roi.Roi(dst_arr, x, y, size_x, size_y).clip()

        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)

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



        barr = roi.Roi(base_arr, bcenter_x, bcenter_y, size_x, size_y).clip()

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



        axs[1][0].imshow(bytescale(darr.astype("f")/barr.astype("f")), cmap="Greys_r", alpha=.6)

        axs[1][0].axhline(y=barr.shape[1]/2, color="red", linestyle="-", alpha=.5)

        axs[1][0].axvline(x=barr.shape[1]/2, color="red", linestyle="-", alpha=.5)

        axs[1][0].set_title("overlap")



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

        plt.show()



    return sample, line, dist, metric, temp_corrmap





def geom_match_classic(base_cube,

                       input_cube,

                       bcenter_x,
@@ -1197,7 +1391,7 @@ def subpixel_register_measure(measureid, 
                      'status': ''}



        try:

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

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

                                                        source.sample, source.line,

                                                        template_kwargs=subpixel_template_kwargs)

        except Exception as e:

bin/acn_cd.py

deleted100755 → 0
+0 −176
Original line number Diff line number Diff line 
#!/usr/bin/env python



import sys

import os

import argparse

from argparse import RawTextHelpFormatter

import yaml

import tempfile



from plio.io.io_gdal import GeoDataset, array_to_raster



from autocnet.graph.network import CandidateGraph

from autocnet.cg import change_detection as cd

from autocnet.examples import get_path

from autocnet.graph.network import CandidateGraph

from autocnet.graph.edge import Edge

from autocnet.spatial.isis import point_info

from autocnet.utils import hirise

from autocnet.utils.utils import bytescale

from autocnet.examples import get_path

from shapely.geometry import Point, MultiPoint

import numpy as np



from affine import Affine



import pysis

from pysis.exceptions import ProcessError

from pysis import isis



import warnings

warnings.simplefilter("ignore")





_cd_functions_ = {

    "okb" : cd.okubogar_detector,

    "okbm" : cd.okbm_detector,

    "blob" : cd.blob_detector

}



def poly_pixel_to_latlon(poly, affine, coord_transform):

    poly_type = type(poly)



    if poly_type == MultiPoint:

        x = [p.x for p in poly]

        y = [p.y for p in poly]

    elif poly_type == Point:

        x,y = poly.xy

    else:

        x,y = poly.exterior.coords.xy



    lonlats = []

    for xval,yval in zip(x,y):

        lon, lat = Affine.from_gdal(*affine) * (xval, yval)

        lon, lat, _ = coord_transform.TransformPoint(lon, lat)

        lonlats.append([lon, lat])



    if poly_type == Point:

        return Point(lonlats[0][0], lonlats[0][1])



    return poly_type(lonlats)







if __name__ == "__main__":

    cd_function_help_string = ("Change detection algorithm to use.\n"

                               "Okubogar method (okb). Simple method which produces an overlay image of change hotspots (i.e. a 2d histogram image of detected change density). Largely based on a method created by Chris Okubo and Brendon Bogar. Histogram step was added for readability, image1, image2 -> image subtraction/ratio -> feature extraction -> feature histogram.\n\n"

                               "Okubogar modified method (okbm). Experimental feature based change detection algorithm that expands on okobubogar to allow for programmatic change detection."

    )



    parser = argparse.ArgumentParser(description="Registers two image and runs a change detection algorithm on the pair of images. WARNING: Runs bundle adjust with update=yes, make sure you are using copies.")

    parser.add_argument('before', action='store', help='Path to image 1, generally the "before image"')

    parser.add_argument('after', action='store', help='Path to image 2, generally the "after image"')

    parser.add_argument('out', action='store', help='Output image path, csv with geometries are also written as a side cart file as a csv.')

    parser.add_argument('--algorithm', '-a', action='store', choices=_cd_functions_.keys(), help=cd_function_help_string, default='okb')

    parser.add_argument('--config', '-c', action='store', default=get_path('cd_config.yml'), help='path to json or yaml file containing parameters for change detection algorithms')

    parser.add_argument('--map','-m',  action='store', help='path to ISIS map file, determines the projection of the two registered images', default=os.path.join(os.environ["ISISROOT"], "appdata", "templates", "maps", "equirectangular.map"))

    parser.add_argument('--register','-r', action="store_true", default=False, help='Whether or not to register the two images, reccomended to set to false if the two images have been registered before.')

    parser.add_argument('--write-registered-cubes','-w', default=False, action="store_true", help='Pass this flag id you want to write out the projected cubes to disk. Useful if you want to run multiple cd algorithms without having to rerun the registration step.')



    args = parser.parse_args()



    with open(args.config) as f:

        config = yaml.load(f, Loader=yaml.FullLoader)



    # temp path for temp files

    dirpath = tempfile.mkdtemp()



    if args.register:

        # Point to the adjacency Graph

        adjacency = {args.before: [args.after], args.after: [args.before]}

        cg = CandidateGraph.from_adjacency(adjacency)



        # Apply SIFT to extract features

        cg.extract_features(extractor_method='vlfeat', extractor_parameters={"edge_thresh":20})

        cg.match()



        # Apply outlier detection

        cg.apply_func_to_edges(Edge.symmetry_check)

        cg.apply_func_to_edges(Edge.ratio_check)

        cg.minimum_spanning_tree()



        # Compute a homography and apply RANSAC

        cg.apply_func_to_edges(Edge.compute_fundamental_matrix, clean_keys=['ratio', 'symmetry'])



        # Generate ISIS compatible control network

        cg.generate_control_network(clean_keys=["fundamental"])



        # write cnet out to temp file, run it through bundle adjust.

        cnet_path = os.path.join(dirpath, "cnet.net")

        filelist_path = os.path.join(dirpath, "cnet.lis")



        cg.to_isis(cnet_path)



        try:

            output = isis.jigsaw(fromlist=filelist_path, cnet=cnet_path, onet=cnet_path, update="yes", **config['jigsaw'])

            print(output.decode())

        except ProcessError as e:

            print(e.stdout.decode('utf-8'))

            print(e.stderr)

            exit()



        if args.write_registered_cubes:

            before_proj = os.path.splitext(args.before)[0] + ".proj.cub"

            after_proj = os.path.splitext(args.after)[0] + ".proj.cub"

        else: # use the temp directory

            before_proj = os.path.join(dirpath, "before.cub")

            after_proj = os.path.join(dirpath, "after.cub")



        try:

            isis.cam2map(from_=args.before, to=before_proj, map=args.map)

            isis.cam2map(from_=args.after, to=after_proj, patchsize=8, map=before_proj, matchmap=True, warpalgorithm="REVERSEPATCH")

        except ProcessError as e:

            print(e.stderr)

            exit()



        args.before = before_proj

        args.after = after_proj



    before_proj_geo = GeoDataset(args.before)

    after_proj_geo = GeoDataset(args.after)



    if args.algorithm == 'blob':

        # Requires sub solar azmith

        ssa_path = "/tmp/ssa.cub" # os.path.join(dirpath, 'ssa.cub')

        try:

            isis.phocube(from_=args.before, to=ssa_path, subsolargroundazimuth=True)

        except ProcessError as e:

            print(e.stderr)

        print("created: ", ssa_path)

        ssa = GeoDataset(ssa_path).read_array(6)

        ret = _cd_functions_[args.algorithm.strip()](before_proj_geo, after_proj_geo, ssa, **config.get(args.algorithm, {}))

    else:

        ret = _cd_functions_[args.algorithm.strip()](before_proj_geo, after_proj_geo, **config.get(args.algorithm, {}))



    # for now, write out raster files assuming okb

    # make it match one of the projected images

    match_srs = before_proj_geo.dataset.GetProjection()

    match_gt = before_proj_geo.geotransform

    match_coord_trans = before_proj_geo.coordinate_transformation



    if os.path.splitext(args.out)[1] == '':

        args.out = args.out + ".tif"



    print(f"Writing {args.out}")

    array_to_raster(ret[1], args.out, projection=match_srs, geotransform=match_gt, outformat="GTiff")



    print(f"Writing {os.path.splitext(args.out)[0]+'.csv'}")



    geodf = ret[0]

    geodf["geometry"] = geodf['geometry']

    geodf['latlon_geometry'] = geodf['geometry'].apply(lambda x: poly_pixel_to_latlon(x, match_gt, match_coord_trans ))

    geodf['geometry'] = [g.wkt for g in geodf['geometry']]

    geodf['latlon_geometry'] = [g.wkt for g in geodf['latlon_geometry']]



    geodf.to_csv(os.path.splitext(args.out)[0]+".csv", index=False)