Commit 02fd6e4a authored by Kelvin Rodriguez's avatar Kelvin Rodriguez Committed by GitHub
Browse files

added hirise coreg + updated to_isis to match network version (#408) 

* added hirise coreg + updated to_isis to match network version

* updated matcher

* faster match, hirise update

* added okubo method

* things workin
parent fda712fa

autocnet/cg/change_detection.py

0 → 100644
+66 −0
Original line number Diff line number Diff line 


import numpy as np

from matplotlib.path import Path

from shapely.geometry import Point, MultiPoint

import geopandas as gpd



import cv2

from sklearn.cluster import  OPTICS



from autocnet.utils.utils import bytescale

from autocnet.matcher.cpu_extractor import extract_features



def image_diff(arr1, arr2):

    diff = arr1-arr2

    diff[np.isnan(diff)] = 0



    bdiff = bytescale(diff)

    return bdiff





def okubogar_detector(image1, image2, nbins=50, extractor_method="orb", extractor_kwargs={"nfeatures": 2000, "scaleFactor": 1.1, "nlevels": 1}, image_func=image_diff):

    arr1 = image1.read_array()

    arr2 = image2.read_array()

    arr1[arr1 == arr1.min()] = np.nan

    arr2[arr2 == arr2.min()] = np.nan



    bdiff = image_func(arr1, arr2)



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



    optics = OPTICS(min_samples=10, max_eps=20,  eps=.3, p=2, xi=.5).fit(list(zip(x,y)))



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

    classes["label"] = optics.labels_

    classes["point"] = points

    class_groups = classes.groupby("label").groups



    polys = []

    weights = []



    # array of x,y pairs

    xv, yv = np.mgrid[0:bdiff.shape[1], 0:bdiff.shape[0]]



    for label, indices in class_groups.items():

        if label == -1:

            continue



        points = classes.loc[indices]["point"]

        poly = MultiPoint(points.__array__()).convex_hull

        xmin, ymin, xmax, ymax = np.asarray(poly.bounds).astype("uint64")

        xv, yv = np.mgrid[xmin:xmax, ymin:ymax]

        xv = xv.flatten()

        yv = yv.flatten()



        points = np.vstack((xv,yv)).T.astype("uint64")



        mask = Path(np.asarray(poly.exterior.xy).T.astype("uint64")).contains_points(points).reshape(int(ymax-ymin), int(xmax-xmin))

        weight = bdiff[ymin:ymax,xmin:xmax].mean()



        polys.append(poly)

        weights.append(weight)



    return polys, weights

autocnet/graph/network.py

+38 −18
Original line number Diff line number Diff line
@@ -22,7 +22,7 @@ import shapely.ops 


import pyproj



from plio.io.io_controlnetwork import from_isis

from plio.io.io_controlnetwork import to_isis, from_isis

from plio.io import io_hdf, io_json

from plio.utils import utils as io_utils

from plio.io.io_gdal import GeoDataset
@@ -43,6 +43,7 @@ from autocnet.io.db.connection import new_connection, Parent 
from autocnet.vis.graph_view import plot_graph, cluster_plot

from autocnet.control import control

from autocnet.spatial.overlap import compute_overlaps_sql

from autocnet.spatial.isis import point_info



#np.warnings.filterwarnings('ignore')
@@ -1253,27 +1254,46 @@ class CandidateGraph(nx.Graph): 
        return self.controlnetwork.groupby('point_id').apply(lambda g: g if len(g) > 1 else None)





    def to_isis(self, outname, serials, olist, *args, **kwargs):  # pragma: no cover

    def to_isis(self, outname, flistpath=None, target="Mars"):  # pragma: no cover

        """

        Write the control network out to the ISIS3 control network format.

        """

        df = self.controlnetwork



        if self.validate_points().any() == True:

            warnings.warn(

                'Control Network is not ISIS3 compliant.  Please run the validate_points method on the control network.')

            return

        serials = [generate_serial_number(self.nodes[id_]["data"]["image_path"]) for id_ in df["image_index"]]



        #create columns in the dataframe; zeros ensure plio (/protobuf) will

        #ignore unless populated with alternate values

        df['aprioriX'] = 0

        df['aprioriY'] = 0

        df['aprioriZ'] = 0

        df['adjustedX'] = 0

        df['adjustedY'] = 0

        df['adjustedZ'] = 0

        df['type'] = 3

        df['measureType'] = 2



        df["serialnumber"] = serials



        # Apply the subpixel shift

        self.controlnetwork.x += self.controlnetwork.x_off

        self.controlnetwork.y += self.controlnetwork.y_off

        #only populate the new columns for ground points. Otherwise, isis will

        #recalculate the control point lat/lon from control measures which where

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

            for j, row in group.iterrows():

                row['aprioriX'] = apriori_geom[0]

                row['aprioriY'] = apriori_geom[1]

                row['aprioriZ'] = apriori_geom[2]

                df.iloc[row.name] = row



        to_isis(outname + '.net', self.controlnetwork.query('valid == True'),

                serials, *args, **kwargs)

        write_filelist(olist, outname + '.lis')

        if flistpath is None:

            flistpath = os.path.splitext(outname)[0] + '.lis'



        # Back out the subpixel shift

        self.controlnetwork.x -= self.controlnetwork.x_off

        self.controlnetwork.y -= self.controlnetwork.y_off

        df = df.rename(columns={'image_index':'image_id','point_id':'id', 'type' : 'pointType',

             'x':'sample', 'y':'line'})

        cnet.to_isis(df, outname, targetname=target)

        cnet.write_filelist(self.files, path=flistpath)



    def to_bal(self):

        """
@@ -1699,9 +1719,9 @@ WHERE 
        # Get the camera objects to manually join. Keeps the caller from

        # having to remember to bring cameras as well.

        ids = [i[0] for i in sourceimages]

        cameras = sourcesession.query(Cameras).filter(Cameras.image_id.in_(ids)).all()

        for c in cameras:

            destinationsession.merge(c)

        #cameras = sourcesession.query(Cameras).filter(Cameras.image_id.in_(ids)).all()

        #for c in cameras:

        #    destinationsession.merge(c)



        destinationsession.commit()

        destinationsession.close()

autocnet/matcher/cpu_matcher.py

+11 −0
Original line number Diff line number Diff line
@@ -3,6 +3,8 @@ import warnings 
import numpy as np

import pandas as pd



import cv2



FLANN_INDEX_KDTREE = 1  # Algorithm to set centers,

DEFAULT_FLANN_PARAMETERS = dict(algorithm=FLANN_INDEX_KDTREE, trees=3)
@@ -84,9 +86,18 @@ def match(edge, k=2, **kwargs): 
    # Swap the indices since mono_matches is generic and source/destin are

    # swapped

    mono_matches(edge.destination, edge.source, aidx=bidx, bidx=aidx)



    source_keypoints = edge.source.keypoints[['x', 'y']]

    source_keypoints.rename(columns={'x': 'source_x', 'y': 'source_y'}, inplace=True)

    edge.matches = edge.matches.join(source_keypoints, 'source_idx')



    destination_keypoints = edge.destination.keypoints[['x', 'y']]

    destination_keypoints.rename(columns={'x': 'destination_x', 'y': 'destination_y'}, inplace=True)

    edge.matches = edge.matches.join(destination_keypoints, 'destination_idx')

    edge.matches.sort_values(by=['distance'])







class FlannMatcher(object):

    """

    A wrapper to the OpenCV Flann based matcher class that adds

autocnet/utils/hirise.py

0 → 100644
+114 −0
Original line number Diff line number Diff line 
import os

from glob import glob

import geopandas as gpd



from pysis import isis

from pysis.exceptions import ProcessError



import plio

from plio.io.io_gdal import GeoDataset



from shapely import wkt

import numpy as np

import pvl

from shapely.geometry import Point, MultiPolygon



def segment_hirise(directory, offset=300):

    images = glob(os.path.join(directory, "*RED*.stitched.norm.cub"))

    for image in images:

        label = pvl.loads(isis.catlab(from_=image))



        dims = label["IsisCube"]["Core"]["Dimensions"]

        nlines, nsamples = dims["Lines"], dims["Samples"]

        print("Lines, Samples: ", nlines, nsamples)



        starts = np.arange(1, nlines, nsamples)

        stops = np.append(np.arange(starts[1], nlines, nsamples), [nlines])



        starts[1:] -= offset

        stops[:-1] += offset



        segments = np.asarray([starts, stops]).T



        for i, seg in enumerate(segments):

            start, stop = seg

            output = os.path.splitext(image)[0] + f".{start}_{stop}" + ".cub"

            print("Writing:", output)

            isis.crop(from_=image, to=output, line=start, nlines=stop-start, sample=1, nsamples=nsamples)

            isis.footprintinit(from_=output)



    return load_segments(directory)





def load_segments(directory):

    images = glob(os.path.join(directory, "*RED*.*_*.cub"))

    objs = [GeoDataset(image) for image in images]

    footprints = [o.footprint for o in objs]

    footprints = [wkt.loads(f.ExportToWkt()) for f in footprints]

    return gpd.GeoDataFrame(data=np.asarray([images, objs, footprints]).T, columns=["path", "image", "footprint"], geometry="footprint")





def ingest_hirise(directory):



    l = glob(os.path.join(directory, "*RED*.IMG"))

    l = [os.path.splitext(i)[0] for i in l]

    print(l)

    cube_name = "_".join(os.path.splitext(os.path.basename(l[0]))[0].split("_")[:-2])



    print("Cube Name:", cube_name)



    print(f"Running hi2isis on {l}")

    for i,cube in enumerate(l):

        print(f"{i+1}/{len(l)}")

        try:

            isis.hi2isis(from_=f'{cube}.IMG', to=f"{cube}.cub")

            print(f"finished {cube}")

        except ProcessError as e:

            print(e.stderr)

            return



    print(f"running spiceinit on {l}")

    for i,cube in enumerate(l):

        print(f"{i+1}/{len(l)}")

        try:

            isis.spiceinit(from_=f'{cube}.cub')

        except ProcessError as e:

            print(e.stderr)

            return



    print(f"running hical on {l}")

    for i,cube in enumerate(l):

        print(f"{i}/{len(l)}")

        try:

            isis.hical(from_=f'{cube}.cub', to=f'{cube}.cal.cub')

        except ProcessError as e:

            print(e.stderr)

            return



    cal_list_0 = sorted(glob(os.path.join(directory, "*0.cal*")))

    cal_list_1 = sorted(glob(os.path.join(directory, "*1.cal*")))

    print(f"Channel 0 images: {cal_list_0}")

    print(f"Channel 1 images: {cal_list_1}")



    for i,cubes in enumerate(zip(cal_list_0, cal_list_1)):

        print(f"{i+1}/{len(cal_list_0)}")

        c0, c1 = cubes

        output ="_".join(c0.split("_")[:-1])

        try:

            isis.histitch(from1=c0, from2=c1, to=f"{output}.stitched.cub")

        except ProcessError as e:

            print(e.stderr)

            return



    stitch_list = glob(os.path.join(directory, "*stitched*"))

    for cube in stitch_list:

        output = os.path.splitext(cube)[0] + ".norm.cub"

        try:

            isis.cubenorm(from_=cube, to=output)

        except ProcessError as e:

            print(e.stderr)

            return

autocnet/vis/graph_view.py

+1 −2
Original line number Diff line number Diff line
@@ -9,8 +9,7 @@ from skimage.transform import resize 
def downsample(array, amount):

    return resize(array,

                      (int(array.shape[0] / amount),

                      int(array.shape[1] / amount)),

                      interp='lanczos')

                      int(array.shape[1] / amount)))



def plot_graph(graph, ax=None, cmap='Spectral', labels=False, font_size=12, clusters=None, **kwargs):

    """