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

Necessary fixes discovered during grounded product creation (#452) 

* bulk commit, may need some cleaning

* removing prints, finishing doc string

* actually finished doc string

* addressing comments

* addressing last comment in first round

* fixing tests pt1

* fix test pt2

* test fix pt 3
parent 204c02a9

autocnet/cg/cg.py

+14 −25
Original line number Diff line number Diff line
@@ -7,6 +7,7 @@ import networkx as nx 
import geopandas as gpd

import ogr



from skimage import transform as tf

from scipy.spatial import ConvexHull

from scipy.spatial import Voronoi

import shapely.geometry
@@ -456,43 +457,31 @@ def distribute_points_new(geom, nspts, ewpts, Session): 
    valid : list

            of point coordinates in the form [(x1,y1), (x2,y2), ..., (xn, yn)]

    """

    geom_coords = np.column_stack(geom.exterior.xy)



    coords = np.array(list(zip(*geom.envelope.exterior.xy))[:-1])



    ll = coords[0]

    lr = coords[1]

    ur = coords[2]

    ul = coords[3]



    # Find the points nearest the ur and ll aligned // with eastern side of ground_poly

    elon, elat = find_side('east', Session)

    ur_actual = np.array(two_point_extrapolate(ur[1], elat, elon))[::-1]

    lr_actual = np.array(two_point_extrapolate(lr[1],elat, elon))[::-1]



    wlon, wlat = find_side('west', Session)

    ul_actual = np.array(two_point_extrapolate(ul[1], wlat, wlon))[::-1]

    ll_actual = np.array(two_point_extrapolate(ll[1], wlat, wlon))[::-1]

    rr_coords = np.array(list(zip(*geom.minimum_rotated_rectangle.exterior.xy))[:-1])

    w = min([i[0] for i in rr_coords])

    s = min([i[1] for i in rr_coords])

    swid = nearest([w, s], rr_coords)

    rr_coords = np.vstack([rr_coords[swid:], rr_coords[0:swid]]) # reorder to match envelope/coords order



    dt = (ur_actual-ul_actual)*0.025 #some offset to make sure endpoints are within geom

    db = (lr_actual-ll_actual)*0.025

    newtop = create_points_along_line(ul_actual+dt, ur_actual-dt, ewpts)

    newbot = create_points_along_line(ll_actual+db, lr_actual-db, ewpts)

    x = np.linspace(ul[0], ur[0], ewpts+2)[1:-1]

    y = np.linspace(ul[1], ll[1], nspts+2)[1:-1]



    points = []

    for i in range(len(newtop)):

        top = newtop[i]

        bot = newbot[i]

    grid = np.transpose([np.tile(x, len(y)), np.repeat(y, len(x))])



        line_of_points = create_points_along_line(top, bot, nspts)

        points.append(line_of_points)



    if len(points) < 1:

    if len(grid) < 1:

        return []



    points = np.vstack(points)

    affine = tf.estimate_transform('affine', coords, rr_coords)

    rr_grid = affine(grid)



    # Perform a spatial intersection check to eject points that are not valid

    valid = [p for p in points if xy_in_polygon(p[0], p[1], geom)]

    valid = [p for p in rr_grid if xy_in_polygon(p[0], p[1], geom)]

    return valid



def distribute_points_in_geom(geom, method="classic",

autocnet/graph/network.py

+2 −1
Original line number Diff line number Diff line
@@ -1756,7 +1756,7 @@ ORDER BY measures."pointid", measures."id"; 
        cnet.to_isis(df, path, targetname=target)

        cnet.write_filelist(fpaths, path=flistpath)



    def update_from_jigsaw(self, session, path):

    def update_from_jigsaw(self,path):

        """

        Updates the measures table in the database with data from

        a jigsaw bundle adjust
@@ -1784,6 +1784,7 @@ ORDER BY measures."pointid", measures."id"; 
        DROP TABLE temp_measures;

        """



        with self.session_scope() as session:

            session.execute(sql)

            session.commit()

autocnet/io/db/tests/test_model.py

+11 −9
Original line number Diff line number Diff line
@@ -160,18 +160,20 @@ def test_measures_exists(tables): 
                                                                                     'adjustedCovar': [[]],

                                                                                     'apriorisample': [0],

                                                                                     'aprioriline': [0]}))

def test_jigsaw_append(mockFunc, session, measure_data, point_data, image_data, ncg):

def test_jigsaw_append(mockFunc, measure_data, point_data, image_data, ncg):

    with ncg.session_scope() as session:

        model.Images.create(session, **image_data)

        model.Points.create(session, **point_data)

        model.Measures.create(session, **measure_data)

    resp = session.query(model.Measures).filter(model.Measures.id == 1).first()

    assert resp.liner == None

    assert resp.sampler == None



    ncg.update_from_jigsaw(session, '/Some/Path/To/An/ISISNetwork.cnet')

    resp = session.query(model.Measures).filter(model.Measures.id == 1).first()

    assert resp.liner == 0.1

    assert resp.sampler == 0.1

        resp1 = session.query(model.Measures).filter(model.Measures.id == 1).first()

        assert resp1.liner == None

        assert resp1.sampler == None



    ncg.update_from_jigsaw('/Some/Path/To/An/ISISNetwork.cnet')

    with ncg.session_scope() as session:

        resp2 = session.query(model.Measures).filter(model.Measures.id == 1).first()

        assert resp2.liner == 0.1

        assert resp2.sampler == 0.1



def test_null_footprint(session):

    i = model.Images.create(session, geom=None,

autocnet/matcher/cross_instrument_matcher.py

+84 −28
Original line number Diff line number Diff line
@@ -12,10 +12,8 @@ import os.path 
import socket

from ctypes.util import find_library



import numpy as np

import pandas as pd

import scipy

from matplotlib import pyplot as plt

from sqlalchemy import (Boolean, Column, Float, ForeignKey, Integer,

                        LargeBinary, String, UniqueConstraint, create_engine,

                        event, orm, pool)
@@ -38,7 +36,6 @@ from geoalchemy2 import functions 
from knoten import csm



from plio.io.io_controlnetwork import from_isis, to_isis

from plio.io.io_gdal import GeoDataset



from shapely import wkt

from shapely.geometry.multipolygon import MultiPolygon
@@ -47,17 +44,18 @@ from shapely.geometry import Point 
from plurmy import Slurm



from autocnet.io.db.model import Images, Points, Measures, JsonEncoder

from autocnet.cg.cg import distribute_points_in_geom

from autocnet.cg.cg import distribute_points_in_geom, xy_in_polygon

from autocnet.io.db.connection import new_connection

from autocnet.spatial import isis

from autocnet.transformation.spatial import reproject

from autocnet.matcher.cpu_extractor import extract_most_interesting

from autocnet.transformation import roi

from autocnet.matcher.subpixel import geom_match

from autocnet.utils.utils import bytescale



import warnings



def generate_ground_points(Session, ground_mosaic, nspts_func=lambda x: int(round(x,1)*1), ewpts_func=lambda x: int(round(x,1)*4)):

def generate_ground_points(Session, ground_mosaic, nspts_func=lambda x: int(round(x,1)*1), ewpts_func=lambda x: int(round(x,1)*4), size=(100,100)):

    """



    Parameters
@@ -75,6 +73,10 @@ def generate_ground_points(Session, ground_mosaic, nspts_func=lambda x: int(roun 
    ewpts_func       : func

                       describes distribution of points along the east-west

                       edge of an overlap.



    size             : tuple of int

                       (size_x, size_y) maximum distances on either access point

                       can move when attempting to find an interesting feature.

    """



    if isinstance(ground_mosaic, str):
@@ -87,34 +89,36 @@ def generate_ground_points(Session, ground_mosaic, nspts_func=lambda x: int(roun 
    fp_poly = wkt.loads(session.query(functions.ST_AsText(functions.ST_Union(Images.geom))).one()[0])

    session.close()



    fp_poly_bounds = list(fp_poly.bounds)



    # just hard code queries to the mars database as it exists for now



    coords = distribute_points_in_geom(fp_poly, nspts_func=nspts_func, ewpts_func=ewpts_func, method="new")

    coords = distribute_points_in_geom(fp_poly, nspts_func=nspts_func, ewpts_func=ewpts_func, method="new", Session=Session)

    coords = np.asarray(coords)



    records = []

    old_coord_list = []

    coord_list = []

    lines = []

    samples = []

    newlines = []

    newsamples = []



    # throw out points not intersecting the ground reference images

    print('points to lay down: ', len(coords))

    for i, coord in enumerate(coords):

        # res = ground_session.execute(formated_sql)

        p = Point(*coord)

        print(f'point {i}'),



        linessamples = isis.point_info(ground_mosaic.file_name, p.x, p.y, 'ground')

        sample = linessamples.get('Sample')

        line = linessamples.get('Line')

        line = linessamples[0].get('Line')

        sample = linessamples[0].get('Sample')



        oldpoint = isis.point_info(ground_mosaic.file_name, sample, line, 'image')

        op = Point(oldpoint[0].get('PositiveEast360Longitude'),

                   oldpoint[0].get('PlanetocentricLatitude'))



        # hardcoded for themis for now

        size = 200



        image = roi.Roi(ground_mosaic, sample, line, size_x=size, size_y=size, dtype="uint64")

        image_roi = image.clip()

        image = roi.Roi(ground_mosaic, sample, line, size_x=size[0], size_y=size[1])

        image_roi = image.clip(dtype="uint64")



        interesting = extract_most_interesting(bytescale(image),  extractor_parameters={'nfeatures':30})

        interesting = extract_most_interesting(bytescale(image_roi),  extractor_parameters={'nfeatures':30})



        # kps are in the image space with upper left origin, so convert to

        # center origin and then convert back into full image space
@@ -123,21 +127,31 @@ def generate_ground_points(Session, ground_mosaic, nspts_func=lambda x: int(roun 
        newline = top_y + interesting.y



        newpoint = isis.point_info(ground_mosaic.file_name, newsample, newline, 'image')

        p = Point(newpoint.get('PositiveEast360Longitude'),

                  newpoint.get('PlanetocentricLatitude'))

        p = Point(newpoint[0].get('PositiveEast360Longitude'),

                  newpoint[0].get('PlanetocentricLatitude'))



        if not (xy_in_polygon(p.x, p.y, fp_poly)):

                print('Interesting point not in mosaic area, ignore')

                continue



        old_coord_list.append(op)

        lines.append(line)

        samples.append(sample)

        coord_list.append(p)

        lines.append(newline)

        samples.append(newsample)

        newlines.append(newline)

        newsamples.append(newsample)





    # start building the cnet

    ground_cnet = pd.DataFrame()

    ground_cnet["path"] = [ground_mosaic.file_name]*len(coord_list)

    ground_cnet["pointid"] = list(range(len(coord_list)))

    ground_cnet["original point"] = old_coord_list

    ground_cnet["point"] = coord_list

    ground_cnet['line'] = lines

    ground_cnet['sample'] = samples

    ground_cnet['original_line'] = lines

    ground_cnet['original_sample'] = samples

    ground_cnet['line'] = newlines

    ground_cnet['sample'] = newsamples

    ground_cnet = gpd.GeoDataFrame(ground_cnet, geometry='point')

    return ground_cnet, fp_poly, coord_list
@@ -151,12 +165,18 @@ def propagate_point(Session, 
                    paths,

                    lines,

                    samples,

                    size_x=40,

                    size_y=40,

                    template_kwargs={'image_size': (39, 39), 'template_size': (21, 21)},

                    verbose=False):

    """



    """

    engine = Session.get_bind()

    images = gpd.GeoDataFrame.from_postgis(f"select * from images where ST_Intersects(geom, ST_SetSRID(ST_Point({lon}, {lat}), {config['spatial']['latitudinal_srid']}))", engine, geom_col="geom")

    session = Session()

    engine = session.get_bind()

    string = f"select * from images where ST_Intersects(geom, ST_SetSRID(ST_Point({lon}, {lat}), {config['spatial']['latitudinal_srid']}))"

    images = pd.read_sql(string, engine)

    session.close()



    image_measures = pd.DataFrame(zip(paths, lines, samples), columns=["path", "line", "sample"])

    measure = image_measures.iloc[0]
@@ -177,7 +197,10 @@ def propagate_point(Session, 
            sx, sy = m["sample"], m["line"]



            try:

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

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

                        size_x=size_x, size_y=size_y, \

                        template_kwargs=template_kwargs, \

                        verbose=verbose)

            except Exception as e:

                raise Exception(e)

                match_results.append(e)
@@ -204,6 +227,7 @@ def propagate_point(Session, 
          print("Winning CORR: ", best_results[3], "Themis Pixel shift: ", best_results[4])

          print("Themis Image: ", best_results[6], "CTX image:", best_results[7])

          print("Themis S,L: ", f"{sx},{sy}", "CTX S,L: ", f"{sample},{line}")

          print('\n')



        # hardcoded for now

        if best_results[3] == None or best_results[3] < 0.7:
@@ -231,8 +255,37 @@ def propagate_point(Session, 


    return new_measures



def propagate_control_network(Session, config, dem, base_cnet, verbose=False):

def propagate_control_network(Session, config, dem, base_cnet,

        size_x=40, size_y=40,

        template_kwargs={'image_size': (39,39), 'template_size': (21,21)}, verbose=False):

    """

    Parameters

    ----------

    Session   : sqlalchemy session maker

                session maker associated with the database you want to propagate to



    config    : dict

                configuation file associated with database you want to propagate to

                In the form: {'username':'somename',

                              'password':'somepassword',

                              'host':'somehost',

                              'pgbouncer_port':6543,

                              'name':'somename'}



    dem       : plio.io.io_gdal.GeoDataset

                Digital elevation model of target body



    base_cnet : pd.DataFrame

                Dataframe representing the points you want to propagate. Must contain line, sample, path.



    verbose   : boolean

                Increase the level of print outs/plots recieved during propagation





    Output

    ------

    ground   : pd.DataFrame

               Dataframe containing successfully propagated points



    """

    warnings.warn('This function is not well tested. No tests currently exists \
@@ -261,6 +314,9 @@ def propagate_control_network(Session, config, dem, base_cnet, verbose=False): 
                                      measures["path"],

                                      measures["line"],

                                      measures["sample"],

                                      size_x,

                                      size_y,

                                      template_kwargs,

                                      verbose=verbose)



        constrained_net.extend(gp_measures)

autocnet/matcher/subpixel.py

+24 −20
Original line number Diff line number Diff line
@@ -151,7 +151,7 @@ def clip_roi(img, center_x, center_y, size_x=200, size_y=200, dtype="uint64"): 


def subpixel_phase(sx, sy, dx, dy,

                   s_img, d_img,

                   image_size=(251, 251),

                   image_size=(51, 51),

                   **kwargs):

    """

    Apply the spectral domain matcher to a search and template image. To
@@ -446,31 +446,30 @@ def geom_match(base_cube, input_cube, bcenter_x, bcenter_y, size_x=60, size_y=60 
    if base_starty < 0:

        raise Exception(f"Window: {base_starty} < 0, center: {bcenter_x},{bcenter_y}")



    # specifically not putting this in a try except, because this should never fail,

    # want to throw error if there is one

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



    match_points = [(base_startx,base_starty),

    base_corners = [(base_startx,base_starty),

                    (base_startx,base_stopy),

                    (base_stopx,base_stopy),

                    (base_stopx,base_starty)]



    cube_points = []

    for x,y in match_points:

    dst_corners = []

    for x,y in base_corners:

        try:

            lat, lon = spatial.isis.image_to_ground(base_cube.file_name, x, y)

            cube_points.append(spatial.isis.ground_to_image(input_cube.file_name, lon, lat)[::-1])

            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([*match_points])

    base_gcps = np.array([*base_corners])

    base_gcps[:,0] -= base_startx

    base_gcps[:,1] -= base_starty



    dst_gcps = np.array([*cube_points])

    dst_gcps = np.array([*dst_corners])

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

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

    stop_x = dst_gcps[:,0].max()
@@ -487,7 +486,7 @@ def geom_match(base_cube, input_cube, bcenter_x, bcenter_y, size_x=60, size_y=60 
            "Real" : "float64"

    }



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

    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)
@@ -512,25 +511,33 @@ def geom_match(base_cube, input_cube, bcenter_x, bcenter_y, size_x=60, size_y=60 
    restemplate = subpixel_template(size_x, size_y, size_x, size_y, bytescale(base_arr), bytescale(dst_arr), **template_kwargs)



    if phase_kwargs:

        resphase = subpixel_phase(size_x, size_y, restemplate[0], restemplate[1], base_arr, dst_arr, **phase_kwargs)

        _,_,maxcorr, temp_corrmap = restemplate

        sample_template, line_template = affine([restemplate[0], restemplate[1]])[0]

        sample_template += start_x

        line_template += start_y

        dist_temp = np.linalg.norm([center_x-sample_template, center_y-line_template])



        resphase = subpixel_phase(size_x, size_y, restemplate[0], restemplate[1], base_arr, dst_arr, **phase_kwargs)

        x,y,(perror, pdiff) = resphase

        if x is None or y is None:

            return None, None, None, None, None

        temp_dist = np.linalg.norm([size_x-restemplate[0], size_y-restemplate[1]])

        phase_dist = np.linalg.norm([restemplate[0]-resphase[0], restemplate[1]-resphase[1]])

        dist = (temp_dist, phase_dist)

        metric = (restemplate[2], perror, pdiff)

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

        sample += start_x

        line += start_y

        phase_dist = np.linalg.norm([sample_template-sample, line_template-line])



        dist = (dist_temp, dist_phase)

        metric = (maxcorr, perror, pdiff)

    else:

        x,y,maxcorr,temp_corrmap = restemplate

        if x is None or y is None:

            return None, None, None, None, None

        metric = maxcorr

        dist = np.linalg.norm([size_x/2-x, size_y/2-y])



        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)
@@ -547,7 +554,6 @@ def geom_match(base_cube, input_cube, bcenter_x, bcenter_y, size_x=60, size_y=60 
      pcm = axs[2].imshow(temp_corrmap**2, interpolation=None, cmap="coolwarm")

      plt.show()



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

    return sample, line, dist, metric, temp_corrmap
@@ -720,14 +726,12 @@ def subpixel_register_point(pointid, iterative_phase_kwargs={}, 
                continue



            if iterative_phase_kwargs:

                print('subpixel_register_point -> PHASE MEASURE WRITE OUT')

                measure.template_metric = metric[0]

                measure.template_shift = dist[0]

                measure.phase_error = metric[1]

                measure.phase_diff = metric[2]

                measure.phase_shift = dist[1]

            else:

                print('subpixel_register_point -> NO PHASE MEASURE WRITE OUT')

                measure.template_metric = metric

                measure.template_shift = dist