Added change detection bin script (#459)

from skimage.feature import blob_log, blob_doh

from math import sqrt, atan2, pi

import matplotlib

from matplotlib import pyplot as plt

from shapely import wkt

from shapely.geometry import Point, MultiPoint

import pandas as pd

import geopandas as gpd

import pysis

from autocnet.utils.utils import bytescale

from autocnet.matcher.cpu_extractor import extract_features

def image_diff(arr1, arr2):

     arr1 = arr1.astype("float32")

     arr2 = arr2.astype("float32")

     arr1[arr1 == 0] = np.nan

     arr2[arr2 == 0] = np.nan

     isis_null = pysis.specialpixels.SPECIAL_PIXELS['Real']['Null']

     arr1[arr1 == isis_null] = np.nan

     arr2[arr2 == isis_null] = np.nan

     diff = arr1-arr2

     diff[np.isnan(diff)] = 0

     return diff

def image_ratio(arr1, arr2):

     arr1 = arr1.astype("float32")

     arr2 = arr2.astype("float32")

     arr1[arr1 == 0] = np.nan

     arr2[arr2 == 0] = np.nan

     isis_null = pysis.specialpixels.SPECIAL_PIXELS['Real']['Null']

     arr1[arr1 == isis_null] = np.nan

     arr2[arr2 == isis_null] = np.nan

     ratio = arr1/arr2

     ratio[np.isnan(ratio)] = 0

     return ratio

def image_diff_sq(arr1, arr2):

     return image_diff(arr1, arr2)**2

func_map = {

    "diff" : image_diff,

    "diff_sq": image_diff_sq,

    "ratio" : image_ratio

}

def okubogar_detector(image1, image2, nbins=50, extractor_method="orb", image_func=image_diff,

                      extractor_kwargs={"nfeatures": 2000, "scaleFactor": 1.1, "nlevels": 1}):

     """

     Largely based on a method created by Chris Okubo and Brendon Bogar. Histogram step

     was added for readability.

     image1

           \

             image subtraction/ratio -> feature_extraction -> feature_histogram

           /

     image2

     image1, image2 -> image subtraction/ratio -> feature extraction -> feature histogram

     TODO: Paper/abstract might exist, cite

     extractor_kwargs : dict

                        A dictionary containing OpenCV SIFT parameters names and values.

     Returns

     -------

     : pd.DataFrame

       Dataframe containing polygon results as wkt

     : np.array

       Numpy array image, the image used to compute change

     See Also

     --------

     if isinstance(image2, GeoDataset):

         image2 = image2.read_array()

     if isinstance(image_func, str):

         try:

             image_func = func_map[image_func]

         except KeyError as e:

             raise Exception(f"{image_func} is not a valid method, available image functions: {func_map.keys()}")

     image1[image1 == image1.min()] = 0

     image2[image2 == image2.min()] = 0

     arr1 = bytescale(image1)

     keys, descriptors = extract_features(bdiff, extractor_method, extractor_parameters=extractor_kwargs)

     x,y = keys["x"], keys["y"]

     points = [Point(xval, yval) for xval,yval in zip(x,y)]

     heatmap, xedges, yedges = np.histogram2d(y, x, bins=nbins, range=[[0, bdiff.shape[0]], [0, bdiff.shape[1]]])

     heatmap = cv2.resize(heatmap, dsize=(bdiff.shape[1], bdiff.shape[0]), interpolation=cv2.INTER_NEAREST)

     return points, heatmap, bdiff

     #square image to improve signal to noise ratio

     heatmap = heatmap**2

     keys = gpd.GeoDataFrame(keys, geometry= gpd.points_from_xy(keys.x, keys.y))

     return keys, heatmap, bdiff

def okbm_detector(image1, image2, nbins=50, extractor_method="orb",  image_func=image_diff,

def okbm_detector(image1, image2, extractor_method="orb",  image_func=image_diff,

                 extractor_kwargs={"nfeatures": 2000, "scaleFactor": 1.1, "nlevels": 1},

                 cluster_params={"min_samples": 10, "max_eps": 10, "eps": .5, "xi":.5}):

                 cluster_kwargs={"min_samples": 10, "max_eps": 10, "eps": .5, "xi":.5}):

     """

     okobubogar modified detector, experimental feature based change detection algorithmthat expands on okobubogar to allow for

     programmatic change detection. Returns detected feature changes as weighted polygons.

                  >>> image1, image2 = random.random((50,50)), random.random((50,50))

                  >>> results = okbm_detector(image1, image2, image_func=lambda im1, im2: im1/im2)

     nbins : int

            number of bins to use in the 2d histogram

     extractor_method : {'orb', 'sift', 'fast', 'surf', 'vl_sift'}

               The detector method to be used.  Note that vl_sift requires that

               vlfeat and cyvlfeat dependencies be installed.

     extractor_kwargs : dict

                        A dictionary containing OpenCV SIFT parameters names and values.

     cluster_params : dict

     cluster_kwargs : dict

                      A dictionary containing sklearn.cluster.OPTICS parameters

     """

     if isinstance(image1, GeoDataset):

         image1 = image1.read_array()

     if isinstance(image2, GeoDataset):

         image2 = image2.read_array()

     image1[image1 == image1.min()] = 0

     image2[image2 == image2.min()] = 0

     arr1 = bytescale(image1)

     arr2 = bytescale(image2)

     Returns

     -------

     : pd.DataFrame

       Dataframe containing polygon results as wkt

     bdiff = image_func(arr1, arr2)

     : np.array

       Numpy array image, the image used to compute change

     keys, descriptors = extract_features(bdiff, extractor_method, extractor_parameters=extractor_kwargs)

     x,y = keys["x"], keys["y"]

     """

     keys, _, bdiff = okubogar_detector(image1, image2, 10, extractor_method, image_func, extractor_kwargs)

     x,y = keys['x'], keys['y']

     points = [Point(xval, yval) for xval,yval in zip(x,y)]

     optics = OPTICS(**cluster_params).fit(list(zip(x,y)))

     optics = OPTICS(**cluster_kwargs).fit(list(zip(x,y)))

     classes = gpd.GeoDataFrame(columns=["label", "point"], geometry="point")

     classes["label"] = optics.labels_

         polys.append(poly)

         weights.append(weight)

     return polys, weights, bdiff

     results = gpd.GeoDataFrame(geometry=polys)

     results['weight'] = weights

     return results, bdiff

def blob_detector(image1, image2, sub_solar_azimuth, image_func=image_diff_sq,

     Returns

     -------

     changes : np.ndarray

               A numpy array containing candidate change points in the form (y,x,radius)

     : pd.DataFrame

       A pandas dataframe containing a points of changed areas

     bdiff : np.ndarray

     : np.ndarray

       A numpy array containing the image upon which the change detection

       algorithm operates, i.e. the image resulting from image_func.

     if isinstance(image2, GeoDataset):

         image2 = image2.read_array()

     if isinstance(image_func, str):

         try:

             image_func = func_map[image_func]

         except KeyError as e:

             raise Exception(f"{image_func} is not a valid method, available image functions: {func_map.keys()}")

     bdiff = image_func(image1,image2)

     bdiff = bytescale(bdiff)

     blobs_log_inv = blob_log(inv, min_sigma=min_sigma, max_sigma=max_sigma,

                              num_sigma=num_sigma, threshold=threshold, overlap=overlap,

                              log_scale=log_scale, exclude_border=exclude_border)

     # Compute radii in the 3rd column.  Radii are appx equal to sqrt2 * sigma

     blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)

     blobs_log_inv[:, 2] = blobs_log_inv[:, 2] * sqrt(2)

     if not len(blobs_log) or not len(blobs_log_inv):

         raise Exception("No blobs detected")

     # Create a KDTree to facilitate nearest neighbor search

     tree = cKDTree(blobs_log)

     # Query the kdtree to find neighboring points

     _, idx_log = tree.query(blobs_log_inv, k=n_neighbors,

                                    distance_upper_bound=dist_upper_bound)

     # Nearest neighbors

     neighbors = [blobs_log[j] for j in [i[i!=len(blobs_log)]for i in idx_log] if j.size > 0]

     changes = []

     polys = []

     for idx, pt1 in enumerate(close_points):

         for pt2 in neighbors[idx]:

             try:

             except IndexError as e:

                 azimuth = sub_solar_azimuth

             if is_azimuth_colinear(pt1, pt2, azimuth, angle_tolerance, subtractive):

                 changes.append([pt1,pt2])

     changes = np.array(changes)

                 if subtractive:

                     polys.append(Point(pt1[1], pt1[0]))

                 else:

                     polys.append(Point(pt2[1], pt2[0]))

     changes = gpd.GeoDataFrame(geometry=polys)

     return changes, bdiff

okb: 

  nbins: 200 

  extractor_method: "orb"

  image_func: "diff"

  extractor_kwargs: 

    nfeatures: 500

    scaleFactor: 1.1 

    nlevels: 1

okbm:

  extractor_method: "orb"

  image_func: "diff"

  extractor_kwargs: 

    nfeatures: 2000

    scaleFactor: 1.1

    nlevels: 1

  cluster_kwargs: 

    min_samples: 10 

    max_eps: 10 

    eps: .5 

    xi: .5

blob: 

  image_func: "diff_sq"

  subtractive: False

  min_sigma: .45 

  max_sigma: 30

  num_sigma: 10 

  threshold: .4

  overlap: .5

  log_scale: False 

  exclude_border: False 

  n_neighbors: 3 

  dist_upper_bound: 5 

  angle_tolerance: 20   

jigsaw:

  # Parameters that seem to have the best results on Hirise pairs

  radius: 'yes'

  errorpropagation: 'yes'

  outlier_rejection: 'yes'

  point_longitude_sigma: 10 

  point_latitude_sigma: 10 

  point_radius_sigma: 2 

  maxits: 30

  sigma0: .1

  camsolve: "velocities"

  twist: "yes"

  overexisting: "yes"

  spsolve: "None"

        #"massaged" by the phase and template matcher.

        for i, group in df.groupby('point_id'):

            zero_group = group.iloc[0]

            apriori_geom = np.array(point_info(self.nodes[zero_group.image_index]['data'].geodata.file_name, zero_group.x, zero_group.y, 'image')['GroundPoint']['BodyFixedCoordinate'].value) * 1000

            apriori_geom = np.array(point_info(self.nodes[zero_group.image_index]['data'].geodata.file_name, zero_group.x, zero_group.y, 'image')['BodyFixedCoordinate'].value) * 1000

            for j, row in group.iterrows():

                row['aprioriX'] = apriori_geom[0]

                row['aprioriY'] = apriori_geom[1]

#!/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)