Commit 7bd9c74c authored by Jesse Mapel's avatar Jesse Mapel Committed by Kelvin Rodriguez
Browse files

cluster overlap point seeding (#297) 

* Fixed install of new script

* more work on the overlap script

* In progress commit

* Now working

* Updates for new geometry handling in DB model

* Moved cluster overlaps into library

* Removed network stuff from place_points_in_overlap

* first pass at overlap test

* Fixed test and comments

* more updates for tests and comments

* Mocked out more

* In progress testing

* Finally test is working

* Updated tests and script

* fixed not passning dem in test

* Removed duplicate logic from serial place points
parent ab118e24

autocnet/spatial/overlap.py

+156 −70
Original line number Diff line number Diff line 
import warnings

from autocnet import config

from autocnet.cg import cg as compgeom

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

from autocnet.matcher.subpixel import iterative_phase

from autocnet import Session, engine

import json



import csmapi

import numpy as np

from redis import StrictRedis

import pyproj

import shapely

import sqlalchemy

from plio.io.io_gdal import GeoDataset



def place_points_in_overlaps(cg, size_threshold=0.0007, reference=None,

                             iterative_phase_kwargs={'size':71}):

from autocnet import config, Session, engine

from autocnet.cg import cg as compgeom

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

from autocnet.matcher.subpixel import iterative_phase

from plurmy import Slurm

import csmapi



# SQL query to decompose pairwise overlaps

compute_overlaps_sql = """

WITH intersectiongeom AS

(SELECT geom AS geom FROM ST_Dump((

   SELECT ST_Polygonize(the_geom) AS the_geom FROM (

     SELECT ST_Union(the_geom) AS the_geom FROM (

	   SELECT ST_ExteriorRing((ST_DUMP(footprint_latlon)).geom) AS the_geom

	     FROM images WHERE images.footprint_latlon IS NOT NULL) AS lines

	) AS noded_lines))),

iid AS (

 SELECT images.id, intersectiongeom.geom AS geom

		FROM images, intersectiongeom

		WHERE images.footprint_latlon is NOT NULL AND

		ST_INTERSECTS(intersectiongeom.geom, images.footprint_latlon) AND

		ST_AREA(ST_INTERSECTION(intersectiongeom.geom, images.footprint_latlon)) > 0.000001

)

INSERT INTO overlay(intersections, geom) SELECT row.intersections, row.geom FROM

(SELECT iid.geom, array_agg(iid.id) AS intersections

  FROM iid GROUP BY iid.geom) AS row WHERE array_length(intersections, 1) > 1;

"""

    Given a geometry, place points into the geometry by back-projecing using

    a sensor model.compgeom



    The DEM specified in the config file will be used to calculate height point elevations.

def place_points_in_overlaps(cg, size_threshold=0.0007,

                             iterative_phase_kwargs={'size':71}):

    """

    Place points in all of the overlap geometries by back-projecing using

    sensor models.



    TODO: This shoucompgeomn once that package is stable.

    The DEM specified in the config file will be used to calculate point elevations.



    Parameters

    ----------
@@ -30,9 +51,8 @@ def place_points_in_overlaps(cg, size_threshold=0.0007, reference=None, 
    size_threshold : float

                     overlaps with area <= this threshold are ignored



    reference : int

                the i.d. of a reference node to use when placing points. If not

                speficied, this is the node with the lowest id

    iterative_phase_kwargs : dict

        Dictionary of keyword arguments for the iterative phase matcher function

    """

    if not Session:

        warnings.warn('This function requires a database connection configured via an autocnet config file.')
@@ -40,49 +60,138 @@ def place_points_in_overlaps(cg, size_threshold=0.0007, reference=None, 


    points = []

    session = Session()

    srid = config['spatial']['srid']

    semi_major = config['spatial']['semimajor_rad']

    semi_minor = config['spatial']['semiminor_rad']

    ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)

    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)

    if 'dem' in config['spatial']:

        dem = config['spatial']['dem']

        gd = GeoDataset(dem)

    else:

        gd = None



    # TODO: This should be a passable query where we can subset.

    for o in session.query(Overlay).\

             filter(sqlalchemy.func.ST_Area(Overlay.geom) >= size_threshold):



        valid = compgeom.distribute_points_in_geom(o.geom)

        if not valid:

            continue



             filter(sqlalchemy.func.ST_Area(Overlay.geom) >= size_threshold).\

             filter(sqlalchemy.func.array_length(Overlay.intersections, 1) > 1):

        overlaps = o.intersections



        if overlaps == None:

            continue

        nodes = [cg.node[id] for id in overlaps]

        points.extend(place_points_in_overlap(nodes, o.geom, dem=gd,

                                              iterative_phase_kwargs=iterative_phase_kwargs))



        if reference is None:

            source = overlaps[0]

        else:

            source = reference

        overlaps.remove(source)

        source = cg.node[source]['data']

        source_camera = source.camera



    session.add_all(points)

    session.commit()



