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

Explicit pointid, remove catch all try/excepts, adding isis point projection error handling (#439) 

* initial commit

* remove extra stuff

* removing explicit point id

* remove nullable requirement on overlapid

ground points are not guaranteed to be in overlaps, so remove this requirement
parent 5ed12c32

autocnet/io/db/model.py

+1 −0
Original line number Diff line number Diff line
@@ -274,6 +274,7 @@ class Points(BaseMixin, Base): 
    id = Column(Integer, primary_key=True, autoincrement=True)

    _pointtype = Column("pointType", IntEnum(PointType), nullable=False)  # 2, 3, 4 - Could be an enum in the future, map str to int in a decorator

    identifier = Column(String, unique=True)

    overlapid = Column(Integer, ForeignKey('overlay.id'))

    _geom = Column("geom", Geometry('POINT', srid=latitudinal_srid, dimension=2, spatial_index=True))

    cam_type = Column(String)

    ignore = Column("pointIgnore", Boolean, default=False)

autocnet/matcher/subpixel.py

+14 −11
Original line number Diff line number Diff line
@@ -11,6 +11,7 @@ from skimage import transform as tf 
from matplotlib import pyplot as plt



from plio.io.io_gdal import GeoDataset

from pysis.exceptions import ProcessError



from autocnet import Session, config

from autocnet.matcher.naive_template import pattern_match, pattern_match_autoreg
@@ -79,6 +80,9 @@ def check_image_size(imagesize): 
    imagesize : tuple

                in the form (size_x, size_y)

    """

    if isinstance(imagesize, int):

        imagesize = (imagesize, imagesize)



    x = imagesize[0]

    y = imagesize[1]
@@ -346,7 +350,7 @@ def subpixel_ciratefi(sx, sy, dx, dy, s_img, d_img, search_size=251, template_si 
    dy += (y_offset + t_roi.ayr)

    return dx, dy, strength



def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=(251, 251), reduction=11, convergence_threshold=1.0, max_dist=50, **kwargs):

def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=(51, 51), reduction=11, convergence_threshold=1.0, max_dist=50, **kwargs):

    """

    Iteratively apply a subpixel phase matcher to source (s_img) and destination (d_img)

    images. The size parameter is used to set the initial search space. The algorithm
@@ -396,10 +400,7 @@ def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=(251, 251), reduction=11, 
    dline = dy



    while True:

        try:

        shifted_dx, shifted_dy, metrics = subpixel_phase(sx, sy, dx, dy, s_img, d_img, image_size=size, **kwargs)

        except:

            return None, None, None



        # Compute the amount of move the matcher introduced

        delta_dx = abs(shifted_dx - dx)
@@ -449,8 +450,10 @@ 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

    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)

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



    match_points = [(base_startx,base_starty),

                    (base_startx,base_stopy),
@@ -459,12 +462,12 @@ def geom_match(base_cube, input_cube, bcenter_x, bcenter_y, size_x=60, size_y=60 


    cube_points = []

    for x,y in match_points:

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



    input_cube_extents = input_cube.raster_size

    for x,y in cube_points:

        if x < 0 or y < 0 or x > input_cube_extents[0] or y > input_cube_extents[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])
@@ -541,7 +544,7 @@ def geom_match(base_cube, input_cube, bcenter_x, bcenter_y, size_x=60, size_y=60 
      pcm = axs[2].imshow(corrmap**2, interpolation=None, cmap="coolwarm")

      plt.show()



    dist = np.linalg.norm([center_x-x, center_y-y])

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

    return sample, line, dist, maxcorr, corrmap
@@ -668,7 +671,7 @@ def subpixel_register_point(pointid, iterative_phase_kwargs={}, subpixel_templat 
            continue



        # Update the measure

        if new_template_x:

        if new_x:

            measure.sample = new_x

            measure.line = new_y

            measure.weight = cost

autocnet/spatial/isis.py

+14 −0
Original line number Diff line number Diff line
@@ -85,10 +85,18 @@ def point_info(cube_path, x, y, point_type, allow_outside=False): 
         dictres = []

         if len(x) > 1 and len(y) > 1:

            for r in pvlres:

                if r['GroundPoint']['Error'] is not None:

                    raise ProcessError(returncode=1, cmd=['pysis.campt()'], stdout=r, stderr=r['GroundPoint']['Error'])

                    return

                else:

                    # convert all pixels to PLIO pixels from ISIS

                    r[1]["Sample"] -= .5

                    r[1]["Line"] -= .5

                    dictres.append(dict(r[1]))

         else:

            if pvlres['GroundPoint']['Error'] is not None:

                raise ProcessError(returncode=1, cmd=['pysis.campt()'], stdout=pvlres, stderr=pvlres['GroundPoint']['Error'])

                return

            else:

                pvlres["GroundPoint"]["Sample"] -= .5

                pvlres["GroundPoint"]["Line"] -= .5
@@ -109,7 +117,10 @@ def image_to_ground(cube_path, sample, line, lattype="PlanetocentricLatitude", l 
           1-D array of longitudes or single floating point longitude



    """

    try:

        res = point_info(cube_path, sample, line, "image")

    except ProcessError as e:

        raise ProcessError(returncode=e.returncode, cmd=e.cmd, stdout=e.stdout, stderr=e.stderr)



    try:

        if isinstance(res, list):
@@ -138,7 +149,10 @@ def ground_to_image(cube_path, lon, lat): 
              array of samples or single dloating point sample



    """

    try:

        res = point_info(cube_path, lon, lat, "ground")

    except ProcessError as e:

        raise ProcessError(returncode=e.returncode, cmd=e.cmd, stdout=e.stdout, stderr=e.stderr)



    try:

        if isinstance(res, list):

autocnet/spatial/overlap.py

+23 −6
Original line number Diff line number Diff line
@@ -7,6 +7,7 @@ import pyproj 
import shapely

import sqlalchemy

from plio.io.io_gdal import GeoDataset

from pysis.exceptions import ProcessError



from autocnet import config, dem, Session

from autocnet.cg import cg as compgeom
@@ -19,7 +20,6 @@ from autocnet.transformation import roi 
from plurmy import Slurm

import csmapi





# SQL query to decompose pairwise overlaps

compute_overlaps_sql = """

WITH intersectiongeom AS
@@ -64,8 +64,9 @@ def place_points_in_overlaps(nodes, size_threshold=0.0007, 
        if overlaps == None:

            continue



        oid = o.id

        overlapnodes = [nodes[id]["data"] for id in overlaps]

        points.extend(place_points_in_overlap(overlapnodes, o.geom, cam_type=cam_type,

        points.extend(place_points_in_overlap(oid, overlapnodes, o.geom, cam_type=cam_type,

                                              distribute_points_kwargs=distribute_points_kwargs))

    Points.bulkadd(points)
@@ -127,7 +128,7 @@ def cluster_place_points_in_overlaps(size_threshold=0.0007, 
    submitter.submit(array='1-{}%24'.format(job_counter), chunksize=chunksize, exclude=exclude)

    return job_counter



def place_points_in_overlap(nodes, geom, cam_type="csm",

def place_points_in_overlap(oid, nodes, geom, cam_type="csm",

                            size=71,

                            distribute_points_kwargs={}):

    """
@@ -190,7 +191,12 @@ def place_points_in_overlap(nodes, geom, cam_type="csm", 
            # Convert to geocentric lon, lat

            geocent_lon, geocent_lat, _ = reproject([x, y, z],

                                                    semi_major, semi_major, 'geocent', 'latlon')

            try:

                line, sample = isis.ground_to_image(node["image_path"], geocent_lon ,geocent_lat)

            except ProcessError as e:

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

                    print(f'point ({geocent_lon}, {geocent_lat}) does not project to reference image {node["image_path"]}')

                    continue

        if cam_type == "csm":

            # The CSM conversion makes the LLA/ECEF conversion explicit

            gnd = csmapi.EcefCoord(x, y, z)
@@ -203,7 +209,6 @@ def place_points_in_overlap(nodes, geom, cam_type="csm", 
        try:

            interesting = extract_most_interesting(image)

        except:

            warnings.warn('Could not find an interesting feature around point')

            continue



        # kps are in the image space with upper left origin and the roi
@@ -215,11 +220,17 @@ def place_points_in_overlap(nodes, geom, cam_type="csm", 


        # Get the updated lat/lon from the feature in the node

        if cam_type == "isis":

            try:

                p = isis.point_info(node["image_path"], newsample, newline, point_type="image")

            except ProcessError as e:

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

                    print(node["image_path"])

                    print(f'interesting point ({newsample}, {newline}) does not project back to ground')

                    continue

            try:

                x, y, z = p["BodyFixedCoordinate"].value

            except:

                x,y,x = ["BodyFixedCoordinate"]

                x, y, z = ["BodyFixedCoordinate"]



            if getattr(p["BodyFixedCoordinate"], "units", "None").lower() == "km":

                x = x * 1000
@@ -255,7 +266,8 @@ def place_points_in_overlap(nodes, geom, cam_type="csm", 
                                                                 'geocent', 'latlon')



        point_geom = shapely.geometry.Point(x, y, z)

        point = Points(apriori=point_geom,

        point = Points(overlapid=oid,

                       apriori=point_geom,

                       adjusted=point_geom,

                       pointtype=2, # Would be 3 or 4 for ground

                       cam_type=cam_type)
@@ -269,7 +281,12 @@ def place_points_in_overlap(nodes, geom, cam_type="csm", 
                image_coord = node.camera.groundToImage(gnd)

                sample, line = image_coord.samp, image_coord.line

            if cam_type == "isis":

                try:

                    line, sample = isis.ground_to_image(node["image_path"], geocent_lon, geocent_lat)

                except ProcessError as e:

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

                        print(f'interesting point ({geocent_lon},{geocent_lat}) does not project to image {node["image_path"]}')

                        continue



            point.measures.append(Measures(sample=sample,

                                           line=line,

bin/acn_overlaps

+1 −1
Original line number Diff line number Diff line
@@ -39,7 +39,7 @@ def main(msg, config): 
        nodes.append(NetworkNode(node_id=id, image_path=res.path))

    session.close()

    

    points = place_points_in_overlap(nodes, geom, cam_type=msg["cam_type"],

    points = place_points_in_overlap(oid, nodes, geom, cam_type=msg["cam_type"],

                                     distribute_points_kwargs=msg['distribute_points_kwargs'])



    print('Adding {} points to the database.'.format(len(points)))