Unverified Commit e386d11c authored by jlaura's avatar jlaura Committed by GitHub
Browse files

Jlaura dev (#418)



* Adds subpixel measure registration where a measure

* adds passing extractor parameters to ferature extractors

* Adds an ORB extractor to placing points in overlaps

* setting default camtype to isis for tests

* fixes test failure

* swaps phase template

* think I missed a typo and save

* branches are better - fies bad import

* fixing two changes in one PR

* fixing poor form syntax issues

* Fixes the overlap to support CSM

* Modularizing to make testable

* Mocks added

* Param for size

* Mocks ficxed

* added if csm/isis

* updated tests

Co-authored-by: default avatarjlaura <jlaura@usgs.gov>
parent 2275751f
Loading
Loading
Loading
Loading
+33 −1
Original line number Diff line number Diff line
@@ -59,7 +59,8 @@ def extract_features(array, extractor_method='sift', extractor_parameters={}):
    if  extractor_method == 'vlfeat':
        keypoint_objs, descriptors  = vl.sift.sift(array,
                                                   compute_descriptor=True,
                                                   float_descriptors=True)
                                                   float_descriptors=True,
                                                   **extractor_parameters)
        # Swap columns for value style access, vl_feat returns y, x
        keypoint_objs[:, 0], keypoint_objs[:, 1] = keypoint_objs[:, 1], keypoint_objs[:, 0].copy()
        keypoints = pd.DataFrame(keypoint_objs, columns=['x', 'y', 'size', 'angle'])
@@ -86,3 +87,34 @@ def extract_features(array, extractor_method='sift', extractor_parameters={}):
            descriptors = descriptors.astype(np.float32)

    return keypoints, descriptors

def extract_most_interesting(image, extractor_method='orb', extractor_parameters={'nfeatures':10}):
    """
    Given an image, extract the most interesting feature. Interesting is defined
    as the feature descriptor that has the maximum variance. By default, this func
    finds 10 features in the image and then selects the best.

    Parameters
    ----------
    image : ndarray
            of DN values
    
    extractor_method : str
                       Any valid, autocnet extractor. Default (orb)

    exctractor_parameters : dict
                            of extractor parameters passed through to the feature extractor

    Returns
    -------
     : pd.series
       The keypoints row with the higest variance. The row has 'x' and 'y' columns to 
       get the location.
    """
    kps, desc = extract_features(image,
                                 extractor_method=extractor_method,
                                 extractor_parameters=extractor_parameters)
    
    # Naively assume that the maximum variance is the most unique feature
    vari = np.var(desc, axis=1)
    return kps.iloc[np.argmax(vari)] 
 No newline at end of file
+110 −4
Original line number Diff line number Diff line
@@ -7,7 +7,7 @@ from redis import StrictRedis
from plurmy import Slurm