def cluster_place_points_in_overlaps(size_threshold=0.0007,

                                     iterative_phase_kwargs={'size':71},

                                     walltime='00:10:00'):

    """

    Place points in all of the overlap geometries by back-projecing using

    sensor models. This method uses the cluster to process all of the overlaps

    in parallel. See place_points_in_overlap and acn_overlaps.



    The DEM specified in the config file will be used to calculate point elevations.



    Parameters

    ----------

    size_threshold : float

        overlaps with area <= this threshold are ignored



    iterative_phase_kwargs : dict

        Dictionary of keyword arguments for the iterative phase matcher function



    walltime : str

        Cluster job wall time as a string HH:MM:SS

    """

    if not Session:

        warnings.warn('This function requires a database connection configured via an autocnet config file.')

        return



    # Get all of the overlaps over the size threshold

    session = Session()

    overlaps = session.query(Overlay.id, Overlay.geom, Overlay.intersections).\

                       filter(sqlalchemy.func.ST_Area(Overlay.geom) >= size_threshold).\

                       filter(sqlalchemy.func.array_length(Overlay.intersections, 1) > 1)

    session.close()



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

    for overlap in overlaps:

        msg = {'id' : overlap.id,

               'iterative_phase_kwargs' : iterative_phase_kwargs,

               'walltime' : walltime}

        rqueue.rpush(queuename, json.dumps(msg))

    job_counter = len([*overlaps]) + 1



    # Submit the jobs

    submitter = Slurm('acn_overlaps',

                 mem_per_cpu=config['cluster']['processing_memory'],

                 time=walltime,

                 partition=config['cluster']['queue'],

                 output=config['cluster']['cluster_log_dir']+'/slurm-%A_%a.out')

    submitter.submit(array='1-{}'.format(job_counter))

    return job_counter



def place_points_in_overlap(nodes, geom, dem=None,

                            iterative_phase_kwargs={'size':71}):

    """

    Place points into an overlap geometry by back-projecing using sensor models.



    Parameters

    ----------

    nodes : list of Nodes

        The CandidateGraph nodes of all the images that intersect the overlap



    geom : geometry

        The geometry of the overlap region



    dem : GeoDataset

         The DEM used to compute point elevations. An elevation of 0 is is used

         if no DEM is passed in.



    iterative_phase_kwargs : dict

        Dictionary of keyword arguments for the iterative phase matcher function



    Returns

    -------

    points : list of Points

        The list of points seeded in the overlap

    """

    points = []

    semi_major = config['spatial']['semimajor_rad']

    semi_minor = config['spatial']['semiminor_rad']

    ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)

    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)



    valid = compgeom.distribute_points_in_geom(geom)

    if not valid:

        raise ValueError('Failed to distribute points in overlap')



    # Grab the source image. This is just the node with the lowest ID, nothing smart.

    source = nodes[0]

    nodes.remove(source)

    source_camera = source.camera

    for v in valid:

            point = Points(geom=shapely.geometry.Point(*v),

        geom = shapely.geometry.Point(v[0], v[1])

        point = Points(geom=geom,

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



        # Calculate the height, the distance (in meters) above or

        # below the aeroid (meters above or below the BCBF spheroid).

            px, py = gd.latlon_to_pixel(v[1], v[0])

            height = gd.read_array(1, [px, py, 1, 1])[0][0]

        if dem is None:

            height = 0

        else:

            px, py = dem.latlon_to_pixel(v[1], v[0])

            height = dem.read_array(1, [px, py, 1, 1])[0][0]



        # Get the BCEF coordinate from the lon, lat

            x, y, z = pyproj.transform(lla, ecef, v[0], v[1], height)

        x, y, z = pyproj.transform(lla, ecef, v[0], v[1], height)  # -3000 working well in elysium, need aeroid

        gnd = csmapi.EcefCoord(x, y, z)



            # Grab the source image. This is just the node with the lowest ID, nothing smart.

        sic = source_camera.groundToImage(gnd)

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

                                       line=sic.line,
@@ -91,40 +200,17 @@ def place_points_in_overlaps(cg, size_threshold=0.0007, reference=None, 
                                       measuretype=3))





            for i, d in enumerate(overlaps):

                destination = cg.node[d]['data']

                destination_camera = destination.camera

                dic = destination_camera.groundToImage(gnd)

                dx, dy, metrics = iterative_phase(sic.samp, sic.line, dic.samp, dic.line,

                                                  source.geodata, destination.geodata,

        for i, dest in enumerate(nodes):

            dic = dest.camera.groundToImage(gnd)

            dx, dy, _ = iterative_phase(sic.samp, sic.line, dic.samp, dic.line,

                                        source.geodata, dest.geodata,

                                        **iterative_phase_kwargs)

            if dx is not None or dy is not None:

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

                                               line=dy,

                                                   imageid=destination['node_id'],

                                                   serial=destination.isis_serial,

                                               imageid=dest['node_id'],

                                               serial=dest.isis_serial,

                                               measuretype=3))

        if len(point.measures) >= 2:

            points.append(point)

    session.add_all(points)

    session.commit()



compute_overlaps_sql = """

WITH intersectiongeom AS

(SELECT geom AS geom FROM ST_Dump((

   SELECT ST_Polygonize(the_geom) AS the_geom FROM (

     SELECT ST_Union(the_geom) AS the_geom FROM (

	   SELECT ST_ExteriorRing((ST_DUMP(footprint_latlon)).geom) AS the_geom

	     FROM images WHERE images.footprint_latlon IS NOT NULL) AS lines

	) AS noded_lines))),

iid AS (

 SELECT images.id, intersectiongeom.geom AS geom

		FROM images, intersectiongeom

		WHERE images.footprint_latlon is NOT NULL AND

		ST_INTERSECTS(intersectiongeom.geom, images.footprint_latlon) AND

		ST_AREA(ST_INTERSECTION(intersectiongeom.geom, images.footprint_latlon)) > 0.000001

)

INSERT INTO overlay(intersections, geom) SELECT row.intersections, row.geom FROM 

(SELECT iid.geom, array_agg(iid.id) AS intersections

  FROM iid GROUP BY iid.geom) AS row WHERE array_length(intersections, 1) > 1;

"""
 
 No newline at end of file
 
    return points

autocnet/spatial/tests/__init__.py

0 → 100644
+0 −0

Empty file added.

autocnet/spatial/tests/test_overlap_points.py

0 → 100644
+60 −0
Original line number Diff line number Diff line 
import pytest

from unittest.mock import MagicMock, patch

from shapely.geometry import Polygon

from autocnet.spatial.overlap import place_points_in_overlap

import csmapi



@patch('autocnet.spatial.overlap.iterative_phase', return_value=(0, 1, 2))