from autocnet import Session, config
from autocnet.matcher import naive_template
from autocnet.matcher.naive_template import pattern_match
from autocnet.matcher import ciratefi
from autocnet.io.db.model import Measures, Points, Images, JsonEncoder
from autocnet.graph.node import NetworkNode
@@ -212,6 +212,7 @@ def subpixel_template(sx, sy, dx, dy, s_img, d_img, image_size=(251, 251), templ
    See Also
    --------
    autocnet.matcher.naive_template.pattern_match : for the kwargs that can be passed to the matcher
    autocnet.matcher.naive_template.pattern_match_autoreg : for the jwargs that can be passed to the autoreg style matcher
    """

    image_size = check_image_size(image_size)
@@ -223,7 +224,7 @@ def subpixel_template(sx, sy, dx, dy, s_img, d_img, image_size=(251, 251), templ
    if (s_image is None) or (d_template is None):
        return None, None, None

    shift_x, shift_y, metrics = naive_template.pattern_match(d_template, s_image, **kwargs)
    shift_x, shift_y, metrics = pattern_match(d_template, s_image, **kwargs)

    dx = (dx - shift_x + dxr)
    dy = (dy - shift_y + dyr)
@@ -362,7 +363,6 @@ def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=251, reduction=11, conver
        # Apply the shift to d_search and compute the new correspondence location
        dx += shift_x  # The implementation already applies the dxr, dyr shifts
        dy += shift_y 

        # Break if the solution has converged
        size = (size[0] - reduction, size[1] - reduction)

@@ -375,6 +375,69 @@ def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=251, reduction=11, conver
            break
    return dx, dy, metrics

def subpixel_register_measure(measureid, iterative_phase_kwargs={}, subpixel_template_kwargs={},
                            cost_func=lambda x,y: 1/x**2 * y, threshold=0.005):

    session = Session()
    
    # Setup the measure that is going to be matched
    destination = session.query(Measures).filter(Measures.id == measureid).one()
    destinationid = destination.imageid
    res = session.query(Images).filter(Images.id == destinationid).one()
    destination_node = NetworkNode(node_id=destinationid, image_path=res.path)

    # Get the point id and set up the reference measure
    pointid = destination.pointid
    source = session.query(Measures).filter(Measures.pointid==pointid).order_by(Measures.id).first()
    sourceid = source.imageid
    res = session.query(Images).filter(Images.id == sourceid).one()
    source_node = NetworkNode(node_id=sourceid, image_path=res.path)

    new_template_x, new_template_y, template_metric = subpixel_template(source.sample,
                                                            source.line,
                                                            destination.sample,
                                                            destination.line,
                                                            source_node.geodata,
                                                            destination_node.geodata,
                                                            **subpixel_template_kwargs)
    if new_template_x == None:
        destination.ignore = True # Unable to template match
        return

    new_phase_x, new_phase_y, phase_metrics = iterative_phase(source.sample,
                                                                source.line,
                                                                new_template_x,
                                                                new_template_y,
                                                                source_node.geodata,
                                                                destination_node.geodata,
                                                                **iterative_phase_kwargs)
    if new_phase_x == None:
        destination.ignore = True # Unable to phase match
        return

    dist = np.linalg.norm([new_phase_x-new_template_x, new_phase_y-new_template_y])
    cost = cost_func(dist, template_metric)

    if cost <= threshold:
        destination.ignore = True # Threshold criteria not met
        return

    # Update the measure
    if new_template_x:
        destination.sample = new_template_x
        destination.line = new_template_y
        destination.weight = cost

    # In case this is a second run, set the ignore to False if this
    # measures passed. Also, set the source measure back to ignore=False
    destination.ignore = False
    source.ignore = False

    session.commit()
    session.close()



def subpixel_register_point(pointid, iterative_phase_kwargs={}, subpixel_template_kwargs={},
                            cost_func=lambda x,y: 1/x**2 * y, threshold=0.005):

@@ -571,3 +634,46 @@ def cluster_subpixel_register_points(iterative_phase_kwargs={'size': 251},
                 output=config['cluster']['cluster_log_dir']+f'/autocnet.subpixel_register-%j')
    submitter.submit(array='1-{}'.format(job_counter), chunksize=chunksize, exclude=exclude)
    return job_counter

def cluster_subpixel_register_measures(iterative_phase_kwargs={'size': 251},
                                     subpixel_template_kwargs={'image_size':(251,251)},
                                     cost_kwargs={},
                                     threshold=0.005,
                                     filters={},
                                     walltime='00:10:00',
                                     chunksize=1000,
                                     exclude=None):
    # Setup the redis queue
    rqueue = StrictRedis(host=config['redis']['host'],
                        port=config['redis']['port'],
                        db=0)

    # Push the job messages onto the queue
    queuename = config['redis']['processing_queue']

    session = Session()
    query = session.query(Measures)
    for attr, value in filters.items():
        query = query.filter(getattr(Measures, attr)==value)
    res = query.all()
    for i, measure in enumerate(res):
        msg = {'id' : measure.id,
               'iterative_phase_kwargs' : iterative_phase_kwargs,
               'subpixel_template_kwargs' : subpixel_template_kwargs,
               'threshold':threshold,
               'cost_kwargs': cost_kwargs,
               'walltime' : walltime}
        rqueue.rpush(queuename, json.dumps(msg, cls=JsonEncoder))
    session.close()

    job_counter = i + 1

    # Submit the jobs
    submitter = Slurm('acn_subpixel_measure',
                 job_name='subpixel_register_measure',
                 mem_per_cpu=config['cluster']['processing_memory'],
                 time=walltime,
                 partition=config['cluster']['queue'],
                 output=config['cluster']['cluster_log_dir']+f'/autocnet.subpixel_register-%j')
    submitter.submit(array='1-{}'.format(job_counter), chunksize=chunksize, exclude=exclude)
    return job_counter
+54 −3
Original line number Diff line number Diff line
@@ -2,6 +2,7 @@ import warnings
import json

from redis import StrictRedis
import numpy as np
import pyproj
import shapely
import sqlalchemy
@@ -11,7 +12,10 @@ from autocnet import config, dem, Session
from autocnet.cg import cg as compgeom
from autocnet.io.db.model import Images, Measures, Overlay, Points, JsonEncoder
from autocnet.spatial import isis
from autocnet.matcher.subpixel import clip_roi
from autocnet.matcher.cpu_extractor import extract_most_interesting
from autocnet.transformation.spatial import reproject

from plurmy import Slurm
import csmapi

@@ -124,6 +128,7 @@ def cluster_place_points_in_overlaps(size_threshold=0.0007,
    return job_counter

def place_points_in_overlap(nodes, geom, cam_type="csm",
                            size=71,
                            distribute_points_kwargs={}):
    """
    Place points into an overlap geometry by back-projecing using sensor models.
@@ -141,6 +146,10 @@ def place_points_in_overlap(nodes, geom, cam_type="csm",
               options: {"csm", "isis"}
               Pick what kind of camera model implementation to use

    size : int
           The size of the window used to extractor features to find an
           interesting feature to which the point is shifted.

    Returns
    -------
    points : list of Points
@@ -162,6 +171,7 @@ def place_points_in_overlap(nodes, geom, cam_type="csm",
    for v in valid:
        lon = v[0]
        lat = v[1]

        # Calculate the height, the distance (in meters) above or
        # below the aeroid (meters above or below the BCBF spheroid).
        if dem is None:
@@ -170,6 +180,47 @@ def place_points_in_overlap(nodes, geom, cam_type="csm",
            px, py = dem.latlon_to_pixel(lat, lon)
            height = dem.read_array(1, [px, py, 1, 1])[0][0]

        # Need to get the first node and then convert from lat/lon to image space
        node = nodes[0]
        if cam_type == "isis":
            line, sample = isis.ground_to_image(node["image_path"], lon ,lat)
        if cam_type == "csm":
            # The CSM conversion makes the LLA/ECEF conversion explicit
            x, y, z = reproject([lon, lat, height],
                                 semi_major, semi_minor,
                                 'latlon', 'geocent')
            gnd = csmapi.EcefCoord(x, y, z)
            image_coord = node.camera.groundToImage(gnd)
            sample, line = image_coord.samp, image_coord.line

        # Extract ORB features in a sub-image around the desired point
        image, _, _ = clip_roi(node.geodata, sample, line, size_x=size, size_y=size)
        interesting = extract_most_interesting(image)

        # kps are in the image space with upper left origin, so convert to
        # center origin and then convert back into full image space
        newsample = sample + (interesting.x - size)
        newline = line + (interesting.y - size)

        # Get the updated lat/lon from the feature in the node
        if cam_type == "isis":
            p = isis.point_info(node["image_path"], newsample, newline, pointtype="image")
            x, y, z = p["GroundPoint"]["BodyFixedCoordinate"].value
        elif cam_type == "csm":
            image_coord = csmapi.ImageCoord(newline, newsample)
            pcoord = node.camera.imageToGround(image_coord)
            # Get the BCEF coordinate from the lon, lat
            lon, lat, _ = reproject([pcoord.x, pcoord.y, pcoord.z], semi_major, semi_minor,
                            'geocent', 'latlon')

            # Get the new DEM height
            if dem is None:
                height = 0
            else:
                px, py = dem.latlon_to_pixel(lat, lon)
                height = dem.read_array(1, [px, py, 1, 1])[0][0]


            # Get the BCEF coordinate from the lon, lat
            x, y, z = reproject([lon, lat, height], semi_major, semi_minor,
                                'latlon', 'geocent')
+18 −6
Original line number Diff line number Diff line
from collections import namedtuple
import pytest
from unittest.mock import MagicMock, patch
import numpy as np
from shapely.geometry import Polygon

from autocnet.spatial.overlap import place_points_in_overlap, place_points_in_overlaps
from autocnet.graph.node import Node
import csmapi

MockKeypoints = namedtuple('Keypoints', ['x', 'y'])
mockkeypoints = MockKeypoints(0,0)

@patch('autocnet.cg.cg.distribute_points_in_geom', return_value=[(0, 0), (5, 5), (10, 10)])
def test_place_points_in_overlap(point_distributer):
@patch('autocnet.spatial.overlap.clip_roi', return_value=np.zeros((3,3)))
@patch('autocnet.spatial.overlap.extract_most_interesting', return_value=mockkeypoints)
def test_place_points_in_overlap(point_distributer, clip_roi, extractor):
    # Mock setup
    first_node = MagicMock()
    first_node.camera = MagicMock()
    first_node.camera.groundToImage.return_value = csmapi.ImageCoord(1.0, 0.0)
    first_node.camera.imageToGround.return_value = csmapi.EcefCoord(1.0,1.0,1.0)
    first_node.isis_serial = '1'
    first_node.__getitem__.return_value = 1
    second_node = MagicMock()
    second_node.camera = MagicMock()
    second_node.camera.groundToImage.return_value = csmapi.ImageCoord(1.0, 1.0)
    second_node.camera.imagetoground.return_value = csmapi.EcefCoord(1.0,1.0,1.0)
    second_node.isis_serial = '2'
    second_node.__getitem__.return_value = 2
    third_node = MagicMock()
    third_node.camera = MagicMock()
    third_node.camera.groundToImage.return_value = csmapi.ImageCoord(0.0, 1.0)
    third_node.camera.imageToGround.return_value = csmapi.EcefCoord(1.0,1.0,1.0)
    third_node.isis_serial = '3'
    third_node.__getitem__.return_value = 3
    fourth_node = MagicMock()
    fourth_node.camera = MagicMock()
    fourth_node.camera.groundToImage.return_value = csmapi.ImageCoord(0.0, 0.0)
    first_node.camera.imageToGround.return_value = csmapi.EcefCoord(1.0,0,1.0)
    fourth_node.isis_serial = '4'
    fourth_node.__getitem__.return_value = 4

    # Actual function being tested
    points = place_points_in_overlap([first_node, second_node, third_node, fourth_node],
                                      Polygon([(0, 0), (0, 10), (10, 10), (10, 0)]))
                                    Polygon([(0, 0), (0, 10), (10, 10), (10, 0)]),
                                    cam_type='csm')

    # Check the function output
    assert len(points) == 3
@@ -43,18 +55,18 @@ def test_place_points_in_overlap(point_distributer):
        assert measure_serials == ['1', '2', '3', '4']

    # Check the mocks
    point_distributer.assert_called_with(Polygon([(0, 0), (0, 10), (10, 10), (10, 0)]))
    np.testing.assert_array_equal(point_distributer.call_args[0][0], np.array([0.0, 0.0, 0.0]))
    first_node.camera.groundToImage.assert_called()
    second_node.camera.groundToImage.assert_called()
    third_node.camera.groundToImage.assert_called()
    fourth_node.camera.groundToImage.assert_called()


class MockOverlap():
    intersections = [0,1]
    geom = Polygon([(0,0),(0,5),(5,5),(5,0),(0,0)])



@patch('autocnet.io.db.model.Overlay.overlapping_larger_than', return_value=[MockOverlap()]*3)
@patch('autocnet.io.db.model.Points.bulkadd')
@pytest.mark.parametrize("distributekwargs",[
@@ -63,8 +75,8 @@ class MockOverlap():
def test_place_points_in_overlaps(overlapper, adder, distributekwargs):
    nodes = [{"id": 0, "data": Node()}, {"id": 1, "data": Node()}]
    with patch('autocnet.spatial.overlap.place_points_in_overlap', return_value=[1,2,3]) as ppio:
        place_points_in_overlaps(nodes,distribute_points_kwargs=distributekwargs)
        place_points_in_overlaps(nodes, cam_type="isis", distribute_points_kwargs=distributekwargs)
        ppio.assert_called_with([Node(), Node()],
                                Polygon([(0,0),(0,5),(5,5),(5,0),(0,0)]),
                                distribute_points_kwargs=distributekwargs,
                                cam_type='csm')
                                cam_type='isis')
+49 −0
Original line number Diff line number Diff line
#!/usr/bin/env python

import copy
import os
import json
import sys
import warnings

from redis import StrictRedis
import yaml

from autocnet.io.db.redis_queue import pop_computetime_push, finalize
from autocnet.io.db.model import Overlay, Images, Measures
from autocnet.matcher.subpixel import subpixel_register_measure
from autocnet.graph.node import NetworkNode
from autocnet import Session, dem

from plio.io.io_gdal import GeoDataset

#Load the config file
try:
    with open(os.environ['autocnet_config'], 'r') as f:
        config = yaml.safe_load(f)
except:
    print("The 'autocnet_config' environment variable is not set.")
    sys.exit(1)

def main(msg, config):
    mid = msg['id']
    print(f'Subpixel registering measure {mid}')
    subpixel_register_measure(mid,
                              iterative_phase_kwargs=msg['iterative_phase_kwargs'],
                              subpixel_template_kwargs=msg['subpixel_template_kwargs'],
                              threshold=msg['threshold'],
                              **msg['cost_kwargs'])
    print('Subpixel registration complete.')

if __name__ == '__main__':
    conf = config['redis']
    queue = StrictRedis(host=conf['host'], port=conf['port'], db=0)

    msg = pop_computetime_push(queue,
                               conf['processing_queue'],
                               conf['working_queue'])
    if msg is None:
        warnings.warn('Expected to process a cluster job, but the message queue is empty.')
        sys.exit()

    main(msg, config)
Loading