@patch('autocnet.cg.cg.distribute_points_in_geom', return_value=[(0, 0), (5, 5), (10, 10)])

def test_place_points_in_overlap(point_distributer, phase_matcher):

    # Mock setup

    first_node = MagicMock()

    first_node.camera = MagicMock()

    first_node.camera.groundToImage.return_value = csmapi.ImageCoord(1.0, 0.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.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.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)

    fourth_node.isis_serial = '4'

    fourth_node.__getitem__.return_value = 4

    dem = MagicMock()

    dem.latlon_to_pixel.return_value = (1.0, 1.0)

    dem.read_array.return_value = [[0.0]]



    # 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)]), dem)



    # Check the function output

    assert len(points) == 3

    for point in points:

        measure_ids = [measure.imageid for measure in point.measures]

        measure_serials = [measure.serial for measure in point.measures]

        assert measure_ids == [1, 2, 3, 4]

        assert measure_serials == ['1', '2', '3', '4']



    # Check the mocks

    point_distributer.assert_called_with(Polygon([(0, 0), (0, 10), (10, 10), (10, 0)]))

    dem.latlon_to_pixel.assert_called()

    dem.read_array.assert_called()

    first_node.camera.groundToImage.assert_called()

    second_node.camera.groundToImage.assert_called()

    third_node.camera.groundToImage.assert_called()

    fourth_node.camera.groundToImage.assert_called()

    phase_matcher.assert_any_call(0.0, 1.0, 1.0, 1.0,

                                  first_node.geodata, second_node.geodata, size=71)

    phase_matcher.assert_any_call(0.0, 1.0, 1.0, 0.0,

                                  first_node.geodata, third_node.geodata, size=71)

    phase_matcher.assert_any_call(0.0, 1.0, 0.0, 0.0,

                                  first_node.geodata, fourth_node.geodata, size=71)

bin/acn_overlaps

100755 → 100644
+37 −5
Original line number Diff line number Diff line
@@ -8,8 +8,11 @@ import sys 
from redis import StrictRedis

import yaml



from autocnet.io.db import connection

from autocnet.io.db.redis_queue import pop_computetime_push, finalize

from autocnet.io.db.model import Overlay, Images

from autocnet.spatial.overlap import place_points_in_overlap

from autocnet.graph.node import NetworkNode

from autocnet import Session



#Load the config file

try:
@@ -17,10 +20,39 @@ try: 
        config = yaml.load(f)

except:

    print("The 'autocnet_config' environment variable is not set.")

    sys.exit()

    sys.exit(1)



def main():

    print('hello')

def main(msg, config):

    session = Session()

    id = msg['id']

    res = session.query(Overlay).filter(Overlay.id == msg['id'])

    if res is None:

        print('Could not find overlap with ID', id)

        sys.exit(1)

    overlap = res.first()

    geom = overlap.geom

    nodes = []

    for id in overlap.intersections:

        res = session.query(Images).filter(Images.id == id).first()

        nodes.append(NetworkNode(node_id=id, image_path=res.path))

    if 'dem' in config['spatial']:

        dem = config['spatial']['dem']

        gd = GeoDataset(dem)

    else:

        gd = None

    print('Placing points in overlap', id)

    points = place_points_in_overlap(nodes, overlap.geom, gd,

                                     msg['iterative_phase_kwargs'])

    session.add_all(points)

    session.commit()

    session.close()



if __name__ == '__main__':

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



    main(msg, config)

config/test_config.yml

+1 −1
Original line number Diff line number Diff line
@@ -45,6 +45,7 @@ spatial: 
    semimajor_rad: 3396190  # in meters

    semiminor_rad: 3376200  # in meters

    proj4_str: '+proj:longlat +a:3396190 +b:3376200 +no_defs'

    dem: '/test/dem.img'



### Working Directories ###

directories:
@@ -77,4 +78,3 @@ algorithms: 
          initial_y_size: 500

          corr_x_size: 40

          corr_y_size: 40